Implemented light tree sampling, for better sampling of many lights.

src/assembly.rs

@@ -3,7 +3,7 @@ use std::collections::HashMap;
 use math::{Matrix4x4, Vector};
 use lerp::lerp_slice;
 use bvh::BVH;
-use light_accel::{LightAccel, LightArray};
+use light_accel::{LightAccel, LightTree};
 use boundable::Boundable;
 use surface::{Surface, SurfaceIntersection};
 use light::LightSource;
@@ -28,7 +28,7 @@ pub struct Assembly {
     pub object_accel: BVH,
 
     // Light accel
-    pub light_accel: LightArray,
+    pub light_accel: LightTree,
 }
 
 impl Assembly {
@@ -40,8 +40,9 @@ impl Assembly {
                          time: f32,
                          intr: &SurfaceIntersection)
                          -> Option<(SpectralSample, Vector, f32)> {
-        if let &SurfaceIntersection::Hit { pos, .. } = intr {
+        if let &SurfaceIntersection::Hit { pos, incoming, nor, closure, .. } = intr {
-            if let Some((light_i, sel_pdf, _)) = self.light_accel.select(n) {
+            if let Some((light_i, sel_pdf, _)) = self.light_accel
+                .select(incoming, pos, nor, closure.as_surface_closure(), time, n) {
                 let inst = self.light_instances[light_i];
                 match inst.instance_type {
                     InstanceType::Object => {
@@ -198,20 +199,42 @@ impl AssemblyBuilder {
                                              |inst| &bbs[bis[inst.id]..bis[inst.id + 1]]);
         println!("Assembly BVH Depth: {}", object_accel.tree_depth());
 
-        // Build light accel
+        // Get list of instances that are for light sources.
-        let mut light_instances = self.instances.clone();
+        // TODO: include assemblies that themselves contain light sources.
-        let light_accel = LightArray::new(&mut light_instances[..], |inst| {
+        let mut light_instances: Vec<_> = self.instances
+            .iter()
+            .filter(|inst| {
                 match inst.instance_type {
                     InstanceType::Object => {
                         if let Object::Light(_) = self.objects[inst.data_index] {
-                        Some((&bbs[bis[inst.id]..bis[inst.id + 1]], 1.0))
+                            true
                         } else {
-                        None
+                            false
                         }
                     }
 
-                _ => None,
+                    _ => false,
                 }
+            })
+            .map(|&a| a)
+            .collect();
+
+        // Build light accel
+        let light_accel = LightTree::from_objects(&mut light_instances[..], |inst| {
+            let bounds = &bbs[bis[inst.id]..bis[inst.id + 1]];
+            let energy = match inst.instance_type {
+                InstanceType::Object => {
+                    if let Object::Light(ref light) = self.objects[inst.data_index] {
+                        light.approximate_energy()
+                    } else {
+                        0.0
+                    }
+                }
+
+                // TODO: handle assemblies.
+                _ => 0.0,
+            };
+            (bounds, energy)
         });
 
         Assembly {

src/bbox.rs

@@ -85,6 +85,14 @@ impl BBox {
 
         ((x * y) + (y * z) + (z * x)) * 2.0
     }
+
+    pub fn center(&self) -> Point {
+        self.min.lerp(self.max, 0.5)
+    }
+
+    pub fn diagonal(&self) -> f32 {
+        (self.max - self.min).length()
+    }
 }
 
 

src/light/mod.rs

@@ -68,11 +68,16 @@ pub trait LightSource: Boundable + Debug + Sync {
                 -> SpectralSample;
 
 
-
     /// Returns whether the light has a delta distribution.
     ///
     /// If a light has no chance of a ray hitting it through random process
     /// then it is a delta light source.  For example, point light sources,
     /// lights that only emit in a single direction, etc.
     fn is_delta(&self) -> bool;
+
+
+    /// Returns an approximation of the total energy emitted by the light
+    /// source.  Note that this does not need to be exact: it is used for
+    /// importance sampling.
+    fn approximate_energy(&self) -> f32;
 }

src/light/rectangle_light.rs

@@ -156,6 +156,12 @@ impl LightSource for RectangleLight {
     fn is_delta(&self) -> bool {
         false
     }
+
+    fn approximate_energy(&self) -> f32 {
+        let color: XYZ = self.colors.iter().fold(XYZ::new(0.0, 0.0, 0.0), |a, &b| a + b) /
+                         self.colors.len() as f32;
+        color.y
+    }
 }
 
 impl Boundable for RectangleLight {

src/light/sphere_light.rs

@@ -155,6 +155,12 @@ impl LightSource for SphereLight {
     fn is_delta(&self) -> bool {
         false
     }
+
+    fn approximate_energy(&self) -> f32 {
+        let color: XYZ = self.colors.iter().fold(XYZ::new(0.0, 0.0, 0.0), |a, &b| a + b) /
+                         self.colors.len() as f32;
+        color.y
+    }
 }
 
 impl Boundable for SphereLight {

src/light_accel/light_tree.rs

@@ -0,0 +1,200 @@
+use bbox::BBox;
+use sah::sah_split;
+use lerp::lerp_slice;
+use algorithm::{partition, merge_slices_append};
+use math::{Vector, Point, Normal};
+use shading::surface_closure::SurfaceClosure;
+
+use super::LightAccel;
+
+#[derive(Debug)]
+pub struct LightTree {
+    nodes: Vec<Node>,
+    bounds: Vec<BBox>,
+    depth: usize,
+    bounds_cache: Vec<BBox>,
+}
+
+#[derive(Debug)]
+struct Node {
+    is_leaf: bool,
+    bounds_range: (usize, usize),
+    energy: f32,
+    child_index: usize,
+}
+
+impl LightTree {
+    pub fn from_objects<'a, T, F>(objects: &mut [T], info_getter: F) -> LightTree
+        where F: 'a + Fn(&T) -> (&'a [BBox], f32)
+    {
+        let mut tree = LightTree {
+            nodes: Vec::new(),
+            bounds: Vec::new(),
+            depth: 0,
+            bounds_cache: Vec::new(),
+        };
+
+        tree.recursive_build(0, 0, objects, &info_getter);
+        tree.bounds_cache.clear();
+        tree.bounds_cache.shrink_to_fit();
+
+        tree
+    }
+
+
+    fn recursive_build<'a, T, F>(&mut self,
+                                 offset: usize,
+                                 depth: usize,
+                                 objects: &mut [T],
+                                 info_getter: &F)
+                                 -> (usize, (usize, usize))
+        where F: 'a + Fn(&T) -> (&'a [BBox], f32)
+    {
+        let me_index = self.nodes.len();
+
+        if objects.len() == 0 {
+            return (0, (0, 0));
+        } else if objects.len() == 1 {
+            // Leaf node
+            let bi = self.bounds.len();
+            let (obj_bounds, energy) = info_getter(&objects[0]);
+            self.bounds.extend(obj_bounds);
+            self.nodes.push(Node {
+                is_leaf: true,
+                bounds_range: (bi, self.bounds.len()),
+                energy: energy,
+                child_index: offset,
+            });
+
+            if self.depth < depth {
+                self.depth = depth;
+            }
+
+            return (me_index, (bi, self.bounds.len()));
+        } else {
+            // Not a leaf node
+            self.nodes.push(Node {
+                is_leaf: false,
+                bounds_range: (0, 0),
+                energy: 0.0,
+                child_index: 0,
+            });
+
+            // Get combined object bounds
+            let bounds = {
+                let mut bb = BBox::new();
+                for obj in &objects[..] {
+                    bb |= lerp_slice(info_getter(obj).0, 0.5);
+                }
+                bb
+            };
+
+
+
+            // Partition objects.
+            let (split_index, split_axis) = sah_split(objects, &|obj_ref| info_getter(obj_ref).0);
+
+            // Create child nodes
+            let (_, c1_bounds) =
+                self.recursive_build(offset, depth + 1, &mut objects[..split_index], info_getter);
+            let (c2_index, c2_bounds) = self.recursive_build(offset + split_index,
+                                                             depth + 1,
+                                                             &mut objects[split_index..],
+                                                             info_getter);
+
+            // Determine bounds
+            // TODO: do merging without the temporary vec.
+            let bi = self.bounds.len();
+            let mut merged = Vec::new();
+            merge_slices_append(&self.bounds[c1_bounds.0..c1_bounds.1],
+                                &self.bounds[c2_bounds.0..c2_bounds.1],
+                                &mut merged,
+                                |b1, b2| *b1 | *b2);
+            self.bounds.extend(merged.drain(0..));
+
+            // Set node
+            let energy = self.nodes[me_index + 1].energy + self.nodes[c2_index].energy;
+            self.nodes[me_index] = Node {
+                is_leaf: false,
+                bounds_range: (bi, self.bounds.len()),
+                energy: energy,
+                child_index: c2_index,
+            };
+
+            return (me_index, (bi, self.bounds.len()));
+        }
+    }
+}
+
+
+impl LightAccel for LightTree {
+    fn select(&self,
+              inc: Vector,
+              pos: Point,
+              nor: Normal,
+              sc: &SurfaceClosure,
+              time: f32,
+              n: f32)
+              -> Option<(usize, f32, f32)> {
+        if self.nodes.len() == 0 {
+            return None;
+        }
+
+        let mut node_index = 0;
+        let mut tot_prob = 1.0;
+        let mut n = n;
+
+        // Calculates the selection probability for a node
+        let node_prob = |node_ref: &Node| {
+            let bounds = &self.bounds[node_ref.bounds_range.0..node_ref.bounds_range.1];
+            let bbox = lerp_slice(bounds, time);
+            let d = bbox.center() - pos;
+            let dist2 = d.length2();
+            let r = bbox.diagonal() * 0.5;
+            let r2 = r * r;
+            let inv_surface_area = 1.0 / r2;
+
+            let cos_theta_max = if dist2 <= r2 {
+                -1.0
+            } else {
+                let sin_theta_max2 = (r2 / dist2).min(1.0);
+                (1.0 - sin_theta_max2).sqrt()
+            };
+
+            // Get the approximate amount of light contribution from the
+            // composite light source.
+            let approx_contrib = sc.estimate_eval_over_solid_angle(inc, d, nor, cos_theta_max);
+
+            node_ref.energy * inv_surface_area * approx_contrib
+        };
+
+        // Traverse down the tree, keeping track of the relative probabilities
+        while !self.nodes[node_index].is_leaf {
+            // Calculate the relative probabilities of the two children
+            let (p1, p2) = {
+                let p1 = node_prob(&self.nodes[node_index + 1]);
+                let p2 = node_prob(&self.nodes[self.nodes[node_index].child_index]);
+                let total = p1 + p2;
+
+                if total <= 0.0 {
+                    (0.5, 0.5)
+                } else {
+                    (p1 / total, p2 / total)
+                }
+            };
+
+            if n <= p1 {
+                tot_prob *= p1;
+                node_index = node_index + 1;
+                n /= p1;
+            } else {
+                tot_prob *= p2;
+                node_index = self.nodes[node_index].child_index;
+                n = (n - p1) / p2;
+            }
+        }
+
+        // Found our light!
+        Some((self.nodes[node_index].child_index, tot_prob, n))
+    }
+}

src/light_accel/mod.rs

@@ -1,8 +1,21 @@
+mod light_tree;
+
+use math::{Vector, Point, Normal};
 use bbox::BBox;
+use shading::surface_closure::SurfaceClosure;
+
+pub use self::light_tree::LightTree;
 
 pub trait LightAccel {
     /// Returns (index_of_light, selection_pdf, whittled_n)
-    fn select(&self, n: f32) -> Option<(usize, f32, f32)>;
+    fn select(&self,
+              inc: Vector,
+              pos: Point,
+              nor: Normal,
+              sc: &SurfaceClosure,
+              time: f32,
+              n: f32)
+              -> Option<(usize, f32, f32)>;
 }
 
 #[derive(Debug, Clone)]
@@ -28,7 +41,16 @@ impl LightArray {
 }
 
 impl LightAccel for LightArray {
-    fn select(&self, n: f32) -> Option<(usize, f32, f32)> {
+    fn select(&self,
+              inc: Vector,
+              pos: Point,
+              nor: Normal,
+              sc: &SurfaceClosure,
+              time: f32,
+              n: f32)
+              -> Option<(usize, f32, f32)> {
+        let _ = (inc, pos, nor, sc, time); // Not using these, silence warnings
+
         assert!(n >= 0.0 && n <= 1.0);
 
         if self.indices.len() == 0 {

src/shading/surface_closure.rs

@@ -3,6 +3,7 @@ use color::{XYZ, SpectralSample, Color};
 use sampling::cosine_sample_hemisphere;
 use std::f32::consts::PI as PI_32;
 const INV_PI: f32 = 1.0 / PI_32;
+const H_PI: f32 = PI_32 / 2.0;
 
 #[derive(Debug, Copy, Clone)]
 pub enum SurfaceClosureUnion {
@@ -57,6 +58,18 @@ pub trait SurfaceClosure {
     /// out: The outgoing light direction.
     /// nor: The surface normal of the reflecting/transmitting surface point.
     fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal) -> f32;
+
+    /// Returns an estimate of the sum total energy that evaluate() would return
+    /// when 'out' is evaluated over a circular solid angle.
+    /// This is used for importance sampling, so does not need to be exact,
+    /// but it does need to be non-zero anywhere that an exact solution would
+    /// be non-zero.
+    fn estimate_eval_over_solid_angle(&self,
+                                      inc: Vector,
+                                      out: Vector,
+                                      nor: Normal,
+                                      cos_theta: f32)
+                                      -> f32;
 }
 
 
@@ -174,6 +187,18 @@ impl SurfaceClosure for EmitClosure {
 
         1.0
     }
+
+    fn estimate_eval_over_solid_angle(&self,
+                                      inc: Vector,
+                                      out: Vector,
+                                      nor: Normal,
+                                      cos_theta: f32)
+                                      -> f32 {
+        let _ = (inc, out, nor, cos_theta); // Not using these, silence warning
+
+        // TODO: what to do here?
+        unimplemented!()
+    }
 }
 
 
@@ -241,4 +266,63 @@ impl SurfaceClosure for LambertClosure {
 
         dot(nn, v).max(0.0) * INV_PI
     }
+
+    fn estimate_eval_over_solid_angle(&self,
+                                      inc: Vector,
+                                      out: Vector,
+                                      nor: Normal,
+                                      cos_theta: f32)
+                                      -> f32 {
+        assert!(cos_theta >= -1.0 && cos_theta <= 1.0);
+
+        // Analytically calculates lambert shading from a uniform light source
+        // subtending a circular solid angle.
+        // Only works for solid angle subtending equal to or less than a hemisphere.
+        //
+        // Formula taken from "Area Light Sources for Real-Time Graphics"
+        // by John M. Snyder
+        fn sphere_lambert(nlcos: f32, rcos: f32) -> f32 {
+            assert!(nlcos >= -1.0 && nlcos <= 1.0);
+            assert!(rcos >= 0.0 && rcos <= 1.0);
+
+            let nlsin: f32 = (1.0 - (nlcos * nlcos)).sqrt();
+            let rsin2: f32 = 1.0 - (rcos * rcos);
+            let rsin: f32 = rsin2.sqrt();
+            let ysin: f32 = rcos / nlsin;
+            let ycos2: f32 = 1.0 - (ysin * ysin);
+            let ycos: f32 = ycos2.sqrt();
+
+            let g: f32 = (-2.0 * nlsin * rcos * ycos) + H_PI - ysin.asin() + (ysin * ycos);
+            let h: f32 = nlcos * ((ycos * (rsin2 - ycos2).sqrt()) + (rsin2 * (ycos / rsin).asin()));
+
+            let nl: f32 = nlcos.acos();
+            let r: f32 = rcos.acos();
+
+            if nl < (H_PI - r) {
+                nlcos * rsin2
+            } else if nl < H_PI {
+                (nlcos * rsin2) + g - h
+            } else if nl < (H_PI + r) {
+                (g + h) * INV_PI
+            } else {
+                0.0
+            }
+        }
+
+        if cos_theta < 0.0 {
+            return 1.0;
+        } else {
+            let v = out.normalized();
+            let nn = if dot(nor.into_vector(), inc) <= 0.0 {
+                    nor.normalized()
+                } else {
+                    -nor.normalized()
+                }
+                .into_vector();
+
+            let cos_nv = dot(nn, v);
+
+            return sphere_lambert(cos_nv, cos_theta);
+        }
+    }
 }

todo.txt

@@ -1,4 +1,6 @@
 //- Implement support for multiple light sources in a scene.
+//- Implement light tree.
+- Implement proper light source selection for hierarchical instancing.
 - Implement shader parsing, and use in rendering.
 - Add gui display of rendering in progress.
 - Unit tests for scene parsing.