Add Convergent Source

pytissueoptics/rayscattering/__init__.py

@@ -18,7 +18,7 @@
 from .opencl import CONFIG, disableOpenCL, hardwareAccelerationIsAvailable
 from .photon import Photon
 from .scatteringScene import ScatteringScene
-from .source import DirectionalSource, DivergentSource, IsotropicPointSource, PencilPointSource
+from .source import DirectionalSource, DivergentSource, ConvergentSource, IsotropicPointSource, PencilPointSource
 from .statistics import Stats
 
 __all__ = [
@@ -28,6 +28,7 @@
     "IsotropicPointSource",
     "DirectionalSource",
     "DivergentSource",
+    "ConvergentSource",
     "EnergyLogger",
     "EnergyType",
     "ScatteringScene",

pytissueoptics/rayscattering/source.py

@@ -335,3 +335,42 @@ def _getInitialDirections(self):
     @property
     def _hashComponents(self) -> tuple:
         return self._position, self._direction, self._diameter, self._divergence
+
+
+class ConvergentSource(DirectionalSource):
+    def __init__(
+        self,
+        position: Vector,
+        direction: Vector,
+        diameter: float,
+        focalLength: float,
+        N: int,
+        useHardwareAcceleration: bool = True,
+        displaySize: float = 0.1,
+        seed: Optional[int] = None,
+    ):
+        if focalLength <= 0:
+            raise ValueError("The focal length of a convergent source must be positive.")
+
+        self._focalLength = focalLength
+
+        super().__init__(
+            position=position,
+            direction=direction,
+            diameter=diameter,
+            N=N,
+            useHardwareAcceleration=useHardwareAcceleration,
+            displaySize=displaySize,
+            seed=seed,
+        )
+
+    def getInitialPositionsAndDirections(self) -> Tuple[np.ndarray, np.ndarray]:
+        positions = self._getInitialPositions()
+        focalPoint = self._position + self._direction * self._focalLength
+        directions = focalPoint.array - positions
+        directions /= np.linalg.norm(directions, axis=1, keepdims=True)
+        return positions, directions
+
+    @property
+    def _hashComponents(self) -> tuple:
+        return self._position, self._direction, self._diameter, self._focalLength

pytissueoptics/rayscattering/tests/testSource.py

@@ -8,7 +8,13 @@
 from pytissueoptics.rayscattering import EnergyLogger, PencilPointSource, Photon
 from pytissueoptics.rayscattering.materials import ScatteringMaterial
 from pytissueoptics.rayscattering.scatteringScene import ScatteringScene
-from pytissueoptics.rayscattering.source import DirectionalSource, DivergentSource, IsotropicPointSource, Source
+from pytissueoptics.rayscattering.source import (
+    DirectionalSource,
+    DivergentSource,
+    IsotropicPointSource,
+    Source,
+    ConvergentSource,
+)
 from pytissueoptics.scene.geometry import Environment, Vector
 from pytissueoptics.scene.logger import Logger
 from pytissueoptics.scene.solids import Solid
@@ -252,3 +258,118 @@ def testGivenTwoDivergentSourcesThatDifferInDivergence_shouldNotHaveSameHash(sel
             position=Vector(), direction=sourceDirection, diameter=1, divergence=divergence2, N=1
         )
         self.assertNotEqual(hash(divergentSource1), hash(divergentSource2))
+
+
+class TestConvergentSource(unittest.TestCase):
+    def testGivenNegativeOrZeroFocalLength_shouldRaiseValueError(self):
+        with self.assertRaises(ValueError):
+            ConvergentSource(
+                position=Vector(),
+                direction=Vector(0, 0, 1),
+                focalLength=0,
+                diameter=1,
+                N=1,
+                useHardwareAcceleration=False,
+            )
+        with self.assertRaises(ValueError):
+            ConvergentSource(
+                position=Vector(),
+                direction=Vector(0, 0, 1),
+                focalLength=-1,
+                diameter=1,
+                N=1,
+                useHardwareAcceleration=False,
+            )
+
+    def testShouldHavePhotonsPointingTowardTheFocalPoint(self):
+        np.random.seed(0)
+        position = Vector(0, 0, 0)
+        direction = Vector(0, 0, 1)
+        focalLength = 5.0
+        diameter = 2.0
+        source = ConvergentSource(
+            position=position,
+            direction=direction,
+            focalLength=focalLength,
+            diameter=diameter,
+            N=10,
+            useHardwareAcceleration=False,
+        )
+
+        focalPoint = position + direction * focalLength
+        for photon in source.photons:
+            expectedDirection = focalPoint - photon.position
+            expectedDirection.normalize()
+            self.assertEqual(expectedDirection, photon.direction)
+
+    def testGivenInfiniteFocalLength_shouldHavePhotonsAllPointingInTheSourceDirection(self):
+        sourceDirection = Vector(1, 0, 0)
+        source = ConvergentSource(
+            position=Vector(),
+            direction=sourceDirection,
+            focalLength=1e10,
+            diameter=1.0,
+            N=10,
+            useHardwareAcceleration=False,
+        )
+        for photon in source.photons:
+            self.assertEqual(sourceDirection, photon.direction)
+
+    def testShouldHavePhotonsUniformlyPositionedInsideTheSourceDiameter(self):
+        np.random.seed(0)
+        sourcePosition = Vector(3, 3, 0)
+        sourceDiameter = 2.0
+        source = ConvergentSource(
+            position=sourcePosition,
+            direction=Vector(0, 1, 0),
+            focalLength=5.0,
+            diameter=sourceDiameter,
+            N=10,
+            useHardwareAcceleration=False,
+        )
+        for photon in source.photons:
+            self.assertTrue(np.isclose(photon.position.y, sourcePosition.y))
+            self.assertTrue(
+                sourcePosition.x - sourceDiameter / 2 <= photon.position.x <= sourcePosition.x + sourceDiameter / 2
+            )
+            self.assertTrue(
+                sourcePosition.z - sourceDiameter / 2 <= photon.position.z <= sourcePosition.z + sourceDiameter / 2
+            )
+
+    def testGivenTwoConvergentSourcesWithSamePropertiesExceptPhotonCount_shouldHaveSameHash(self):
+        source1 = ConvergentSource(
+            position=Vector(),
+            direction=Vector(0, 0, 1),
+            focalLength=5.0,
+            diameter=1.0,
+            N=1,
+            useHardwareAcceleration=False,
+        )
+        source2 = ConvergentSource(
+            position=Vector(),
+            direction=Vector(0, 0, 1),
+            focalLength=5.0,
+            diameter=1.0,
+            N=2,
+            useHardwareAcceleration=False,
+        )
+        self.assertEqual(hash(source1), hash(source2))
+
+    def testGivenTwoConvergentSourcesThatDifferInFocalLength_shouldNotHaveSameHash(self):
+        source1 = ConvergentSource(
+            position=Vector(),
+            direction=Vector(0, 0, 1),
+            focalLength=5.0,
+            diameter=1.0,
+            N=1,
+            useHardwareAcceleration=False,
+        )
+        source2 = ConvergentSource(
+            position=Vector(),
+            direction=Vector(0, 0, 1),
+            focalLength=10.0,
+            diameter=1.0,
+            N=1,
+            useHardwareAcceleration=False,
+        )
+        self.assertNotEqual(hash(source1), hash(source2))