Beginnings of global illumination.

There are still some obvious bugs in it that I haven't tracked down,
so the renders aren't correct yet.

src/color/mod.rs

@@ -93,6 +93,24 @@ impl AddAssign for SpectralSample {
     }
 }
 
+impl Mul for SpectralSample {
+    type Output = SpectralSample;
+    fn mul(self, rhs: SpectralSample) -> Self::Output {
+        assert!(self.hero_wavelength == rhs.hero_wavelength);
+        SpectralSample {
+            e: self.e * rhs.e,
+            hero_wavelength: self.hero_wavelength,
+        }
+    }
+}
+
+impl MulAssign for SpectralSample {
+    fn mul_assign(&mut self, rhs: SpectralSample) {
+        assert!(self.hero_wavelength == rhs.hero_wavelength);
+        self.e = self.e * rhs.e;
+    }
+}
+
 impl Mul<f32> for SpectralSample {
     type Output = SpectralSample;
     fn mul(self, rhs: f32) -> Self::Output {

src/main.rs

@@ -27,6 +27,7 @@ mod assembly;
 mod halton;
 mod sampling;
 mod color;
+mod shading;
 
 use std::mem;
 use std::io;

src/math/mod.rs

@@ -55,7 +55,25 @@ pub fn coordinate_system_from_vector(v: Vector) -> (Vector, Vector, Vector) {
     (v, v2, v3)
 }
 
+/// Simple mapping of a vector that exists in a z-up space to
+/// the space of another vector who's direction is considered
+/// z-up for the purpose.
+/// Obviously this doesn't care about the direction _around_
+/// the z-up, although it will be sufficiently consistent for
+/// isotropic sampling purposes.
+///
+/// from: The vector we're transforming.
+/// toz: The vector whose space we are transforming "from" into.
+///
+/// Returns he transformed vector.
+pub fn zup_to_vec(from: Vector, toz: Vector) -> Vector {
+    let (toz, tox, toy) = coordinate_system_from_vector(toz.normalized());
 
+    // Use simple linear algebra to convert the "from"
+    // vector to a space composed of tox, toy, and toz
+    // as the x, y, and z axes.
+    (tox * from[0]) + (toy * from[1]) + (toz * from[2])
+}
 
 /// The logit function, scaled to approximate the probit function.
 ///

src/renderer.rs

@@ -13,11 +13,12 @@ use ray::Ray;
 use assembly::Object;
 use tracer::Tracer;
 use halton;
-use math::{Matrix4x4, dot, fast_logit};
+use math::{Matrix4x4, fast_logit};
 use image::Image;
 use surface;
 use scene::Scene;
 use color::{Color, XYZ, SpectralSample, map_0_1_to_wavelength};
+use shading::surface_closure::{SurfaceClosure, LambertClosure};
 
 #[derive(Debug)]
 pub struct Renderer {
@@ -182,6 +183,7 @@ pub struct LightPath {
     round: u32,
     time: f32,
     wavelength: f32,
+    interaction: surface::SurfaceIntersection,
     light_attenuation: SpectralSample,
     pending_color_addition: SpectralSample,
     color: SpectralSample,
@@ -203,6 +205,7 @@ impl LightPath {
             round: 0,
             time: time,
             wavelength: wavelength,
+            interaction: surface::SurfaceIntersection::Miss,
             light_attenuation: SpectralSample::from_value(1.0, wavelength),
             pending_color_addition: SpectralSample::new(wavelength),
             color: SpectralSample::new(wavelength),
@@ -222,77 +225,114 @@ impl LightPath {
     }
 
     fn next(&mut self, scene: &Scene, isect: &surface::SurfaceIntersection, ray: &mut Ray) -> bool {
-        match self.round {
+        self.round += 1;
-            // Result of camera rays, prepare light rays
-            0 => {
-                self.round += 1;
-                if let &surface::SurfaceIntersection::Hit { t: _,
-                                                            pos,
-                                                            nor,
-                                                            local_space: _,
-                                                            uv: _ } = isect {
-                    // Hit something!  Do lighting!
-                    if scene.root.light_accel.len() > 0 {
-                        // Get the light and the mapping to its local space
-                        let (light, space) = {
-                            let l1 = &scene.root.objects[scene.root.light_accel[0].data_index];
-                            let light = if let &Object::Light(ref light) = l1 {
-                                light
-                            } else {
-                                panic!()
-                            };
-                            let space = if let Some((start, end)) = scene.root.light_accel[0]
-                                .transform_indices {
-                                lerp_slice(&scene.root.xforms[start..end], self.time)
-                            } else {
-                                Matrix4x4::new()
-                            };
-                            (light, space)
-                        };
 
+        // Result of shading ray, prepare light ray
+        if self.round % 2 == 1 {
+            if let &surface::SurfaceIntersection::Hit { t: _,
+                                                        incoming: _,
+                                                        pos,
+                                                        nor,
+                                                        local_space: _,
+                                                        uv: _ } = isect {
+                // Hit something!  Do the stuff
+                self.interaction = *isect; // Store interaction for use in next phase
+
+                // Prepare light ray
+                if scene.root.light_accel.len() > 0 {
+                    // Get the light and the mapping to its local space
+                    let (light, space) = {
+                        let l1 = &scene.root.objects[scene.root.light_accel[0].data_index];
+                        let light = if let &Object::Light(ref light) = l1 {
+                            light
+                        } else {
+                            panic!()
+                        };
+                        let space = if let Some((start, end)) = scene.root.light_accel[0]
+                            .transform_indices {
+                            lerp_slice(&scene.root.xforms[start..end], self.time)
+                        } else {
+                            Matrix4x4::new()
+                        };
+                        (light, space)
+                    };
+
+                    // Sample the light
+                    let (light_color, shadow_vec, light_pdf) = {
                         let lu = self.next_lds_samp();
                         let lv = self.next_lds_samp();
-                        // TODO: store incident light info and pdf, and use them properly
+                        light.sample(&space, pos, lu, lv, self.wavelength, self.time)
-                        let (light_color, shadow_vec, light_pdf) =
+                    };
-                            light.sample(&space, pos, lu, lv, self.wavelength, self.time);
 
-                        let rnor = if dot(nor.into_vector(), ray.dir) > 0.0 {
+                    // Calculate and store the light that will be contributed
-                            -nor.into_vector().normalized()
+                    // to the film plane if the light is not in shadow.
-                        } else {
+                    self.pending_color_addition = {
-                            nor.into_vector().normalized()
+                        let material = LambertClosure::new(XYZ::new(0.8, 0.8, 0.8));
-                        };
+                        let la = material.evaluate(ray.dir, shadow_vec, nor, self.wavelength);
-                        let la = dot(rnor, shadow_vec.normalized()).max(0.0);
+                        light_color * la * self.light_attenuation / light_pdf
-                        // self.light_attenuation = SpectralSample::from_value(la);
+                    };
-                        self.pending_color_addition = light_color * la / light_pdf;
-                        *ray = Ray::new(pos + rnor * 0.0001,
-                                        shadow_vec - rnor * 0.0001,
-                                        self.time,
-                                        true);
 
-                        return true;
+                    // Calculate the shadow ray for testing if the light is
-                    } else {
+                    // in shadow or not.
-                        return false;
+                    // TODO: use proper ray offsets for avoiding self-shadowing
-                    }
+                    // rather than this hacky stupid stuff.
+                    *ray = Ray::new(pos + shadow_vec.normalized() * 0.0001,
+                                    shadow_vec,
+                                    self.time,
+                                    true);
+
+                    return true;
                 } else {
-                    // Didn't hit anything, so background color
-                    let xyz = XYZ::new(0.02, 0.02, 0.02);
-                    self.color += xyz.to_spectral_sample(self.wavelength);
                     return false;
                 }
-
+            } else {
-            }
+                // Didn't hit anything, so background color
-
+                let xyz = XYZ::new(0.0, 0.0, 0.0);
-            // Result of light rays
+                self.color += xyz.to_spectral_sample(self.wavelength);
-            1 => {
-                self.round += 1;
-                if let &surface::SurfaceIntersection::Miss = isect {
-                    self.color += self.pending_color_addition;
-                }
                 return false;
             }
+        }
+        // Result of light ray, prepare shading ray
+        else if self.round % 2 == 0 {
+            // If the light was not in shadow, add it's light to the film
+            // plane.
+            if let &surface::SurfaceIntersection::Miss = isect {
+                self.color += self.pending_color_addition;
+            }
 
+            // Calculate bounced lighting!
+            if self.round < 6 {
+                if let surface::SurfaceIntersection::Hit { t: _,
+                                                           pos,
+                                                           incoming,
+                                                           nor,
+                                                           local_space: _,
+                                                           uv: _ } = self.interaction {
+                    // Sample material
+                    let (dir, filter, pdf) = {
+                        let material = LambertClosure::new(XYZ::new(0.8, 0.8, 0.8));
+                        let u = self.next_lds_samp();
+                        let v = self.next_lds_samp();
+                        material.sample(incoming, nor, (u, v), self.wavelength)
+                    };
+
+                    // Account for the additional light attenuation from
+                    // this bounce
+                    self.light_attenuation *= filter / pdf;
+
+                    // Calculate the ray for this bounce
+                    *ray = Ray::new(pos + dir.normalized() * 0.0001, dir, self.time, false);
+
+                    return true;
+                } else {
+                    return false;
+                }
+            } else {
+                return false;
+            }
+        } else {
             // TODO
-            _ => unimplemented!(),
+            unimplemented!()
         }
     }
 }

src/shading/mod.rs

@@ -0,0 +1 @@
+pub mod surface_closure;

src/shading/surface_closure.rs

@@ -0,0 +1,229 @@
+use math::{Vector, Normal, dot, zup_to_vec};
+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;
+
+/// Trait for surface closures.
+pub trait SurfaceClosure: Copy {
+    /// Returns whether the closure has a delta distribution or not.
+    fn is_delta(&self) -> bool;
+
+    /// Given an incoming ray and sample values, generates an outgoing ray and
+    /// color filter.
+    ///
+    /// inc: Incoming light direction.
+    /// nor: The surface normal at the surface point.
+    /// uv:  The sampling values.
+    /// wavelength: The wavelength of light to sample at.
+    ///
+    /// Returns a tuple with the generated outgoing light direction, color filter, and pdf.
+    fn sample(&self,
+              inc: Vector,
+              nor: Normal,
+              uv: (f32, f32),
+              wavelength: f32)
+              -> (Vector, SpectralSample, f32);
+
+    /// Evaluates the closure for the given incoming and outgoing rays.
+    ///
+    /// inc: The incoming light direction.
+    /// out: The outgoing light direction.
+    /// nor: The surface normal of the reflecting/transmitting surface point.
+    /// wavelength: The wavelength of light to evaluate for.
+    ///
+    /// Returns the resulting filter color.
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, wavelength: f32) -> SpectralSample;
+
+    /// Returns the pdf for the given 'in' direction producing the given 'out'
+    /// direction with the given differential geometry.
+    ///
+    /// inc: The incoming light direction.
+    /// 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;
+}
+
+
+/// Utility function that calculates the fresnel reflection factor of a given
+/// incoming ray against a surface with the given ior outside/inside ratio.
+///
+/// ior_ratio: The ratio of the outside material ior (probably 1.0 for air)
+///            over the inside ior.
+/// c: The cosine of the angle between the incoming light and the
+///    surface's normal.  Probably calculated e.g. with a normalized
+///    dot product.
+#[allow(dead_code)]
+fn dielectric_fresnel(ior_ratio: f32, c: f32) -> f32 {
+    let g = (ior_ratio - 1.0 + (c * c)).sqrt();
+
+    let f1 = g - c;
+    let f2 = g + c;
+    let f3 = (f1 * f1) / (f2 * f2);
+
+    let f4 = (c * f2) - 1.0;
+    let f5 = (c * f1) + 1.0;
+    let f6 = 1.0 + ((f4 * f4) / (f5 * f5));
+
+    return 0.5 * f3 * f6;
+}
+
+
+/// Schlick's approximation of the fresnel reflection factor.
+///
+/// Same interface as dielectric_fresnel(), above.
+#[allow(dead_code)]
+fn schlick_fresnel(ior_ratio: f32, c: f32) -> f32 {
+    let f1 = (1.0 - ior_ratio) / (1.0 + ior_ratio);
+    let f2 = f1 * f1;
+    let c1 = 1.0 - c;
+    let c2 = c1 * c1;
+    return f2 + ((1.0 - f2) * c1 * c2 * c2);
+}
+
+
+/// Utility function that calculates the fresnel reflection factor of a given
+/// incoming ray against a surface with the given normal-reflectance factor.
+///
+/// frensel_fac: The ratio of light reflected back if the ray were to
+///              hit the surface head-on (perpendicular to the surface).
+/// c The cosine of the angle between the incoming light and the
+///   surface's normal.  Probably calculated e.g. with a normalized
+///   dot product.
+#[allow(dead_code)]
+fn dielectric_fresnel_from_fac(fresnel_fac: f32, c: f32) -> f32 {
+    let tmp1 = fresnel_fac.sqrt() - 1.0;
+
+    // Protect against divide by zero.
+    if tmp1.abs() < 0.000001 {
+        return 1.0;
+    }
+
+    // Find the ior ratio
+    let tmp2 = (-2.0 / tmp1) - 1.0;
+    let ior_ratio = tmp2 * tmp2;
+
+    // Calculate fresnel factor
+    return dielectric_fresnel(ior_ratio, c);
+}
+
+
+/// Schlick's approximation version of dielectric_fresnel_from_fac() above.
+#[allow(dead_code)]
+fn schlick_fresnel_from_fac(frensel_fac: f32, c: f32) -> f32 {
+    let c1 = 1.0 - c;
+    let c2 = c1 * c1;
+    return frensel_fac + ((1.0 - frensel_fac) * c1 * c2 * c2);
+}
+
+
+/// Emit closure.
+///
+/// NOTE: this needs to be handled specially by the integrator!  It does not
+/// behave like a standard closure!
+#[derive(Debug, Copy, Clone)]
+pub struct EmitClosure {
+    col: XYZ,
+}
+
+impl EmitClosure {
+    pub fn emitted_color(&self, wavelength: f32) -> SpectralSample {
+        self.col.to_spectral_sample(wavelength)
+    }
+}
+
+impl SurfaceClosure for EmitClosure {
+    fn is_delta(&self) -> bool {
+        false
+    }
+
+    fn sample(&self,
+              inc: Vector,
+              nor: Normal,
+              uv: (f32, f32),
+              wavelength: f32)
+              -> (Vector, SpectralSample, f32) {
+        let _ = (inc, nor, uv); // Not using these, silence warning
+
+        (Vector::new(0.0, 0.0, 0.0), SpectralSample::new(wavelength), 1.0)
+    }
+
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, wavelength: f32) -> SpectralSample {
+        let _ = (inc, out, nor); // Not using these, silence warning
+
+        SpectralSample::new(wavelength)
+    }
+
+    fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal) -> f32 {
+        let _ = (inc, out, nor); // Not using these, silence warning
+
+        1.0
+    }
+}
+
+
+/// Lambertian surface closure
+#[derive(Debug, Copy, Clone)]
+pub struct LambertClosure {
+    col: XYZ,
+}
+
+impl LambertClosure {
+    pub fn new(col: XYZ) -> LambertClosure {
+        LambertClosure { col: col }
+    }
+}
+
+impl SurfaceClosure for LambertClosure {
+    fn is_delta(&self) -> bool {
+        false
+    }
+
+    fn sample(&self,
+              inc: Vector,
+              nor: Normal,
+              uv: (f32, f32),
+              wavelength: f32)
+              -> (Vector, SpectralSample, f32) {
+        let nn = if dot(nor.into_vector().normalized(), inc.normalized()) <= 0.0 {
+                nor.normalized()
+            } else {
+                -nor.normalized()
+            }
+            .into_vector();
+
+        // Generate a random ray direction in the hemisphere
+        // of the surface.
+        let dir = cosine_sample_hemisphere(uv.0, uv.1);
+        let pdf = dir[2] * INV_PI;
+        let out = zup_to_vec(dir, nn);
+        let filter = self.evaluate(inc, out, nor, wavelength);
+
+        (dir, filter, pdf)
+    }
+
+    fn evaluate(&self, inc: Vector, out: Vector, nor: Normal, wavelength: f32) -> SpectralSample {
+        let v = out.normalized();
+        let nn = if dot(nor.into_vector(), inc) <= 0.0 {
+                nor.normalized()
+            } else {
+                -nor.normalized()
+            }
+            .into_vector();
+        let fac = dot(nn, v).max(0.0) * INV_PI;
+
+        self.col.to_spectral_sample(wavelength) * fac
+    }
+
+    fn sample_pdf(&self, inc: Vector, out: Vector, nor: Normal) -> f32 {
+        let v = out.normalized();
+        let nn = if dot(nor.into_vector(), inc) <= 0.0 {
+                nor.normalized()
+            } else {
+                -nor.normalized()
+            }
+            .into_vector();
+
+        dot(nn, v).max(0.0) * INV_PI
+    }
+}

src/surface/mod.rs

@@ -5,17 +5,18 @@ use std::fmt::Debug;
 pub mod triangle_mesh;
 
 use ray::{Ray, AccelRay};
-use math::{Point, Normal, Matrix4x4};
+use math::{Point, Vector, Normal, Matrix4x4};
 use boundable::Boundable;
 
 
-#[derive(Debug, Clone)]
+#[derive(Debug, Copy, Clone)]
 pub enum SurfaceIntersection {
     Miss,
     Occlude,
     Hit {
         t: f32,
         pos: Point,
+        incoming: Vector,
         nor: Normal,
         local_space: Matrix4x4,
         uv: (f32, f32),

src/surface/triangle_mesh.rs

@@ -85,6 +85,7 @@ impl Surface for TriangleMesh {
                             isects[r.id as usize] = SurfaceIntersection::Hit {
                                 t: t,
                                 pos: wr.orig + (wr.dir * t),
+                                incoming: wr.dir,
                                 nor: cross(tri.0 - tri.1, tri.0 - tri.2).into_normal(),
                                 local_space: mat_space,
                                 uv: (tri_u, tri_v),