Multiple importance sampling is now basically working.

Intersecting rectangular lights still isn't implemented, and there
are likely bugs in the MIS implementation, but it's more-or-less
working!

src/light/rectangle_light.rs

@@ -192,7 +192,8 @@ impl<'a> Surface for RectangleLight<'a> {
         space: &[Matrix4x4],
     ) {
         let _ = (accel_rays, wrays, isects, shader, space);
-        unimplemented!()
+        // TODO
+        // unimplemented!()
     }
 }
 

src/light/sphere_light.rs

@@ -15,6 +15,11 @@ use surface::{Surface, SurfaceIntersection, SurfaceIntersectionData};
 
 use super::SurfaceLight;
 
+// TODO: use proper error bounds for sample generation to avoid self-shadowing
+// instead of these fudge factors.
+const SAMPLE_RADIUS_EXPAND_FACTOR: f32 = 1.001;
+const SAMPLE_RADIUS_SHRINK_FACTOR: f32 = 0.99;
+
 // TODO: handle case where radius = 0.0.
 
 #[derive(Copy, Clone, Debug)]
@@ -140,14 +145,23 @@ impl<'a> SurfaceLight for SphereLight<'a> {
             );
 
             // Calculate the final values and return everything.
-            let shadow_vec = ((x * sample.x()) + (y * sample.y()) + (z * sample.z())) *
+            let shadow_vec = {
-                space.inverse();
+                let sample_vec = (x * sample.x()) + (y * sample.y()) + (z * sample.z());
+                let sample_point = (arr + sample_vec).into_vector().normalized() * radius as f32;
+                let adjusted_sample_point = sample_point * SAMPLE_RADIUS_EXPAND_FACTOR;
+                (adjusted_sample_point.into_point() - arr) * space.inverse()
+            };
             let pdf = uniform_sample_cone_pdf(cos_theta_max);
             let spectral_sample = (col * surface_area_inv as f32).to_spectral_sample(wavelength);
             return (spectral_sample, shadow_vec, pdf as f32);
         } else {
             // If we're inside the sphere, there's light from every direction.
-            let shadow_vec = uniform_sample_sphere(u, v) * space.inverse();
+            let shadow_vec = {
+                let sample_vec = uniform_sample_sphere(u, v);
+                let sample_point = (arr + sample_vec).into_vector().normalized() * radius as f32;
+                let adjusted_sample_point = sample_point * SAMPLE_RADIUS_SHRINK_FACTOR;
+                (adjusted_sample_point.into_point() - arr) * space.inverse()
+            };
             let pdf = 1.0 / (4.0 * PI_64);
             let spectral_sample = (col * surface_area_inv as f32).to_spectral_sample(wavelength);
             return (spectral_sample, shadow_vec, pdf as f32);

src/main.rs

@@ -44,6 +44,7 @@ mod image;
 mod lerp;
 mod light;
 mod math;
+mod mis;
 mod parse;
 mod ray;
 mod renderer;

src/mis.rs

@@ -0,0 +1,13 @@
+#![allow(dead_code)]
+
+pub fn balance_heuristic(a: f32, b: f32) -> f32 {
+    let mis_fac = a / (a + b);
+    a / mis_fac
+}
+
+pub fn power_heuristic(a: f32, b: f32) -> f32 {
+    let a2 = a * a;
+    let b2 = b * b;
+    let mis_fac = a2 / (a2 + b2);
+    a / mis_fac
+}

src/renderer.rs

@@ -19,6 +19,7 @@ use hash::hash_u32;
 use hilbert;
 use image::Image;
 use math::{fast_logit, upper_power_of_two};
+use mis::power_heuristic;
 use ray::Ray;
 use scene::Scene;
 use surface;
@@ -372,8 +373,9 @@ pub struct LightPath {
     wavelength: f32,
 
     next_bounce_ray: Option<Ray>,
-    next_attentuation_fac: Float4,
+    next_attenuation_fac: Float4,
 
+    closure_sample_pdf: f32,
     light_attenuation: Float4,
     pending_color_addition: Float4,
     color: Float4,
@@ -401,8 +403,9 @@ impl LightPath {
                 wavelength: wavelength,
 
                 next_bounce_ray: None,
-                next_attentuation_fac: Float4::splat(1.0),
+                next_attenuation_fac: Float4::splat(1.0),
 
+                closure_sample_pdf: 1.0,
                 light_attenuation: Float4::splat(1.0),
                 pending_color_addition: Float4::splat(0.0),
                 color: Float4::splat(0.0),
@@ -449,10 +452,20 @@ impl LightPath {
                     // - Terminate the path.
                     use shading::surface_closure::SurfaceClosureUnion;
                     if let &SurfaceClosureUnion::EmitClosure(ref clsr) = closure {
-                        self.color += clsr.emitted_color().e * self.light_attenuation;
+                        if let LightPathEvent::CameraRay = self.event {
+                            self.color += clsr.emitted_color().e;
+                        } else {
+                            let mis_pdf =
+                                power_heuristic(self.closure_sample_pdf, idata.sample_pdf);
+                            self.color += clsr.emitted_color().e * self.light_attenuation / mis_pdf;
+                        };
+
                         return false;
                     }
 
+                    // Roll the previous closure pdf into the attentuation
+                    self.light_attenuation /= self.closure_sample_pdf;
+
                     // Prepare light ray
                     let light_n = self.next_lds_samp();
                     let light_uvw = (
@@ -482,11 +495,20 @@ impl LightPath {
                                 material.evaluate(ray.dir, shadow_vec, idata.nor, idata.nor_g);
 
                             if attenuation.e.h_max() > 0.0 {
+                                // Calculate MIS
+                                let closure_pdf = material.sample_pdf(
+                                    ray.dir,
+                                    shadow_vec,
+                                    idata.nor,
+                                    idata.nor_g,
+                                );
+                                let light_mis_pdf = power_heuristic(light_pdf, closure_pdf);
+
                                 // Calculate and store the light that will be contributed
                                 // to the film plane if the light is not in shadow.
                                 self.pending_color_addition = light_color.e * attenuation.e *
                                     self.light_attenuation /
-                                    (light_pdf * light_sel_pdf);
+                                    (light_mis_pdf * light_sel_pdf);
 
                                 // Calculate the shadow ray for testing if the light is
                                 // in shadow or not.
@@ -536,7 +558,8 @@ impl LightPath {
                         if (pdf > 0.0) && (filter.e.h_max() > 0.0) {
                             // Account for the additional light attenuation from
                             // this bounce
-                            self.next_attentuation_fac = filter.e / pdf;
+                            self.next_attenuation_fac = filter.e;
+                            self.closure_sample_pdf = pdf;
 
                             // Calculate the ray for this bounce
                             let offset_pos = robust_ray_origin(
@@ -564,7 +587,7 @@ impl LightPath {
                     } else if do_bounce {
                         *ray = self.next_bounce_ray.unwrap();
                         self.event = LightPathEvent::BounceRay;
-                        self.light_attenuation *= self.next_attentuation_fac;
+                        self.light_attenuation *= self.next_attenuation_fac;
                         return true;
                     } else {
                         return false;
@@ -575,7 +598,8 @@ impl LightPath {
                         .world
                         .background_color
                         .to_spectral_sample(self.wavelength)
-                        .e * self.light_attenuation;
+                        .e * self.light_attenuation /
+                        self.closure_sample_pdf;
                     return false;
                 }
             }
@@ -592,7 +616,7 @@ impl LightPath {
                 // Set up for the next bounce, if any
                 if let Some(ref nbr) = self.next_bounce_ray {
                     *ray = *nbr;
-                    self.light_attenuation *= self.next_attentuation_fac;
+                    self.light_attenuation *= self.next_attenuation_fac;
                     self.event = LightPathEvent::BounceRay;
                     return true;
                 } else {

src/tracer.rs

@@ -185,7 +185,7 @@ impl<'a> TracerInner<'a> {
     ) {
         match *obj {
             Object::Surface(surface) => {
-                let unassigned_shader = SimpleSurfaceShader::Lambert {
+                let unassigned_shader = SimpleSurfaceShader::Emit {
                     color: XYZ::from_tuple(rec709_to_xyz((1.0, 0.0, 1.0))),
                 };
                 let shader = surface_shader.unwrap_or(&unassigned_shader);
@@ -199,8 +199,19 @@ impl<'a> TracerInner<'a> {
                 );
             }
 
-            Object::SurfaceLight(_) => {
+            Object::SurfaceLight(surface) => {
-                // TODO
+                // Lights don't use shaders
+                let bogus_shader = SimpleSurfaceShader::Emit {
+                    color: XYZ::from_tuple(rec709_to_xyz((1.0, 0.0, 1.0))),
+                };
+
+                surface.intersect_rays(
+                    rays,
+                    wrays,
+                    &mut self.isects,
+                    &bogus_shader,
+                    self.xform_stack.top(),
+                );
             }
         }
     }