Fix intersection mechanic around floating-point limit.

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
 

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
 

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;
@@ -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;

pytissueoptics/rayscattering/opencl/src/intersection.c

@@ -1,17 +1,19 @@
-__constant float EPS = 0.00001f;
-__constant float EPS_CORRECTION = 0.0005f;
-__constant float EPS_PARALLEL = 0.00001f;
-__constant float EPS_SIDE = 0.000001f;
+__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;
 
 struct Intersection {
     uint exists;
-    uint isTooClose;
     float distance;
     float3 position;
     float3 normal;
     uint surfaceID;
     uint polygonID;
     float distanceLeft;
+    bool isSmooth;
+    float3 rawNormal;
 };
 
 typedef struct Intersection Intersection;
@@ -113,11 +115,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) {
@@ -149,16 +157,15 @@ void _sortSolidCandidates(Scene *scene, uint gid) {
 
 struct HitPoint {
     bool exists;
-    bool isTooClose;
+    float distance;
     float3 position;
 };
 
 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;
 
     float3 edgeA = v2 - v1;
     float3 edgeB = v3 - v1;
@@ -179,56 +186,119 @@ HitPoint _getTriangleIntersection(Ray ray, float3 v1, float3 v2, float3 v3) {
 
     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;
     }
 
     float t = dot(edgeB, qVector) * invDet;
+    hitPoint.distance = t;
+    hitPoint.position = ray.origin + t * ray.direction;
+
+    // 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.0f) {
+    if (t >= 0 && ray.length >= t){
+        hitPoint.exists = true;
         return hitPoint;
     }
 
-    if (t > (ray.length + EPS)) {
-        // No Intersection, it's too far away
+    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) {
+        // Forward catch.
+        hitPoint.exists = true;
+        return hitPoint;
+    }
+
+    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;
-    } 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;
+
+    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);
+            HitPoint hitPoint = _getTriangleIntersection(ray, vertices[vertexIDs[0]].position, vertices[vertexIDs[1]].position, vertices[vertexIDs[2]].position, triangles[p].normal);
+
             if (!hitPoint.exists) {
                 continue;
             }
-            float distance = length(hitPoint.position - ray.origin);
-            if (distance < intersection.distance) {
+
+            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;
+            }
+
+            if (fabs(hitPoint.distance) < fabs(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;
                 intersection.polygonID = p;
             }
         }
     }
+
+    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;
 }
 
@@ -284,34 +354,37 @@ 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) {
     if (!intersection->exists) {
         return;
     }
 
+    intersection->isSmooth = false;
     if (scene->surfaces[intersection->surfaceID].toSmooth) {
         setSmoothNormal(intersection, scene->triangles, scene->vertices, ray);
     }
     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;
     closestIntersection.exists = false;
-    closestIntersection.isTooClose = false;
     closestIntersection.distance = INFINITY;
     if (scene->nSolids == 0) {
         return closestIntersection;
@@ -323,14 +396,16 @@ Intersection findIntersection(Ray ray, Scene *scene, uint gid) {
             // 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;
         }
 
         uint solidID = scene->solidCandidates[boxGID].solidID;
-        Intersection intersection = _findClosestPolygonIntersection(ray, solidID, scene->solids, scene->surfaces, scene->triangles, scene->vertices);
-        if (intersection.exists  && intersection.distance < closestIntersection.distance) {
+        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 +420,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);
 }