perf

README.md

@@ -82,6 +82,18 @@ cargo build --release
 
 It uses all cores via rayon, shows a progress bar and the region brightness table, but does not write `preview.png` or `render_metrics.png`.
 
+Beyond the Python version it adds:
+
+- **Russian roulette** - after 3 bounces, dim paths are terminated randomly (survivors are re-weighted, so the image stays unbiased). About 2.5x more image quality per second of rendering in this scene.
+- **Cheap cosine-weighted diffuse sampling** - two uniform random numbers and an orthonormal basis instead of Gaussian samples.
+- **Adaptive sampling** (optional 6th argument, off by default) - after 32 samples per pixel, keeps sampling only pixels whose 3x3 neighbourhood still has a display-space standard error above the threshold, up to `samples_per_pixel`:
+
+  ```bash
+  ./target/release/pathtracer 800 600 4096 15 cornell.png 0.03
+  ```
+
+  In this scene noise is spread fairly evenly (caustic fireflies reach almost every pixel), so adaptive sampling roughly breaks even with uniform sampling and can darken the image slightly at loose thresholds. It pays off in scenes with large easy regions.
+
 ## Scene Description
 
 The scene is the classic **Cornell Box**:

rust/Cargo.lock

@@ -201,12 +201,6 @@ version = "0.2.189"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2"
 
-[[package]]
-name = "libm"
-version = "0.2.16"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "b6d2cec3eae94f9f509c767b45932f1ada8350c4bdb85af2fcab4a3c14807981"
-
 [[package]]
 name = "miniz_oxide"
 version = "0.8.9"
@@ -244,7 +238,6 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
 dependencies = [
  "autocfg",
- "libm",
 ]
 
 [[package]]
@@ -266,7 +259,6 @@ dependencies = [
  "image",
  "indicatif",
  "rand",
- "rand_distr",
  "rayon",
 ]
 
@@ -358,16 +350,6 @@ dependencies = [
  "getrandom",
 ]
 
-[[package]]
-name = "rand_distr"
-version = "0.4.3"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31"
-dependencies = [
- "num-traits",
- "rand",
-]
-
 [[package]]
 name = "rayon"
 version = "1.12.0"

rust/Cargo.toml

@@ -6,7 +6,6 @@ description = "Cornell Box path tracer (Rust port of raytracer.py)"
 
 [dependencies]
 rand = { version = "0.8", features = ["small_rng"] }
-rand_distr = "0.4"
 rayon = "1"
 indicatif = "0.17"
 image = { version = "0.25", default-features = false, features = ["png"] }

rust/src/main.rs

@@ -67,88 +67,251 @@ fn sample_box_light<R: Rng>(scene: &Scene, light: &BoxLight, rec: &HitRecord, al
     light.emission.mul_elem(brdf) * (geometry * total_area)
 }
 
-/// `count_emission` is false after a diffuse bounce: that light was already
-/// gathered by next event estimation, so hitting it again would double count.
-fn trace<R: Rng>(scene: &Scene, ray: &Ray, depth: u32, max_depth: u32, rng: &mut R, count_emission: bool) -> Vec3 {
-    if depth >= max_depth {
-        return Vec3::ZERO;
-    }
+/// Bounces before Russian roulette may terminate a path.
+const RR_MIN_DEPTH: u32 = 3;
 
-    let Some(rec) = scene.hit(ray, EPSILON, f64::INFINITY) else {
-        let t = 0.5 * (ray.direction.y + 1.0);
-        return Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t;
-    };
+fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> Vec3 {
+    let mut radiance = Vec3::ZERO;
+    // Product of the attenuations along the path so far (divided by the RR
+    // survival probabilities), i.e. how much of any light found here reaches the camera.
+    let mut throughput = Vec3::ONE;
+    // False after a diffuse bounce: that light was already gathered by next
+    // event estimation, so hitting it again would double count.
+    let mut count_emission = true;
+    let mut ray = *primary;
 
-    if let Material::Light { emission } = rec.material {
-        return if count_emission { emission } else { Vec3::ZERO };
-    }
+    for depth in 0..max_depth {
+        let Some(rec) = scene.hit(&ray, EPSILON, f64::INFINITY) else {
+            let t = 0.5 * (ray.direction.y + 1.0);
+            let sky = Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t;
+            radiance += throughput.mul_elem(sky);
+            break;
+        };
 
-    let Some((attenuation, scattered)) = rec.material.scatter(ray, &rec, rng) else {
-        return Vec3::ZERO;
-    };
-
-    if let Material::Diffuse { albedo } = rec.material {
-        let mut direct = Vec3::ZERO;
-        for light in &scene.lights {
-            direct += sample_box_light(scene, light, &rec, albedo, rng);
+        if let Material::Light { emission } = rec.material {
+            if count_emission {
+                radiance += throughput.mul_elem(emission);
+            }
+            break;
+        }
+
+        let Some((attenuation, scattered)) = rec.material.scatter(&ray, &rec, rng) else {
+            break;
+        };
+
+        if let Material::Diffuse { albedo } = rec.material {
+            for light in &scene.lights {
+                radiance += throughput.mul_elem(sample_box_light(scene, light, &rec, albedo, rng));
+            }
+            count_emission = false;
+        } else {
+            count_emission = true;
+        }
+
+        throughput = throughput.mul_elem(attenuation);
+        ray = scattered;
+
+        // Russian roulette: end dim paths early, boosting survivors so the
+        // estimate stays unbiased. Glass/mirror paths (throughput ~1) survive.
+        if depth + 1 >= RR_MIN_DEPTH {
+            let survive = throughput.x.max(throughput.y).max(throughput.z).min(0.95);
+            if rng.gen::<f64>() >= survive {
+                break;
+            }
+            throughput = throughput / survive;
         }
-        let indirect = trace(scene, &scattered, depth + 1, max_depth, rng, false);
-        return direct + attenuation.mul_elem(indirect);
     }
 
-    attenuation.mul_elem(trace(scene, &scattered, depth + 1, max_depth, rng, true))
+    radiance
 }
 
 struct Settings {
     width: usize,
     height: usize,
+    /// Maximum samples per pixel (adaptive sampling may stop earlier).
     spp: u32,
     max_depth: u32,
     output: String,
+    /// Target standard error of a pixel in display (gamma) space, 0..1.
+    /// 0 disables adaptive sampling.
+    noise_threshold: f64,
 }
 
-fn render(scene: &Scene, s: &Settings) -> Vec<Vec3> {
-    let (width, height) = (s.width, s.height);
-    let cam = Camera::new(
-        Vec3::new(278.0, 278.0, -800.0),
-        Vec3::new(278.0, 278.0, 278.0),
-        Vec3::new(0.0, 1.0, 0.0),
-        40.0,
-        width as f64 / height as f64,
-        1078.0,
-    );
+/// Adaptive sampling: every pixel first takes this many samples...
+const ADAPTIVE_MIN_SPP: u32 = 32;
+/// ...then later passes grow each still-noisy pixel's count by this fraction.
+const ADAPTIVE_GROWTH: f64 = 0.5;
 
-    let pb = ProgressBar::new(height as u64);
-    pb.set_style(
-        ProgressStyle::with_template(
-            "{spinner:.green} Rendering [{bar:50.bright_green/blue}] {percent:>3}% {pos}/{len} rows  {elapsed_precise}  ETA {eta_precise}",
-        )
-        .unwrap()
-        .progress_chars("=> "),
-    );
+fn luminance(c: Vec3) -> f64 {
+    0.2126 * c.x + 0.7152 * c.y + 0.0722 * c.z
+}
 
-    let mut image = vec![Vec3::ZERO; width * height];
-    image
-        .par_chunks_mut(width)
-        .enumerate()
-        .for_each(|(j, row)| {
-            for (i, px) in row.iter_mut().enumerate() {
-                let mut rng = SmallRng::seed_from_u64((j * width + i) as u64 ^ 0x9E37_79B9_7F4A_7C15);
-                let mut color = Vec3::ZERO;
-                for _ in 0..s.spp {
-                    let u = (i as f64 + rng.gen::<f64>()) / width as f64;
-                    let v = ((height - 1 - j) as f64 + rng.gen::<f64>()) / height as f64;
-                    let ray = cam.get_ray(u, v);
-                    color += trace(scene, &ray, 0, s.max_depth, &mut rng, true);
+/// Running mean and variance of a pixel's luminance (Welford's algorithm).
+#[derive(Default)]
+struct PixelStats {
+    n: u32,
+    mean: f64,
+    m2: f64,
+}
+
+impl PixelStats {
+    fn add(&mut self, x: f64) {
+        self.n += 1;
+        let delta = x - self.mean;
+        self.mean += delta / self.n as f64;
+        self.m2 += delta * (x - self.mean);
+    }
+
+    /// Standard error of the pixel mean after gamma (sqrt) encoding, using
+    /// d(sqrt(x))/dx = 1 / (2 sqrt(x)) to map it from linear space.
+    fn display_error(&self) -> f64 {
+        if self.n < 2 {
+            return f64::INFINITY;
+        }
+        let variance = self.m2 / (self.n - 1) as f64;
+        let std_error = (variance / self.n as f64).sqrt();
+        std_error / (2.0 * self.mean.max(1e-4).sqrt())
+    }
+}
+
+struct Pixel {
+    sum: Vec3,
+    stats: PixelStats,
+    rng: SmallRng,
+}
+
+struct Renderer<'a> {
+    scene: &'a Scene,
+    cam: Camera,
+    s: &'a Settings,
+}
+
+impl Renderer<'_> {
+    fn sample(&self, px: &mut Pixel, i: usize, j: usize, count: u32) {
+        let (width, height) = (self.s.width, self.s.height);
+        for _ in 0..count {
+            let u = (i as f64 + px.rng.gen::<f64>()) / width as f64;
+            let v = ((height - 1 - j) as f64 + px.rng.gen::<f64>()) / height as f64;
+            let sample = trace(self.scene, &self.cam.get_ray(u, v), self.s.max_depth, &mut px.rng);
+            px.sum += sample;
+            px.stats.add(luminance(sample));
+        }
+    }
+
+    /// Per-pixel noise, taken as the maximum over each 3x3 neighbourhood. A
+    /// pixel's own estimate is unreliable while it has only seen a few samples
+    /// (a pixel that has not yet hit a rare bright path looks converged), so
+    /// a neighbour that did keeps it sampling.
+    fn neighbourhood_error(&self, pixels: &[Pixel]) -> Vec<f64> {
+        let (width, height) = (self.s.width, self.s.height);
+        let own: Vec<f64> = pixels.iter().map(|p| p.stats.display_error()).collect();
+        (0..width * height)
+            .into_par_iter()
+            .map(|idx| {
+                let (i, j) = (idx % width, idx / width);
+                let mut worst = 0.0_f64;
+                for nj in j.saturating_sub(1)..(j + 2).min(height) {
+                    for ni in i.saturating_sub(1)..(i + 2).min(width) {
+                        worst = worst.max(own[nj * width + ni]);
+                    }
                 }
-                color /= s.spp as f64;
-                *px = color.map(|c| c.max(0.0).sqrt().min(1.0));
-            }
-            pb.inc(1);
-        });
+                worst
+            })
+            .collect()
+    }
+}
+
+fn bar(len: u64, template: &str) -> ProgressBar {
+    let pb = ProgressBar::new(len);
+    pb.set_style(
+        ProgressStyle::with_template(template)
+            .unwrap()
+            .progress_chars("=> "),
+    );
+    pb
+}
+
+/// Returns the gamma-encoded image and the total number of samples taken.
+fn render(scene: &Scene, s: &Settings) -> (Vec<Vec3>, u64) {
+    let (width, height) = (s.width, s.height);
+    let r = Renderer {
+        scene,
+        cam: Camera::new(
+            Vec3::new(278.0, 278.0, -800.0),
+            Vec3::new(278.0, 278.0, 278.0),
+            Vec3::new(0.0, 1.0, 0.0),
+            40.0,
+            width as f64 / height as f64,
+            1078.0,
+        ),
+        s,
+    };
+
+    let mut pixels: Vec<Pixel> = (0..width * height)
+        .map(|idx| Pixel {
+            sum: Vec3::ZERO,
+            stats: PixelStats::default(),
+            rng: SmallRng::seed_from_u64(idx as u64 ^ 0x9E37_79B9_7F4A_7C15),
+        })
+        .collect();
+
+    // Base pass: every pixel gets the full budget, or the adaptive minimum.
+    let adaptive = s.noise_threshold > 0.0;
+    let base_spp = if adaptive { s.spp.min(ADAPTIVE_MIN_SPP) } else { s.spp };
+    let pb = bar(
+        height as u64,
+        "{spinner:.green} Rendering [{bar:50.bright_green/blue}] {percent:>3}% {pos}/{len} rows  {elapsed_precise}  ETA {eta_precise}",
+    );
+    pixels.par_chunks_mut(width).enumerate().for_each(|(j, row)| {
+        for (i, px) in row.iter_mut().enumerate() {
+            r.sample(px, i, j, base_spp);
+        }
+        pb.inc(1);
+    });
     pb.finish();
 
-    image
+    // Adaptive passes: keep sampling pixels whose neighbourhood is still noisy.
+    if adaptive {
+        let pb = bar(
+            (width * height) as u64,
+            "{spinner:.green} Refining  [{bar:50.bright_green/blue}] {percent:>3}% {pos}/{len} pixels converged  {msg}",
+        );
+        for pass in 1.. {
+            let error = r.neighbourhood_error(&pixels);
+            let active: Vec<bool> = pixels
+                .iter()
+                .zip(&error)
+                .map(|(p, &e)| p.stats.n < s.spp && e >= s.noise_threshold)
+                .collect();
+            let remaining = active.iter().filter(|&&a| a).count();
+            pb.set_position((width * height - remaining) as u64);
+            pb.set_message(format!("pass {pass}"));
+            if remaining == 0 {
+                break;
+            }
+            pixels
+                .par_chunks_mut(width)
+                .zip(active.par_chunks(width))
+                .enumerate()
+                .for_each(|(j, (row, row_active))| {
+                    for (i, (px, &a)) in row.iter_mut().zip(row_active).enumerate() {
+                        if a {
+                            let n = px.stats.n;
+                            let extra = ((n as f64 * ADAPTIVE_GROWTH).ceil() as u32).clamp(1, s.spp - n);
+                            r.sample(px, i, j, extra);
+                        }
+                    }
+                });
+        }
+        pb.finish();
+    }
+
+    let total_samples = pixels.iter().map(|p| p.stats.n as u64).sum();
+    let image = pixels
+        .iter()
+        .map(|p| (p.sum / p.stats.n as f64).map(|c| c.max(0.0).sqrt().min(1.0)))
+        .collect();
+    (image, total_samples)
 }
 
 fn save_png(image: &[Vec3], width: usize, height: usize, path: &str) -> image::ImageResult<()> {
@@ -196,22 +359,23 @@ fn print_region_table(image: &[Vec3], width: usize, height: usize) {
 fn main() {
     let args: Vec<String> = std::env::args().collect();
     let arg = |i: usize, default: &str| args.get(i).cloned().unwrap_or_else(|| default.to_string());
-    let parse = |i: usize, default: &str| -> usize {
-        arg(i, default).parse().unwrap_or_else(|_| {
-            eprintln!("usage: pathtracer [width] [height] [samples_per_pixel] [max_bounces] [output]");
+    fn parse<T: std::str::FromStr>(s: String) -> T {
+        s.parse().unwrap_or_else(|_| {
+            eprintln!("usage: pathtracer [width] [height] [samples_per_pixel] [max_bounces] [output] [noise_threshold]");
             std::process::exit(2);
         })
-    };
+    }
 
     let settings = Settings {
-        width: parse(1, "400"),
-        height: parse(2, "300"),
-        spp: parse(3, "50") as u32,
-        max_depth: parse(4, "15") as u32,
+        width: parse(arg(1, "400")),
+        height: parse(arg(2, "300")),
+        spp: parse(arg(3, "50")),
+        max_depth: parse(arg(4, "15")),
         output: arg(5, "cornell_box.png"),
+        noise_threshold: parse(arg(6, "0")),
     };
-    if settings.width == 0 || settings.height == 0 || settings.spp == 0 {
-        eprintln!("width, height and samples_per_pixel must be >= 1");
+    if settings.width == 0 || settings.height == 0 || settings.spp == 0 || settings.noise_threshold < 0.0 {
+        eprintln!("width, height and samples_per_pixel must be >= 1, noise_threshold >= 0");
         std::process::exit(2);
     }
 
@@ -219,16 +383,22 @@ fn main() {
     println!("Features: reflections, refraction/caustics, multiple bounces, NEE");
     println!();
     println!("Resolution:        {}x{}", settings.width, settings.height);
-    println!("Samples per pixel: {}", settings.spp);
-    println!("Max bounces:       {}", settings.max_depth);
+    println!("Samples per pixel: {} (max)", settings.spp);
+    println!("Max bounces:       {} (Russian roulette after {})", settings.max_depth, RR_MIN_DEPTH);
+    if settings.noise_threshold > 0.0 {
+        println!("Adaptive sampling: stop at noise < {}", settings.noise_threshold);
+    } else {
+        println!("Adaptive sampling: off");
+    }
     println!("Threads:           {}", rayon::current_num_threads());
     println!("Output:            {}", settings.output);
     println!();
 
     let scene = cornell_box();
     let start = Instant::now();
-    let image = render(&scene, &settings);
+    let (image, total_samples) = render(&scene, &settings);
     let elapsed = start.elapsed().as_secs_f64();
+    let pixels = (settings.width * settings.height) as f64;
 
     if let Err(e) = save_png(&image, settings.width, settings.height, &settings.output) {
         eprintln!("failed to write {}: {e}", settings.output);
@@ -238,9 +408,10 @@ fn main() {
     println!();
     println!("Rendered {}", settings.output);
     println!(
-        "Total time: {:.2}s | Speed: {:.2} MPix*samp/s",
+        "Total time: {:.2}s | Speed: {:.2} MPix*samp/s | Avg samples/pixel: {:.1}",
         elapsed,
-        (settings.width * settings.height) as f64 * settings.spp as f64 / elapsed / 1e6
+        total_samples as f64 / elapsed / 1e6,
+        total_samples as f64 / pixels
     );
     println!();
     print_region_table(&image, settings.width, settings.height);

rust/src/scene.rs

@@ -1,22 +1,25 @@
 use rand::Rng;
-use rand_distr::StandardNormal;
 
 use crate::vec3::{reflect, refract, schlick, Vec3};
 
 pub const EPSILON: f64 = 1e-6;
 
-pub fn random_unit_vector<R: Rng>(rng: &mut R) -> Vec3 {
-    loop {
-        let v = Vec3::new(
-            rng.sample(StandardNormal),
-            rng.sample(StandardNormal),
-            rng.sample(StandardNormal),
-        );
-        let n = v.length();
-        if n > 1e-12 {
-            return v / n;
-        }
-    }
+/// Cosine-weighted direction in the hemisphere around the unit vector `n`.
+/// Uses the branchless orthonormal basis of Duff et al. 2017 ("Building an
+/// Orthonormal Basis, Revisited"); the result is already unit length.
+pub fn cosine_hemisphere<R: Rng>(n: Vec3, rng: &mut R) -> Vec3 {
+    let r1: f64 = rng.gen();
+    let r2: f64 = rng.gen();
+    let phi = 2.0 * std::f64::consts::PI * r1;
+    let r = r2.sqrt();
+    let (x, y, z) = (r * phi.cos(), r * phi.sin(), (1.0 - r2).max(0.0).sqrt());
+
+    let sign = 1.0_f64.copysign(n.z);
+    let a = -1.0 / (sign + n.z);
+    let b = n.x * n.y * a;
+    let t = Vec3::new(1.0 + sign * n.x * n.x * a, sign * b, -sign * n.x);
+    let bt = Vec3::new(b, sign + n.y * n.y * a, -n.y);
+    t * x + bt * y + n * z
 }
 
 #[derive(Clone, Copy, Debug)]
@@ -33,6 +36,11 @@ impl Ray {
         }
     }
 
+    /// For directions already known to be unit length.
+    pub fn from_unit(origin: Vec3, direction: Vec3) -> Self {
+        Ray { origin, direction }
+    }
+
     pub fn at(&self, t: f64) -> Vec3 {
         self.origin + t * self.direction
     }
@@ -71,13 +79,10 @@ impl Material {
     pub fn scatter<R: Rng>(&self, ray_in: &Ray, rec: &HitRecord, rng: &mut R) -> Option<(Vec3, Ray)> {
         match *self {
             Material::Diffuse { albedo } => {
-                // normal + uniform unit vector is cosine-distributed, so the
-                // Lambertian BRDF * cos / pdf reduces to just the albedo.
-                let mut target = rec.normal + random_unit_vector(rng);
-                if target.length() < 1e-8 {
-                    target = rec.normal;
-                }
-                Some((albedo, Ray::new(rec.point, target)))
+                // Cosine-weighted sampling: the Lambertian BRDF * cos / pdf
+                // reduces to just the albedo.
+                let direction = cosine_hemisphere(rec.normal, rng);
+                Some((albedo, Ray::from_unit(rec.point, direction)))
             }
             Material::Mirror { albedo } => {
                 let reflected = reflect(ray_in.direction, rec.normal);