Generalized the BVH to... pretty much anything.

The BVH is now generic over any kind of data.  The building
function takes in a closure that can bound the given data type
in 3d space, and the rest just works.

src/bvh.rs

@@ -7,8 +7,10 @@ use triangle;
 use algorithm::partition;
 
 #[derive(Debug)]
-pub struct BVH {
+pub struct BVH<'a, T: 'a> {
     nodes: Vec<BVHNode>,
+    objects: &'a [T],
+    depth: usize,
 }
 
 #[derive(Debug)]
@@ -20,36 +22,51 @@ enum BVHNode {
 
     Leaf {
         bounds: BBox,
-        triangle: (Point, Point, Point),
+        object_index: usize,
     },
 }
 
-impl BVH {
+impl<'a, T> BVH<'a, T> {
-    pub fn from_triangles(triangles: &mut [(Point, Point, Point)]) -> BVH {
+    pub fn from_objects<F>(objects: &'a mut [T], bounder: F) -> BVH<'a, T>
-        let mut bvh = BVH { nodes: Vec::new() };
+        where F: Fn(&T) -> BBox
+    {
+        let mut bvh = BVH {
+            nodes: Vec::new(),
+            objects: &[],
+            depth: 0,
+        };
 
-        bvh.recursive_build(triangles);
+        bvh.recursive_build(0, 0, objects, &bounder);
+        bvh.objects = objects;
+
+        println!("BVH Depth: {}", bvh.depth);
 
         bvh
     }
 
-    // Recursively builds the BVH starting at the given node with the given
+
-    // first and last primitive indices (in bag).
+    fn recursive_build<F>(&mut self,
-    fn recursive_build(&mut self, triangles: &mut [(Point, Point, Point)]) -> usize {
+                          offset: usize,
+                          depth: usize,
+                          objects: &mut [T],
+                          bounder: &F)
+                          -> usize
+        where F: Fn(&T) -> BBox
+    {
         let me = self.nodes.len();
 
-        if triangles.len() == 1 {
+        if objects.len() == 0 {
+            return 0;
+        } else if objects.len() == 1 {
             // Leaf node
-            let tri = triangles[0];
-
             self.nodes.push(BVHNode::Leaf {
-                bounds: {
+                bounds: bounder(&objects[0]),
-                    let minimum = tri.0.min(tri.1.min(tri.2));
+                object_index: offset,
-                    let maximum = tri.0.max(tri.1.max(tri.2));
-                    BBox::from_points(minimum, maximum)
-                },
-                triangle: tri,
             });
+
+            if self.depth < depth {
+                self.depth = depth;
+            }
         } else {
             // Not a leaf node
             self.nodes.push(BVHNode::Internal {
@@ -58,18 +75,10 @@ impl BVH {
             });
 
             // Determine which axis to split on
-            fn tri_bounds(tri: (Point, Point, Point)) -> BBox {
-                let minimum = tri.0.min(tri.1.min(tri.2));
-                let maximum = tri.0.max(tri.1.max(tri.2));
-                BBox {
-                    min: minimum,
-                    max: maximum,
-                }
-            }
             let bounds = {
                 let mut bounds = BBox::new();
-                for tri in &triangles[..] {
+                for obj in objects.iter() {
-                    bounds = bounds | tri_bounds(*tri);
+                    bounds = bounds | bounder(obj);
                 }
                 bounds
             };
@@ -87,10 +96,10 @@ impl BVH {
             };
             let split_pos = (bounds.min[split_axis] + bounds.max[split_axis]) * 0.5;
 
-            // Partition triangles based on split
+            // Partition objects based on split
             let split_index = {
-                let mut split_i = partition(triangles, |tri| {
+                let mut split_i = partition(objects, |obj| {
-                    let tb = tri_bounds(*tri);
+                    let tb = bounder(obj);
                     let centroid = (tb.min[split_axis] + tb.max[split_axis]) * 0.5;
                     centroid < split_pos
                 });
@@ -102,8 +111,11 @@ impl BVH {
             };
 
             // Create child nodes
-            self.recursive_build(&mut triangles[..split_index]);
+            self.recursive_build(offset, depth + 1, &mut objects[..split_index], bounder);
-            let child2_index = self.recursive_build(&mut triangles[split_index..]);
+            let child2_index = self.recursive_build(offset + split_index,
+                                                    depth + 1,
+                                                    &mut objects[split_index..],
+                                                    bounder);
 
             // Set node
             self.nodes[me] = BVHNode::Internal {
@@ -117,11 +129,18 @@ impl BVH {
 }
 
 
-pub fn intersect_bvh(bvh: &BVH, ray: &Ray) -> bool {
+pub fn intersect_bvh(bvh: &BVH<(Point, Point, Point)>, ray: &mut Ray) -> Option<(f32, f32, f32)> {
-    let mut i_stack = [0; 64];
+    if bvh.nodes.len() == 0 {
-    let mut stack_ptr = 0;
+        return None;
+    }
 
-    loop {
+    let mut i_stack = [0; 65];
+    let mut stack_ptr: usize = 1;
+    let mut hit = false;
+    let mut u = 0.0;
+    let mut v = 0.0;
+
+    while stack_ptr > 0 {
         match bvh.nodes[i_stack[stack_ptr]] {
             BVHNode::Internal { bounds, second_child_index } => {
                 if bounds.intersect_ray(ray) {
@@ -129,25 +148,29 @@ pub fn intersect_bvh(bvh: &BVH, ray: &Ray) -> bool {
                     i_stack[stack_ptr + 1] = second_child_index;
                     stack_ptr += 1;
                 } else {
-                    if stack_ptr == 0 {
-                        break;
-                    }
                     stack_ptr -= 1;
                 }
             }
 
-            BVHNode::Leaf{bounds: _, triangle: tri} => {
+            BVHNode::Leaf { bounds: _, object_index } => {
-                if let Some(_) = triangle::intersect_ray(ray, tri) {
+                if let Some((t, tri_u, tri_v)) =
-                    return true;
+                       triangle::intersect_ray(ray, bvh.objects[object_index]) {
+                    if t < ray.max_t {
+                        hit = true;
+                        ray.max_t = t;
+                        u = tri_u;
+                        v = tri_v;
+                    }
+                }
+
+                stack_ptr -= 1;
+            }
+        }
+    }
+
+    if hit {
+        return Some((ray.max_t, u, v));
     } else {
-                    if stack_ptr == 0 {
+        return None;
-                        break;
     }
-                    stack_ptr -= 1;
-                }
-            }
-        }
-    }
-
-    return false;
 }

src/main.rs

@@ -20,6 +20,7 @@ use docopt::Docopt;
 use image::Image;
 use math::{Point, Vector, fast_logit};
 use ray::Ray;
+use bbox::BBox;
 
 // ----------------------------------------------------------------
 
@@ -74,18 +75,32 @@ fn main() {
     // Generate a scene of triangles
     let mut triangles = {
         let mut triangles = Vec::new();
-        for x in 0..10 {
+        let xres = 512;
-            for y in 0..10 {
+        let yres = 512;
-                let cx = x as f32 * 32.0;
+        let xinc = 512.0 / (xres as f32);
-                let cy = y as f32 * 32.0;
+        let yinc = 512.0 / (yres as f32);
+        for x in 0..xres {
+            for y in 0..yres {
+                let cx = x as f32 * xinc;
+                let cy = y as f32 * yinc;
                 triangles.push((Point::new(cx, cy, 1.0),
-                                Point::new(cx + 32.0, cy, 1.0),
+                                Point::new(cx + xinc, cy, 1.1),
-                                Point::new(cx, cy + 32.0, 1.0)));
+                                Point::new(cx, cy + yinc, 1.2)));
+                triangles.push((Point::new(cx + xinc, cy + yinc, 1.0),
+                                Point::new(cx, cy + yinc, 1.1),
+                                Point::new(cx + xinc, cy, 1.2)));
             }
         }
         triangles
     };
-    let scene = bvh::BVH::from_triangles(&mut triangles[..]);
+    let scene = bvh::BVH::from_objects(&mut triangles[..], |tri| {
+        let minimum = tri.0.min(tri.1.min(tri.2));
+        let maximum = tri.0.max(tri.1.max(tri.2));
+        BBox {
+            min: minimum,
+            max: maximum,
+        }
+    });
     println!("Scene built.");
 
     // Write output image of ray-traced triangle
@@ -96,23 +111,29 @@ fn main() {
             let mut g = 0.0;
             let mut b = 0.0;
             let offset = hash_u32(((x as u32) << 16) ^ (y as u32), 0);
-            const SAMPLES: usize = 64;
+            const SAMPLES: usize = 16;
             for si in 0..SAMPLES {
-                let ray = Ray::new(Point::new(x as f32 +
+                let mut ray = Ray::new(Point::new(x as f32 +
-                                              fast_logit(halton::sample(0, offset + si as u32),
+                                                  fast_logit(halton::sample(0,
+                                                                            offset + si as u32),
                                                              1.5),
                                                   y as f32 +
-                                              fast_logit(halton::sample(3, offset + si as u32),
+                                                  fast_logit(halton::sample(3,
+                                                                            offset + si as u32),
                                                              1.5),
                                                   0.0),
                                        Vector::new(0.0, 0.0, 1.0));
-                if bvh::intersect_bvh(&scene, &ray) {
+                if let Some((_, u, v)) = bvh::intersect_bvh(&scene, &mut ray) {
-                    r += 1.0;
+                    r += u;
-                    g += 1.0;
+                    g += v;
-                    b += 1.0;
+                    b += (1.0 - u - v).max(0.0);
-                    // r += u;
+                    // r += 1.0;
-                    // g += v;
+                    // g += 1.0;
-                    // b += (1.0 - u - v).max(0.0);
+                    // b += 1.0;
+                } else {
+                    r += 0.1;
+                    g += 0.1;
+                    b += 0.1;
                 }
             }
             r *= 255.0 / SAMPLES as f32;

src/math/mod.rs

@@ -40,14 +40,20 @@ pub fn fast_ln(x: f32) -> f32 {
     return y - 87.989971088;
 }
 
-pub fn fast_logit(p: f32, width: f32) -> f32 {
-    let n = 0.001 + (p * 0.998);
 
-    fast_ln((n / (1.0 - n))) * width * (0.6266 / 4.0)
-}
 
+/// The logit function, scaled to approximate the probit function.
+///
+/// We use this as a close approximation to the gaussian inverse CDF,
+/// since the gaussian inverse CDF (probit) has no analytic formula.
 pub fn logit(p: f32, width: f32) -> f32 {
     let n = 0.001 + (p * 0.998);
 
     (n / (1.0 - n)).ln() * width * (0.6266 / 4.0)
 }
+
+pub fn fast_logit(p: f32, width: f32) -> f32 {
+    let n = 0.001 + (p * 0.998);
+
+    fast_ln((n / (1.0 - n))) * width * (0.6266 / 4.0)
+}