rtiow: add commented out failing test.

This seems to be faster than using both HT buddies.

rtiow/Cargo.lock

@@ -812,6 +812,16 @@ dependencies = [
  "subtle",
 ]
 
+[[package]]
+name = "ctor"
+version = "0.1.26"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6d2301688392eb071b0bf1a37be05c469d3cc4dbbd95df672fe28ab021e6a096"
+dependencies = [
+ "quote 1.0.23",
+ "syn 1.0.107",
+]
+
 [[package]]
 name = "ctr"
 version = "0.6.0"
@@ -885,6 +895,12 @@ dependencies = [
  "byteorder",
 ]
 
+[[package]]
+name = "diff"
+version = "0.1.13"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "56254986775e3233ffa9c4d7d3faaf6d36a2c09d30b20687e9f88bc8bafc16c8"
+
 [[package]]
 name = "digest"
 version = "0.8.1"
@@ -1910,6 +1926,15 @@ version = "6.4.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "9b7820b9daea5457c9f21c69448905d723fbd21136ccf521748f23fd49e723ee"
 
+[[package]]
+name = "output_vt100"
+version = "0.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "628223faebab4e3e40667ee0b2336d34a5b960ff60ea743ddfdbcf7770bcfb66"
+dependencies = [
+ "winapi 0.3.9",
+]
+
 [[package]]
 name = "owning_ref"
 version = "0.4.1"
@@ -2076,6 +2101,18 @@ version = "0.2.17"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de"
 
+[[package]]
+name = "pretty_assertions"
+version = "1.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a25e9bcb20aa780fd0bb16b72403a9064d6b3f22f026946029acb941a50af755"
+dependencies = [
+ "ctor",
+ "diff",
+ "output_vt100",
+ "yansi",
+]
+
 [[package]]
 name = "proc-macro-error"
 version = "1.0.4"
@@ -2434,6 +2471,7 @@ dependencies = [
  "lazy_static 1.4.0",
  "log 0.4.17",
  "num_cpus",
+ "pretty_assertions",
  "rand 0.8.5",
  "serde",
  "serde_derive",
@@ -3695,3 +3733,9 @@ dependencies = [
  "winapi 0.2.8",
  "winapi-build",
 ]
+
+[[package]]
+name = "yansi"
+version = "0.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "09041cd90cf85f7f8b2df60c646f853b7f535ce68f85244eb6731cf89fa498ec"

rtiow/renderer/Cargo.toml

@@ -33,6 +33,7 @@ strum_macros = "0.24.3"
 
 [dev-dependencies]
 criterion = "0.4"
+pretty_assertions = "1.3.0"
 
 [features]
 profile = ["cpuprofiler"]

rtiow/renderer/src/aabb.rs

@@ -2,7 +2,7 @@ use std::fmt;
 
 use crate::{ray::Ray, vec3::Vec3};
 
-#[derive(Default, Debug, Copy, Clone, PartialEq)]
+#[derive(Default, Copy, Clone, PartialEq)]
 pub struct AABB {
     bounds: [Vec3; 2],
 }
@@ -29,13 +29,25 @@ fn max(x: f32, y: f32) -> f32 {
     }
 }
 
+impl fmt::Debug for AABB {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(
+            f,
+            "AABB: <{} - {}> Vol: {}",
+            self.bounds[0],
+            self.bounds[1],
+            self.volume()
+        )
+    }
+}
+
 impl AABB {
     pub fn new<V: Into<Vec3>>(min: V, max: V) -> AABB {
         let min: Vec3 = min.into();
         let max: Vec3 = max.into();
-        assert!(min.x < max.x);
-        assert!(min.y < max.y);
-        assert!(min.z < max.z);
+        assert!(min.x <= max.x);
+        assert!(min.y <= max.y);
+        assert!(min.z <= max.z);
         AABB { bounds: [min, max] }
     }
 
@@ -167,6 +179,9 @@ impl AABB {
             let tmin = max(max(vmin4.0, max(vmin4.1, vmin4.2)), t_min);
             tmin <= tmax
         }
+        #[cfg(target_arch = "aarch64")]
+        // TODO(wathiede): add NEON implementation.
+        self.hit2(r, t_min, t_max)
     }
 
     pub fn hit_fast(&self, r: Ray, _t_min: f32, _t_max: f32) -> bool {

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,13 +221,14 @@ where
         if !node.aabb.hit(r, t_min, t_max) {
             return None;
         }
+        //dbg!(&self);
 
         if node.is_leaf() {
-            for tri in &self.triangles {
-                if let Some(RayTriangleResult { t, p }) =
-                    ray_triangle_intersect_moller_trumbore(r, tri)
-                {
-                    //if let Some(RayTriangleResult { t, p }) = ray_triangle_intersect_geometric(r, tri) {
+            return self
+                .triangles
+                .iter()
+                .map(|tri| {
+                    if let Some(RayTriangleResult { t, p }) = intersect_tri(r, tri) {
                         // We don't support UV (yet?).
                         let uv = (0.5, 0.5);
                         let v0 = tri.verts[0];
@@ -184,15 +238,20 @@ where
                         let v0v1 = v1 - v0;
                         let v0v2 = v2 - v0;
                         let normal = cross(v0v1, v0v2).unit_vector();
-                    return Some(HitRecord {
+                        //println!("hit triangle {tri:?}");
+                        Some(HitRecord {
                             t,
                             uv,
                             p,
                             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,32 +441,78 @@ 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 = [
-            Ray::new([-1., 1., 1.], [1., 0., 0.], 0.),
-            Ray::new([21., 1., 1.], [-1., 0., 0.], 0.),
-            Ray::new([1., -1., 1.], [0., 1., 0.], 0.),
-            Ray::new([1., 21., 1.], [0., -1., 0.], 0.),
-            Ray::new([1., 1., -1.], [0., 0., 1.], 0.),
-            Ray::new([1., 1., 21.], [0., 0., -1.], 0.),
-            // TODO more
-        ];
+        //dbg!(&s);
+        let mut rays: Vec<_> = (1..20)
+            .flat_map(|y| {
+                (1..20).flat_map(move |x| {
+                    let x = x as f32;
+                    let y = y as f32;
+                    vec![
+                        // Outward in angle
+                        Ray::new([-1., x, y], [1., 0., 0.], 0.),
+                        Ray::new([21., x, y], [-1., 0., 0.], 0.),
+                        Ray::new([x, -1., y], [0., 1., 0.], 0.),
+                        Ray::new([x, 21., y], [0., -1., 0.], 0.),
+                        Ray::new([x, y, -1.], [0., 0., 1.], 0.),
+                        Ray::new([x, y, 21.], [0., 0., -1.], 0.),
+                        // Inward out (
+                        Ray::new([x, y, 10.], [1., 0., 0.], 0.),
+                        Ray::new([x, y, 10.], [-1., 0., 0.], 0.),
+                        Ray::new([x, 10., y], [1., 0., 0.], 0.),
+                        Ray::new([x, 10., y], [-1., 0., 0.], 0.),
+                        Ray::new([10., x, y], [1., 0., 0.], 0.),
+                        Ray::new([10., x, y], [-1., 0., 0.], 0.),
+                        Ray::new([x, y, 10.], [0., 1., 0.], 0.),
+                        Ray::new([x, y, 10.], [0., -1., 0.], 0.),
+                        Ray::new([x, 10., y], [0., 1., 0.], 0.),
+                        Ray::new([x, 10., y], [0., -1., 0.], 0.),
+                        Ray::new([10., x, y], [0., 1., 0.], 0.),
+                        Ray::new([10., x, y], [0., -1., 0.], 0.),
+                        Ray::new([x, y, 10.], [0., 0., 1.], 0.),
+                        Ray::new([x, y, 10.], [0., 0., -1.], 0.),
+                        Ray::new([x, 10., y], [0., 0., 1.], 0.),
+                        Ray::new([x, 10., y], [0., 0., -1.], 0.),
+                        Ray::new([10., x, y], [0., 0., 1.], 0.),
+                        Ray::new([10., x, y], [0., 0., -1.], 0.),
+                    ]
+                    .into_iter()
+                })
+            })
+            .collect();
+        // These currently differ between STL and cuboid.
+        if false {
+            // Outward in at an angle.
+            let sqrt2 = 2f32.sqrt();
+            rays.extend(vec![
+                Ray::new([-1., 10., 10.], [sqrt2, sqrt2, 0.], 0.),
+                Ray::new([-1., 10., 10.], [sqrt2, -sqrt2, 0.], 0.),
+                Ray::new([-1., 10., 10.], [sqrt2, 0., sqrt2], 0.),
+                Ray::new([-1., 10., 10.], [sqrt2, 0., -sqrt2], 0.),
+            ]);
+        }
 
-        for (i, r) in rays.into_iter().enumerate() {
-            let c_hit = c.hit(r, 0., f32::MAX).expect("c_hit missed");
-            let s_hit = s.hit(r, 0., f32::MAX).expect("s_hit missed");
-            assert_eq!(
-                c_hit.t, s_hit.t,
-                "{i}: [t] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
+        // TODO(wathiede): proptest this, it's still not perfectly equal when rendering.
+
+        for r in rays.into_iter() {
+            let c_hit = c
+                .hit(r, 0., f32::MAX)
+                .expect(&format!("c_hit missed {r:#?}"));
+            let s_hit = s
+                .hit(r, 0., f32::MAX)
+                .expect(&format!("s_hit missed {r:#?}"));
+            assert!(
+                (c_hit.t - s_hit.t).abs() < EPSILON,
+                "{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:?}"
             );
         }
     }

rtiow/renderer/src/material.rs

@@ -251,6 +251,20 @@ where
     }
 }
 
+#[derive(Debug)]
+pub struct DebugMaterial {}
+
+impl Material for DebugMaterial {
+    fn scatter(&self, _r_in: &Ray, rec: &HitRecord) -> ScatterResponse {
+        let dir = Vec3::new(0., -1., -1.).unit_vector();
+        ScatterResponse {
+            scattered: Ray::new(rec.p, dir, 0.),
+            attenutation: [1., 1., 1.].into(),
+            reflected: true,
+        }
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;

rtiow/renderer/src/renderer.rs

@@ -501,7 +501,8 @@ fn render_worker(
 }
 
 pub fn render(scene: Scene, output_dir: &Path) -> std::result::Result<(), std::io::Error> {
-    let num_threads = scene.num_threads.unwrap_or_else(num_cpus::get);
+    // Default to half the cores to disable hyperthreading.
+    let num_threads = scene.num_threads.unwrap_or_else(|| num_cpus::get() / 2);
     let (pixel_req_tx, pixel_req_rx) = sync_channel(2 * num_threads);
     let (pixel_resp_tx, pixel_resp_rx) = sync_channel(2 * num_threads);
 

rtiow/vec3/src/vec3.rs

@@ -102,7 +102,7 @@ impl From<[f32; 3]> for Vec3 {
 
 impl fmt::Display for Vec3 {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        write!(f, "{} {} {}", self.x, self.y, self.z)
+        write!(f, "{:>6.2} {:>6.2} {:>6.2}", self.x, self.y, self.z)
     }
 }