Cycles: more closely match some math and intersection operations in Embree

This helps with debugging, and gives a slightly closer match between CPU
and CUDA/HIP/Metal renders when it comes to ray tracing precision.

1 files changed, 16 insertions, 19 deletions

diff --git a/intern/cycles/util/math_intersect.h b/intern/cycles/util/math_intersect.h
index c5b1cd51030..3e5891b2507 100644
--- a/intern/cycles/util/math_intersect.h
+++ b/intern/cycles/util/math_intersect.h
@@ -105,10 +105,10 @@ ccl_device bool ray_disk_intersect(float3 ray_P,
   return false;
 }
 
-ccl_device_forceinline bool ray_triangle_intersect(float3 ray_P,
-                                                   float3 ray_dir,
-                                                   float ray_tmin,
-                                                   float ray_tmax,
+ccl_device_forceinline bool ray_triangle_intersect(const float3 ray_P,
+                                                   const float3 ray_D,
+                                                   const float ray_tmin,
+                                                   const float ray_tmax,
                                                    const float3 tri_a,
                                                    const float3 tri_b,
                                                    const float3 tri_c,
@@ -116,14 +116,13 @@ ccl_device_forceinline bool ray_triangle_intersect(float3 ray_P,
                                                    ccl_private float *isect_v,
                                                    ccl_private float *isect_t)
 {
-#define dot3(a, b) dot(a, b)
-  const float3 P = ray_P;
-  const float3 dir = ray_dir;
+  /* This implementation matches the Plücker coordinates triangle intersection
+   * in Embree. */
 
   /* Calculate vertices relative to ray origin. */
-  const float3 v0 = tri_c - P;
-  const float3 v1 = tri_a - P;
-  const float3 v2 = tri_b - P;
+  const float3 v0 = tri_c - ray_P;
+  const float3 v1 = tri_a - ray_P;
+  const float3 v2 = tri_b - ray_P;
 
   /* Calculate triangle edges. */
   const float3 e0 = v2 - v0;
@@ -131,29 +130,29 @@ ccl_device_forceinline bool ray_triangle_intersect(float3 ray_P,
   const float3 e2 = v1 - v2;
 
   /* Perform edge tests. */
-  const float U = dot(cross(v2 + v0, e0), ray_dir);
-  const float V = dot(cross(v0 + v1, e1), ray_dir);
-  const float W = dot(cross(v1 + v2, e2), ray_dir);
+  const float U = dot(cross(v2 + v0, e0), ray_D);
+  const float V = dot(cross(v0 + v1, e1), ray_D);
+  const float W = dot(cross(v1 + v2, e2), ray_D);
 
+  const float eps = FLT_EPSILON * fabsf(U + V + W);
   const float minUVW = min(U, min(V, W));
   const float maxUVW = max(U, max(V, W));
 
-  if (minUVW < 0.0f && maxUVW > 0.0f) {
+  if (!(minUVW >= -eps || maxUVW <= eps)) {
     return false;
   }
 
   /* Calculate geometry normal and denominator. */
   const float3 Ng1 = cross(e1, e0);
-  // const Vec3vfM Ng1 = stable_triangle_normal(e2,e1,e0);
   const float3 Ng = Ng1 + Ng1;
-  const float den = dot3(Ng, dir);
+  const float den = dot(Ng, ray_D);
   /* Avoid division by 0. */
   if (UNLIKELY(den == 0.0f)) {
     return false;
   }
 
   /* Perform depth test. */
-  const float T = dot3(v0, Ng);
+  const float T = dot(v0, Ng);
   const float t = T / den;
   if (!(t >= ray_tmin && t <= ray_tmax)) {
     return false;
@@ -163,8 +162,6 @@ ccl_device_forceinline bool ray_triangle_intersect(float3 ray_P,
   *isect_v = V / den;
   *isect_t = t;
   return true;
-
-#undef dot3
 }
 
 /* Tests for an intersection between a ray and a quad defined by