rtiow: add support to send rendering output to https://github.com/Tom94/tev

Add debugging hit DebugHit.
Move some color helpers from mandelbrot to colors mod.
Implement Texture trait for [f32;3]

rtiow/Cargo.toml

@@ -1,13 +1,15 @@
 [workspace]
+resolver = "2"
 
 members = [
   "noise_explorer",
   "noise_explorer_warp",
   "renderer",
   "tracer",
+  "vec3",
 ]
-
-[profile]
 [profile.release]
 debug = true
 
+[profile.dev]
+opt-level = 3 

rtiow/noise_explorer/Cargo.toml

@@ -2,19 +2,19 @@
 name = "noise_explorer"
 version = "0.1.0"
 authors = ["Bill Thiede <git@xinu.tv>"]
-edition = "2018"
+edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-actix-web = "0.7.8"
-askama = "0.7.1"
-image = "0.22.3"
-log = "0.4.5"
-rand = "0.8.3"
+actix-web = "0.7.19"
+askama = "0.7.2"
+image = "0.22.5"
+log = "0.4.17"
+rand = "0.8.5"
 rand_xorshift = "0.3.0"
 renderer = { path = "../renderer" }
-serde = "1.0.79"
-serde_derive = "1.0.79"
-stderrlog = "0.4.1"
-structopt = "0.2.10"
+serde = "1.0.152"
+serde_derive = "1.0.152"
+stderrlog = "0.4.3"
+structopt = "0.2.18"

rtiow/noise_explorer_warp/Cargo.toml

@@ -2,18 +2,18 @@
 name = "noise_explorer_warp"
 version = "0.1.0"
 authors = ["Bill Thiede <git@xinu.tv>"]
-edition = "2018"
+edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-askama = "0.7.1"
-image = "0.22.3"
-log = "0.4.5"
-rand = "0.5.5"
+askama = "0.7.2"
+image = "0.22.5"
+log = "0.4.17"
+rand = "0.5.6"
 renderer = { path = "../renderer" }
-serde = "1.0.79"
-serde_derive = "1.0.79"
-stderrlog = "0.4.1"
-structopt = "0.2.10"
-warp = "0.1.20"
+serde = "1.0.152"
+serde_derive = "1.0.152"
+stderrlog = "0.4.3"
+structopt = "0.2.18"
+warp = "0.1.23"

rtiow/renderer/Cargo.toml

@@ -2,29 +2,41 @@
 name = "renderer"
 version = "0.1.0"
 authors = ["Bill Thiede <git@xinu.tv>"]
-edition = "2018"
+edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 [[bench]]
 harness = false
 name = "spheres"
+[[bench]]
+harness = false
+name = "aabb"
 
 [dependencies]
 chrono = "*"
 core_affinity = "0.5"
 cpuprofiler = { version = "0.0.3", optional = true }
-image = "0.22.3"
-lazy_static = "1.1.0"
-log = "0.4.5"
-num_cpus = "1.8.0"
-rand = "0.8.3"
-serde = "1.0.79"
-serde_derive = "1.0.79"
-serde_json = "1.0.41"
-structopt = "0.2.10"
+image = "0.22.5"
+lazy_static = "1.4.0"
+log = "0.4.17"
+num_cpus = "1.15.0"
+rand = "0.8.5"
+serde = "1.0.152"
+serde_derive = "1.0.152"
+serde_json = "1.0.93"
+structopt = "0.2.18"
+vec3 = {path = "../vec3"}
+stl = {path = "../../../stl"}
+strum = { version = "0.24.1", features = ["derive"] }
+strum_macros = "0.24.3"
+thiserror = "1.0.38"
+tev_client = "0.5.2"
+#stl = {git = "https://git-private.z.xinu.tv/wathiede/stl"}
 
 [dev-dependencies]
-criterion = "0.2"
+criterion = "0.4"
+pretty_assertions = "1.3.0"
+proptest = "1.1.0"
 
 [features]
 profile = ["cpuprofiler"]

rtiow/renderer/benches/aabb.rs

@@ -0,0 +1,47 @@
+use criterion::*;
+use renderer::{aabb::AABB, ray::Ray};
+
+fn bench(c: &mut Criterion) {
+    let bb = AABB::new([1., -1., -1.], [3., 1., 1.]);
+    let r_hit = Ray::new([0., 0., 0.], [1., 0., 0.], 0.);
+    let r_miss = Ray::new([0., 0., 0.], [-1., 0., 0.], 0.);
+    let t_min = 0.001;
+    let t_max = f32::MAX;
+
+    let mut group = c.benchmark_group("aabb");
+    group.throughput(Throughput::Elements(1));
+    group.bench_with_input(BenchmarkId::new("hit_naive", "r_hit"), &r_hit, |b, r| {
+        b.iter(|| bb.hit_naive(*r, t_min, t_max))
+    });
+    group.bench_with_input(BenchmarkId::new("hit2", "r_hit"), &r_hit, |b, r| {
+        b.iter(|| bb.hit2(*r, t_min, t_max))
+    });
+    //group.bench_with_input(BenchmarkId::new("hit_precompute", "r_hit"), &r_hit, |b, r| { b.iter(|| bb.hit_precompute(*r, t_min, t_max)) });
+    group.bench_with_input(BenchmarkId::new("hit_fast", "r_hit"), &r_hit, |b, r| {
+        b.iter(|| bb.hit_fast(*r, t_min, t_max))
+    });
+    #[cfg(target_arch = "x86_64")]
+    group.bench_with_input(BenchmarkId::new("hit_simd", "r_hit"), &r_hit, |b, r| {
+        b.iter(|| bb.hit_simd(*r, t_min, t_max))
+    });
+
+    group.bench_with_input(BenchmarkId::new("hit_naive", "r_miss"), &r_miss, |b, r| {
+        b.iter(|| bb.hit_naive(*r, t_min, t_max))
+    });
+    group.bench_with_input(BenchmarkId::new("hit2", "r_miss"), &r_miss, |b, r| {
+        b.iter(|| bb.hit2(*r, t_min, t_max))
+    });
+    //group.bench_with_input(BenchmarkId::new("hit_precompute", "r_miss"), &r_miss, |b, r| { b.iter(|| bb.hit_precompute(*r, t_min, t_max)) });
+    group.bench_with_input(BenchmarkId::new("hit_fast", "r_miss"), &r_miss, |b, r| {
+        b.iter(|| bb.hit_fast(*r, t_min, t_max))
+    });
+    #[cfg(target_arch = "x86_64")]
+    group.bench_with_input(BenchmarkId::new("hit_simd", "r_miss"), &r_miss, |b, r| {
+        b.iter(|| bb.hit_simd(*r, t_min, t_max))
+    });
+
+    group.finish();
+}
+
+criterion_group!(benches, bench);
+criterion_main!(benches);

rtiow/renderer/benches/spheres.rs

@@ -1,9 +1,6 @@
-#[macro_use]
-extern crate criterion;
+use criterion::*;
 
-use criterion::{Criterion, ParameterizedBenchmark};
-
-use rtiow::{
+use renderer::{
     hitable::Hit, material::Lambertian, ray::Ray, sphere::Sphere, texture::ConstantTexture,
     vec3::Vec3,
 };
@@ -20,14 +17,15 @@ fn criterion_benchmark(c: &mut Criterion) {
         // Miss
         Ray::new([0., 0., -2.], [0., 0., -1.], 0.),
     ];
-    c.bench(
-        "sphere",
-        ParameterizedBenchmark::new(
-            "Sphere",
-            move |b, r| b.iter(|| sphere.hit(*r, 0., 1.)),
-            rays,
-        ),
-    );
+    let mut group = c.benchmark_group("sphere");
+    group.throughput(Throughput::Elements(1));
+    group.bench_with_input(BenchmarkId::new("Sphere", "hit"), &rays[0], |b, r| {
+        b.iter(|| sphere.hit(*r, 0., 1.))
+    });
+    group.bench_with_input(BenchmarkId::new("Sphere", "miss"), &rays[1], |b, r| {
+        b.iter(|| sphere.hit(*r, 0., 1.))
+    });
+    group.finish()
 }
 
 criterion_group!(benches, criterion_benchmark);

rtiow/renderer/src/aabb.rs

@@ -2,7 +2,7 @@ use std::fmt;
 
 use crate::{ray::Ray, vec3::Vec3};
 
-#[derive(Debug, Clone, PartialEq)]
+#[derive(Default, Copy, Clone, PartialEq)]
 pub struct AABB {
     bounds: [Vec3; 2],
 }
@@ -29,15 +29,65 @@ 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: {} Area: {}",
+            self.bounds[0],
+            self.bounds[1],
+            self.volume(),
+            self.area()
+        )
+    }
+}
+
 impl AABB {
-    pub fn new(min: Vec3, max: Vec3) -> AABB {
+    // Create AABB with min = f32::MAX and max = -f32::MAX. It is expected the caller will use grow to create
+    // a vaild AABB.
+    pub fn infinite() -> AABB {
+        AABB {
+            bounds: [Vec3::from(f32::MAX), Vec3::from(f32::MIN)],
+        }
+    }
+    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);
         AABB { bounds: [min, max] }
     }
 
+    pub fn area(&self) -> f32 {
+        if self.max().x == f32::MIN || self.min().x == f32::MAX {
+            return 0.;
+        }
+        let e = self.max() - self.min();
+        let v = e.x * e.y + e.y * e.z + e.z * e.x;
+        if v.is_finite() {
+            v
+        } else {
+            0.
+        }
+    }
     pub fn volume(&self) -> f32 {
-        (self.min().x - self.max().x).abs()
+        if self.max().x == f32::MIN || self.min().x == f32::MAX {
+            return 0.;
+        }
+        let v = (self.min().x - self.max().x).abs()
             * (self.min().y - self.max().y).abs()
-            * (self.min().z - self.max().z).abs()
+            * (self.min().z - self.max().z).abs();
+        if v.is_finite() {
+            v
+        } else {
+            0.
+        }
+    }
+
+    pub fn grow(&mut self, v: Vec3) {
+        self.bounds[0] = vec3::min(self.bounds[0], v);
+        self.bounds[1] = vec3::max(self.bounds[1], v);
     }
 
     pub fn longest_axis(&self) -> usize {
@@ -61,11 +111,39 @@ impl AABB {
     pub fn min(&self) -> Vec3 {
         self.bounds[0]
     }
+
     pub fn max(&self) -> Vec3 {
         self.bounds[1]
     }
 
-    pub fn hit2(&self, r: Ray, t_min: f32, t_max: f32) -> bool {
+    // TODO(wathiede): implement branchless https://tavianator.com/cgit/dimension.git/tree/libdimension/bvh/bvh.c#n194
+
+    pub fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> bool {
+        self.hit_simd(r, t_min, t_max).is_some()
+    }
+
+    pub fn hit_distance(&self, r: Ray, t_min: f32, t_max: f32) -> Option<f32> {
+        self.hit_simd(r, t_min, t_max)
+    }
+
+    pub fn hit_naive(&self, r: Ray, t_min: f32, t_max: f32) -> Option<f32> {
+        let mut t_min = t_min;
+        let mut t_max = t_max;
+        for axis in 0..3 {
+            let t0 = ((self.min()[axis] - r.origin[axis]) * r.inv_direction[axis])
+                .min((self.max()[axis] - r.origin[axis]) * r.inv_direction[axis]);
+            let t1 = ((self.min()[axis] - r.origin[axis]) * r.inv_direction[axis])
+                .max((self.max()[axis] - r.origin[axis]) * r.inv_direction[axis]);
+            t_min = t0.max(t_min);
+            t_max = t1.min(t_max);
+            if t_max <= t_min {
+                return None;
+            }
+        }
+        Some(t_min)
+    }
+
+    pub fn hit2(&self, r: Ray, t_min: f32, t_max: f32) -> Option<f32> {
         let mut t_min = t_min;
         let mut t_max = t_max;
         for axis in 0..3 {
@@ -81,13 +159,13 @@ impl AABB {
             t_min = max(t0, t_min);
             t_max = min(t1, t_max);
             if t_max <= t_min {
-                return false;
+                return None;
             }
         }
-        true
+        Some(t_min)
     }
 
-    pub fn hit_precompute(&self, r: Ray, t0: f32, t1: f32) -> bool {
+    pub fn hit_precompute(&self, r: Ray, t0: f32, t1: f32) -> Option<f32> {
         // TODO(wathiede): this has bugs.
         let mut t_min = (self.bounds[r.sign[0]].x - r.origin.x) * r.inv_direction.x;
         let mut t_max = (self.bounds[1 - r.sign[0]].x - r.origin.x) * r.inv_direction.x;
@@ -95,7 +173,7 @@ impl AABB {
         let t_y_max = (self.bounds[1 - r.sign[0]].y - r.origin.y) * r.inv_direction.y;
 
         if t_min > t_y_max || t_y_min > t_max {
-            return false;
+            return None;
         }
 
         if t_y_min > t_min {
@@ -108,7 +186,7 @@ impl AABB {
         let t_z_min = (self.bounds[r.sign[2]].z - r.origin.z) * r.inv_direction.z;
         let t_z_max = (self.bounds[1 - r.sign[2]].z - r.origin.z) * r.inv_direction.z;
         if t_min > t_z_max || t_z_min > t_max {
-            return false;
+            return None;
         }
         if t_z_min > t_min {
             t_min = t_z_min;
@@ -116,27 +194,42 @@ impl AABB {
         if t_z_max < t_max {
             t_max = t_z_max;
         }
-        t_min < t1 && t_max > t0
+        if t_min < t1 && t_max > t0 {
+            Some(t_min)
+        } else {
+            None
+        }
     }
 
-    pub fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> bool {
-        let mut t_min = t_min;
-        let mut t_max = t_max;
-        for axis in 0..3 {
-            let t0 = ((self.min()[axis] - r.origin[axis]) * r.inv_direction[axis])
-                .min((self.max()[axis] - r.origin[axis]) * r.inv_direction[axis]);
-            let t1 = ((self.min()[axis] - r.origin[axis]) * r.inv_direction[axis])
-                .max((self.max()[axis] - r.origin[axis]) * r.inv_direction[axis]);
-            t_min = t0.max(t_min);
-            t_max = t1.min(t_max);
-            if t_max <= t_min {
-                return false;
+    pub fn hit_simd(&self, r: Ray, t_min: f32, t_max: f32) -> Option<f32> {
+        #[cfg(target_arch = "x86_64")]
+        unsafe {
+            use std::arch::x86_64::*;
+            let o4 = _mm_set_ps(0., r.origin.z, r.origin.y, r.origin.x);
+            let d4 = _mm_set_ps(0., r.inv_direction.z, r.inv_direction.y, r.inv_direction.x);
+            let bmin4 = _mm_set_ps(0., self.min().z, self.min().y, self.min().x);
+            let bmax4 = _mm_set_ps(0., self.max().z, self.max().y, self.max().x);
+            let mask4 = _mm_cmpeq_ps(_mm_setzero_ps(), _mm_set_ps(1., 0., 0., 0.));
+            let t1 = _mm_mul_ps(_mm_sub_ps(_mm_and_ps(bmin4, mask4), o4), d4);
+            let t2 = _mm_mul_ps(_mm_sub_ps(_mm_and_ps(bmax4, mask4), o4), d4);
+            let vmax4 = _mm_max_ps(t1, t2);
+            let vmin4 = _mm_min_ps(t1, t2);
+            let vmax4: (f32, f32, f32, f32) = std::mem::transmute(vmax4);
+            let vmin4: (f32, f32, f32, f32) = std::mem::transmute(vmin4);
+            let tmax = min(min(vmax4.0, min(vmax4.1, vmax4.2)), t_max);
+            let tmin = max(max(vmin4.0, max(vmin4.1, vmin4.2)), t_min);
+            if tmin <= tmax {
+                Some(tmin)
+            } else {
+                None
             }
         }
-        true
+        #[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 {
+    pub fn hit_fast(&self, r: Ray, _t_min: f32, _t_max: f32) -> Option<f32> {
         let b = self;
         let mut t1 = (b.min()[0] - r.origin[0]) * 1. / r.direction[0];
         let mut t2 = (b.max()[0] - r.origin[0]) * 1. / r.direction[0];
@@ -152,7 +245,11 @@ impl AABB {
             tmax = tmax.min(t1.max(t2));
         }
 
-        tmax > tmin.max(0.0)
+        if tmax > tmin.max(0.0) {
+            Some(tmin)
+        } else {
+            None
+        }
     }
 }
 
@@ -169,3 +266,112 @@ pub fn surrounding_box(box0: &AABB, box1: &AABB) -> AABB {
     );
     AABB::new(min, max)
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn infinite() {
+        let bb = AABB::infinite();
+        assert_eq!(bb.area(), 0.);
+        assert_eq!(bb.volume(), 0.);
+    }
+
+    macro_rules! hit_test {
+        ($($name:ident,)*) => {
+            $(
+                mod $name {
+                    use super::*;
+                    const T_MIN: f32 = 0.001;
+                    const T_MAX: f32 = f32::MAX;
+
+                    fn test_bb() -> AABB {
+                        AABB::new([-1., -1., -1.], [1., 1., 1.])
+                    }
+
+                    #[test]
+                    fn hit_front() {
+                        let bb = test_bb();
+                        let r = Ray::new([-2., 0., 0.], [1., 0., 0.], 0.5);
+                        assert!(bb.$name(r, T_MIN, T_MAX).is_some());
+                    }
+                    #[test]
+                    fn hit_back() {
+                        let bb = test_bb();
+                        let r = Ray::new([2., 0., 0.], [-1., 0., 0.], 0.5);
+                        assert!(bb.$name(r, T_MIN, T_MAX).is_some());
+                    }
+                    #[test]
+                    fn hit_top() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., 2., 0.], [0., -1., 0.], 0.5);
+                        assert!(bb.$name(r, T_MIN, T_MAX).is_some());
+                    }
+                    #[test]
+                    fn hit_bottom() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., -2., 0.], [0., 1., 0.], 0.5);
+                        assert!(bb.$name(r, T_MIN, T_MAX).is_some());
+                    }
+                    #[test]
+                    fn hit_left() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., 0., -2.], [0., 0., 1.], 0.5);
+                        assert!(bb.$name(r, T_MIN, T_MAX).is_some());
+                    }
+                    #[test]
+                    fn hit_right() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., 0., 2.], [0., 0., -1.], 0.5);
+                        assert!(bb.$name(r, T_MIN, T_MAX).is_some());
+                    }
+
+                    #[test]
+                    fn miss_front() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., 1.1, -1.1], [0., -1., 0.], 0.5);
+                        assert_eq!(bb.$name(r, T_MIN, T_MAX), None);
+                    }
+                    #[test]
+                    fn miss_back() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., 1.1, 1.1], [0., -1., 0.], 0.5);
+                        assert_eq!(bb.$name(r, T_MIN, T_MAX), None);
+                    }
+                    #[test]
+                    fn miss_top() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., 1.1, -1.1], [0., 0., 1.], 0.5);
+                        assert_eq!(bb.$name(r, T_MIN, T_MAX), None);
+                    }
+                    #[test]
+                    fn miss_bottom() {
+                        let bb = test_bb();
+                        let r = Ray::new([0., -1.1, -1.1], [0., 0., 1.], 0.5);
+                        assert_eq!(bb.$name(r, T_MIN, T_MAX), None);
+                    }
+                    #[test]
+                    fn miss_left() {
+                        let bb = test_bb();
+                        let r = Ray::new([-1.1, 0., -1.1], [0., 0., 1.], 0.5);
+                        assert_eq!(bb.$name(r, T_MIN, T_MAX), None);
+                    }
+                    #[test]
+                    fn miss_right() {
+                        let bb = test_bb();
+                        let r = Ray::new([1.1, 0., -1.1], [0., 0., 1.], 0.5);
+                        assert_eq!(bb.$name(r, T_MIN, T_MAX), None);
+                    }
+                }
+            )*
+        }
+    }
+
+    hit_test! {
+        hit_naive,
+        hit_simd,
+        hit_fast,
+        hit2,
+    }
+}

rtiow/renderer/src/bvh.rs

@@ -9,6 +9,7 @@ use crate::{
     ray::Ray,
 };
 
+#[derive(Debug)]
 enum BVHNode {
     Leaf(Box<dyn Hit>),
     Branch {
@@ -178,6 +179,7 @@ impl Hit for BVHNode {
     }
 }
 
+#[derive(Debug)]
 pub struct BVH {
     root: BVHNode,
 }

rtiow/renderer/src/bvh_triangles.rs

@@ -0,0 +1,711 @@
+/// 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 log::info;
+use stl::STL;
+
+use crate::{
+    aabb::AABB,
+    hitable::{Hit, HitRecord},
+    material::Material,
+    ray::Ray,
+    vec3::{cross, dot, Vec3},
+};
+
+#[derive(Debug, PartialEq)]
+struct BVHNode {
+    aabb: AABB,
+    // When tri_count==0, left_first holds the left child's index in bvh_nodes. When >0 left_first
+    // holds the index for the first triangle in triangles.
+    left_first: u32,
+    tri_count: u32,
+}
+
+impl BVHNode {
+    fn is_leaf(&self) -> bool {
+        self.tri_count > 0
+    }
+
+    fn cost(&self) -> f32 {
+        let area = self.aabb.area();
+        self.tri_count as f32 * area
+    }
+}
+
+#[derive(Clone, 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
+        )
+    }
+}
+
+pub struct BVHTriangles<M>
+where
+    M: Material,
+{
+    pub triangles: Vec<Triangle>,
+    triangle_index: Vec<usize>,
+    material: M,
+    bvh_nodes: Vec<BVHNode>,
+}
+
+#[derive(Debug, Default)]
+struct Bin {
+    bounds: AABB,
+    tri_count: usize,
+}
+
+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.left_first..(n.left_first + n.tri_count) {
+                    if f.alternate() {
+                        write!(f, "\t")?;
+                    }
+                    write!(f, "{:?} ", self.triangles[t_idx as usize])?;
+                    if f.alternate() {
+                        writeln!(f)?;
+                    }
+                }
+            }
+        }
+        Ok(())
+    }
+}
+
+#[derive(Debug)]
+struct SplitCost {
+    pos: f32,
+    axis: usize,
+    cost: f32,
+}
+
+const ROOT_NODE_IDX: usize = 0;
+impl<M> BVHTriangles<M>
+where
+    M: Material,
+{
+    pub fn new(stl: &STL, material: M, scale_factor: f32) -> BVHTriangles<M> {
+        let now = std::time::Instant::now();
+
+        assert_eq!(std::mem::size_of::<BVHNode>(), 32);
+        let div3 = 1. / 3.;
+        let triangles: Vec<_> = stl
+            .triangles
+            .iter()
+            .map(|t| {
+                let v0 = t.verts[0] * scale_factor;
+                let v1 = t.verts[1] * scale_factor;
+                let v2 = t.verts[2] * scale_factor;
+                let centroid = (v0 + v1 + v2) * div3;
+                Triangle {
+                    centroid,
+                    verts: [v0, v1, v2],
+                }
+            })
+            .collect();
+        let triangle_index = (0..triangles.len()).collect();
+
+        let n = 2 * triangles.len() - 2;
+        let bvh_nodes = Vec::with_capacity(n);
+        let mut bvh = BVHTriangles {
+            triangles,
+            triangle_index,
+            bvh_nodes,
+            material,
+        };
+        bvh.build_bvh();
+        info!(
+            "BVHTriangles build time {:0.3}s",
+            now.elapsed().as_secs_f32()
+        );
+        struct Stats {
+            nodes: usize,
+            leafs: usize,
+            min_tris: u32,
+            max_tris: u32,
+        }
+        impl Default for Stats {
+            fn default() -> Self {
+                Stats {
+                    nodes: 0,
+                    leafs: 0,
+                    min_tris: u32::MAX,
+                    max_tris: u32::MIN,
+                }
+            }
+        }
+        let stats = bvh.bvh_nodes.iter().fold(Stats::default(), |mut stats, n| {
+            stats.nodes += 1;
+            if n.is_leaf() {
+                stats.leafs += 1;
+                stats.min_tris = n.tri_count.min(stats.min_tris);
+                stats.max_tris = n.tri_count.max(stats.max_tris);
+            }
+            stats
+        });
+        info!("BVHTriangles build stats:");
+        info!("  Nodes:   {}", stats.nodes);
+        info!("  Leaves:  {}", stats.leafs);
+        info!("  Tris:    {}", bvh.triangles.len());
+        info!("  Min Tri: {}", stats.min_tris);
+        info!("  Max Tri: {}", stats.max_tris);
+        info!("  Avg Tri: {}", bvh.triangles.len() / stats.leafs);
+        info!("  Predict: {}", bvh.bvh_nodes.capacity());
+        info!("  Actual:  {}", bvh.bvh_nodes.len());
+        info!(
+            "  Savings:  {} bytes",
+            (bvh.bvh_nodes.capacity() - bvh.bvh_nodes.len()) * std::mem::size_of::<BVHNode>()
+        );
+        bvh.bvh_nodes.shrink_to_fit();
+        //dbg!(&bvh);
+        bvh
+    }
+
+    fn build_bvh(&mut self) {
+        // assign all triangles to root node
+        let root = BVHNode {
+            aabb: AABB::default(),
+            left_first: 0,
+            tri_count: self.triangles.len() as u32,
+        };
+        self.bvh_nodes.push(root);
+        self.update_node_bounds(ROOT_NODE_IDX);
+        // subdivide recursively
+        self.subdivide(ROOT_NODE_IDX);
+    }
+    fn update_node_bounds(&mut self, node_idx: usize) {
+        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.left_first..(node.left_first + node.tri_count) {
+            let leaf_tri_idx = self.triangle_index[i as usize];
+            let leaf_tri = &self.triangles[leaf_tri_idx];
+            aabb_min = vec3::min(aabb_min, leaf_tri.verts[0]);
+            aabb_min = vec3::min(aabb_min, leaf_tri.verts[1]);
+            aabb_min = vec3::min(aabb_min, leaf_tri.verts[2]);
+            aabb_max = vec3::max(aabb_max, leaf_tri.verts[0]);
+            aabb_max = vec3::max(aabb_max, leaf_tri.verts[1]);
+            aabb_max = vec3::max(aabb_max, leaf_tri.verts[2]);
+        }
+        node.aabb = AABB::new(aabb_min, aabb_max);
+    }
+
+    fn find_best_split_plane(&self, node: &BVHNode) -> SplitCost {
+        let mut best = SplitCost {
+            pos: 0.,
+            cost: f32::MAX,
+            axis: usize::MAX,
+        };
+
+        for axis in 0..3 {
+            let mut bounds_min = f32::MAX;
+            let mut bounds_max = f32::MIN;
+
+            for i in 0..node.tri_count {
+                let triangle = &self.triangles[self.triangle_index[(node.left_first + i) as usize]];
+                bounds_min = bounds_min.min(triangle.centroid[axis]);
+                bounds_max = bounds_max.max(triangle.centroid[axis]);
+            }
+
+            if bounds_min == bounds_max {
+                continue;
+            }
+
+            const NUM_BINS: usize = 8;
+            let mut bins: Vec<_> = (0..NUM_BINS)
+                .map(|_| Bin {
+                    bounds: AABB::infinite(),
+                    tri_count: 0,
+                })
+                .collect();
+            let scale = bins.len() as f32 / (bounds_max - bounds_min);
+            // populate the bins
+            for i in 0..node.tri_count {
+                let triangle = &self.triangles[self.triangle_index[(node.left_first + i) as usize]];
+                let bin_idx = (triangle.centroid[axis] - bounds_min) * scale;
+                let bin_idx = (bins.len() - 1).min(bin_idx as usize);
+
+                bins[bin_idx].tri_count += 1;
+                bins[bin_idx].bounds.grow(triangle.verts[0]);
+                bins[bin_idx].bounds.grow(triangle.verts[1]);
+                bins[bin_idx].bounds.grow(triangle.verts[2]);
+            }
+
+            // gather data for the 7 planes between the 8 bins
+            let mut left_area: Vec<_> = (0..bins.len() - 1).map(|_| 0.).collect();
+            let mut right_area: Vec<_> = (0..bins.len() - 1).map(|_| 0.).collect();
+            let mut left_count: Vec<_> = (0..bins.len() - 1).map(|_| 0).collect();
+            let mut right_count: Vec<_> = (0..bins.len() - 1).map(|_| 0).collect();
+            let mut left_box = AABB::infinite();
+            let mut right_box = AABB::infinite();
+            let mut left_sum = 0;
+            let mut right_sum = 0;
+            for i in 0..(bins.len() - 1) {
+                left_sum += bins[i].tri_count;
+                left_count[i] = left_sum;
+                left_box.grow(bins[i].bounds.min());
+                left_box.grow(bins[i].bounds.max());
+                left_area[i] = left_box.area();
+
+                right_sum += bins[bins.len() - 1 - i].tri_count;
+                right_count[bins.len() - 2 - i] = right_sum;
+                right_box.grow(bins[bins.len() - 1 - i].bounds.min());
+                right_box.grow(bins[bins.len() - 1 - i].bounds.max());
+                right_area[bins.len() - 2 - i] = right_box.area();
+            }
+
+            let scale = (bounds_max - bounds_min) / bins.len() as f32;
+            // calculate SAH cost for the 7 planes
+            for i in 0..(bins.len() - 1) {
+                let plane_cost =
+                    left_count[i] as f32 * left_area[i] + right_count[i] as f32 * right_area[i];
+                if plane_cost < best.cost {
+                    best.axis = axis;
+                    best.pos = bounds_min + scale * (i as f32 + 1.);
+                    best.cost = plane_cost;
+                }
+            }
+        }
+        best
+    }
+
+    fn subdivide(&mut self, idx: usize) {
+        let (left_first, tri_count, left_count, i) = {
+            let node = &self.bvh_nodes[idx];
+            let split = self.find_best_split_plane(&node);
+            let no_split_cost = node.cost();
+            // Stop subdividing if it isn't getting any better.
+            if split.cost >= no_split_cost {
+                return;
+            }
+
+            // Split the group in two halves.
+            let mut i = node.left_first as isize;
+            let mut j = i + node.tri_count as isize - 1;
+            while i <= j {
+                if self.triangles[self.triangle_index[i as usize]].centroid[split.axis] < split.pos
+                {
+                    i += 1;
+                } else {
+                    self.triangles.swap(
+                        self.triangle_index[i as usize],
+                        self.triangle_index[j as usize],
+                    );
+                    j -= 1;
+                }
+            }
+
+            // Create child nodes for each half.
+            let left_count = i as u32 - node.left_first;
+            if left_count == 0 || left_count == node.tri_count {
+                return;
+            }
+            (node.left_first, node.tri_count, left_count, i)
+        };
+
+        // create child nodes
+        let left_child_idx = self.bvh_nodes.len() as u32;
+        let right_child_idx = left_child_idx + 1 as u32;
+        let left = BVHNode {
+            aabb: AABB::default(),
+            left_first,
+            tri_count: left_count as u32,
+        };
+        let right = BVHNode {
+            aabb: AABB::default(),
+            left_first: i as u32,
+            tri_count: (tri_count - left_count) as u32,
+        };
+        self.bvh_nodes.push(left);
+        self.bvh_nodes.push(right);
+        let node = &mut self.bvh_nodes[idx];
+        node.left_first = left_child_idx;
+        node.tri_count = 0;
+
+        // Recurse
+        self.update_node_bounds(left_child_idx as usize);
+        self.subdivide(left_child_idx as usize);
+        self.update_node_bounds(right_child_idx as usize);
+        self.subdivide(right_child_idx as usize);
+    }
+
+    fn intersect_bvh(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        let mut node = &self.bvh_nodes[ROOT_NODE_IDX];
+        let mut stack = Vec::with_capacity(2);
+        let mut nearest = None;
+        loop {
+            if node.is_leaf() {
+                let canditate = (node.left_first..(node.left_first + node.tri_count))
+                    // Map from idx to Triangle
+                    .map(|idx| &self.triangles[self.triangle_index[idx as usize]])
+                    // Try to hit all triangles for this node, filtering out misses.
+                    .filter_map(|tri| intersect_tri(r, &tri))
+                    // Find the nearest hit (if any).
+                    .min_by(|a, b| a.t.partial_cmp(&b.t).unwrap());
+
+                // merge candidate with nearest.
+                nearest = match (&canditate, &nearest) {
+                    (Some(_), None) => canditate,
+                    (None, Some(_)) => nearest,
+                    (Some(c), Some(n)) => {
+                        //info!("merging {c:#?} and {n:#?}");
+                        if c.t < n.t {
+                            canditate
+                        } else {
+                            nearest
+                        }
+                    }
+                    (None, None) => None,
+                };
+                if stack.is_empty() {
+                    break;
+                }
+                node = stack.pop().unwrap();
+                continue;
+            }
+
+            let child1 = &self.bvh_nodes[node.left_first as usize];
+            let child2 = &self.bvh_nodes[node.left_first as usize + 1];
+            let dist1 = child1.aabb.hit_distance(r, t_min, t_max);
+            let dist2 = child2.aabb.hit_distance(r, t_min, t_max);
+            // Swap c1/c2 & d1/d2 based on d1/d2.
+            let (child1, child2, dist1, dist2) = match (dist1, dist2) {
+                (Some(d1), Some(d2)) if d1 > d2 => (child2, child1, dist2, dist1),
+                (None, Some(_)) => (child2, child1, dist2, dist1),
+                _ => (child1, child2, dist1, dist2),
+            };
+
+            // dist1/child1 should now be the nearest hit.
+
+            // If we missed dist1/child1, then we implicitly missed dist2/child2, so pop a child
+            // from the stack or exit the function.
+            if dist1.is_none() {
+                if stack.is_empty() {
+                    break;
+                }
+                node = stack.pop().unwrap();
+            } else {
+                // We hit child1, so process it next.
+                node = child1;
+                // If we also hit child2 save it on the stack so we can process it later.
+                if dist2.is_some() {
+                    stack.push(child2);
+                }
+            }
+        }
+        nearest
+            .and_then(|rtr: RayTriangleResult| Some(rtr.hit_record_with_material(&self.material)))
+    }
+}
+
+/// Based on https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection.html
+fn intersect_tri(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.abs() < 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.point_at_parameter(t),
+            tri: tri.clone(),
+        });
+    }
+    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, t_min, t_max)
+    }
+
+    fn bounding_box(&self, _t_min: f32, _t_max: f32) -> Option<AABB> {
+        Some(self.bvh_nodes[ROOT_NODE_IDX].aabb)
+    }
+}
+
+#[derive(Debug)]
+struct RayTriangleResult {
+    t: f32,
+    p: Vec3,
+    tri: Triangle,
+}
+
+impl RayTriangleResult {
+    fn hit_record_with_material<'m>(self, material: &'m dyn Material) -> HitRecord<'m> {
+        // We don't support UV (yet?).
+        let uv = (0.5, 0.5);
+        let v0 = self.tri.verts[0];
+        let v1 = self.tri.verts[1];
+        let v2 = self.tri.verts[2];
+
+        let v0v1 = v1 - v0;
+        let v0v2 = v2 - v0;
+        let normal = cross(v0v1, v0v2).unit_vector();
+        //println!("hit triangle {tri:?}");
+        HitRecord {
+            t: self.t,
+            uv,
+            p: self.p,
+            normal,
+            material,
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::{
+        bvh_triangles::BVHTriangles,
+        cuboid::Cuboid,
+        hitable::Hit,
+        material::{Dielectric, Lambertian},
+        ray::Ray,
+        texture::ConstantTexture,
+    };
+    use pretty_assertions::assert_eq;
+    use proptest::prelude::*;
+    use std::{
+        io::{BufReader, Cursor},
+        sync::Arc,
+    };
+    use stl::STL;
+
+    #[test]
+    fn compare_cuboid() {
+        let c = Cuboid::new(
+            [0., 0., 0.].into(),
+            [20., 20., 20.].into(),
+            Arc::new(Dielectric::new(1.5)),
+        );
+        let stl_cube = STL::parse(
+            BufReader::new(Cursor::new(include_bytes!("../stls/cube.stl"))),
+            false,
+        )
+        .expect("failed to parse cube");
+
+        let s = BVHTriangles::new(
+            &stl_cube,
+            Dielectric::new(1.5),
+            //Metal::new(Vec3::new(0.8, 0.8, 0.8), 0.2),
+            //Lambertian::new(ConstantTexture::new(Vec3::new(1.0, 0.2, 0.2))),
+        );
+        //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();
+        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 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,
+                "{r:?}: [p] c_hit: {c_hit:#?}, s_hit: {s_hit:#?}"
+            );
+            assert_eq!(
+                c_hit.normal, s_hit.normal,
+                "{r:?}: [normal] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
+            );
+        }
+    }
+
+    proptest! {
+        // TODO(wathiede): make this work.
+        //#[test]
+        fn compare_cuboid_proptest(
+            ox in -20.0f32..40.0,
+            oy in -20.0f32..40.0,
+            oz in -20.0f32..40.0,
+            dx in -1.0f32..1.0,
+            dy in -1.0f32..1.0,
+            dz in -1.0f32..1.0,
+        ) {
+            let r = Ray::new([ox,oy,oz].into(), Vec3::new(dx, dy, dz).unit_vector(), 0.5);
+            let c = Cuboid::new(
+                [0., 0., 0.].into(),
+                [20., 20., 20.].into(),
+                Arc::new(Dielectric::new(1.5)),
+            );
+            let stl_cube = STL::parse(
+                BufReader::new(Cursor::new(include_bytes!("../stls/cube.stl"))),
+                false,
+            )
+                .expect("failed to parse cube");
+
+            let s = BVHTriangles::new(
+                &stl_cube,
+                Dielectric::new(1.5),
+                //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 c_hit = c .hit(r, 0., f32::MAX);
+            let s_hit = s .hit(r, 0., f32::MAX);
+
+            match (c_hit, s_hit) {
+                (Some(_), None)=>assert!(false, "hit cuboid but not STL"),
+                (None, Some(_))=>assert!(false, "hit STL but not cuboid"),
+                (Some(c_hit), Some(s_hit))=> {
+                    assert!(
+                        (c_hit.t - s_hit.t).abs() < 0.00001,
+                        "{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,
+                        "{r:?}: [p] c_hit: {c_hit:#?}, s_hit: {s_hit:#?}"
+                    );
+                    assert_eq!(
+                        c_hit.normal, s_hit.normal,
+                        "{r:?}: [normal] c_hit: {c_hit:?}, s_hit: {s_hit:?}"
+                    );
+                },
+                // It's okay if they both miss.
+                (None,None)=>(),
+            }
+        }
+    }
+}

rtiow/renderer/src/camera.rs

@@ -18,6 +18,7 @@ fn random_in_unit_disk() -> Vec3 {
     }
 }
 
+#[derive(Debug)]
 pub struct Camera {
     origin: Vec3,
     lower_left_corner: Vec3,

rtiow/renderer/src/colors.rs

@@ -0,0 +1,46 @@
+use rand::{self, Rng};
+
+use crate::vec3::Vec3;
+
+// HSV values in [0..1]
+// returns [r, g, b] values from 0 to 255
+//From https://martin.ankerl.com/2009/12/09/how-to-create-random-colors-programmatically/
+pub fn hsv_to_rgb(h: f32, s: f32, v: f32) -> Vec3 {
+    let h_i = (h * 6.) as i32;
+    let f = h * 6. - h_i as f32;
+    let p = v * (1. - s);
+    let q = v * (1. - f * s);
+    let t = v * (1. - (1. - f) * s);
+    match h_i {
+        0 => Vec3::new(v, t, p),
+        1 => Vec3::new(q, v, p),
+        2 => Vec3::new(p, v, t),
+        3 => Vec3::new(p, q, v),
+        4 => Vec3::new(t, p, v),
+        5 => Vec3::new(v, p, q),
+        _ => panic!("Unknown H value {}", h_i),
+    }
+}
+
+pub fn generate_rainbow(num: usize) -> Vec<Vec3> {
+    (0..num)
+        .map(|n| {
+            let h = n as f32 / num as f32;
+            hsv_to_rgb(h, 0.99, 0.99)
+        })
+        .collect()
+}
+pub fn generate_palette(num: usize) -> Vec<Vec3> {
+    let mut rng = rand::thread_rng();
+    let mut random = || rng.gen();
+    // use golden ratio
+    let golden_ratio_conjugate = 0.618_034;
+    let mut h = random();
+    (0..num)
+        .map(|_| {
+            h += golden_ratio_conjugate;
+            h %= 1.0;
+            hsv_to_rgb(h, 0.99, 0.99)
+        })
+        .collect()
+}

rtiow/renderer/src/constant_medium.rs

@@ -11,6 +11,7 @@ use crate::{
     vec3::Vec3,
 };
 
+#[derive(Debug)]
 pub struct ConstantMedium<H, T>
 where
     H: Hit,

rtiow/renderer/src/cuboid.rs

@@ -11,6 +11,7 @@ use crate::{
     vec3::Vec3,
 };
 
+#[derive(Debug)]
 pub struct Cuboid {
     p_min: Vec3,
     p_max: Vec3,
@@ -21,6 +22,9 @@ impl Cuboid {
     // This clippy doesn't work right with Arc.
     #[allow(clippy::needless_pass_by_value)]
     pub fn new(p_min: Vec3, p_max: Vec3, material: Arc<dyn Material>) -> Cuboid {
+        assert!(p_min.x <= p_max.x);
+        assert!(p_min.y <= p_max.y);
+        assert!(p_min.z <= p_max.z);
         Cuboid {
             p_min,
             p_max,

rtiow/renderer/src/debug_hit.rs

@@ -0,0 +1,44 @@
+use crate::{
+    aabb::AABB,
+    hitable::{Hit, HitRecord},
+    material::Lambertian,
+    ray::Ray,
+};
+
+#[derive(Debug)]
+pub struct DebugHit<H>
+where
+    H: Hit,
+{
+    hitable: H,
+    material: Lambertian<[f32; 3]>,
+}
+
+impl<H> DebugHit<H>
+where
+    H: Hit,
+{
+    pub fn new(hitable: H) -> DebugHit<H>
+where {
+        DebugHit {
+            hitable,
+            material: Lambertian::new([0.2, 0.2, 0.2]),
+        }
+    }
+}
+
+impl<H> Hit for DebugHit<H>
+where
+    H: Hit,
+{
+    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        if let Some(hit) = self.hitable.hit(r, t_min, t_max) {
+            return Some(HitRecord { t: hit.t, ..hit });
+        }
+        None
+    }
+
+    fn bounding_box(&self, t_min: f32, t_max: f32) -> Option<AABB> {
+        self.hitable.bounding_box(t_min, t_max)
+    }
+}

rtiow/renderer/src/flip_normals.rs

@@ -4,6 +4,7 @@ use crate::{
     ray::Ray,
 };
 
+#[derive(Debug)]
 pub struct FlipNormals<H>
 where
     H: Hit,

rtiow/renderer/src/glowybox.rs

@@ -0,0 +1,139 @@
+use std::sync::Arc;
+
+use crate::{
+    aabb::AABB,
+    cuboid::Cuboid,
+    hitable::{Hit, HitRecord},
+    material::Material,
+    ray::Ray,
+    vec3::Vec3,
+};
+
+#[derive(Debug)]
+pub struct Glowybox {
+    p_min: Vec3,
+    p_max: Vec3,
+    main: Cuboid,
+    edges: [Cuboid; 12],
+}
+
+impl Glowybox {
+    // This clippy doesn't work right with Arc.
+    #[allow(clippy::needless_pass_by_value)]
+    pub fn new(
+        p_min: Vec3,
+        p_max: Vec3,
+        edge_thickness: f32,
+        main_material: Arc<dyn Material>,
+        edge_material: Arc<dyn Material>,
+    ) -> Glowybox {
+        assert!(p_min.x < p_max.x);
+        assert!(p_min.y < p_max.y);
+        assert!(p_min.z < p_max.z);
+        let main = Cuboid::new(p_min, p_max, main_material);
+        // Top edges
+        let ht = edge_thickness / 2.;
+        let edges = [
+            // Top edges
+            Cuboid::new(
+                [p_min.x - ht, p_max.y - ht, p_min.z - ht].into(),
+                [p_min.x + ht, p_max.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_min.x - ht, p_max.y - ht, p_min.z - ht].into(),
+                [p_max.x + ht, p_max.y + ht, p_min.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_max.x - ht, p_max.y - ht, p_min.z - ht].into(),
+                [p_max.x + ht, p_max.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_min.x - ht, p_max.y - ht, p_max.z - ht].into(),
+                [p_max.x + ht, p_max.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            // Bottom edges
+            Cuboid::new(
+                [p_min.x - ht, p_min.y - ht, p_min.z - ht].into(),
+                [p_min.x + ht, p_min.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_min.x - ht, p_min.y - ht, p_min.z - ht].into(),
+                [p_max.x + ht, p_min.y + ht, p_min.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_max.x - ht, p_min.y - ht, p_min.z - ht].into(),
+                [p_max.x + ht, p_min.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_min.x - ht, p_min.y - ht, p_max.z - ht].into(),
+                [p_max.x + ht, p_min.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            // Middle edges
+            Cuboid::new(
+                [p_min.x - ht, p_min.y - ht, p_min.z - ht].into(),
+                [p_min.x + ht, p_max.y + ht, p_min.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_min.x - ht, p_min.y - ht, p_max.z - ht].into(),
+                [p_min.x + ht, p_max.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_max.x - ht, p_min.y - ht, p_min.z - ht].into(),
+                [p_max.x + ht, p_max.y + ht, p_min.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+            Cuboid::new(
+                [p_max.x - ht, p_min.y - ht, p_max.z - ht].into(),
+                [p_max.x + ht, p_max.y + ht, p_max.z + ht].into(),
+                Arc::clone(&edge_material),
+            ),
+        ];
+        let p_min = [p_min.x - ht, p_min.y - ht, p_min.z - ht].into();
+        let p_max = [p_max.x + ht, p_max.y + ht, p_max.z + ht].into();
+        Glowybox {
+            p_min,
+            p_max,
+            main,
+            edges,
+        }
+    }
+}
+
+impl Hit for Glowybox {
+    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        let mut edge_hit = None;
+        for edge in &self.edges {
+            if let Some(hit) = edge.hit(r, t_min, t_max) {
+                edge_hit = Some(hit);
+                break;
+            }
+        }
+        let main_hit = self.main.hit(r, t_min, t_max);
+        match (edge_hit, main_hit) {
+            (Some(ehit), Some(mhit)) => {
+                if mhit.t < ehit.t {
+                    Some(mhit)
+                } else {
+                    Some(ehit)
+                }
+            }
+            (Some(ehit), None) => Some(ehit),
+            (None, Some(mhit)) => Some(mhit),
+            _ => None,
+        }
+    }
+
+    fn bounding_box(&self, _t_min: f32, _t_max: f32) -> Option<AABB> {
+        Some(AABB::new(self.p_min, self.p_max))
+    }
+}

rtiow/renderer/src/hitable.rs

@@ -1,7 +1,8 @@
-use std::sync::Arc;
+use std::{fmt::Debug, sync::Arc};
 
 use crate::{aabb::AABB, material::Material, ray::Ray, vec3::Vec3};
 
+#[derive(Debug)]
 pub struct HitRecord<'m> {
     pub t: f32,
     pub uv: (f32, f32),
@@ -10,7 +11,7 @@ pub struct HitRecord<'m> {
     pub material: &'m dyn Material,
 }
 
-pub trait Hit: Send + Sync {
+pub trait Hit: Send + Sync + Debug {
     fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord>;
     fn bounding_box(&self, t_min: f32, t_max: f32) -> Option<AABB>;
 }
@@ -23,3 +24,12 @@ impl Hit for Arc<dyn Hit> {
         (**self).bounding_box(t_min, t_max)
     }
 }
+
+impl Hit for Box<dyn Hit> {
+    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        (**self).hit(r, t_min, t_max)
+    }
+    fn bounding_box(&self, t_min: f32, t_max: f32) -> Option<AABB> {
+        (**self).bounding_box(t_min, t_max)
+    }
+}

rtiow/renderer/src/hitable_list.rs

@@ -7,7 +7,7 @@ use crate::{
     vec3::Vec3,
 };
 
-#[derive(Default)]
+#[derive(Debug, Default)]
 pub struct HitableList {
     list: Vec<Box<dyn Hit>>,
 }

rtiow/renderer/src/kdtree.rs

@@ -5,9 +5,11 @@ use log::info;
 use crate::{
     aabb::{surrounding_box, AABB},
     hitable::{Hit, HitRecord},
+    hitable_list::HitableList,
     ray::Ray,
 };
 
+#[derive(Debug)]
 pub enum KDTree {
     Leaf(Box<dyn Hit>),
     Branch {
@@ -102,6 +104,14 @@ impl KDTree {
             }),
             _ => panic!("Unreachable"),
         };
+        //info!("left_half {:?}", left_half);
+        //info!("right_half {:?}", right_half);
+        if left_half.is_empty() {
+            return KDTree::Leaf(Box::new(HitableList::new(right_half)));
+        };
+        if right_half.is_empty() {
+            return KDTree::Leaf(Box::new(HitableList::new(left_half)));
+        };
         KDTree::Branch {
             left: Box::new(KDTree::new(left_half, t_min, t_max)),
             right: Box::new(KDTree::new(right_half, t_min, t_max)),

rtiow/renderer/src/lib.rs

@@ -1,9 +1,13 @@
 pub mod aabb;
 pub mod bvh;
+pub mod bvh_triangles;
 pub mod camera;
+pub mod colors;
 pub mod constant_medium;
 pub mod cuboid;
+pub mod debug_hit;
 pub mod flip_normals;
+pub mod glowybox;
 pub mod hitable;
 pub mod hitable_list;
 pub mod human;
@@ -12,12 +16,15 @@ pub mod material;
 pub mod moving_sphere;
 pub mod noise;
 pub mod output;
+pub mod parser;
 pub mod ray;
 pub mod rect;
 pub mod renderer;
 pub mod rotate;
+pub mod scale;
 pub mod scenes;
 pub mod sphere;
 pub mod texture;
 pub mod translate;
+pub mod triangles;
 pub mod vec3;

rtiow/renderer/src/material.rs

@@ -1,4 +1,4 @@
-use std::sync::Arc;
+use std::{fmt::Debug, sync::Arc};
 
 use rand::{self, Rng};
 
@@ -20,14 +20,14 @@ fn random_in_unit_sphere() -> Vec3 {
     }
 }
 
-#[derive(Default)]
+#[derive(Default, Debug)]
 pub struct ScatterResponse {
     pub scattered: Ray,
     pub attenutation: Vec3,
     pub reflected: bool,
 }
 
-pub trait Material: Send + Sync {
+pub trait Material: Send + Sync + Debug {
     fn scatter(&self, r_in: &Ray, rec: &HitRecord) -> ScatterResponse;
     fn emitted(&self, _u: f32, _v: f32, _p: Vec3) -> Vec3 {
         Vec3::new(0., 0., 0.)
@@ -52,6 +52,7 @@ impl Material for Box<dyn Material> {
     }
 }
 
+#[derive(Clone, Debug)]
 pub struct Isotropic<T>
 where
     T: Texture,
@@ -82,6 +83,7 @@ where
     }
 }
 
+#[derive(Clone, Debug)]
 pub struct Lambertian<T>
 where
     T: Texture,
@@ -113,6 +115,7 @@ where
     }
 }
 
+#[derive(Clone, Debug)]
 pub struct Metal {
     albedo: Vec3,
     fuzzy: f32,
@@ -167,6 +170,7 @@ fn schlick(cosine: f32, ref_idx: f32) -> f32 {
     r0 + (1. - r0) * (1. - cosine).powf(5.)
 }
 
+#[derive(Clone, Debug)]
 pub struct Dielectric {
     ref_idx: f32,
 }
@@ -213,6 +217,7 @@ impl Material for Dielectric {
     }
 }
 
+#[derive(Debug)]
 pub struct DiffuseLight<T>
 where
     T: Texture,
@@ -246,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/moving_sphere.rs

@@ -7,6 +7,7 @@ use crate::{
     vec3::{dot, Vec3},
 };
 
+#[derive(Debug)]
 pub struct MovingSphere<M>
 where
     M: Material,

rtiow/renderer/src/noise/lode.rs

@@ -5,6 +5,7 @@ use rand;
 use crate::{noise::NoiseSource, vec3::Vec3};
 
 const NOISE_SIZE: usize = 128;
+#[derive(Debug)]
 pub struct Lode {
     // Using fixed array causes stack overflow.
     noise: Vec<Vec<Vec<f32>>>, //[[[f32; NOISE_SIZE]; NOISE_SIZE]; NOISE_SIZE],

rtiow/renderer/src/noise/mod.rs

@@ -1,13 +1,13 @@
 pub mod lode;
 pub mod perlin;
 
-use std::f32::consts::PI;
+use std::{f32::consts::PI, fmt::Debug};
 
 use serde_derive::Deserialize;
 
 use crate::vec3::Vec3;
 
-pub trait NoiseSource: Send + Sync {
+pub trait NoiseSource: Send + Sync + Debug {
     /// value returns noise on the interval [0., 1.).
     fn value(&self, p: Vec3) -> f32;
 

rtiow/renderer/src/noise/perlin.rs

@@ -8,6 +8,7 @@ use crate::{
     vec3::{dot, Vec3},
 };
 
+#[derive(Debug)]
 pub struct Perlin {
     ran_vec: Vec<Vec3>,
     perm_x: Vec<usize>,

rtiow/renderer/src/output.rs

@@ -2,6 +2,7 @@ use std::{
     collections::HashMap,
     fs::File,
     io::BufWriter,
+    net::TcpStream,
     path::Path,
     sync::{Arc, Mutex},
     time,
@@ -9,9 +10,9 @@ use std::{
 
 use chrono::Local;
 use image;
-use lazy_static::lazy_static;
 use log::info;
 use serde_derive::Serialize;
+use tev_client::{PacketCreateImage, PacketUpdateImage, TevClient};
 
 use crate::{renderer::Scene, vec3::Vec3};
 
@@ -24,10 +25,6 @@ pub const ADAPTIVE_DEPTH: &str = "adaptive_depth";
 // Grey scale showing rays cast per pixel.
 pub const RAYS_PER_PIXEL: &str = "rays_per_pixel";
 
-lazy_static! {
-    static ref DEBUGGER: Arc<Mutex<Debugger>> = Arc::new(Mutex::new(Debugger::new()));
-}
-
 #[derive(Serialize)]
 struct ImageMetadata {
     name: String,
@@ -74,43 +71,176 @@ impl Image {
     }
 }
 
-struct Debugger {
-    images: HashMap<String, (ImageType, Image)>,
+pub struct OutputManager {
+    images: Arc<Mutex<HashMap<String, (ImageType, Image)>>>,
+    tev_client: Option<Arc<Mutex<TevClient>>>,
 }
 
-impl Debugger {
-    fn new() -> Debugger {
-        Debugger {
-            images: HashMap::new(),
-        }
+impl OutputManager {
+    pub fn new(tev_addr: &Option<String>) -> std::io::Result<OutputManager> {
+        let tev_client = if let Some(addr) = tev_addr {
+            Some(Arc::new(Mutex::new(TevClient::wrap(TcpStream::connect(
+                addr,
+            )?))))
+        } else {
+            None
+        };
+        Ok(OutputManager {
+            images: Arc::new(Mutex::new(HashMap::new())),
+            tev_client,
+        })
     }
-}
 
-pub fn register_image(name: String, dimensions: (usize, usize), it: ImageType) {
-    let mut debugger = DEBUGGER.lock().unwrap();
-    debugger
-        .images
-        .insert(name, (it, Image::new(dimensions.0, dimensions.1)));
-}
+    pub fn register_image(&self, name: String, dimensions: (usize, usize), it: ImageType) {
+        let mut images = self.images.lock().unwrap();
+        images.insert(name.clone(), (it, Image::new(dimensions.0, dimensions.1)));
+        self.tev_client.clone().map(|c| {
+            c.lock().unwrap().send(PacketCreateImage {
+                image_name: &name,
+                grab_focus: false,
+                width: dimensions.0 as u32,
+                height: dimensions.1 as u32,
+                channel_names: &["R", "G", "B"],
+            })
+        });
+    }
 
-pub fn set_pixel(name: &str, x: usize, y: usize, pixel: Vec3) {
-    let mut debugger = DEBUGGER.lock().unwrap();
-    let (_it, img) = debugger
-        .images
-        .get_mut(name)
-        .unwrap_or_else(|| panic!("couldn't find image named '{}'", name));
-    let y_inv = img.h - y - 1;
-    img.put_pixel(x, y_inv, pixel);
-}
+    pub fn set_pixel(&self, name: &str, x: usize, y: usize, pixel: Vec3) {
+        let mut images = self.images.lock().unwrap();
+        let (_it, img) = images
+            .get_mut(name)
+            .unwrap_or_else(|| panic!("couldn't find image named '{}'", name));
+        let y_inv = img.h - y - 1;
+        img.put_pixel(x, y_inv, pixel);
+        self.tev_client.clone().map(|c| {
+            c.lock().unwrap().send(PacketUpdateImage {
+                image_name: &name,
+                grab_focus: false,
+                channel_names: &["R", "G", "B"],
+                channel_offsets: &[0, 1, 2],
+                channel_strides: &[0, 0, 0],
+                x: x as u32,
+                y: y_inv as u32,
+                width: 1,
+                height: 1,
+                data: &[pixel.x, pixel.y, pixel.z],
+            })
+        });
+    }
 
-pub fn set_pixel_grey(name: &str, x: usize, y: usize, grey: f32) {
-    let mut debugger = DEBUGGER.lock().unwrap();
-    let (_it, img) = debugger
-        .images
-        .get_mut(name)
-        .unwrap_or_else(|| panic!("couldn't find image named '{}'", name));
-    let y_inv = img.h - y - 1;
-    img.put_pixel(x, y_inv, [grey, grey, grey].into());
+    pub fn set_pixel_grey(&self, name: &str, x: usize, y: usize, grey: f32) {
+        let mut images = self.images.lock().unwrap();
+        let (_it, img) = images
+            .get_mut(name)
+            .unwrap_or_else(|| panic!("couldn't find image named '{}'", name));
+        let y_inv = img.h - y - 1;
+        img.put_pixel(x, y_inv, [grey, grey, grey].into());
+    }
+    pub fn write_images<P: AsRef<Path>>(
+        &self,
+        scene: &Scene,
+        render_time: time::Duration,
+        output_dir: P,
+    ) -> std::io::Result<()> {
+        let output_dir: &Path = output_dir.as_ref();
+        let now = Local::now();
+        let images = self.images.lock().unwrap();
+        // Write out images in consistent order.
+        let mut names = images.keys().collect::<Vec<_>>();
+        names.sort();
+        let mut image_metadata = Vec::new();
+        for name in &names {
+            let (it, img) = images.get(*name).unwrap();
+            let image = format!("{}.png", name);
+            let binary = format!("{}.json", name);
+            let ratio = img.w as f32 / img.h as f32;
+            let size = (img.w, img.h);
+            let image_path = output_dir.join(&image);
+            let binary_path = output_dir.join(&binary);
+            image_metadata.push(ImageMetadata {
+                name: name.to_string(),
+                image,
+                binary,
+                ratio,
+                size,
+                format: *it,
+            });
+            info!("Saving {}", image_path.to_string_lossy());
+            match it {
+                ImageType::RGB01 => {
+                    let mut out_img = image::RgbImage::new(img.w as u32, img.h as u32);
+                    out_img
+                        .enumerate_pixels_mut()
+                        .enumerate()
+                        .for_each(|(i, (_x, _y, p))| {
+                            let pixel = img.pix[i];
+                            *p = image::Rgb([
+                                (pixel[0] * 255.).min(255.) as u8,
+                                (pixel[1] * 255.).min(255.) as u8,
+                                (pixel[2] * 255.).min(255.) as u8,
+                            ])
+                        });
+                    out_img.save(image_path)?;
+                }
+                ImageType::Grey01 => {
+                    let mut out_img = image::GrayImage::new(img.w as u32, img.h as u32);
+                    out_img
+                        .enumerate_pixels_mut()
+                        .enumerate()
+                        .for_each(|(i, (_x, _y, p))| {
+                            let pixel = img.pix[i];
+                            *p = image::Luma([(pixel[0] * 255.).min(255.) as u8])
+                        });
+                    out_img.save(image_path)?;
+                }
+                ImageType::GreyNormalized => {
+                    let mut out_img = image::GrayImage::new(img.w as u32, img.h as u32);
+
+                    let max_val = img.pix.iter().map(|v| v.x).fold(0., f32::max);
+                    out_img
+                        .enumerate_pixels_mut()
+                        .enumerate()
+                        .for_each(|(i, (_x, _y, p))| {
+                            let pixel = img.pix[i];
+                            *p = image::Luma([(pixel[0] / max_val * 255.).min(255.) as u8])
+                        });
+                    out_img.save(image_path)?;
+                }
+            };
+            info!("Saving {}", binary_path.to_string_lossy());
+            let f = File::create(output_dir.join(binary_path))?;
+            let f = BufWriter::new(f);
+            match it {
+                ImageType::RGB01 => {
+                    serde_json::ser::to_writer(
+                        f,
+                        &img.pix
+                            .iter()
+                            .map(|v| [v.x, v.y, v.z])
+                            .collect::<Vec<[f32; 3]>>(),
+                    )?;
+                }
+                ImageType::Grey01 | ImageType::GreyNormalized => {
+                    serde_json::ser::to_writer(
+                        f,
+                        &img.pix.iter().map(|v| v.x).collect::<Vec<f32>>(),
+                    )?;
+                }
+            };
+        }
+        let f = File::create(output_dir.join("data.json"))?;
+        let f = BufWriter::new(f);
+        serde_json::ser::to_writer(
+            f,
+            &Data {
+                timestamp: now.timestamp(),
+                render_time_seconds: render_time.as_secs_f32(),
+                scene,
+                image_metadata,
+            },
+        )?;
+        Ok(())
+    }
 }
 
 trait ImageSaver {
@@ -118,105 +248,3 @@ trait ImageSaver {
     where
         Q: AsRef<Path> + Sized;
 }
-
-pub fn write_images<P: AsRef<Path>>(
-    scene: &Scene,
-    render_time: time::Duration,
-    output_dir: P,
-) -> std::io::Result<()> {
-    let output_dir: &Path = output_dir.as_ref();
-    let debugger = DEBUGGER.lock().unwrap();
-    let now = Local::now();
-    // Write out images in consistent order.
-    let mut names = debugger.images.keys().collect::<Vec<_>>();
-    names.sort();
-    let mut image_metadata = Vec::new();
-    for name in &names {
-        let (it, img) = debugger.images.get(*name).unwrap();
-        let image = format!("{}.png", name);
-        let binary = format!("{}.json", name);
-        let ratio = img.w as f32 / img.h as f32;
-        let size = (img.w, img.h);
-        let image_path = output_dir.join(&image);
-        let binary_path = output_dir.join(&binary);
-        image_metadata.push(ImageMetadata {
-            name: name.to_string(),
-            image,
-            binary,
-            ratio,
-            size,
-            format: *it,
-        });
-        info!("Saving {}", image_path.to_string_lossy());
-        match it {
-            ImageType::RGB01 => {
-                let mut out_img = image::RgbImage::new(img.w as u32, img.h as u32);
-                out_img
-                    .enumerate_pixels_mut()
-                    .enumerate()
-                    .for_each(|(i, (_x, _y, p))| {
-                        let pixel = img.pix[i];
-                        *p = image::Rgb([
-                            (pixel[0] * 255.).min(255.) as u8,
-                            (pixel[1] * 255.).min(255.) as u8,
-                            (pixel[2] * 255.).min(255.) as u8,
-                        ])
-                    });
-                out_img.save(image_path)?;
-            }
-            ImageType::Grey01 => {
-                let mut out_img = image::GrayImage::new(img.w as u32, img.h as u32);
-                out_img
-                    .enumerate_pixels_mut()
-                    .enumerate()
-                    .for_each(|(i, (_x, _y, p))| {
-                        let pixel = img.pix[i];
-                        *p = image::Luma([(pixel[0] * 255.).min(255.) as u8])
-                    });
-                out_img.save(image_path)?;
-            }
-            ImageType::GreyNormalized => {
-                let mut out_img = image::GrayImage::new(img.w as u32, img.h as u32);
-
-                let max_val = img.pix.iter().map(|v| v.x).fold(0., f32::max);
-                out_img
-                    .enumerate_pixels_mut()
-                    .enumerate()
-                    .for_each(|(i, (_x, _y, p))| {
-                        let pixel = img.pix[i];
-                        *p = image::Luma([(pixel[0] / max_val * 255.).min(255.) as u8])
-                    });
-                out_img.save(image_path)?;
-            }
-        };
-        info!("Saving {}", binary_path.to_string_lossy());
-        let f = File::create(output_dir.join(binary_path))?;
-        let f = BufWriter::new(f);
-        match it {
-            ImageType::RGB01 => {
-                serde_json::ser::to_writer(
-                    f,
-                    &img.pix
-                        .iter()
-                        .map(|v| [v.x, v.y, v.z])
-                        .collect::<Vec<[f32; 3]>>(),
-                )?;
-            }
-            ImageType::Grey01 | ImageType::GreyNormalized => {
-                serde_json::ser::to_writer(f, &img.pix.iter().map(|v| v.x).collect::<Vec<f32>>())?;
-            }
-        };
-    }
-    let f = File::create(output_dir.join("data.json"))?;
-    let f = BufWriter::new(f);
-    serde_json::ser::to_writer(
-        f,
-        &Data {
-            timestamp: now.timestamp(),
-            render_time_seconds: render_time.as_secs_f32(),
-            scene,
-            image_metadata,
-        },
-    )?;
-    Ok(())
-}

rtiow/renderer/src/parser.rs

@@ -0,0 +1,206 @@
+use crate::{
+    bvh_triangles::BVHTriangles,
+    camera::Camera,
+    cuboid::Cuboid,
+    hitable::Hit,
+    hitable_list::HitableList,
+    material::{Dielectric, DiffuseLight, Isotropic, Lambertian, Material, Metal},
+    renderer::Scene,
+    sphere::Sphere,
+    texture::{EnvMap, Texture},
+};
+use chrono::IsoWeek;
+use serde::Deserialize;
+use std::{collections::HashMap, fs::File, io::BufReader, path::PathBuf, sync::Arc};
+use stl::STL;
+use thiserror::Error;
+use vec3::Vec3;
+
+#[derive(Debug, Deserialize)]
+pub struct Config {
+    scene: SceneConfig,
+    camera: CameraConfig,
+    materials: Vec<MaterialConfig>,
+    hitables: Vec<HitableConfig>,
+    envmap: Option<EnvMapConfig>,
+}
+
+#[derive(Error, Debug)]
+pub enum ConfigError {
+    #[error("failed to load image")]
+    ImageError(#[from] image::ImageError),
+    #[error("failed to parser STL")]
+    STLError(#[from] stl::ParseError),
+    #[error("I/O error")]
+    IOError(#[from] std::io::Error),
+    #[error("duplication material named '{0}'")]
+    DuplicateMaterial(String),
+    #[error("unkown material named '{0}'")]
+    UnknownMaterial(String),
+}
+
+impl TryFrom<Config> for Scene {
+    type Error = ConfigError;
+
+    fn try_from(c: Config) -> Result<Scene, Self::Error> {
+        let mut materials = HashMap::new();
+        for mc in c.materials {
+            let v: Arc<dyn Material> = match mc.material {
+                Materials::Metal { albedo, fuzzy } => Arc::new(Metal::new(albedo, fuzzy)),
+                Materials::Dielectric { ref_idx } => Arc::new(Dielectric::new(ref_idx)),
+                Materials::DiffuseLight { texture } => Arc::new(DiffuseLight::new(texture)),
+                Materials::Isotropic { texture } => Arc::new(Isotropic::new(texture)),
+                Materials::Lambertian { texture } => Arc::new(Lambertian::new(texture)),
+            };
+            if materials.insert(mc.name.clone(), v).is_some() {
+                return Err(ConfigError::DuplicateMaterial(mc.name));
+            }
+        }
+
+        let hitables: Result<Vec<Box<dyn Hit>>, Self::Error> = c
+            .hitables
+            .into_iter()
+            .map(|hc| -> Result<Box<dyn Hit>, Self::Error> {
+                match hc.hitable {
+                    Hitables::Sphere { center, radius } => Ok(Box::new(Sphere::new(
+                        center,
+                        radius,
+                        Arc::clone(
+                            materials
+                                .get(&hc.material_name)
+                                .ok_or(ConfigError::UnknownMaterial(hc.material_name))?,
+                        ),
+                    ))),
+                    Hitables::Cuboid { min, max } => Ok(Box::new(Cuboid::new(
+                        min.into(),
+                        max.into(),
+                        Arc::clone(
+                            materials
+                                .get(&hc.material_name)
+                                .ok_or(ConfigError::UnknownMaterial(hc.material_name))?,
+                        ),
+                    ))),
+                    Hitables::STL { path, scale } => {
+                        let r = BufReader::new(File::open(path)?);
+                        let stl = STL::parse(r, false)?;
+                        Ok(Box::new(BVHTriangles::new(
+                            &stl,
+                            Arc::clone(
+                                materials
+                                    .get(&hc.material_name)
+                                    .ok_or(ConfigError::UnknownMaterial(hc.material_name))?,
+                            ),
+                            scale.unwrap_or(1.),
+                        )))
+                    }
+                }
+            })
+            .collect();
+        let hitables = hitables?;
+
+        let world: Box<dyn Hit> = Box::new(HitableList::new(hitables));
+        let mut env_map: Option<EnvMap> = None;
+
+        if let Some(em) = c.envmap {
+            let im = image::open(em.path)?.into_rgb();
+            env_map = Some(EnvMap::new(im));
+        };
+
+        let camera = make_camera(&c.camera, c.scene.width, c.scene.height);
+
+        let scene = Scene {
+            world,
+            camera,
+            env_map,
+            subsamples: c.scene.subsamples.unwrap_or(8),
+            adaptive_subsampling: c.scene.adaptive_subsampling,
+            num_threads: c.scene.num_threads,
+            width: c.scene.width,
+            height: c.scene.height,
+            global_illumination: c.scene.global_illumination.unwrap_or(true),
+        };
+        Ok(scene)
+    }
+}
+
+fn make_camera(cfg: &CameraConfig, width: usize, height: usize) -> Camera {
+    Camera::new(
+        cfg.lookfrom.into(),
+        cfg.lookat.into(),
+        Vec3::new(0., 1., 0.),
+        cfg.fov,
+        width as f32 / height as f32,
+        cfg.aperture,
+        cfg.focus_dist,
+        cfg.time_min,
+        cfg.time_max,
+    )
+}
+
+#[derive(Debug, Deserialize)]
+struct SceneConfig {
+    subsamples: Option<usize>,
+    adaptive_subsampling: Option<f32>,
+    num_threads: Option<usize>,
+    width: usize,
+    height: usize,
+    global_illumination: Option<bool>,
+}
+
+#[derive(Debug, Deserialize)]
+#[serde(tag = "type")]
+struct HitableConfig {
+    material_name: String,
+    #[serde(flatten)]
+    hitable: Hitables,
+}
+
+#[derive(Debug, Deserialize)]
+#[serde(tag = "type")]
+enum Hitables {
+    #[serde(rename = "sphere")]
+    Sphere { center: [f32; 3], radius: f32 },
+    #[serde(rename = "cuboid")]
+    Cuboid { min: [f32; 3], max: [f32; 3] },
+    #[serde(rename = "stl")]
+    STL { path: PathBuf, scale: Option<f32> },
+}
+
+#[derive(Debug, Deserialize)]
+struct MaterialConfig {
+    name: String,
+    #[serde(flatten)]
+    material: Materials,
+}
+
+#[derive(Debug, Deserialize)]
+#[serde(tag = "type")]
+enum Materials {
+    #[serde(rename = "metal")]
+    Metal { albedo: [f32; 3], fuzzy: f32 },
+    #[serde(rename = "dielectric")]
+    Dielectric { ref_idx: f32 },
+    // TODO(wathiede): these all take Textures, for now, only support RGB
+    #[serde(rename = "diffuse_light")]
+    DiffuseLight { texture: [f32; 3] },
+    #[serde(rename = "isotropic")]
+    Isotropic { texture: [f32; 3] },
+    #[serde(rename = "lambertian")]
+    Lambertian { texture: [f32; 3] },
+}
+
+#[derive(Debug, Deserialize)]
+pub struct CameraConfig {
+    lookfrom: [f32; 3],
+    lookat: [f32; 3],
+    fov: f32,
+    aperture: f32,
+    focus_dist: f32,
+    time_min: f32,
+    time_max: f32,
+}
+
+#[derive(Debug, Deserialize)]
+struct EnvMapConfig {
+    path: PathBuf,
+}

rtiow/renderer/src/ray.rs

@@ -5,6 +5,9 @@ pub struct Ray {
     pub origin: Vec3,
     pub direction: Vec3,
     pub time: f32,
+
+    // Precache 1/direction, a single ray intersects multiple AABB's, and divides are more
+    // expensive than multiplies.
     pub inv_direction: Vec3,
     pub sign: [usize; 3],
 }
@@ -14,7 +17,7 @@ impl Ray {
     where
         V: Into<Vec3>,
     {
-        let direction = direction.into();
+        let direction: Vec3 = direction.into();
         let origin = origin.into();
         let inv = 1. / direction;
         Ray {

rtiow/renderer/src/rect.rs

@@ -8,6 +8,7 @@ use crate::{
     vec3::Vec3,
 };
 
+#[derive(Debug)]
 pub struct XYRect<M>
 where
     M: Material,
@@ -69,6 +70,7 @@ where
     }
 }
 
+#[derive(Debug)]
 pub struct XZRect<M>
 where
     M: Material,
@@ -130,6 +132,7 @@ where
     }
 }
 
+#[derive(Debug)]
 pub struct YZRect<M>
 where
     M: Material,

rtiow/renderer/src/renderer.rs

@@ -1,4 +1,5 @@
 use std::{
+    collections::HashMap,
     fmt,
     ops::{AddAssign, Range},
     path::{Path, PathBuf},
@@ -7,7 +8,8 @@ use std::{
         mpsc::{sync_channel, Receiver, SyncSender},
         Arc, Mutex,
     },
-    thread, time,
+    thread,
+    time::{Duration, Instant},
 };
 
 use core_affinity;
@@ -21,16 +23,19 @@ use crate::{
     camera::Camera,
     hitable::Hit,
     human,
-    material::Lambertian,
+    material::{Lambertian, Material},
     output,
+    output::OutputManager,
     ray::Ray,
     scenes,
     sphere::Sphere,
     texture::{ConstantTexture, EnvMap},
     vec3::Vec3,
 };
+use strum::{EnumString, EnumVariantNames};
 
-#[derive(Debug)]
+#[derive(Debug, EnumString, EnumVariantNames, strum::Display)]
+#[strum(serialize_all = "snake_case")]
 pub enum Model {
     BVH,
     Bench,
@@ -41,8 +46,11 @@ pub enum Model {
     Mandelbrot,
     PerlinDebug,
     Spheramid,
+    Stltest,
     Test,
+    Tron,
     Tutorial,
+    Dragon,
 }
 
 impl Model {
@@ -53,11 +61,14 @@ impl Model {
             Model::Book => scenes::book::new(opt),
             Model::CornellBox => scenes::cornell_box::new(opt),
             Model::CornellSmoke => scenes::cornell_smoke::new(opt),
+            Model::Dragon => scenes::dragon::new(opt),
             Model::Final => scenes::final_scene::new(opt),
             Model::Mandelbrot => scenes::mandelbrot::new(opt),
             Model::PerlinDebug => scenes::perlin_debug::new(opt),
             Model::Spheramid => scenes::spheramid::new(opt),
+            Model::Stltest => scenes::stltest::new(opt),
             Model::Test => scenes::test::new(opt),
+            Model::Tron => scenes::tron::new(opt),
             Model::Tutorial => scenes::tutorial::new(opt),
         }
     }
@@ -72,44 +83,6 @@ impl fmt::Display for ModelParseError {
     }
 }
 
-impl str::FromStr for Model {
-    type Err = ModelParseError;
-    fn from_str(s: &str) -> std::result::Result<Self, Self::Err> {
-        match s {
-            "bench" => Ok(Model::Bench),
-            "book" => Ok(Model::Book),
-            "bvh" => Ok(Model::BVH),
-            "cornell_box" => Ok(Model::CornellBox),
-            "cornell_smoke" => Ok(Model::CornellSmoke),
-            "final" => Ok(Model::Final),
-            "mandelbrot" => Ok(Model::Mandelbrot),
-            "perlin_debug" => Ok(Model::PerlinDebug),
-            "spheramid" => Ok(Model::Spheramid),
-            "test" => Ok(Model::Test),
-            "tutorial" => Ok(Model::Tutorial),
-            _ => Err(ModelParseError(s.to_owned())),
-        }
-    }
-}
-
-impl std::string::ToString for Model {
-    fn to_string(&self) -> String {
-        match self {
-            Model::BVH => "bvh".to_string(),
-            Model::Bench => "bench".to_string(),
-            Model::Book => "book".to_string(),
-            Model::CornellBox => "cornell_box".to_string(),
-            Model::CornellSmoke => "cornell_smoke".to_string(),
-            Model::Final => "final".to_string(),
-            Model::Mandelbrot => "mandelbrot".to_string(),
-            Model::PerlinDebug => "perlin_debug".to_string(),
-            Model::Spheramid => "spheramid".to_string(),
-            Model::Test => "test".to_string(),
-            Model::Tutorial => "tutorial".to_string(),
-        }
-    }
-}
-
 #[derive(Debug, StructOpt)]
 #[structopt(name = "tracer", about = "An experimental ray tracer.")]
 pub struct Opt {
@@ -125,19 +98,29 @@ pub struct Opt {
     /// Sub-samples per pixel
     #[structopt(short = "s", long = "subsample", default_value = "8")]
     pub subsamples: usize,
-    /// Select scene to render, one of: "bench", "book", "tutorial", "bvh", "test", "cornell_box",
-    /// "cornell_smoke", "perlin_debug", "final"
-    #[structopt(long = "model", default_value = "book")]
-    pub model: Model,
+    /// Select scene to render.
+    #[structopt(long = "model")]
+    pub model: Option<Model>,
+    /// Toml config describing scene.
+    #[structopt(long = "config")]
+    pub config: Option<PathBuf>,
     /// Path to store pprof profile data, i.e. /tmp/cpuprofile.pprof
     #[structopt(long = "pprof", parse(from_os_str))]
     pub pprof: Option<PathBuf>,
     /// Use acceleration data structure, may be BVH or kd-tree depending on scene.
     #[structopt(long = "use_accel")]
     pub use_accel: bool,
+    /// Host:port of running tev instance.
+    #[structopt(long = "tev_addr")]
+    pub tev_addr: Option<String>,
 
     /// Output directory
-    #[structopt(parse(from_os_str), default_value = "/tmp/tracer")]
+    #[structopt(
+        short = "o",
+        long = "output",
+        parse(from_os_str),
+        default_value = "/tmp/tracer"
+    )]
     pub output: PathBuf,
 }
 
@@ -149,18 +132,20 @@ pub fn opt_hash(opt: &Opt) -> String {
         opt.height,
         opt.subsamples,
         opt.pprof.is_some(),
-        opt.model.to_string(),
+        opt.model.as_ref().unwrap().to_string(),
         opt.use_accel,
         opt.output.display().to_string().replace('/', "_")
     )
 }
 
 // TODO(wathiede): implement the skips and then the renderer could use json as an input file type.
-#[derive(Serialize)]
+#[derive(Debug, Serialize)]
 #[serde(rename_all = "camelCase")]
 pub struct Scene {
     #[serde(skip)]
     pub world: Box<dyn Hit>,
+    //#[serde(skip)]
+    //pub materials: HashMap<String, Box<dyn Material>>,
     #[serde(skip)]
     pub camera: Camera,
     pub subsamples: usize,
@@ -265,6 +250,7 @@ fn trace_pixel_adaptive(
     x_range: Range<f32>,
     y_range: Range<f32>,
     scene: &Scene,
+    output: &OutputManager,
 ) -> (Vec3, usize) {
     let w = scene.width as f32;
     let h = scene.height as f32;
@@ -279,7 +265,7 @@ fn trace_pixel_adaptive(
         &scene.env_map,
     );
     if depth == 0 {
-        output::set_pixel(output::ADAPTIVE_DEPTH, x, y, [1., 0., 0.].into());
+        output.set_pixel(output::ADAPTIVE_DEPTH, x, y, [1., 0., 0.].into());
         return (center, rays);
     }
     // t = top
@@ -321,6 +307,7 @@ fn trace_pixel_adaptive(
             x_range.start..x_mid,
             y_range.start..y_mid,
             scene,
+            output,
         );
         let tr = trace_pixel_adaptive(
             depth - 1,
@@ -330,6 +317,7 @@ fn trace_pixel_adaptive(
             x_mid..x_range.end,
             y_range.start..y_mid,
             scene,
+            output,
         );
         let bl = trace_pixel_adaptive(
             depth - 1,
@@ -339,6 +327,7 @@ fn trace_pixel_adaptive(
             x_range.start..x_mid,
             y_mid..y_range.end,
             scene,
+            output,
         );
         let br = trace_pixel_adaptive(
             depth - 1,
@@ -348,13 +337,14 @@ fn trace_pixel_adaptive(
             x_mid..x_range.end,
             y_mid..y_range.end,
             scene,
+            output,
         );
         let pixel = (tl.0 + tr.0 + bl.0 + br.0) / 4.;
         let rays = tl.1 + tr.1 + bl.1 + br.1;
         (pixel, rays)
     } else {
         if depth == MAX_ADAPTIVE_DEPTH {
-            output::set_pixel(output::ADAPTIVE_DEPTH, x, y, [0., 1., 0.].into());
+            output.set_pixel(output::ADAPTIVE_DEPTH, x, y, [0., 1., 0.].into());
         }
         (corners, rays)
     }
@@ -390,27 +380,37 @@ impl AddAssign for RenderStats {
 }
 
 fn progress(
+    start_time: Instant,
     last_stat: &RenderStats,
     current_stat: &RenderStats,
-    time_diff: time::Duration,
+    time_diff: Duration,
     pixel_total: usize,
 ) -> String {
     let human = human::Formatter::new();
     let pixel_diff = current_stat.pixels - last_stat.pixels;
     let ray_diff = current_stat.rays - last_stat.rays;
+    let now = Instant::now();
+    let start_diff = now - start_time;
+    let ratio = current_stat.pixels as f32 / pixel_total as f32;
+    let percent = ratio * 100.;
+    let elapsed = start_diff.as_secs_f32();
+    let total = elapsed * (1. / ratio);
+    let eta = total - elapsed;
     format!(
-        "{:7} / {:7}pixels ({:2}%) {:7}pixels/s {:7}rays/s",
+        "{:7} / {:7}pixels ({:2.0}%) {:7}pixels/s {:7}rays/s eta {:.0}s",
         human.format(current_stat.pixels as f64),
         human.format(pixel_total as f64),
-        100 * current_stat.pixels / pixel_total,
+        percent,
         human.format(pixel_diff as f64 / time_diff.as_secs_f64()),
-        human.format(ray_diff as f64 / time_diff.as_secs_f64())
+        human.format(ray_diff as f64 / time_diff.as_secs_f64()),
+        eta
     )
 }
 
 enum Request {
     Pixel { x: usize, y: usize },
     Line { width: usize, y: usize },
+    // TODO(wathiede): add Cohort that does 4x4 or 8x8 pixel chunks.
 }
 
 enum Response {
@@ -427,7 +427,7 @@ enum Response {
     },
 }
 
-fn render_pixel(scene: &Scene, x: usize, y: usize) -> (Vec3, usize) {
+fn render_pixel(scene: &Scene, x: usize, y: usize, output: &OutputManager) -> (Vec3, usize) {
     let (pixel, rays) = if let Some(threshold) = scene.adaptive_subsampling {
         trace_pixel_adaptive(
             MAX_ADAPTIVE_DEPTH,
@@ -437,6 +437,7 @@ fn render_pixel(scene: &Scene, x: usize, y: usize) -> (Vec3, usize) {
             0.0..1.0,
             0.0..1.0,
             scene,
+            output,
         )
     } else {
         let (pixel, rays) = (0..scene.subsamples)
@@ -445,7 +446,7 @@ fn render_pixel(scene: &Scene, x: usize, y: usize) -> (Vec3, usize) {
                 ([0., 0., 0.].into(), 0),
                 |(p1, r1): (Vec3, usize), (p2, r2): (Vec3, usize)| ((p1 + p2), (r1 + r2)),
             );
-        output::set_pixel_grey(output::RAYS_PER_PIXEL, x, y, rays as f32);
+        output.set_pixel_grey(output::RAYS_PER_PIXEL, x, y, rays as f32);
         (pixel / scene.subsamples as f32, rays)
     };
     // Gamma correct, use gamma 2 correction, which is 1/gamma where gamma=2 which is 1/2 or
@@ -461,6 +462,7 @@ fn render_worker(
     scene: &Scene,
     input_chan: Arc<Mutex<Receiver<Request>>>,
     output_chan: &SyncSender<Response>,
+    output: &OutputManager,
 ) {
     loop {
         let job = { input_chan.lock().unwrap().recv() };
@@ -475,7 +477,7 @@ fn render_worker(
                     let batch = false;
                     if batch {
                         let (pixels, rays): (Vec<Vec3>, Vec<usize>) = (0..width)
-                            .map(|x| render_pixel(scene, x, y))
+                            .map(|x| render_pixel(scene, x, y, output))
                             .collect::<Vec<(_, _)>>()
                             .into_iter()
                             .unzip();
@@ -492,7 +494,7 @@ fn render_worker(
                             .expect("failed to send pixel response");
                     } else {
                         (0..width).for_each(|x| {
-                            let (pixel, rays) = render_pixel(scene, x, y);
+                            let (pixel, rays) = render_pixel(scene, x, y, output);
                             output_chan
                                 .send(Response::Pixel {
                                     x,
@@ -506,7 +508,7 @@ fn render_worker(
                 }
                 Request::Pixel { x, y } => {
                     trace!("tid {} x {} y {}", tid, x, y);
-                    let (pixel, rays) = render_pixel(scene, x, y);
+                    let (pixel, rays) = render_pixel(scene, x, y, output);
                     output_chan
                         .send(Response::Pixel {
                             x,
@@ -521,8 +523,19 @@ 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);
+/*
+lazy_static! {
+    static ref DEBUGGER: Arc<Mutex<OutputManager>> = Arc::new(Mutex::new(OutputManager::new()));
+}
+*/
+
+pub fn render(
+    scene: Scene,
+    output_dir: &Path,
+    tev_addr: &Option<String>,
+) -> std::result::Result<(), std::io::Error> {
+    // 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);
 
@@ -536,20 +549,23 @@ pub fn render(scene: Scene, output_dir: &Path) -> std::result::Result<(), std::i
     } else {
         core_ids
     };
+    let output = output::OutputManager::new(tev_addr)?;
+    let output = Arc::new(output);
+
     info!("Creating {} render threads", core_ids.len());
-    output::register_image(
+    output.register_image(
         output::MAIN_IMAGE.to_string(),
         (scene.width, scene.height),
         output::ImageType::RGB01,
     );
     if scene.adaptive_subsampling.is_some() {
-        output::register_image(
+        output.register_image(
             output::ADAPTIVE_DEPTH.to_string(),
             (scene.width, scene.height),
             output::ImageType::RGB01,
         );
     }
-    output::register_image(
+    output.register_image(
         output::RAYS_PER_PIXEL.to_string(),
         (scene.width, scene.height),
         output::ImageType::GreyNormalized,
@@ -563,19 +579,19 @@ pub fn render(scene: Scene, output_dir: &Path) -> std::result::Result<(), std::i
             let s = sync::Arc::clone(&scene);
             let pixel_req_rx = pixel_req_rx.clone();
             let pixel_resp_tx = pixel_resp_tx.clone();
+            let output = sync::Arc::clone(&output);
             thread::spawn(move || {
                 core_affinity::set_for_current(id);
-                render_worker(i, &s, pixel_req_rx, &pixel_resp_tx);
+                render_worker(i, &s, pixel_req_rx, &pixel_resp_tx, &output);
             })
         })
         .collect::<Vec<_>>();
     drop(pixel_req_rx);
     drop(pixel_resp_tx);
 
-    let start_time = time::Instant::now();
     let (w, h) = (scene.width, scene.height);
     handles.push(thread::spawn(move || {
-        let batch_line_requests = true;
+        let batch_line_requests = false;
         if batch_line_requests {
             for y in 0..h {
                 pixel_req_tx
@@ -597,29 +613,36 @@ pub fn render(scene: Scene, output_dir: &Path) -> std::result::Result<(), std::i
     info!("Rendering with {} subsamples", scene.subsamples);
 
     let pixel_total = scene.width * scene.height;
-    let mut last_time = time::Instant::now();
+    let mut last_time = Instant::now();
     let mut last_stat: RenderStats = Default::default();
     let mut current_stat: RenderStats = Default::default();
+    let render_start_time = Instant::now();
     for resp in pixel_resp_rx {
         match resp {
             Response::Pixel { x, y, pixel, rs } => {
                 current_stat += rs;
-                output::set_pixel(output::MAIN_IMAGE, x, y, pixel);
+                output.set_pixel(output::MAIN_IMAGE, x, y, pixel);
             }
             Response::Line { y, pixels, rs } => {
                 current_stat += rs;
                 for (x, pixel) in pixels.iter().enumerate() {
-                    output::set_pixel(output::MAIN_IMAGE, x, y, *pixel);
+                    output.set_pixel(output::MAIN_IMAGE, x, y, *pixel);
                 }
             }
         }
 
-        let now = time::Instant::now();
+        let now = Instant::now();
         let time_diff = now - last_time;
-        if time_diff > time::Duration::from_secs(1) {
-            info!(
+        if time_diff > Duration::from_secs(5) {
+            println!(
                 "{}",
-                progress(&last_stat, &current_stat, time_diff, pixel_total)
+                progress(
+                    render_start_time,
+                    &last_stat,
+                    &current_stat,
+                    time_diff,
+                    pixel_total
+                )
             );
             last_stat = current_stat;
             last_time = now;
@@ -628,12 +651,18 @@ pub fn render(scene: Scene, output_dir: &Path) -> std::result::Result<(), std::i
     for thr in handles {
         thr.join().expect("thread join");
     }
-    let time_diff = time::Instant::now() - start_time;
-    info!(
-        "Runtime {} seconds {}",
+    let time_diff = Instant::now() - render_start_time;
+    println!(
+        "Render {} seconds {}",
         time_diff.as_secs_f32(),
-        progress(&Default::default(), &current_stat, time_diff, pixel_total)
+        progress(
+            render_start_time,
+            &Default::default(),
+            &current_stat,
+            time_diff,
+            pixel_total
+        )
     );
 
-    output::write_images(&scene, time_diff, output_dir)
+    output.write_images(&scene, time_diff, output_dir)
 }

rtiow/renderer/src/rotate.rs

@@ -7,6 +7,7 @@ use crate::{
     vec3::Vec3,
 };
 
+#[derive(Debug)]
 pub struct RotateY<H>
 where
     H: Hit,

rtiow/renderer/src/scale.rs

@@ -0,0 +1,47 @@
+use crate::{
+    aabb::AABB,
+    hitable::{Hit, HitRecord},
+    ray::Ray,
+};
+
+#[derive(Debug)]
+pub struct Scale<H>
+where
+    H: Hit,
+{
+    hitable: H,
+    scale: f32,
+}
+
+impl<H> Scale<H>
+where
+    H: Hit,
+{
+    pub fn new(hitable: H, scale: f32) -> Scale<H> {
+        Scale { hitable, scale }
+    }
+}
+
+impl<H> Hit for Scale<H>
+where
+    H: Hit,
+{
+    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        let moved_r = Ray::new(r.origin / self.scale, r.direction, r.time);
+        if let Some(rec) = self.hitable.hit(moved_r, t_min, t_max) {
+            return Some(HitRecord {
+                p: rec.p * self.scale,
+                t: rec.t * self.scale,
+                ..rec
+            });
+        }
+        None
+    }
+
+    fn bounding_box(&self, t_min: f32, t_max: f32) -> Option<AABB> {
+        if let Some(bbox) = self.hitable.bounding_box(t_min, t_max) {
+            return Some(AABB::new(bbox.min() * self.scale, bbox.max() * self.scale));
+        }
+        None
+    }
+}

rtiow/renderer/src/scenes/book.rs

@@ -52,6 +52,7 @@ pub fn new(opt: &Opt) -> Scene {
         height: opt.height,
         global_illumination: true,
         env_map: Some(EnvMap::new(skybox)),
+        ..Default::default()
     }
 }
 

rtiow/renderer/src/scenes/dragon.rs

@@ -0,0 +1,116 @@
+use std::{
+    f32::consts::PI,
+    io::{BufReader, Cursor},
+};
+
+use stl::STL;
+
+use crate::{
+    bvh_triangles::BVHTriangles,
+    camera::Camera,
+    colors::generate_rainbow,
+    hitable::Hit,
+    hitable_list::HitableList,
+    kdtree::KDTree,
+    material::{Lambertian, Metal},
+    renderer::{Opt, Scene},
+    rotate::RotateY,
+    scale::Scale,
+    sphere::Sphere,
+    texture::{ConstantTexture, EnvMap},
+    translate::Translate,
+    vec3::Vec3,
+};
+
+pub fn new(opt: &Opt) -> Scene {
+    let lookfrom = Vec3::new(0., 80., 80.);
+    let lookat = Vec3::new(0., 0., 0.);
+    let dist_to_focus = 10.0;
+    let aperture = 0.0;
+    let time_min = 0.;
+    let time_max = 1.;
+    let camera = Camera::new(
+        lookfrom,
+        lookat,
+        Vec3::new(0., 1., 0.),
+        45.,
+        opt.width as f32 / opt.height as f32,
+        aperture,
+        dist_to_focus,
+        time_min,
+        time_max,
+    );
+    //let dragon_material = Dielectric::new(1.5);
+    let dragon_material = Metal::new(Vec3::new(0.6, 0.6, 0.6), 0.0);
+    //let dragon_material = Lambertian::new(ConstantTexture::new(Vec3::new(1.0, 1.0, 0.2)));
+
+    let ground_color = if opt.use_accel {
+        ConstantTexture::new(Vec3::new(1.0, 0.4, 0.4))
+    } else {
+        ConstantTexture::new(Vec3::new(0.4, 0.4, 0.4))
+    };
+
+    let stl_cube = STL::parse(
+        BufReader::new(Cursor::new(include_bytes!("../../stls/dragon.stl"))),
+        false,
+    )
+    .expect("failed to parse cube");
+    let _light_size = 50.;
+    let _light_height = 200.;
+    let sphere_radius = 5.;
+    let circle_radius = 40.;
+    let num_spheres = 16;
+    let palette = generate_rainbow(num_spheres);
+    let spheres: Vec<Box<dyn Hit>> = (0..num_spheres)
+        .map(|i| (i, i as f32, num_spheres as f32))
+        .map(|(idx, idx_f, n)| (idx, idx_f * 2. * PI / n))
+        .map(|(idx, rad)| -> Box<dyn Hit> {
+            let x = circle_radius * rad.cos();
+            let y = 4. * sphere_radius;
+            let z = circle_radius * rad.sin();
+            let c = palette[idx];
+            Box::new(Sphere::new(
+                [x, y, z],
+                sphere_radius,
+                Lambertian::new(ConstantTexture::new(c)),
+            ))
+        })
+        .collect();
+    let mut objects: Vec<Box<dyn Hit>> = vec![
+        Box::new(Sphere::new(
+            Vec3::new(0., 0.1, -0.5),
+            0.1,
+            Lambertian::new([0., 1., 1.]),
+        )),
+        // Earth sized sphere
+        //Box::new(Sphere::new( Vec3::new(0., -10000., 0.), 10000., Lambertian::new(ground_color),)),
+        // STL Mesh
+        Box::new(crate::debug_hit::DebugHit::new(RotateY::new(
+            Translate::new(
+                BVHTriangles::new(&stl_cube, dragon_material, 250.),
+                [0., -10., 0.],
+            ),
+            180.,
+        ))),
+    ];
+    objects.extend(spheres);
+
+    let world: Box<dyn Hit> = if opt.use_accel {
+        Box::new(KDTree::new(objects, time_min, time_max))
+    } else {
+        Box::new(HitableList::new(objects))
+    };
+    let skybox_bytes = include_bytes!("../../images/envmap.jpg");
+    let skybox = image::load_from_memory(skybox_bytes).unwrap().to_rgb();
+    Scene {
+        camera,
+        world,
+        subsamples: opt.subsamples,
+        num_threads: opt.num_threads,
+        width: opt.width,
+        height: opt.height,
+        global_illumination: true,
+        env_map: Some(EnvMap::new(skybox)),
+        ..Default::default()
+    }
+}

rtiow/renderer/src/scenes/mod.rs

@@ -3,9 +3,12 @@ pub mod book;
 pub mod bvh;
 pub mod cornell_box;
 pub mod cornell_smoke;
+pub mod dragon;
 pub mod final_scene;
 pub mod mandelbrot;
 pub mod perlin_debug;
 pub mod spheramid;
+pub mod stltest;
 pub mod test;
+pub mod tron;
 pub mod tutorial;

rtiow/renderer/src/scenes/spheramid.rs

@@ -112,5 +112,6 @@ pub fn new(opt: &Opt) -> Scene {
         height: opt.height,
         global_illumination: true,
         env_map: Some(EnvMap::new(skybox)),
+        ..Default::default()
     }
 }

rtiow/renderer/src/scenes/stltest.rs

@@ -0,0 +1,122 @@
+use std::{
+    io::{BufReader, Cursor},
+    sync::Arc,
+};
+
+use stl::STL;
+
+use crate::{
+    bvh_triangles::BVHTriangles,
+    camera::Camera,
+    cuboid::Cuboid,
+    hitable::Hit,
+    hitable_list::HitableList,
+    kdtree::KDTree,
+    material::{Dielectric, Lambertian, Metal},
+    renderer::{Opt, Scene},
+    sphere::Sphere,
+    texture::{ConstantTexture, EnvMap},
+    translate::Translate,
+    vec3::Vec3,
+};
+
+pub fn new(opt: &Opt) -> Scene {
+    let lookfrom = Vec3::new(0., 40., -100.);
+    let lookat = Vec3::new(0., 10., 0.);
+    let dist_to_focus = 10.0;
+    let aperture = 0.0;
+    let time_min = 0.;
+    let time_max = 1.;
+    let camera = Camera::new(
+        lookfrom,
+        lookat,
+        Vec3::new(0., 1., 0.),
+        45.,
+        opt.width as f32 / opt.height as f32,
+        aperture,
+        dist_to_focus,
+        time_min,
+        time_max,
+    );
+    let ground_color = if opt.use_accel {
+        ConstantTexture::new(Vec3::new(1.0, 0.4, 0.4))
+    } else {
+        ConstantTexture::new(Vec3::new(0.4, 0.4, 0.4))
+    };
+
+    let glass = Dielectric::new(1.5);
+    let metal = Metal::new(Vec3::new(0.8, 0.8, 0.8), 0.2);
+    let red = Lambertian::new(ConstantTexture::new(Vec3::new(1.0, 0.2, 0.2)));
+
+    //let box_material = glass;
+    let _ = glass;
+    let _ = metal;
+    let _ = red;
+
+    let stl_cube = STL::parse(
+        BufReader::new(Cursor::new(include_bytes!("../../stls/cube.stl"))),
+        false,
+    )
+    .expect("failed to parse cube");
+    let objects: Vec<Box<dyn Hit>> = vec![
+        // Light from above - white
+        // Earth sized sphere
+        Box::new(Sphere::new(
+            Vec3::new(0., -10000., 0.),
+            10000.,
+            Lambertian::new(ground_color),
+        )),
+        Box::new(Sphere::new(
+            Vec3::new(0., 20., -40.),
+            10.,
+            Lambertian::new(ConstantTexture::new(Vec3::new(1., 0.2, 1.))),
+        )),
+        // Blue sphere
+        Box::new(Sphere::new(
+            Vec3::new(0., 20., 40.),
+            20.,
+            Lambertian::new(ConstantTexture::new(Vec3::new(0.2, 0.2, 1.))),
+        )),
+        Box::new(Sphere::new(
+            Vec3::new(40., 20., 40.),
+            20.,
+            //Metal::new(Vec3::new(0.8, 0.8, 0.8), 0.2),
+            Lambertian::new(ConstantTexture::new(Vec3::new(0.2, 1.0, 0.2))),
+        )),
+        Box::new(Sphere::new(
+            Vec3::new(-40., 20., 40.),
+            20.,
+            //Metal::new(Vec3::new(0.8, 0.8, 0.8), 0.2),
+            Lambertian::new(ConstantTexture::new(Vec3::new(1.0, 0.2, 0.2))),
+        )),
+        // STL Mesh
+        Box::new(Translate::new(
+            BVHTriangles::new(&stl_cube, glass, 1.),
+            [0., 10., 0.],
+        )),
+        //Box::new(BVHTriangles::new(&stl_cube, box_material.clone())),
+        Box::new(Cuboid::new(
+            [-20., 0., 0.].into(),
+            [0., 20., 20.].into(),
+            Arc::new(red),
+        )),
+    ];
+    let world: Box<dyn Hit> = if opt.use_accel {
+        Box::new(KDTree::new(objects, time_min, time_max))
+    } else {
+        Box::new(HitableList::new(objects))
+    };
+    let skybox_bytes = include_bytes!("../../images/envmap.jpg");
+    let skybox = image::load_from_memory(skybox_bytes).unwrap().to_rgb();
+    Scene {
+        camera,
+        world,
+        subsamples: opt.subsamples,
+        num_threads: opt.num_threads,
+        width: opt.width,
+        height: opt.height,
+        global_illumination: true,
+        env_map: Some(EnvMap::new(skybox)),
+        ..Default::default()
+    }
+}

rtiow/renderer/src/scenes/tron.rs

@@ -0,0 +1,124 @@
+use std::sync::Arc;
+
+use log::info;
+
+use crate::{
+    camera::Camera,
+    glowybox::Glowybox,
+    hitable::Hit,
+    hitable_list::HitableList,
+    kdtree::KDTree,
+    material::{DiffuseLight, Lambertian},
+    rect::{XYRect, XZRect, YZRect},
+    renderer::{Opt, Scene},
+    sphere::Sphere,
+    texture::ConstantTexture,
+    vec3::Vec3,
+};
+
+pub fn new(opt: &Opt) -> Scene {
+    let lookfrom = Vec3::new(20., 20., 20.);
+    let lookat = Vec3::new(0., 1., 0.);
+    let dist_to_focus = 10.0;
+    let aperture = 0.0;
+    let time_min = 0.;
+    let time_max = 1.;
+    let camera = Camera::new(
+        lookfrom,
+        lookat,
+        Vec3::new(0., 1., 0.),
+        20.,
+        opt.width as f32 / opt.height as f32,
+        aperture,
+        dist_to_focus,
+        time_min,
+        time_max,
+    );
+    let ground_color = if opt.use_accel {
+        ConstantTexture::new(Vec3::new(1.0, 0.4, 0.4))
+    } else {
+        ConstantTexture::new(Vec3::new(0.4, 1.0, 0.4))
+    };
+
+    let lights: Vec<Box<dyn Hit>> = vec![
+        Box::new(XZRect::new(
+            -100.,
+            100.,
+            -100.,
+            1000.,
+            60.,
+            DiffuseLight::new(ConstantTexture::new(Vec3::new(1., 1., 1.))),
+        )),
+        Box::new(YZRect::new(
+            1.,
+            3.,
+            -1.,
+            1.,
+            4.,
+            DiffuseLight::new(ConstantTexture::new(Vec3::new(4., 0., 4.))),
+        )),
+        Box::new(YZRect::new(
+            1.,
+            3.,
+            -1.,
+            1.,
+            -4.,
+            DiffuseLight::new(ConstantTexture::new(Vec3::new(0., 4., 0.))),
+        )),
+        Box::new(XZRect::new(
+            -1.,
+            1.,
+            -1.,
+            1.,
+            6.,
+            DiffuseLight::new(ConstantTexture::new(Vec3::new(4., 4., 0.))),
+        )),
+        Box::new(XYRect::new(
+            -1.,
+            1.,
+            1.,
+            3.,
+            -4.,
+            DiffuseLight::new(ConstantTexture::new(Vec3::new(0., 0., 4.))),
+        )),
+        Box::new(XYRect::new(
+            -1.,
+            1.,
+            1.,
+            3.,
+            4.,
+            DiffuseLight::new(ConstantTexture::new(Vec3::new(0., 4., 4.))),
+        )),
+    ];
+    let mut objects: Vec<Box<dyn Hit>> = vec![
+        // Earth sized sphere
+        Box::new(Sphere::new(
+            Vec3::new(0., -1010., 0.),
+            1000.,
+            Lambertian::new(ground_color),
+        )),
+        Box::new(Glowybox::new(
+            [-4., -4., -4.].into(),
+            [4., 4., 4.].into(),
+            0.01,
+            Arc::new(Lambertian::new(ConstantTexture::new([0., 0., 0.]))),
+            Arc::new(DiffuseLight::new(ConstantTexture::new([100., 0., 0.]))),
+        )),
+    ];
+    objects.extend(lights);
+    info!("objects {:?}", objects);
+    let world: Box<dyn Hit> = if opt.use_accel {
+        Box::new(KDTree::new(objects, time_min, time_max))
+    } else {
+        Box::new(HitableList::new(objects))
+    };
+    Scene {
+        camera,
+        world,
+        subsamples: opt.subsamples,
+        num_threads: opt.num_threads,
+        width: opt.width,
+        height: opt.height,
+        ..Default::default()
+    }
+}

rtiow/renderer/src/sphere.rs

@@ -8,6 +8,7 @@ use crate::{
     vec3::{dot, Vec3},
 };
 
+#[derive(Debug)]
 pub struct Sphere<M>
 where
     M: Material,

rtiow/renderer/src/texture/checker.rs

@@ -1,5 +1,6 @@
 use crate::{texture::Texture, vec3::Vec3};
 
+#[derive(Debug)]
 pub struct CheckerTexture<T>
 where
     T: Texture,

rtiow/renderer/src/texture/constant.rs

@@ -1,6 +1,6 @@
 use crate::{texture::Texture, vec3::Vec3};
 
-#[derive(Debug, PartialEq)]
+#[derive(Clone, Debug, PartialEq)]
 pub struct ConstantTexture {
     color: Vec3,
 }

rtiow/renderer/src/texture/mandelbrot.rs

@@ -1,47 +1,12 @@
 #![allow(clippy::many_single_char_names)]
-use rand::{self, Rng};
 
-use crate::{texture::Texture, vec3::Vec3};
+use crate::{colors::generate_palette, texture::Texture, vec3::Vec3};
 
+#[derive(Debug)]
 pub struct Mandelbrot {
     palette: Vec<Vec3>,
 }
 
-// HSV values in [0..1]
-// returns [r, g, b] values from 0 to 255
-//From https://martin.ankerl.com/2009/12/09/how-to-create-random-colors-programmatically/
-fn hsv_to_rgb(h: f32, s: f32, v: f32) -> Vec3 {
-    let h_i = (h * 6.) as i32;
-    let f = h * 6. - h_i as f32;
-    let p = v * (1. - s);
-    let q = v * (1. - f * s);
-    let t = v * (1. - (1. - f) * s);
-    match h_i {
-        0 => Vec3::new(v, t, p),
-        1 => Vec3::new(q, v, p),
-        2 => Vec3::new(p, v, t),
-        3 => Vec3::new(p, q, v),
-        4 => Vec3::new(t, p, v),
-        5 => Vec3::new(v, p, q),
-        _ => panic!("Unknown H value {}", h_i),
-    }
-}
-
-fn generate_palette(num: usize) -> Vec<Vec3> {
-    let mut rng = rand::thread_rng();
-    let mut random = || rng.gen();
-    // use golden ratio
-    let golden_ratio_conjugate = 0.618_034;
-    let mut h = random();
-    (0..num)
-        .map(|_| {
-            h += golden_ratio_conjugate;
-            h %= 1.0;
-            hsv_to_rgb(h, 0.99, 0.99)
-        })
-        .collect()
-}
-
 impl Default for Mandelbrot {
     fn default() -> Self {
         Mandelbrot {

rtiow/renderer/src/texture/mod.rs

@@ -9,11 +9,11 @@ pub use crate::texture::{
     mandelbrot::Mandelbrot, noise::NoiseTexture,
 };
 
-use std::sync::Arc;
+use std::{fmt::Debug, sync::Arc};
 
 use crate::vec3::Vec3;
 
-pub trait Texture: Send + Sync {
+pub trait Texture: Send + Sync + Debug {
     fn value(&self, u: f32, v: f32, p: Vec3) -> Vec3;
 }
 
@@ -29,6 +29,12 @@ impl Texture for Box<dyn Texture> {
     }
 }
 
+impl Texture for [f32; 3] {
+    fn value(&self, u: f32, v: f32, p: Vec3) -> Vec3 {
+        (*self).into()
+    }
+}
+
 #[cfg(test)]
 mod test {
     use super::*;

rtiow/renderer/src/texture/noise.rs

@@ -4,6 +4,7 @@ use crate::{
     vec3::Vec3,
 };
 
+#[derive(Debug)]
 pub struct NoiseTexture<N>
 where
     N: NoiseSource,

rtiow/renderer/src/translate.rs

@@ -5,6 +5,7 @@ use crate::{
     vec3::Vec3,
 };
 
+#[derive(Debug)]
 pub struct Translate<H>
 where
     H: Hit,

rtiow/renderer/src/triangles.rs

@@ -0,0 +1,168 @@
+use std::f32::EPSILON;
+
+use stl::STL;
+
+use crate::{
+    aabb::AABB,
+    hitable::{Hit, HitRecord},
+    material::Material,
+    ray::Ray,
+    vec3::{cross, dot, Vec3},
+};
+
+#[derive(Debug)]
+pub struct Triangle {
+    normal: Vec3,
+    verts: [Vec3; 3],
+    // Precomputed data
+}
+
+#[derive(Debug)]
+pub struct Triangles<M>
+where
+    M: Material,
+{
+    pub triangles: Vec<Triangle>,
+    pub bbox: AABB,
+    material: M,
+}
+
+impl<M> Triangles<M>
+where
+    M: Material,
+{
+    pub fn new(stl: &STL, material: M, scale_factor: f32) -> Triangles<M> {
+        let triangles: Vec<_> = stl
+            .triangles
+            .iter()
+            .map(|t| {
+                let v0 = t.verts[0] * scale_factor;
+                let v1 = t.verts[1] * scale_factor;
+                let v2 = t.verts[2] * scale_factor;
+                Triangle {
+                    normal: t.normal,
+                    verts: [v0, v1, v2],
+                }
+            })
+            .collect();
+        let (min, max) = triangles.iter().fold(
+            (Vec3::from(f32::MAX), Vec3::from(f32::MIN)),
+            |(min, max), t| {
+                let t_min_x = t.verts[0].x.min(t.verts[1].x.min(t.verts[2].x));
+                let t_min_y = t.verts[0].y.min(t.verts[1].y.min(t.verts[2].y));
+                let t_min_z = t.verts[0].z.min(t.verts[1].z.min(t.verts[2].z));
+                let t_max_x = t.verts[0].x.max(t.verts[1].x.max(t.verts[2].x));
+                let t_max_y = t.verts[0].y.max(t.verts[1].y.max(t.verts[2].y));
+                let t_max_z = t.verts[0].z.max(t.verts[1].z.max(t.verts[2].z));
+
+                (
+                    Vec3::from([min.x.min(t_min_x), min.y.min(t_min_y), min.z.min(t_min_z)]),
+                    Vec3::from([max.x.max(t_max_x), max.y.max(t_max_y), max.z.max(t_max_z)]),
+                )
+            },
+        );
+        let bbox = AABB::new(min, max);
+        Triangles {
+            triangles,
+            bbox,
+            material,
+        }
+    }
+}
+
+/// 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> {
+    // #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
+}
+
+impl<M> Hit for Triangles<M>
+where
+    M: Material,
+{
+    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 {
+            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) {
+                // We don't support UV (yet?).
+                let uv = (0.5, 0.5);
+                return Some(HitRecord {
+                    t,
+                    uv,
+                    p,
+                    normal: tri.normal,
+                    material: &self.material,
+                });
+            }
+        }
+        None
+    }
+
+    fn bounding_box(&self, _t_min: f32, _t_max: f32) -> Option<AABB> {
+        Some(self.bbox)
+    }
+}
+
+struct RayTriangleResult {
+    t: f32,
+    p: Vec3,
+}

rtiow/renderer/src/vec3.rs

@@ -1,274 +1 @@
-use std::{
-    convert::From,
-    fmt,
-    num::ParseFloatError,
-    ops::{Add, Div, Index, Mul, Neg, Sub},
-    str,
-};
-
-use serde_derive::{Deserialize, Serialize};
-
-#[derive(Default, Debug, Deserialize, Serialize, PartialEq, Copy, Clone)]
-pub struct Vec3 {
-    pub x: f32,
-    pub y: f32,
-    pub z: f32,
-}
-
-pub fn cross(v1: Vec3, v2: Vec3) -> Vec3 {
-    Vec3 {
-        x: v1.y * v2.z - v1.z * v2.y,
-        y: v1.z * v2.x - v1.x * v2.z,
-        z: v1.x * v2.y - v1.y * v2.x,
-    }
-}
-
-pub fn dot(v1: Vec3, v2: Vec3) -> f32 {
-    v1.x * v2.x + v1.y * v2.y + v1.z * v2.z
-}
-
-impl Vec3 {
-    pub fn new(x: f32, y: f32, z: f32) -> Vec3 {
-        Vec3 { x, y, z }
-    }
-
-    pub fn min(self) -> f32 {
-        self.x.min(self.y).min(self.z)
-    }
-
-    pub fn max(self) -> f32 {
-        self.x.max(self.y).max(self.z)
-    }
-
-    pub fn length(self) -> f32 {
-        self.squared_length().sqrt()
-    }
-
-    pub fn squared_length(self) -> f32 {
-        self.x * self.x + self.y * self.y + self.z * self.z
-    }
-
-    pub fn unit_vector(self) -> Vec3 {
-        self / self.length()
-    }
-
-    pub fn make_unit_vector(&mut self) {
-        *self = self.unit_vector();
-    }
-}
-
-impl From<[f32; 3]> for Vec3 {
-    fn from(v: [f32; 3]) -> Self {
-        Vec3 {
-            x: v[0],
-            y: v[1],
-            z: v[2],
-        }
-    }
-}
-
-impl fmt::Display for Vec3 {
-    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        write!(f, "{} {} {}", self.x, self.y, self.z)
-    }
-}
-
-impl str::FromStr for Vec3 {
-    type Err = ParseFloatError;
-    fn from_str(s: &str) -> Result<Self, Self::Err> {
-        let coords: Vec<&str> = s.split(' ').collect();
-        Ok(Vec3 {
-            x: coords[0].parse::<f32>()?,
-            y: coords[1].parse::<f32>()?,
-            z: coords[2].parse::<f32>()?,
-        })
-    }
-}
-
-impl Add<f32> for Vec3 {
-    type Output = Vec3;
-
-    fn add(self, r: f32) -> Vec3 {
-        Vec3 {
-            x: self.x + r,
-            y: self.y + r,
-            z: self.z + r,
-        }
-    }
-}
-
-impl Add for Vec3 {
-    type Output = Vec3;
-
-    fn add(self, r: Vec3) -> Vec3 {
-        Vec3 {
-            x: self.x + r.x,
-            y: self.y + r.y,
-            z: self.z + r.z,
-        }
-    }
-}
-
-impl Div<Vec3> for f32 {
-    type Output = Vec3;
-
-    fn div(self, r: Vec3) -> Vec3 {
-        Vec3 {
-            x: self / r.x,
-            y: self / r.y,
-            z: self / r.z,
-        }
-    }
-}
-
-impl Div<f32> for Vec3 {
-    type Output = Vec3;
-
-    fn div(self, r: f32) -> Vec3 {
-        Vec3 {
-            x: self.x / r,
-            y: self.y / r,
-            z: self.z / r,
-        }
-    }
-}
-
-impl Index<usize> for Vec3 {
-    type Output = f32;
-    fn index(&self, idx: usize) -> &f32 {
-        match idx {
-            0 => &self.x,
-            1 => &self.y,
-            2 => &self.z,
-            _ => panic!("idx {} out of range for vec3", idx),
-        }
-    }
-}
-
-impl Mul for Vec3 {
-    type Output = Vec3;
-
-    fn mul(self, r: Vec3) -> Vec3 {
-        Vec3 {
-            x: self.x * r.x,
-            y: self.y * r.y,
-            z: self.z * r.z,
-        }
-    }
-}
-
-impl Mul<Vec3> for f32 {
-    type Output = Vec3;
-
-    fn mul(self, v: Vec3) -> Vec3 {
-        Vec3 {
-            x: v.x * self,
-            y: v.y * self,
-            z: v.z * self,
-        }
-    }
-}
-
-impl Mul<f32> for Vec3 {
-    type Output = Vec3;
-
-    fn mul(self, r: f32) -> Vec3 {
-        Vec3 {
-            x: self.x * r,
-            y: self.y * r,
-            z: self.z * r,
-        }
-    }
-}
-
-impl Neg for Vec3 {
-    type Output = Vec3;
-
-    fn neg(self) -> Vec3 {
-        -1. * self
-    }
-}
-
-impl Sub for Vec3 {
-    type Output = Vec3;
-
-    fn sub(self, r: Vec3) -> Vec3 {
-        Vec3 {
-            x: self.x - r.x,
-            y: self.y - r.y,
-            z: self.z - r.z,
-        }
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use std::{f32::consts::PI, str::FromStr};
-
-    use super::*;
-
-    #[test]
-    fn vec_add() {
-        let v0 = Vec3::new(1., 2., 3.);
-        let v1 = Vec3::new(1., 1., 1.);
-        assert_eq!(v0 + v1, Vec3::new(2., 3., 4.));
-    }
-
-    #[test]
-    fn vec_div() {
-        let v0 = Vec3::new(1., 2., 4.);
-        assert_eq!(v0 / 2., Vec3::new(0.5, 1., 2.));
-    }
-
-    #[test]
-    fn vec_mul() {
-        let v0 = Vec3::new(1., 2., 4.);
-        assert_eq!(v0 * 0.5, Vec3::new(0.5, 1., 2.));
-        assert_eq!(v0 * v0, Vec3::new(1., 4., 16.));
-    }
-
-    #[test]
-    fn vec_sub() {
-        let v0 = Vec3::new(1., 2., 3.);
-        let v1 = Vec3::new(1., 1., 1.);
-        assert_eq!(v0 - v1, Vec3::new(0., 1., 2.));
-    }
-
-    #[test]
-    fn vec_idx() {
-        let v0 = Vec3::new(1., 2., 3.);
-        assert_eq!(v0[2], 3.);
-    }
-
-    #[test]
-    fn vec_display() {
-        let v0 = Vec3::new(1., 2., 3.);
-        assert_eq!(format!("{}", v0), "1 2 3".to_owned());
-    }
-
-    #[test]
-    fn vec_from_str() {
-        assert_eq!(Vec3::from_str("1. 2. 3.").unwrap(), Vec3::new(1., 2., 3.));
-    }
-
-    #[test]
-    fn vec_str_roundtrip() {
-        let v = Vec3::from_str("1 2 3").unwrap();
-        let s = format!("{}", v);
-        assert_eq!(v, Vec3::from_str(&s).unwrap());
-    }
-
-    #[test]
-    fn vec_dot() {
-        let v0 = Vec3::new(1., 0., 0.);
-        let v1 = Vec3::new(-1., 0., 0.);
-        assert_eq!(dot(v0, v1), v0.length() * v1.length() * PI.cos());
-    }
-
-    #[test]
-    fn vec_cross() {
-        let v0 = Vec3::new(1., 0., 0.);
-        let v1 = Vec3::new(0., 1., 0.);
-        assert_eq!(cross(v0, v1), Vec3::new(0., 0., 1.));
-    }
-}
+pub use vec3::*;

rtiow/renderer/stls/dragon.stl

@@ -0,0 +1 @@
+stanford_dragon.stl

rtiow/renderer/stls/stanford_dragon-lowres.stl

@@ -0,0 +1 @@
+/home/wathiede/3dprint/stl/stanford_dragon-lowres.stl

rtiow/rustfmt.toml

@@ -1,2 +0,0 @@
-imports_granularity = "Crate"
-format_code_in_doc_comments = true

rtiow/tracer/Cargo.toml

@@ -2,12 +2,16 @@
 name = "tracer"
 version = "0.1.0"
 authors = ["Bill Thiede <git@xinu.tv>"]
-edition = "2018"
+edition = "2021"
+default-run = "tracer"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-log = "0.4.5"
+anyhow = "1.0.69"
+log = "0.4.17"
 renderer = { path = "../renderer" }
-stderrlog = "0.4.1"
-structopt = "0.2.10"
+stderrlog = "0.4.3"
+structopt = "0.2.18"
+strum = "0.24.1"
+toml = "0.7.2"

rtiow/tracer/configs/test.toml

@@ -0,0 +1,74 @@
+[scene]
+width = 768
+height = 512
+subsamples = 100
+
+[camera]
+lookfrom = [0.0, 50.0, 100.0]
+lookat = [0.0, 10.0, 0.0]
+fov = 45
+aperture = 0.0
+focus_dist = 10.0
+time_min = 0.0
+time_max = 1.0
+
+[[materials]]
+name = "light1"
+type = "isotropic"
+texture = [20, 10, 10]
+[[materials]]
+name = "yellow"
+type = "isotropic"
+texture = [1, 1, 0]
+[[materials]]
+name = "magenta"
+type = "lambertian"
+texture = [1, 0, 1]
+[[materials]]
+name = "green"
+type = "diffuse_light"
+texture = [0, 1, 0]
+[[materials]]
+name = "metal"
+type = "metal"
+albedo = [1, 1, 1]
+fuzzy = 0
+[[materials]]
+name = "glass"
+type = "dielectric"
+ref_idx = 1.5
+
+
+[[hitables]]
+type = "sphere"
+center = [-30.0, 0.0, 0.0]
+radius = 10
+material_name = "yellow"
+[[hitables]]
+type = "sphere"
+center = [30.0, 0.0, 0.0]
+radius = 10
+material_name = "green"
+[[hitables]]
+type = "sphere"
+center = [0.0, -10.0, 0.0]
+radius = 10
+material_name = "metal"
+[[hitables]]
+type = "sphere"
+center = [0.0, 0.0, -30.0]
+radius = 10
+material_name = "magenta"
+[[hitables]]
+type = "stl"
+path = "/net/nasx.h.xinu.tv/x/3dprint/stl/stanford_dragon.stl"
+scale = 200
+material_name = "glass"
+#[[hitables]]
+#type = "sphere"
+#center = [0.0, 50.0, -100.0]
+#radius = 10
+#material_name = "light1"
+
+[envmap]
+path = "/home/wathiede/src/xinu.tv/raytracers/rtiow/renderer/images/52681723945_e1d94d3df9_6k.jpg"

rtiow/tracer/src/main.rs

@@ -1,12 +1,17 @@
 #![warn(unused_extern_crates)]
-use std::fs;
+use std::{fs, time::Instant};
 
+use anyhow::Result;
 #[cfg(feature = "profile")]
 use cpuprofiler::PROFILER;
 use log::info;
 use structopt::StructOpt;
 
-use renderer::renderer::{render, Opt};
+use renderer::{
+    parser::Config,
+    renderer::{render, Model, Opt},
+};
+use strum::VariantNames;
 
 #[cfg(not(feature = "profile"))]
 struct MockTimer;
@@ -34,15 +39,36 @@ impl MockProfiler {
 #[cfg(not(feature = "profile"))]
 static PROFILER: MockProfiler = MockProfiler {};
 
-fn main() -> Result<(), std::io::Error> {
+fn main() -> Result<()> {
+    let start_time = Instant::now();
     stderrlog::new()
         .verbosity(3)
         .timestamp(stderrlog::Timestamp::Millisecond)
         .init()
         .unwrap();
     let opt = Opt::from_args();
-    info!("{:?}", opt);
-    let scene = opt.model.scene(&opt);
+    if opt.model.is_none() && opt.config.is_none() {
+        eprintln!(
+            "--config <path> or --model should be one of {:?}",
+            Model::VARIANTS
+        );
+        return Ok(());
+    }
+    if opt.model.is_some() && opt.config.is_some() {
+        eprintln!("only specify one of --config or --model");
+        return Ok(());
+    }
+    info!("{:#?}", opt);
+    let scene = match (&opt.model, &opt.config) {
+        (Some(model), None) => model.scene(&opt),
+        (None, Some(config)) => {
+            let s = std::fs::read_to_string(config)?;
+            let cfg: Config = toml::from_str(&s)?;
+            println!("{:#?}", cfg);
+            cfg.try_into()?
+        }
+        _ => unreachable!(),
+    };
     fs::create_dir_all(&opt.output)?;
     if opt.pprof.is_some() && !cfg!(feature = "profile") {
         panic!("profiling disabled at compile time, but -pprof specified");
@@ -54,11 +80,13 @@ fn main() -> Result<(), std::io::Error> {
             .start(pprof_path.to_str().unwrap().as_bytes())
             .unwrap();
     }
-    let res = render(scene, &opt.output);
+    let res = render(scene, &opt.output, &opt.tev_addr);
     if let Some(pprof_path) = &opt.pprof {
         info!("Saving pprof to {}", pprof_path.to_string_lossy());
         PROFILER.lock().unwrap().stop().unwrap();
     }
 
-    res
+    let time_diff = Instant::now() - start_time;
+    info!("Total runtime {} seconds", time_diff.as_secs_f32());
+    Ok(res?)
 }

rtiow/vec3/Cargo.toml

@@ -0,0 +1,10 @@
+[package]
+name = "vec3"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+serde = { version = "1.0.152", features = ["serde_derive"] }
+serde_derive = "1.0.152"

rtiow/vec3/src/lib.rs

@@ -0,0 +1,2 @@
+mod vec3;
+pub use vec3::*;

rtiow/vec3/src/vec3.rs

@@ -0,0 +1,330 @@
+use std::{
+    convert::From,
+    fmt,
+    num::ParseFloatError,
+    ops::{Add, Div, Index, Mul, Neg, Sub},
+    str,
+};
+
+use serde::{Deserialize, Serialize};
+
+#[derive(Default, Debug, Deserialize, Serialize, PartialEq, Copy, Clone)]
+pub struct Vec3 {
+    pub x: f32,
+    pub y: f32,
+    pub z: f32,
+}
+
+// Return a `Vec3` with the lowest of each component in v1 or v2.
+#[inline]
+pub fn min(v1: Vec3, v2: Vec3) -> Vec3 {
+    Vec3 {
+        x: v1.x.min(v2.x),
+        y: v1.y.min(v2.y),
+        z: v1.z.min(v2.z),
+    }
+}
+// Return a `Vec3` with the greatest of each component in v1 or v2.
+#[inline]
+pub fn max(v1: Vec3, v2: Vec3) -> Vec3 {
+    Vec3 {
+        x: v1.x.max(v2.x),
+        y: v1.y.max(v2.y),
+        z: v1.z.max(v2.z),
+    }
+}
+
+#[inline]
+pub fn cross(v1: Vec3, v2: Vec3) -> Vec3 {
+    Vec3 {
+        x: v1.y * v2.z - v1.z * v2.y,
+        y: v1.z * v2.x - v1.x * v2.z,
+        z: v1.x * v2.y - v1.y * v2.x,
+    }
+}
+
+#[inline]
+pub fn dot(v1: Vec3, v2: Vec3) -> f32 {
+    v1.x * v2.x + v1.y * v2.y + v1.z * v2.z
+}
+
+impl Vec3 {
+    pub fn new(x: f32, y: f32, z: f32) -> Vec3 {
+        Vec3 { x, y, z }
+    }
+
+    #[inline]
+    pub fn min(self) -> f32 {
+        self.x.min(self.y).min(self.z)
+    }
+
+    #[inline]
+    pub fn max(self) -> f32 {
+        self.x.max(self.y).max(self.z)
+    }
+
+    #[inline]
+    pub fn length(self) -> f32 {
+        self.squared_length().sqrt()
+    }
+
+    #[inline]
+    pub fn squared_length(self) -> f32 {
+        self.x * self.x + self.y * self.y + self.z * self.z
+    }
+
+    #[inline]
+    pub fn unit_vector(self) -> Vec3 {
+        self / self.length()
+    }
+
+    #[inline]
+    pub fn make_unit_vector(&mut self) {
+        *self = self.unit_vector();
+    }
+}
+
+impl From<f32> for Vec3 {
+    fn from(v: f32) -> Self {
+        Vec3 { x: v, y: v, z: v }
+    }
+}
+
+impl From<[f32; 3]> for Vec3 {
+    fn from(v: [f32; 3]) -> Self {
+        Vec3 {
+            x: v[0],
+            y: v[1],
+            z: v[2],
+        }
+    }
+}
+
+impl fmt::Display for Vec3 {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(f, "{:>6.2} {:>6.2} {:>6.2}", self.x, self.y, self.z)
+    }
+}
+
+impl str::FromStr for Vec3 {
+    type Err = ParseFloatError;
+    fn from_str(s: &str) -> Result<Self, Self::Err> {
+        let coords: Vec<&str> = s.split(' ').collect();
+        Ok(Vec3 {
+            x: coords[0].parse::<f32>()?,
+            y: coords[1].parse::<f32>()?,
+            z: coords[2].parse::<f32>()?,
+        })
+    }
+}
+
+impl Add<f32> for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn add(self, r: f32) -> Vec3 {
+        Vec3 {
+            x: self.x + r,
+            y: self.y + r,
+            z: self.z + r,
+        }
+    }
+}
+
+impl Add for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn add(self, r: Vec3) -> Vec3 {
+        Vec3 {
+            x: self.x + r.x,
+            y: self.y + r.y,
+            z: self.z + r.z,
+        }
+    }
+}
+
+impl Div<Vec3> for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn div(self, r: Vec3) -> Vec3 {
+        Vec3 {
+            x: self.x / r.x,
+            y: self.y / r.y,
+            z: self.z / r.z,
+        }
+    }
+}
+
+impl Div<Vec3> for f32 {
+    type Output = Vec3;
+
+    #[inline]
+    fn div(self, r: Vec3) -> Vec3 {
+        Vec3 {
+            x: self / r.x,
+            y: self / r.y,
+            z: self / r.z,
+        }
+    }
+}
+
+impl Div<f32> for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn div(self, r: f32) -> Vec3 {
+        Vec3 {
+            x: self.x / r,
+            y: self.y / r,
+            z: self.z / r,
+        }
+    }
+}
+
+impl Index<usize> for Vec3 {
+    type Output = f32;
+    #[inline]
+    fn index(&self, idx: usize) -> &f32 {
+        match idx {
+            0 => &self.x,
+            1 => &self.y,
+            2 => &self.z,
+            _ => panic!("idx {} out of range for vec3", idx),
+        }
+    }
+}
+
+impl Mul for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn mul(self, r: Vec3) -> Vec3 {
+        Vec3 {
+            x: self.x * r.x,
+            y: self.y * r.y,
+            z: self.z * r.z,
+        }
+    }
+}
+
+impl Mul<Vec3> for f32 {
+    type Output = Vec3;
+
+    #[inline]
+    fn mul(self, v: Vec3) -> Vec3 {
+        Vec3 {
+            x: v.x * self,
+            y: v.y * self,
+            z: v.z * self,
+        }
+    }
+}
+
+impl Mul<f32> for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn mul(self, r: f32) -> Vec3 {
+        Vec3 {
+            x: self.x * r,
+            y: self.y * r,
+            z: self.z * r,
+        }
+    }
+}
+
+impl Neg for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn neg(self) -> Vec3 {
+        -1. * self
+    }
+}
+
+impl Sub for Vec3 {
+    type Output = Vec3;
+
+    #[inline]
+    fn sub(self, r: Vec3) -> Vec3 {
+        Vec3 {
+            x: self.x - r.x,
+            y: self.y - r.y,
+            z: self.z - r.z,
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use std::{f32::consts::PI, str::FromStr};
+
+    use super::*;
+
+    #[test]
+    fn vec_add() {
+        let v0 = Vec3::new(1., 2., 3.);
+        let v1 = Vec3::new(1., 1., 1.);
+        assert_eq!(v0 + v1, Vec3::new(2., 3., 4.));
+    }
+
+    #[test]
+    fn vec_div() {
+        let v0 = Vec3::new(1., 2., 4.);
+        assert_eq!(v0 / 2., Vec3::new(0.5, 1., 2.));
+    }
+
+    #[test]
+    fn vec_mul() {
+        let v0 = Vec3::new(1., 2., 4.);
+        assert_eq!(v0 * 0.5, Vec3::new(0.5, 1., 2.));
+        assert_eq!(v0 * v0, Vec3::new(1., 4., 16.));
+    }
+
+    #[test]
+    fn vec_sub() {
+        let v0 = Vec3::new(1., 2., 3.);
+        let v1 = Vec3::new(1., 1., 1.);
+        assert_eq!(v0 - v1, Vec3::new(0., 1., 2.));
+    }
+
+    #[test]
+    fn vec_idx() {
+        let v0 = Vec3::new(1., 2., 3.);
+        assert_eq!(v0[2], 3.);
+    }
+
+    #[test]
+    fn vec_display() {
+        let v0 = Vec3::new(1., 2., 3.);
+        assert_eq!(format!("{}", v0), "1 2 3".to_owned());
+    }
+
+    #[test]
+    fn vec_from_str() {
+        assert_eq!(Vec3::from_str("1. 2. 3.").unwrap(), Vec3::new(1., 2., 3.));
+    }
+
+    #[test]
+    fn vec_str_roundtrip() {
+        let v = Vec3::from_str("1 2 3").unwrap();
+        let s = format!("{}", v);
+        assert_eq!(v, Vec3::from_str(&s).unwrap());
+    }
+
+    #[test]
+    fn vec_dot() {
+        let v0 = Vec3::new(1., 0., 0.);
+        let v1 = Vec3::new(-1., 0., 0.);
+        assert_eq!(dot(v0, v1), v0.length() * v1.length() * PI.cos());
+    }
+
+    #[test]
+    fn vec_cross() {
+        let v0 = Vec3::new(1., 0., 0.);
+        let v1 = Vec3::new(0., 1., 0.);
+        assert_eq!(cross(v0, v1), Vec3::new(0., 0., 1.));
+    }
+}

zigrtiow/dev.sh

@@ -0,0 +1,2 @@
+fd zig | entr -c bash -c "zig build test && zig build -Drelease-safe=true &&
+  ./zig-out/bin/zigrtiow | convert ppm:- /tmp/output.png"

zigrtiow/src/camera.zig

@@ -1,33 +1,70 @@
 const vec = @import("./vec.zig");
+const std = @import("std");
 
 const Ray = @import("./ray.zig").Ray;
 
 const Point3 = vec.Point3;
 const Vec3 = vec.Vec3;
+const random_in_unit_disk = vec.random_in_unit_disk;
+const math = std.math;
+const pi = math.pi;
+const tan = math.tan;
+
+// Utility Functions
+
+fn degrees_to_radians(degrees: f32) f32 {
+    return degrees * pi / 180.0;
+}
 
 pub const Camera = struct {
     origin: Point3,
     lower_left_corner: Point3,
     horizontal: Vec3,
     vertical: Vec3,
+    u: Vec3,
+    v: Vec3,
+    w: Vec3,
+    lens_radius: f32,
 
-    pub fn init() Camera {
-        const aspect_ratio = 16.0 / 9.0;
-        const viewport_height = 2.0;
+    pub fn init(
+        lookfrom: Point3,
+        lookat: Point3,
+        vup: Vec3,
+        vfov: f32, // vertical field-of-view in degrees
+        aspect_ratio: f32,
+        aperature: f32,
+        focus_dist: f32,
+    ) Camera {
+        const theta = degrees_to_radians(vfov);
+        const h = tan(theta / 2);
+        const viewport_height = 2 * h;
         const viewport_width = aspect_ratio * viewport_height;
-        const focal_length = 1.0;
-        const origin = Point3.init(0, 0, 0);
-        const horizontal = Vec3.init(viewport_width, 0.0, 0.0);
-        const vertical = Vec3.init(0.0, viewport_height, 0.0);
+
+        const w = lookfrom.sub(lookat).unit();
+        const u = vup.cross(w).unit();
+        const v = w.cross(u);
+
+        const origin = lookfrom;
+        const horizontal = u.scale(focus_dist * viewport_width);
+        const vertical = v.scale(focus_dist * viewport_height);
+        const lower_left_corner = origin.sub(horizontal.scale(0.5)).sub(vertical.scale(0.5)).sub(w.scale(focus_dist));
+
+        const lens_radius = aperature / 2;
 
         return Camera{
             .origin = origin,
             .horizontal = horizontal,
             .vertical = vertical,
-            .lower_left_corner = origin.sub(horizontal.scale(0.5)).sub(vertical.scale(0.5)).sub(Vec3.init(0, 0, focal_length)),
+            .lower_left_corner = lower_left_corner,
+            .u = u,
+            .v = v,
+            .w = w,
+            .lens_radius = lens_radius,
         };
     }
-    pub fn get_ray(camera: Camera, u: f32, v: f32) Ray {
-        return Ray.init(camera.origin, camera.lower_left_corner.add(camera.horizontal.scale(u).add(camera.vertical.scale(v).sub(camera.origin))));
+    pub fn get_ray(camera: Camera, s: f32, t: f32) Ray {
+        const rd = random_in_unit_disk().scale(camera.lens_radius);
+        const offset = camera.u.scale(rd.x()).add(camera.v.scale(rd.y()));
+        return Ray.init(camera.origin.add(offset), camera.lower_left_corner.add(camera.horizontal.scale(s).add(camera.vertical.scale(t).sub(camera.origin).sub(offset))));
     }
 };

zigrtiow/src/color.zig

@@ -8,11 +8,11 @@ pub fn write_color(c: Color, samples_per_pixel: isize) anyerror!void {
     var g = c.y();
     var b = c.z();
 
-    // Divide the color by the number of samples.
-    var scale = 1.0 / @intToFloat(f32, samples_per_pixel);
-    r *= scale;
-    g *= scale;
-    b *= scale;
+    // Divide the color by the number of samples and gamma-correct for gamma=2.0.
+    const scale = 1.0 / @intToFloat(f32, samples_per_pixel);
+    r = @sqrt(scale * r);
+    g = @sqrt(scale * g);
+    b = @sqrt(scale * b);
 
     // Write the translated [0,255] value of each color component.
     var ir = @floatToInt(u8, 256 * clamp(r, 0, 0.999));

zigrtiow/src/hittable.zig

@@ -1,6 +1,8 @@
 const vec = @import("./vec.zig");
 const ray = @import("./ray.zig");
 
+const Material = @import("./material.zig").Material;
+
 const Sphere = @import("./sphere.zig").Sphere;
 const HittableList = @import("./hittable_list.zig").HittableList;
 
@@ -25,6 +27,7 @@ pub const HitRecord = struct {
     normal: Vec3,
     t: f32,
     front_face: bool,
+    material: Material,
 
     pub fn set_face_normal(hr: *HitRecord, r: Ray, outward_normal: Vec3) void {
         hr.front_face = r.direction().dot(outward_normal) < 0;

zigrtiow/src/main.zig

@@ -6,72 +6,204 @@ const vec = @import("./vec.zig");
 const sphere = @import("./sphere.zig");
 const hittable = @import("./hittable.zig");
 const hittable_list = @import("./hittable_list.zig");
+const material = @import("./material.zig");
 
+const Thread = std.Thread;
+const Queue = std.atomic.Queue;
 const Camera = camera.Camera;
 const Color = vec.Color;
 const Hittable = hittable.Hittable;
 const HittableList = hittable_list.HittableList;
+const Labertian = material.Labertian;
+const Metal = material.Metal;
+const Dielectric = material.Dielectric;
+const Material = material.Material;
 const Point3 = vec.Point3;
 const Ray = ray.Ray;
 const Sphere = sphere.Sphere;
 const Vec3 = vec.Vec3;
-const random_in_unit_sphere = vec.random_in_unit_sphere;
 const info = std.log.info;
+const random_in_hemisphere = vec.random_in_hemisphere;
 const write_color = color.write_color;
+const math = std.math;
+const pi = math.pi;
+const cos = math.cos;
+const create_labertian = material.create_labertian;
+const create_metal = material.create_metal;
+const create_dielectric = material.create_dielectric;
 
 fn ray_color(r: Ray, world: Hittable, depth: isize) Color {
     // If we've exceeded the ray bounce limit, no more light is gathered.
     if (depth <= 0) return Color.init(0, 0, 0);
 
-    var hit = world.hit(r, 0, std.math.inf(f32));
-    if (hit) |rec| {
-        const target = rec.p.add(rec.normal.add(random_in_unit_sphere()));
-        return ray_color(Ray.init(rec.p, target.sub(rec.p)), world, depth - 1);
+    if (world.hit(r, 0.0001, std.math.inf(f32))) |hit_rec| {
+        if (hit_rec.material.scatter(r, hit_rec)) |scatter_rec| {
+            return scatter_rec.attenuation.mul(ray_color(scatter_rec.scattered, world, depth - 1));
+        }
+        return Color.init(0, 0, 0);
     }
     var unit_direction = r.direction().unit();
     const t = 0.5 * (unit_direction.y() + 1);
     return Color.init(1, 1, 1).scale(1 - t).add(Color.init(0.5, 0.7, 1.0).scale(t));
 }
 
+const Task = struct {
+    cam: Camera,
+    samples_per_pixel: usize,
+    max_depth: isize,
+    row_idx: usize,
+    image_height: isize,
+    image_width: isize,
+    row_pixels: []Color,
+    world: Hittable,
+};
+
+fn render_row(q: *Queue(Task)) void {
+    var prng = std.rand.DefaultPrng.init(0);
+    const rand = prng.random();
+    while (true) {
+        if (q.get()) |node| {
+            const task = node.data;
+            info("Rendering row: {}", .{task.row_idx});
+            const j = task.row_idx;
+            var i: usize = 0;
+            while (i < task.image_width) : (i += 1) {
+                var pixel_color = Color.init(0, 0, 0);
+                var s: isize = 0;
+                while (s < task.samples_per_pixel) : (s += 1) {
+                    const u = (@intToFloat(f32, i) + rand.float(f32)) / @intToFloat(f32, task.image_width - 1);
+                    const v = (@intToFloat(f32, j) + rand.float(f32)) / @intToFloat(f32, task.image_height - 1);
+                    const r = task.cam.get_ray(u, v);
+                    pixel_color = pixel_color.add(ray_color(r, task.world, task.max_depth));
+                }
+                task.row_pixels[i] = pixel_color;
+            }
+        } else {
+            return;
+        }
+    }
+}
+fn random_scene() anyerror!Hittable {
+    var prng = std.rand.DefaultPrng.init(0);
+    const rand = prng.random();
+    var world = HittableList.init();
+    const ground_material = create_labertian(Color.init(0.5, 0.5, 0.5));
+    try world.add(Hittable{ .sphere = Sphere.init(Point3.init(0, -1000, 0), 1000, ground_material) });
+
+    var a: isize = -11;
+    while (a < 11) : (a += 1) {
+        var b: isize = -11;
+        while (b < 11) : (b += 1) {
+            const choose_mat = rand.float(f32);
+            const center = Point3.init(@intToFloat(f32, a) + 0.9 * rand.float(f32), 0.2, @intToFloat(f32, b) + 0.9 * rand.float(f32));
+            if (center.sub(Point3.init(4, 0.2, 0)).length() > 0.9) {
+                if (choose_mat < 0.8) {
+                    const albedo = Color.random_min_max(0, 1).mul(Color.random_min_max(0, 1));
+                    const sphere_material = create_labertian(albedo);
+
+                    try world.add(Hittable{ .sphere = Sphere.init(center, 0.2, sphere_material) });
+                } else if (choose_mat < 0.95) {
+                    const albedo = Color.init(
+                        0.5 + rand.float(f32) / 2.0,
+                        0.5 + rand.float(f32) / 2.0,
+                        0.5 + rand.float(f32) / 2.0,
+                    );
+                    const fuzz = rand.float(f32) / 2.0;
+                    const sphere_material = create_metal(albedo, fuzz);
+
+                    try world.add(Hittable{ .sphere = Sphere.init(center, 0.2, sphere_material) });
+                } else {
+                    const sphere_material = create_dielectric(1.5);
+
+                    try world.add(Hittable{ .sphere = Sphere.init(center, 0.2, sphere_material) });
+                }
+            }
+        }
+    }
+
+    const material1 = create_dielectric(1.5);
+    try world.add(Hittable{ .sphere = Sphere.init(Point3.init(0, 1, 0), 1.0, material1) });
+
+    const material2 = create_labertian(Color.init(0.4, 0.2, 0.1));
+    try world.add(Hittable{ .sphere = Sphere.init(Point3.init(-4, 1, 0), 1.0, material2) });
+
+    const material3 = create_metal(Color.init(0.7, 0.6, 0.5), 0.0);
+    try world.add(Hittable{ .sphere = Sphere.init(Point3.init(4, 1, 0), 1.0, material3) });
+
+    return Hittable{ .hittable_list = world };
+}
+
 pub fn main() anyerror!void {
     // Image
-    const aspect_ratio: f32 = 16.0 / 9.0;
-    const image_width = 400;
+    const aspect_ratio: f32 = 3.0 / 2.0;
+    const image_width = 1200;
     const image_height = @floatToInt(isize, @intToFloat(f32, image_width) / aspect_ratio);
-    const samples_per_pixel = 100;
+    const samples_per_pixel = 500;
     const max_depth = 50;
 
     // World
-    var tmp_world = HittableList.init();
-    try tmp_world.add(Hittable{ .sphere = Sphere.init(Point3.init(0, 0, -1), 0.5) });
-    try tmp_world.add(Hittable{ .sphere = Sphere.init(Point3.init(0, -100.5, -1), 100) });
-    const world = Hittable{ .hittable_list = tmp_world };
+    const world = try random_scene();
 
     // Camera
-    const cam = Camera.init();
+    const lookfrom = Point3.init(13, 2, 3);
+    const lookat = Point3.init(0, 0, 0);
+    const vup = Vec3.init(0, 1, 0);
+    const dist_to_focus = 10;
+    const aperature = 0.1;
+
+    const cam = Camera.init(lookfrom, lookat, vup, 20, aspect_ratio, aperature, dist_to_focus);
 
     // Render
+    const Node = Queue(Task).Node;
+    const allocator = std.heap.page_allocator;
+    const cpus = try Thread.getCpuCount();
+    var pixels = try allocator.create([image_width * image_height]Color);
+    var q = Queue(Task).init();
+    var threads: std.ArrayList(std.Thread) = std.ArrayList(std.Thread).init(allocator);
+    var row: usize = 0;
+    while (row < image_height) : (row += 1) {
+        const node = try allocator.create(Node);
+        const start: usize = row * image_width;
+        const end: usize = (row + 1) * image_width;
+        node.* = .{
+            .prev = undefined,
+            .next = undefined,
+            .data = Task{
+                .cam = cam,
+                .samples_per_pixel = samples_per_pixel,
+                .max_depth = max_depth,
+                .row_idx = row,
+                .image_width = image_width,
+                .image_height = image_height,
+                .row_pixels = pixels[start..end],
+                .world = world,
+            },
+        };
+        q.put(node);
+    }
+
+    var t: usize = 0;
+    while (t < cpus) : (t += 1) {
+        try threads.append(try Thread.spawn(.{}, render_row, .{&q}));
+    }
+
+    // Wait for all threads to finish.
+    for (threads.items) |thread| {
+        Thread.join(thread);
+    }
+
     const stdout = std.io.getStdOut();
 
     info("Writing image {d}x{d}", .{ image_width, image_height });
     try stdout.writer().print("P3\n{d} {d}\n255\n", .{ image_width, image_height });
 
-    var prng = std.rand.DefaultPrng.init(0);
-    const rand = prng.random();
-    var j: isize = image_height - 1;
-    while (j >= 0) : (j -= 1) {
-        info("Scanlines remaining: {d}", .{j});
-
-        var i: isize = 0;
+    var j: usize = 0;
+    while (j < image_height) : (j += 1) {
+        var i: usize = 0;
         while (i < image_width) : (i += 1) {
-            var pixel_color = Color.init(0, 0, 0);
-            var s: isize = 0;
-            while (s < samples_per_pixel) : (s += 1) {
-                const u = (@intToFloat(f32, i) + rand.float(f32)) / @intToFloat(f32, image_width - 1);
-                const v = (@intToFloat(f32, j) + rand.float(f32)) / @intToFloat(f32, image_height - 1);
-                const r = cam.get_ray(u, v);
-                pixel_color = pixel_color.add(ray_color(r, world, max_depth));
-            }
+            // Flip image upside down.
+            const j_idx = @as(usize, image_height) - j - 1;
+            const pixel_color = pixels[i + j_idx * image_width];
             try write_color(pixel_color, samples_per_pixel);
         }
     }

zigrtiow/src/material.zig

@@ -0,0 +1,116 @@
+const std = @import("std");
+const vec = @import("./vec.zig");
+const Ray = @import("./ray.zig").Ray;
+const HitRecord = @import("./hittable.zig").HitRecord;
+
+const Vec3 = vec.Vec3;
+const Color = vec.Color;
+const random_unit_vector = vec.random_unit_vector;
+const random_in_unit_sphere = vec.random_in_unit_sphere;
+const pow = std.math.pow;
+
+pub fn create_labertian(albedo: Vec3) Material {
+    return Material{ .labertian = Labertian{ .albedo = albedo } };
+}
+
+pub fn create_metal(albedo: Vec3, fuzz: f32) Material {
+    return Material{ .metal = Metal{ .albedo = albedo, .fuzz = fuzz } };
+}
+
+pub fn create_dielectric(ir: f32) Material {
+    return Material{ .dielectric = Dielectric{ .ir = ir } };
+}
+
+pub const Material = union(enum) {
+    labertian: Labertian,
+    metal: Metal,
+    dielectric: Dielectric,
+
+    pub fn scatter(material: Material, r_in: Ray, rec: HitRecord) ?ScatterRec {
+        return switch (material) {
+            .labertian => |labertian| labertian.scatter(r_in, rec),
+            .metal => |metal| metal.scatter(r_in, rec),
+            .dielectric => |dielectric| dielectric.scatter(r_in, rec),
+        };
+    }
+};
+
+pub const ScatterRec = struct { attenuation: Color, scattered: Ray };
+
+fn reflect(v: Vec3, n: Vec3) Vec3 {
+    return v.sub(n.scale(v.dot(n) * 2));
+}
+fn refract(uv: Vec3, n: Vec3, etai_over_etat: f32) Vec3 {
+    const cos_theta = @minimum(-uv.dot(n), 1.0);
+
+    const r_out_perp = uv.add(n.scale(cos_theta)).scale(etai_over_etat);
+    const r_out_parallel = n.scale(-@sqrt(@fabs(1.0 - r_out_perp.length_squared())));
+    return r_out_perp.add(r_out_parallel);
+}
+
+pub const Labertian = struct {
+    albedo: Vec3,
+    pub fn scatter(labertian: Labertian, r_in: Ray, rec: HitRecord) ?ScatterRec {
+        _ = r_in;
+        var scatter_direction = rec.normal.add(random_unit_vector());
+        if (scatter_direction.near_zero()) {
+            scatter_direction = rec.normal;
+        }
+        return ScatterRec{
+            .scattered = Ray.init(rec.p, scatter_direction),
+            .attenuation = labertian.albedo,
+        };
+    }
+};
+
+pub const Metal = struct {
+    albedo: Vec3,
+    fuzz: f32,
+
+    pub fn scatter(metal: Metal, r_in: Ray, rec: HitRecord) ?ScatterRec {
+        const reflected = reflect(r_in.direction().unit(), rec.normal);
+        const scattered = Ray.init(rec.p, reflected.add(random_in_unit_sphere().scale(metal.fuzz)));
+        const attenuation = metal.albedo;
+        if (scattered.direction().dot(rec.normal) > 0) {
+            return ScatterRec{
+                .scattered = scattered,
+                .attenuation = attenuation,
+            };
+        } else {
+            return null;
+        }
+    }
+};
+
+var prng = std.rand.DefaultPrng.init(0);
+const rand = prng.random();
+pub const Dielectric = struct {
+    ir: f32, // Index of Refraction
+
+    pub fn scatter(dielectric: Dielectric, r_in: Ray, rec: HitRecord) ?ScatterRec {
+        const attenuation = Color.init(1, 1, 1);
+        const refraction_ratio = if (rec.front_face) 1.0 / dielectric.ir else dielectric.ir;
+
+        const unit_direction = r_in.direction().unit();
+        const cos_theta = @minimum(-unit_direction.dot(rec.normal), 1.0);
+        const sin_theta = @sqrt(1 - cos_theta * cos_theta);
+
+        const cannot_refract = refraction_ratio * sin_theta > 1;
+        const direction = if (cannot_refract or reflectance(cos_theta, refraction_ratio) > rand.float(f32))
+            reflect(unit_direction, rec.normal)
+        else
+            refract(unit_direction, rec.normal, refraction_ratio);
+
+        const scattered = Ray.init(rec.p, direction);
+        return ScatterRec{
+            .attenuation = attenuation,
+            .scattered = scattered,
+        };
+    }
+    fn reflectance(cosine: f32, ref_idx: f32) f32 {
+        // Use Schlick's approximation for reflectance.
+        var r0 = (1 - ref_idx) / (1 + ref_idx);
+        r0 = r0 * r0;
+        return r0 + (1 - r0) * pow(f32, (1 - cosine), 5);
+    }
+};

zigrtiow/src/sphere.zig

@@ -2,6 +2,8 @@ const vec = @import("./vec.zig");
 const ray = @import("./ray.zig");
 const hittable = @import("./hittable.zig");
 
+const Material = @import("./material.zig").Material;
+
 const Vec3 = vec.Vec3;
 const Point3 = vec.Point3;
 const Ray = ray.Ray;
@@ -10,11 +12,13 @@ const HitRecord = hittable.HitRecord;
 pub const Sphere = struct {
     center: Point3,
     radius: f32,
+    material: Material,
 
-    pub fn init(center: Point3, radius: f32) Sphere {
+    pub fn init(center: Point3, radius: f32, material: Material) Sphere {
         return Sphere{
             .center = center,
             .radius = radius,
+            .material = material,
         };
     }
     pub fn hit(sphere: Sphere, r: Ray, t_min: f32, t_max: f32) ?HitRecord {
@@ -44,6 +48,7 @@ pub const Sphere = struct {
             .p = p,
             .normal = Vec3.init(0, 0, 0),
             .front_face = false,
+            .material = sphere.material,
         };
         hr.set_face_normal(r, outward_normal);
         return hr;

zigrtiow/src/vec.zig

@@ -34,6 +34,9 @@ pub const Vec3 = struct {
     pub fn sub(lhs: Vec3, rhs: Vec3) Vec3 {
         return Vec3{ .v = lhs.v - rhs.v };
     }
+    pub fn mul(lhs: Vec3, rhs: Vec3) Vec3 {
+        return Vec3{ .v = lhs.v * rhs.v };
+    }
     pub fn scale(vec: Vec3, t: f32) Vec3 {
         return Vec3{ .v = vec.v * @splat(3, t) };
     }
@@ -44,23 +47,50 @@ pub const Vec3 = struct {
         const t = u.v * v.v;
         return t[0] + t[1] + t[2];
     }
-    /// Return value in unit cube from -1, 1
-    pub fn random() Vec3 {
+    pub fn cross(u: Vec3, v: Vec3) Vec3 {
+        return Vec3.init(u.v[1] * v.v[2] - u.v[2] * v.v[1], u.v[2] * v.v[0] - u.v[0] * v.v[2], u.v[0] * v.v[1] - u.v[1] * v.v[0]);
+    }
+    /// Return value in unit cube from [min, max).
+    pub fn random_min_max(min: f32, max: f32) Vec3 {
         return Vec3.init(
-            rand.float(f32) * 2 - 0.5,
-            rand.float(f32) * 2 - 0.5,
-            rand.float(f32) * 2 - 0.5,
+            min + (max - min) * rand.float(f32),
+            min + (max - min) * rand.float(f32),
+            min + (max - min) * rand.float(f32),
         );
     }
+    pub fn near_zero(vec: Vec3) bool {
+        // Return true if the vector is close to zero in all dimensions.
+        const s = 1e-8;
+        const e = vec.v;
+        return (@fabs(e[0]) < s) and (@fabs(e[1]) < s) and (@fabs(e[2]) < s);
+    }
 };
 
+pub fn random_in_unit_disk() Vec3 {
+    while (true) {
+        const p = Vec3.init(rand.float(f32) * 2 - 1, rand.float(f32) * 2 - 1, 0);
+        if (p.length_squared() >= 1) continue;
+        return p;
+    }
+}
 pub fn random_in_unit_sphere() Vec3 {
     while (true) {
-        const p = Vec3.random();
+        const p = Vec3.random_min_max(-1, 1);
         if (p.length_squared() > 1) continue;
         return p;
     }
 }
+pub fn random_unit_vector() Vec3 {
+    return random_in_unit_sphere().unit();
+}
+pub fn random_in_hemisphere(normal: Vec3) Vec3 {
+    const in_unit_sphere = random_in_unit_sphere();
+    if (in_unit_sphere.dot(normal) > 0.0) { // In the same hemisphere as the normal
+        return in_unit_sphere;
+    } else {
+        return in_unit_sphere.scale(-1);
+    }
+}
 
 pub const Color = Vec3;
 pub const Point3 = Vec3;

zigrtiow/test/threads.zig

@@ -0,0 +1,56 @@
+const std = @import("std");
+const Thread = std.Thread;
+const info = std.log.info;
+const Queue = std.atomic.Queue;
+
+const Task = struct {
+    row_idx: usize,
+    row_pixels: []u8,
+};
+
+fn worker(q: *Queue(Task)) void {
+    while (true) {
+        if (q.get()) |node| {
+            const task = node.data;
+            info("starting thread: {}", .{task.row_idx});
+        } else {
+            return;
+        }
+    }
+}
+
+pub fn main() anyerror!void {
+    const allocator = std.heap.page_allocator;
+    const cpus = try Thread.getCpuCount();
+    const width = 512;
+    const height = 256;
+    var pixels: [width * height * 3]u8 = undefined;
+
+    var threads: std.ArrayList(std.Thread) = std.ArrayList(std.Thread).init(allocator);
+    var q = Queue(Task).init();
+    const Node = Queue(Task).Node;
+    var row: usize = 0;
+    while (row < height) : (row += 1) {
+        const node = try allocator.create(Node);
+        node.* = .{
+            .prev = undefined,
+            .next = undefined,
+            .data = Task{
+                .row_idx = row,
+                .row_pixels = pixels[row * width * 3 .. (row + 1) * width * 3],
+            },
+        };
+        q.put(node);
+    }
+
+    var t: usize = 0;
+    while (t < cpus) : (t += 1) {
+        try threads.append(try Thread.spawn(.{}, worker, .{&q}));
+    }
+
+    // Wait for all threads to finish.
+    for (threads.items) |thread| {
+        Thread.join(thread);
+    }
+    info("main done", .{});
+}