From 4c86338df97da95cc59ad9e816410e26dd32e61e Mon Sep 17 00:00:00 2001 From: jlbegin Date: Sun, 22 Jun 2025 14:14:05 -0400 Subject: [PATCH 01/14] add rectangle and circle solids --- pytissueoptics/scene/solids/circle.py | 27 ++++++++++++++++++ pytissueoptics/scene/solids/rectangle.py | 35 ++++++++++++++++++++++++ pytissueoptics/scene/solids/solid.py | 4 ++- 3 files changed, 65 insertions(+), 1 deletion(-) create mode 100644 pytissueoptics/scene/solids/circle.py create mode 100644 pytissueoptics/scene/solids/rectangle.py diff --git a/pytissueoptics/scene/solids/circle.py b/pytissueoptics/scene/solids/circle.py new file mode 100644 index 00000000..e42e6fd6 --- /dev/null +++ b/pytissueoptics/scene/solids/circle.py @@ -0,0 +1,27 @@ +from ..geometry import primitives, Vector +from .cylinder import Cylinder + + +class Circle(Cylinder): + def __init__(self, radius: float, orientation: Vector = Vector(0, 0, 1), position: Vector = Vector(0, 0, 0), + u: int = 32, s: int = 1, primitive: str = primitives.DEFAULT): + super().__init__(radius, length=0, u=u, v=1, s=s, position=position, label="circle", primitive=primitive, smooth=False) + self.orient(towards=orientation) + + def _computeTriangleMesh(self): + # A circle is a cylinder of length 0 with a single (back) surface pointing towards the positive z axis before orientation. + lateralLayers = self._computeSectionVertices(lateralSteps=[0], radialSteps=[self._s]) + backLayers = self._computeSectionVertices(lateralSteps=[0], radialSteps=list(range(self._s - 1, 0, -1))) + backLayers.insert(0, lateralLayers[-1]) + backLayers.append([self._backCenter]) + self._vertices.append(self._backCenter) + self._surfaces.add("surface", self._getSurfaceTriangles(backLayers)) + + def _geometryParams(self) -> dict: + return { + "radius": self._radius, + "position": self._position, + "orientation": self._orientation, + "u": self._u, + "s": self._s + } diff --git a/pytissueoptics/scene/solids/rectangle.py b/pytissueoptics/scene/solids/rectangle.py new file mode 100644 index 00000000..03ba56cb --- /dev/null +++ b/pytissueoptics/scene/solids/rectangle.py @@ -0,0 +1,35 @@ +from ..geometry import primitives, Vertex, Triangle, Quad, Vector +from .solid import Solid + + +class Rectangle(Solid): + def __init__(self, a: float, b: float, orientation = Vector(0, 0, 1), position = Vector(0, 0, 0), + primitive: str = primitives.DEFAULT): + self._a = a + self._b = b + + vertices = [ # Mesh points towards (0, 0, 1) by default, as expected by the orient() method. + Vertex(-a / 2, -b / 2, 0), + Vertex(a / 2, -b / 2, 0), + Vertex(a / 2, b / 2, 0), + Vertex(-a / 2, b / 2, 0) + ] + + super().__init__(vertices, position, label="rectangle", primitive=primitive) + self.orient(towards=orientation) + + def _computeTriangleMesh(self): + V = self._vertices + self._surfaces.add("surface", [Triangle(V[0], V[1], V[2]), Triangle(V[2], V[3], V[0])]) + + def _computeQuadMesh(self): + V = self._vertices + self._surfaces.add("surface", [Quad(V[0], V[1], V[2], V[3])]) + + def _geometryParams(self) -> dict: + return { + "a": self._a, + "b": self._b, + "position": self._position, + "orientation": self._orientation + } diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py index 007b9fd0..e71a91e5 100644 --- a/pytissueoptics/scene/solids/solid.py +++ b/pytissueoptics/scene/solids/solid.py @@ -348,4 +348,6 @@ def geometryExport(self) -> dict[str, Any]: def _geometryParams(self) -> dict: """To be implemented by Solid subclasses to detail other geometry parameters.""" - return {} + raise NotImplementedError( + f"Geometry parameters not implemented for Solids of type {type(self).__name__}" + ) From 97c9c0adf329f12081e294bff06a69458c580898 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Fri, 10 Oct 2025 16:58:43 -0400 Subject: [PATCH 02/14] force vtk 9.4 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 4ba54b80..17b1bc35 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,7 +27,7 @@ dependencies = [ "Pygments", "siphash24; python_version < '3.13'", "pyopencl", - "vtk>=9.4", + "vtk==9.4", "mayavi-dev", "pyqt5", ] From c9f7e07d03fa913e67ac1597e512605961b17a56 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Fri, 10 Oct 2025 17:00:51 -0400 Subject: [PATCH 03/14] fix GetAxisAngleBetween for opposite case --- pytissueoptics/scene/geometry/utils.py | 8 +++++++- pytissueoptics/scene/tests/geometry/testUtils.py | 16 +++++++++++----- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/pytissueoptics/scene/geometry/utils.py b/pytissueoptics/scene/geometry/utils.py index 0d61d4c7..6bdaf8ca 100644 --- a/pytissueoptics/scene/geometry/utils.py +++ b/pytissueoptics/scene/geometry/utils.py @@ -55,7 +55,13 @@ def getAxisAngleBetween(fromDirection: Vector, toDirection: Vector) -> Tuple[Vec dot = fromDirection.dot(toDirection) dot = max(min(dot, 1), -1) angle = np.arccos(dot) - axis = fromDirection.cross(toDirection) + + # Special case: vectors are opposite (antiparallel) + if np.isclose(dot, -1): + axis = fromDirection.getAnyOrthogonal() + else: + axis = fromDirection.cross(toDirection) + axis.normalize() axis = axis.array return Vector(*axis), angle diff --git a/pytissueoptics/scene/tests/geometry/testUtils.py b/pytissueoptics/scene/tests/geometry/testUtils.py index ffa2c41b..f88d06cf 100644 --- a/pytissueoptics/scene/tests/geometry/testUtils.py +++ b/pytissueoptics/scene/tests/geometry/testUtils.py @@ -94,13 +94,12 @@ def testGivenTwoVectorsOnABaseAxisPlaneWithLessThan180DegreesBetween_shouldRetur testCases = self._generateAxisAlignedTestCases() for testCase in testCases: axis, angle = getAxisAngleBetween(testCase.v1, testCase.v2) - self.assertTrue(np.allclose(testCase.expectedAxis.array, axis.array)) self.assertTrue(np.isclose(testCase.expectedAngle, angle)) + self.assertTrue(np.isclose(1, axis.getNorm())) - axis, angle = getAxisAngleBetween(testCase.v2, testCase.v1) - testCase.expectedAxis *= -1 - self.assertTrue(np.allclose(testCase.expectedAxis.array, axis.array)) - self.assertTrue(np.isclose(testCase.expectedAngle, angle)) + # For non-degenerate cases (angle not 0 or 180), check axis direction + if not np.isclose(testCase.expectedAngle % np.pi, 0): + self.assertTrue(np.allclose(testCase.expectedAxis.array, axis.array)) def _generateAxisAlignedTestCases(self): testCases = [] @@ -126,3 +125,10 @@ def testGivenTwoVectorsWithSameOrientation_shouldReturnZeroAngle(self): v2 = Vector(1, 0, 0) _, angle = getAxisAngleBetween(v1, v2) self.assertTrue(np.isclose(0, angle)) + + def testGivenTwoOppositeVectors_shouldReturnAnyAxisAnd180Degrees(self): + v1 = Vector(1, 0, 0) + v2 = Vector(-1, 0, 0) + axis, angle = getAxisAngleBetween(v1, v2) + self.assertTrue(np.isclose(np.pi, angle)) + self.assertTrue(np.isclose(1, axis.getNorm())) From a24b22343a818d31235ca16668ec7269141aabc4 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Fri, 10 Oct 2025 17:05:45 -0400 Subject: [PATCH 04/14] define detector through Solid.asDetector(halfAngle) --- pytissueoptics/scene/__init__.py | 4 +++ pytissueoptics/scene/solids/__init__.py | 4 +++ pytissueoptics/scene/solids/circle.py | 4 +-- pytissueoptics/scene/solids/rectangle.py | 4 +-- pytissueoptics/scene/solids/solid.py | 31 ++++++++++++++++++++++++ 5 files changed, 43 insertions(+), 4 deletions(-) diff --git a/pytissueoptics/scene/__init__.py b/pytissueoptics/scene/__init__.py index c3edf9b6..bc7a54cf 100644 --- a/pytissueoptics/scene/__init__.py +++ b/pytissueoptics/scene/__init__.py @@ -14,6 +14,8 @@ Sphere, SymmetricLens, ThickLens, + Rectangle, + Circle, ) from .viewer import ViewPointStyle, get3DViewer @@ -37,4 +39,6 @@ "SymmetricLens", "PlanoConvexLens", "PlanoConcaveLens", + "Rectangle", + "Circle", ] diff --git a/pytissueoptics/scene/solids/__init__.py b/pytissueoptics/scene/solids/__init__.py index 9e63b939..bb39b960 100644 --- a/pytissueoptics/scene/solids/__init__.py +++ b/pytissueoptics/scene/solids/__init__.py @@ -7,6 +7,8 @@ from .solid import Solid from .solidFactory import SolidFactory from .sphere import Sphere +from .circle import Circle +from .rectangle import Rectangle __all__ = [ "Solid", @@ -21,4 +23,6 @@ "PlanoConcaveLens", "ThickLens", "SymmetricLens", + "Rectangle", + "Circle", ] diff --git a/pytissueoptics/scene/solids/circle.py b/pytissueoptics/scene/solids/circle.py index e42e6fd6..26b4ae32 100644 --- a/pytissueoptics/scene/solids/circle.py +++ b/pytissueoptics/scene/solids/circle.py @@ -4,8 +4,8 @@ class Circle(Cylinder): def __init__(self, radius: float, orientation: Vector = Vector(0, 0, 1), position: Vector = Vector(0, 0, 0), - u: int = 32, s: int = 1, primitive: str = primitives.DEFAULT): - super().__init__(radius, length=0, u=u, v=1, s=s, position=position, label="circle", primitive=primitive, smooth=False) + u: int = 32, s: int = 1, primitive: str = primitives.DEFAULT, label="circle"): + super().__init__(radius, length=0, u=u, v=1, s=s, position=position, label=label, primitive=primitive, smooth=False) self.orient(towards=orientation) def _computeTriangleMesh(self): diff --git a/pytissueoptics/scene/solids/rectangle.py b/pytissueoptics/scene/solids/rectangle.py index 03ba56cb..3158a259 100644 --- a/pytissueoptics/scene/solids/rectangle.py +++ b/pytissueoptics/scene/solids/rectangle.py @@ -4,7 +4,7 @@ class Rectangle(Solid): def __init__(self, a: float, b: float, orientation = Vector(0, 0, 1), position = Vector(0, 0, 0), - primitive: str = primitives.DEFAULT): + primitive: str = primitives.DEFAULT, label="rectangle"): self._a = a self._b = b @@ -15,7 +15,7 @@ def __init__(self, a: float, b: float, orientation = Vector(0, 0, 1), position = Vertex(-a / 2, b / 2, 0) ] - super().__init__(vertices, position, label="rectangle", primitive=primitive) + super().__init__(vertices, position, label=label, primitive=primitive) self.orient(towards=orientation) def _computeTriangleMesh(self): diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py index e71a91e5..f2cde0b7 100644 --- a/pytissueoptics/scene/solids/solid.py +++ b/pytissueoptics/scene/solids/solid.py @@ -42,6 +42,7 @@ def __init__( self._bbox = None self._label = label self._layerLabels = {} + self._detectorAcceptanceCosine = None if not self._surfaces: self._computeMesh() @@ -80,6 +81,36 @@ def primitive(self) -> str: def bbox(self) -> BoundingBox: return self._bbox + def asDetector(self, halfAngle: float) -> "Solid": + """Treat this solid as a detector with a given half angle in radians. + + Detectors will fully absorb a photon when incident within the half angle, else the photon will go through + unaffected. Enabling this makes any previous material assigned to the solid irrelevant. Note that using a half + angle above pi/2 will detect rays coming from the back as well. + """ + if halfAngle < 0 or halfAngle > np.pi: + raise ValueError("Detector half angle must be between 0 and pi radians.") + + if self._material: + warnings.warn( + f"Solid '{self._label}' is now treated as a detector. Its material '{self._material}' will be ignored.", + RuntimeWarning, + ) + self._material = None + + self._detectorAcceptanceCosine = np.cos(halfAngle) + return self + + @property + def isDetector(self) -> bool: + return self._detectorAcceptanceCosine is not None + + @property + def detectorAcceptanceCosine(self) -> float: + if self._detectorAcceptanceCosine is None: + raise RuntimeError("This solid is not a detector. Set the half angle first with asDetector().") + return self._detectorAcceptanceCosine + def getBoundingBox(self) -> BoundingBox: return self.bbox From 2625bcccf23a2c430caab9e970ceccb3aa3401c7 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Fri, 10 Oct 2025 17:18:38 -0400 Subject: [PATCH 05/14] (python) Implement detectors in propagation logic. Pure sink (inside NA) without tracing for now. Pass-through outside NA. --- pytissueoptics/rayscattering/photon.py | 34 ++++++- .../rayscattering/scatteringScene.py | 2 +- .../rayscattering/statistics/statistics.py | 14 ++- .../rayscattering/tests/testPhoton.py | 89 +++++++++++++++++++ .../scene/intersection/intersectionFinder.py | 35 +++++--- pytissueoptics/scene/solids/solid.py | 2 +- 6 files changed, 158 insertions(+), 18 deletions(-) diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py index 523d5497..9a8bcd61 100644 --- a/pytissueoptics/rayscattering/photon.py +++ b/pytissueoptics/rayscattering/photon.py @@ -29,6 +29,8 @@ def __init__(self, position: Vector, direction: Vector): self._intersectionFinder: Optional[IntersectionFinder] = None self._logger: Optional[Logger] = None + self._lastIntersectedDetector: Optional[str] = None # Prevent intersecting with a detector when passing through it. + @property def isAlive(self) -> bool: return self._weight > 0 @@ -83,10 +85,23 @@ def step(self, distance=0) -> float: distance = 0 intersection = self._getIntersection(distance) + self._lastIntersectedDetector = None # Reset ignored intersection label after each step. if intersection: self.moveTo(intersection.position) - distanceLeft = self.reflectOrRefract(intersection) + isDetector = intersection.insideEnvironment.solid.isDetector + + if isDetector: + if self.detectOrIgnore(intersection): + return 0 # Skip unnecessary vertex check if detected. + + # Prevent re-intersecting with the same detector when passing through it. + self._lastIntersectedDetector = intersection.insideEnvironment.solid.getLabel() + + # skipping vertex check for now. TODO? Unclear stepCorrection and implications. + return intersection.distanceLeft + else: + distanceLeft = self.reflectOrRefract(intersection) # Check if intersection lies too close to a vertex. for vertex in intersection.polygon.vertices: @@ -118,7 +133,22 @@ def _getIntersection(self, distance) -> Optional[Intersection]: return None stepRay = Ray(self._position, self._direction, distance) - return self._intersectionFinder.findIntersection(stepRay, self.solidLabel) + return self._intersectionFinder.findIntersection(stepRay, self.solidLabel, self._lastIntersectedDetector) + + def detectOrIgnore(self, intersection: Intersection) -> bool: + # If the incidence angle is within the numerical aperture, absorb photon. + cosIncidence = -intersection.normal.dot(self._direction) + cosDetector = intersection.insideEnvironment.solid.detectorAcceptanceCosine + if cosIncidence >= cosDetector: + self._detectAndLog(intersection.insideEnvironment.solidLabel) + return True + return False + + def _detectAndLog(self, solidLabel: str): + if self._logger: + key = InteractionKey(solidLabel) + self._logger.logDataPoint(self._weight, self._position, key) + self._weight = 0 def reflectOrRefract(self, intersection: Intersection): fresnelIntersection = self._getFresnelIntersection(intersection) diff --git a/pytissueoptics/rayscattering/scatteringScene.py b/pytissueoptics/rayscattering/scatteringScene.py index 2cb05cbe..17cd5ce2 100644 --- a/pytissueoptics/rayscattering/scatteringScene.py +++ b/pytissueoptics/rayscattering/scatteringScene.py @@ -14,7 +14,7 @@ def __init__(self, solids: List[Solid], worldMaterial=ScatteringMaterial(), igno def add(self, solid: Solid, position: Vector = None): polygonSample = solid.getPolygons()[0] - if not isinstance(polygonSample.insideEnvironment.material, ScatteringMaterial): + if not isinstance(polygonSample.insideEnvironment.material, ScatteringMaterial) and not solid.isDetector: raise Exception( f"Solid '{solid.getLabel()}' has no ScatteringMaterial defined. " f"This is required for any RayScatteringScene. " diff --git a/pytissueoptics/rayscattering/statistics/statistics.py b/pytissueoptics/rayscattering/statistics/statistics.py index 5416535b..354d743d 100644 --- a/pytissueoptics/rayscattering/statistics/statistics.py +++ b/pytissueoptics/rayscattering/statistics/statistics.py @@ -1,3 +1,4 @@ +import math import os from dataclasses import dataclass from typing import Dict @@ -88,6 +89,11 @@ def _reportSolid(self, solidLabel: str): solidStats = self._solidStatsMap[solidLabel] reportString = "Report of solid '{}'\n".format(solidLabel) + if solidStats.absorbance == math.inf: + # Detectors don't log energy input. + reportString += " Detected {:.2f}% of total power\n".format(solidStats.totalAbsorbance) + return reportString + if solidStats.absorbance is None: reportString += " Absorbance: N/A ({:.2f}% of total power)\n".format(solidStats.totalAbsorbance) reportString += " Absorbance + Transmittance: N/A\n" @@ -110,12 +116,12 @@ def getAbsorbance(self, solidLabel: str, useTotalEnergy=False) -> float: return self._getAbsorbanceFromViews(solidLabel, useTotalEnergy) points = self._getPointCloud(solidLabel).solidPoints energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount() - return 100 * self._sumEnergy(points) / energyInput + return 100 * self._sumEnergy(points) / energyInput if energyInput else math.inf def _getAbsorbanceFromViews(self, solidLabel: str, useTotalEnergy=False) -> float: energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount() absorbedEnergy = self._getAbsorbedEnergyFromViews(solidLabel) - return 100 * absorbedEnergy / energyInput + return 100 * absorbedEnergy / energyInput if energyInput else math.inf def _getAbsorbedEnergyFromViews(self, solidLabel: str) -> float: for view in self._logger.views: @@ -207,7 +213,7 @@ def getTransmittance(self, solidLabel: str, surfaceLabel: str = None, useTotalEn points = self._getPointCloud(solidLabel, surfaceLabel).leavingSurfacePoints energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount() - return 100 * self._sumEnergy(points) / energyInput + return 100 * self._sumEnergy(points) / energyInput if energyInput else math.inf def _getTransmittanceFromViews(self, solidLabel: str, surfaceLabel: str = None, useTotalEnergy=False): if surfaceLabel is None: @@ -216,7 +222,7 @@ def _getTransmittanceFromViews(self, solidLabel: str, surfaceLabel: str = None, energyLeaving = self._getSurfaceEnergyFromViews(solidLabel, surfaceLabel, leaving=True) energyInput = self.getEnergyInput(solidLabel) if not useTotalEnergy else self.getPhotonCount() - return 100 * energyLeaving / energyInput + return 100 * energyLeaving / energyInput if energyInput else math.inf @staticmethod def _sumEnergy(points: np.ndarray): diff --git a/pytissueoptics/rayscattering/tests/testPhoton.py b/pytissueoptics/rayscattering/tests/testPhoton.py index 592701ac..974b3220 100644 --- a/pytissueoptics/rayscattering/tests/testPhoton.py +++ b/pytissueoptics/rayscattering/tests/testPhoton.py @@ -31,8 +31,10 @@ def setUp(self): self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy()) self.solidInside = mock(Solid) + self.solidInside.isDetector = False when(self.solidInside).getLabel().thenReturn(self.SOLID_INSIDE_LABEL) self.solidOutside = mock(Solid) + self.solidOutside.isDetector = False when(self.solidOutside).getLabel().thenReturn(self.SOLID_OUTSIDE_LABEL) def testShouldBeInTheGivenState(self): @@ -421,6 +423,93 @@ def testWhenInteractWithWeightAtFloatLimit_shouldKillPhoton(self): self.photon.interact() self.assertFalse(self.photon.isAlive) + def testWhenIntersectingDetectorInsideNA_shouldKillPhoton(self): + # Mock detector intersection + halfAngle = math.pi / 8 # 22.5 degrees + detector = mock(Solid) + detector.isDetector = True + detector.detectorAcceptanceCosine = math.cos(halfAngle) + when(detector).getLabel().thenReturn("detector") + self.solidInside = self.solidOutside = detector + + # Setup photon going towards detector within NA. + photonIncidence = halfAngle + self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence)) + self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy()) + + # Setup intersection with detector. + distance = 8 + intersectionFinder = self._createIntersectionFinder(distance) + self.photon.setContext( + Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder + ) + + self.photon.step(distance + 2) + + expectedIntersectionPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance + self.assertVectorEqual(expectedIntersectionPosition, self.photon.position) + self.assertFalse(self.photon.isAlive) + + def testWhenIntersectingDetectorOutsideNA_shouldNotKillPhoton(self): + # Mock detector intersection + halfAngle = math.pi / 8 # 22.5 degrees + detector = mock(Solid) + detector.isDetector = True + detector.detectorAcceptanceCosine = math.cos(halfAngle) + when(detector).getLabel().thenReturn("detector") + self.solidInside = self.solidOutside = detector + + # Setup photon going towards detector just outside NA. + photonIncidence = halfAngle + math.pi / 100 # Just outside NA + self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence)) + self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy()) + + # Setup intersection with detector. + distance = 8 + intersectionFinder = self._createIntersectionFinder(distance) + self.photon.setContext( + Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder + ) + + self.photon.step(distance + 2) + + expectedIntersectionPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance + self.assertVectorEqual(expectedIntersectionPosition, self.photon.position) + self.assertTrue(self.photon.isAlive) + self.assertEqual(self.photon.solidLabel, WORLD_LABEL) # Still in world + + def testGivenALogger_whenDetected_shouldLogWeightAtThisPositionWithDetectorLabel(self): + # Mock detector intersection + detectorLabel = "detector" + halfAngle = math.pi / 8 # 22.5 degrees + detector = mock(Solid) + detector.isDetector = True + detector.detectorAcceptanceCosine = math.cos(halfAngle) + when(detector).getLabel().thenReturn(detectorLabel) + self.solidInside = self.solidOutside = detector + + # Setup photon going towards detector within NA. + photonIncidence = halfAngle # Within NA + self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence)) + self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy()) + initial_weight = self.photon.weight + + # Setup intersection with detector. + distance = 8 + intersectionFinder = self._createIntersectionFinder(distance) + logger = self._createLogger() + self.photon.setContext( + Environment(ScatteringMaterial()), + intersectionFinder=intersectionFinder, + logger=logger, + ) + + self.photon.step(distance + 2) + + self.assertFalse(self.photon.isAlive) + intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance + verify(logger).logDataPoint(initial_weight, intersectionPoint, InteractionKey(detectorLabel)) + def assertVectorEqual(self, v1, v2): self.assertAlmostEqual(v1.x, v2.x, places=7) self.assertAlmostEqual(v1.y, v2.y, places=7) diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py index 2be19ec5..2a9af76c 100644 --- a/pytissueoptics/scene/intersection/intersectionFinder.py +++ b/pytissueoptics/scene/intersection/intersectionFinder.py @@ -35,11 +35,13 @@ def __init__(self, scene: Scene): self._polygonIntersect = MollerTrumboreIntersect() self._boxIntersect = GemsBoxIntersect() - def findIntersection(self, ray: Ray, currentSolidLabel: Optional[str]) -> Optional[Intersection]: + def findIntersection( + self, ray: Ray, currentSolidLabel: Optional[str], ignoreLabel: Optional[str] = None + ) -> Optional[Intersection]: raise NotImplementedError def _findClosestPolygonIntersection( - self, ray: Ray, polygons: List[Polygon], currentSolidLabel: str + self, ray: Ray, polygons: List[Polygon], currentSolidLabel: str, ignoreLabel: Optional[str] = None ) -> Optional[Intersection]: closestIntersection = Intersection(sys.maxsize) minSameSolidDistance = -sys.maxsize @@ -52,6 +54,9 @@ def _findClosestPolygonIntersection( if insideLabel != currentSolidLabel and outsideLabel != currentSolidLabel: continue + if ignoreLabel and insideLabel == ignoreLabel: + continue + intersectionPoint = self._polygonIntersect.getIntersection(ray, polygon) if not intersectionPoint: continue @@ -107,7 +112,9 @@ def _composeIntersection(ray: Ray, intersection: Intersection) -> Optional[Inter class SimpleIntersectionFinder(IntersectionFinder): - def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]: + def findIntersection( + self, ray: Ray, currentSolidLabel: str, ignoreLabel: Optional[str] = None + ) -> Optional[Intersection]: """ Find the closest intersection between a ray and the scene. @@ -119,7 +126,7 @@ def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Interse If the solid bbox distance is greater than the closest intersection distance, then we can stop testing since they are ordered by distance. Note that bbox distance is zero if the ray starts inside. """ - bboxIntersections = self._findBBoxIntersectingSolids(ray, currentSolidLabel) + bboxIntersections = self._findBBoxIntersectingSolids(ray, currentSolidLabel, ignoreLabel) bboxIntersections.sort(key=lambda x: x[0]) closestDistance = sys.maxsize @@ -134,12 +141,16 @@ def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Interse return self._composeIntersection(ray, closestIntersection) - def _findBBoxIntersectingSolids(self, ray: Ray, currentSolidLabel: str) -> Optional[List[Tuple[float, Solid]]]: + def _findBBoxIntersectingSolids( + self, ray: Ray, currentSolidLabel: str, ignoreLabel: Optional[str] = None + ) -> Optional[List[Tuple[float, Solid]]]: """We need to handle the special case where ray starts inside bbox. The Box Intersect will not compute the intersection for this case and will instead return ray.origin. When that happens, distance will be 0, and we continue to check for possibly other contained solids.""" solidCandidates = [] for solid in self._scene.solids: + if ignoreLabel and solid.getLabel() == ignoreLabel: + continue if solid.getLabel() == currentSolidLabel: bboxIntersectionPoint = ray.origin else: @@ -158,16 +169,20 @@ def __init__(self, scene: Scene, constructor=NoSplitThreeAxesConstructor(), maxD self._scene.getBoundingBox(), self._scene.getPolygons(), constructor, maxDepth, minLeafSize ) - def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]: + def findIntersection( + self, ray: Ray, currentSolidLabel: str, ignoreLabel: Optional[str] = None + ) -> Optional[Intersection]: self._currentSolidLabel = currentSolidLabel - intersection = self._findIntersection(ray, self._partition.root) + intersection = self._findIntersection(ray, self._partition.root, ignoreLabel=ignoreLabel) return self._composeIntersection(ray, intersection) - def _findIntersection(self, ray: Ray, node: Node, closestDistance=sys.maxsize) -> Optional[Intersection]: + def _findIntersection( + self, ray: Ray, node: Node, closestDistance=sys.maxsize, ignoreLabel: Optional[str] = None + ) -> Optional[Intersection]: # todo: test if this node search is still compatible with the new core intersection logic # which can now require testing a polygon slightly behind the ray origin. if node.isLeaf: - intersection = self._findClosestPolygonIntersection(ray, node.polygons, self._currentSolidLabel) + intersection = self._findClosestPolygonIntersection(ray, node.polygons, self._currentSolidLabel, ignoreLabel) return intersection if not self._nodeIsWorthExploring(ray, node, closestDistance): @@ -175,7 +190,7 @@ def _findIntersection(self, ray: Ray, node: Node, closestDistance=sys.maxsize) - closestIntersection = None for child in node.children: - intersection = self._findIntersection(ray, child, closestDistance) + intersection = self._findIntersection(ray, child, closestDistance, ignoreLabel) if intersection is None: continue if intersection.distance < closestDistance: diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py index f2cde0b7..23ac28f8 100644 --- a/pytissueoptics/scene/solids/solid.py +++ b/pytissueoptics/scene/solids/solid.py @@ -81,7 +81,7 @@ def primitive(self) -> str: def bbox(self) -> BoundingBox: return self._bbox - def asDetector(self, halfAngle: float) -> "Solid": + def asDetector(self, halfAngle: float = np.pi / 2) -> "Solid": """Treat this solid as a detector with a given half angle in radians. Detectors will fully absorb a photon when incident within the half angle, else the photon will go through From 8ab7dbd29a7c8ab9945fe3034a5b8d1862c09a06 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Fri, 10 Oct 2025 18:15:38 -0400 Subject: [PATCH 06/14] (opencl) Implement detectors in propagation logic. --- .../energyLogging/energyLogger.py | 4 +- .../rayscattering/opencl/CLScene.py | 23 ++++++---- .../rayscattering/opencl/buffers/photonCL.py | 4 ++ .../rayscattering/opencl/buffers/surfaceCL.py | 6 +++ .../rayscattering/opencl/src/intersection.c | 13 ++++-- .../rayscattering/opencl/src/propagation.c | 46 +++++++++++++++++-- .../rayscattering/opencl/utils/CLKeyLog.py | 4 +- .../rayscattering/statistics/statistics.py | 8 ++-- .../tests/energyLogging/testEnergyLogger.py | 4 +- .../tests/opencl/src/testCLPhoton.py | 18 ++++---- .../tests/opencl/testCLKeyLog.py | 12 ++--- pytissueoptics/scene/scene/scene.py | 2 + 12 files changed, 103 insertions(+), 41 deletions(-) diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py index 98fda268..c24859b1 100644 --- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py +++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py @@ -13,7 +13,7 @@ from pytissueoptics.scene.logger.listArrayContainer import ListArrayContainer from pytissueoptics.scene.logger.logger import DataType, InteractionData, InteractionKey, Logger -from ..opencl.CLScene import NO_SOLID_LABEL +from ..opencl.CLScene import WORLD_SOLID_LABEL from .energyType import EnergyType @@ -344,7 +344,7 @@ def export(self, exportPath: str): filepath = f"{exportPath}.csv" with open(filepath, "w") as file: file.write("energy,x,y,z,solid_index,surface_index\n") - self._writeKeyData(file, InteractionKey(NO_SOLID_LABEL), -1, -1) + self._writeKeyData(file, InteractionKey(WORLD_SOLID_LABEL), -1, -1) for i, solidLabel in enumerate(solidLabels): self._writeKeyData(file, InteractionKey(solidLabel), i, -1) for j, surfaceLabel in enumerate(self._scene.getSurfaceLabels(solidLabel)): diff --git a/pytissueoptics/rayscattering/opencl/CLScene.py b/pytissueoptics/rayscattering/opencl/CLScene.py index 22793a3d..01de2db2 100644 --- a/pytissueoptics/rayscattering/opencl/CLScene.py +++ b/pytissueoptics/rayscattering/opencl/CLScene.py @@ -12,10 +12,10 @@ from pytissueoptics.rayscattering.scatteringScene import ScatteringScene NO_LOG_ID = 0 -NO_SOLID_ID = -1 +WORLD_SOLID_ID = -1 NO_SURFACE_ID = -1 FIRST_SOLID_ID = 1 -NO_SOLID_LABEL = "world" +WORLD_SOLID_LABEL = "world" class CLScene: @@ -40,32 +40,35 @@ def __init__(self, scene: ScatteringScene, nWorkUnits: int): self.vertices = VertexCL(self._vertices) def getMaterialID(self, material): + if material is None: + # Detector case. Set dummy value (not used). + return 0 return self._sceneMaterials.index(material) def getSolidID(self, solid): if solid is None: - return NO_SOLID_ID + return WORLD_SOLID_ID return self._solidLabels.index(solid.getLabel()) + FIRST_SOLID_ID def getSolidLabel(self, solidID): - if solidID == NO_SOLID_ID: - return NO_SOLID_LABEL + if solidID == WORLD_SOLID_ID: + return WORLD_SOLID_LABEL return self._solidLabels[solidID - FIRST_SOLID_ID] def getSolidIDs(self) -> List[int]: solidIDs = list(self._surfaceLabels.keys()) - solidIDs.insert(0, NO_SOLID_ID) + solidIDs.insert(0, WORLD_SOLID_ID) return solidIDs def getSurfaceIDs(self, solidID): - if solidID == NO_SOLID_ID: + if solidID == WORLD_SOLID_ID: return [NO_SURFACE_ID] surfaceIDs = list(self._surfaceLabels[solidID].keys()) surfaceIDs.insert(0, NO_SURFACE_ID) return surfaceIDs def getSurfaceLabel(self, solidID, surfaceID): - if solidID == NO_SOLID_ID: + if solidID == WORLD_SOLID_ID: return None if surfaceID == NO_SURFACE_ID: return None @@ -100,6 +103,8 @@ def _compileSurface(self, polygonRef, firstPolygonID, lastPolygonID): insideSolidID = self.getSolidID(insideEnvironment.solid) outsideSolidID = self.getSolidID(outsideEnvironment.solid) toSmooth = polygonRef.toSmooth + isDetector = insideEnvironment.solid.isDetector if insideEnvironment.solid else False + detectorCosine = insideEnvironment.solid.detectorAcceptanceCosine if isDetector else 0.0 self._surfacesInfo.append( SurfaceCLInfo( @@ -110,6 +115,8 @@ def _compileSurface(self, polygonRef, firstPolygonID, lastPolygonID): insideSolidID, outsideSolidID, toSmooth, + isDetector, + detectorCosine, ) ) diff --git a/pytissueoptics/rayscattering/opencl/buffers/photonCL.py b/pytissueoptics/rayscattering/opencl/buffers/photonCL.py index dc6775c7..8757f465 100644 --- a/pytissueoptics/rayscattering/opencl/buffers/photonCL.py +++ b/pytissueoptics/rayscattering/opencl/buffers/photonCL.py @@ -3,6 +3,8 @@ from .CLObject import CLObject, cl +NULL_SOLID_ID = 0 + class PhotonCL(CLObject): STRUCT_NAME = "Photon" @@ -14,6 +16,7 @@ class PhotonCL(CLObject): ("weight", cl.cltypes.float), ("materialID", cl.cltypes.uint), ("solidID", cl.cltypes.int), + ("lastIntersectedDetectorID", cl.cltypes.int), ] ) @@ -36,4 +39,5 @@ def _getInitialHostBuffer(self) -> np.ndarray: buffer["weight"] = self._weight buffer["materialID"] = self._materialID buffer["solidID"] = self._solidID + buffer["lastIntersectedDetectorID"] = NULL_SOLID_ID return buffer diff --git a/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py b/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py index 265aef05..c1346fc1 100644 --- a/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py +++ b/pytissueoptics/rayscattering/opencl/buffers/surfaceCL.py @@ -14,6 +14,8 @@ ("insideSolidID", int), ("outsideSolidID", int), ("toSmooth", bool), + ("isDetector", bool), + ("detectorCosine", float), ], ) @@ -29,6 +31,8 @@ class SurfaceCL(CLObject): ("insideSolidID", cl.cltypes.int), ("outsideSolidID", cl.cltypes.int), ("toSmooth", cl.cltypes.uint), + ("isDetector", cl.cltypes.uint), + ("detectorCosine", cl.cltypes.float), ] ) @@ -47,4 +51,6 @@ def _getInitialHostBuffer(self) -> np.ndarray: buffer[i]["insideSolidID"] = np.int32(surfaceInfo.insideSolidID) buffer[i]["outsideSolidID"] = np.int32(surfaceInfo.outsideSolidID) buffer[i]["toSmooth"] = np.uint32(surfaceInfo.toSmooth) + buffer[i]["isDetector"] = np.uint32(surfaceInfo.isDetector) + buffer[i]["detectorCosine"] = np.float32(surfaceInfo.detectorCosine) return buffer diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c index 2d62b5fd..fecfa67e 100644 --- a/pytissueoptics/rayscattering/opencl/src/intersection.c +++ b/pytissueoptics/rayscattering/opencl/src/intersection.c @@ -115,13 +115,18 @@ GemsBoxIntersection _getBBoxIntersection(Ray ray, float3 minCornerVector, float3 return intersection; } -void _findBBoxIntersectingSolids(Ray ray, Scene *scene, uint gid, uint photonSolidID) { +void _findBBoxIntersectingSolids(Ray ray, Scene *scene, uint gid, uint photonSolidID, uint ignoreSolidID) { for (uint i = 0; i < scene->nSolids; i++) { uint boxGID = gid * scene->nSolids + i; uint solidID = i + 1; scene->solidCandidates[boxGID].solidID = solidID; + if (solidID == ignoreSolidID) { + scene->solidCandidates[boxGID].distance = -1; + continue; + } + if (solidID == photonSolidID) { scene->solidCandidates[boxGID].distance = 0; continue; @@ -375,12 +380,12 @@ void _composeIntersection(Intersection *intersection, Ray *ray, Scene *scene) { intersection->distanceLeft = ray->length - intersection->distance; } -Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidID) { +Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidID, uint ignoreSolidID) { /* OpenCL implementation of the Python module SimpleIntersectionFinder See the Python module documentation for more details. */ - _findBBoxIntersectingSolids(ray, scene, gid, photonSolidID); + _findBBoxIntersectingSolids(ray, scene, gid, photonSolidID, ignoreSolidID); _sortSolidCandidates(scene, gid); Intersection closestIntersection; @@ -420,7 +425,7 @@ __kernel void findIntersections(__global Ray *rays, uint nSolids, __global Solid __global Triangle *triangles, __global Vertex *vertices, __global SolidCandidate *solidCandidates, __global Intersection *intersections) { uint gid = get_global_id(0); Scene scene = {nSolids, solids, surfaces, triangles, vertices, solidCandidates}; - intersections[gid] = findIntersection(rays[gid], &scene, gid, -1); + intersections[gid] = findIntersection(rays[gid], &scene, gid, -1, 0); } diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c index 1a471e0e..17120218 100644 --- a/pytissueoptics/rayscattering/opencl/src/propagation.c +++ b/pytissueoptics/rayscattering/opencl/src/propagation.c @@ -4,7 +4,8 @@ #include "intersection.c" #include "fresnel.c" -__constant int NO_SOLID_ID = -1; +__constant int NULL_SOLID_ID = 0; +__constant int WORLD_SOLID_ID = -1; __constant int NO_SURFACE_ID = -1; __constant float MIN_ANGLE = 0.0001f; @@ -86,7 +87,7 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo (*logIndex)++; int outsideSolidID = surfaces[intersection->surfaceID].outsideSolidID; - if (outsideSolidID == NO_SOLID_ID){ + if (outsideSolidID == WORLD_SOLID_ID){ return; } logID++; @@ -99,6 +100,29 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo (*logIndex)++; } +bool detectOrIgnore(Intersection *intersection, __global Photon *photons, __global Surface *surfaces, + __global DataPoint *logger, uint *logIndex, uint gid, uint photonID){ + // If the incidence angle is within the numerical aperture, absorb photon. + float cosIncidence = -1 * dot(intersection->normal, photons[photonID].direction); + float cosDetector = surfaces[intersection->surfaceID].detectorCosine; + + if (cosIncidence < cosDetector){ + return false; // Outside NA, ignore. + } + + logger[*logIndex].x = photons[photonID].position.x; + logger[*logIndex].y = photons[photonID].position.y; + logger[*logIndex].z = photons[photonID].position.z; + logger[*logIndex].solidID = surfaces[intersection->surfaceID].insideSolidID; + logger[*logIndex].delta_weight = photons[photonID].weight; + logger[*logIndex].surfaceID = NO_SURFACE_ID; + (*logIndex)++; + + // Absorb photon. + photons[photonID].weight = 0; + return true; +} + float reflectOrRefract(Intersection *intersection, __global Photon *photons, __constant Material *materials, __global Surface *surfaces, __global DataPoint *logger, uint *logIndex, __global uint *seeds, uint gid, uint photonID){ FresnelIntersection fresnelIntersection = computeFresnelIntersection(photons[photonID].direction, intersection, @@ -156,13 +180,27 @@ float propagateStep(float distance, __global Photon *photons, __constant Materia } Ray stepRay = {photons[photonID].position, photons[photonID].direction, distance}; - Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID); + Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID, photons[photonID].lastIntersectedDetectorID); + + photons[photonID].lastIntersectedDetectorID = NULL_SOLID_ID; // Reset ignored detector ID. float distanceLeft = 0; if (intersection.exists){ moveTo(intersection.position, photons, photonID); - distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID); + if (scene->surfaces[intersection.surfaceID].isDetector) { + if (detectOrIgnore(&intersection, photons, scene->surfaces, logger, logIndex, gid, photonID)) {; + return 0; // Skip unnecessary vertex check if detected. + } + + // Prevent re-intersecting with the same detector when passing through it. + photons[photonID].lastIntersectedDetectorID = scene->surfaces[intersection.surfaceID].insideSolidID; + + // Skipping vertex check for now. + return intersection.distanceLeft; + } else { + distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID); + } // Check if intersection lies too close to a vertex. int closeToVertexID = -1; diff --git a/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py b/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py index fcf22250..96ca3915 100644 --- a/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py +++ b/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py @@ -3,7 +3,7 @@ import numpy as np -from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SOLID_LABEL, CLScene +from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_LABEL, CLScene from pytissueoptics.scene.logger import InteractionKey, Logger SOLID_ID_COL = 4 @@ -42,7 +42,7 @@ def _extractNoKeyLog(self): noInteractionIndices = np.where(self._log[:, SOLID_ID_COL] == NO_LOG_ID)[0] self._log = self._log[:, :4] self._log = np.delete(self._log, noInteractionIndices, axis=0) - self._keyLog[InteractionKey(NO_SOLID_LABEL, None)] = self._log + self._keyLog[InteractionKey(WORLD_SOLID_LABEL, None)] = self._log def _sortLocal(self): """Sorts the log locally by solidID and surfaceID, diff --git a/pytissueoptics/rayscattering/statistics/statistics.py b/pytissueoptics/rayscattering/statistics/statistics.py index 354d743d..4e59a874 100644 --- a/pytissueoptics/rayscattering/statistics/statistics.py +++ b/pytissueoptics/rayscattering/statistics/statistics.py @@ -8,7 +8,7 @@ from pytissueoptics.rayscattering import utils from pytissueoptics.rayscattering.display.views.defaultViews import View2DProjection from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloud, PointCloudFactory -from pytissueoptics.rayscattering.opencl.CLScene import NO_SOLID_LABEL +from pytissueoptics.rayscattering.opencl.CLScene import WORLD_SOLID_LABEL @dataclass @@ -50,11 +50,11 @@ def report(self, solidLabel: str = None, saveToFile: str = None, verbose=True): def _computeStats(self, solidLabel: str = None): solidLabels = [solidLabel] - if solidLabel is None or utils.labelsEqual(solidLabel, NO_SOLID_LABEL): + if solidLabel is None or utils.labelsEqual(solidLabel, WORLD_SOLID_LABEL): solidLabels = self._logger.getSeenSolidLabels() for solidLabel in solidLabels: - if solidLabel == NO_SOLID_LABEL: + if solidLabel == WORLD_SOLID_LABEL: continue try: absorbance = self.getAbsorbance(solidLabel) @@ -70,7 +70,7 @@ def _computeStats(self, solidLabel: str = None): def _makeReport(self, solidLabel: str = None, reportString: str = ""): if solidLabel: - if solidLabel == NO_SOLID_LABEL: + if solidLabel == WORLD_SOLID_LABEL: reportString += self._reportWorld(solidLabel) else: reportString += self._reportSolid(solidLabel) diff --git a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py index 6b599cde..4c41c438 100644 --- a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py +++ b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py @@ -18,7 +18,7 @@ ) from pytissueoptics.rayscattering.energyLogging import EnergyLogger from pytissueoptics.rayscattering.materials import ScatteringMaterial -from pytissueoptics.rayscattering.opencl.CLScene import NO_SOLID_LABEL +from pytissueoptics.rayscattering.opencl.CLScene import WORLD_SOLID_LABEL from pytissueoptics.rayscattering.samples import PhantomTissue from pytissueoptics.rayscattering.scatteringScene import ScatteringScene from pytissueoptics.scene.geometry import Vector @@ -299,7 +299,7 @@ def testWhenExport_shouldExport3DDataPointsToFile(self): self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("middleLayer")) self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("frontLayer", "interface1")) self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere")) - self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(NO_SOLID_LABEL)) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL)) with tempfile.TemporaryDirectory() as tempDir: filePath = os.path.join(tempDir, "test_sim") diff --git a/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py b/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py index 64c65846..3784b2d1 100644 --- a/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py +++ b/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py @@ -21,7 +21,7 @@ VertexCL, ) from pytissueoptics.rayscattering.opencl.CLProgram import CLProgram -from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SOLID_ID, NO_SURFACE_ID, CLScene +from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_ID, NO_SURFACE_ID, CLScene from pytissueoptics.rayscattering.opencl.config.CLConfig import OPENCL_SOURCE_DIR from pytissueoptics.scene.geometry import Vertex @@ -172,7 +172,7 @@ def testWhenLogIntersectionWithPhotonLeavingIntoWorld_shouldLogOnlyOneIntersecti self.INITIAL_DIRECTION.normalize() intersectionNormal = Vector(0, 0, 1) insideSolidID = self.INITIAL_SOLID_ID - outsideSolidID = NO_SOLID_ID + outsideSolidID = WORLD_SOLID_ID surfaceID = 0 surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)]) @@ -196,10 +196,10 @@ def testWhenLogIntersectionWithPhotonEnteringFromWorld_shouldLogOnlyOneIntersect ): self.INITIAL_DIRECTION = Vector(0, 1, -1) self.INITIAL_DIRECTION.normalize() - self.INITIAL_SOLID_ID = NO_SOLID_ID + self.INITIAL_SOLID_ID = WORLD_SOLID_ID intersectionNormal = Vector(0, 0, 1) insideSolidID = 9 - outsideSolidID = NO_SOLID_ID + outsideSolidID = WORLD_SOLID_ID surfaceID = 0 surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)]) @@ -255,13 +255,13 @@ def testWhenInteract_shouldLogInteraction(self): def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self): # => Photon is trying to enter solid intersectionNormal = Vector(0, 1, 0) - self.INITIAL_SOLID_ID = NO_SOLID_ID + self.INITIAL_SOLID_ID = WORLD_SOLID_ID self.INITIAL_DIRECTION = Vector(1, -1, 0) self.INITIAL_DIRECTION.normalize() self._mockFresnelIntersection(isReflected=True, incidencePlane=Vector(0, 0, 1), angleDeflection=np.pi / 2) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=NO_SOLID_ID, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=WORLD_SOLID_ID, toSmooth=False)]) photonResult = self._photonFunc( "reflectOrRefract", intersectionNormal, @@ -279,7 +279,7 @@ def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self) self._assertVectorAlmostEqual(expectedDirection, photonResult.direction, places=6) expectedPosition = self.INITIAL_POSITION self._assertVectorAlmostEqual(expectedPosition, photonResult.position) - self.assertEqual(NO_SOLID_ID, photonResult.solidID) + self.assertEqual(WORLD_SOLID_ID, photonResult.solidID) dataPoint = self._getDataPointResult(logger) self.assertEqual(NO_LOG_ID, dataPoint.solidID) @@ -287,7 +287,7 @@ def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self) def testWhenReflectOrRefractWithRefractingIntersection_shouldRefractPhoton(self): # => Photon enters solid intersectionNormal = Vector(0, 1, 0) - self.INITIAL_SOLID_ID = NO_SOLID_ID + self.INITIAL_SOLID_ID = WORLD_SOLID_ID self.INITIAL_DIRECTION = Vector(1, -1, 0) self.INITIAL_DIRECTION.normalize() insideSolidID = 9 @@ -301,7 +301,7 @@ def testWhenReflectOrRefractWithRefractingIntersection_shouldRefractPhoton(self) ) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID=NO_SOLID_ID, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID=WORLD_SOLID_ID, toSmooth=False)]) photonResult = self._photonFunc( "reflectOrRefract", intersectionNormal, diff --git a/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py b/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py index 9e27708c..ca75b726 100644 --- a/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py +++ b/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py @@ -5,7 +5,7 @@ from pytissueoptics import Cube, EnergyLogger, ScatteringMaterial, ScatteringScene, Sphere from pytissueoptics.rayscattering.opencl import OPENCL_OK -from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SOLID_ID, NO_SOLID_LABEL, NO_SURFACE_ID, CLScene +from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_ID, WORLD_SOLID_LABEL, NO_SURFACE_ID, CLScene from pytissueoptics.rayscattering.opencl.utils import CLKeyLog from pytissueoptics.scene.logger import InteractionKey @@ -29,7 +29,7 @@ def _createTestLog(self, sceneCL): [ [0, 0, 0, 0, cubeID, NO_SURFACE_ID], [2, 0, 0, 0, sphereID, NO_SURFACE_ID], - [3, 0, 0, 0, NO_SOLID_ID, NO_SURFACE_ID], + [3, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], [4, 9, 9, 9, NO_LOG_ID, 9], [5, 0, 0, 0, cubeID, cubeSurfaceIDs[1]], [6, 0, 0, 0, sphereID, sphereSurfaceIDs[1]], @@ -54,7 +54,7 @@ def testGivenCLKeyLog_whenTransferToSceneLogger_shouldLogDataWithInteractionKeys expectedValuesWithKeys = [ (2, InteractionKey(self.sphere.getLabel())), - (3, InteractionKey(NO_SOLID_LABEL)), + (3, InteractionKey(WORLD_SOLID_LABEL)), (5, InteractionKey(self.cube.getLabel(), self.cube.surfaceLabels[0])), (6, InteractionKey(self.sphere.getLabel(), self.sphere.surfaceLabels[0])), ] @@ -67,9 +67,9 @@ def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWit sceneCL = CLScene(self.scene, nWorkUnits=10) log = np.array( [ - [1, 0, 0, 0, NO_SOLID_ID, NO_SURFACE_ID], + [1, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], [2, 0, 0, 0, NO_LOG_ID, 99], - [3, 0, 0, 0, NO_SOLID_ID, NO_SURFACE_ID], + [3, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], ] ) clKeyLog = CLKeyLog(log, sceneCL) @@ -80,4 +80,4 @@ def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWit verify(sceneLogger, times=1).logDataPointArray(...) expectedWorldData = arg_that(lambda arg: np.array_equal(arg, np.array([[1, 0, 0, 0], [3, 0, 0, 0]]))) - verify(sceneLogger).logDataPointArray(expectedWorldData, InteractionKey(NO_SOLID_LABEL)) + verify(sceneLogger).logDataPointArray(expectedWorldData, InteractionKey(WORLD_SOLID_LABEL)) diff --git a/pytissueoptics/scene/scene/scene.py b/pytissueoptics/scene/scene/scene.py index 7e7b4085..25ad9d1e 100644 --- a/pytissueoptics/scene/scene/scene.py +++ b/pytissueoptics/scene/scene/scene.py @@ -142,6 +142,8 @@ def getPolygons(self) -> List[Polygon]: def getMaterials(self) -> list: materials = [self._worldMaterial] for solid in self._solids: + if solid.isDetector: + continue surfaceLabels = solid.surfaceLabels for surfaceLabel in surfaceLabels: material = solid.getPolygons(surfaceLabel)[0].insideEnvironment.material From cd0554bf99e3b48c200dafca44459979724fa315 Mon Sep 17 00:00:00 2001 From: JLBegin Date: Sat, 11 Oct 2025 11:47:23 -0400 Subject: [PATCH 07/14] fix cuboid hash --- pytissueoptics/scene/solids/cuboid.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/pytissueoptics/scene/solids/cuboid.py b/pytissueoptics/scene/solids/cuboid.py index 609efd84..29dfa204 100644 --- a/pytissueoptics/scene/solids/cuboid.py +++ b/pytissueoptics/scene/solids/cuboid.py @@ -154,3 +154,15 @@ def _getLayerPolygons(self, layerLabel: str) -> List[Polygon]: def _geometryParams(self) -> dict: return {"shape": self.shape} + + def __hash__(self): + baseHash = super().__hash__() + if not self.isStack(): + return baseHash + + materials = set() + for surfaceLabel in self.surfaceLabels: + material = self.getPolygons(surfaceLabel)[0].insideEnvironment.material + materials.add(material) + materialsHash = hash(tuple(materials)) + return hash((baseHash, materialsHash)) From 3f26ff5c036f674f44153466f5f8a728397e4ec5 Mon Sep 17 00:00:00 2001 From: JLBegin Date: Sat, 11 Oct 2025 11:48:14 -0400 Subject: [PATCH 08/14] fix photons per batch when workUnits > N --- pytissueoptics/rayscattering/opencl/utils/CLParameters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytissueoptics/rayscattering/opencl/utils/CLParameters.py b/pytissueoptics/rayscattering/opencl/utils/CLParameters.py index 00883007..21ba3a43 100644 --- a/pytissueoptics/rayscattering/opencl/utils/CLParameters.py +++ b/pytissueoptics/rayscattering/opencl/utils/CLParameters.py @@ -12,8 +12,8 @@ def __init__(self, N, AVG_IT_PER_PHOTON): nBatch = 1 / CONFIG.BATCH_LOAD_FACTOR avgPhotonsPerBatch = int(np.ceil(N / min(nBatch, CONFIG.N_WORK_UNITS))) self._maxLoggerMemory = self._calculateAverageBatchMemorySize(avgPhotonsPerBatch, AVG_IT_PER_PHOTON) - self._maxPhotonsPerBatch = min(2 * avgPhotonsPerBatch, N) self._workItemAmount = CONFIG.N_WORK_UNITS + self.maxPhotonsPerBatch = min(2 * avgPhotonsPerBatch, N) self._assertEnoughRAM() From ed1ce904de069ab80abaed243a838f3d54979954 Mon Sep 17 00:00:00 2001 From: JLBegin Date: Sat, 11 Oct 2025 15:30:29 -0400 Subject: [PATCH 09/14] track photonIDs to support tracing --- .../rayscattering/opencl/CLPhotons.py | 1 + .../rayscattering/opencl/buffers/dataPointCL.py | 1 + .../rayscattering/opencl/buffers/photonCL.py | 7 ++++++- .../rayscattering/opencl/src/propagation.c | 4 ++++ .../rayscattering/opencl/utils/CLKeyLog.py | 17 +++++++++-------- pytissueoptics/rayscattering/photon.py | 11 ++++++----- pytissueoptics/rayscattering/source.py | 2 +- pytissueoptics/scene/logger/logger.py | 16 ++++++++++------ 8 files changed, 38 insertions(+), 21 deletions(-) diff --git a/pytissueoptics/rayscattering/opencl/CLPhotons.py b/pytissueoptics/rayscattering/opencl/CLPhotons.py index 7e793901..43ed16d4 100644 --- a/pytissueoptics/rayscattering/opencl/CLPhotons.py +++ b/pytissueoptics/rayscattering/opencl/CLPhotons.py @@ -54,6 +54,7 @@ def propagate(self, IPP: float, verbose: bool = False): self._directions[params.maxPhotonsPerBatch :], materialID=scene.getMaterialID(self._initialMaterial), solidID=scene.getSolidID(self._initialSolid), + startID=params.maxPhotonsPerBatch, ) photonPool.make(program.device) seeds = SeedCL(params.maxPhotonsPerBatch) diff --git a/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py b/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py index 5a3cfeb4..cc153956 100644 --- a/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py +++ b/pytissueoptics/rayscattering/opencl/buffers/dataPointCL.py @@ -12,6 +12,7 @@ class DataPointCL(CLObject): ("x", cl.cltypes.float), ("y", cl.cltypes.float), ("z", cl.cltypes.float), + ("photonID", cl.cltypes.uint), ("solidID", cl.cltypes.int), ("surfaceID", cl.cltypes.int), ] diff --git a/pytissueoptics/rayscattering/opencl/buffers/photonCL.py b/pytissueoptics/rayscattering/opencl/buffers/photonCL.py index 8757f465..9a74782a 100644 --- a/pytissueoptics/rayscattering/opencl/buffers/photonCL.py +++ b/pytissueoptics/rayscattering/opencl/buffers/photonCL.py @@ -17,16 +17,20 @@ class PhotonCL(CLObject): ("materialID", cl.cltypes.uint), ("solidID", cl.cltypes.int), ("lastIntersectedDetectorID", cl.cltypes.int), + ("ID", cl.cltypes.uint), ] ) - def __init__(self, positions: np.ndarray, directions: np.ndarray, materialID: int, solidID: int, weight=1.0): + def __init__( + self, positions: np.ndarray, directions: np.ndarray, materialID: int, solidID: int, weight=1.0, startID=0 + ): self._positions = positions self._directions = directions self._N = positions.shape[0] self._materialID = materialID self._solidID = solidID self._weight = weight + self._startID = startID super().__init__() @@ -40,4 +44,5 @@ def _getInitialHostBuffer(self) -> np.ndarray: buffer["materialID"] = self._materialID buffer["solidID"] = self._solidID buffer["lastIntersectedDetectorID"] = NULL_SOLID_ID + buffer["ID"] = np.arange(self._startID, self._startID + self._N, dtype=np.uint32) return buffer diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c index 17120218..d4e296bb 100644 --- a/pytissueoptics/rayscattering/opencl/src/propagation.c +++ b/pytissueoptics/rayscattering/opencl/src/propagation.c @@ -37,6 +37,7 @@ void interact(__global Photon *photons, __constant Material *materials, __global logger[logIndex].delta_weight = delta_weight; logger[logIndex].solidID = photons[photonID].solidID; logger[logIndex].surfaceID = NO_SURFACE_ID; + logger[logIndex].photonID = photons[photonID].ID; } void scatter(__global Photon *photons, __constant Material *materials, __global uint *seeds, __global DataPoint *logger, @@ -80,6 +81,7 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo logger[logID].z = photons[photonID].position.z; logger[logID].surfaceID = intersection->surfaceID; logger[logID].solidID = surfaces[intersection->surfaceID].insideSolidID; + logger[logID].photonID = photons[photonID].ID; bool isLeavingSurface = dot(photons[photonID].direction, intersection->normal) > 0; int sign = isLeavingSurface ? 1 : -1; @@ -97,6 +99,7 @@ void logIntersection(Intersection *intersection, __global Photon *photons, __glo logger[logID].surfaceID = intersection->surfaceID; logger[logID].solidID = outsideSolidID; logger[logID].delta_weight = -sign * photons[photonID].weight; + logger[logID].photonID = photons[photonID].ID; (*logIndex)++; } @@ -116,6 +119,7 @@ bool detectOrIgnore(Intersection *intersection, __global Photon *photons, __glob logger[*logIndex].solidID = surfaces[intersection->surfaceID].insideSolidID; logger[*logIndex].delta_weight = photons[photonID].weight; logger[*logIndex].surfaceID = NO_SURFACE_ID; + logger[*logIndex].photonID = photons[photonID].ID; (*logIndex)++; // Absorb photon. diff --git a/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py b/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py index 96ca3915..7a2aeb57 100644 --- a/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py +++ b/pytissueoptics/rayscattering/opencl/utils/CLKeyLog.py @@ -6,15 +6,16 @@ from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_LABEL, CLScene from pytissueoptics.scene.logger import InteractionKey, Logger -SOLID_ID_COL = 4 -SURFACE_ID_COL = 5 +SOLID_ID_COL = 5 +SURFACE_ID_COL = 6 class CLKeyLog: - """Parses a DataPointCL array of shape (N, 6) where each point is of the form - (weight, x, y, z, solidID, surfaceID) to extract a dictionary of InteractionKey - and their corresponding datapoint array of the form (weight, x, y, z). The - translation from IDs to their corresponding labels is done using the given CLScene.""" + """Parses a DataPointCL array of shape (N, 7) where each point is of the form + (weight, x, y, z, photonID, solidID, surfaceID) to extract a dictionary of InteractionKey + and their corresponding datapoint array of the form (weight, x, y, z, photonID). The + translation from IDs to their corresponding labels is done using the given CLScene. + """ def __init__(self, log: np.ndarray, sceneCL: CLScene): self._log = log @@ -40,7 +41,7 @@ def _extractKeyLog(self): def _extractNoKeyLog(self): noInteractionIndices = np.where(self._log[:, SOLID_ID_COL] == NO_LOG_ID)[0] - self._log = self._log[:, :4] + self._log = self._log[:, :5] self._log = np.delete(self._log, noInteractionIndices, axis=0) self._keyLog[InteractionKey(WORLD_SOLID_LABEL, None)] = self._log @@ -95,7 +96,7 @@ def _sort(log: np.ndarray, columns: List[int]): def _merge(self): """Merges the local batches into a single key log with unique interaction keys.""" - self._log = self._log[:, :4] + self._log = self._log[:, :5] for i, batchKeyIndices in enumerate(self._keyIndices): batchStartIndex = i * self._batchSize diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py index 9a8bcd61..a7456f42 100644 --- a/pytissueoptics/rayscattering/photon.py +++ b/pytissueoptics/rayscattering/photon.py @@ -16,11 +16,12 @@ class Photon: - def __init__(self, position: Vector, direction: Vector): + def __init__(self, position: Vector, direction: Vector, ID: int = 0): self._position = position self._direction = direction self._weight = 1 self._environment: Environment = None + self._ID = ID self._er = self._direction.getAnyOrthogonal() self._hasContext = False @@ -147,7 +148,7 @@ def detectOrIgnore(self, intersection: Intersection) -> bool: def _detectAndLog(self, solidLabel: str): if self._logger: key = InteractionKey(solidLabel) - self._logger.logDataPoint(self._weight, self._position, key) + self._logger.logDataPoint(self._weight, self._position, key, self._ID) self._weight = 0 def reflectOrRefract(self, intersection: Intersection): @@ -242,16 +243,16 @@ def _logIntersection(self, intersection: Intersection): key = InteractionKey(solidLabelA, intersection.surfaceLabel) isLeavingSurface = self._direction.dot(intersection.normal) > 0 sign = 1 if isLeavingSurface else -1 - self._logger.logDataPoint(sign * self._weight, self._position, key) + self._logger.logDataPoint(sign * self._weight, self._position, key, self._ID) solidB = intersection.outsideEnvironment.solid if solidB is None: return solidLabelB = solidB.getLabel() key = InteractionKey(solidLabelB, intersection.surfaceLabel) - self._logger.logDataPoint(-sign * self._weight, self._position, key) + self._logger.logDataPoint(-sign * self._weight, self._position, key, self._ID) def _logWeightDecrease(self, delta): if self._logger: key = InteractionKey(self.solidLabel) - self._logger.logDataPoint(delta, self._position, key) + self._logger.logDataPoint(delta, self._position, key, self._ID) diff --git a/pytissueoptics/rayscattering/source.py b/pytissueoptics/rayscattering/source.py index 59c0c1f1..a1b0f289 100644 --- a/pytissueoptics/rayscattering/source.py +++ b/pytissueoptics/rayscattering/source.py @@ -138,7 +138,7 @@ def _loadPhotons(self): def _loadPhotonsCPU(self): positions, directions = self.getInitialPositionsAndDirections() for i in range(self._N): - self._photons.append(Photon(Vector(*positions[i]), Vector(*directions[i]))) + self._photons.append(Photon(Vector(*positions[i]), Vector(*directions[i]), ID=i)) def _loadPhotonsOpenCL(self): positions, directions = self.getInitialPositionsAndDirections() diff --git a/pytissueoptics/scene/logger/logger.py b/pytissueoptics/scene/logger/logger.py index 2e813795..3b5c8d08 100644 --- a/pytissueoptics/scene/logger/logger.py +++ b/pytissueoptics/scene/logger/logger.py @@ -78,8 +78,11 @@ def getStoredSurfaceLabels(self, solidLabel: str) -> List[str]: def logPoint(self, point: Vector, key: InteractionKey = None): self._appendData([point.x, point.y, point.z], DataType.POINT, key) - def logDataPoint(self, value: float, position: Vector, key: InteractionKey): - self._appendData([value, position.x, position.y, position.z], DataType.DATA_POINT, key) + def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: Optional[int] = None): + dataPoint = [value, position.x, position.y, position.z] + if ID is not None: + dataPoint.append(ID) + self._appendData(dataPoint, DataType.DATA_POINT, key) def logSegment(self, start: Vector, end: Vector, key: InteractionKey = None): self._appendData([start.x, start.y, start.z, end.x, end.y, end.z], DataType.SEGMENT, key) @@ -90,8 +93,9 @@ def logPointArray(self, array: np.ndarray, key: InteractionKey = None): self._appendData(array, DataType.POINT, key) def logDataPointArray(self, array: np.ndarray, key: InteractionKey): - """'array' must be of shape (n, 4) where the second axis is (value, x, y, z)""" - assert array.shape[1] == 4 and array.ndim == 2, "Data point array must be of shape (n, 4)" + """'array' must be of shape (n, 4) or (n, 5) where the second axis is (value, x, y, z) + or (value, x, y, z, photonID). The photonID column is optional.""" + assert array.shape[1] in [4, 5] and array.ndim == 2, "Data point array must be of shape (n, 4) or (n, 5)" self._appendData(array, DataType.DATA_POINT, key) def logSegmentArray(self, array: np.ndarray, key: InteractionKey = None): @@ -127,8 +131,8 @@ def getPoints(self, key: InteractionKey = None) -> np.ndarray: return self._getData(DataType.POINT, key) def getRawDataPoints(self, key: InteractionKey = None) -> np.ndarray: - """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) where - the second axis is (value, x, y, z).""" + """All raw 3D data points recorded for this InteractionKey (not binned). Array of shape (n, 4) or (n, 5) where + the second axis is (value, x, y, z) or (value, x, y, z, photonID) if photon IDs were logged.""" return self._getData(DataType.DATA_POINT, key) def getSegments(self, key: InteractionKey = None) -> np.ndarray: From 126b34a679ae2d72d51a86647d381a8ba3d301b1 Mon Sep 17 00:00:00 2001 From: JLBegin Date: Sat, 11 Oct 2025 19:27:50 -0400 Subject: [PATCH 10/14] implement `detectedBy` filter for logger and displays --- pytissueoptics/examples/rayscattering/ex06.py | 46 ++++++++++ .../rayscattering/display/viewer.py | 13 +-- .../display/views/defaultViews.py | 42 ++++++++- .../rayscattering/display/views/view2D.py | 11 +++ .../energyLogging/energyLogger.py | 85 ++++++++++++++++--- pytissueoptics/scene/logger/logger.py | 2 +- 6 files changed, 178 insertions(+), 21 deletions(-) create mode 100644 pytissueoptics/examples/rayscattering/ex06.py diff --git a/pytissueoptics/examples/rayscattering/ex06.py b/pytissueoptics/examples/rayscattering/ex06.py new file mode 100644 index 00000000..c0f8bceb --- /dev/null +++ b/pytissueoptics/examples/rayscattering/ex06.py @@ -0,0 +1,46 @@ +import numpy as np + +import env # noqa: F401 + +from pytissueoptics import * # noqa: F403 + +TITLE = "Detectors" + +DESCRIPTION = """ Sampling volume simulation with a directional source and a circle detector. """ + + +def exampleCode(): + N = 1000000 if hardwareAccelerationIsAvailable() else 1000 + + material = ScatteringMaterial(mu_s=10, mu_a=1, g=0.98, n=1.0) + + cube = Cuboid( + a=3, b=3, c=3, position=Vector(0, 0, 0), material=material, label="cube" + ) + detector = Circle( + radius=0.25, orientation=Vector(0, 0, -1), position=Vector(0, 0.5, 1.501), label="detector", + ).asDetector(halfAngle=np.pi / 4) + scene = ScatteringScene([cube, detector]) + + logger = EnergyLogger(scene) + source = DirectionalSource( + position=Vector(0, -0.5, -2), direction=Vector(0, 0, 1), N=N, diameter=0.5 + ) + + source.propagate(scene, logger) + + viewer = Viewer(scene, source, logger) + viewer.reportStats() + + # Either filter the whole logger + logger.filter(detectedBy="detector") + viewer.show2D(View2DProjectionX()) + viewer.show3D() + + # ... or keep the logger intact and filter individual views + # viewer.show2D(View2DProjectionX(detectedBy="detector")) + # viewer.show3D(pointCloudStyle=PointCloudStyle(detectedBy="detector")) + + +if __name__ == "__main__": + exampleCode() diff --git a/pytissueoptics/rayscattering/display/viewer.py b/pytissueoptics/rayscattering/display/viewer.py index adf529cb..fc793f3e 100644 --- a/pytissueoptics/rayscattering/display/viewer.py +++ b/pytissueoptics/rayscattering/display/viewer.py @@ -66,6 +66,7 @@ def __init__( showSurfacePointsLeaving: bool = True, showSurfacePointsEntering: bool = False, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, showPointsAsSpheres: bool = False, pointSize: float = 0.15, scaleWithValue: bool = True, @@ -82,6 +83,7 @@ def __init__( self.showSurfacePointsLeaving = showSurfacePointsLeaving self.showSurfacePointsEntering = showSurfacePointsEntering self.energyType = energyType + self.detectedBy = detectedBy self.showPointsAsSpheres = showPointsAsSpheres self.pointSize = pointSize @@ -103,7 +105,6 @@ def __init__(self, scene: ScatteringScene, source: Source, logger: EnergyLogger) self._logger = logger self._viewer3D = None - self._pointCloudFactory = PointCloudFactory(logger) self._profileFactory = ProfileFactory(scene, logger) def listViews(self): @@ -169,7 +170,7 @@ def show3DVolumeSlicer( f"Consider using a larger binSize or tighter limits." ) - points = self._pointCloudFactory.getPointCloudOfSolids(energyType=energyType).solidPoints + points = PointCloudFactory(self._logger).getPointCloudOfSolids(energyType=energyType).solidPoints try: hist, _ = np.histogramdd(points[:, 1:], bins=bins, weights=points[:, 0], range=limits) except MemoryError: @@ -215,8 +216,10 @@ def reportStats(self, solidLabel: str = None, saveToFile: str = None, verbose=Tr stats.report(solidLabel=solidLabel, saveToFile=saveToFile, verbose=verbose) def _addPointCloud(self, style: PointCloudStyle): - pointCloud = self._pointCloudFactory.getPointCloud( - solidLabel=style.solidLabel, surfaceLabel=style.surfaceLabel, energyType=style.energyType + logger = self._logger if style.detectedBy is None else self._logger.getFiltered(style.detectedBy) + + pointCloud = PointCloudFactory(logger).getPointCloud( + solidLabel=style.solidLabel, surfaceLabel=style.surfaceLabel, energyType=style.energyType, ) self._drawPointCloudOfSolids(pointCloud, style) @@ -241,7 +244,7 @@ def _drawPointCloudOfSurfaces(self, pointCloud: PointCloud, style: PointCloudSty if pointCloud.surfacePoints is None: return - surfacePoints = np.empty((0, 4)) + surfacePoints = np.empty((0, 5)) if style.showSurfacePointsLeaving: surfacePoints = np.vstack((surfacePoints, pointCloud.leavingSurfacePoints)) diff --git a/pytissueoptics/rayscattering/display/views/defaultViews.py b/pytissueoptics/rayscattering/display/views/defaultViews.py index a18092eb..6a2e17ee 100644 --- a/pytissueoptics/rayscattering/display/views/defaultViews.py +++ b/pytissueoptics/rayscattering/display/views/defaultViews.py @@ -21,6 +21,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( projectionDirection, @@ -29,6 +30,7 @@ def __init__( limits=limits, binSize=binSize, energyType=energyType, + detectedBy=detectedBy, ) def _filter(self, dataPoints: np.ndarray) -> np.ndarray: @@ -42,9 +44,15 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( - *DEFAULT_X_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType + *DEFAULT_X_VIEW_DIRECTIONS, + solidLabel=solidLabel, + limits=limits, + binSize=binSize, + energyType=energyType, + detectedBy=detectedBy, ) @@ -55,9 +63,15 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( - *DEFAULT_Y_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType + *DEFAULT_Y_VIEW_DIRECTIONS, + solidLabel=solidLabel, + limits=limits, + binSize=binSize, + energyType=energyType, + detectedBy=detectedBy, ) @@ -68,9 +82,15 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( - *DEFAULT_Z_VIEW_DIRECTIONS, solidLabel=solidLabel, limits=limits, binSize=binSize, energyType=energyType + *DEFAULT_Z_VIEW_DIRECTIONS, + solidLabel=solidLabel, + limits=limits, + binSize=binSize, + energyType=energyType, + detectedBy=detectedBy, ) @@ -84,6 +104,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + detectedBy: str = None, ): super().__init__( projectionDirection, @@ -93,6 +114,7 @@ def __init__( surfaceEnergyLeaving=surfaceEnergyLeaving, limits=limits, binSize=binSize, + detectedBy=detectedBy, ) def _filter(self, dataPoints: np.ndarray) -> np.ndarray: @@ -116,6 +138,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + detectedBy: str = None, ): super().__init__( *DEFAULT_X_VIEW_DIRECTIONS, @@ -124,6 +147,7 @@ def __init__( surfaceEnergyLeaving=surfaceEnergyLeaving, limits=limits, binSize=binSize, + detectedBy=detectedBy, ) @@ -135,6 +159,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + detectedBy: str = None, ): super().__init__( *DEFAULT_Y_VIEW_DIRECTIONS, @@ -143,6 +168,7 @@ def __init__( surfaceEnergyLeaving=surfaceEnergyLeaving, limits=limits, binSize=binSize, + detectedBy=detectedBy, ) @@ -154,6 +180,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, + detectedBy: str = None, ): super().__init__( *DEFAULT_Z_VIEW_DIRECTIONS, @@ -162,6 +189,7 @@ def __init__( surfaceEnergyLeaving=surfaceEnergyLeaving, limits=limits, binSize=binSize, + detectedBy=detectedBy, ) @@ -176,6 +204,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( projectionDirection, @@ -186,6 +215,7 @@ def __init__( limits=limits, binSize=binSize, energyType=energyType, + detectedBy=detectedBy, ) def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float, float, float]): @@ -210,6 +240,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( *DEFAULT_X_VIEW_DIRECTIONS, @@ -219,6 +250,7 @@ def __init__( limits=limits, binSize=binSize, energyType=energyType, + detectedBy=detectedBy, ) @@ -231,6 +263,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( *DEFAULT_Y_VIEW_DIRECTIONS, @@ -240,6 +273,7 @@ def __init__( limits=limits, binSize=binSize, energyType=energyType, + detectedBy=detectedBy, ) @@ -252,6 +286,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): super().__init__( *DEFAULT_Z_VIEW_DIRECTIONS, @@ -261,4 +296,5 @@ def __init__( limits=limits, binSize=binSize, energyType=energyType, + detectedBy=detectedBy, ) diff --git a/pytissueoptics/rayscattering/display/views/view2D.py b/pytissueoptics/rayscattering/display/views/view2D.py index e61e6de2..88f86769 100644 --- a/pytissueoptics/rayscattering/display/views/view2D.py +++ b/pytissueoptics/rayscattering/display/views/view2D.py @@ -50,6 +50,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, + detectedBy: str = None, ): """ The 2D view plane is obtained by looking towards the 'projectionDirection'. The 'horizontalDirection' @@ -78,6 +79,7 @@ def __init__( self._position = position self._thickness = thickness self._energyType = energyType + self._detectedBy = detectedBy limits = [sorted(limit) for limit in limits] if limits else [None, None] self._limitsU, self._limitsV = limits @@ -103,6 +105,7 @@ def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float limits = [self._limitsU, self._limitsV] self._binsU, self._binsV = [int((limit[1] - limit[0]) / bin_) for limit, bin_ in zip(limits, self._binSize)] + self._binsU, self._binsV = max(1, self._binsU), max(1, self._binsV) if self._verticalIsNegative: self._limitsV = self._limitsV[::-1] @@ -121,6 +124,10 @@ def setContext(self, limits3D: List[Tuple[float, float]], binSize3D: Tuple[float def energyType(self) -> EnergyType: return self._energyType + @property + def detectedBy(self) -> Optional[str]: + return self._detectedBy + def extractData(self, dataPoints: np.ndarray): """ Used internally by Logger2D to store 3D datapoints into this 2D view. @@ -245,6 +252,8 @@ def isEqualTo(self, other: "View2D") -> bool: return False if self._energyType != other._energyType: return False + if self._detectedBy != other._detectedBy: + return False return True def isContainedBy(self, other: "View2D") -> bool: @@ -260,6 +269,8 @@ def isContainedBy(self, other: "View2D") -> bool: return False if self._thickness != other._thickness: return False + if self._detectedBy != other._detectedBy: + return False # TODO: change/remove the following once the algorithm can extract a view contained inside a bigger view. isTransposed = self.axisU != other.axisU diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py index c24859b1..d21bd765 100644 --- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py +++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py @@ -1,7 +1,9 @@ +from __future__ import annotations + import json import os import pickle -from typing import List, Optional, TextIO, Union +from typing import List, Optional, TextIO, Union, Dict import numpy as np @@ -76,7 +78,7 @@ def addView(self, view: View2D) -> bool: return True if self.has3D: - self._compileViews([view]) + self._compileViews([view], detectedBy=view.detectedBy) self._views.append(view) return True @@ -233,11 +235,24 @@ def logDataPointArray(self, array: np.ndarray, key: InteractionKey): self._compileViews(self._views) self._delete3DData() - def logDataPoint(self, value: float, position: Vector, key: InteractionKey): - self.logDataPointArray(np.array([[value, *position.array]]), key) + def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: Optional[int] = None): + dataPoint = [value, *position.array] + if ID is not None: + dataPoint.append(ID) + self.logDataPointArray(np.array([dataPoint]), key) + + def _compileViews(self, views: List[View2D], detectedBy: str = None): + if detectedBy is None: + dataPerInteraction = self._data + if any(view.detectedBy for view in views): + utils.warn( + "Ignoring the detectedBy property of a view. Can only use detectedBy when adding a view after the " + "simulation with keep3D=True." + ) + else: + dataPerInteraction = self.getFiltered(detectedBy)._data - def _compileViews(self, views: List[View2D]): - for key, data in self._data.items(): + for key, data in dataPerInteraction.items(): datapointsContainer: Optional[ListArrayContainer] = data.dataPoints if datapointsContainer is None or len(datapointsContainer) == 0: continue @@ -304,6 +319,47 @@ def getDataPoints(self, key: InteractionKey, energyType=EnergyType.DEPOSITION) - return self._getData(DataType.DATA_POINT, key) + def filter(self, detectedBy: str) -> None: + """Keeps only the data points that were detected by the given solid label.""" + if not self._keep3D: + utils.warn("Cannot filter a logger that has discarded the 3D data.") + return + + filteredPhotonIDs = self._getPhotonIDs(InteractionKey(detectedBy)) + self._data = self._getDataForPhotons(filteredPhotonIDs) + self._outdatedViews = set(self._views) + + def getFiltered(self, detectedBy: str) -> EnergyLogger: + filteredLogger = EnergyLogger(self._scene, views=[]) + filteredPhotonIDs = self._getPhotonIDs(InteractionKey(detectedBy)) + filteredLogger._data = self._getDataForPhotons(filteredPhotonIDs) + return filteredLogger + + def _getPhotonIDs(self, key: InteractionKey = None) -> np.ndarray: + """Get all unique photon IDs that interacted with the given interaction key.""" + data = self.getRawDataPoints(key) + if data is None or data.shape[1] < 5: + return np.array([]) + return np.unique(data[:, 4].astype(int)) + + def _getDataForPhotons(self, photonIDs: np.ndarray) -> Dict[InteractionKey, InteractionData]: + keyToData: Dict[InteractionKey, InteractionData] = {} + photonIDs = np.asarray(photonIDs, dtype=np.uint32) + for key, interactionData in self._data.items(): + points: Optional[ListArrayContainer] = interactionData.dataPoints + if points is None: + continue + data = points.getData() + if data.shape[1] < 5: + continue + mask = np.isin(data[:, 4].astype(np.uint32), photonIDs) + filteredData = data[mask] + if filteredData.size > 0: + container = ListArrayContainer() + container.append(filteredData) + keyToData[key] = InteractionData(dataPoints=container) + return keyToData + def _fluenceTransform(self, key: InteractionKey, data: Optional[np.ndarray]) -> Optional[np.ndarray]: # Converts volumetric data to fluence rate when needed. if not key.volumetric or data is None: @@ -317,7 +373,7 @@ def export(self, exportPath: str): Export the raw 3D data points to a CSV file, along with the scene information to a JSON file. The data file .csv will be comma-delimited and will contain the following columns: - - energy, x, y, z, solid_index, surface_index + - energy, x, y, z, photon_index, solid_index, surface_index Two types of interactions are logged: scattering and surface crossings. In the first case, the energy will be the delta energy deposited at the point and the surface index will be -1. In the second case, the energy @@ -343,7 +399,7 @@ def export(self, exportPath: str): print("Exporting raw data to file...") filepath = f"{exportPath}.csv" with open(filepath, "w") as file: - file.write("energy,x,y,z,solid_index,surface_index\n") + file.write("energy,x,y,z,photon_index,solid_index,surface_index\n") self._writeKeyData(file, InteractionKey(WORLD_SOLID_LABEL), -1, -1) for i, solidLabel in enumerate(solidLabels): self._writeKeyData(file, InteractionKey(solidLabel), i, -1) @@ -382,7 +438,12 @@ def export(self, exportPath: str): def _writeKeyData(self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int): if key not in self._data or self._data[key].dataPoints is None: return - dataArray = self._data[key].dataPoints.getData().astype(str) - dataArray = np.hstack((dataArray, np.full((dataArray.shape[0], 2), str(solidIndex)))) - dataArray[:, 5] = str(surfaceIndex) - file.write("\n".join([",".join(row) for row in dataArray]) + "\n") + dataArray = self._data[key].dataPoints.getData() + n_rows = dataArray.shape[0] + + energy_xyz = dataArray[:, :4].astype(str) + photon_ids = dataArray[:, 4].astype(np.uint32).astype(str) + solid_indices = np.full(n_rows, solidIndex, dtype=int).astype(str) + surface_indices = np.full(n_rows, surfaceIndex, dtype=int).astype(str) + output = np.column_stack([energy_xyz, photon_ids, solid_indices, surface_indices]) + file.write("\n".join([",".join(row) for row in output]) + "\n") diff --git a/pytissueoptics/scene/logger/logger.py b/pytissueoptics/scene/logger/logger.py index 3b5c8d08..8f701bb0 100644 --- a/pytissueoptics/scene/logger/logger.py +++ b/pytissueoptics/scene/logger/logger.py @@ -79,7 +79,7 @@ def logPoint(self, point: Vector, key: InteractionKey = None): self._appendData([point.x, point.y, point.z], DataType.POINT, key) def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: Optional[int] = None): - dataPoint = [value, position.x, position.y, position.z] + dataPoint = [value, *position.array] if ID is not None: dataPoint.append(ID) self._appendData(dataPoint, DataType.DATA_POINT, key) From 7ced77ece20b51caf51f788ae49223b9cd250c16 Mon Sep 17 00:00:00 2001 From: JLBegin Date: Sat, 11 Oct 2025 20:12:17 -0400 Subject: [PATCH 11/14] optimize data export decrease RAM usage --- .../energyLogging/energyLogger.py | 51 +++++++++++-------- 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py index d21bd765..49445991 100644 --- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py +++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py @@ -368,11 +368,11 @@ def _fluenceTransform(self, key: InteractionKey, data: Optional[np.ndarray]) -> data[:, 0] = data[:, 0] / self._scene.getMaterial(key.solidLabel).mu_a return data - def export(self, exportPath: str): + def export(self, exportName: str): """ Export the raw 3D data points to a CSV file, along with the scene information to a JSON file. - The data file .csv will be comma-delimited and will contain the following columns: + The data file .csv will be comma-delimited and will contain the following columns: - energy, x, y, z, photon_index, solid_index, surface_index Two types of interactions are logged: scattering and surface crossings. In the first case, the energy will be @@ -380,7 +380,7 @@ def export(self, exportPath: str): will be the total photon energy when crossing the surface, either as positive if leaving the surface (along the normal) or as negative if entering the surface. - The scene information will be saved in a JSON file named .json, which includes details for each solid + The scene information will be saved in a JSON file named .json, which includes details for each solid index and surface index, such as their labels, materials, and geometry. The world information is also exported as solid index -1. """ @@ -397,7 +397,7 @@ def export(self, exportPath: str): solidLabels.sort() print("Exporting raw data to file...") - filepath = f"{exportPath}.csv" + filepath = f"{exportName}.csv" with open(filepath, "w") as file: file.write("energy,x,y,z,photon_index,solid_index,surface_index\n") self._writeKeyData(file, InteractionKey(WORLD_SOLID_LABEL), -1, -1) @@ -407,6 +407,31 @@ def export(self, exportPath: str): self._writeKeyData(file, InteractionKey(solidLabel, surfaceLabel), i, j) print(f"Exported data points to {filepath}") + self._exportSceneInfo(f"{exportName}.json", solidLabels) + + def _writeKeyData( + self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int + ): + if key not in self._data or self._data[key].dataPoints is None: + return + + dataArray = self._data[key].dataPoints.getData() + n_rows = dataArray.shape[0] + + output = np.empty((n_rows, 7), dtype=np.float64) + output[:, :4] = dataArray[:, :4] + output[:, 4] = dataArray[:, 4].astype(np.uint32) + output[:, 5] = solidIndex + output[:, 6] = surfaceIndex + + np.savetxt( + file, + output, + delimiter=",", + fmt=["%.8e", "%.8e", "%.8e", "%.8e", "%d", "%d", "%d"], + ) + + def _exportSceneInfo(self, filepath: str, solidLabels: List[str]): sceneInfo = {} material = self._scene.getWorldEnvironment().material sceneInfo["-1"] = {"label": "world", "material": material.__dict__ if material else None} @@ -430,20 +455,6 @@ def export(self, exportPath: str): "surfaces": surfaces, } - sceneFilepath = f"{exportPath}.json" - with open(sceneFilepath, "w") as file: + with open(filepath, "w") as file: json.dump(sceneInfo, file, indent=4) - print(f"Exported scene information to {sceneFilepath}") - - def _writeKeyData(self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int): - if key not in self._data or self._data[key].dataPoints is None: - return - dataArray = self._data[key].dataPoints.getData() - n_rows = dataArray.shape[0] - - energy_xyz = dataArray[:, :4].astype(str) - photon_ids = dataArray[:, 4].astype(np.uint32).astype(str) - solid_indices = np.full(n_rows, solidIndex, dtype=int).astype(str) - surface_indices = np.full(n_rows, surfaceIndex, dtype=int).astype(str) - output = np.column_stack([energy_xyz, photon_ids, solid_indices, surface_indices]) - file.write("\n".join([",".join(row) for row in output]) + "\n") + print(f"Exported scene information to {filepath}") From db2eff486bb5594c4b5a88e32f8d3e5c92122cea Mon Sep 17 00:00:00 2001 From: JLBegin Date: Sat, 11 Oct 2025 20:51:55 -0400 Subject: [PATCH 12/14] fix tests and lint --- pytissueoptics/examples/rayscattering/ex06.py | 16 ++-- .../rayscattering/display/viewer.py | 4 +- .../energyLogging/energyLogger.py | 6 +- pytissueoptics/rayscattering/photon.py | 4 +- .../rayscattering/tests/display/testViewer.py | 80 ++++++++----------- .../tests/energyLogging/testEnergyLogger.py | 21 ++--- .../tests/opencl/src/testCLFresnel.py | 2 + .../tests/opencl/src/testCLPhoton.py | 32 ++++---- .../tests/opencl/testCLKeyLog.py | 34 ++++---- .../rayscattering/tests/testPhoton.py | 24 +++--- pytissueoptics/scene/__init__.py | 4 +- .../scene/intersection/intersectionFinder.py | 4 +- pytissueoptics/scene/solids/__init__.py | 4 +- pytissueoptics/scene/solids/circle.py | 20 +++-- pytissueoptics/scene/solids/rectangle.py | 22 ++--- pytissueoptics/scene/solids/solid.py | 4 +- 16 files changed, 140 insertions(+), 141 deletions(-) diff --git a/pytissueoptics/examples/rayscattering/ex06.py b/pytissueoptics/examples/rayscattering/ex06.py index c0f8bceb..46b3c815 100644 --- a/pytissueoptics/examples/rayscattering/ex06.py +++ b/pytissueoptics/examples/rayscattering/ex06.py @@ -1,5 +1,3 @@ -import numpy as np - import env # noqa: F401 from pytissueoptics import * # noqa: F403 @@ -10,22 +8,22 @@ def exampleCode(): + import numpy as np N = 1000000 if hardwareAccelerationIsAvailable() else 1000 material = ScatteringMaterial(mu_s=10, mu_a=1, g=0.98, n=1.0) - cube = Cuboid( - a=3, b=3, c=3, position=Vector(0, 0, 0), material=material, label="cube" - ) + cube = Cuboid(a=3, b=3, c=3, position=Vector(0, 0, 0), material=material, label="cube") detector = Circle( - radius=0.25, orientation=Vector(0, 0, -1), position=Vector(0, 0.5, 1.501), label="detector", + radius=0.25, + orientation=Vector(0, 0, -1), + position=Vector(0, 0.5, 1.501), + label="detector", ).asDetector(halfAngle=np.pi / 4) scene = ScatteringScene([cube, detector]) logger = EnergyLogger(scene) - source = DirectionalSource( - position=Vector(0, -0.5, -2), direction=Vector(0, 0, 1), N=N, diameter=0.5 - ) + source = DirectionalSource(position=Vector(0, -0.5, -2), direction=Vector(0, 0, 1), N=N, diameter=0.5) source.propagate(scene, logger) diff --git a/pytissueoptics/rayscattering/display/viewer.py b/pytissueoptics/rayscattering/display/viewer.py index fc793f3e..276c908e 100644 --- a/pytissueoptics/rayscattering/display/viewer.py +++ b/pytissueoptics/rayscattering/display/viewer.py @@ -219,7 +219,9 @@ def _addPointCloud(self, style: PointCloudStyle): logger = self._logger if style.detectedBy is None else self._logger.getFiltered(style.detectedBy) pointCloud = PointCloudFactory(logger).getPointCloud( - solidLabel=style.solidLabel, surfaceLabel=style.surfaceLabel, energyType=style.energyType, + solidLabel=style.solidLabel, + surfaceLabel=style.surfaceLabel, + energyType=style.energyType, ) self._drawPointCloudOfSolids(pointCloud, style) diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py index 49445991..5cd023d2 100644 --- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py +++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py @@ -3,7 +3,7 @@ import json import os import pickle -from typing import List, Optional, TextIO, Union, Dict +from typing import Dict, List, Optional, TextIO, Union import numpy as np @@ -409,9 +409,7 @@ def export(self, exportName: str): self._exportSceneInfo(f"{exportName}.json", solidLabels) - def _writeKeyData( - self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int - ): + def _writeKeyData(self, file: TextIO, key: InteractionKey, solidIndex: int, surfaceIndex: int): if key not in self._data or self._data[key].dataPoints is None: return diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py index a7456f42..38c6a860 100644 --- a/pytissueoptics/rayscattering/photon.py +++ b/pytissueoptics/rayscattering/photon.py @@ -30,7 +30,7 @@ def __init__(self, position: Vector, direction: Vector, ID: int = 0): self._intersectionFinder: Optional[IntersectionFinder] = None self._logger: Optional[Logger] = None - self._lastIntersectedDetector: Optional[str] = None # Prevent intersecting with a detector when passing through it. + self._lastIntersectedDetector: Optional[str] = None @property def isAlive(self) -> bool: @@ -99,7 +99,7 @@ def step(self, distance=0) -> float: # Prevent re-intersecting with the same detector when passing through it. self._lastIntersectedDetector = intersection.insideEnvironment.solid.getLabel() - # skipping vertex check for now. TODO? Unclear stepCorrection and implications. + # skipping vertex check if detector return intersection.distanceLeft else: distanceLeft = self.reflectOrRefract(intersection) diff --git a/pytissueoptics/rayscattering/tests/display/testViewer.py b/pytissueoptics/rayscattering/tests/display/testViewer.py index 468c4531..608276e2 100644 --- a/pytissueoptics/rayscattering/tests/display/testViewer.py +++ b/pytissueoptics/rayscattering/tests/display/testViewer.py @@ -8,7 +8,7 @@ from pytissueoptics.rayscattering.display.profiles import Profile1D, ProfileFactory from pytissueoptics.rayscattering.display.viewer import PointCloudStyle, Viewer, Visibility from pytissueoptics.rayscattering.display.views import View2D -from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloud, PointCloudFactory +from pytissueoptics.rayscattering.energyLogging import EnergyLogger, PointCloud from pytissueoptics.rayscattering.scatteringScene import ScatteringScene from pytissueoptics.rayscattering.source import Source from pytissueoptics.scene.geometry import BoundingBox @@ -112,14 +112,11 @@ def testWhenShow3DWithSource_shouldDisplaySource(self): self.mock3DViewer.show.assert_called_once() def testWhenShow3DWithDefaultPointCloud_shouldDisplayPointCloudOfSolidsAndSurfaceLeaving(self): - mockPointCloudFactory = mock(PointCloudFactory) aPointCloud = PointCloud( - solidPoints=np.array([[0.5, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0], [-1, 0, 0, 0]]) + solidPoints=np.array([[0.5, 0, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0, 0], [-1, 0, 0, 0, 0]]) ) - when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud) - self.viewer._pointCloudFactory = mockPointCloudFactory - - self.viewer.show3D(visibility=Visibility.POINT_CLOUD) + with self._mockPointCloud(aPointCloud): + self.viewer.show3D(visibility=Visibility.POINT_CLOUD) self.mock3DViewer.addDataPoints.assert_called() addedSolidPoints = self.mock3DViewer.addDataPoints.call_args_list[0][0][0] @@ -130,25 +127,21 @@ def testWhenShow3DWithDefaultPointCloud_shouldDisplayPointCloudOfSolidsAndSurfac self.mock3DViewer.show.assert_called_once() def testGivenNoData_whenShow3DWithPointCloud_shouldNotDisplayPointCloud(self): - mockPointCloudFactory = mock(PointCloudFactory) aPointCloud = PointCloud() - when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud) - self.viewer._pointCloudFactory = mockPointCloudFactory - - self.viewer.show3D(visibility=Visibility.POINT_CLOUD) + with self._mockPointCloud(aPointCloud): + self.viewer.show3D(visibility=Visibility.POINT_CLOUD) self.mock3DViewer.addDataPoints.assert_not_called() self.mock3DViewer.show.assert_called_once() def testWhenShow3DWithSurfacePointCloud_shouldOnlyDisplaySurfacePoints(self): - mockPointCloudFactory = mock(PointCloudFactory) aPointCloud = PointCloud( - solidPoints=np.array([[0.5, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0], [-1, 0, 0, 0]]) + solidPoints=np.array([[0.5, 0, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0, 0], [-1, 0, 0, 0, 0]]) ) - when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud) - self.viewer._pointCloudFactory = mockPointCloudFactory - - self.viewer.show3D(visibility=Visibility.POINT_CLOUD, pointCloudStyle=PointCloudStyle(showSolidPoints=False)) + with self._mockPointCloud(aPointCloud): + self.viewer.show3D( + visibility=Visibility.POINT_CLOUD, pointCloudStyle=PointCloudStyle(showSolidPoints=False) + ) self.mock3DViewer.addDataPoints.assert_called_once() self.assertTrue( @@ -157,19 +150,16 @@ def testWhenShow3DWithSurfacePointCloud_shouldOnlyDisplaySurfacePoints(self): self.mock3DViewer.show.assert_called_once() def testWhenShow3DWithEnteringSurfacePointCloud_shouldOnlyDisplayEnteringSurfacePoints(self): - mockPointCloudFactory = mock(PointCloudFactory) aPointCloud = PointCloud( - solidPoints=np.array([[0.5, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0], [-1, 1, 1, 1]]) - ) - when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud) - self.viewer._pointCloudFactory = mockPointCloudFactory - - self.viewer.show3D( - visibility=Visibility.POINT_CLOUD, - pointCloudStyle=PointCloudStyle( - showSolidPoints=False, showSurfacePointsLeaving=False, showSurfacePointsEntering=True - ), + solidPoints=np.array([[0.5, 0, 0, 0, 0]]), surfacePoints=np.array([[1, 0, 0, 0, 0], [-1, 1, 1, 0, 1]]) ) + with self._mockPointCloud(aPointCloud): + self.viewer.show3D( + visibility=Visibility.POINT_CLOUD, + pointCloudStyle=PointCloudStyle( + showSolidPoints=False, showSurfacePointsLeaving=False, showSurfacePointsEntering=True + ), + ) self.mock3DViewer.addDataPoints.assert_called_once() self.assertTrue( @@ -206,16 +196,13 @@ def testWhenShow3DWithViewsIndexList_shouldAdd2DImageOfTheseViewsInThe3DDisplay( self.mock3DViewer.show.assert_called_once() def testGiven3DLogger_whenShow3DDefault_shouldDisplayEverythingExceptViews(self): - mockPointCloudFactory = mock(PointCloudFactory) aPointCloud = PointCloud() - when(mockPointCloudFactory).getPointCloud(...).thenReturn(aPointCloud) - self.viewer._pointCloudFactory = mockPointCloudFactory - self.viewer.show3D() + with self._mockPointCloud(aPointCloud): + self.viewer.show3D() verify(self.source, times=1).addToViewer(...) verify(self.scene, times=1).addToViewer(...) - verify(mockPointCloudFactory, times=1).getPointCloud(...) self.mock3DViewer.addImage.assert_not_called() self.mock3DViewer.show.assert_called_once() @@ -223,15 +210,11 @@ def testGiven2DLogger_whenShow3DDefault_shouldDisplayEverythingExceptPointCloud( self._givenLoggerWithXSceneView() self.logger.has3D = False - mockPointCloudFactory = mock(PointCloudFactory) - when(mockPointCloudFactory).getPointCloud(...).thenReturn() - self.viewer._pointCloudFactory = mockPointCloudFactory - - self.viewer.show3D() + with self._mockPointCloud(PointCloud()): + self.viewer.show3D() verify(self.source, times=1).addToViewer(...) verify(self.scene, times=1).addToViewer(...) - verify(mockPointCloudFactory, times=0).getPointCloud(...) self.mock3DViewer.addImage.assert_called() self.mock3DViewer.show.assert_called_once() @@ -239,16 +222,12 @@ def testGiven2DLogger_whenShow3DWithDefault3DVisibility_shouldWarnAndDisplayDefa self._givenLoggerWithXSceneView() self.logger.has3D = False - mockPointCloudFactory = mock(PointCloudFactory) - when(mockPointCloudFactory).getPointCloud(...).thenReturn() - self.viewer._pointCloudFactory = mockPointCloudFactory - - with self.assertWarns(UserWarning): - self.viewer.show3D(visibility=Visibility.DEFAULT_3D) + with self._mockPointCloud(PointCloud()): + with self.assertWarns(UserWarning): + self.viewer.show3D(visibility=Visibility.DEFAULT_3D) verify(self.source, times=1).addToViewer(...) verify(self.scene, times=1).addToViewer(...) - verify(mockPointCloudFactory, times=0).getPointCloud(...) self.mock3DViewer.addImage.assert_called() self.mock3DViewer.show.assert_called_once() @@ -260,3 +239,10 @@ def _givenLoggerWithXSceneView(self): self.logger.views = [sceneView] when(self.logger).updateView(sceneView).thenReturn() when(self.scene).getBoundingBox().thenReturn(BoundingBox([-2, 2], [-2, 2], [0, 5])) + + @staticmethod + def _mockPointCloud(pointCloud: PointCloud): + return patch( + "pytissueoptics.rayscattering.display.viewer.PointCloudFactory.getPointCloud", + return_value=pointCloud, + ) diff --git a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py index 4c41c438..7b4a6f2d 100644 --- a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py +++ b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py @@ -296,10 +296,10 @@ def testWhenExport_shouldExport3DDataPointsToFile(self): self.logger = EnergyLogger(scene) # Log entering surface event, world scattering event and scattering event in both solids. - self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("middleLayer")) - self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("frontLayer", "interface1")) - self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere")) - self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL)) + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("middleLayer"), ID=1) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("frontLayer", "interface1"), ID=0) + self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=0) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0) with tempfile.TemporaryDirectory() as tempDir: filePath = os.path.join(tempDir, "test_sim") @@ -309,12 +309,15 @@ def testWhenExport_shouldExport3DDataPointsToFile(self): with open(filePath + ".csv", "r") as f: lines = f.readlines() + def parse_line(line): + return list(map(float, line.strip().split(","))) + self.assertEqual(5, len(lines)) - self.assertEqual("energy,x,y,z,solid_index,surface_index\n", lines[0]) - self.assertEqual("0.2,0.0,0.0,0.0,-1,-1\n", lines[1]) - self.assertEqual("-0.9,0.5,1.0,0.75,1,5\n", lines[2]) - self.assertEqual("0.1,0.7,0.8,0.8,2,-1\n", lines[3]) - self.assertEqual("0.4,0.0,5.0,0.0,3,-1\n", lines[4]) + self.assertEqual("energy,x,y,z,photon_index,solid_index,surface_index\n", lines[0]) + self.assertEqual(parse_line(lines[1]), [0.2, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0]) + self.assertEqual(parse_line(lines[2]), [-0.9, 0.5, 1.0, 0.75, 0.0, 1.0, 5.0]) + self.assertEqual(parse_line(lines[3]), [0.1, 0.7, 0.8, 0.8, 1.0, 2.0, -1.0]) + self.assertEqual(parse_line(lines[4]), [0.4, 0.0, 5.0, 0.0, 0.0, 3.0, -1.0]) def testWhenExport_shouldExportMetadataToFile(self): scene = PhantomTissue(worldMaterial=ScatteringMaterial(0.1, 0.1, 0.99)) diff --git a/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py b/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py index 0c34f98b..afcaef63 100644 --- a/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py +++ b/pytissueoptics/rayscattering/tests/opencl/src/testCLFresnel.py @@ -55,6 +55,8 @@ def _setUpWith(self, n1=1.0, n2=1.5, normal=Vector(0, 0, 1)): self.INSIDE_SOLID_ID, self.OUTSIDE_SOLID_ID, False, + False, + 0, ) ] ) diff --git a/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py b/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py index 3784b2d1..92dccd66 100644 --- a/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py +++ b/pytissueoptics/rayscattering/tests/opencl/src/testCLPhoton.py @@ -21,7 +21,7 @@ VertexCL, ) from pytissueoptics.rayscattering.opencl.CLProgram import CLProgram -from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_ID, NO_SURFACE_ID, CLScene +from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, NO_SURFACE_ID, WORLD_SOLID_ID, CLScene from pytissueoptics.rayscattering.opencl.config.CLConfig import OPENCL_SOURCE_DIR from pytissueoptics.scene.geometry import Vertex @@ -175,7 +175,7 @@ def testWhenLogIntersectionWithPhotonLeavingIntoWorld_shouldLogOnlyOneIntersecti outsideSolidID = WORLD_SOLID_ID surfaceID = 0 - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False, False, 0)]) logger = DataPointCL(2) self._photonFunc("logIntersection", intersectionNormal, 0, surfaces, logger, 0) @@ -202,7 +202,7 @@ def testWhenLogIntersectionWithPhotonEnteringFromWorld_shouldLogOnlyOneIntersect outsideSolidID = WORLD_SOLID_ID surfaceID = 0 - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False, False, 0)]) logger = DataPointCL(2) self._photonFunc("logIntersection", intersectionNormal, surfaceID, surfaces, logger, 0) @@ -225,7 +225,7 @@ def testWhenLogIntersectionBetweenTwoSolids_shouldLogIntersectionOnEachSolidWith insideSolidID = self.INITIAL_SOLID_ID outsideSolidID = self.INITIAL_SOLID_ID + 11 - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID, False, False, 0)]) logger = DataPointCL(2) self._photonFunc("logIntersection", intersectionNormal, 0, surfaces, logger, 0) @@ -261,7 +261,7 @@ def testWhenReflectOrRefractWithReflectingIntersection_shouldReflectPhoton(self) self._mockFresnelIntersection(isReflected=True, incidencePlane=Vector(0, 0, 1), angleDeflection=np.pi / 2) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=WORLD_SOLID_ID, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, WORLD_SOLID_ID, False, False, 0)]) photonResult = self._photonFunc( "reflectOrRefract", intersectionNormal, @@ -301,7 +301,7 @@ def testWhenReflectOrRefractWithRefractingIntersection_shouldRefractPhoton(self) ) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, outsideSolidID=WORLD_SOLID_ID, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID, WORLD_SOLID_ID, False, False, 0)]) photonResult = self._photonFunc( "reflectOrRefract", intersectionNormal, @@ -334,7 +334,7 @@ def testWhenStepToInfinityWithNoIntersection_shouldKillPhoton(self): self._mockFindIntersection(exists=False) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)]) photonResult = self._photonFunc( "propagateStep", stepDistance, @@ -354,7 +354,7 @@ def testWhenStepWithNoIntersection_shouldMovePhotonAndScatter(self): self._mockFindIntersection(exists=False) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)]) photonResult = self._photonFunc( "propagateStep", stepDistance, @@ -381,9 +381,7 @@ def testWhenStepWithIntersectionTooCloseToVertex_shouldMovePhotonAwayABitAndScat intersectionPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * stepDistance logger = DataPointCL(2) - surfaces = SurfaceCL( - [SurfaceCLInfo(0, 0, 0, 0, insideSolidID=nextSolidID, outsideSolidID=self.INITIAL_SOLID_ID, toSmooth=False)] - ) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, nextSolidID, self.INITIAL_SOLID_ID, False, False, 0)]) triangles = TriangleCL([TriangleCLInfo([0, 1, 2], normal)]) # Put vertex slightly after the intersection point vertex = Vertex(intersectionPosition.x, intersectionPosition.y, intersectionPosition.z - 2e-7) @@ -410,7 +408,7 @@ def testWhenStepWithNoDistance_shouldStepWithANewScatteringDistance(self): self._mockFindIntersection(exists=False) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)]) photonResult = self._photonFunc( "propagateStep", stepDistance, @@ -432,7 +430,7 @@ def testWhenStepWithIntersectionReflecting_shouldMovePhotonToIntersection(self): self._mockFresnelIntersection(isReflected=True) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)]) triangles = TriangleCL([TriangleCLInfo([0, 1, 2], Vector(0, 0, 1))]) vertices = VertexCL([Vertex(0, 0, 0)] * 3) photonResult = self._photonFunc( @@ -457,7 +455,7 @@ def testWhenStepWithIntersectionRefracting_shouldMovePhotonToIntersection(self): self._mockFresnelIntersection(isReflected=False) logger = DataPointCL(2) - surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, insideSolidID=9, outsideSolidID=10, toSmooth=False)]) + surfaces = SurfaceCL([SurfaceCLInfo(0, 0, 0, 0, 9, 10, False, False, 0)]) triangles = TriangleCL([TriangleCLInfo([0, 1, 2], Vector(0, 0, 1))]) vertices = VertexCL([Vertex(0, 0, 0)] * 3) photonResult = self._photonFunc( @@ -590,7 +588,7 @@ def _getPhotonResult(self, photonBuffer: PhotonCL): def _getDataPointResult(self, dataPointBuffer: DataPointCL, i=0): data = self.program.getData(dataPointBuffer)[i] - return DataPointResult(deltaWeight=data[0], position=Vector(*data[1:4]), solidID=data[4], surfaceID=data[5]) + return DataPointResult(deltaWeight=data[0], position=Vector(*data[1:4]), solidID=data[5], surfaceID=data[6]) def _assertVectorAlmostEqual(self, v1: Vector, v2: Vector, places=7): self.assertAlmostEqual(v1.x, v2.x, places=places) @@ -644,9 +642,7 @@ def _mockFresnelIntersection( def _mockFindIntersection(self, exists=True, distance=8.0, normal=Vector(0, 0, 1), surfaceID=0, distanceLeft=2): expectedPosition = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance - intersectionCall = ( - """Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID);""" - ) + intersectionCall = """Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID, photons[photonID].lastIntersectedDetectorID);""" px, py, pz = expectedPosition.array nx, ny, nz = normal.array mockCall = """Intersection intersection; diff --git a/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py b/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py index ca75b726..5f81db63 100644 --- a/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py +++ b/pytissueoptics/rayscattering/tests/opencl/testCLKeyLog.py @@ -5,7 +5,13 @@ from pytissueoptics import Cube, EnergyLogger, ScatteringMaterial, ScatteringScene, Sphere from pytissueoptics.rayscattering.opencl import OPENCL_OK -from pytissueoptics.rayscattering.opencl.CLScene import NO_LOG_ID, WORLD_SOLID_ID, WORLD_SOLID_LABEL, NO_SURFACE_ID, CLScene +from pytissueoptics.rayscattering.opencl.CLScene import ( + NO_LOG_ID, + NO_SURFACE_ID, + WORLD_SOLID_ID, + WORLD_SOLID_LABEL, + CLScene, +) from pytissueoptics.rayscattering.opencl.utils import CLKeyLog from pytissueoptics.scene.logger import InteractionKey @@ -27,13 +33,13 @@ def _createTestLog(self, sceneCL): log = np.array( [ - [0, 0, 0, 0, cubeID, NO_SURFACE_ID], - [2, 0, 0, 0, sphereID, NO_SURFACE_ID], - [3, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], - [4, 9, 9, 9, NO_LOG_ID, 9], - [5, 0, 0, 0, cubeID, cubeSurfaceIDs[1]], - [6, 0, 0, 0, sphereID, sphereSurfaceIDs[1]], - [1, 0, 0, 0, cubeID, NO_SURFACE_ID], + [0, 0, 0, 0, 0, cubeID, NO_SURFACE_ID], + [2, 0, 0, 0, 0, sphereID, NO_SURFACE_ID], + [3, 0, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], + [4, 9, 9, 9, 0, NO_LOG_ID, 9], + [5, 0, 0, 0, 0, cubeID, cubeSurfaceIDs[1]], + [6, 0, 0, 0, 0, sphereID, sphereSurfaceIDs[1]], + [1, 0, 0, 0, 0, cubeID, NO_SURFACE_ID], ] ) return log @@ -49,7 +55,7 @@ def testGivenCLKeyLog_whenTransferToSceneLogger_shouldLogDataWithInteractionKeys verify(sceneLogger, times=5).logDataPointArray(...) - expectedCubeData = arg_that(lambda arg: np.array_equal(arg, np.array([[0, 0, 0, 0], [1, 0, 0, 0]]))) + expectedCubeData = arg_that(lambda arg: np.array_equal(arg, np.array([[0, 0, 0, 0, 0], [1, 0, 0, 0, 0]]))) verify(sceneLogger).logDataPointArray(expectedCubeData, InteractionKey(self.cube.getLabel())) expectedValuesWithKeys = [ @@ -59,7 +65,7 @@ def testGivenCLKeyLog_whenTransferToSceneLogger_shouldLogDataWithInteractionKeys (6, InteractionKey(self.sphere.getLabel(), self.sphere.surfaceLabels[0])), ] for value, expectedKey in expectedValuesWithKeys: - expectedData = arg_that(lambda arg: np.array_equal(arg, np.array([[value, 0, 0, 0]]))) + expectedData = arg_that(lambda arg: np.array_equal(arg, np.array([[value, 0, 0, 0, 0]]))) verify(sceneLogger).logDataPointArray(expectedData, expectedKey) def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWithInteractionKeys(self): @@ -67,9 +73,9 @@ def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWit sceneCL = CLScene(self.scene, nWorkUnits=10) log = np.array( [ - [1, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], - [2, 0, 0, 0, NO_LOG_ID, 99], - [3, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], + [1, 0, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], + [2, 0, 0, 0, 0, NO_LOG_ID, 99], + [3, 0, 0, 0, 0, WORLD_SOLID_ID, NO_SURFACE_ID], ] ) clKeyLog = CLKeyLog(log, sceneCL) @@ -79,5 +85,5 @@ def testGivenCLKeyLogForInfiniteScene_whenTransferToSceneLogger_shouldLogDataWit clKeyLog.toSceneLogger(sceneLogger) verify(sceneLogger, times=1).logDataPointArray(...) - expectedWorldData = arg_that(lambda arg: np.array_equal(arg, np.array([[1, 0, 0, 0], [3, 0, 0, 0]]))) + expectedWorldData = arg_that(lambda arg: np.array_equal(arg, np.array([[1, 0, 0, 0, 0], [3, 0, 0, 0, 0]]))) verify(sceneLogger).logDataPointArray(expectedWorldData, InteractionKey(WORLD_SOLID_LABEL)) diff --git a/pytissueoptics/rayscattering/tests/testPhoton.py b/pytissueoptics/rayscattering/tests/testPhoton.py index 974b3220..90f96817 100644 --- a/pytissueoptics/rayscattering/tests/testPhoton.py +++ b/pytissueoptics/rayscattering/tests/testPhoton.py @@ -312,7 +312,7 @@ def testGivenALogger_whenSteppingOutsideASolid_shouldLogPositiveWeightOnSurfaceI intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance interactionKey = InteractionKey(self.SOLID_INSIDE_LABEL, self.SURFACE_LABEL) - verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey) + verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey, 0) def testGivenALogger_whenSteppingOutsideASolid_shouldLogNegativeWeightOnOtherSolidSurfaceIntersection(self): distance = 8 @@ -320,7 +320,7 @@ def testGivenALogger_whenSteppingOutsideASolid_shouldLogNegativeWeightOnOtherSol intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance interactionKey = InteractionKey(self.SOLID_OUTSIDE_LABEL, self.SURFACE_LABEL) - verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey) + verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey, 0) def givenALogger_whenSteppingInsideASolidAt(self, distance): logger = self._createLogger() @@ -340,7 +340,7 @@ def testGivenALogger_whenSteppingInsideASurface_shouldLogNegativeWeightAtInterse intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance interactionKey = InteractionKey(self.SOLID_INSIDE_LABEL, self.SURFACE_LABEL) - verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey) + verify(logger).logDataPoint(-self.photon.weight, intersectionPoint, interactionKey, 0) def testGivenALogger_whenSteppingInsideASurface_shouldLogPositiveWeightOnOtherSolidSurfaceIntersection(self): distance = 8 @@ -348,7 +348,7 @@ def testGivenALogger_whenSteppingInsideASurface_shouldLogPositiveWeightOnOtherSo intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance interactionKey = InteractionKey(self.SOLID_OUTSIDE_LABEL, self.SURFACE_LABEL) - verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey) + verify(logger).logDataPoint(self.photon.weight, intersectionPoint, interactionKey, 0) def testGivenALogger_whenScatter_shouldLogWeightLossAtThisPositionWithSolidLabel(self): solidInside = mock(Solid) @@ -359,7 +359,7 @@ def testGivenALogger_whenScatter_shouldLogWeightLossAtThisPositionWithSolidLabel self.photon.scatter() weightLoss = self.photon.material.getAlbedo() - verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(self.SOLID_INSIDE_LABEL)) + verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(self.SOLID_INSIDE_LABEL), 0) def testGivenALogger_whenScatterInWorldMaterial_shouldLogWeightLossAtThisPositionWithWorldLabel(self): logger = self._createLogger() @@ -368,7 +368,7 @@ def testGivenALogger_whenScatterInWorldMaterial_shouldLogWeightLossAtThisPositio self.photon.scatter() weightLoss = self.photon.material.getAlbedo() - verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(WORLD_LABEL)) + verify(logger).logDataPoint(weightLoss, self.INITIAL_POSITION, InteractionKey(WORLD_LABEL), 0) def givenALoggerAndNoSolidOutside_whenSteppingInsideASolidAt(self, distance): self.solidOutside = None @@ -440,9 +440,7 @@ def testWhenIntersectingDetectorInsideNA_shouldKillPhoton(self): # Setup intersection with detector. distance = 8 intersectionFinder = self._createIntersectionFinder(distance) - self.photon.setContext( - Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder - ) + self.photon.setContext(Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder) self.photon.step(distance + 2) @@ -467,9 +465,7 @@ def testWhenIntersectingDetectorOutsideNA_shouldNotKillPhoton(self): # Setup intersection with detector. distance = 8 intersectionFinder = self._createIntersectionFinder(distance) - self.photon.setContext( - Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder - ) + self.photon.setContext(Environment(ScatteringMaterial()), intersectionFinder=intersectionFinder) self.photon.step(distance + 2) @@ -491,7 +487,7 @@ def testGivenALogger_whenDetected_shouldLogWeightAtThisPositionWithDetectorLabel # Setup photon going towards detector within NA. photonIncidence = halfAngle # Within NA self.INITIAL_DIRECTION = Vector(0, math.sin(photonIncidence), -math.cos(photonIncidence)) - self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy()) + self.photon = Photon(self.INITIAL_POSITION.copy(), self.INITIAL_DIRECTION.copy(), ID=42) initial_weight = self.photon.weight # Setup intersection with detector. @@ -508,7 +504,7 @@ def testGivenALogger_whenDetected_shouldLogWeightAtThisPositionWithDetectorLabel self.assertFalse(self.photon.isAlive) intersectionPoint = self.INITIAL_POSITION + self.INITIAL_DIRECTION * distance - verify(logger).logDataPoint(initial_weight, intersectionPoint, InteractionKey(detectorLabel)) + verify(logger).logDataPoint(initial_weight, intersectionPoint, InteractionKey(detectorLabel), 42) def assertVectorEqual(self, v1, v2): self.assertAlmostEqual(v1.x, v2.x, places=7) diff --git a/pytissueoptics/scene/__init__.py b/pytissueoptics/scene/__init__.py index bc7a54cf..7d9ea20c 100644 --- a/pytissueoptics/scene/__init__.py +++ b/pytissueoptics/scene/__init__.py @@ -4,6 +4,7 @@ from .material import RefractiveMaterial from .scene import Scene from .solids import ( + Circle, Cone, Cube, Cuboid, @@ -11,11 +12,10 @@ Ellipsoid, PlanoConcaveLens, PlanoConvexLens, + Rectangle, Sphere, SymmetricLens, ThickLens, - Rectangle, - Circle, ) from .viewer import ViewPointStyle, get3DViewer diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py index 2a9af76c..f66a98bf 100644 --- a/pytissueoptics/scene/intersection/intersectionFinder.py +++ b/pytissueoptics/scene/intersection/intersectionFinder.py @@ -182,7 +182,9 @@ def _findIntersection( # todo: test if this node search is still compatible with the new core intersection logic # which can now require testing a polygon slightly behind the ray origin. if node.isLeaf: - intersection = self._findClosestPolygonIntersection(ray, node.polygons, self._currentSolidLabel, ignoreLabel) + intersection = self._findClosestPolygonIntersection( + ray, node.polygons, self._currentSolidLabel, ignoreLabel + ) return intersection if not self._nodeIsWorthExploring(ray, node, closestDistance): diff --git a/pytissueoptics/scene/solids/__init__.py b/pytissueoptics/scene/solids/__init__.py index bb39b960..b4e464f8 100644 --- a/pytissueoptics/scene/solids/__init__.py +++ b/pytissueoptics/scene/solids/__init__.py @@ -1,14 +1,14 @@ +from .circle import Circle from .cone import Cone from .cube import Cube from .cuboid import Cuboid from .cylinder import Cylinder from .ellipsoid import Ellipsoid from .lens import PlanoConcaveLens, PlanoConvexLens, SymmetricLens, ThickLens +from .rectangle import Rectangle from .solid import Solid from .solidFactory import SolidFactory from .sphere import Sphere -from .circle import Circle -from .rectangle import Rectangle __all__ = [ "Solid", diff --git a/pytissueoptics/scene/solids/circle.py b/pytissueoptics/scene/solids/circle.py index 26b4ae32..7cfba17e 100644 --- a/pytissueoptics/scene/solids/circle.py +++ b/pytissueoptics/scene/solids/circle.py @@ -1,11 +1,21 @@ -from ..geometry import primitives, Vector +from ..geometry import Vector, primitives from .cylinder import Cylinder class Circle(Cylinder): - def __init__(self, radius: float, orientation: Vector = Vector(0, 0, 1), position: Vector = Vector(0, 0, 0), - u: int = 32, s: int = 1, primitive: str = primitives.DEFAULT, label="circle"): - super().__init__(radius, length=0, u=u, v=1, s=s, position=position, label=label, primitive=primitive, smooth=False) + def __init__( + self, + radius: float, + orientation: Vector = Vector(0, 0, 1), + position: Vector = Vector(0, 0, 0), + u: int = 32, + s: int = 1, + primitive: str = primitives.DEFAULT, + label="circle", + ): + super().__init__( + radius, length=0, u=u, v=1, s=s, position=position, label=label, primitive=primitive, smooth=False + ) self.orient(towards=orientation) def _computeTriangleMesh(self): @@ -23,5 +33,5 @@ def _geometryParams(self) -> dict: "position": self._position, "orientation": self._orientation, "u": self._u, - "s": self._s + "s": self._s, } diff --git a/pytissueoptics/scene/solids/rectangle.py b/pytissueoptics/scene/solids/rectangle.py index 3158a259..ad7c753d 100644 --- a/pytissueoptics/scene/solids/rectangle.py +++ b/pytissueoptics/scene/solids/rectangle.py @@ -1,10 +1,17 @@ -from ..geometry import primitives, Vertex, Triangle, Quad, Vector +from ..geometry import Quad, Triangle, Vector, Vertex, primitives from .solid import Solid class Rectangle(Solid): - def __init__(self, a: float, b: float, orientation = Vector(0, 0, 1), position = Vector(0, 0, 0), - primitive: str = primitives.DEFAULT, label="rectangle"): + def __init__( + self, + a: float, + b: float, + orientation=Vector(0, 0, 1), + position=Vector(0, 0, 0), + primitive: str = primitives.DEFAULT, + label="rectangle", + ): self._a = a self._b = b @@ -12,7 +19,7 @@ def __init__(self, a: float, b: float, orientation = Vector(0, 0, 1), position = Vertex(-a / 2, -b / 2, 0), Vertex(a / 2, -b / 2, 0), Vertex(a / 2, b / 2, 0), - Vertex(-a / 2, b / 2, 0) + Vertex(-a / 2, b / 2, 0), ] super().__init__(vertices, position, label=label, primitive=primitive) @@ -27,9 +34,4 @@ def _computeQuadMesh(self): self._surfaces.add("surface", [Quad(V[0], V[1], V[2], V[3])]) def _geometryParams(self) -> dict: - return { - "a": self._a, - "b": self._b, - "position": self._position, - "orientation": self._orientation - } + return {"a": self._a, "b": self._b, "position": self._position, "orientation": self._orientation} diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py index 23ac28f8..70408e3f 100644 --- a/pytissueoptics/scene/solids/solid.py +++ b/pytissueoptics/scene/solids/solid.py @@ -379,6 +379,4 @@ def geometryExport(self) -> dict[str, Any]: def _geometryParams(self) -> dict: """To be implemented by Solid subclasses to detail other geometry parameters.""" - raise NotImplementedError( - f"Geometry parameters not implemented for Solids of type {type(self).__name__}" - ) + raise NotImplementedError(f"Geometry parameters not implemented for Solids of type {type(self).__name__}") From a61b5cc5db553c6471d84be551dbaf5b378388b0 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Sun, 12 Oct 2025 11:14:19 -0400 Subject: [PATCH 13/14] support filtering multiple detectors --- pytissueoptics/examples/rayscattering/ex06.py | 26 +++++++---- .../rayscattering/display/viewer.py | 2 +- .../display/views/defaultViews.py | 24 +++++----- .../rayscattering/display/views/view2D.py | 4 +- .../energyLogging/energyLogger.py | 44 ++++++++++++++----- 5 files changed, 67 insertions(+), 33 deletions(-) diff --git a/pytissueoptics/examples/rayscattering/ex06.py b/pytissueoptics/examples/rayscattering/ex06.py index 46b3c815..16d57071 100644 --- a/pytissueoptics/examples/rayscattering/ex06.py +++ b/pytissueoptics/examples/rayscattering/ex06.py @@ -9,21 +9,31 @@ def exampleCode(): import numpy as np + N = 1000000 if hardwareAccelerationIsAvailable() else 1000 material = ScatteringMaterial(mu_s=10, mu_a=1, g=0.98, n=1.0) cube = Cuboid(a=3, b=3, c=3, position=Vector(0, 0, 0), material=material, label="cube") - detector = Circle( + + detectorA = Circle( radius=0.25, orientation=Vector(0, 0, -1), - position=Vector(0, 0.5, 1.501), - label="detector", + position=Vector(0, 1, 1.501), + label="detectorA", ).asDetector(halfAngle=np.pi / 4) - scene = ScatteringScene([cube, detector]) + + detectorB = Circle( + radius=0.25, + orientation=Vector(0, 0, -1), + position=Vector(0, -1, 1.501), + label="detectorB", + ).asDetector(halfAngle=np.pi / 4) + + scene = ScatteringScene([cube, detectorA, detectorB]) logger = EnergyLogger(scene) - source = DirectionalSource(position=Vector(0, -0.5, -2), direction=Vector(0, 0, 1), N=N, diameter=0.5) + source = DirectionalSource(position=Vector(0, 0, -2), direction=Vector(0, 0, 1), N=N, diameter=0.5) source.propagate(scene, logger) @@ -31,13 +41,13 @@ def exampleCode(): viewer.reportStats() # Either filter the whole logger - logger.filter(detectedBy="detector") + logger.filter(detectedBy=["detectorA", "detectorB"]) viewer.show2D(View2DProjectionX()) viewer.show3D() # ... or keep the logger intact and filter individual views - # viewer.show2D(View2DProjectionX(detectedBy="detector")) - # viewer.show3D(pointCloudStyle=PointCloudStyle(detectedBy="detector")) + # viewer.show2D(View2DProjectionX(detectedBy="detectorA")) + # viewer.show3D(pointCloudStyle=PointCloudStyle(detectedBy=["detectorA", "detectorB"])) if __name__ == "__main__": diff --git a/pytissueoptics/rayscattering/display/viewer.py b/pytissueoptics/rayscattering/display/viewer.py index 276c908e..ba3b398f 100644 --- a/pytissueoptics/rayscattering/display/viewer.py +++ b/pytissueoptics/rayscattering/display/viewer.py @@ -66,7 +66,7 @@ def __init__( showSurfacePointsLeaving: bool = True, showSurfacePointsEntering: bool = False, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, showPointsAsSpheres: bool = False, pointSize: float = 0.15, scaleWithValue: bool = True, diff --git a/pytissueoptics/rayscattering/display/views/defaultViews.py b/pytissueoptics/rayscattering/display/views/defaultViews.py index 6a2e17ee..18f0ebfa 100644 --- a/pytissueoptics/rayscattering/display/views/defaultViews.py +++ b/pytissueoptics/rayscattering/display/views/defaultViews.py @@ -21,7 +21,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( projectionDirection, @@ -44,7 +44,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_X_VIEW_DIRECTIONS, @@ -63,7 +63,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_Y_VIEW_DIRECTIONS, @@ -82,7 +82,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_Z_VIEW_DIRECTIONS, @@ -104,7 +104,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( projectionDirection, @@ -138,7 +138,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_X_VIEW_DIRECTIONS, @@ -159,7 +159,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_Y_VIEW_DIRECTIONS, @@ -180,7 +180,7 @@ def __init__( surfaceEnergyLeaving: bool = True, limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_Z_VIEW_DIRECTIONS, @@ -204,7 +204,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( projectionDirection, @@ -240,7 +240,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_X_VIEW_DIRECTIONS, @@ -263,7 +263,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_Y_VIEW_DIRECTIONS, @@ -286,7 +286,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): super().__init__( *DEFAULT_Z_VIEW_DIRECTIONS, diff --git a/pytissueoptics/rayscattering/display/views/view2D.py b/pytissueoptics/rayscattering/display/views/view2D.py index 88f86769..025d5c8a 100644 --- a/pytissueoptics/rayscattering/display/views/view2D.py +++ b/pytissueoptics/rayscattering/display/views/view2D.py @@ -50,7 +50,7 @@ def __init__( limits: Tuple[Tuple[float, float], Tuple[float, float]] = None, binSize: Union[float, Tuple[int, int]] = None, energyType=EnergyType.DEPOSITION, - detectedBy: str = None, + detectedBy: Union[str, List[str]] = None, ): """ The 2D view plane is obtained by looking towards the 'projectionDirection'. The 'horizontalDirection' @@ -125,7 +125,7 @@ def energyType(self) -> EnergyType: return self._energyType @property - def detectedBy(self) -> Optional[str]: + def detectedBy(self) -> Union[str, List[str], None]: return self._detectedBy def extractData(self, dataPoints: np.ndarray): diff --git a/pytissueoptics/rayscattering/energyLogging/energyLogger.py b/pytissueoptics/rayscattering/energyLogging/energyLogger.py index 5cd023d2..db43b650 100644 --- a/pytissueoptics/rayscattering/energyLogging/energyLogger.py +++ b/pytissueoptics/rayscattering/energyLogging/energyLogger.py @@ -241,7 +241,7 @@ def logDataPoint(self, value: float, position: Vector, key: InteractionKey, ID: dataPoint.append(ID) self.logDataPointArray(np.array([dataPoint]), key) - def _compileViews(self, views: List[View2D], detectedBy: str = None): + def _compileViews(self, views: List[View2D], detectedBy: Union[str, List[str]] = None): if detectedBy is None: dataPerInteraction = self._data if any(view.detectedBy for view in views): @@ -319,28 +319,52 @@ def getDataPoints(self, key: InteractionKey, energyType=EnergyType.DEPOSITION) - return self._getData(DataType.DATA_POINT, key) - def filter(self, detectedBy: str) -> None: - """Keeps only the data points that were detected by the given solid label.""" + def filter(self, detectedBy: Union[str, List[str]]) -> None: + """Keeps only the data points from photons detected by one of the specified detector(s).""" if not self._keep3D: utils.warn("Cannot filter a logger that has discarded the 3D data.") return - filteredPhotonIDs = self._getPhotonIDs(InteractionKey(detectedBy)) + filteredPhotonIDs = self._getDetectedPhotonIDs(detectedBy) self._data = self._getDataForPhotons(filteredPhotonIDs) self._outdatedViews = set(self._views) - def getFiltered(self, detectedBy: str) -> EnergyLogger: + def getFiltered(self, detectedBy: Union[str, List[str]]) -> "EnergyLogger": + """ + Returns a new logger with only data from photons detected by one of the specified detector(s). + """ filteredLogger = EnergyLogger(self._scene, views=[]) - filteredPhotonIDs = self._getPhotonIDs(InteractionKey(detectedBy)) + filteredPhotonIDs = self._getDetectedPhotonIDs(detectedBy) filteredLogger._data = self._getDataForPhotons(filteredPhotonIDs) return filteredLogger - def _getPhotonIDs(self, key: InteractionKey = None) -> np.ndarray: - """Get all unique photon IDs that interacted with the given interaction key.""" + def _getDetectedPhotonIDs(self, detectedBy: Union[str, List[str]]) -> np.ndarray: + """Helper to get photon IDs detected by one of the specified detector(s).""" + detector_labels = [detectedBy] if isinstance(detectedBy, str) else detectedBy + detector_keys = [InteractionKey(label) for label in detector_labels] + photonIDs = self._getPhotonIDs(detector_keys) + if len(photonIDs) == 0: + utils.warn(f"No photons detected by: {detectedBy}") + return photonIDs + + def _getPhotonIDs(self, key: Union[InteractionKey, List[InteractionKey], None] = None) -> np.ndarray: + """Get all unique photon IDs that interacted with one of the given interaction key(s).""" + if isinstance(key, list): + all_photon_ids = [] + for k in key: + data = self.getRawDataPoints(k) + if data is not None and data.shape[1] >= 5: + all_photon_ids.append(data[:, 4].astype(np.uint32)) + + if len(all_photon_ids) == 0: + return np.array([], dtype=np.uint32) + combined = np.concatenate(all_photon_ids) + return np.unique(combined) + data = self.getRawDataPoints(key) if data is None or data.shape[1] < 5: - return np.array([]) - return np.unique(data[:, 4].astype(int)) + return np.array([], dtype=np.uint32) + return np.unique(data[:, 4].astype(np.uint32)) def _getDataForPhotons(self, photonIDs: np.ndarray) -> Dict[InteractionKey, InteractionData]: keyToData: Dict[InteractionKey, InteractionData] = {} From ec601c5f7a669a6a38a307d3e315ab7433851626 Mon Sep 17 00:00:00 2001 From: jlbegin Date: Sun, 12 Oct 2025 13:44:43 -0400 Subject: [PATCH 14/14] improve test coverage --- pyproject.toml | 13 +++- pytissueoptics/examples/rayscattering/ex06.py | 2 +- pytissueoptics/examples/scene/example2.py | 6 +- .../rayscattering/scatteringScene.py | 3 + .../tests/energyLogging/testEnergyLogger.py | 77 +++++++++++++++++++ .../tests/testScatteringScene.py | 10 +++ pytissueoptics/scene/solids/circle.py | 7 ++ pytissueoptics/scene/solids/rectangle.py | 4 + pytissueoptics/scene/solids/solid.py | 4 + .../intersection/testIntersectionFinder.py | 7 ++ .../scene/tests/solids/testCircle.py | 37 +++++++++ .../scene/tests/solids/testRectangle.py | 37 +++++++++ .../scene/tests/solids/testSolid.py | 18 +++++ 13 files changed, 220 insertions(+), 5 deletions(-) create mode 100644 pytissueoptics/scene/tests/solids/testCircle.py create mode 100644 pytissueoptics/scene/tests/solids/testRectangle.py diff --git a/pyproject.toml b/pyproject.toml index 17b1bc35..3b39c0ae 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,7 +27,7 @@ dependencies = [ "Pygments", "siphash24; python_version < '3.13'", "pyopencl", - "vtk==9.4", + "vtk>=9.4,<9.5", "mayavi-dev", "pyqt5", ] @@ -67,3 +67,14 @@ ignore = [ [tool.ruff.format] quote-style = "double" indent-style = "space" + +[tool.coverage.run] +source = ["pytissueoptics"] +omit = [ + "*/tests/*", + "*/test*.py", + "*/examples/*", +] + +[tool.coverage.report] +show_missing = true diff --git a/pytissueoptics/examples/rayscattering/ex06.py b/pytissueoptics/examples/rayscattering/ex06.py index 16d57071..9f1a74ae 100644 --- a/pytissueoptics/examples/rayscattering/ex06.py +++ b/pytissueoptics/examples/rayscattering/ex06.py @@ -4,7 +4,7 @@ TITLE = "Detectors" -DESCRIPTION = """ Sampling volume simulation with a directional source and a circle detector. """ +DESCRIPTION = """ Sampling volume simulation with a directional source and two circle detectors. """ def exampleCode(): diff --git a/pytissueoptics/examples/scene/example2.py b/pytissueoptics/examples/scene/example2.py index b58415ee..bb3d3117 100644 --- a/pytissueoptics/examples/scene/example2.py +++ b/pytissueoptics/examples/scene/example2.py @@ -9,9 +9,9 @@ def exampleCode(): from pytissueoptics.scene import Cuboid, Vector, get3DViewer - cuboid1 = Cuboid(1, 1, 1, position=Vector(2, 0, 0)) - cuboid2 = Cuboid(2, 1, 1, position=Vector(0, 2, 0)) - cuboid3 = Cuboid(3, 1, 1, position=Vector(0, 0, 2)) + cuboid1 = Cuboid(1, 1, 1, position=Vector(2, 0, 0), label="1") + cuboid2 = Cuboid(2, 1, 1, position=Vector(0, 2, 0), label="2") + cuboid3 = Cuboid(3, 1, 1, position=Vector(0, 0, 2), label="3") viewer = get3DViewer() viewer.add(cuboid1, cuboid2, cuboid3, representation="wireframe", lineWidth=5) diff --git a/pytissueoptics/rayscattering/scatteringScene.py b/pytissueoptics/rayscattering/scatteringScene.py index 17cd5ce2..b068d90a 100644 --- a/pytissueoptics/rayscattering/scatteringScene.py +++ b/pytissueoptics/rayscattering/scatteringScene.py @@ -13,6 +13,9 @@ def __init__(self, solids: List[Solid], worldMaterial=ScatteringMaterial(), igno super().__init__(solids, worldMaterial=worldMaterial, ignoreIntersections=ignoreIntersections) def add(self, solid: Solid, position: Vector = None): + if solid.isFlat and not solid.isDetector: + raise Exception(f"Solid '{solid.getLabel()}' is flat. Flat solids must be used with asDetector().") + polygonSample = solid.getPolygons()[0] if not isinstance(polygonSample.insideEnvironment.material, ScatteringMaterial) and not solid.isDetector: raise Exception( diff --git a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py index 7b4a6f2d..446d6e0d 100644 --- a/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py +++ b/pytissueoptics/rayscattering/tests/energyLogging/testEnergyLogger.py @@ -177,6 +177,17 @@ def testGiven2DLoggerWithData_whenShowWithCustomViewNotContainedByExistingViews_ self.assertEqual(initialNumberOfViews, len(self.logger.views)) + def testGiven3DLoggerWithData_whenAddViewWithDetectedBy_shouldInitializeViewWithFilteredData(self): + # Photon 0 + self.logger.logDataPoint(0.9, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0) + # Photon 1 interacted with "sphere". + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1) + self.logger.logDataPoint(0.4, Vector(0.5, 0.5, 0.5), InteractionKey("sphere"), ID=1) + + newView = View2DProjectionX(detectedBy="sphere") + self.logger.addView(newView) + self.assertAlmostEqual(0.1 + 0.4, newView.getSum()) + def testWhenSave_shouldSaveLoggerToFile(self): with tempfile.TemporaryDirectory() as tempDir: filePath = os.path.join(tempDir, "test.log") @@ -352,3 +363,69 @@ def testWhenExport_shouldExportMetadataToFile(self): self.assertEqual(expectedLayerGeometry, sceneInfo["0"]["geometry"]) self.assertEqual(16, len(sceneInfo["0"]["surfaces"])) self.assertEqual({"label": "interface0", "normal": [0.0, 0.0, -1.0]}, sceneInfo["0"]["surfaces"]["14"]) + + def testWhenFilterWithOneLabel_shouldOnlyKeepPhotonsThatInteractedWithTheGivenLabel(self): + # Photon 0 + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=0) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0) + # Photon 1 interacted with "sphere". + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=1) + self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=1) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=1) + + self.logger.filter(detectedBy="sphere") + + data = self.logger.getRawDataPoints() + + # We should only have the data related to photon 1. + self.assertEqual(4, data.shape[0]) + self.assertTrue(np.array_equal(np.unique(data[:, 4]), [1])) + + def testWhenFilterWithMultipleLabels_shouldKeepPhotonsThatInteractedWithAnyOfTheGivenLabels(self): + # Photon 0 + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=0) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0) + # Photon 1 interacted with "sphere". + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=1) + self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=1) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=1) + # Photon 2 interacted with "otherSphere". + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=2) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=2) + self.logger.logDataPoint(0.4, Vector(3, 5, 0), InteractionKey("otherSphere"), ID=2) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=2) + + self.logger.filter(detectedBy=["sphere", "otherSphere"]) + + data = self.logger.getRawDataPoints() + + # We should only have the data related to photon 1 and photon 2. + self.assertEqual(8, data.shape[0]) + self.assertTrue(np.array_equal(np.unique(data[:, 4]), [1, 2])) + + def testWhenGetFiltered_shouldReturnANewFilteredLogger(self): + # Photon 0 + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=0) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=0) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=0) + # Photon 1 interacted with "sphere". + self.logger.logDataPoint(0.1, Vector(0.7, 0.8, 0.8), InteractionKey("cube"), ID=1) + self.logger.logDataPoint(-0.9, Vector(0.5, 1.0, 0.75), InteractionKey("cube", "cube_top"), ID=1) + self.logger.logDataPoint(0.4, Vector(0, 5, 0), InteractionKey("sphere"), ID=1) + self.logger.logDataPoint(0.2, Vector(0, 0, 0), InteractionKey(WORLD_SOLID_LABEL), ID=1) + + filteredLogger = self.logger.getFiltered(detectedBy="sphere") + + data = filteredLogger.getRawDataPoints() + + # We should only have the data related to photon 1. + self.assertEqual(4, data.shape[0]) + self.assertTrue(np.array_equal(np.unique(data[:, 4]), [1])) + + # Should not modify the original logger. + originalData = self.logger.getRawDataPoints() + self.assertEqual(7, originalData.shape[0]) diff --git a/pytissueoptics/rayscattering/tests/testScatteringScene.py b/pytissueoptics/rayscattering/tests/testScatteringScene.py index e2a9a4e4..739eedbc 100644 --- a/pytissueoptics/rayscattering/tests/testScatteringScene.py +++ b/pytissueoptics/rayscattering/tests/testScatteringScene.py @@ -4,6 +4,7 @@ from mockito import mock, verify, when +from pytissueoptics import Circle from pytissueoptics.rayscattering.materials import ScatteringMaterial from pytissueoptics.rayscattering.scatteringScene import ScatteringScene from pytissueoptics.scene.solids import Cuboid @@ -54,3 +55,12 @@ def testShouldHaveIPPEstimationUsingMeanAlbedoInInfiniteMedium(self): estimation = scene.getEstimatedIPP(weightThreshold) expectedEstimation = -math.log(weightThreshold) / meanAlbedo self.assertAlmostEqual(expectedEstimation, estimation, places=7) + + def testCannotAddFlatSolidThatIsNotADetector(self): + flatSolid = Circle(radius=1) + with self.assertRaises(Exception): + ScatteringScene([flatSolid]) + + def testCanAddFlatDetectors(self): + flatSolid = Circle(radius=1).asDetector() + ScatteringScene([flatSolid]) diff --git a/pytissueoptics/scene/solids/circle.py b/pytissueoptics/scene/solids/circle.py index 7cfba17e..c35b9b44 100644 --- a/pytissueoptics/scene/solids/circle.py +++ b/pytissueoptics/scene/solids/circle.py @@ -35,3 +35,10 @@ def _geometryParams(self) -> dict: "u": self._u, "s": self._s, } + + @property + def isFlat(self) -> bool: + return True + + def __hash__(self): + return hash((self._radius, self._orientation, self._position)) diff --git a/pytissueoptics/scene/solids/rectangle.py b/pytissueoptics/scene/solids/rectangle.py index ad7c753d..736f151d 100644 --- a/pytissueoptics/scene/solids/rectangle.py +++ b/pytissueoptics/scene/solids/rectangle.py @@ -35,3 +35,7 @@ def _computeQuadMesh(self): def _geometryParams(self) -> dict: return {"a": self._a, "b": self._b, "position": self._position, "orientation": self._orientation} + + @property + def isFlat(self) -> bool: + return True diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py index 70408e3f..53874bdc 100644 --- a/pytissueoptics/scene/solids/solid.py +++ b/pytissueoptics/scene/solids/solid.py @@ -380,3 +380,7 @@ def geometryExport(self) -> dict[str, Any]: def _geometryParams(self) -> dict: """To be implemented by Solid subclasses to detail other geometry parameters.""" raise NotImplementedError(f"Geometry parameters not implemented for Solids of type {type(self).__name__}") + + @property + def isFlat(self) -> bool: + return False diff --git a/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py b/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py index 50e86bfe..c502a1be 100644 --- a/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py +++ b/pytissueoptics/scene/tests/intersection/testIntersectionFinder.py @@ -227,6 +227,13 @@ def testGivenSmoothSolid_shouldNotSmoothTheIntersectionNormalIfItChangesTheSignO self.assertIsNotNone(intersection) self.assertNotEqual(smoothNormal, intersection.normal) + def testGivenIgnoreLabel_shouldNotIntersectWithIt(self): + ray = Ray(origin=Vector(0, 0.5, 0), direction=Vector(0, 0, 1)) + solid = Cube(2, position=Vector(0, 0, 5), label="ignoreMe") + + intersection = self.getIntersectionFinder([solid]).findIntersection(ray, WORLD_LABEL, ignoreLabel="ignoreMe") + self.assertIsNone(intersection) + def assertVectorEqual(self, expected, actual): self.assertEqual(expected.x, actual.x) self.assertEqual(expected.y, actual.y) diff --git a/pytissueoptics/scene/tests/solids/testCircle.py b/pytissueoptics/scene/tests/solids/testCircle.py new file mode 100644 index 00000000..c3629e65 --- /dev/null +++ b/pytissueoptics/scene/tests/solids/testCircle.py @@ -0,0 +1,37 @@ +import unittest + +from pytissueoptics.scene.geometry import Vector +from pytissueoptics.scene.solids import Circle + + +class TestCircle(unittest.TestCase): + def testShouldBeFlat(self): + circle = Circle(radius=1) + self.assertTrue(circle.isFlat) + + def testShouldHaveDefaultZOrientation(self): + circle = Circle(radius=1) + polygons = circle.getPolygons() + normals = [p.normal for p in polygons] + for n in normals: + self.assertEqual(n, Vector(0, 0, 1)) + + def testGivenAnOrientation_shouldOrientTheCircle(self): + orientation = Vector(1, 1, 0) + orientation.normalize() + circle = Circle(radius=1, orientation=orientation) + polygons = circle.getPolygons() + normals = [p.normal for p in polygons] + for n in normals: + self.assertEqual(n, orientation) + + def testShouldHaveASingleSurface(self): + circle = Circle(radius=1) + assert len(circle.surfaceLabels) == 1 + + def testShouldHaveDifferentHashWhenOrientationChanges(self): + circle1 = Circle(radius=1, orientation=Vector(0, 0, 1)) + circle2 = Circle(radius=1, orientation=Vector(1, 1, 0)) + self.assertNotEqual(hash(circle1), hash(circle2)) + circle3 = Circle(radius=1, orientation=Vector(0, 0, 1)) + self.assertEqual(hash(circle1), hash(circle3)) diff --git a/pytissueoptics/scene/tests/solids/testRectangle.py b/pytissueoptics/scene/tests/solids/testRectangle.py new file mode 100644 index 00000000..55a6ca11 --- /dev/null +++ b/pytissueoptics/scene/tests/solids/testRectangle.py @@ -0,0 +1,37 @@ +import unittest + +from pytissueoptics.scene.geometry import Vector +from pytissueoptics.scene.solids import Rectangle + + +class TestRectangle(unittest.TestCase): + def testShouldBeFlat(self): + rectangle = Rectangle(2, 2) + self.assertTrue(rectangle.isFlat) + + def testShouldHaveDefaultZOrientation(self): + rectangle = Rectangle(2, 2) + polygons = rectangle.getPolygons() + normals = [p.normal for p in polygons] + for n in normals: + self.assertEqual(n, Vector(0, 0, 1)) + + def testGivenAnOrientation_shouldOrientTheRectangle(self): + orientation = Vector(1, 1, 0) + orientation.normalize() + rectangle = Rectangle(2, 2, orientation=orientation) + polygons = rectangle.getPolygons() + normals = [p.normal for p in polygons] + for n in normals: + self.assertEqual(n, orientation) + + def testShouldHaveASingleSurface(self): + rectangle = Rectangle(2, 2) + assert len(rectangle.surfaceLabels) == 1 + + def testShouldHaveDifferentHashWhenOrientationChanges(self): + rectangle1 = Rectangle(2, 2, orientation=Vector(0, 0, 1)) + rectangle2 = Rectangle(2, 2, orientation=Vector(1, 1, 0)) + self.assertNotEqual(hash(rectangle1), hash(rectangle2)) + rectangle3 = Rectangle(2, 2, orientation=Vector(0, 0, 1)) + self.assertEqual(hash(rectangle1), hash(rectangle3)) diff --git a/pytissueoptics/scene/tests/solids/testSolid.py b/pytissueoptics/scene/tests/solids/testSolid.py index 3dac4eb8..e5329e99 100644 --- a/pytissueoptics/scene/tests/solids/testSolid.py +++ b/pytissueoptics/scene/tests/solids/testSolid.py @@ -257,6 +257,24 @@ def testWhenCheckIfContainsAVertexPartiallyInside_shouldWarnAndReturnFalse(self) with self.assertWarns(RuntimeWarning): self.assertFalse(self.solid.contains(Vertex(2, 2, -0.75))) + def testShouldNotBeFlat(self): + self.assertFalse(self.solid.isFlat) + + def shouldHaveNoDetectorProperties(self): + with self.assertRaises(RuntimeError): + _ = self.solid.detectorAcceptanceCosine + + def testWhenConvertToDetector_shouldSetDetectorCosine(self): + halfAngle = math.pi / 8 + self.solid.asDetector(halfAngle) + expectedAcceptanceCosine = math.cos(halfAngle) + self.assertEqual(expectedAcceptanceCosine, self.solid.detectorAcceptanceCosine) + + def testWhenConvertToDetector_shouldRemoveMaterial(self): + self.assertIsNotNone(self.solid.getEnvironment().material) + self.solid.asDetector() + self.assertIsNone(self.solid.getEnvironment().material) + @staticmethod def createPolygonMock() -> Polygon: polygon = mock(Polygon)