Working basic triangle intersection.

rtiow/renderer/src/scenes/stltest.rs

@@ -11,13 +11,14 @@ use crate::{
     renderer::{Opt, Scene},
     sphere::Sphere,
     texture::ConstantTexture,
+    translate::Translate,
     triangles::Triangles,
     vec3::Vec3,
 };
 
 pub fn new(opt: &Opt) -> Scene {
-    let lookfrom = Vec3::new(-100., 200., 200.);
+    let lookfrom = Vec3::new(-20., 100., -100.);
-    let lookat = Vec3::new(0., 1., 0.);
+    let lookat = Vec3::new(0., 10., 0.);
     let dist_to_focus = 10.0;
     let aperture = 0.0;
     let time_min = 0.;
@@ -46,12 +47,6 @@ pub fn new(opt: &Opt) -> Scene {
     .expect("failed to parse cube");
     let light_size = 50.;
     let light_height = 200.;
-    let tri_mesh = Triangles::new(
-        &stl_cube,
-        Lambertian::new(ConstantTexture::new(Vec3::new(1., 1., 1.))),
-        1.,
-    );
-    dbg!(tri_mesh.triangles.len(), tri_mesh.bbox);
     let objects: Vec<Box<dyn Hit>> = vec![
         // Light from above
         Box::new(XZRect::new(
@@ -98,7 +93,15 @@ pub fn new(opt: &Opt) -> Scene {
             20.,
             Metal::new(Vec3::new(0.8, 0.8, 0.8), 0.2),
         )),
-        Box::new(tri_mesh),
+        // STL Mesh
+        Box::new(Translate::new(
+            Triangles::new(
+                &stl_cube,
+                Lambertian::new(ConstantTexture::new(Vec3::new(1., 1.0, 1.0))),
+                1.,
+            ),
+            [-10., 0., 0.],
+        )),
     ];
     let world: Box<dyn Hit> = if opt.use_accel {
         Box::new(KDTree::new(objects, time_min, time_max))

rtiow/renderer/src/triangles.rs

@@ -1,3 +1,5 @@
+use std::f32::EPSILON;
+
 use stl::STL;
 
 use crate::{
@@ -13,11 +15,6 @@ pub struct Triangle {
     normal: Vec3,
     verts: [Vec3; 3],
     // Precomputed data
-    // TODO(wathiede): precompute `d` on load.
-    d: f32,
-    edge0: Vec3,
-    edge1: Vec3,
-    edge2: Vec3,
 }
 
 #[derive(Debug)]
@@ -42,13 +39,17 @@ where
                 let v0 = t.verts[0] * scale_factor;
                 let v1 = t.verts[1] * scale_factor;
                 let v2 = t.verts[2] * scale_factor;
+                assert_eq!(
+                    t.normal,
+                    cross(v1 - v0, v2 - v0).unit_vector(),
+                    "v1 {} v2 {} v3 {}",
+                    v0,
+                    v1,
+                    v2
+                );
                 Triangle {
                     normal: t.normal,
                     verts: [v0, v1, v2],
-                    d: -dot(t.normal, t.verts[0]),
-                    edge0: v1 - v0,
-                    edge1: v2 - v1,
-                    edge2: v0 - v2,
                 }
             })
             .collect();
@@ -77,62 +78,172 @@ where
     }
 }
 
+/// Based on https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection.html
+fn ray_triangle_intersect_moller_trumbore1(r: Ray, tri: &Triangle) -> Option<RayTriangleResult> {
+    // #ifdef MOLLER_TRUMBORE
+    //     Vec3f v0v1 = v1 - v0;
+    //     Vec3f v0v2 = v2 - v0;
+    //     Vec3f pvec = dir.crossProduct(v0v2);
+    //     float det = v0v1.dotProduct(pvec);
+    // #ifdef CULLING
+    //     // if the determinant is negative, the triangle is 'back facing'
+    //     // if the determinant is close to 0, the ray misses the triangle
+    //     if (det < kEpsilon) return false;
+    // #else
+    //     // ray and triangle are parallel if det is close to 0
+    //     if (fabs(det) < kEpsilon) return false;
+    // #endif
+    //     float invDet = 1 / det;
+    //
+    //     Vec3f tvec = orig - v0;
+    //     u = tvec.dotProduct(pvec) * invDet;
+    //     if (u < 0 || u > 1) return false;
+    //
+    //     Vec3f qvec = tvec.crossProduct(v0v1);
+    //     v = dir.dotProduct(qvec) * invDet;
+    //     if (v < 0 || u + v > 1) return false;
+    //
+    //     t = v0v2.dotProduct(qvec) * invDet;
+    //
+    let v0 = tri.verts[0];
+    let v1 = tri.verts[1];
+    let v2 = tri.verts[2];
+
+    let v0v1 = v1 - v0;
+    let v0v2 = v2 - v0;
+    let p = cross(r.direction, v0v2);
+    let det = dot(v0v1, p);
+    if det < EPSILON {
+        return None;
+    }
+
+    let inv_det = 1. / det;
+
+    let t = r.origin - v0;
+    let u = dot(t, p) * inv_det;
+    if u < 0. || u > 1. {
+        return None;
+    }
+
+    let q = cross(t, v0v1);
+    let v = dot(r.direction, q) * inv_det;
+    if v < 0. || u + v > 1. {
+        return None;
+    }
+    let t = dot(v0v2, q) * inv_det;
+
+    if t > EPSILON {
+        return Some(RayTriangleResult {
+            t,
+            p: r.origin + r.direction * t,
+        });
+    }
+    None
+}
+
+// From https://en.wikipedia.org/wiki/M%C3%B6ller%E2%80%93Trumbore_intersection_algorithm
+fn ray_triangle_intersect_moller_trumbore2(r: Ray, tri: &Triangle) -> Option<RayTriangleResult> {
+    let v0 = tri.verts[0];
+    let v1 = tri.verts[1];
+    let v2 = tri.verts[2];
+
+    let edge1 = v1 - v0;
+    let edge2 = v2 - v0;
+    let h = cross(r.direction, edge2);
+    let a = dot(edge1, h);
+    if a < EPSILON {
+        return None;
+    }
+
+    let f = 1. / a;
+    let s = r.origin - v0;
+    let u = f * dot(s, h);
+    if u < 0. || u > 1. {
+        return None;
+    }
+    let q = cross(s, edge1);
+    let v = f * dot(r.direction, q);
+    if v < 0. || u + v > 1. {
+        return None;
+    }
+    // At this stage we can compute t to find out where the intersection point is on the line.
+    let t = f * dot(edge2, q);
+    // ray intersection
+    if t > EPSILON {
+        return Some(RayTriangleResult {
+            t,
+            p: r.origin + r.direction * t,
+        });
+    }
+    // This means that there is a line intersection but not a ray intersection.
+    None
+}
+/// Based on https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/ray-triangle-intersection-geometric-solution.html
+fn ray_triangle_intersect_geometric(r: Ray, tri: &Triangle) -> Option<RayTriangleResult> {
+    let n = tri.normal;
+
+    // close enough to parallel, won't hit
+    let n_dot_dir = dot(n, r.direction);
+    if n_dot_dir < EPSILON {
+        return None;
+    }
+
+    let v0 = tri.verts[0];
+    let v1 = tri.verts[1];
+    let v2 = tri.verts[2];
+    let d = -dot(n, v0);
+    let t = -(dot(n, r.origin) + d) / n_dot_dir;
+    // check if the triangle is behind the ray
+    if t < 0. {
+        return None;
+    }
+
+    let p = r.origin + t * r.direction;
+
+    let edge0 = v1 - v0;
+    let edge1 = v2 - v1;
+    let edge2 = v0 - v2;
+    let vp0 = p - v0;
+    let c = cross(edge0, vp0);
+    if dot(n, c) < 0. {
+        return None;
+    }
+
+    let vp1 = p - v1;
+    let c = cross(edge1, vp1);
+    if dot(n, c) < 0. {
+        return None;
+    }
+
+    let vp2 = p - v2;
+    let c = cross(edge2, vp2);
+    if dot(n, c) < 0. {
+        return None;
+    }
+    Some(RayTriangleResult { t, p })
+}
+
 impl<M> Hit for Triangles<M>
 where
     M: Material,
 {
-    // Based on https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/ray-triangle-intersection-geometric-solution.html
     fn hit(&self, r: Ray, _t_min: f32, _t_max: f32) -> Option<HitRecord> {
         // TODO(wathiede): add an acceleration structure to more cheaply skip some triangles.
         for tri in &self.triangles {
-            let n = tri.normal;
+            if let Some(RayTriangleResult { t, p }) =
-
+                ray_triangle_intersect_moller_trumbore1(r, tri)
-            // close enough to parallel, won't hit
+            {
-            // TODO(wathiede): verify EPSILON
+                //if let Some(RayTriangleResult { t, p }) = ray_triangle_intersect_geometric(r, tri) {
-            const EPSILON: f32 = 0.00001;
+                // We don't support UV (yet?).
-            let n_dot_dir = dot(n, r.direction);
+                let uv = (0.5, 0.5);
-            if n_dot_dir < EPSILON {
+                return Some(HitRecord {
-                continue;
+                    t,
+                    uv,
+                    p,
+                    normal: tri.normal,
+                    material: &self.material,
+                });
             }
-
-            let t = -(dot(n, r.origin) + tri.d) / n_dot_dir;
-            // check if the triangle is behind the ray
-            if t < 0. {
-                continue;
-            }
-
-            let p = r.origin + t * r.direction;
-            let v0 = tri.verts[0];
-            let v1 = tri.verts[1];
-            let v2 = tri.verts[2];
-
-            let vp0 = p - v0;
-            let c = cross(tri.edge0, vp0);
-            if dot(n, c) < 0. {
-                continue;
-            }
-
-            let vp1 = p - v1;
-            let c = cross(tri.edge1, vp1);
-            if dot(n, c) < 0. {
-                continue;
-            }
-
-            let vp2 = p - v2;
-            let c = cross(tri.edge2, vp2);
-            if dot(n, c) < 0. {
-                continue;
-            }
-
-            // We don't support UV (yet?).
-            let uv = (0., 0.);
-            return Some(HitRecord {
-                t,
-                uv,
-                p,
-                normal: n,
-                material: &self.material,
-            });
         }
         None
     }
@@ -141,3 +252,8 @@ where
         Some(self.bbox)
     }
 }
+
+struct RayTriangleResult {
+    t: f32,
+    p: Vec3,
+}