From 3e7b142cd879ae889c72881a80a019643a7f6eae Mon Sep 17 00:00:00 2001 From: Nathan Vegdahl Date: Sun, 31 Jul 2016 11:14:33 -0700 Subject: [PATCH] Implemented light tree sampling, for better sampling of many lights. --- src/assembly.rs | 47 ++++++-- src/bbox.rs | 8 ++ src/light/mod.rs | 7 +- src/light/rectangle_light.rs | 6 + src/light/sphere_light.rs | 6 + src/light_accel/light_tree.rs | 200 +++++++++++++++++++++++++++++++++ src/light_accel/mod.rs | 26 ++++- src/shading/surface_closure.rs | 84 ++++++++++++++ todo.txt | 2 + 9 files changed, 371 insertions(+), 15 deletions(-) create mode 100644 src/light_accel/light_tree.rs diff --git a/src/assembly.rs b/src/assembly.rs index f595ea8..2728742 100644 --- a/src/assembly.rs +++ b/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 Some((light_i, sel_pdf, _)) = self.light_accel.select(n) { + if let &SurfaceIntersection::Hit { pos, incoming, nor, closure, .. } = intr { + 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()); + // Get list of instances that are for light sources. + // TODO: include assemblies that themselves contain light sources. + 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] { + true + } else { + false + } + } + + _ => false, + } + }) + .map(|&a| a) + .collect(); + // Build light accel - let mut light_instances = self.instances.clone(); - let light_accel = LightArray::new(&mut light_instances[..], |inst| { - match inst.instance_type { + 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(_) = self.objects[inst.data_index] { - Some((&bbs[bis[inst.id]..bis[inst.id + 1]], 1.0)) + if let Object::Light(ref light) = self.objects[inst.data_index] { + light.approximate_energy() } else { - None + 0.0 } } - _ => None, - } + // TODO: handle assemblies. + _ => 0.0, + }; + (bounds, energy) }); Assembly { diff --git a/src/bbox.rs b/src/bbox.rs index 01e028f..aa070e7 100644 --- a/src/bbox.rs +++ b/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() + } } diff --git a/src/light/mod.rs b/src/light/mod.rs index 0011ef4..21df3a4 100644 --- a/src/light/mod.rs +++ b/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; } diff --git a/src/light/rectangle_light.rs b/src/light/rectangle_light.rs index 49fe252..5e0fc4e 100644 --- a/src/light/rectangle_light.rs +++ b/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 { diff --git a/src/light/sphere_light.rs b/src/light/sphere_light.rs index 6f4a90c..9cc4c51 100644 --- a/src/light/sphere_light.rs +++ b/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 { diff --git a/src/light_accel/light_tree.rs b/src/light_accel/light_tree.rs new file mode 100644 index 0000000..211dfbc --- /dev/null +++ b/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, + bounds: Vec, + depth: usize, + bounds_cache: Vec, +} + +#[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)) + } +} diff --git a/src/light_accel/mod.rs b/src/light_accel/mod.rs index 633427d..a26916e 100644 --- a/src/light_accel/mod.rs +++ b/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 { diff --git a/src/shading/surface_closure.rs b/src/shading/surface_closure.rs index a4b22cf..ddd05da 100644 --- a/src/shading/surface_closure.rs +++ b/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); + } + } } diff --git a/todo.txt b/todo.txt index 680ed6f..7d02ef4 100644 --- a/todo.txt +++ b/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. \ No newline at end of file