rtiow: better debugging, testing and fix some BVHTriangles bugs.

rtiow/renderer/src/bvh_triangles.rs

@@ -1,6 +1,7 @@
 /// Implementation based on blog post @
 /// https://jacco.ompf2.com/2022/04/13/how-to-build-a-bvh-part-1-basics/
 use std::f32::EPSILON;
+use std::fmt;
 
 use stl::STL;
 
@@ -12,7 +13,7 @@ use crate::{
     vec3::{cross, dot, Vec3},
 };
 
-#[derive(Debug)]
+#[derive(Debug, PartialEq)]
 struct BVHNode {
     aabb: AABB,
     left_child: usize,
@@ -26,13 +27,21 @@ impl BVHNode {
     }
 }
 
-#[derive(Debug)]
+#[derive(PartialEq)]
 pub struct Triangle {
     centroid: Vec3,
     verts: [Vec3; 3],
 }
+impl fmt::Debug for Triangle {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(
+            f,
+            "Tri: <{}, {}, {}> @ {}",
+            self.verts[0], self.verts[1], self.verts[2], self.centroid
+        )
+    }
+}
 
-#[derive(Debug)]
 pub struct BVHTriangles<M>
 where
     M: Material,
@@ -42,12 +51,56 @@ where
     bvh_nodes: Vec<BVHNode>,
 }
 
+impl<M> fmt::Debug for BVHTriangles<M>
+where
+    M: Material,
+{
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(f, "{} triangles", self.triangles.len())?;
+        if f.alternate() {
+            writeln!(f)?;
+        }
+        for (i, t) in self.triangles.iter().enumerate() {
+            if f.alternate() {
+                write!(f, "\t")?;
+            }
+            write!(f, "{i:>3} {:?} ", t)?;
+            if f.alternate() {
+                writeln!(f)?;
+            }
+        }
+        if f.alternate() {
+            writeln!(f)?;
+        }
+        for (i, n) in self.bvh_nodes.iter().enumerate() {
+            write!(f, "N[{i}] {n:?}")?;
+            if f.alternate() {
+                writeln!(f)?;
+            }
+            let n = &self.bvh_nodes[i];
+            if n.is_leaf() {
+                for t_idx in n.first_prim..(n.first_prim + n.prim_count) {
+                    if f.alternate() {
+                        write!(f, "\t")?;
+                    }
+                    write!(f, "{:?} ", self.triangles[t_idx])?;
+                    if f.alternate() {
+                        writeln!(f)?;
+                    }
+                }
+            }
+        }
+        Ok(())
+    }
+}
+
 const ROOT_NODE_IDX: usize = 0;
 impl<M> BVHTriangles<M>
 where
     M: Material,
 {
     pub fn new(stl: &STL, material: M) -> BVHTriangles<M> {
+        let div3 = 1. / 3.;
         let triangles: Vec<_> = stl
             .triangles
             .iter()
@@ -55,7 +108,7 @@ where
                 let v0 = t.verts[0];
                 let v1 = t.verts[1];
                 let v2 = t.verts[2];
-                let centroid = (v0 + v1 + v2) * 0.3333;
+                let centroid = (v0 + v1 + v2) * div3;
                 Triangle {
                     centroid,
                     verts: [v0, v1, v2],
@@ -80,7 +133,7 @@ where
             aabb: AABB::default(),
             left_child: 0,
             first_prim: 0,
-            prim_count: self.triangles.len() - 1,
+            prim_count: self.triangles.len(),
         };
         self.bvh_nodes.push(root);
         self.update_node_bounds(ROOT_NODE_IDX);
@@ -91,7 +144,7 @@ where
         let node = &mut self.bvh_nodes[node_idx];
         let mut aabb_min: Vec3 = f32::MAX.into();
         let mut aabb_max: Vec3 = f32::MIN.into();
-        for i in node.first_prim..node.prim_count {
+        for i in node.first_prim..(node.first_prim + node.prim_count) {
             let leaf_tri = &self.triangles[i];
             aabb_min = vec3::min(aabb_min, leaf_tri.verts[0]);
             aabb_min = vec3::min(aabb_min, leaf_tri.verts[1]);
@@ -111,7 +164,7 @@ where
             let node = &self.bvh_nodes[idx];
 
             // Compute split plane and position.
-            let extent = node.aabb.min() - node.aabb.max();
+            let extent = node.aabb.max() - node.aabb.min();
             let axis = node.aabb.longest_axis();
             let split_pos = node.aabb.min()[axis] + extent[axis] * 0.5;
 
@@ -168,31 +221,37 @@ where
         if !node.aabb.hit(r, t_min, t_max) {
             return None;
         }
+        //dbg!(&self);
 
         if node.is_leaf() {
-            for tri in &self.triangles {
+            return self
-                if let Some(RayTriangleResult { t, p }) =
+                .triangles
-                    ray_triangle_intersect_moller_trumbore(r, tri)
+                .iter()
-                {
+                .map(|tri| {
-                    //if let Some(RayTriangleResult { t, p }) = ray_triangle_intersect_geometric(r, tri) {
+                    if let Some(RayTriangleResult { t, p }) = intersect_tri(r, tri) {
-                    // We don't support UV (yet?).
+                        // We don't support UV (yet?).
-                    let uv = (0.5, 0.5);
+                        let uv = (0.5, 0.5);
-                    let v0 = tri.verts[0];
+                        let v0 = tri.verts[0];
-                    let v1 = tri.verts[1];
+                        let v1 = tri.verts[1];
-                    let v2 = tri.verts[2];
+                        let v2 = tri.verts[2];
 
-                    let v0v1 = v1 - v0;
+                        let v0v1 = v1 - v0;
-                    let v0v2 = v2 - v0;
+                        let v0v2 = v2 - v0;
-                    let normal = cross(v0v1, v0v2).unit_vector();
+                        let normal = cross(v0v1, v0v2).unit_vector();
-                    return Some(HitRecord {
+                        //println!("hit triangle {tri:?}");
-                        t,
+                        Some(HitRecord {
-                        uv,
+                            t,
-                        p,
+                            uv,
-                        normal,
+                            p,
-                        material: &self.material,
+                            normal,
-                    });
+                            material: &self.material,
-                }
+                        })
-            }
+                    } else {
+                        None
+                    }
+                })
+                .filter_map(|hr| hr)
+                .min_by(|a, b| a.t.partial_cmp(&b.t).unwrap());
         } else {
             let r1 = self.intersect_bvh(r, node.left_child, t_min, t_max);
             let r2 = self.intersect_bvh(r, node.left_child + 1, t_min, t_max);
@@ -209,26 +268,8 @@ where
     }
 }
 
-impl<M> Hit for BVHTriangles<M>
-where
-    M: Material,
-{
-    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
-        self.intersect_bvh(r, 0, t_min, t_max)
-    }
-
-    fn bounding_box(&self, _t_min: f32, _t_max: f32) -> Option<AABB> {
-        Some(self.bvh_nodes[0].aabb)
-    }
-}
-
-struct RayTriangleResult {
-    t: f32,
-    p: Vec3,
-}
-///
 /// 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_trumbore(r: Ray, tri: &Triangle) -> Option<RayTriangleResult> {
+fn intersect_tri(r: Ray, tri: &Triangle) -> Option<RayTriangleResult> {
     // #ifdef MOLLER_TRUMBORE
     //     Vec3f v0v1 = v1 - v0;
     //     Vec3f v0v2 = v2 - v0;
@@ -262,7 +303,7 @@ fn ray_triangle_intersect_moller_trumbore(r: Ray, tri: &Triangle) -> Option<RayT
     let v0v2 = v2 - v0;
     let p = cross(r.direction, v0v2);
     let det = dot(v0v1, p);
-    if det < EPSILON {
+    if det.abs() < EPSILON {
         return None;
     }
 
@@ -290,17 +331,97 @@ fn ray_triangle_intersect_moller_trumbore(r: Ray, tri: &Triangle) -> Option<RayT
     None
 }
 
+impl<M> Hit for BVHTriangles<M>
+where
+    M: Material,
+{
+    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        self.intersect_bvh(r, 0, t_min, t_max)
+    }
+
+    fn bounding_box(&self, _t_min: f32, _t_max: f32) -> Option<AABB> {
+        Some(self.bvh_nodes[0].aabb)
+    }
+}
+
+struct RayTriangleResult {
+    t: f32,
+    p: Vec3,
+}
+
 #[cfg(test)]
 mod tests {
+    use super::*;
     use crate::{
-        bvh_triangles::BVHTriangles, cuboid::Cuboid, hitable::Hit, material::Dielectric, ray::Ray,
+        bvh_triangles::BVHTriangles,
+        cuboid::Cuboid,
+        hitable::Hit,
+        material::{Dielectric, Lambertian},
+        ray::Ray,
+        texture::ConstantTexture,
     };
+    //use pretty_assertions::assert_eq;
     use std::{
         io::{BufReader, Cursor},
         sync::Arc,
     };
     use stl::STL;
 
+    /*
+        #[test]
+        fn build_bvh() {
+            let stl_triangles: Vec<_> = (0..4)
+                .flat_map(|y| {
+                    (0..2).map(move |x| {
+                        let x = x as f32;
+                        let y = y as f32;
+                        stl::Triangle {
+                            normal: [1., 0., 0.].into(),
+                            verts: [
+                                [2. * x + 0., 2. * y + 0., 0.].into(),
+                                [2. * x + 1., 2. * y + 0., 0.].into(),
+                                [2. * x + 1., 2. * y + 1., 0.].into(),
+                            ],
+                            attr: 0,
+                        }
+                    })
+                })
+                .collect();
+            let stl = STL {
+                header: [0; 80],
+                triangles: stl_triangles,
+            };
+            /*
+            let mut bvh_triangles: Vec<_> = stl_triangles
+            .iter()
+            .map(|tri| {
+            let div3 = 1. / 3.;
+            let v0 = tri.verts[0];
+            let v1 = tri.verts[1];
+            let v2 = tri.verts[2];
+            let centroid = (v0 + v1 + v2) * div3;
+            Triangle {
+            centroid,
+            verts: tri.verts,
+            }
+            })
+            .collect();
+            bvh_triangles.sort_by(|a, b| a.centroid.y.partial_cmp(&b.centroid.y).unwrap());
+            let material = Lambertian::new(ConstantTexture::new([0., 0., 0.]));
+            let bvh_nodes = Default::default();
+            let want = BVHTriangles {
+            triangles: bvh_triangles,
+            bvh_nodes,
+            material,
+            };
+            */
+            let material = Lambertian::new(ConstantTexture::new([0., 0., 0.]));
+            let bvh = BVHTriangles::new(&stl, material);
+            dbg!(&bvh);
+            assert_eq!(bvh.bvh_nodes.len(), 2 * bvh.triangles.len() - 2);
+        }
+    */
+
     #[test]
     fn compare_cuboid() {
         let c = Cuboid::new(
@@ -320,7 +441,8 @@ mod tests {
             //Metal::new(Vec3::new(0.8, 0.8, 0.8), 0.2),
             //Lambertian::new(ConstantTexture::new(Vec3::new(1.0, 0.2, 0.2))),
         );
-        let rays: Vec<_> = (1..20)
+        //dbg!(&s);
+        let mut rays: Vec<_> = (1..20)
             .flat_map(|y| {
                 (1..20).flat_map(move |x| {
                     let x = x as f32;
@@ -337,27 +459,28 @@ mod tests {
                 })
             })
             .collect();
+        rays.push(Ray::new([10., 10., 10.], [1., 0., 0.], 0.));
 
-        for (i, r) in rays.into_iter().enumerate() {
+        for r in rays.into_iter() {
             let c_hit = c
                 .hit(r, 0., f32::MAX)
-                .expect(&format!("c_hit missed {i}: {r:?}"));
+                .expect(&format!("c_hit missed {r:#?}"));
             let s_hit = s
                 .hit(r, 0., f32::MAX)
-                .expect(&format!("s_hit missed {i}: {r:?}"));
+                .expect(&format!("s_hit missed {r:#?}"));
             assert_eq!(
                 c_hit.t, s_hit.t,
-                "{i}: [t] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
+                "{r:?} [t] c_hit: {c_hit:#?}, s_hit: {s_hit:#?}"
             );
             // uv isn't valid for BVHTriangles.
             // assert_eq!( c_hit.uv, s_hit.uv, "{i}: [uv] c_hit: {c_hit:?}, s_hit: {s_hit:?}");
             assert_eq!(
                 c_hit.p, s_hit.p,
-                "{i}: [p] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
+                "{r:?}: [p] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
             );
             assert_eq!(
                 c_hit.normal, s_hit.normal,
-                "{i}: [normal] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
+                "{r:?}: [normal] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
             );
         }
     }