From 5f26732fbcc1c4216f83ee9980f67b867a85794b Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Fri, 8 Nov 2024 21:11:34 -0500
Subject: [PATCH 01/23] use math.random and tweak equality condition to fit
 mcml

---
 pytissueoptics/rayscattering/fresnel.py                  | 2 +-
 .../rayscattering/materials/scatteringMaterial.py        | 9 +++++----
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/pytissueoptics/rayscattering/fresnel.py b/pytissueoptics/rayscattering/fresnel.py
index 76d3f465..572048a0 100644
--- a/pytissueoptics/rayscattering/fresnel.py
+++ b/pytissueoptics/rayscattering/fresnel.py
@@ -56,7 +56,7 @@ def _create(self, nextEnvironment, incidencePlane) -> FresnelIntersection:
 
     def _getIsReflected(self) -> bool:
         R = self._getReflectionCoefficient()
-        if random.random() < R:
+        if random.random() <= R:
             return True
         return False
 
diff --git a/pytissueoptics/rayscattering/materials/scatteringMaterial.py b/pytissueoptics/rayscattering/materials/scatteringMaterial.py
index 8da70d73..6fca0f71 100644
--- a/pytissueoptics/rayscattering/materials/scatteringMaterial.py
+++ b/pytissueoptics/rayscattering/materials/scatteringMaterial.py
@@ -1,4 +1,5 @@
 import math
+import random
 
 import numpy as np
 
@@ -32,16 +33,16 @@ def getScatteringDistance(self):
 
         rnd = 0
         while rnd == 0:
-            rnd = np.random.random()
+            rnd = random.random()
         return -np.log(rnd) / self.mu_t
 
     def getScatteringAngles(self):
-        phi = np.random.random() * 2 * np.pi
+        phi = random.random() * 2 * np.pi
         g = self.g
         if g == 0:
-            cost = 2 * np.random.random() - 1
+            cost = 2 * random.random() - 1
         else:
-            temp = (1 - g * g) / (1 - g + 2 * g * np.random.random())
+            temp = (1 - g * g) / (1 - g + 2 * g * random.random())
             cost = (1 + g * g - temp * temp) / (2 * g)
         return np.arccos(cost), phi
 

From bdace21b3abaec2885019ebfe3f83b19c7453305 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Fri, 8 Nov 2024 21:21:14 -0500
Subject: [PATCH 02/23] switch to MCML spin formula

---
 pytissueoptics/rayscattering/photon.py  |  3 +-
 pytissueoptics/scene/geometry/vector.py | 38 +++++++++++++++++--------
 2 files changed, 27 insertions(+), 14 deletions(-)

diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index f4172c2f..30faf269 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -168,8 +168,7 @@ def scatter(self):
         self.interact()
 
     def scatterBy(self, theta, phi):
-        self._er.rotateAround(self._direction, phi)
-        self._direction.rotateAround(self._er, theta)
+        self._direction.spin(theta, phi)
 
     def interact(self):
         delta = self._weight * self.material.getAlbedo()
diff --git a/pytissueoptics/scene/geometry/vector.py b/pytissueoptics/scene/geometry/vector.py
index a91193dc..df834d7c 100644
--- a/pytissueoptics/scene/geometry/vector.py
+++ b/pytissueoptics/scene/geometry/vector.py
@@ -1,5 +1,7 @@
 import math
 
+COS_ZERO = 1 - 1e-12
+
 
 class Vector:
     """
@@ -96,22 +98,34 @@ def update(self, x: float, y: float, z: float):
     def copy(self) -> 'Vector':
         return Vector(self._x, self._y, self._z)
 
+    def spin(self, theta: float, phi: float):
+        """ In spherical coordinates. Taken directly from MCML code. """
+        cosp = math.cos(phi)
+        if phi < math.pi:
+            sinp = math.sqrt(1 - cosp*cosp)
+        else:
+            sinp = -math.sqrt(1 - cosp*cosp)
+
+        cost = math.cos(theta)
+        sint = math.sqrt(1 - cost*cost)
+
+        if abs(self._z) > COS_ZERO:
+            ux = sint*cosp
+            uy = sint*sinp
+            uz = cost * (1 if self._z >= 0 else -1)
+        else:
+            temp = math.sqrt(1 - self._z*self._z)
+            ux = sint*(self._x*self._z*cosp - self._y*sinp) / temp + self._x*cost
+            uy = sint*(self._y*self._z*cosp + self._x*sinp) / temp + self._y*cost
+            uz = -sint*cosp*temp + self._z*cost
+
+        self.update(ux, uy, uz)
+
     def rotateAround(self, unitAxis: 'Vector', theta: float):
         """
         Rotate the vector around `unitAxis` by `theta` radians. Assumes the axis to be a unit vector.
+        Uses Rodrigues' rotation formula.
         """
-        # This is the most expensive (and most common)
-        # operation when performing Monte Carlo in tissue
-        # (15% of time spent here). It is difficult to optimize without
-        # making it even less readable than it currently is
-        # http://en.wikipedia.org/wiki/Rotation_matrix
-        #
-        # Several options were tried in the past such as
-        # external not-so-portable C library, unreadable
-        # shortcuts, sine and cosine lookup tables, etc...
-        # and the performance gain was minimal (<20%).
-        # For now, this is the best, most readable solution.
-
         cost = math.cos(theta)
         sint = math.sin(theta)
         one_cost = 1 - cost

From 06168513ac9fe217e951c1e51eee6502488f8cc1 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Fri, 8 Nov 2024 21:29:56 -0500
Subject: [PATCH 03/23] opencl: tweak fresnel equality

---
 pytissueoptics/rayscattering/opencl/src/fresnel.c | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/fresnel.c b/pytissueoptics/rayscattering/opencl/src/fresnel.c
index c8021616..1decf6f0 100644
--- a/pytissueoptics/rayscattering/opencl/src/fresnel.c
+++ b/pytissueoptics/rayscattering/opencl/src/fresnel.c
@@ -40,7 +40,7 @@ float _getReflectionCoefficient(float n1, float n2, float thetaIn) {
 bool _getIsReflected(float nIn, float nOut, float thetaIn, __global uint *seeds, uint gid) {
     float R = _getReflectionCoefficient(nIn, nOut, thetaIn);
     float randomFloat = getRandomFloatValue(seeds, gid);
-    if (R > randomFloat) {
+    if (R >= randomFloat) {
         return true;
     }
     return false;

From 9d7bf7e7fc636670a6d2d7aed979337edc029e65 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Fri, 8 Nov 2024 22:31:29 -0500
Subject: [PATCH 04/23] opencl: switch to mcml rotate

---
 .../rayscattering/opencl/src/propagation.c    |  3 +--
 .../opencl/src/vectorOperators.c              | 23 +++++++++++++++++++
 2 files changed, 24 insertions(+), 2 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index 5e0db101..b6d1a9a5 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -12,8 +12,7 @@ void moveBy(float distance, __global Photon *photons, uint photonID){
 }
 
 void scatterBy(float phi, float theta, __global Photon *photons, uint photonID){
-    rotateAroundAxisGlobal(&photons[photonID].er, &photons[photonID].direction, phi);
-    rotateAroundAxisGlobal(&photons[photonID].direction, &photons[photonID].er, theta);
+    spin(&photons[photonID].direction, theta, phi);
 }
 
 void decreaseWeightBy(float delta_weight, __global Photon *photons, uint photonID){
diff --git a/pytissueoptics/rayscattering/opencl/src/vectorOperators.c b/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
index 32dec7c2..efb0358b 100644
--- a/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
+++ b/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
@@ -1,3 +1,4 @@
+__constant float COS_ZERO = 0.999999f;
 
 void normalizeVectorLocal(float3 *vector){
     float length = sqrt(vector->x * vector->x + vector->y * vector->y + vector->z * vector->z);
@@ -67,6 +68,28 @@ void rotateAround(__global float3 *mainVector, float3 *axisVector, float theta){
     mainVector->z = z;
 }
 
+void spin(__global float3 *mainVector, float theta, float phi){
+    float cosp = cos(phi);
+    float sinp = phi < M_PI ? sqrt(1 - cosp*cosp) : -sqrt(1 - cosp*cosp);
+    float cost = cos(theta);
+    float sint = sqrt(1 - cost*cost);
+    float ux, uy, uz;
+    if (fabs(mainVector->z) > COS_ZERO){
+        ux = sint*cosp;
+        uy = sint*sinp;
+        uz = cost * (mainVector->z >= 0 ? 1 : -1);
+    }
+    else{
+        float temp = sqrt(1 - mainVector->z*mainVector->z);
+        ux = sint*(mainVector->x*mainVector->z*cosp - mainVector->y*sinp) / temp + mainVector->x*cost;
+        uy = sint*(mainVector->y*mainVector->z*cosp + mainVector->x*sinp) / temp + mainVector->y*cost;
+        uz = -sint*cosp*temp + mainVector->z*cost;
+    }
+    mainVector->x = ux;
+    mainVector->y = uy;
+    mainVector->z = uz;
+}
+
 float3 getAnyOrthogonalGlobal(__global float3 *vector){
     if (fabs(vector->z) < fabs(vector->x)){
         return (float3)(vector->y, -vector->x, 0.0f);

From 3739521ac9ea83e3b27cf269351f65af7d2f7496 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Fri, 8 Nov 2024 22:40:48 -0500
Subject: [PATCH 05/23] implement (CPU) new float limit intersections

might not work with fast intersector for now (which is not used in opencl)
---
 pytissueoptics/rayscattering/photon.py        | 23 ++-----------
 .../rayscattering/tests/testPhoton.py         |  5 ++-
 .../scene/intersection/intersectionFinder.py  | 34 +++++++++++++------
 .../intersection/mollerTrumboreIntersect.py   | 33 +++++++++++-------
 .../intersection/testPolygonIntersect.py      |  6 ++--
 5 files changed, 51 insertions(+), 50 deletions(-)

diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index 30faf269..0b4ba59b 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -7,7 +7,6 @@
 from pytissueoptics.scene.geometry import Environment, Vector
 from pytissueoptics.scene.intersection import Ray
 from pytissueoptics.scene.intersection.intersectionFinder import IntersectionFinder, Intersection
-from pytissueoptics.scene.intersection.mollerTrumboreIntersect import EPS_CORRECTION
 from pytissueoptics.scene.logger import Logger, InteractionKey
 
 WORLD_LABEL = "world"
@@ -78,7 +77,7 @@ def step(self, distance=0) -> float:
 
         intersection = self._getIntersection(distance)
 
-        if intersection and not intersection.isTooClose:
+        if intersection:
             self.moveBy(intersection.distance)
             distanceLeft = self.reflectOrRefract(intersection)
         else:
@@ -89,15 +88,6 @@ def step(self, distance=0) -> float:
             self.moveBy(distance)
             distanceLeft = 0
 
-            if intersection and intersection.isTooClose:
-                # Photon will land too close to the surface, so we need to move it away from the surface.
-                stepSign = 1
-                solidTowardsNormal = intersection.outsideEnvironment.solid
-                if solidTowardsNormal != self._environment.solid:
-                    stepSign = -1
-                stepCorrection = intersection.normal * stepSign * EPS_CORRECTION
-                self._position += stepCorrection
-
             self.scatter()
 
         return distanceLeft
@@ -107,7 +97,7 @@ def _getIntersection(self, distance) -> Optional[Intersection]:
             return None
 
         stepRay = Ray(self._position, self._direction, distance)
-        return self._intersectionFinder.findIntersection(stepRay)
+        return self._intersectionFinder.findIntersection(stepRay, self.solidLabel)
 
     def reflectOrRefract(self, intersection: Intersection):
         fresnelIntersection = self._getFresnelIntersection(intersection)
@@ -137,15 +127,6 @@ def reflectOrRefract(self, intersection: Intersection):
 
             self._environment = fresnelIntersection.nextEnvironment
 
-        # Move away from intersecting surface by a small amount
-        stepCorrection = intersection.normal * stepSign * EPS_CORRECTION
-        self._position += stepCorrection
-
-        # Remove this distance correction from the distance left, but set to zero if the result is negative.
-        intersection.distanceLeft -= EPS_CORRECTION
-        if intersection.distanceLeft < 0:
-            intersection.distanceLeft = 0
-
         return intersection.distanceLeft
 
     def _getFresnelIntersection(self, intersection: Intersection) -> FresnelIntersection:
diff --git a/pytissueoptics/rayscattering/tests/testPhoton.py b/pytissueoptics/rayscattering/tests/testPhoton.py
index 1e91e8b6..c5821215 100644
--- a/pytissueoptics/rayscattering/tests/testPhoton.py
+++ b/pytissueoptics/rayscattering/tests/testPhoton.py
@@ -17,7 +17,7 @@
 from pytissueoptics.scene.solids import Solid
 
 
-EPS = MollerTrumboreIntersect.EPS
+EPS = MollerTrumboreIntersect.EPS_CATCH
 
 
 class TestPhoton(unittest.TestCase):
@@ -230,7 +230,7 @@ def testWhenStepWithReflectingIntersection_shouldReturnDistanceLeft(self):
 
         distanceLeft = self.photon.step(totalDistance)
 
-        expectedDistanceLeft = totalDistance - intersectionDistance - EPS_CORRECTION
+        expectedDistanceLeft = totalDistance - intersectionDistance
         self.assertAlmostEqual(expectedDistanceLeft, distanceLeft)
 
     def testWhenStepWithRefractingIntersection_shouldUpdatePhotonMaterialToNextMaterial(self):
@@ -267,7 +267,6 @@ def testWhenStepWithRefractingIntersection_shouldReturnDistanceLeftInNextMateria
         distanceLeft = self.photon.step(scatteringDistance)
 
         expectedDistanceLeft = (scatteringDistance - intersectionDistance) * material.mu_t / nextMaterial.mu_t
-        expectedDistanceLeft -= EPS_CORRECTION
         self.assertAlmostEqual(expectedDistanceLeft, distanceLeft)
 
     def testWhenStepWithRefractingIntersectionToVacuum_shouldReturnInfiniteDistanceLeft(self):
diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py
index 6705bd02..b89b623f 100644
--- a/pytissueoptics/scene/intersection/intersectionFinder.py
+++ b/pytissueoptics/scene/intersection/intersectionFinder.py
@@ -32,19 +32,28 @@ def __init__(self, scene: Scene):
         self._polygonIntersect = MollerTrumboreIntersect()
         self._boxIntersect = GemsBoxIntersect()
 
-    def findIntersection(self, ray: Ray) -> Optional[Intersection]:
+    def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]:
         raise NotImplementedError
 
-    def _findClosestPolygonIntersection(self, ray: Ray, polygons: List[Polygon]) -> Optional[Intersection]:
+    def _findClosestPolygonIntersection(
+        self, ray: Ray, polygons: List[Polygon], currentSolidLabel: str
+    ) -> Optional[Intersection]:
         closestIntersection = Intersection(sys.maxsize)
         for polygon in polygons:
+            # Skip intersection test if the ray is heading towards its current solid
+            # (possible because of the epsilon catch zone in our Moller Trumbore intersect).
+            isGoingInside = ray.direction.dot(polygon.normal) < 0
+            nextSolid = polygon.insideEnvironment.solid if isGoingInside else polygon.outsideEnvironment.solid
+            nextLabel = 'world' if nextSolid is None else nextSolid.getLabel()
+            if nextLabel == currentSolidLabel:
+                continue
+
             intersectionPoint = self._polygonIntersect.getIntersection(ray, polygon)
             if not intersectionPoint:
                 continue
             distance = (intersectionPoint - ray.origin).getNorm()
             if distance < closestIntersection.distance:
                 closestIntersection = Intersection(distance, intersectionPoint, polygon)
-                closestIntersection.isTooClose = ray.isTooClose
 
         if closestIntersection.distance == sys.maxsize:
             return None
@@ -73,7 +82,7 @@ def _composeIntersection(ray: Ray, intersection: Intersection) -> Optional[Inter
 
 
 class SimpleIntersectionFinder(IntersectionFinder):
-    def findIntersection(self, ray: Ray) -> Optional[Intersection]:
+    def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]:
         """
         Find the closest intersection between a ray and the scene.
 
@@ -88,7 +97,7 @@ def findIntersection(self, ray: Ray) -> Optional[Intersection]:
             meaningful way, so in this case we need to test all of them (candidate distance is zero) before
             returning the closest intersection found.
         """
-        bboxIntersections = self._findBBoxIntersectingSolids(ray)
+        bboxIntersections = self._findBBoxIntersectingSolids(ray, currentSolidLabel)
         bboxIntersections.sort(key=lambda x: x[0])
 
         closestDistance = sys.maxsize
@@ -97,20 +106,23 @@ def findIntersection(self, ray: Ray) -> Optional[Intersection]:
             contained = distance == 0
             if not contained and closestIntersection:
                 break
-            intersection = self._findClosestPolygonIntersection(ray, solid.getPolygons())
+            intersection = self._findClosestPolygonIntersection(ray, solid.getPolygons(), currentSolidLabel)
             if intersection and intersection.distance < closestDistance:
                 closestDistance = intersection.distance
                 closestIntersection = intersection
 
         return self._composeIntersection(ray, closestIntersection)
 
-    def _findBBoxIntersectingSolids(self, ray) -> Optional[List[Tuple[float, Solid]]]:
+    def _findBBoxIntersectingSolids(self, ray: Ray, currentSolidLabel: str) -> 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:
-            bboxIntersectionPoint = self._boxIntersect.getIntersection(ray, solid.bbox)
+            if solid.getLabel() == currentSolidLabel:
+                bboxIntersectionPoint = ray.origin
+            else:
+                bboxIntersectionPoint = self._boxIntersect.getIntersection(ray, solid.bbox)
             if not bboxIntersectionPoint:
                 continue
             distance = (bboxIntersectionPoint - ray.origin).getNorm()
@@ -124,13 +136,15 @@ def __init__(self, scene: Scene, constructor=NoSplitThreeAxesConstructor(), maxD
         self._partition = SpacePartition(self._scene.getBoundingBox(), self._scene.getPolygons(), constructor,
                                          maxDepth, minLeafSize)
 
-    def findIntersection(self, ray: Ray) -> Optional[Intersection]:
+    def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Intersection]:
+        self._currentSolidLabel = currentSolidLabel
         intersection = self._findIntersection(ray, self._partition.root)
         return self._composeIntersection(ray, intersection)
 
     def _findIntersection(self, ray: Ray, node: Node, closestDistance=sys.maxsize) -> Optional[Intersection]:
+        # TODO: implement a way to test triangles that are close behind the ray origin.
         if node.isLeaf:
-            intersection = self._findClosestPolygonIntersection(ray, node.polygons)
+            intersection = self._findClosestPolygonIntersection(ray, node.polygons, self._currentSolidLabel)
             return intersection
 
         if not self._nodeIsWorthExploring(ray, node, closestDistance):
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index 999f619f..40cd34e3 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -8,7 +8,7 @@
 
 
 class MollerTrumboreIntersect:
-    EPS = 0.00001
+    EPS_CATCH = 0.000001
     EPS_PARALLEL = 0.00001
     EPS_SIDE = 0.000001
 
@@ -44,6 +44,8 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
         tVector = ray.origin - v1
         u = tVector.dot(pVector) * inverseDeterminant
         if u < -self.EPS_SIDE or u > 1.:
+            # EPS_SIDE is used to make the triangle a bit larger than it is
+            # to be sure a ray could not sneak between two triangles.
             return None
 
         qVector = tVector.cross(edgeA)
@@ -51,21 +53,26 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
         if v < -self.EPS_SIDE or u + v > 1.:
             return None
 
+        # Distance to intersection point
         t = edgeB.dot(qVector) * inverseDeterminant
-        if t < 0.:
-            return None
-
-        if ray.length is None:
+        if t > 0 and (ray.length > t or ray.length is None):
+            # Case 1: Trivial case. Intersects.
             return ray.origin + ray.direction * t
 
-        if t > (ray.length + self.EPS):
-            # No intersection, it's too far away
-            return None
-        elif t > ray.length:
-            # Just a bit too far away. There is no intersection, but we cannot accept ray to land here.
-            ray.isTooClose = True
-
-        return ray.origin + ray.direction * t
+        # Next we need to check if the intersection is inside the epsilon catch zone (forward or backward).
+        # Note that this mechanic only works when same-solid intersections are ignored before calling this function.
+        dt = t - ray.length
+        dt_T = abs(triangle.normal.dot(ray.direction) * dt)
+        if t > ray.length and dt_T < self.EPS_CATCH:
+            # Case 2: Forward epsilon catch. Ray ends close to the triangle, so we intersect at the ray's end.
+            return ray.origin + ray.direction * ray.length
+        if t < 0 and dt_T < self.EPS_CATCH:
+            # Case 3: Backward epsilon catch. Ray starts close to the triangle, so we intersect at the origin.
+            # This requires the intersector to always test triangles (or at least, close ones) of the origin solid.
+            return ray.origin
+
+        # Case 4: No intersection.
+        return None
 
     def _getQuadIntersection(self, ray: Ray, quad: Quad) -> Optional[Vector]:
         v1, v2, v3, v4 = quad.vertices
diff --git a/pytissueoptics/scene/tests/intersection/testPolygonIntersect.py b/pytissueoptics/scene/tests/intersection/testPolygonIntersect.py
index 129a6edc..cfe2de5c 100644
--- a/pytissueoptics/scene/tests/intersection/testPolygonIntersect.py
+++ b/pytissueoptics/scene/tests/intersection/testPolygonIntersect.py
@@ -75,7 +75,7 @@ def testGivenRayLandsInEpsilonRegionBeforePolygon_shouldStillReturnIntersectionP
         rayOrigin = Vector(0.25, 0.25, 2)
         rayDirection = Vector(0, 0, -1)
         rayDirection.normalize()
-        ray = Ray(rayOrigin, rayDirection, length=2 - MollerTrumboreIntersect.EPS / 2)
+        ray = Ray(rayOrigin, rayDirection, length=2 - MollerTrumboreIntersect.EPS_CATCH / 2)
 
         for poly in ([self.triangle, self.quad, self.polygon]):
             intersection = self.intersectStrategy.getIntersection(ray, poly)
@@ -89,7 +89,7 @@ def testGivenRayLandsInEpsilonRegionBeforePolygon_shouldLabelTheRayAsTooCloseFor
         rayOrigin = Vector(0.25, 0.25, 2)
         rayDirection = Vector(0, 0, -1)
         rayDirection.normalize()
-        ray = Ray(rayOrigin, rayDirection, length=2 - MollerTrumboreIntersect.EPS * 0.9)
+        ray = Ray(rayOrigin, rayDirection, length=2 - MollerTrumboreIntersect.EPS_CATCH * 0.9)
 
         for poly in ([self.triangle, self.quad, self.polygon]):
             _ = self.intersectStrategy.getIntersection(ray, poly)
@@ -99,7 +99,7 @@ def testGivenRayLandsBeforeEpsilonRegionOfPolygon_shouldReturnNone(self):
         rayOrigin = Vector(0.25, 0.25, 2)
         rayDirection = Vector(0, 0, -1)
         rayDirection.normalize()
-        ray = Ray(rayOrigin, rayDirection, length=2 - MollerTrumboreIntersect.EPS * 1.1)
+        ray = Ray(rayOrigin, rayDirection, length=2 - MollerTrumboreIntersect.EPS_CATCH * 1.1)
 
         for poly in ([self.triangle, self.quad, self.polygon]):
             intersection = self.intersectStrategy.getIntersection(ray, poly)

From 4c77a909bc25cf12bbbdd096150c1d3ae491cbf7 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Fri, 8 Nov 2024 23:36:41 -0500
Subject: [PATCH 06/23] opencl: new float limit intersections

---
 .../rayscattering/opencl/src/intersection.c   | 56 ++++++++++++-------
 .../rayscattering/opencl/src/propagation.c    | 22 +-------
 2 files changed, 38 insertions(+), 40 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index ffeaa803..1659e9b1 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -1,5 +1,4 @@
-__constant float EPS = 0.00001f;
-__constant float EPS_CORRECTION = 0.0005f;
+__constant float EPS_CATCH = 0.000001f;
 __constant float EPS_PARALLEL = 0.00001f;
 __constant float EPS_SIDE = 0.000001f;
 
@@ -113,11 +112,17 @@ GemsBoxIntersection _getBBoxIntersection(Ray ray, float3 minCornerVector, float3
     return intersection;
 }
 
-void _findBBoxIntersectingSolids(Ray ray, Scene *scene, uint gid){
+void _findBBoxIntersectingSolids(Ray ray, Scene *scene, uint gid, uint photonSolidID) {
 
     for (uint i = 0; i < scene->nSolids; i++) {
         uint boxGID = gid * scene->nSolids + i;
-        scene->solidCandidates[boxGID].solidID = i + 1;
+        uint solidID = i + 1;
+        scene->solidCandidates[boxGID].solidID = solidID;
+
+        if (solidID == photonSolidID) {
+            scene->solidCandidates[boxGID].distance = 0;
+            continue;
+        }
 
         GemsBoxIntersection gemsIntersection = _getBBoxIntersection(ray, scene->solids[i].bbox_min, scene->solids[i].bbox_max);
         if (gemsIntersection.rayIsInside) {
@@ -155,7 +160,7 @@ struct HitPoint {
 
 typedef struct HitPoint HitPoint;
 
-HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3) {
+HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, float3 normal) {
     HitPoint hitPoint;
     hitPoint.exists = false;
     hitPoint.isTooClose = false;
@@ -185,35 +190,46 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3) {
 
     float t = dot(edgeB, qVector) * invDet;
 
-    if (t < 0.0f) {
+    if (t > 0 && ray.length > t){
+        hitPoint.exists = true;
+        hitPoint.position = ray.origin + t * ray.direction;
         return hitPoint;
     }
 
-    if (t > (ray.length + EPS)) {
-        // No Intersection, it's too far away
+    float dt = t - ray.length;
+    float dt_T = fabs(dot(normal, ray.direction) * dt);
+    if (t > ray.length && dt_T < EPS_CATCH) {
+        hitPoint.exists = true;
+        hitPoint.position = ray.origin + ray.length * ray.direction;
+        return hitPoint;
+    }
+    if (t < 0 && dt_T < EPS_CATCH) {
+        hitPoint.exists = true;
+        hitPoint.position = ray.origin;
         return hitPoint;
-    } else if (t > ray.length) {
-        // Just a bit too far away. There is no intersection, but we cannot accept photon to land here.
-        // hitPoint will also be returned with exists = true in order to consider this event and possibly process it if it was the closest "hit".
-        hitPoint.isTooClose = true;
     }
 
-    hitPoint.exists = true;
-    hitPoint.position = ray.origin + t * ray.direction;
     return hitPoint;
 }
 
 Intersection _findClosestPolygonIntersection(Ray ray, uint solidID,
                                             __global Solid *solids, __global Surface *surfaces,
-                                            __global Triangle *triangles, __global Vertex *vertices) {
+                                            __global Triangle *triangles, __global Vertex *vertices,
+                                            uint photonSolidID) {
     Intersection intersection;
     intersection.exists = false;
     intersection.isTooClose = false;
     intersection.distance = INFINITY;
     for (uint s = solids[solidID-1].firstSurfaceID; s <= solids[solidID-1].lastSurfaceID; s++) {
         for (uint p = surfaces[s].firstPolygonID; p <= surfaces[s].lastPolygonID; p++) {
+            bool isGoingInside = dot(ray.direction, triangles[p].normal) < 0;
+            uint nextSolidID = isGoingInside ? surfaces[s].insideSolidID : surfaces[s].outsideSolidID;
+            if (nextSolidID == photonSolidID) {
+                continue;
+            }
+
             uint vertexIDs[3] = {triangles[p].vertexIDs[0], triangles[p].vertexIDs[1], triangles[p].vertexIDs[2]};
-            HitPoint hitPoint = _getTriangleIntersection(ray, vertices[vertexIDs[0]].position, vertices[vertexIDs[1]].position, vertices[vertexIDs[2]].position);
+            HitPoint hitPoint = _getTriangleIntersection(ray, vertices[vertexIDs[0]].position, vertices[vertexIDs[1]].position, vertices[vertexIDs[2]].position, triangles[p].normal);
             if (!hitPoint.exists) {
                 continue;
             }
@@ -301,12 +317,12 @@ void _composeIntersection(Intersection *intersection, Ray *ray, Scene *scene) {
     intersection->distanceLeft = ray->length - intersection->distance;
 }
 
-Intersection findIntersection(Ray ray, Scene *scene, uint gid) {
+Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidID) {
     /*
     OpenCL implementation of the Python module SimpleIntersectionFinder
     See the Python module documentation for more details.
     */
-    _findBBoxIntersectingSolids(ray, scene, gid);
+    _findBBoxIntersectingSolids(ray, scene, gid, photonSolidID);
     _sortSolidCandidates(scene, gid);
 
     Intersection closestIntersection;
@@ -329,7 +345,7 @@ Intersection findIntersection(Ray ray, Scene *scene, uint gid) {
         }
 
         uint solidID = scene->solidCandidates[boxGID].solidID;
-        Intersection intersection = _findClosestPolygonIntersection(ray, solidID, scene->solids, scene->surfaces, scene->triangles, scene->vertices);
+        Intersection intersection = _findClosestPolygonIntersection(ray, solidID, scene->solids, scene->surfaces, scene->triangles, scene->vertices, photonSolidID);
         if (intersection.exists  && intersection.distance < closestIntersection.distance) {
             closestIntersection = intersection;
         }
@@ -345,7 +361,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);
+    intersections[gid] = findIntersection(rays[gid], &scene, gid, 0);
 }
 
 
diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index b6d1a9a5..b40bf1c2 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -125,14 +125,6 @@ float reflectOrRefract(Intersection *intersection, __global Photon *photons, __c
         photons[photonID].solidID = fresnelIntersection.nextSolidID;
     }
 
-    float3 stepCorrection = stepSign * intersection->normal * EPS_CORRECTION;
-    photons[photonID].position += stepCorrection;
-
-    intersection->distanceLeft -= EPS_CORRECTION;
-    if (intersection->distanceLeft < 0) {
-        intersection->distanceLeft = 0;
-    }
-
     return intersection->distanceLeft;
 }
 
@@ -146,11 +138,11 @@ float propagateStep(float distance, __global Photon *photons, __constant Materia
     }
 
     Ray stepRay = {photons[photonID].position, photons[photonID].direction, distance};
-    Intersection intersection = findIntersection(stepRay, scene, gid);
+    Intersection intersection = findIntersection(stepRay, scene, gid, photons[photonID].solidID);
 
     float distanceLeft = 0;
 
-    if (intersection.exists && !intersection.isTooClose){
+    if (intersection.exists){
         moveBy(intersection.distance, photons, photonID);
         distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID);
     } else {
@@ -161,16 +153,6 @@ float propagateStep(float distance, __global Photon *photons, __constant Materia
 
         moveBy(distance, photons, photonID);
 
-        if (intersection.isTooClose){
-            int stepSign = 1;
-            int solidIDTowardsNormal = scene->surfaces[intersection.surfaceID].outsideSolidID;
-            if (solidIDTowardsNormal != photons[photonID].solidID) {
-                stepSign = -1;
-            }
-            float3 stepCorrection = stepSign * intersection.normal * EPS_CORRECTION;
-            photons[photonID].position += stepCorrection;
-        }
-
         scatter(photons, materials, seeds, logger, logIndex, gid, photonID);
     }
 

From 60ddc84fb06a0f28b75d2c0da4bd0e4d2581133a Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sun, 10 Nov 2024 12:55:34 -0500
Subject: [PATCH 07/23] fix spin: normalize direction

---
 pytissueoptics/rayscattering/opencl/src/intersection.c       | 4 ++--
 pytissueoptics/rayscattering/opencl/src/vectorOperators.c    | 1 +
 pytissueoptics/scene/geometry/vector.py                      | 1 +
 pytissueoptics/scene/intersection/mollerTrumboreIntersect.py | 2 +-
 4 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 1659e9b1..13cc7ce8 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -190,7 +190,7 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
 
     float t = dot(edgeB, qVector) * invDet;
 
-    if (t > 0 && ray.length > t){
+    if (t > 0 && ray.length >= t){
         hitPoint.exists = true;
         hitPoint.position = ray.origin + t * ray.direction;
         return hitPoint;
@@ -346,7 +346,7 @@ Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidI
 
         uint solidID = scene->solidCandidates[boxGID].solidID;
         Intersection intersection = _findClosestPolygonIntersection(ray, solidID, scene->solids, scene->surfaces, scene->triangles, scene->vertices, photonSolidID);
-        if (intersection.exists  && intersection.distance < closestIntersection.distance) {
+        if (intersection.exists && intersection.distance < closestIntersection.distance) {
             closestIntersection = intersection;
         }
     }
diff --git a/pytissueoptics/rayscattering/opencl/src/vectorOperators.c b/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
index efb0358b..554dc14c 100644
--- a/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
+++ b/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
@@ -88,6 +88,7 @@ void spin(__global float3 *mainVector, float theta, float phi){
     mainVector->x = ux;
     mainVector->y = uy;
     mainVector->z = uz;
+    normalizeVectorGlobal(mainVector);
 }
 
 float3 getAnyOrthogonalGlobal(__global float3 *vector){
diff --git a/pytissueoptics/scene/geometry/vector.py b/pytissueoptics/scene/geometry/vector.py
index df834d7c..41dfaa52 100644
--- a/pytissueoptics/scene/geometry/vector.py
+++ b/pytissueoptics/scene/geometry/vector.py
@@ -120,6 +120,7 @@ def spin(self, theta: float, phi: float):
             uz = -sint*cosp*temp + self._z*cost
 
         self.update(ux, uy, uz)
+        self.normalize()
 
     def rotateAround(self, unitAxis: 'Vector', theta: float):
         """
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index 40cd34e3..ac02bba7 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -55,7 +55,7 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
 
         # Distance to intersection point
         t = edgeB.dot(qVector) * inverseDeterminant
-        if t > 0 and (ray.length > t or ray.length is None):
+        if t > 0 and (ray.length >= t or ray.length is None):
             # Case 1: Trivial case. Intersects.
             return ray.origin + ray.direction * t
 

From 842b1f3eeaacf20fd880569c1e270bf47b96e461 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sun, 17 Nov 2024 10:51:22 -0500
Subject: [PATCH 08/23] fix triangle intersect dt sign when hitpoint behind
 origin and store distance

---
 .../rayscattering/opencl/src/intersection.c   | 23 +++++++++++--------
 .../intersection/mollerTrumboreIntersect.py   | 11 +++++----
 2 files changed, 21 insertions(+), 13 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 13cc7ce8..e2c022f5 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -1,10 +1,9 @@
-__constant float EPS_CATCH = 0.000001f;
-__constant float EPS_PARALLEL = 0.00001f;
+__constant float EPS_CATCH = 0.00001f;
+__constant float EPS_PARALLEL = 0.000001f;
 __constant float EPS_SIDE = 0.000001f;
 
 struct Intersection {
     uint exists;
-    uint isTooClose;
     float distance;
     float3 position;
     float3 normal;
@@ -154,7 +153,7 @@ void _sortSolidCandidates(Scene *scene, uint gid) {
 
 struct HitPoint {
     bool exists;
-    bool isTooClose;
+    float distance;
     float3 position;
 };
 
@@ -192,19 +191,27 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
 
     if (t > 0 && ray.length >= t){
         hitPoint.exists = true;
+        hitPoint.distance = t;
         hitPoint.position = ray.origin + t * ray.direction;
         return hitPoint;
     }
 
-    float dt = t - ray.length;
+    float dt;
+    if (t <= 0) {
+        dt = t;
+    } else {
+        dt = t - ray.length;
+    }
     float dt_T = fabs(dot(normal, ray.direction) * dt);
     if (t > ray.length && dt_T < EPS_CATCH) {
         hitPoint.exists = true;
+        hitPoint.distance = ray.length;
         hitPoint.position = ray.origin + ray.length * ray.direction;
         return hitPoint;
     }
     if (t < 0 && dt_T < EPS_CATCH) {
         hitPoint.exists = true;
+        hitPoint.distance = 0;
         hitPoint.position = ray.origin;
         return hitPoint;
     }
@@ -233,11 +240,9 @@ Intersection _findClosestPolygonIntersection(Ray ray, uint solidID,
             if (!hitPoint.exists) {
                 continue;
             }
-            float distance = length(hitPoint.position - ray.origin);
-            if (distance < intersection.distance) {
+            if (hitPoint.distance < intersection.distance) {
                 intersection.exists = true;
-                intersection.isTooClose = hitPoint.isTooClose;
-                intersection.distance = distance;
+                intersection.distance = hitPoint.distance;
                 intersection.position = hitPoint.position;
                 intersection.normal = triangles[p].normal;
                 intersection.surfaceID = s;
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index ac02bba7..82cdcd0a 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -4,12 +4,12 @@
 from pytissueoptics.scene.intersection import Ray
 
 
-EPS_CORRECTION = 0.0005
+EPS_CORRECTION = 0.0  # TODO: remove
 
 
 class MollerTrumboreIntersect:
-    EPS_CATCH = 0.000001
-    EPS_PARALLEL = 0.00001
+    EPS_CATCH = 0.00001
+    EPS_PARALLEL = 0.000001
     EPS_SIDE = 0.000001
 
     def getIntersection(self, ray: Ray, polygon: Union[Triangle, Quad, Polygon]) -> Optional[Vector]:
@@ -61,7 +61,10 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
 
         # Next we need to check if the intersection is inside the epsilon catch zone (forward or backward).
         # Note that this mechanic only works when same-solid intersections are ignored before calling this function.
-        dt = t - ray.length
+        if t <= 0:
+            dt = t
+        else:
+            dt = t - ray.length
         dt_T = abs(triangle.normal.dot(ray.direction) * dt)
         if t > ray.length and dt_T < self.EPS_CATCH:
             # Case 2: Forward epsilon catch. Ray ends close to the triangle, so we intersect at the ray's end.

From c7dc7a699b222001f5f398faebc404883edcd43e Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sun, 17 Nov 2024 11:07:43 -0500
Subject: [PATCH 09/23] fix smooth fresnel edge case: limit deflection angles

---
 .../rayscattering/opencl/src/intersection.c   | 10 ++++++--
 .../rayscattering/opencl/src/propagation.c    | 24 ++++++++++++------
 pytissueoptics/rayscattering/photon.py        | 25 +++++++++++++------
 .../scene/intersection/intersectionFinder.py  | 21 ++++++++++------
 4 files changed, 54 insertions(+), 26 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index e2c022f5..f0e8a3ed 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -10,6 +10,8 @@ struct Intersection {
     uint surfaceID;
     uint polygonID;
     float distanceLeft;
+    bool isSmooth;
+    float3 rawNormal;
 };
 
 typedef struct Intersection Intersection;
@@ -305,10 +307,13 @@ void setSmoothNormal(Intersection *intersection, __global Triangle *triangles, _
     // Not accounting for this can lead to a photon slightly going inside another solid mesh, but being considered as leaving the other solid (during FresnelIntersection calculations).
     // Which would result in the wrong next environment being set as well as the wrong step correction being applied after refraction.
     if (dot(newNormal, ray->direction) * dot(intersection->normal, ray->direction) < 0) {
+        intersection->isSmooth = false;
         return;
     }
-    intersection->normal = newNormal;
-    intersection->normal = normalize(intersection->normal);
+    intersection->normal = normalize(newNormal);
+
+    intersection->isSmooth = true;
+    intersection->rawNormal = triangles[intersection->polygonID].normal;
 }
 
 void _composeIntersection(Intersection *intersection, Ray *ray, Scene *scene) {
@@ -316,6 +321,7 @@ void _composeIntersection(Intersection *intersection, Ray *ray, Scene *scene) {
         return;
     }
 
+    intersection->isSmooth = false;
     if (scene->surfaces[intersection->surfaceID].toSmooth) {
         setSmoothNormal(intersection, scene->triangles, scene->vertices, ray);
     }
diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index b40bf1c2..a3845941 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -6,6 +6,7 @@
 
 __constant int NO_SOLID_ID = -1;
 __constant int NO_SURFACE_ID = -1;
+__constant float MIN_ANGLE = 0.0001f;
 
 void moveBy(float distance, __global Photon *photons, uint photonID){
     photons[photonID].position += (distance * photons[photonID].direction);
@@ -96,20 +97,27 @@ float reflectOrRefract(Intersection *intersection, __global Photon *photons, __c
         __global Surface *surfaces, __global DataPoint *logger, uint *logIndex, __global uint *seeds, uint gid, uint photonID){
     FresnelIntersection fresnelIntersection = computeFresnelIntersection(photons[photonID].direction, intersection,
                                                                          materials, surfaces, seeds, gid);
-    int stepSign = 1;
-    int solidIDTowardsNormal = surfaces[intersection->surfaceID].outsideSolidID;
-    if (solidIDTowardsNormal != photons[photonID].solidID) {
-        stepSign = -1;
-    }
-    if (!fresnelIntersection.isReflected) {
-        stepSign *= -1;
-    }
 
     if (fresnelIntersection.isReflected) {
+        if (intersection->isSmooth) {
+            // Prevent reflection from crossing the raw surface.
+            float smoothAngle = acos(dot(intersection->normal, intersection->rawNormal));
+            float minDeflectionAngle = smoothAngle + fabs(fresnelIntersection.angleDeflection) / 2 + MIN_ANGLE;
+            if (fabs(fresnelIntersection.angleDeflection) < minDeflectionAngle) {
+                fresnelIntersection.angleDeflection = sign(fresnelIntersection.angleDeflection) * minDeflectionAngle;
+            }
+        }
         reflect(&fresnelIntersection, photons, photonID);
     }
     else {
         logIntersection(intersection, photons, surfaces, logger, logIndex, photonID);
+        if (intersection->isSmooth) {
+            // Prevent refraction from not crossing the raw surface.
+            float maxDeflectionAngle = fabs(M_PI_F / 2 - acos(dot(intersection->rawNormal, photons[photonID].direction))) - MIN_ANGLE;
+            if (fabs(fresnelIntersection.angleDeflection) > maxDeflectionAngle) {
+                fresnelIntersection.angleDeflection = sign(fresnelIntersection.angleDeflection) * maxDeflectionAngle;
+            }
+        }
         refract(&fresnelIntersection, photons, photonID);
 
         float mut1 = materials[photons[photonID].materialID].mu_t;
diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index 0b4ba59b..07438f9f 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -2,6 +2,8 @@
 import random
 from typing import Optional
 
+import numpy as np
+
 from pytissueoptics.rayscattering.fresnel import FresnelIntersect, FresnelIntersection
 from pytissueoptics.rayscattering.materials import ScatteringMaterial
 from pytissueoptics.scene.geometry import Environment, Vector
@@ -11,6 +13,7 @@
 
 WORLD_LABEL = "world"
 WEIGHT_THRESHOLD = 1e-4
+MIN_ANGLE = 0.0001
 
 
 class Photon:
@@ -102,18 +105,24 @@ def _getIntersection(self, distance) -> Optional[Intersection]:
     def reflectOrRefract(self, intersection: Intersection):
         fresnelIntersection = self._getFresnelIntersection(intersection)
 
-        # Determine required step sign to move away from intersecting surface
-        stepSign = 1
-        solidTowardsNormal = intersection.outsideEnvironment.solid
-        if solidTowardsNormal != self._environment.solid:
-            stepSign = -1
-        if not fresnelIntersection.isReflected:
-            stepSign *= -1
-
         if fresnelIntersection.isReflected:
+            if intersection.isSmooth:
+                # Prevent reflection from crossing the raw surface.
+                smoothAngle = math.acos(intersection.normal.dot(intersection.polygon.normal))
+                minDeflectionAngle = smoothAngle + abs(fresnelIntersection.angleDeflection) / 2 + MIN_ANGLE
+                if abs(fresnelIntersection.angleDeflection) < minDeflectionAngle:
+                    fresnelIntersection.angleDeflection = minDeflectionAngle * np.sign(fresnelIntersection.angleDeflection)
+
             self.reflect(fresnelIntersection)
         else:
             self._logIntersection(intersection)
+
+            if intersection.isSmooth:
+                # Prevent refraction from not crossing the raw surface.
+                maxDeflectionAngle = abs(np.pi / 2 - math.acos(intersection.polygon.normal.dot(self._direction))) - MIN_ANGLE
+                if abs(fresnelIntersection.angleDeflection) > maxDeflectionAngle:
+                    fresnelIntersection.angleDeflection = maxDeflectionAngle * np.sign(fresnelIntersection.angleDeflection)
+
             self.refract(fresnelIntersection)
 
             mut1 = self.material.mu_t
diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py
index b89b623f..a973bdc6 100644
--- a/pytissueoptics/scene/intersection/intersectionFinder.py
+++ b/pytissueoptics/scene/intersection/intersectionFinder.py
@@ -23,7 +23,8 @@ class Intersection:
     outsideEnvironment: Environment = None
     surfaceLabel: str = None
     distanceLeft: float = None
-    isTooClose: bool = False
+    isSmooth: bool = False
+    rawNormal: Vector = None
 
 
 class IntersectionFinder:
@@ -64,20 +65,24 @@ def _composeIntersection(ray: Ray, intersection: Intersection) -> Optional[Inter
         if not intersection:
             return None
 
-        smoothNormal = shader.getSmoothNormal(intersection.polygon, intersection.position)
-
-        # If the resulting smooth normal changes the sign of the dot product with the ray direction, do not smooth.
-        if smoothNormal.dot(ray.direction) * intersection.polygon.normal.dot(ray.direction) < 0:
-            smoothNormal = intersection.polygon.normal
-
-        intersection.normal = smoothNormal
         intersection.insideEnvironment = intersection.polygon.insideEnvironment
         intersection.outsideEnvironment = intersection.polygon.outsideEnvironment
         intersection.surfaceLabel = intersection.polygon.surfaceLabel
+        intersection.rawNormal = intersection.polygon.normal
 
         if ray.length is not None:
             intersection.distanceLeft = ray.length - intersection.distance
 
+        smoothNormal = shader.getSmoothNormal(intersection.polygon, intersection.position)
+
+        # If the resulting smooth normal changes the sign of the dot product with the ray direction, do not smooth.
+        if smoothNormal.dot(ray.direction) * intersection.polygon.normal.dot(ray.direction) < 0:
+            intersection.normal = intersection.polygon.normal
+            intersection.isSmooth = False
+            return intersection
+
+        intersection.normal = smoothNormal
+        intersection.isSmooth = True
         return intersection
 
 

From 8fd0f975e082fd26be63513fe487ab61998521f3 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Wed, 20 Nov 2024 20:54:58 -0500
Subject: [PATCH 10/23] fix backward catch: ignore if ray is not over triangle

---
 .../rayscattering/opencl/src/intersection.c   | 15 ++++++++++-
 .../intersection/mollerTrumboreIntersect.py   | 26 ++++++++++++++-----
 2 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index f0e8a3ed..718b177d 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -211,7 +211,20 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
         hitPoint.position = ray.origin + ray.length * ray.direction;
         return hitPoint;
     }
-    if (t < 0 && dt_T < EPS_CATCH) {
+    if (t <= 0 && dt_T < EPS_CATCH) {
+        // Create a test ray to compute if the origin lies inside the triangle (from the normal).
+        pVector = cross(normal, edgeB);
+        det = dot(edgeA, pVector);
+        invDet = 1.0f / det;
+        u = dot(tVector, pVector) * invDet;
+        if (u < 0.0f || u > 1.0f) {
+            return hitPoint;
+        }
+
+        v = dot(normal, qVector) * invDet;
+        if (v < 0.0f || u + v > 1.0f) {
+            return hitPoint;
+        }
         hitPoint.exists = true;
         hitPoint.distance = 0;
         hitPoint.position = ray.origin;
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index 82cdcd0a..306d4464 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -24,11 +24,13 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
         """ Möller–Trumbore ray-triangle 3D intersection algorithm.
         Added epsilon zones to avoid numerical errors in the OpenCL implementation.
         Modified to support rays with finite length:
-            A. If the intersection is too far away, return None.
-            B. If the intersection is just a bit too far away, there is no intersection, but we cannot accept
-                the ray to land there (it would be too close to the surface and might yield intersection errors
-                later on). Therefore, we tag the intersection as 'tooClose' and use that later to move the ray's
-                landing position a bit away from this surface.
+            A. If the intersection is too far away, do not intersect.
+            B. (Forward catch) If the intersection is just a bit too far away, such that the resulting shortest distance
+                to surface is under epsilon, we must intersect already to prevent floating point errors in the next
+                intersection search.
+            C. (Backward catch) If, after a forward catch, the photon attempts to scatter back (inside epsilon region)
+                before actually crossing the surface, we must trigger an intersection event (at the origin). Ignore if
+                the origin does not lie over the triangle (meaning that the photon already crossed another surface).
         """
         v1, v2, v3 = triangle.vertices
         edgeA = v2 - v1
@@ -69,9 +71,21 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
         if t > ray.length and dt_T < self.EPS_CATCH:
             # Case 2: Forward epsilon catch. Ray ends close to the triangle, so we intersect at the ray's end.
             return ray.origin + ray.direction * ray.length
-        if t < 0 and dt_T < self.EPS_CATCH:
+        if t <= 0 and dt_T < self.EPS_CATCH:
             # Case 3: Backward epsilon catch. Ray starts close to the triangle, so we intersect at the origin.
             # This requires the intersector to always test triangles (or at least, close ones) of the origin solid.
+
+            # Do not catch if the origin lies outside the triangle (intersection test using triangle normal).
+            pVector = triangle.normal.cross(edgeB)
+            determinant = edgeA.dot(pVector)
+            inverseDeterminant = 1. / determinant
+            u = tVector.dot(pVector) * inverseDeterminant
+            if u < 0 or u > 1:
+                return None
+            v = triangle.normal.dot(qVector) * inverseDeterminant
+            if v < 0 or u + v > 1:
+                return None
+
             return ray.origin
 
         # Case 4: No intersection.

From acdcccdb36d016fb8566d8923d76f1eed6a07075 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Wed, 20 Nov 2024 20:58:31 -0500
Subject: [PATCH 11/23] revert spin to use original rotate around method

---
 pytissueoptics/rayscattering/photon.py  |  3 +-
 pytissueoptics/scene/geometry/vector.py | 38 ++++++++-----------------
 2 files changed, 14 insertions(+), 27 deletions(-)

diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index 07438f9f..baa34f47 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -158,7 +158,8 @@ def scatter(self):
         self.interact()
 
     def scatterBy(self, theta, phi):
-        self._direction.spin(theta, phi)
+        self._er.rotateAround(self._direction, phi)
+        self._direction.rotateAround(self._er, theta)
 
     def interact(self):
         delta = self._weight * self.material.getAlbedo()
diff --git a/pytissueoptics/scene/geometry/vector.py b/pytissueoptics/scene/geometry/vector.py
index 41dfaa52..25e90d77 100644
--- a/pytissueoptics/scene/geometry/vector.py
+++ b/pytissueoptics/scene/geometry/vector.py
@@ -1,7 +1,5 @@
 import math
 
-COS_ZERO = 1 - 1e-12
-
 
 class Vector:
     """
@@ -98,35 +96,23 @@ def update(self, x: float, y: float, z: float):
     def copy(self) -> 'Vector':
         return Vector(self._x, self._y, self._z)
 
-    def spin(self, theta: float, phi: float):
-        """ In spherical coordinates. Taken directly from MCML code. """
-        cosp = math.cos(phi)
-        if phi < math.pi:
-            sinp = math.sqrt(1 - cosp*cosp)
-        else:
-            sinp = -math.sqrt(1 - cosp*cosp)
-
-        cost = math.cos(theta)
-        sint = math.sqrt(1 - cost*cost)
-
-        if abs(self._z) > COS_ZERO:
-            ux = sint*cosp
-            uy = sint*sinp
-            uz = cost * (1 if self._z >= 0 else -1)
-        else:
-            temp = math.sqrt(1 - self._z*self._z)
-            ux = sint*(self._x*self._z*cosp - self._y*sinp) / temp + self._x*cost
-            uy = sint*(self._y*self._z*cosp + self._x*sinp) / temp + self._y*cost
-            uz = -sint*cosp*temp + self._z*cost
-
-        self.update(ux, uy, uz)
-        self.normalize()
-
     def rotateAround(self, unitAxis: 'Vector', theta: float):
         """
         Rotate the vector around `unitAxis` by `theta` radians. Assumes the axis to be a unit vector.
         Uses Rodrigues' rotation formula.
         """
+        # This is the most expensive (and most common)
+        # operation when performing Monte Carlo in tissue
+        # (15% of time spent here). It is difficult to optimize without
+        # making it even less readable than it currently is
+        # http://en.wikipedia.org/wiki/Rotation_matrix
+        #
+        # Several options were tried in the past such as
+        # external not-so-portable C library, unreadable
+        # shortcuts, sine and cosine lookup tables, etc...
+        # and the performance gain was minimal (<20%).
+        # For now, this is the best, most readable solution.
+
         cost = math.cos(theta)
         sint = math.sin(theta)
         one_cost = 1 - cost

From 1ea0203a0a20def7cdf525e4b7d1b031ffc0e728 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Wed, 20 Nov 2024 20:59:45 -0500
Subject: [PATCH 12/23] (opencl) revert spin method as well

---
 .../rayscattering/opencl/src/propagation.c    |  3 ++-
 .../opencl/src/vectorOperators.c              | 24 -------------------
 2 files changed, 2 insertions(+), 25 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index a3845941..5f976230 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -13,7 +13,8 @@ void moveBy(float distance, __global Photon *photons, uint photonID){
 }
 
 void scatterBy(float phi, float theta, __global Photon *photons, uint photonID){
-    spin(&photons[photonID].direction, theta, phi);
+    rotateAroundAxisGlobal(&photons[photonID].er, &photons[photonID].direction, phi);
+    rotateAroundAxisGlobal(&photons[photonID].direction, &photons[photonID].er, theta);
 }
 
 void decreaseWeightBy(float delta_weight, __global Photon *photons, uint photonID){
diff --git a/pytissueoptics/rayscattering/opencl/src/vectorOperators.c b/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
index 554dc14c..32dec7c2 100644
--- a/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
+++ b/pytissueoptics/rayscattering/opencl/src/vectorOperators.c
@@ -1,4 +1,3 @@
-__constant float COS_ZERO = 0.999999f;
 
 void normalizeVectorLocal(float3 *vector){
     float length = sqrt(vector->x * vector->x + vector->y * vector->y + vector->z * vector->z);
@@ -68,29 +67,6 @@ void rotateAround(__global float3 *mainVector, float3 *axisVector, float theta){
     mainVector->z = z;
 }
 
-void spin(__global float3 *mainVector, float theta, float phi){
-    float cosp = cos(phi);
-    float sinp = phi < M_PI ? sqrt(1 - cosp*cosp) : -sqrt(1 - cosp*cosp);
-    float cost = cos(theta);
-    float sint = sqrt(1 - cost*cost);
-    float ux, uy, uz;
-    if (fabs(mainVector->z) > COS_ZERO){
-        ux = sint*cosp;
-        uy = sint*sinp;
-        uz = cost * (mainVector->z >= 0 ? 1 : -1);
-    }
-    else{
-        float temp = sqrt(1 - mainVector->z*mainVector->z);
-        ux = sint*(mainVector->x*mainVector->z*cosp - mainVector->y*sinp) / temp + mainVector->x*cost;
-        uy = sint*(mainVector->y*mainVector->z*cosp + mainVector->x*sinp) / temp + mainVector->y*cost;
-        uz = -sint*cosp*temp + mainVector->z*cost;
-    }
-    mainVector->x = ux;
-    mainVector->y = uy;
-    mainVector->z = uz;
-    normalizeVectorGlobal(mainVector);
-}
-
 float3 getAnyOrthogonalGlobal(__global float3 *vector){
     if (fabs(vector->z) < fabs(vector->x)){
         return (float3)(vector->y, -vector->x, 0.0f);

From 4780ffb41874cda7776d7ed2cf60f7b45dc78c0e Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Wed, 20 Nov 2024 21:08:51 -0500
Subject: [PATCH 13/23] fix rotateAround by recalculating orthogonal vector

---
 pytissueoptics/rayscattering/opencl/src/propagation.c | 1 +
 pytissueoptics/rayscattering/photon.py                | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index 5f976230..925752ab 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -15,6 +15,7 @@ void moveBy(float distance, __global Photon *photons, uint photonID){
 void scatterBy(float phi, float theta, __global Photon *photons, uint photonID){
     rotateAroundAxisGlobal(&photons[photonID].er, &photons[photonID].direction, phi);
     rotateAroundAxisGlobal(&photons[photonID].direction, &photons[photonID].er, theta);
+    photons[photonID].er = getAnyOrthogonalGlobal(&photons[photonID].direction);
 }
 
 void decreaseWeightBy(float delta_weight, __global Photon *photons, uint photonID){
diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index baa34f47..2f84b440 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -24,7 +24,6 @@ def __init__(self, position: Vector, direction: Vector):
         self._environment: Environment = None
 
         self._er = self._direction.getAnyOrthogonal()
-        self._er.normalize()
         self._hasContext = False
         self._fresnelIntersect: FresnelIntersect = None
 
@@ -160,6 +159,7 @@ def scatter(self):
     def scatterBy(self, theta, phi):
         self._er.rotateAround(self._direction, phi)
         self._direction.rotateAround(self._er, theta)
+        self._er = self._direction.getAnyOrthogonal()
 
     def interact(self):
         delta = self._weight * self.material.getAlbedo()

From dccbed51f56a502def72829bce670ff9a49cd15d Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Wed, 20 Nov 2024 21:13:22 -0500
Subject: [PATCH 14/23] remove deprecated field isTooClose

---
 pytissueoptics/rayscattering/opencl/src/fresnel.c      | 1 -
 pytissueoptics/rayscattering/opencl/src/intersection.c | 3 ---
 pytissueoptics/scene/intersection/ray.py               | 2 --
 3 files changed, 6 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/fresnel.c b/pytissueoptics/rayscattering/opencl/src/fresnel.c
index 1decf6f0..1ada0389 100644
--- a/pytissueoptics/rayscattering/opencl/src/fresnel.c
+++ b/pytissueoptics/rayscattering/opencl/src/fresnel.c
@@ -107,7 +107,6 @@ FresnelIntersection computeFresnelIntersection(float3 rayDirection, Intersection
 Intersection getLocalIntersection(__global Intersection *intersections, uint gid) {
     Intersection intersection;
     intersection.exists = intersections[gid].exists;
-    intersection.isTooClose = intersections[gid].isTooClose;
     intersection.distance = intersections[gid].distance;
     intersection.position = intersections[gid].position;
     intersection.normal = intersections[gid].normal;
diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 718b177d..a3990ff3 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -164,7 +164,6 @@ typedef struct HitPoint HitPoint;
 HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, float3 normal) {
     HitPoint hitPoint;
     hitPoint.exists = false;
-    hitPoint.isTooClose = false;
 
     float3 edgeA = v2 - v1;
     float3 edgeB = v3 - v1;
@@ -240,7 +239,6 @@ Intersection _findClosestPolygonIntersection(Ray ray, uint solidID,
                                             uint photonSolidID) {
     Intersection intersection;
     intersection.exists = false;
-    intersection.isTooClose = false;
     intersection.distance = INFINITY;
     for (uint s = solids[solidID-1].firstSurfaceID; s <= solids[solidID-1].lastSurfaceID; s++) {
         for (uint p = surfaces[s].firstPolygonID; p <= surfaces[s].lastPolygonID; p++) {
@@ -351,7 +349,6 @@ Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidI
 
     Intersection closestIntersection;
     closestIntersection.exists = false;
-    closestIntersection.isTooClose = false;
     closestIntersection.distance = INFINITY;
     if (scene->nSolids == 0) {
         return closestIntersection;
diff --git a/pytissueoptics/scene/intersection/ray.py b/pytissueoptics/scene/intersection/ray.py
index 1f336d8f..ec80c3af 100644
--- a/pytissueoptics/scene/intersection/ray.py
+++ b/pytissueoptics/scene/intersection/ray.py
@@ -8,8 +8,6 @@ def __init__(self, origin: Vector, direction: Vector, length: float = None):
         self._direction.normalize()
         self._length = length
 
-        self.isTooClose = False
-
     @property
     def origin(self):
         return self._origin

From 8b2d6a751e354e75cd6ccea83a6f4fdecde28b4f Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 11:34:26 -0500
Subject: [PATCH 15/23] fix solid intersection test skipper based on bbox

simply break when bbox distance is greater than current current intersection distance.
---
 pytissueoptics/rayscattering/opencl/src/intersection.c |  6 ++++--
 .../scene/intersection/intersectionFinder.py           | 10 +++-------
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index a3990ff3..89fdf20f 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -360,8 +360,10 @@ Intersection findIntersection(Ray ray, Scene *scene, uint gid, uint photonSolidI
             // Default buffer value -1 means that there is no intersection with this solid
             continue;
         }
-        bool contained = scene->solidCandidates[boxGID].distance == 0;
-        if (!contained && closestIntersection.exists) {
+
+        if (scene->solidCandidates[boxGID].distance > closestIntersection.distance) {
+            // The solid candidates are sorted by distance, so we can break early if the BBox distance
+            // is greater than the closest intersection found so far.
             break;
         }
 
diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py
index a973bdc6..c17bfe12 100644
--- a/pytissueoptics/scene/intersection/intersectionFinder.py
+++ b/pytissueoptics/scene/intersection/intersectionFinder.py
@@ -96,11 +96,8 @@ def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Interse
                 possible, and we simply set the bbox intersection distance to zero.
         2. Sort these solid candidates by bbox intersection distance.
         3. For each solid, find the closest polygon intersection.
-            3.1 If a polygon intersection is found for a given solid candidate, then we return this intersection point
-            without testing the other solid candidates since they are ordered by distance.
-            N.B.: Except for the case of multiple contained solids where it is not possible to order them in a
-            meaningful way, so in this case we need to test all of them (candidate distance is zero) before
-            returning the closest intersection found.
+            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.sort(key=lambda x: x[0])
@@ -108,8 +105,7 @@ def findIntersection(self, ray: Ray, currentSolidLabel: str) -> Optional[Interse
         closestDistance = sys.maxsize
         closestIntersection = None
         for i, (distance, solid) in enumerate(bboxIntersections):
-            contained = distance == 0
-            if not contained and closestIntersection:
+            if distance > closestDistance:
                 break
             intersection = self._findClosestPolygonIntersection(ray, solid.getPolygons(), currentSolidLabel)
             if intersection and intersection.distance < closestDistance:

From 8dbe3cb40b114038727e7393d8cd6dae60b00b87 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 12:03:14 -0500
Subject: [PATCH 16/23] fix forward catch: move photon to hitpoint and support
 negative distance remainder.

---
 .../rayscattering/opencl/src/intersection.c      |  8 ++++----
 .../rayscattering/opencl/src/propagation.c       | 16 ++++++++++++----
 pytissueoptics/rayscattering/photon.py           | 12 +++++++++---
 .../intersection/mollerTrumboreIntersect.py      |  8 +++++---
 4 files changed, 30 insertions(+), 14 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 89fdf20f..6fa645d1 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -189,11 +189,11 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
     }
 
     float t = dot(edgeB, qVector) * invDet;
+    hitPoint.distance = t;
+    hitPoint.position = ray.origin + t * ray.direction;
 
     if (t > 0 && ray.length >= t){
         hitPoint.exists = true;
-        hitPoint.distance = t;
-        hitPoint.position = ray.origin + t * ray.direction;
         return hitPoint;
     }
 
@@ -204,10 +204,10 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
         dt = t - ray.length;
     }
     float dt_T = fabs(dot(normal, ray.direction) * dt);
+
     if (t > ray.length && dt_T < EPS_CATCH) {
+        // Forward catch.
         hitPoint.exists = true;
-        hitPoint.distance = ray.length;
-        hitPoint.position = ray.origin + ray.length * ray.direction;
         return hitPoint;
     }
     if (t <= 0 && dt_T < EPS_CATCH) {
diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index 925752ab..e67a7a39 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -12,6 +12,10 @@ void moveBy(float distance, __global Photon *photons, uint photonID){
     photons[photonID].position += (distance * photons[photonID].direction);
 }
 
+void moveTo(float3 position, __global Photon *photons, uint photonID){
+    photons[photonID].position = position;
+}
+
 void scatterBy(float phi, float theta, __global Photon *photons, uint photonID){
     rotateAroundAxisGlobal(&photons[photonID].er, &photons[photonID].direction, phi);
     rotateAroundAxisGlobal(&photons[photonID].direction, &photons[photonID].er, theta);
@@ -141,10 +145,14 @@ float reflectOrRefract(Intersection *intersection, __global Photon *photons, __c
 float propagateStep(float distance, __global Photon *photons, __constant Material *materials, Scene *scene,
                     __global uint *seeds, __global DataPoint *logger, uint *logIndex, uint gid, uint photonID){
 
-    if (distance == 0) {
+    if (distance <= 0) {
         float mu_t = materials[photons[photonID].materialID].mu_t;
         float randomNumber = getRandomFloatValue(seeds, gid);
-        distance = getScatteringDistance(mu_t, randomNumber);
+        distance += getScatteringDistance(mu_t, randomNumber);
+        if (distance < 0){
+            // Not really possible until mu_t is very high (> 1000) and intense smoothing is applied (order-1 spheres).
+            distance = 0;
+        }
     }
 
     Ray stepRay = {photons[photonID].position, photons[photonID].direction, distance};
@@ -153,7 +161,7 @@ float propagateStep(float distance, __global Photon *photons, __constant Materia
     float distanceLeft = 0;
 
     if (intersection.exists){
-        moveBy(intersection.distance, photons, photonID);
+        moveTo(intersection.position, photons, photonID);
         distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID);
     } else {
         if (distance == INFINITY){
@@ -197,7 +205,7 @@ __kernel void propagate(uint maxPhotons, uint maxInteractions, float weightThres
             distance = propagateStep(distance, photons, materials, &scene,
                                      seeds, logger, &logIndex, gid, currentPhotonIndex);
             roulette(weightThreshold, photons, seeds, gid, currentPhotonIndex);
-            }
+        }
         photonCount++;
     }
 }
diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index 2f84b440..dde06f66 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -74,13 +74,16 @@ def propagate(self):
             self.roulette()
 
     def step(self, distance=0) -> float:
-        if distance == 0:
-            distance = self.material.getScatteringDistance()
+        if distance <= 0:
+            distance += self.material.getScatteringDistance()
+            if distance < 0:
+                # Not really possible until mu_t is very high (> 1000) and intense smoothing is applied (order-1 spheres).
+                distance = 0
 
         intersection = self._getIntersection(distance)
 
         if intersection:
-            self.moveBy(intersection.distance)
+            self.moveTo(intersection.position)
             distanceLeft = self.reflectOrRefract(intersection)
         else:
             if math.isinf(distance):
@@ -143,6 +146,9 @@ def _getFresnelIntersection(self, intersection: Intersection) -> FresnelIntersec
     def moveBy(self, distance):
         self._position += self._direction * distance
 
+    def moveTo(self, position: Vector):
+        self._position = position
+
     def reflect(self, fresnelIntersection: FresnelIntersection):
         self._direction.rotateAround(fresnelIntersection.incidencePlane,
                                      fresnelIntersection.angleDeflection)
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index 306d4464..3974a05c 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -57,9 +57,11 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
 
         # Distance to intersection point
         t = edgeB.dot(qVector) * inverseDeterminant
+        hitPoint = ray.origin + ray.direction * t
+
         if t > 0 and (ray.length >= t or ray.length is None):
             # Case 1: Trivial case. Intersects.
-            return ray.origin + ray.direction * t
+            return hitPoint
 
         # Next we need to check if the intersection is inside the epsilon catch zone (forward or backward).
         # Note that this mechanic only works when same-solid intersections are ignored before calling this function.
@@ -69,8 +71,8 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
             dt = t - ray.length
         dt_T = abs(triangle.normal.dot(ray.direction) * dt)
         if t > ray.length and dt_T < self.EPS_CATCH:
-            # Case 2: Forward epsilon catch. Ray ends close to the triangle, so we intersect at the ray's end.
-            return ray.origin + ray.direction * ray.length
+            # Case 2: Forward epsilon catch. Ray ends too close to the triangle, so we intersect.
+            return hitPoint
         if t <= 0 and dt_T < self.EPS_CATCH:
             # Case 3: Backward epsilon catch. Ray starts close to the triangle, so we intersect at the origin.
             # This requires the intersector to always test triangles (or at least, close ones) of the origin solid.

From accb44204e4f932a9d9699d02c2e821349fec3cd Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 12:13:38 -0500
Subject: [PATCH 17/23] fix side-limit intersections: use larger virtual
 triangles then move inside true surface

---
 .../rayscattering/opencl/src/intersection.c   | 28 +++++++++++++++----
 .../intersection/mollerTrumboreIntersect.py   | 27 ++++++++++++------
 2 files changed, 42 insertions(+), 13 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 6fa645d1..a38539bd 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -1,6 +1,7 @@
-__constant float EPS_CATCH = 0.00001f;
-__constant float EPS_PARALLEL = 0.000001f;
-__constant float EPS_SIDE = 0.000001f;
+__constant float EPS_CATCH = 1e-7f;
+__constant float EPS_PARALLEL = 1e-6f;
+__constant float EPS_SIDE = 3e-6f;
+__constant float EPS = 1e-7;
 
 struct Intersection {
     uint exists;
@@ -184,7 +185,7 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
 
     float3 qVector = cross(tVector, edgeA);
     float v = dot(ray.direction, qVector) * invDet;
-    if (v < -EPS_SIDE || u + v > 1.0f) {
+    if (v < -EPS_SIDE || u + v > 1.0f + EPS_SIDE) {
         return hitPoint;
     }
 
@@ -192,7 +193,24 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
     hitPoint.distance = t;
     hitPoint.position = ray.origin + t * ray.direction;
 
-    if (t > 0 && ray.length >= t){
+    // Check if the intersection is slightly outside the triangle.
+    float error = 0;
+    if (u < -EPS){
+        error -= u;
+    }
+    if (v < -EPS){
+        error -= v;
+    }
+    if (u + v > 1.0 + EPS){
+        error += u + v - 1.0;
+    }
+    if (error > 0){
+        // Move the hit point towards the triangle center by this error factor.
+        float3 correction = v1 + v2 + v3 - hitPoint.position * 3;
+        hitPoint.position += 2.0f * error * correction;
+    }
+
+    if (t >= 0 && ray.length >= t){
         hitPoint.exists = true;
         return hitPoint;
     }
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index 3974a05c..c6ccca18 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -4,13 +4,11 @@
 from pytissueoptics.scene.intersection import Ray
 
 
-EPS_CORRECTION = 0.0  # TODO: remove
-
-
 class MollerTrumboreIntersect:
-    EPS_CATCH = 0.00001
-    EPS_PARALLEL = 0.000001
-    EPS_SIDE = 0.000001
+    EPS_CATCH = 1e-7
+    EPS_PARALLEL = 1e-6
+    EPS_SIDE = 3e-6
+    EPS = 1e-7
 
     def getIntersection(self, ray: Ray, polygon: Union[Triangle, Quad, Polygon]) -> Optional[Vector]:
         if isinstance(polygon, Triangle):
@@ -52,14 +50,27 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
 
         qVector = tVector.cross(edgeA)
         v = ray.direction.dot(qVector) * inverseDeterminant
-        if v < -self.EPS_SIDE or u + v > 1.:
+        if v < -self.EPS_SIDE or u + v > 1. + self.EPS_SIDE:
             return None
 
         # Distance to intersection point
         t = edgeB.dot(qVector) * inverseDeterminant
         hitPoint = ray.origin + ray.direction * t
 
-        if t > 0 and (ray.length >= t or ray.length is None):
+        # Check if the intersection is slightly outside the true triangle surface.
+        error = 0
+        if u < -self.EPS:
+            error -= u
+        if v < -self.EPS:
+            error -= v
+        if u + v > 1. + self.EPS:
+            error += u + v - 1.
+        if error > 0:
+            # Move the hit point towards the triangle center by this error factor.
+            correctionDirection = v1 + v2 + v3 - hitPoint * 3
+            hitPoint += correctionDirection * 2 * error
+
+        if t >= 0 and (ray.length >= t or ray.length is None):
             # Case 1: Trivial case. Intersects.
             return hitPoint
 

From 03cb1e5e30f9b172434822a023b3aec0757c8c45 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 12:36:34 -0500
Subject: [PATCH 18/23] fix back catch: tweak eps, use true hitpoint and
 replace locality test with same-solid distance condition

which cancels back catch (and all intersections for that matter) if there is an intersection further away that sends the photon towards its current environment.
---
 .../rayscattering/opencl/src/intersection.c   | 49 +++++++++++--------
 .../scene/intersection/intersectionFinder.py  | 24 ++++++---
 .../intersection/mollerTrumboreIntersect.py   | 21 +++-----
 3 files changed, 52 insertions(+), 42 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index a38539bd..af910f1f 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -1,4 +1,5 @@
 __constant float EPS_CATCH = 1e-7f;
+__constant float EPS_BACK_CATCH = 2e-6f;
 __constant float EPS_PARALLEL = 1e-6f;
 __constant float EPS_SIDE = 3e-6f;
 __constant float EPS = 1e-7;
@@ -228,23 +229,10 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
         hitPoint.exists = true;
         return hitPoint;
     }
-    if (t <= 0 && dt_T < EPS_CATCH) {
-        // Create a test ray to compute if the origin lies inside the triangle (from the normal).
-        pVector = cross(normal, edgeB);
-        det = dot(edgeA, pVector);
-        invDet = 1.0f / det;
-        u = dot(tVector, pVector) * invDet;
-        if (u < 0.0f || u > 1.0f) {
-            return hitPoint;
-        }
 
-        v = dot(normal, qVector) * invDet;
-        if (v < 0.0f || u + v > 1.0f) {
-            return hitPoint;
-        }
+    if (t < 0 && (t > -EPS_BACK_CATCH || dt_T < EPS_CATCH)) {
+        // Backward catch.
         hitPoint.exists = true;
-        hitPoint.distance = 0;
-        hitPoint.position = ray.origin;
         return hitPoint;
     }
 
@@ -258,20 +246,28 @@ Intersection _findClosestPolygonIntersection(Ray ray, uint solidID,
     Intersection intersection;
     intersection.exists = false;
     intersection.distance = INFINITY;
+
+    float minSameSolidDistance = -INFINITY;
+
     for (uint s = solids[solidID-1].firstSurfaceID; s <= solids[solidID-1].lastSurfaceID; s++) {
         for (uint p = surfaces[s].firstPolygonID; p <= surfaces[s].lastPolygonID; p++) {
+            uint vertexIDs[3] = {triangles[p].vertexIDs[0], triangles[p].vertexIDs[1], triangles[p].vertexIDs[2]};
+            HitPoint hitPoint = _getTriangleIntersection(ray, vertices[vertexIDs[0]].position, vertices[vertexIDs[1]].position, vertices[vertexIDs[2]].position, triangles[p].normal);
+
+            if (!hitPoint.exists) {
+                continue;
+            }
+
             bool isGoingInside = dot(ray.direction, triangles[p].normal) < 0;
             uint nextSolidID = isGoingInside ? surfaces[s].insideSolidID : surfaces[s].outsideSolidID;
             if (nextSolidID == photonSolidID) {
+                if (hitPoint.distance > minSameSolidDistance) {
+                    minSameSolidDistance = hitPoint.distance;
+                }
                 continue;
             }
 
-            uint vertexIDs[3] = {triangles[p].vertexIDs[0], triangles[p].vertexIDs[1], triangles[p].vertexIDs[2]};
-            HitPoint hitPoint = _getTriangleIntersection(ray, vertices[vertexIDs[0]].position, vertices[vertexIDs[1]].position, vertices[vertexIDs[2]].position, triangles[p].normal);
-            if (!hitPoint.exists) {
-                continue;
-            }
-            if (hitPoint.distance < intersection.distance) {
+            if (fabs(hitPoint.distance) < fabs(intersection.distance)) {
                 intersection.exists = true;
                 intersection.distance = hitPoint.distance;
                 intersection.position = hitPoint.position;
@@ -281,6 +277,17 @@ Intersection _findClosestPolygonIntersection(Ray ray, uint solidID,
             }
         }
     }
+
+    if (intersection.distance == 0 && minSameSolidDistance == 0){
+        // Cancel back catch. Surface overlap.
+        intersection.exists = false;
+    } else if (intersection.distance < 0) {
+        // Cancel backward catch if the same-solid intersect distance is greater.
+        if (minSameSolidDistance > intersection.distance + 1e-7) {
+            intersection.exists = false;
+        }
+    }
+
     return intersection;
 }
 
diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py
index c17bfe12..153334d3 100644
--- a/pytissueoptics/scene/intersection/intersectionFinder.py
+++ b/pytissueoptics/scene/intersection/intersectionFinder.py
@@ -40,24 +40,36 @@ def _findClosestPolygonIntersection(
         self, ray: Ray, polygons: List[Polygon], currentSolidLabel: str
     ) -> Optional[Intersection]:
         closestIntersection = Intersection(sys.maxsize)
+        minSameSolidDistance = -sys.maxsize
+
         for polygon in polygons:
-            # Skip intersection test if the ray is heading towards its current solid
+            intersectionPoint = self._polygonIntersect.getIntersection(ray, polygon)
+            if not intersectionPoint:
+                continue
+            distance = (intersectionPoint - ray.origin).getNorm()
+
+            # Discard intersection result if the ray is heading towards its current solid
             # (possible because of the epsilon catch zone in our Moller Trumbore intersect).
             isGoingInside = ray.direction.dot(polygon.normal) < 0
             nextSolid = polygon.insideEnvironment.solid if isGoingInside else polygon.outsideEnvironment.solid
             nextLabel = 'world' if nextSolid is None else nextSolid.getLabel()
             if nextLabel == currentSolidLabel:
+                minSameSolidDistance = max(minSameSolidDistance, distance)
                 continue
 
-            intersectionPoint = self._polygonIntersect.getIntersection(ray, polygon)
-            if not intersectionPoint:
-                continue
-            distance = (intersectionPoint - ray.origin).getNorm()
-            if distance < closestIntersection.distance:
+            if abs(distance) < abs(closestIntersection.distance):
                 closestIntersection = Intersection(distance, intersectionPoint, polygon)
 
         if closestIntersection.distance == sys.maxsize:
             return None
+
+        if closestIntersection.distance == 0 and minSameSolidDistance == 0:
+            # Cancel back catch. Surface overlap.
+            return None
+        if closestIntersection.distance < 0 and minSameSolidDistance > closestIntersection.distance + 1e-7:
+            # Cancel back catch if the same-solid intersect distance is greater.
+            return None
+
         return closestIntersection
 
     @staticmethod
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index c6ccca18..bec64342 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -6,6 +6,7 @@
 
 class MollerTrumboreIntersect:
     EPS_CATCH = 1e-7
+    EPS_BACK_CATCH = 2e-6
     EPS_PARALLEL = 1e-6
     EPS_SIDE = 3e-6
     EPS = 1e-7
@@ -81,25 +82,15 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
         else:
             dt = t - ray.length
         dt_T = abs(triangle.normal.dot(ray.direction) * dt)
+
         if t > ray.length and dt_T < self.EPS_CATCH:
             # Case 2: Forward epsilon catch. Ray ends too close to the triangle, so we intersect.
             return hitPoint
-        if t <= 0 and dt_T < self.EPS_CATCH:
-            # Case 3: Backward epsilon catch. Ray starts close to the triangle, so we intersect at the origin.
-            # This requires the intersector to always test triangles (or at least, close ones) of the origin solid.
 
-            # Do not catch if the origin lies outside the triangle (intersection test using triangle normal).
-            pVector = triangle.normal.cross(edgeB)
-            determinant = edgeA.dot(pVector)
-            inverseDeterminant = 1. / determinant
-            u = tVector.dot(pVector) * inverseDeterminant
-            if u < 0 or u > 1:
-                return None
-            v = triangle.normal.dot(qVector) * inverseDeterminant
-            if v < 0 or u + v > 1:
-                return None
-
-            return ray.origin
+        if t < 0 and (t > -self.EPS_BACK_CATCH or dt_T < self.EPS_CATCH):
+            # Case 3: Backward epsilon catch. Ray starts too close to the triangle, so we intersect.
+            # This requires the intersector to always test triangles (or at least, close ones) of the origin solid.
+            return hitPoint
 
         # Case 4: No intersection.
         return None

From f3a7c3910230d9910b59827ee6515e3e7bb08e30 Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 12:42:01 -0500
Subject: [PATCH 19/23] (opencl only) prioritize back catch when ray lies on
 the triangle.

---
 pytissueoptics/rayscattering/opencl/src/intersection.c | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index af910f1f..784bcff1 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -233,6 +233,11 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3, floa
     if (t < 0 && (t > -EPS_BACK_CATCH || dt_T < EPS_CATCH)) {
         // Backward catch.
         hitPoint.exists = true;
+
+        // If ray lies on the triangle, return a distance of 0 to prioritize this intersection.
+        if (dt_T < EPS) {
+            hitPoint.distance = 0;
+        }
         return hitPoint;
     }
 

From 4e73f3f14edd26a13fb972ab267c5cca17451dce Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 12:53:45 -0500
Subject: [PATCH 20/23] prevent solid interfaces to be intersected from outside

---
 .../rayscattering/opencl/src/intersection.c        |  6 ++++++
 pytissueoptics/rayscattering/photon.py             |  5 +----
 pytissueoptics/scene/geometry/polygon.py           |  8 ++++++++
 .../scene/intersection/intersectionFinder.py       | 14 ++++++++++----
 .../scene/intersection/mollerTrumboreIntersect.py  |  6 ++----
 5 files changed, 27 insertions(+), 12 deletions(-)

diff --git a/pytissueoptics/rayscattering/opencl/src/intersection.c b/pytissueoptics/rayscattering/opencl/src/intersection.c
index 784bcff1..a2da8216 100644
--- a/pytissueoptics/rayscattering/opencl/src/intersection.c
+++ b/pytissueoptics/rayscattering/opencl/src/intersection.c
@@ -255,6 +255,12 @@ Intersection _findClosestPolygonIntersection(Ray ray, uint solidID,
     float minSameSolidDistance = -INFINITY;
 
     for (uint s = solids[solidID-1].firstSurfaceID; s <= solids[solidID-1].lastSurfaceID; s++) {
+        // When an interface joins a side surface, an outside photon could try to intersect with the interface
+        //  while this is not allowed. So we skip these tests (where surface environments dont match the photon).
+        if (photonSolidID != surfaces[s].insideSolidID && photonSolidID != surfaces[s].outsideSolidID) {
+            continue;
+        }
+
         for (uint p = surfaces[s].firstPolygonID; p <= surfaces[s].lastPolygonID; p++) {
             uint vertexIDs[3] = {triangles[p].vertexIDs[0], triangles[p].vertexIDs[1], triangles[p].vertexIDs[2]};
             HitPoint hitPoint = _getTriangleIntersection(ray, vertices[vertexIDs[0]].position, vertices[vertexIDs[1]].position, vertices[vertexIDs[2]].position, triangles[p].normal);
diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index dde06f66..6a3300bf 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -11,7 +11,6 @@
 from pytissueoptics.scene.intersection.intersectionFinder import IntersectionFinder, Intersection
 from pytissueoptics.scene.logger import Logger, InteractionKey
 
-WORLD_LABEL = "world"
 WEIGHT_THRESHOLD = 1e-4
 MIN_ANGLE = 0.0001
 
@@ -52,9 +51,7 @@ def material(self) -> ScatteringMaterial:
 
     @property
     def solidLabel(self):
-        if not self._environment.solid:
-            return WORLD_LABEL
-        return self._environment.solid.getLabel()
+        return self._environment.solidLabel
 
     def setContext(self, environment: Environment, intersectionFinder: IntersectionFinder = None, logger: Logger = None,
                    fresnelIntersect=FresnelIntersect()):
diff --git a/pytissueoptics/scene/geometry/polygon.py b/pytissueoptics/scene/geometry/polygon.py
index fdee03c8..34640e30 100644
--- a/pytissueoptics/scene/geometry/polygon.py
+++ b/pytissueoptics/scene/geometry/polygon.py
@@ -4,12 +4,20 @@
 from pytissueoptics.scene.geometry import Vector, Vertex
 from pytissueoptics.scene.geometry import BoundingBox
 
+WORLD_LABEL = "world"
+
 
 @dataclass
 class Environment:
     material: ...
     solid: 'Solid' = None
 
+    @property
+    def solidLabel(self) -> str:
+        if self.solid:
+            return self.solid.getLabel()
+        return WORLD_LABEL
+
 
 class Polygon:
     """
diff --git a/pytissueoptics/scene/intersection/intersectionFinder.py b/pytissueoptics/scene/intersection/intersectionFinder.py
index 153334d3..1f7a68e6 100644
--- a/pytissueoptics/scene/intersection/intersectionFinder.py
+++ b/pytissueoptics/scene/intersection/intersectionFinder.py
@@ -43,17 +43,23 @@ def _findClosestPolygonIntersection(
         minSameSolidDistance = -sys.maxsize
 
         for polygon in polygons:
+            # When an interface joins a side surface, an outside photon could try to intersect with the interface.
+            #  This is not allowed, so we skip these tests (where surface environments dont match the photon).
+            insideLabel = polygon.insideEnvironment.solidLabel
+            outsideLabel = polygon.outsideEnvironment.solidLabel
+            if insideLabel != currentSolidLabel and outsideLabel != currentSolidLabel:
+                continue
+
             intersectionPoint = self._polygonIntersect.getIntersection(ray, polygon)
             if not intersectionPoint:
                 continue
             distance = (intersectionPoint - ray.origin).getNorm()
 
             # Discard intersection result if the ray is heading towards its current solid
-            # (possible because of the epsilon catch zone in our Moller Trumbore intersect).
+            # (possible because of the epsilon catch zone in our Moller-Trumbore intersect).
             isGoingInside = ray.direction.dot(polygon.normal) < 0
-            nextSolid = polygon.insideEnvironment.solid if isGoingInside else polygon.outsideEnvironment.solid
-            nextLabel = 'world' if nextSolid is None else nextSolid.getLabel()
-            if nextLabel == currentSolidLabel:
+            nextSolidLabel = insideLabel if isGoingInside else outsideLabel
+            if nextSolidLabel == currentSolidLabel:
                 minSameSolidDistance = max(minSameSolidDistance, distance)
                 continue
 
diff --git a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
index bec64342..eecfaea2 100644
--- a/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
+++ b/pytissueoptics/scene/intersection/mollerTrumboreIntersect.py
@@ -25,11 +25,9 @@ def _getTriangleIntersection(self, ray: Ray, triangle: Triangle) -> Optional[Vec
         Modified to support rays with finite length:
             A. If the intersection is too far away, do not intersect.
             B. (Forward catch) If the intersection is just a bit too far away, such that the resulting shortest distance
-                to surface is under epsilon, we must intersect already to prevent floating point errors in the next
+                to surface is under epsilon, we must intersect to prevent floating point errors in the next
                 intersection search.
-            C. (Backward catch) If, after a forward catch, the photon attempts to scatter back (inside epsilon region)
-                before actually crossing the surface, we must trigger an intersection event (at the origin). Ignore if
-                the origin does not lie over the triangle (meaning that the photon already crossed another surface).
+            C. (Backward catch) If the photon attempts to scatter back before actually crossing the surface, we must trigger an intersection event.
         """
         v1, v2, v3 = triangle.vertices
         edgeA = v2 - v1

From e32cc47020af3dcb4e801a190658c6950b291ebb Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 12:54:34 -0500
Subject: [PATCH 21/23] always define vertex normals independently of smoothing
 conditions

---
 pytissueoptics/scene/solids/solid.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/pytissueoptics/scene/solids/solid.py b/pytissueoptics/scene/solids/solid.py
index d0e64bfc..fe813194 100644
--- a/pytissueoptics/scene/solids/solid.py
+++ b/pytissueoptics/scene/solids/solid.py
@@ -37,6 +37,7 @@ def __init__(self, vertices: List[Vertex], position: Vector = Vector(0, 0, 0),
         self._resetPolygonsCentroids()
 
         self._smoothing = False
+        self._setVertexNormals()
         if smooth:
             self.smooth()
 
@@ -281,7 +282,10 @@ def smooth(self, surfaceLabel: str = None, reset: bool = True):
         be changed by overwriting the signature with a specific surfaceLabel in
         another solid implementation and calling super().smooth(surfaceLabel).
         """
-        self._smoothing = True
+        self._setVertexNormals(surfaceLabel, smooth=True, reset=reset)
+
+    def _setVertexNormals(self, surfaceLabel: str = None, smooth=False, reset=True):
+        self._smoothing = smooth
         if reset:
             for vertex in self.vertices:
                 vertex.normal = None
@@ -289,7 +293,7 @@ def smooth(self, surfaceLabel: str = None, reset: bool = True):
         polygons = self.getPolygons(surfaceLabel)
 
         for polygon in polygons:
-            polygon.toSmooth = True
+            polygon.toSmooth = smooth
             for vertex in polygon.vertices:
                 if vertex.normal:
                     vertex.normal += polygon.normal

From 20b989b2c9a6e4e1e1c99064c63542ff458bf3de Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 13:08:36 -0500
Subject: [PATCH 22/23] prevent photons to land on a vertex

---
 .../rayscattering/opencl/src/propagation.c    | 22 +++++++++++++++++++
 pytissueoptics/rayscattering/photon.py        | 14 ++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/pytissueoptics/rayscattering/opencl/src/propagation.c b/pytissueoptics/rayscattering/opencl/src/propagation.c
index e67a7a39..d64c11bb 100644
--- a/pytissueoptics/rayscattering/opencl/src/propagation.c
+++ b/pytissueoptics/rayscattering/opencl/src/propagation.c
@@ -163,6 +163,28 @@ float propagateStep(float distance, __global Photon *photons, __constant Materia
     if (intersection.exists){
         moveTo(intersection.position, photons, photonID);
         distanceLeft = reflectOrRefract(&intersection, photons, materials, scene->surfaces, logger, logIndex, seeds, gid, photonID);
+
+        // Check if intersection lies too close to a vertex.
+        int closeToVertexID = -1;
+        for (uint i = 0; i < 3; i++) {
+            uint vertexID = scene->triangles[intersection.polygonID].vertexIDs[i];
+            if (length(intersection.position - scene->vertices[vertexID].position) < 3e-7) {
+                closeToVertexID = vertexID;
+                break;
+            }
+        }
+
+        // If too close to a vertex, move photon away slightly.
+        if (closeToVertexID != -1) {
+            int stepSign = 1;
+            int solidIDTowardsNormal = scene->surfaces[intersection.surfaceID].outsideSolidID;
+            if (solidIDTowardsNormal != photons[photonID].solidID) {
+                stepSign = -1;
+            }
+            float3 stepCorrection = stepSign * scene->vertices[closeToVertexID].normal * EPS_CATCH;
+            photons[photonID].position += stepCorrection;
+        }
+
     } else {
         if (distance == INFINITY){
             photons[photonID].weight = 0;
diff --git a/pytissueoptics/rayscattering/photon.py b/pytissueoptics/rayscattering/photon.py
index 6a3300bf..7945762f 100644
--- a/pytissueoptics/rayscattering/photon.py
+++ b/pytissueoptics/rayscattering/photon.py
@@ -82,6 +82,20 @@ def step(self, distance=0) -> float:
         if intersection:
             self.moveTo(intersection.position)
             distanceLeft = self.reflectOrRefract(intersection)
+
+            # Check if intersection lies too close to a vertex.
+            for vertex in intersection.polygon.vertices:
+                if (intersection.position - vertex).getNorm() > 3e-7:
+                    continue
+                # If too close to a vertex, move photon away slightly.
+                stepSign = 1
+                solidLabelTowardsNormal = intersection.outsideEnvironment.solidLabel
+                if solidLabelTowardsNormal != self.solidLabel:
+                    stepSign = -1
+                stepCorrection = vertex.normal * stepSign * 1e-7
+                self._position += stepCorrection
+                break
+
         else:
             if math.isinf(distance):
                 self._weight = 0

From bc3dc03a309065c4828d8fb9f2eda2b5fdfabc0e Mon Sep 17 00:00:00 2001
From: jlbegin <ludovick.begin@gmail.com>
Date: Sat, 14 Dec 2024 13:25:18 -0500
Subject: [PATCH 23/23] update Source: add proper seed support by skipping IPP
 update

---
 pytissueoptics/rayscattering/source.py | 24 ++++++++++++++++--------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/pytissueoptics/rayscattering/source.py b/pytissueoptics/rayscattering/source.py
index 72d00303..69318f08 100644
--- a/pytissueoptics/rayscattering/source.py
+++ b/pytissueoptics/rayscattering/source.py
@@ -1,4 +1,5 @@
 import hashlib
+import random
 import time
 from typing import List, Union, Optional, Tuple
 import numpy as np
@@ -21,9 +22,14 @@
 
 
 class Source(Displayable):
-    def __init__(self, position: Vector, N: int, useHardwareAcceleration: bool = True, displaySize: float = 0.1):
+    def __init__(self, position: Vector, N: int, useHardwareAcceleration: bool = True, displaySize: float = 0.1, seed: Optional[int] = None):
         self._position = position
         self._N = N
+        self._seed = seed
+        if seed is not None:
+            np.random.seed(seed)
+            random.seed(seed)
+
         self._photons: Union[List[Photon], CLPhotons] = []
         self._environment = None
         self.displaySize = displaySize
@@ -41,7 +47,9 @@ def propagate(self, scene: ScatteringScene, logger: Logger = None, showProgress:
         if self._useHardwareAcceleration:
             IPP = self._getAverageInteractionsPerPhoton(scene)
             self._propagateOpenCL(IPP, scene, logger, showProgress)
-            self._updateIPP(scene, logger)
+            if self._seed is None:
+                # Do not update IPP if the seed is set, since it will alter batch statistics.
+                self._updateIPP(scene, logger)
         else:
             self._propagateCPU(scene, logger, showProgress)
 
@@ -182,7 +190,7 @@ def __hash__(self):
 
 class DirectionalSource(Source):
     def __init__(self, position: Vector, direction: Vector, diameter: float, N: int,
-                 useHardwareAcceleration: bool = True, displaySize: float = 0.1):
+                 useHardwareAcceleration: bool = True, displaySize: float = 0.1, seed: Optional[int] = None):
         self._diameter = diameter
         self._direction = direction
         self._direction.normalize()
@@ -190,7 +198,7 @@ def __init__(self, position: Vector, direction: Vector, diameter: float, N: int,
         self._xAxis.normalize()
         self._yAxis = self._direction.cross(self._xAxis)
         self._yAxis.normalize()
-        super().__init__(position=position, N=N, useHardwareAcceleration=useHardwareAcceleration, displaySize=displaySize)
+        super().__init__(position=position, N=N, useHardwareAcceleration=useHardwareAcceleration, displaySize=displaySize, seed=seed)
 
     def getInitialPositionsAndDirections(self) -> Tuple[np.ndarray, np.ndarray]:
         positions = self._getInitialPositions()
@@ -239,9 +247,9 @@ def _hashComponents(self) -> tuple:
 
 
 class PencilPointSource(DirectionalSource):
-    def __init__(self, position: Vector, direction: Vector, N: int, useHardwareAcceleration: bool = True, displaySize: float = 0.1):
+    def __init__(self, position: Vector, direction: Vector, N: int, useHardwareAcceleration: bool = True, displaySize: float = 0.1, seed: Optional[int] = None):
         super().__init__(position=position, direction=direction, diameter=0, N=N,
-                         useHardwareAcceleration=useHardwareAcceleration, displaySize=displaySize)
+                         useHardwareAcceleration=useHardwareAcceleration, displaySize=displaySize, seed=seed)
 
 
 class IsotropicPointSource(Source):
@@ -258,11 +266,11 @@ def _hashComponents(self) -> tuple:
 
 class DivergentSource(DirectionalSource):
     def __init__(self, position: Vector, direction: Vector, diameter: float, divergence: float, N: int,
-                 useHardwareAcceleration: bool = True, displaySize: float = 0.1):
+                 useHardwareAcceleration: bool = True, displaySize: float = 0.1, seed: Optional[int] = None):
         self._divergence = divergence
 
         super().__init__(position=position, direction=direction, diameter=diameter, N=N,
-                         useHardwareAcceleration=useHardwareAcceleration, displaySize=displaySize)
+                         useHardwareAcceleration=useHardwareAcceleration, displaySize=displaySize, seed=seed)
 
     def _getInitialDirections(self):
         thetaDiameter = np.tan(self._divergence/2) * 2