Prepping for adding lighting.

The part of the renderer responsible for light transport has been
split out into a LightPath struct.  Also moving over to spectral
rendering, although it's a bit silly at the moment.

src/algorithm.rs

@@ -45,3 +45,67 @@ pub fn partition<T, F>(slc: &mut [T], mut pred: F) -> usize
         }
     }
 }
+
+
+/// Partitions two slices in-place in concert based on the given unary
+/// predicate, returning the index of the first element for which the
+/// predicate evaluates false.
+///
+/// Because this runs on two slices at once, they must both be the same
+/// length.
+///
+/// The predicate takes a usize (which will receive the index of the elments
+/// being tested), a mutable reference to an element of the first slice's type,
+/// and a mutable reference to an element of the last slice's type.
+///
+/// The predicate is executed precisely once on every element in
+/// the slices, and is allowed to modify the elements.
+pub fn partition_pair<A, B, F>(slc1: &mut [A], slc2: &mut [B], mut pred: F) -> usize
+    where F: FnMut(usize, &mut A, &mut B) -> bool
+{
+    assert!(slc1.len() == slc2.len());
+
+    // This version uses raw pointers and pointer arithmetic to squeeze more
+    // performance out of the code.
+    unsafe {
+        let mut a1 = slc1.as_mut_ptr();
+        let mut a2 = slc2.as_mut_ptr();
+        let mut b1 = a1.offset(slc1.len() as isize);
+        let mut b2 = a2.offset(slc2.len() as isize);
+        let start = a1 as usize;
+
+        loop {
+            loop {
+                if a1 == b1 {
+                    return ((a1 as usize) - start) / std::mem::size_of::<A>();
+                }
+                if !pred(((a1 as usize) - start) / std::mem::size_of::<A>(),
+                         &mut *a1,
+                         &mut *a2) {
+                    break;
+                }
+                a1 = a1.offset(1);
+                a2 = a2.offset(1);
+            }
+
+            loop {
+                b1 = b1.offset(-1);
+                b2 = b2.offset(-1);
+                if a1 == b1 {
+                    return ((a1 as usize) - start) / std::mem::size_of::<A>();
+                }
+                if pred(((b1 as usize) - start) / std::mem::size_of::<A>(),
+                        &mut *b1,
+                        &mut *b2) {
+                    break;
+                }
+            }
+
+            std::ptr::swap(a1, b1);
+            std::ptr::swap(a2, b2);
+
+            a1 = a1.offset(1);
+            a2 = a2.offset(1);
+        }
+    }
+}

src/color/mod.rs

@@ -11,10 +11,14 @@ const WL_MAX: f32 = 700.0;
 const WL_RANGE: f32 = WL_MAX - WL_MIN;
 const WL_RANGE_Q: f32 = WL_RANGE / 4.0;
 
+pub fn map_0_1_to_wavelength(n: f32) -> f32 {
+    n * WL_RANGE + WL_MIN
+}
+
 pub trait Color {
     fn sample_spectrum(&self, wavelength: f32) -> f32;
 
-    fn get_spectral_sample(&self, hero_wavelength: f32) -> SpectralSample {
+    fn to_spectral_sample(&self, hero_wavelength: f32) -> SpectralSample {
         SpectralSample {
             e: [self.sample_spectrum(nth_wavelength(hero_wavelength, 0)),
                 self.sample_spectrum(nth_wavelength(hero_wavelength, 1)),
@@ -89,6 +93,14 @@ impl XYZ {
         XYZ { x: x, y: y, z: z }
     }
 
+    pub fn from_tuple(xyz: (f32, f32, f32)) -> XYZ {
+        XYZ {
+            x: xyz.0,
+            y: xyz.1,
+            z: xyz.2,
+        }
+    }
+
     pub fn from_wavelength(wavelength: f32, intensity: f32) -> XYZ {
         XYZ {
             x: x_1931(wavelength) * intensity,
@@ -97,6 +109,14 @@ impl XYZ {
         }
     }
 
+    pub fn from_spectral_sample(ss: &SpectralSample) -> XYZ {
+        let xyz0 = XYZ::from_wavelength(ss.wl_n(0), ss.e[0]);
+        let xyz1 = XYZ::from_wavelength(ss.wl_n(1), ss.e[1]);
+        let xyz2 = XYZ::from_wavelength(ss.wl_n(2), ss.e[2]);
+        let xyz3 = XYZ::from_wavelength(ss.wl_n(3), ss.e[3]);
+        (xyz0 + xyz1 + xyz2 + xyz3) * 0.75
+    }
+
     pub fn to_tuple(&self) -> (f32, f32, f32) {
         (self.x, self.y, self.z)
     }

src/light/sphere_light.rs

@@ -92,13 +92,13 @@ impl LightSource for SphereLight {
             // Calculate the final values and return everything.
             let shadow_vec = (x * sample[0]) + (y * sample[1]) + (z * sample[2]);
             let pdf = uniform_sample_cone_pdf(cos_theta_max);
-            let spectral_sample = (col * surface_area_inv as f32).get_spectral_sample(wavelength);
+            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);
             let pdf = 1.0 / (4.0 * PI_64);
-            let spectral_sample = (col * surface_area_inv as f32).get_spectral_sample(wavelength);
+            let spectral_sample = (col * surface_area_inv as f32).to_spectral_sample(wavelength);
             return (spectral_sample, shadow_vec, pdf as f32);
         }
     }
@@ -133,7 +133,7 @@ impl LightSource for SphereLight {
         let radius = lerp_slice(&self.radii, time) as f64;
         let col = lerp_slice(&self.colors, time);
         let surface_area = 4.0 * PI_64 * radius * radius;
-        (col / surface_area as f32).get_spectral_sample(wavelength)
+        (col / surface_area as f32).to_spectral_sample(wavelength)
     }
 
     fn is_delta(&self) -> bool {

src/main.rs

@@ -36,6 +36,7 @@ use std::fs::File;
 use docopt::Docopt;
 
 use ray::{Ray, AccelRay};
+use renderer::LightPath;
 use parse::{parse_scene, DataTree};
 
 // ----------------------------------------------------------------
@@ -95,6 +96,7 @@ fn main() {
 
     println!("Ray size:       {} bytes", mem::size_of::<Ray>());
     println!("AccelRay size:  {} bytes", mem::size_of::<AccelRay>());
+    println!("LightPath size: {} bytes", mem::size_of::<LightPath>());
 
     // Iterate through scenes and render them
     if let DataTree::Internal { ref children, .. } = dt {

src/renderer.rs

@@ -8,13 +8,15 @@ use std::cell::RefCell;
 use scoped_threadpool::Pool;
 use crossbeam::sync::MsQueue;
 
+use algorithm::partition_pair;
+use ray::Ray;
 use tracer::Tracer;
 use halton;
 use math::fast_logit;
 use image::Image;
 use surface;
 use scene::Scene;
-use color::{XYZ, rec709e_to_xyz};
+use color::{Color, XYZ, rec709e_to_xyz, map_0_1_to_wavelength};
 
 #[derive(Debug)]
 pub struct Renderer {
@@ -59,13 +61,13 @@ impl Renderer {
                 let ajq = &all_jobs_queued;
                 let img = ℑ
                 scope.execute(move || {
+                    let mut paths = Vec::new();
                     let mut rays = Vec::new();
-                    let mut pixel_mapping = Vec::new();
                     let mut tracer = Tracer::from_assembly(&self.scene.root);
 
                     loop {
+                        paths.clear();
                         rays.clear();
-                        pixel_mapping.clear();
 
                         // Get bucket, or exit if no more jobs left
                         let bucket: BucketJob;
@@ -85,51 +87,57 @@ impl Renderer {
                             for x in bucket.x..(bucket.x + bucket.w) {
                                 let offset = hash_u32(((x as u32) << 16) ^ (y as u32), 0);
                                 for si in 0..self.spp {
-                                    let mut ray = {
+                                    // Calculate image plane x and y coordinates
+                                    let (img_x, img_y) = {
                                         let filter_x =
-                                            fast_logit(halton::sample(3, offset + si as u32), 1.5) +
+                                            fast_logit(halton::sample(4, offset + si as u32), 1.5) +
                                             0.5;
                                         let filter_y =
-                                            fast_logit(halton::sample(4, offset + si as u32), 1.5) +
+                                            fast_logit(halton::sample(5, offset + si as u32), 1.5) +
                                             0.5;
                                         let samp_x = (filter_x + x as f32) * cmpx;
                                         let samp_y = (filter_y + y as f32) * cmpy;
-
+                                        ((samp_x - 0.5) * x_extent, (0.5 - samp_y) * y_extent)
-                                        self.scene
-                                            .camera
-                                            .generate_ray((samp_x - 0.5) * x_extent,
-                                                          (0.5 - samp_y) * y_extent,
-                                                          halton::sample(0, offset + si as u32),
-                                                          halton::sample(1, offset + si as u32),
-                                                          halton::sample(2, offset + si as u32))
                                     };
+
+                                    // Create the light path and initial ray for this sample
+                                    let (path, ray) =
+                                        LightPath::new(&self.scene,
+                                                       (x, y),
+                                                       (img_x, img_y),
+                                                       (halton::sample(0, offset + si as u32),
+                                                        halton::sample(1, offset + si as u32)),
+                                                       halton::sample(2, offset + si as u32),
+                                                       map_0_1_to_wavelength(
+                                                           halton::sample(3, offset + si as u32)
+                                                       ),
+                                                       offset + si as u32);
+                                    paths.push(path);
                                     rays.push(ray);
-                                    pixel_mapping.push((x, y))
                                 }
                             }
                         }
 
+                        // Trace the paths!
+                        let mut pi = paths.len();
+                        while pi > 0 {
                             // Test rays against scene
                             let isects = tracer.trace(&rays);
 
+                            // Determine next rays to shoot based on result
+                            pi = partition_pair(&mut paths[..pi],
+                                                &mut rays[..pi],
+                                                |i, path, ray| path.next(&isects[i], &mut *ray));
+                        }
+
                         // Calculate color based on ray hits
                         let img = img.lock().unwrap();
                         let mut img = img.borrow_mut();
-                        for (isect, co) in Iterator::zip(isects.iter(), pixel_mapping.iter()) {
+                        for path in paths.iter() {
-                            let mut col = img.get(co.0 as usize, co.1 as usize);
+                            let mut col =
-                            if let &surface::SurfaceIntersection::Hit { t: _,
+                                img.get(path.pixel_co.0 as usize, path.pixel_co.1 as usize);
-                                                                        pos: _,
+                            col += path.color / self.spp as f32;
-                                                                        nor: _,
+                            img.set(path.pixel_co.0 as usize, path.pixel_co.1 as usize, col);
-                                                                        local_space: _,
-                                                                        uv } = isect {
-                                let rgbcol = (uv.0, uv.1, (1.0 - uv.0 - uv.1).max(0.0));
-                                let xyzcol = rec709e_to_xyz(rgbcol);
-                                col += XYZ::new(xyzcol.0, xyzcol.1, xyzcol.2) / self.spp as f32;
-
-                            } else {
-                                col += XYZ::new(0.02, 0.02, 0.02) / self.spp as f32;
-                            }
-                            img.set(co.0 as usize, co.1 as usize, col);
                         }
                     }
                 });
@@ -168,6 +176,61 @@ impl Renderer {
     }
 }
 
+
+#[derive(Debug)]
+pub struct LightPath {
+    pixel_co: (u32, u32),
+    lds_offset: u32,
+    dim_offset: u32,
+    round: u32,
+    time: f32,
+    wavelength: f32,
+    color: XYZ,
+}
+
+impl LightPath {
+    fn new(scene: &Scene,
+           pixel_co: (u32, u32),
+           image_plane_co: (f32, f32),
+           lens_uv: (f32, f32),
+           time: f32,
+           wavelength: f32,
+           lds_offset: u32)
+           -> (LightPath, Ray) {
+        (LightPath {
+            pixel_co: pixel_co,
+            lds_offset: lds_offset,
+            dim_offset: 6,
+            round: 1,
+            time: time,
+            wavelength: wavelength,
+            color: XYZ::new(0.0, 0.0, 0.0),
+        },
+
+         scene.camera.generate_ray(image_plane_co.0,
+                                   image_plane_co.1,
+                                   time,
+                                   lens_uv.0,
+                                   lens_uv.1))
+    }
+
+    fn next(&mut self, isect: &surface::SurfaceIntersection, ray: &mut Ray) -> bool {
+        if let &surface::SurfaceIntersection::Hit { t: _, pos: _, nor: _, local_space: _, uv } =
+               isect {
+            let rgbcol = (uv.0, uv.1, (1.0 - uv.0 - uv.1).max(0.0));
+            let xyz = XYZ::from_tuple(rec709e_to_xyz(rgbcol));
+            self.color += XYZ::from_spectral_sample(&xyz.to_spectral_sample(self.wavelength));
+
+        } else {
+            let xyz = XYZ::new(0.02, 0.02, 0.02);
+            self.color += XYZ::from_spectral_sample(&xyz.to_spectral_sample(self.wavelength));
+        }
+
+        return false;
+    }
+}
+
+
 #[derive(Debug)]
 struct BucketJob {
     x: u32,