Very inefficient first-pass at global illumination

Cargo.toml

@@ -9,6 +9,7 @@ itertools = "0.9"
 obj = "0.9"
 quickcheck = "0.9"
 quickcheck_macros = "0.9"
+rand = "0.7"
 rayon = "1.3"
 sdl2 = "0.32"
 csv = "1.1.3"

src/accumulation_buffer.rs

@@ -5,16 +5,25 @@ use crate::util::{Array2D, Tile};
 #[derive(Clone, Debug)]
 pub struct AccumulationBuffer {
     colour_buffer: Array2D<ColourXyz>,
+    colour_sum_buffer: Array2D<ColourXyz>,
+    colour_bias_buffer: Array2D<ColourXyz>,
     weight_buffer: Array2D<f64>,
+    weight_bias_buffer: Array2D<f64>,
 }
 
 impl AccumulationBuffer {
     pub fn new(height: usize, width: usize) -> AccumulationBuffer {
         let colour_buffer = Array2D::new(width, height);
+        let colour_sum_buffer = Array2D::new(width, height);
+        let colour_bias_buffer = Array2D::new(width, height);
         let weight_buffer = Array2D::new(width, height);
+        let weight_bias_buffer = Array2D::new(width, height);
         AccumulationBuffer {
             colour_buffer,
+            colour_sum_buffer,
+            colour_bias_buffer,
             weight_buffer,
+            weight_bias_buffer,
         }
     }
 
@@ -34,15 +43,20 @@ impl AccumulationBuffer {
 
     pub fn update_pixel(&mut self, row: usize, column: usize, photon: &Photon, weight: f64) {
         let buffer_colour = &mut self.colour_buffer[row][column];
+        let buffer_colour_sum = &mut self.colour_sum_buffer[row][column];
+        let buffer_colour_bias = &mut self.colour_bias_buffer[row][column];
         let buffer_weight = &mut self.weight_buffer[row][column];
-
+        let buffer_weight_bias = &mut self.weight_bias_buffer[row][column];
-        *buffer_colour = blend(
+        let photon_colour = ColourXyz::from_photon(&photon);
-            buffer_colour,
+        let weight_sum_y = weight - *buffer_weight_bias;
-            *buffer_weight,
+        let weight_sum_t = *buffer_weight + weight_sum_y;
-            &ColourXyz::from_photon(&photon),
+        *buffer_weight_bias = (weight_sum_t - *buffer_weight) - weight_sum_y;
-            weight,
+        *buffer_weight = weight_sum_t;
-        );
+        let colour_sum_y = photon_colour.values * weight - buffer_colour_bias.values;
-        *buffer_weight += weight;
+        let colour_sum_t = buffer_colour_sum.values + colour_sum_y;
+        buffer_colour_bias.values = (colour_sum_t - buffer_colour_sum.values) - colour_sum_y;
+        buffer_colour_sum.values = colour_sum_t;
+        buffer_colour.values = buffer_colour_sum.values * (1.0 / *buffer_weight);
     }
 
     pub fn merge_tile(&mut self, tile: &Tile, src: &AccumulationBuffer) {

src/camera.rs

@@ -1,11 +1,8 @@
 use crate::math::Vec3;
 
 use super::accumulation_buffer::AccumulationBuffer;
-use super::colour::{
+use super::colour::Photon;
-    ColourRgbF, NamedColour, Photon, Spectrum, LONGEST_VISIBLE_WAVELENGTH,
+use super::integrators::{Integrator, SimpleRandomIntegrator};
-    SHORTEST_VISIBLE_WAVELENGTH,
-};
-use super::integrators::{DirectionalLight, Integrator, WhittedIntegrator};
 use super::raycasting::Ray;
 use super::sampler::Sampler;
 use super::scene::Scene;
@@ -67,7 +64,7 @@ impl ImageSampler {
     }
 }
 
-const RECURSION_LIMIT: u16 = 32;
+const RECURSION_LIMIT: u16 = 128;
 
 /// Render a rectangular section of the image.
 ///
@@ -101,44 +98,11 @@ pub fn partial_render_scene(
 ) -> AccumulationBuffer {
     let mut output_image_tile = AccumulationBuffer::new(tile.width(), tile.height());
     let image_sampler = ImageSampler::new(width, height, scene.camera_location);
-    let ambient_intensity = 0.0;
+    let integrator = SimpleRandomIntegrator {};
-    let directional_intensity1 = 7.0;
-    let directional_intensity2 = 3.0;
-    let directional_intensity3 = 2.0;
-    let integrator = WhittedIntegrator {
-        ambient_light: Spectrum::reflection_from_linear_rgb(
-            &(ColourRgbF::from_named(NamedColour::White) * ambient_intensity),
-        ),
-        lights: vec![
-            DirectionalLight {
-                direction: Vec3::new(1.0, 1.0, -1.0).normalize(),
-                spectrum: Spectrum::reflection_from_linear_rgb(
-                    &(ColourRgbF::from_named(NamedColour::White) * directional_intensity1),
-                ),
-            },
-            DirectionalLight {
-                direction: Vec3::new(-0.5, 2.0, -0.5).normalize(),
-                spectrum: Spectrum::reflection_from_linear_rgb(
-                    &(ColourRgbF::from_named(NamedColour::White) * directional_intensity2),
-                ),
-            },
-            DirectionalLight {
-                direction: Vec3::new(-3.0, 0.1, -0.5).normalize(),
-                spectrum: Spectrum::reflection_from_linear_rgb(
-                    &(ColourRgbF::from_named(NamedColour::White) * directional_intensity3),
-                ),
-            },
-        ],
-    };
     let sampler = Sampler { scene: &scene };
     for column in 0..tile.width() {
         for row in 0..tile.height() {
-            for wavelength_number in 0..8 {
+            let ray = image_sampler.ray_for_pixel(tile.start_row + row, tile.start_column + column);
-                let wavelength_ratio = wavelength_number as f64 / 8.0;
-                let wavelength = SHORTEST_VISIBLE_WAVELENGTH
-                    + wavelength_ratio * (LONGEST_VISIBLE_WAVELENGTH - SHORTEST_VISIBLE_WAVELENGTH);
-                let ray =
-                    image_sampler.ray_for_pixel(tile.start_row + row, tile.start_column + column);
             let hit = sampler.sample(&ray);
             let photon = match hit {
                 None => Photon {
@@ -148,15 +112,16 @@ pub fn partial_render_scene(
                 Some(intersection_info) => integrator.integrate(
                     &sampler,
                     &intersection_info,
-                        &Photon {
+                    &Photon::random_wavelength(),
-                            wavelength,
-                            intensity: 0.0,
-                        },
                     RECURSION_LIMIT,
                 ),
             };
-                output_image_tile.update_pixel(row, column, &photon, 1.0);
+            output_image_tile.update_pixel(
-            }
+                row,
+                column,
+                &photon.scale_intensity(Photon::random_wavelength_pdf(photon.wavelength)),
+                1.0,
+            );
         }
     }
     output_image_tile

src/colour/photon.rs

@@ -1,3 +1,7 @@
+use crate::colour::{LONGEST_VISIBLE_WAVELENGTH, SHORTEST_VISIBLE_WAVELENGTH};
+
+use rand::random;
+
 /// A quantum of light with a given wavelength and intensity
 #[derive(Clone, Default, Debug)]
 pub struct Photon {
@@ -11,6 +15,18 @@ pub struct Photon {
 }
 
 impl Photon {
+    pub fn random_wavelength() -> Photon {
+        Photon {
+            wavelength: SHORTEST_VISIBLE_WAVELENGTH
+                + (LONGEST_VISIBLE_WAVELENGTH - SHORTEST_VISIBLE_WAVELENGTH) * random::<f64>(),
+            intensity: 0.0,
+        }
+    }
+
+    pub fn random_wavelength_pdf(_wavelength: f64) -> f64 {
+        LONGEST_VISIBLE_WAVELENGTH - SHORTEST_VISIBLE_WAVELENGTH
+    }
+
     pub fn scale_intensity(&self, scale_factor: f64) -> Photon {
         Photon {
             wavelength: self.wavelength,

src/integrators.rs

@@ -55,8 +55,7 @@ impl Integrator for WhittedIntegrator {
                 }
             })
             .chain(
-                info.material
+                [info.material.sample(&(world_to_bsdf_space * info.retro))]
-                    .sample(&(world_to_bsdf_space * info.retro))
                     .iter()
                     .map(|direction| {
                         let world_space_direction = bsdf_to_world_space * direction;
@@ -93,3 +92,41 @@ impl Integrator for WhittedIntegrator {
             })
     }
 }
+
+pub struct SimpleRandomIntegrator {}
+
+impl Integrator for SimpleRandomIntegrator {
+    fn integrate(
+        &self,
+        sampler: &Sampler,
+        info: &IntersectionInfo,
+        photon: &Photon,
+        recursion_limit: u16,
+    ) -> Photon {
+        if recursion_limit == 0 {
+            return Photon {
+                wavelength: 0.0,
+                intensity: 0.0,
+            };
+        }
+        let world_to_bsdf_space =
+            try_change_of_basis_matrix(&info.tangent, &info.cotangent, &info.normal)
+                .expect("Normal, tangent and cotangent don't form a valid basis.");
+        let bsdf_to_world_space = world_to_bsdf_space
+            .try_inverse()
+            .expect("Expected matrix to be invertable.");
+        let w_i = info.material.sample(&(world_to_bsdf_space * info.retro));
+        let world_space_w_i = bsdf_to_world_space * w_i;
+        info.material.bsdf()(
+            &(world_to_bsdf_space * info.retro),
+            &w_i,
+            &match sampler.sample(&Ray::new(info.location, world_space_w_i).bias(0.000_000_1)) {
+                None => photon.set_intensity(world_space_w_i.y()),
+                Some(recursive_hit) => {
+                    self.integrate(&sampler, &recursive_hit, &photon, recursion_limit - 1)
+                }
+            }
+            .scale_intensity(world_space_w_i.dot(&info.normal).abs()),
+        )
+    }
+}

src/main.rs

@@ -16,7 +16,7 @@ use std::time::Duration;
 use vanrijn::accumulation_buffer::AccumulationBuffer;
 use vanrijn::colour::{ColourRgbF, NamedColour, Spectrum};
 use vanrijn::image::{ClampingToneMapper, ImageRgbU8};
-use vanrijn::materials::{LambertianMaterial, PhongMaterial, ReflectiveMaterial};
+use vanrijn::materials::LambertianMaterial;
 use vanrijn::math::Vec3;
 use vanrijn::mesh::load_obj;
 use vanrijn::partial_render_scene;
@@ -122,12 +122,12 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
     println!("Loading object...");
     let mut model_object = load_obj(
         &model_file_path,
-        Arc::new(ReflectiveMaterial {
+        Arc::new(LambertianMaterial {
             colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named(
                 NamedColour::Yellow,
             )),
             diffuse_strength: 0.05,
-            reflection_strength: 0.9,
+            //reflection_strength: 0.9,
         }),
     )?;
     println!("Building BVH...");
@@ -162,24 +162,25 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
                 Box::new(Sphere::new(
                     Vec3::new(-4.25, -0.5, 2.0),
                     1.0,
-                    Arc::new(ReflectiveMaterial {
+                    Arc::new(LambertianMaterial {
                         colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named(
                             NamedColour::Blue,
                         )),
-                        diffuse_strength: 0.01,
+                        diffuse_strength: 0.1,
-                        reflection_strength: 0.99,
+                        //                        diffuse_strength: 0.01,
+                        //                        reflection_strength: 0.99,
                     }),
                 )),
                 Box::new(Sphere::new(
                     Vec3::new(-5.0, 1.5, 1.0),
                     1.0,
-                    Arc::new(PhongMaterial {
+                    Arc::new(LambertianMaterial {
                         colour: Spectrum::reflection_from_linear_rgb(&ColourRgbF::from_named(
                             NamedColour::Red,
                         )),
                         diffuse_strength: 0.05,
-                        smoothness: 100.0,
+                        //smoothness: 100.0,
-                        specular_strength: 1.0,
+                        //specular_strength: 1.0,
                     }),
                 )),
             ]) as Box<dyn Aggregate>,
@@ -195,7 +196,8 @@ pub fn main() -> Result<(), Box<dyn std::error::Error>> {
 
     let worker_boss = std::thread::spawn(move || {
         let end_tx = tile_tx.clone();
-        TileIterator::new(image_width as usize, image_height as usize, 256)
+        TileIterator::new(image_width as usize, image_height as usize, 2048)
+            .cycle()
             .map(move |tile| (tile, tile_tx.clone()))
             .par_bridge()
             .try_for_each(|(tile, tx)| {

src/materials/mod.rs

@@ -2,6 +2,8 @@ use crate::math::Vec3;
 
 use super::colour::Photon;
 
+use rand::distributions::{Open01, OpenClosed01};
+use rand::{thread_rng, Rng};
 use std::fmt::Debug;
 
 pub mod lambertian_material;
@@ -16,7 +18,24 @@ pub use reflective_material::ReflectiveMaterial;
 pub trait Material: Debug + Sync + Send {
     fn bsdf<'a>(&'a self) -> Box<dyn Fn(&Vec3, &Vec3, &Photon) -> Photon + 'a>;
 
-    fn sample(&self, _w_o: &Vec3) -> Vec<Vec3> {
+    fn sample(&self, _w_i: &Vec3) -> Vec3 {
-        vec![]
+        let mut rng = thread_rng();
+        let mut w_o = Vec3::new(
+            2.0 * rng.sample::<f64, _>(Open01) - 1.0,
+            2.0 * rng.sample::<f64, _>(Open01) - 1.0,
+            rng.sample::<f64, _>(OpenClosed01),
+        );
+        while w_o.norm_squared() > 1.0 {
+            w_o = Vec3::new(
+                2.0 * rng.sample::<f64, _>(Open01) - 1.0,
+                2.0 * rng.sample::<f64, _>(Open01) - 1.0,
+                rng.sample::<f64, _>(OpenClosed01),
+            );
+        }
+        w_o
+    }
+
+    fn pdf(&self, _w_i: &Vec3, _w_o: &Vec3) -> f64 {
+        1.0
     }
 }

src/materials/reflective_material.rs

@@ -39,7 +39,7 @@ impl Material for ReflectiveMaterial {
         })
     }
 
-    fn sample(&self, w_o: &Vec3) -> Vec<Vec3> {
+    fn sample(&self, w_o: &Vec3) -> Vec3 {
-        vec![Vec3::new(-w_o.x(), -w_o.y(), w_o.z())]
+        Vec3::new(-w_o.x(), -w_o.y(), w_o.z())
     }
 }

src/math/vec3.rs

@@ -302,6 +302,14 @@ impl MulAssign<Mat3> for Vec3 {
     }
 }
 
+impl Mul<Vec3> for f64 {
+    type Output = Vec3;
+
+    fn mul(self, rhs: Vec3) -> Vec3 {
+        rhs * self
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;

src/util/tile_iterator.rs

@@ -15,6 +15,7 @@ impl Tile {
     }
 }
 
+#[derive(Clone)]
 pub struct TileIterator {
     tile_size: usize,
     total_height: usize,