Sample lights by power

Produces somewhat lower noise in the staircase scene, though the
fireflies are still there.

src/bvh.rs

@@ -31,6 +31,7 @@ pub enum SplitMethod {
 #[derive(Debug, PartialEq)]
 pub struct Bvh {
     pub root: BvhNode,
+    pub bounds: Bounds,
 }
 
 impl Bvh {
@@ -49,7 +50,9 @@ impl Bvh {
             SplitMethod::SAH => BvhNode::from_sah_splitting(&mut primitive_infos),
         };
 
-        Bvh { root }
+        let bounds: Bounds = primitive_infos.iter().map(|i| i.bounds).sum();
+
+        Bvh { root, bounds }
     }
 
     pub fn intersect(&self, ray: &mut Ray) -> Option<PrimitiveIntersection> {

src/color.rs

@@ -160,6 +160,11 @@ impl Mul<Color> for Color {
 
 impl Display for Color {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(f, "({},{},{})", self.r, self.g, self.b)
+        let precision = f.precision().unwrap_or(8);
+        write!(
+            f,
+            "({:.precision$},{:.precision$},{:.precision$})",
+            self.r, self.g, self.b
+        )
     }
 }

src/light.rs

@@ -48,9 +48,6 @@ pub struct LightSample {
     pub shadow_ray: Ray,
 }
 
-// TODO: This should be computed using the scene's bounds
-const WORLD_RADIUS: f64 = 1e6;
-
 impl Light {
     /// Samples the light arriving at a given point from this light source.
     ///
@@ -168,12 +165,50 @@ impl Light {
         }
     }
 
-    pub fn power(self: &Self) -> Color {
+    pub fn power(self: &Self, world_radius: f64) -> Color {
         match &self {
             Light::Point { intensity, .. } => *intensity * 4.0 * PI,
-            Light::Distant { intensity, .. } => *intensity * PI * WORLD_RADIUS * WORLD_RADIUS,
-            Light::Infinite { intensity, .. } => *intensity * PI * WORLD_RADIUS * WORLD_RADIUS,
+            Light::Distant { intensity, .. } => *intensity * PI * world_radius * world_radius,
+            Light::Infinite { intensity, .. } => *intensity * PI * world_radius * world_radius,
             Light::Area { emittance, shape } => *emittance * PI * shape.area(),
         }
     }
 }
+
+/// Samples lights in proportion to their power
+#[derive(Debug, PartialEq)]
+pub struct LightSampler {
+    cdfs: Vec<f64>,
+    total_cdf: f64,
+}
+
+impl LightSampler {
+    pub fn new(lights: &Vec<Arc<Light>>, world_radius: f64) -> Self {
+        let mut total_cdf = 0.0;
+        let mut cdfs = vec![];
+        for light in lights.iter() {
+            let power = light.power(world_radius);
+            let pdf = (power.r + power.g + power.b) / 3.0;
+            total_cdf += pdf;
+            cdfs.push(total_cdf);
+        }
+        Self { cdfs, total_cdf }
+    }
+
+    // TODO: Try the "AliasTable" method to do this in constant time
+    pub fn sample(&self, sample: Sample1d) -> (usize, f64) {
+        let u = sample.take() * self.total_cdf;
+        let idx = self
+            .cdfs
+            .binary_search_by(|probe| probe.total_cmp(&u))
+            .unwrap_or_else(|idx| idx);
+
+        let pdf = if idx > 0 {
+            self.cdfs[idx] - self.cdfs[idx - 1]
+        } else {
+            self.cdfs[idx]
+        };
+
+        (idx, pdf / self.total_cdf)
+    }
+}

src/scene.rs

@@ -6,7 +6,7 @@ use crate::{
     bvh::{Bvh, SplitMethod},
     camera::Camera,
     intersection::PrimitiveIntersection,
-    light::Light,
+    light::{Light, LightSampler},
     primitive::Primitive,
     ray::Ray,
 };
@@ -17,6 +17,7 @@ pub struct Scene {
     pub num_samples: usize,
     pub camera: Camera,
     pub lights: Vec<Arc<Light>>,
+    pub light_sampler: LightSampler,
     bvh: Bvh,
 }
 
@@ -37,11 +38,16 @@ impl Scene {
         );
         let bvh = Bvh::new(primitives, SplitMethod::SAH);
         debug!("BVH constructed in {:?}", start.elapsed());
+
+        let world_radius = bvh.bounds.diagonal().magnitude() * 0.5;
+        let light_sampler = LightSampler::new(&lights, world_radius);
+
         Self {
             max_depth,
             num_samples,
             camera,
             lights,
+            light_sampler,
             bvh,
         }
     }

src/trace.rs

@@ -226,8 +226,7 @@ where
 
         // Estimate the contribution from a path that ends here. We will reuse
         // the path without the terminator in the loop.
-        let light_pdf = 1.0 / scene.lights.len() as f64;
-        let light_index = (path_samples.light_index.take() * scene.lights.len() as f64) as usize;
+        let (light_index, light_pdf) = scene.light_sampler.sample(path_samples.light_index);
         let light = &scene.lights[light_index];
         L +=
             beta * estimate_direct(