WIP multiple importance sampling.

Added method for intersecting finite light sources, and implemented
the method for SphereLight.

src/light/mod.rs

@@ -7,6 +7,8 @@ use std::fmt::Debug;
 use boundable::Boundable;
 use color::SpectralSample;
 use math::{Vector, Point, Matrix4x4};
+use ray::{Ray, AccelRay};
+use surface::SurfaceIntersection;
 
 pub use self::distant_disk_light::DistantDiskLight;
 pub use self::rectangle_light::RectangleLight;
@@ -38,13 +40,17 @@ pub trait LightSource: Boundable + Debug + Sync {
 
 
     /// Calculates the pdf of sampling the given
-    /// sample_dir/sample_u/sample_v from the given point arr.  This is used
+    /// `sample_dir`/`sample_u`/`sample_v` from the given point `arr`.  This is used
     /// primarily to calculate probabilities for multiple importance sampling.
     ///
     /// NOTE: this function CAN assume that sample_dir, sample_u, and sample_v
     /// are a valid sample for the light source (i.e. hits/lies on the light
     /// source).  No guarantees are made about the correctness of the return
     /// value if they are not valid.
+    ///
+    /// TODO: this probably shouldn't be part of the public interface.  In the
+    /// rest of the renderer, the PDF is always calculated by the `sample` and
+    /// and `intersect_rays` methods.
     fn sample_pdf(
         &self,
         space: &Matrix4x4,
@@ -65,6 +71,10 @@ pub trait LightSource: Boundable + Debug + Sync {
     ///     - v: Random parameter V.
     ///     - wavelength: The hero wavelength of light to sample at.
     ///     - time: The time to sample at.
+    ///
+    /// TODO: this probably shouldn't be part of the public interface.  In the
+    /// rest of the renderer, this is handled by the `sample` and
+    /// `intersect_rays` methods.
     fn outgoing(
         &self,
         space: &Matrix4x4,
@@ -88,6 +98,14 @@ pub trait LightSource: Boundable + Debug + Sync {
     /// source.  Note that this does not need to be exact: it is used for
     /// importance sampling.
     fn approximate_energy(&self) -> f32;
+
+    fn intersect_rays(
+        &self,
+        accel_rays: &mut [AccelRay],
+        wrays: &[Ray],
+        isects: &mut [SurfaceIntersection],
+        space: &[Matrix4x4],
+    );
 }
 
 

src/light/rectangle_light.rs

@@ -5,7 +5,9 @@ use boundable::Boundable;
 use color::{XYZ, SpectralSample, Color};
 use lerp::lerp_slice;
 use math::{Vector, Point, Matrix4x4};
+use ray::{Ray, AccelRay};
 use sampling::{spherical_triangle_solid_angle, uniform_sample_spherical_triangle};
+use surface::SurfaceIntersection;
 
 use super::LightSource;
 
@@ -176,6 +178,16 @@ impl<'a> LightSource for RectangleLight<'a> {
         ) / self.colors.len() as f32;
         color.y
     }
+
+    fn intersect_rays(
+        &self,
+        accel_rays: &mut [AccelRay],
+        wrays: &[Ray],
+        isects: &mut [SurfaceIntersection],
+        space: &[Matrix4x4],
+    ) {
+        unimplemented!()
+    }
 }
 
 impl<'a> Boundable for RectangleLight<'a> {

src/light/sphere_light.rs

@@ -6,8 +6,11 @@ use bbox::BBox;
 use boundable::Boundable;
 use color::{XYZ, SpectralSample, Color};
 use lerp::lerp_slice;
-use math::{Vector, Point, Matrix4x4, coordinate_system_from_vector};
+use math::{Vector, Point, Matrix4x4, dot, coordinate_system_from_vector};
+use ray::{Ray, AccelRay};
 use sampling::{uniform_sample_cone, uniform_sample_cone_pdf, uniform_sample_sphere};
+use surface::{SurfaceIntersection, SurfaceIntersectionData};
+use shading::surface_closure::{SurfaceClosureUnion, EmitClosure};
 
 use super::LightSource;
 
@@ -178,6 +181,139 @@ impl<'a> LightSource for SphereLight<'a> {
         ) / self.colors.len() as f32;
         color.y
     }
+
+    fn intersect_rays(
+        &self,
+        accel_rays: &mut [AccelRay],
+        wrays: &[Ray],
+        isects: &mut [SurfaceIntersection],
+        space: &[Matrix4x4],
+    ) {
+        for r in accel_rays.iter_mut() {
+            let wr = &wrays[r.id as usize];
+
+            // Get the transform space
+            let xform = lerp_slice(space, r.time);
+
+            // Get the radius of the sphere at the ray's time
+            let radius = lerp_slice(self.radii, r.time); // Radius of the sphere
+
+            // Get the ray origin and direction in local space
+            let orig = r.orig.into_vector();
+            let dir = wr.dir * xform;
+
+            // Code adapted to Rust from https://github.com/Tecla/Rayito
+            // Ray-sphere intersection can result in either zero, one or two points
+            // of intersection.  It turns into a quadratic equation, so we just find
+            // the solution using the quadratic formula.  Note that there is a
+            // slightly more stable form of it when computing it on a computer, and
+            // we use that method to keep everything accurate.
+
+            // Calculate quadratic coeffs
+            let a = dir.length2();
+            let b = 2.0 * dot(dir, orig);
+            let c = orig.length2() - (radius * radius);
+
+            let discriminant = (b * b) - (4.0 * a * c);
+            if discriminant < 0.0 {
+                // Discriminant less than zero?  No solution => no intersection.
+                continue;
+            }
+            let discriminant = discriminant.sqrt();
+
+            // Compute a more stable form of our param t (t0 = q/a, t1 = c/q)
+            // q = -0.5 * (b - sqrt(b * b - 4.0 * a * c)) if b < 0, or
+            // q = -0.5 * (b + sqrt(b * b - 4.0 * a * c)) if b >= 0
+            let q = if b < 0.0 {
+                -0.5 * (b - discriminant)
+            } else {
+                -0.5 * (b + discriminant)
+            };
+
+            // Get our final parametric values
+            let mut t0 = q / a;
+            let mut t1 = if q != 0.0 { c / q } else { r.max_t };
+
+            // Swap them so they are ordered right
+            if t0 > t1 {
+                use std::mem::swap;
+                swap(&mut t0, &mut t1);
+            }
+
+            // Check our intersection for validity against this ray's extents
+            if t0 > r.max_t || t1 <= 0.0 {
+                // Didn't hit because shere is entirely outside of ray's extents
+                continue;
+            }
+
+            let t = if t0 > 0.0 {
+                t0
+            } else if t1 <= r.max_t {
+                t1
+            } else {
+                // Didn't hit because ray is entirely within the sphere, and
+                // therefore doesn't hit its surface.
+                continue;
+            };
+
+            // We hit the sphere, so calculate intersection info.
+            if r.is_occlusion() {
+                isects[r.id as usize] = SurfaceIntersection::Occlude;
+                r.mark_done();
+            } else {
+                let inv_xform = xform.inverse();
+
+                // Position is calculated from the local-space ray and t, and then
+                // re-projected onto the surface of the sphere.
+                let t_pos = orig + (dir * t);
+                let unit_pos = t_pos.normalized();
+                let pos = (unit_pos * radius * inv_xform).into_point();
+
+                // TODO: proper error bounds.
+                let pos_err = 0.001;
+
+                let normal = unit_pos.into_normal() * inv_xform;
+
+                let intersection_data = SurfaceIntersectionData {
+                    incoming: wr.dir,
+                    t: t,
+                    pos: pos,
+                    pos_err: pos_err,
+                    nor: normal,
+                    nor_g: normal,
+                    uv: (0.0, 0.0), // TODO
+                    local_space: xform,
+                    sample_pdf: self.sample_pdf(
+                        &xform,
+                        wr.orig,
+                        wr.dir,
+                        0.0,
+                        0.0,
+                        wr.wavelength,
+                        r.time,
+                    ),
+                };
+
+                let closure = {
+                    let inv_surface_area = (1.0 / (4.0 * PI_64 * radius as f64 * radius as f64)) as
+                        f32;
+                    let color = lerp_slice(self.colors, r.time).to_spectral_sample(
+                        wr.wavelength,
+                    ) * inv_surface_area;
+                    SurfaceClosureUnion::EmitClosure(EmitClosure::new(color))
+                };
+
+                // Fill in intersection
+                isects[r.id as usize] = SurfaceIntersection::Hit {
+                    intersection_data: intersection_data,
+                    closure: closure,
+                };
+
+                // Set ray's max t
+                r.max_t = t;
+            }
+        }
+    }
 }
 
 impl<'a> Boundable for SphereLight<'a> {

src/surface/mod.rs

@@ -46,4 +46,5 @@ pub struct SurfaceIntersectionData {
     pub local_space: Matrix4x4, // Matrix from global space to local space
     pub t: f32, // Ray t-value at the intersection point
     pub uv: (f32, f32), // 2d surface parameters
+    pub sample_pdf: f32, // The PDF of getting this point by explicitly sampling the surface
 }

src/surface/triangle_mesh.rs

@@ -244,6 +244,7 @@ impl<'a> Surface for TriangleMesh<'a> {
                                     nor_g: geo_normal,
                                     uv: (0.0, 0.0), // TODO
                                     local_space: mat_space,
+                                    sample_pdf: 0.0,
                                 };
 
                                 // Fill in intersection data