From 71bdf203aa208bd738fcdbb151a1897b916da672 Mon Sep 17 00:00:00 2001
From: Nathan Vegdahl <cessen@cessen.com>
Date: Sun, 18 Jun 2017 20:51:53 -0700
Subject: [PATCH] Fixed bug in mesh intersection code.

Very small triangles were being missed because of the
not-so-robust ray-triangle intersection algorithm I was using.

Switched to the algorithm from the paper "Watertight
Ray/Triangle Intersection" by Woop et al.  Happily, the new
algorithm doesn't seem to measurably slow down renders at all.
---
 src/surface/triangle.rs      | 131 +++++++++++++++++++++++++++--------
 src/surface/triangle_mesh.rs |   2 +-
 2 files changed, 103 insertions(+), 30 deletions(-)

diff --git a/src/surface/triangle.rs b/src/surface/triangle.rs
index e935127..0a776eb 100644
--- a/src/surface/triangle.rs
+++ b/src/surface/triangle.rs
@@ -1,40 +1,113 @@
 #![allow(dead_code)]
 
-use std;
-
-use math::{Point, cross, dot};
+use math::Point;
 use ray::Ray;
 
 
-/// Intersects ray with tri, returning (t, u, v), or None if no intersection.
-pub fn intersect_ray(ray: &Ray, tri: (Point, Point, Point)) -> Option<(f32, f32, f32)> {
-    let edge1 = tri.1 - tri.0;
-    let edge2 = tri.2 - tri.0;
-    let pvec = cross(ray.dir, edge2);
-    let det = dot(edge1, pvec);
+/// Intersects `ray` with `tri`, returning `Some((t, b0, b1, b2))`, or `None`
+/// if no intersection.
+///
+/// Returned values:
+///
+/// * `t` is the ray t at the hit point.
+/// * `b0`, `b1`, and `b2` are the barycentric coordinates of the triangle at
+///   the hit point.
+///
+/// Uses the ray-triangle test from the paper "Watertight Ray/Triangle
+/// Intersection" by Woop et al.
+pub fn intersect_ray(ray: &Ray, tri: (Point, Point, Point)) -> Option<(f32, f32, f32, f32)> {
+    // Calculate the permuted dimension indices for the new ray space.
+    let zi = {
+        let xabs = ray.dir.x().abs();
+        let yabs = ray.dir.y().abs();
+        let zabs = ray.dir.z().abs();
 
-    if det <= -std::f32::EPSILON || det >= std::f32::EPSILON {
-        let inv_det = 1.0 / det;
-        let tvec = ray.orig - tri.0;
-        let qvec = cross(tvec, edge1);
-
-        let t = dot(edge2, qvec) * inv_det;
-        if t < 0.0 {
-            return None;
+        if xabs > yabs && xabs > zabs {
+            0
+        } else if yabs > zabs {
+            1
+        } else {
+            2
         }
-
-        let u = dot(tvec, pvec) * inv_det;
-        if u < 0.0 || u > 1.0 {
-            return None;
-        }
-
-        let v = dot(ray.dir, qvec) * inv_det;
-        if v < 0.0 || (u + v) > 1.0 {
-            return None;
-        }
-
-        return Some((t, u, v));
+    };
+    let (xi, yi) = if ray.dir.get_n(zi) >= 0.0 {
+        let xi = if zi == 2 { 0 } else { zi + 1 };
+        let yi = if xi == 2 { 0 } else { xi + 1 };
+        (xi, yi)
     } else {
+        let xi = if zi == 2 { 0 } else { zi + 1 };
+        let yi = if xi == 2 { 0 } else { xi + 1 };
+        (yi, xi)
+    };
+
+    let dir_x = ray.dir.get_n(xi);
+    let dir_y = ray.dir.get_n(yi);
+    let dir_z = ray.dir.get_n(zi);
+
+    // Calculate shear constants.
+    let sx = dir_x / dir_z;
+    let sy = dir_y / dir_z;
+    let sz = 1.0 / dir_z;
+
+    // Calculate vertices in ray space.
+    let p0 = tri.0 - ray.orig;
+    let p1 = tri.1 - ray.orig;
+    let p2 = tri.2 - ray.orig;
+
+    let p0x = p0.get_n(xi) - sx * p0.get_n(zi);
+    let p0y = p0.get_n(yi) - sy * p0.get_n(zi);
+    let p1x = p1.get_n(xi) - sx * p1.get_n(zi);
+    let p1y = p1.get_n(yi) - sy * p1.get_n(zi);
+    let p2x = p2.get_n(xi) - sx * p2.get_n(zi);
+    let p2y = p2.get_n(yi) - sy * p2.get_n(zi);
+
+    // Calculate scaled barycentric coordinates.
+    let mut e0 = (p2x * p1y) - (p2y * p1x);
+    let mut e1 = (p0x * p2y) - (p0y * p2x);
+    let mut e2 = (p1x * p0y) - (p1y * p0x);
+
+    // Fallback to test against edges using double precision.
+    if e0 == 0.0 || e1 == 0.0 || e2 == 0.0 {
+        let p2xp1y = p2x as f64 * p1y as f64;
+        let p2yp1x = p2y as f64 * p1x as f64;
+        e0 = (p2xp1y - p2yp1x) as f32;
+        let p0xp2y = p0x as f64 * p2y as f64;
+        let p0yp2x = p0y as f64 * p2x as f64;
+        e1 = (p0xp2y - p0yp2x) as f32;
+        let p1xp0y = p1x as f64 * p0y as f64;
+        let p1yp0x = p1y as f64 * p0x as f64;
+        e2 = (p1xp0y - p1yp0x) as f32;
+    }
+
+    // Check if the ray hit the triangle.
+    if (e0 < 0.0 || e1 < 0.0 || e2 < 0.0) && (e0 > 0.0 || e1 > 0.0 || e2 > 0.0) {
         return None;
     }
+
+    // Determinant
+    let det = e0 + e1 + e2;
+    if det == 0.0 {
+        return None;
+    }
+
+    // Calculate t of hitpoint.
+    let p0z = sz * p0.get_n(zi);
+    let p1z = sz * p1.get_n(zi);
+    let p2z = sz * p2.get_n(zi);
+    let t = (e0 * p0z) + (e1 * p1z) + (e2 * p2z);
+
+    // Check if the hitpoint t is within ray min/max t.
+    if det >= 0.0 {
+        if t <= 0.0 || t > (ray.max_t * det) {
+            return None;
+        }
+    } else {
+        if -t <= 0.0 || -t > (ray.max_t * -det) {
+            return None;
+        }
+    }
+
+    // Return t and the hitpoint barycentric coordinates.
+    let inv_det = 1.0 / det;
+    Some((t * inv_det, e0 * inv_det, e1 * inv_det, e2 * inv_det))
 }
diff --git a/src/surface/triangle_mesh.rs b/src/surface/triangle_mesh.rs
index aefe708..219693f 100644
--- a/src/surface/triangle_mesh.rs
+++ b/src/surface/triangle_mesh.rs
@@ -92,7 +92,7 @@ impl<'a> Surface for TriangleMesh<'a> {
                         };
                         let mat_inv = mat_space.inverse();
                         let tri = (tri.0 * mat_inv, tri.1 * mat_inv, tri.2 * mat_inv);
-                        if let Some((t, _, _)) = triangle::intersect_ray(wr, tri) {
+                        if let Some((t, _, _, _)) = triangle::intersect_ray(wr, tri) {
                             if t < r.max_t {
                                 if r.is_occlusion() {
                                     isects[r.id as usize] = SurfaceIntersection::Occlude;