added scenes

README.md

@@ -72,17 +72,34 @@ Multiprocessing uses one worker per CPU core by default. Pass `workers=N` to the
 
 ## Rust version
 
-A much faster port lives in [`rust/`](rust/) (renders 800x600 at 100 spp in seconds). It takes the same arguments:
+A much faster port lives in [`rust/`](rust/) (renders 800x600 at 100 spp in seconds). It takes the same arguments, plus an optional `--scene`:
 
 ```bash
 cd rust
 cargo build --release
 ./target/release/pathtracer 800 600 100 15 cornell_hq.png
+./target/release/pathtracer --scene caustics 800 600 1024 15 caustics.png
+./target/release/pathtracer --list-scenes
 ```
 
-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`.
+It uses all cores via rayon, shows a progress bar and (for the `cornell` scene) the region brightness table, but does not write `preview.png` or `render_metrics.png`. Without an output name it writes `<scene>.png`.
 
-Beyond the Python version it adds:
+### Scenes
+
+| Scene | What it shows |
+|-------|---------------|
+| `cornell` (default) | The scene from `raytracer.py` |
+| `caustics` | Clear, amber and blue glass spheres and a mirror under one light, placed so each sphere focuses a sharp caustic into its shadow. Caustics are the hardest light paths for a path tracer, so use ~1000+ spp |
+| `spheres` | Large scene: ~470 random diffuse/mirror/glass spheres under a sky (after *Ray Tracing in One Weekend*). The BVH makes it ~11x faster than testing every object |
+| `night` | The same field at night, lit by 61 glowing spheres. Each shading point picks one light, weighted by power / distance², so many lights cost no more shadow rays than one |
+| `classic` | The classic Cornell box: one-sided ceiling panel light, two boxes, colour bleeding from the walls |
+
+Scenes are defined in [`rust/src/scenes.rs`](rust/src/scenes.rs) with a small builder (`plane`, `sphere`, `cuboid`, `box_light`, `sphere_light`, `rect_light`). Every emitter must be added with a `*_light` method so light sampling knows about it.
+
+### Beyond the Python version
+
+- **BVH** - bounded objects live in a bounding volume hierarchy (median split); planes are tested separately.
+- **Many-light sampling** - box, sphere and one-sided rectangle lights, one picked per shading point in proportion to its estimated contribution.
 
 - **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.
@@ -108,7 +125,7 @@ The scene is the classic **Cornell Box**:
 
 - Position: (278, 278, -800)
 - Look-at: (278, 278, 278) (box center)
-- FOV: 40°
+- FOV: 40° vertical (earlier versions effectively used half of this, zooming the image in 2x)
 - Focus distance: 1078
 
 ## Materials

raytracer.py

@@ -283,8 +283,8 @@ class Camera:
         w = normalize(self.origin - np.asarray(lookat, dtype=np.float64))
         u = normalize(np.cross(vup, w))
         v = np.cross(w, u)
-        self.horizontal = focus_dist * half_width * u
-        self.vertical = focus_dist * half_height * v
+        self.horizontal = 2.0 * focus_dist * half_width * u
+        self.vertical = 2.0 * focus_dist * half_height * v
         self.lower_left_corner = self.origin - self.horizontal / 2.0 - self.vertical / 2.0 - focus_dist * w
         self.aperture = aperture
 

rust/src/bvh.rs

@@ -0,0 +1,185 @@
+//! Bounding volume hierarchy over the scene's bounded objects, so a ray only
+//! tests the objects whose boxes it passes through instead of all of them.
+
+use crate::scene::Ray;
+use crate::vec3::Vec3;
+
+#[derive(Clone, Copy, Debug)]
+pub struct Aabb {
+    pub min: Vec3,
+    pub max: Vec3,
+}
+
+impl Aabb {
+    pub const EMPTY: Aabb = Aabb {
+        min: Vec3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY),
+        max: Vec3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY),
+    };
+
+    pub fn union(self, o: Aabb) -> Aabb {
+        Aabb {
+            min: Vec3::new(self.min.x.min(o.min.x), self.min.y.min(o.min.y), self.min.z.min(o.min.z)),
+            max: Vec3::new(self.max.x.max(o.max.x), self.max.y.max(o.max.y), self.max.z.max(o.max.z)),
+        }
+    }
+
+    pub fn centroid(&self) -> Vec3 {
+        (self.min + self.max) * 0.5
+    }
+
+    /// Slab test. `inv_dir` components may be infinite for axis-parallel rays;
+    /// the resulting NaNs are discarded by f64::min/max.
+    fn hit(&self, origin: Vec3, inv_dir: Vec3, mut t_min: f64, mut t_max: f64) -> bool {
+        for axis in 0..3 {
+            let mut t0 = (self.min[axis] - origin[axis]) * inv_dir[axis];
+            let mut t1 = (self.max[axis] - origin[axis]) * inv_dir[axis];
+            if inv_dir[axis] < 0.0 {
+                std::mem::swap(&mut t0, &mut t1);
+            }
+            t_min = t_min.max(t0);
+            t_max = t_max.min(t1);
+            if t_max < t_min {
+                return false;
+            }
+        }
+        true
+    }
+}
+
+const LEAF_SIZE: usize = 4;
+
+/// Flattened node: a leaf covers `prims[start..start + count]`; an interior
+/// node (count == 0) has its left child at the next index and its right child at `right`.
+struct Node {
+    bounds: Aabb,
+    start: u32,
+    count: u32,
+    right: u32,
+    axis: u8,
+}
+
+pub struct Bvh {
+    nodes: Vec<Node>,
+    /// Object indices, reordered so each leaf's objects are contiguous.
+    prims: Vec<u32>,
+}
+
+impl Bvh {
+    pub fn build(bounds: &[Aabb]) -> Bvh {
+        let mut bvh = Bvh {
+            nodes: Vec::with_capacity(2 * bounds.len()),
+            prims: (0..bounds.len() as u32).collect(),
+        };
+        if !bounds.is_empty() {
+            let centroids: Vec<Vec3> = bounds.iter().map(Aabb::centroid).collect();
+            bvh.build_range(0, bounds.len(), bounds, &centroids);
+        }
+        bvh
+    }
+
+    /// Median split along the axis where the object centroids spread widest.
+    fn build_range(&mut self, start: usize, end: usize, bounds: &[Aabb], centroids: &[Vec3]) -> usize {
+        let prims = &self.prims[start..end];
+        let node_bounds = prims.iter().fold(Aabb::EMPTY, |b, &p| b.union(bounds[p as usize]));
+        let centroid_bounds = prims.iter().fold(Aabb::EMPTY, |b, &p| {
+            let c = centroids[p as usize];
+            b.union(Aabb { min: c, max: c })
+        });
+        let extent = centroid_bounds.max - centroid_bounds.min;
+        let axis = if extent.x >= extent.y && extent.x >= extent.z {
+            0
+        } else if extent.y >= extent.z {
+            1
+        } else {
+            2
+        };
+
+        let idx = self.nodes.len();
+        self.nodes.push(Node {
+            bounds: node_bounds,
+            start: start as u32,
+            count: (end - start) as u32,
+            right: 0,
+            axis: axis as u8,
+        });
+        if end - start <= LEAF_SIZE || extent[axis] <= 0.0 {
+            return idx;
+        }
+
+        let half = (end - start) / 2;
+        self.prims[start..end].select_nth_unstable_by(half, |&a, &b| {
+            centroids[a as usize][axis].total_cmp(&centroids[b as usize][axis])
+        });
+        self.build_range(start, start + half, bounds, centroids);
+        let right = self.build_range(start + half, end, bounds, centroids);
+        self.nodes[idx].count = 0;
+        self.nodes[idx].right = right as u32;
+        idx
+    }
+
+    /// Calls `visit` for each object in a leaf whose box the ray passes through
+    /// within (t_min, t_max). `visit` returns Some(t) to shrink t_max after a
+    /// hit, or Err(()) to stop traversal entirely.
+    fn traverse(
+        &self,
+        r: &Ray,
+        t_min: f64,
+        mut t_max: f64,
+        mut visit: impl FnMut(usize, f64) -> Result<Option<f64>, ()>,
+    ) {
+        if self.nodes.is_empty() {
+            return;
+        }
+        let inv_dir = Vec3::new(1.0 / r.direction.x, 1.0 / r.direction.y, 1.0 / r.direction.z);
+        let mut stack = [0u32; 64];
+        let mut sp = 1;
+        while sp > 0 {
+            sp -= 1;
+            let node = &self.nodes[stack[sp] as usize];
+            if !node.bounds.hit(r.origin, inv_dir, t_min, t_max) {
+                continue;
+            }
+            if node.count > 0 {
+                for &p in &self.prims[node.start as usize..(node.start + node.count) as usize] {
+                    match visit(p as usize, t_max) {
+                        Ok(Some(t)) => t_max = t,
+                        Ok(None) => {}
+                        Err(()) => return,
+                    }
+                }
+            } else {
+                // Push the far child first so the near one is visited first,
+                // which shrinks t_max sooner.
+                let left = stack[sp] + 1;
+                let (near, far) = if r.direction[node.axis as usize] < 0.0 {
+                    (node.right, left)
+                } else {
+                    (left, node.right)
+                };
+                stack[sp] = far;
+                stack[sp + 1] = near;
+                sp += 2;
+            }
+        }
+    }
+
+    /// Closest hit: `hit(object, t_max)` returns the hit distance if the object
+    /// is hit closer than t_max.
+    pub fn closest(&self, r: &Ray, t_min: f64, t_max: f64, mut hit: impl FnMut(usize, f64) -> Option<f64>) {
+        self.traverse(r, t_min, t_max, |p, t_max| Ok(hit(p, t_max)));
+    }
+
+    /// True as soon as `hit(object)` reports any hit.
+    pub fn any(&self, r: &Ray, t_min: f64, t_max: f64, mut hit: impl FnMut(usize) -> bool) -> bool {
+        let mut found = false;
+        self.traverse(r, t_min, t_max, |p, _| {
+            if hit(p) {
+                found = true;
+                Err(())
+            } else {
+                Ok(None)
+            }
+        });
+        found
+    }
+}

rust/src/camera.rs

@@ -15,8 +15,8 @@ impl Camera {
         let w = (lookfrom - lookat).normalize();
         let u = vup.cross(w).normalize();
         let v = w.cross(u);
-        let horizontal = focus_dist * half_width * u;
-        let vertical = focus_dist * half_height * v;
+        let horizontal = 2.0 * focus_dist * half_width * u;
+        let vertical = 2.0 * focus_dist * half_height * v;
         let lower_left_corner = lookfrom - horizontal / 2.0 - vertical / 2.0 - focus_dist * w;
         Camera {
             origin: lookfrom,

rust/src/main.rs

@@ -1,5 +1,7 @@
+mod bvh;
 mod camera;
 mod scene;
+mod scenes;
 mod vec3;
 
 use std::f64::consts::PI;
@@ -11,43 +13,20 @@ use rand::{Rng, SeedableRng};
 use rayon::prelude::*;
 
 use camera::Camera;
-use scene::{cornell_box, BoxLight, HitRecord, Material, Ray, Scene, EPSILON};
+use scene::{HitRecord, Light, Material, Ray, Scene, EPSILON};
 use vec3::Vec3;
 
-const SHADOW_EPSILON: f64 = 1e-3;
+/// Shadow rays stop this fraction short of the light, so they don't hit the
+/// light's own surface. Relative, because scenes range from ~1 to ~1000 units.
+const SHADOW_RELATIVE_EPSILON: f64 = 1e-4;
 
-/// Direct light from one uniformly sampled point on the box light's surface.
-fn sample_box_light<R: Rng>(scene: &Scene, light: &BoxLight, rec: &HitRecord, albedo: Vec3, rng: &mut R) -> Vec3 {
-    let (lo, hi) = (light.min, light.max);
-    let size = hi - lo;
-    let face_areas = [size.y * size.z, size.x * size.z, size.x * size.y];
-    let half_area: f64 = face_areas.iter().sum();
-    let total_area = 2.0 * half_area;
-
-    // Pick a face with probability proportional to its area, then a uniform
-    // point on it: overall pdf is 1 / total_area.
-    let u = rng.gen::<f64>() * half_area;
-    let axis = if u < face_areas[0] {
-        0
-    } else if u < face_areas[0] + face_areas[1] {
-        1
-    } else {
-        2
-    };
-    let positive_side = rng.gen::<f64>() < 0.5;
-    let mut p = [
-        lo.x + size.x * rng.gen::<f64>(),
-        lo.y + size.y * rng.gen::<f64>(),
-        lo.z + size.z * rng.gen::<f64>(),
-    ];
-    p[axis] = if positive_side { hi[axis] } else { lo[axis] };
-    let mut n = [0.0; 3];
-    n[axis] = if positive_side { 1.0 } else { -1.0 };
-    let (lp, light_normal) = (Vec3::new(p[0], p[1], p[2]), Vec3::new(n[0], n[1], n[2]));
+/// Direct light from one uniformly sampled point on `light`'s surface.
+fn sample_direct<R: Rng>(scene: &Scene, light: &Light, rec: &HitRecord, albedo: Vec3, rng: &mut R) -> Vec3 {
+    let (lp, light_normal) = light.sample(rng);
 
     let to_light = lp - rec.point;
     let dist_sq = to_light.dot(to_light);
-    if dist_sq < 1e-8 {
+    if dist_sq < 1e-12 {
         return Vec3::ZERO;
     }
     let dist = dist_sq.sqrt();
@@ -58,13 +37,14 @@ fn sample_box_light<R: Rng>(scene: &Scene, light: &BoxLight, rec: &HitRecord, al
         return Vec3::ZERO;
     }
 
-    if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist - SHADOW_EPSILON) {
+    if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist * (1.0 - SHADOW_RELATIVE_EPSILON)) {
         return Vec3::ZERO;
     }
 
+    // Uniform area sampling has pdf 1 / area.
     let brdf = albedo / PI;
     let geometry = cos_at_rec * cos_at_light / dist_sq;
-    light.emission.mul_elem(brdf) * (geometry * total_area)
+    light.emission.mul_elem(brdf) * (geometry * light.area())
 }
 
 /// Bounces before Russian roulette may terminate a path.
@@ -82,14 +62,13 @@ fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> V
 
     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);
+            radiance += throughput.mul_elem(scene.background.radiance(ray.direction));
             break;
         };
 
         if let Material::Light { emission } = rec.material {
-            if count_emission {
+            // Lights emit from their front side only (matters for one-sided rects).
+            if count_emission && rec.front_face {
                 radiance += throughput.mul_elem(emission);
             }
             break;
@@ -100,8 +79,9 @@ fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> V
         };
 
         if let Material::Diffuse { albedo } = rec.material {
-            for light in &scene.lights {
-                radiance += throughput.mul_elem(sample_box_light(scene, light, &rec, albedo, rng));
+            if let Some((light, pick_prob)) = scene.pick_light(rec.point, rng) {
+                let direct = sample_direct(scene, light, &rec, albedo, rng) / pick_prob;
+                radiance += throughput.mul_elem(direct);
             }
             count_emission = false;
         } else {
@@ -142,10 +122,6 @@ 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;
 
-fn luminance(c: Vec3) -> f64 {
-    0.2126 * c.x + 0.7152 * c.y + 0.0722 * c.z
-}
-
 /// Running mean and variance of a pixel's luminance (Welford's algorithm).
 #[derive(Default)]
 struct PixelStats {
@@ -194,7 +170,7 @@ impl Renderer<'_> {
             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));
+            px.stats.add(sample.luminance());
         }
     }
 
@@ -234,15 +210,16 @@ fn bar(len: u64, template: &str) -> ProgressBar {
 /// 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 view = scene.view;
     let r = Renderer {
         scene,
         cam: Camera::new(
-            Vec3::new(278.0, 278.0, -800.0),
-            Vec3::new(278.0, 278.0, 278.0),
+            view.lookfrom,
+            view.lookat,
             Vec3::new(0.0, 1.0, 0.0),
-            40.0,
+            view.vfov,
             width as f64 / height as f64,
-            1078.0,
+            (view.lookfrom - view.lookat).length(),
         ),
         s,
     };
@@ -356,32 +333,73 @@ fn print_region_table(image: &[Vec3], width: usize, height: usize) {
     }
 }
 
+const USAGE: &str = "usage: pathtracer [--scene NAME] [--list-scenes] \
+[width] [height] [samples_per_pixel] [max_bounces] [output] [noise_threshold]";
+
+fn list_scenes() {
+    println!("Scenes:");
+    for s in scenes::SCENES {
+        println!("  {:<10} {}", s.name, s.description);
+    }
+}
+
 fn main() {
-    let args: Vec<String> = std::env::args().collect();
+    // Pull out the --scene / --list-scenes options; the rest are positional.
+    let mut scene_name = String::from("cornell");
+    let mut args: Vec<String> = Vec::new();
+    let mut raw = std::env::args().skip(1);
+    while let Some(a) = raw.next() {
+        if a == "--list-scenes" {
+            list_scenes();
+            return;
+        } else if a == "--scene" {
+            scene_name = raw.next().unwrap_or_else(|| {
+                eprintln!("{USAGE}");
+                std::process::exit(2);
+            });
+        } else if let Some(name) = a.strip_prefix("--scene=") {
+            scene_name = name.to_string();
+        } else if a == "--help" || a == "-h" {
+            println!("{USAGE}");
+            list_scenes();
+            return;
+        } else {
+            args.push(a);
+        }
+    }
+    let Some(scene_info) = scenes::find(&scene_name) else {
+        eprintln!("unknown scene '{scene_name}'");
+        list_scenes();
+        std::process::exit(2);
+    };
+
     let arg = |i: usize, default: &str| args.get(i).cloned().unwrap_or_else(|| default.to_string());
     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]");
+            eprintln!("{USAGE}");
             std::process::exit(2);
         })
     }
 
     let settings = Settings {
-        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")),
+        width: parse(arg(0, "400")),
+        height: parse(arg(1, "300")),
+        spp: parse(arg(2, "50")),
+        max_depth: parse(arg(3, "15")),
+        output: arg(4, &format!("{scene_name}.png")),
+        noise_threshold: parse(arg(5, "0")),
     };
     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);
     }
 
-    println!("=== Cornell Box Path Tracer (Rust) ===");
-    println!("Features: reflections, refraction/caustics, multiple bounces, NEE");
+    let scene = (scene_info.build)();
+
+    println!("=== Path Tracer (Rust) ===");
+    println!("Features: reflections, refraction/caustics, multiple bounces, NEE, BVH");
     println!();
+    println!("Scene:             {} ({} objects, {} lights)", scene_name, scene.object_count(), scene.light_count());
     println!("Resolution:        {}x{}", settings.width, settings.height);
     println!("Samples per pixel: {} (max)", settings.spp);
     println!("Max bounces:       {} (Russian roulette after {})", settings.max_depth, RR_MIN_DEPTH);
@@ -394,7 +412,6 @@ fn main() {
     println!("Output:            {}", settings.output);
     println!();
 
-    let scene = cornell_box();
     let start = Instant::now();
     let (image, total_samples) = render(&scene, &settings);
     let elapsed = start.elapsed().as_secs_f64();
@@ -413,6 +430,9 @@ fn main() {
         total_samples as f64 / elapsed / 1e6,
         total_samples as f64 / pixels
     );
-    println!();
-    print_region_table(&image, settings.width, settings.height);
+    // The regions (left wall, floor, ...) only mean something for the Cornell box.
+    if scene_name == "cornell" {
+        println!();
+        print_region_table(&image, settings.width, settings.height);
+    }
 }

rust/src/scene.rs

@@ -1,5 +1,6 @@
 use rand::Rng;
 
+use crate::bvh::{Aabb, Bvh};
 use crate::vec3::{reflect, refract, schlick, Vec3};
 
 pub const EPSILON: f64 = 1e-6;
@@ -114,6 +115,10 @@ pub enum Shape {
     Sphere { center: Vec3, radius: f64 },
     Plane { point: Vec3, normal: Vec3 },
     Box { min: Vec3, max: Vec3 },
+    /// Axis-aligned rectangle spanning min..max, flat along the axis of
+    /// `normal` (which must be ±x, ±y or ±z). As a light it emits only on the
+    /// `normal` side.
+    Rect { min: Vec3, max: Vec3, normal: Vec3 },
 }
 
 impl Shape {
@@ -183,6 +188,38 @@ impl Shape {
                 n[hit_axis] = sign;
                 Some((t_lo, Vec3::new(n[0], n[1], n[2])))
             }
+            Shape::Rect { min, max, normal } => {
+                let denom = r.direction.dot(normal);
+                if denom.abs() < 1e-12 {
+                    return None;
+                }
+                let t = (min - r.origin).dot(normal) / denom;
+                if t < t_min || t > t_max {
+                    return None;
+                }
+                let p = r.at(t);
+                let tol = 1e-9;
+                let inside = (0..3).all(|a| p[a] >= min[a] - tol && p[a] <= max[a] + tol);
+                inside.then_some((t, normal))
+            }
+        }
+    }
+}
+
+impl Shape {
+    /// None for unbounded shapes (planes), which are tested outside the BVH.
+    fn bounds(&self) -> Option<Aabb> {
+        match *self {
+            Shape::Sphere { center, radius } => Some(Aabb {
+                min: center - Vec3::splat(radius),
+                max: center + Vec3::splat(radius),
+            }),
+            Shape::Box { min, max } => Some(Aabb { min, max }),
+            Shape::Rect { min, max, .. } => Some(Aabb {
+                min: min - Vec3::splat(1e-9),
+                max: max + Vec3::splat(1e-9),
+            }),
+            Shape::Plane { .. } => None,
         }
     }
 }
@@ -193,79 +230,278 @@ pub struct Object {
     pub material: Material,
 }
 
-/// An axis-aligned box light, sampled explicitly for next event estimation.
+/// An emitter that next event estimation can sample: a box, sphere or rectangle.
 #[derive(Clone, Copy, Debug)]
-pub struct BoxLight {
-    pub min: Vec3,
-    pub max: Vec3,
+pub struct Light {
+    shape: Shape,
     pub emission: Vec3,
+    /// Total emitted power (up to a constant): luminance * area.
+    power: f64,
+    center: Vec3,
+    /// Half the bounding box diagonal: how close a point can be to the centre
+    /// while still outside the light.
+    extent: f64,
+}
+
+impl Light {
+    fn new(shape: Shape, emission: Vec3) -> Self {
+        let bounds = shape.bounds().expect("planes cannot be lights");
+        let mut light = Light {
+            shape,
+            emission,
+            power: 0.0,
+            center: bounds.centroid(),
+            extent: 0.5 * (bounds.max - bounds.min).length(),
+        };
+        light.power = emission.luminance() * light.area();
+        light
+    }
+
+    pub fn area(&self) -> f64 {
+        match self.shape {
+            Shape::Box { min, max } => {
+                let d = max - min;
+                2.0 * (d.x * d.y + d.y * d.z + d.x * d.z)
+            }
+            Shape::Sphere { radius, .. } => 4.0 * std::f64::consts::PI * radius * radius,
+            Shape::Rect { min, max, normal } => {
+                let d = max - min;
+                if normal.x != 0.0 {
+                    d.y * d.z
+                } else if normal.y != 0.0 {
+                    d.x * d.z
+                } else {
+                    d.x * d.y
+                }
+            }
+            Shape::Plane { .. } => unreachable!("planes cannot be lights"),
+        }
+    }
+
+    /// Rough estimate of how much this light contributes at `point`: its
+    /// power over the squared distance, clamped so points near or inside the
+    /// light's bounds don't blow up.
+    fn importance(&self, point: Vec3) -> f64 {
+        let d = point - self.center;
+        self.power / d.dot(d).max(self.extent * self.extent)
+    }
+
+    /// A uniformly distributed point on the surface (pdf 1 / area) and the
+    /// outward normal there.
+    pub fn sample<R: Rng>(&self, rng: &mut R) -> (Vec3, Vec3) {
+        match self.shape {
+            Shape::Box { min, max } => {
+                // Pick a face with probability proportional to its area, then
+                // a uniform point on it.
+                let size = max - min;
+                let face_areas = [size.y * size.z, size.x * size.z, size.x * size.y];
+                let u = rng.gen::<f64>() * (face_areas[0] + face_areas[1] + face_areas[2]);
+                let axis = if u < face_areas[0] {
+                    0
+                } else if u < face_areas[0] + face_areas[1] {
+                    1
+                } else {
+                    2
+                };
+                let positive_side = rng.gen::<f64>() < 0.5;
+                let mut p = [
+                    min.x + size.x * rng.gen::<f64>(),
+                    min.y + size.y * rng.gen::<f64>(),
+                    min.z + size.z * rng.gen::<f64>(),
+                ];
+                p[axis] = if positive_side { max[axis] } else { min[axis] };
+                let mut n = [0.0; 3];
+                n[axis] = if positive_side { 1.0 } else { -1.0 };
+                (Vec3::new(p[0], p[1], p[2]), Vec3::new(n[0], n[1], n[2]))
+            }
+            Shape::Sphere { center, radius } => {
+                let z = 1.0 - 2.0 * rng.gen::<f64>();
+                let r = (1.0 - z * z).max(0.0).sqrt();
+                let phi = 2.0 * std::f64::consts::PI * rng.gen::<f64>();
+                let n = Vec3::new(r * phi.cos(), r * phi.sin(), z);
+                (center + radius * n, n)
+            }
+            Shape::Rect { min, max, normal } => {
+                let size = max - min;
+                let p = Vec3::new(
+                    min.x + size.x * rng.gen::<f64>(),
+                    min.y + size.y * rng.gen::<f64>(),
+                    min.z + size.z * rng.gen::<f64>(),
+                );
+                (p, normal)
+            }
+            Shape::Plane { .. } => unreachable!("planes cannot be lights"),
+        }
+    }
+}
+
+/// What a ray that leaves the scene sees.
+#[derive(Clone, Copy, Debug)]
+pub enum Background {
+    /// White-to-blue gradient, as in raytracer.py.
+    Sky,
+    Color(Vec3),
+}
+
+impl Background {
+    pub fn radiance(&self, direction: Vec3) -> Vec3 {
+        match *self {
+            Background::Sky => {
+                let t = 0.5 * (direction.y + 1.0);
+                Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t
+            }
+            Background::Color(c) => c,
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct View {
+    pub lookfrom: Vec3,
+    pub lookat: Vec3,
+    pub vfov: f64,
+}
+
+pub struct SceneBuilder {
+    objects: Vec<Object>,
+    planes: Vec<Object>,
+    lights: Vec<Light>,
+    background: Background,
+    view: View,
+}
+
+impl SceneBuilder {
+    pub fn new(view: View, background: Background) -> Self {
+        SceneBuilder {
+            objects: Vec::new(),
+            planes: Vec::new(),
+            lights: Vec::new(),
+            background,
+            view,
+        }
+    }
+
+    pub fn add(&mut self, shape: Shape, material: Material) {
+        let obj = Object { shape, material };
+        match shape {
+            Shape::Plane { .. } => self.planes.push(obj),
+            _ => self.objects.push(obj),
+        }
+    }
+
+    pub fn plane(&mut self, point: Vec3, normal: Vec3, material: Material) {
+        self.add(Shape::Plane { point, normal: normal.normalize() }, material);
+    }
+
+    pub fn sphere(&mut self, center: Vec3, radius: f64, material: Material) {
+        self.add(Shape::Sphere { center, radius }, material);
+    }
+
+    pub fn cuboid(&mut self, min: Vec3, max: Vec3, material: Material) {
+        self.add(Shape::Box { min, max }, material);
+    }
+
+    /// Every emitter must be added through a *_light method: after a diffuse
+    /// bounce, trace() only counts light found by sampling these.
+    pub fn box_light(&mut self, min: Vec3, max: Vec3, emission: Vec3) {
+        self.add_light(Shape::Box { min, max }, emission);
+    }
+
+    pub fn sphere_light(&mut self, center: Vec3, radius: f64, emission: Vec3) {
+        self.add_light(Shape::Sphere { center, radius }, emission);
+    }
+
+    /// One-sided rectangular light emitting towards `normal` (±x, ±y or ±z).
+    pub fn rect_light(&mut self, min: Vec3, max: Vec3, normal: Vec3, emission: Vec3) {
+        self.add_light(Shape::Rect { min, max, normal }, emission);
+    }
+
+    fn add_light(&mut self, shape: Shape, emission: Vec3) {
+        self.lights.push(Light::new(shape, emission));
+        self.add(shape, Material::Light { emission });
+    }
+
+    pub fn build(self) -> Scene {
+        let bounds: Vec<Aabb> = self.objects.iter().map(|o| o.shape.bounds().unwrap()).collect();
+        Scene {
+            bvh: Bvh::build(&bounds),
+            objects: self.objects,
+            planes: self.planes,
+            lights: self.lights,
+            background: self.background,
+            view: self.view,
+        }
+    }
 }
 
-#[derive(Default)]
 pub struct Scene {
-    pub objects: Vec<Object>,
-    pub lights: Vec<BoxLight>,
+    objects: Vec<Object>,
+    planes: Vec<Object>,
+    bvh: Bvh,
+    lights: Vec<Light>,
+    pub background: Background,
+    pub view: View,
 }
 
 impl Scene {
-    pub fn add(&mut self, shape: Shape, material: Material) {
-        self.objects.push(Object { shape, material });
+    pub fn object_count(&self) -> usize {
+        self.objects.len() + self.planes.len()
     }
 
-    pub fn add_box_light(&mut self, min: Vec3, max: Vec3, emission: Vec3) {
-        self.lights.push(BoxLight { min, max, emission });
-        self.add(Shape::Box { min, max }, Material::Light { emission });
+    pub fn light_count(&self) -> usize {
+        self.lights.len()
+    }
+
+    /// Picks one light for shading `point`, with probability proportional to
+    /// its estimated contribution there (power / distance^2), so a scene with
+    /// many lights costs one shadow ray per shading point and mostly samples
+    /// nearby lights. Returns the light and the probability of picking it.
+    /// This is O(lights) per call; thousands of lights would want a light tree.
+    pub fn pick_light<R: Rng>(&self, point: Vec3, rng: &mut R) -> Option<(&Light, f64)> {
+        match self.lights.len() {
+            0 => None,
+            1 => Some((&self.lights[0], 1.0)),
+            _ => {
+                let total: f64 = self.lights.iter().map(|l| l.importance(point)).sum();
+                let u = rng.gen::<f64>() * total;
+                let mut acc = 0.0;
+                for light in &self.lights {
+                    let w = light.importance(point);
+                    acc += w;
+                    if u < acc {
+                        return Some((light, w / total));
+                    }
+                }
+                // Rounding left u just past the last bucket.
+                let last = self.lights.last().unwrap();
+                Some((last, last.importance(point) / total))
+            }
+        }
     }
 
     pub fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
         let mut closest: Option<(f64, Vec3, Material)> = None;
         let mut min_t = t_max;
-        for obj in &self.objects {
+        for obj in &self.planes {
             if let Some((t, outward)) = obj.shape.hit(r, t_min, min_t) {
-                if t < min_t {
-                    min_t = t;
-                    closest = Some((t, outward, obj.material));
-                }
+                min_t = t;
+                closest = Some((t, outward, obj.material));
             }
         }
+        self.bvh.closest(r, t_min, min_t, |i, t_max| {
+            let obj = &self.objects[i];
+            let (t, outward) = obj.shape.hit(r, t_min, t_max)?;
+            closest = Some((t, outward, obj.material));
+            Some(t)
+        });
         closest.map(|(t, outward, material)| HitRecord::new(r, t, outward, material))
     }
 
     pub fn blocked(&self, origin: Vec3, direction: Vec3, t_max: f64) -> bool {
-        let r = Ray::new(origin, direction);
-        self.objects
-            .iter()
-            .any(|obj| obj.shape.hit(&r, EPSILON, t_max).is_some())
+        let r = Ray::from_unit(origin, direction);
+        self.planes.iter().any(|obj| obj.shape.hit(&r, EPSILON, t_max).is_some())
+            || self
+                .bvh
+                .any(&r, EPSILON, t_max, |i| self.objects[i].shape.hit(&r, EPSILON, t_max).is_some())
     }
 }
-
-pub fn cornell_box() -> Scene {
-    let mut scene = Scene::default();
-
-    let white = Material::Diffuse { albedo: Vec3::new(0.73, 0.73, 0.73) };
-    let gray = Material::Diffuse { albedo: Vec3::new(0.5, 0.5, 0.5) };
-    let red = Material::Diffuse { albedo: Vec3::new(0.65, 0.05, 0.05) };
-    let green = Material::Diffuse { albedo: Vec3::new(0.12, 0.45, 0.15) };
-    let mirror = Material::Mirror { albedo: Vec3::new(0.95, 0.95, 0.95) };
-    let glass = Material::Dielectric { ref_idx: 1.5, albedo: Vec3::ONE };
-
-    let plane = |p: Vec3, n: Vec3| Shape::Plane { point: p, normal: n.normalize() };
-    scene.add(plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0)), gray);
-    scene.add(plane(Vec3::new(0.0, 500.0, 0.0), Vec3::new(0.0, -1.0, 0.0)), white);
-    scene.add(plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0)), white);
-    scene.add(plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0)), red);
-    scene.add(plane(Vec3::new(500.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0)), green);
-
-    scene.add_box_light(
-        Vec3::new(0.0, 180.0, 180.0),
-        Vec3::new(25.0, 220.0, 320.0),
-        Vec3::splat(50.0),
-    );
-
-    let sphere = |c: Vec3, r: f64| Shape::Sphere { center: c, radius: r };
-    scene.add(sphere(Vec3::new(250.0, 180.0, 250.0), 80.0), glass);
-    scene.add(sphere(Vec3::new(150.0, 180.0, 380.0), 80.0), mirror);
-    scene.add(sphere(Vec3::new(350.0, 160.0, 300.0), 60.0), red);
-
-    scene
-}

rust/src/scenes.rs

@@ -0,0 +1,220 @@
+//! The built-in scenes, selectable with `--scene NAME`.
+
+use rand::rngs::SmallRng;
+use rand::{Rng, SeedableRng};
+
+use crate::scene::{Background, Material, Scene, SceneBuilder, View};
+use crate::vec3::Vec3;
+
+pub struct SceneInfo {
+    pub name: &'static str,
+    pub description: &'static str,
+    pub build: fn() -> Scene,
+}
+
+pub const SCENES: &[SceneInfo] = &[
+    SceneInfo {
+        name: "cornell",
+        description: "Cornell box from raytracer.py: glass, mirror and diffuse spheres, wall light",
+        build: cornell,
+    },
+    SceneInfo {
+        name: "caustics",
+        description: "Clear, amber and blue glass plus a mirror under one light, focusing caustics onto the floor",
+        build: caustics,
+    },
+    SceneInfo {
+        name: "spheres",
+        description: "Large scene: ~480 random spheres under a sky (exercises the BVH)",
+        build: spheres,
+    },
+    SceneInfo {
+        name: "night",
+        description: "The same sphere field at night, lit by dozens of glowing spheres (many-light sampling)",
+        build: night,
+    },
+    SceneInfo {
+        name: "classic",
+        description: "The classic Cornell box: ceiling light, two boxes, colour bleeding",
+        build: classic,
+    },
+];
+
+pub fn find(name: &str) -> Option<&'static SceneInfo> {
+    SCENES.iter().find(|s| s.name == name)
+}
+
+fn diffuse(r: f64, g: f64, b: f64) -> Material {
+    Material::Diffuse { albedo: Vec3::new(r, g, b) }
+}
+
+fn mirror(r: f64, g: f64, b: f64) -> Material {
+    Material::Mirror { albedo: Vec3::new(r, g, b) }
+}
+
+/// Glass; `tint` is applied at every reflection/refraction, so colour
+/// deepens with each pass through the surface.
+fn glass(tint: Vec3) -> Material {
+    Material::Dielectric { ref_idx: 1.5, albedo: tint }
+}
+
+fn cornell() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(278.0, 278.0, -800.0),
+            lookat: Vec3::new(278.0, 278.0, 278.0),
+            vfov: 40.0,
+        },
+        Background::Sky,
+    );
+
+    let white = diffuse(0.73, 0.73, 0.73);
+    b.plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), diffuse(0.5, 0.5, 0.5));
+    b.plane(Vec3::new(0.0, 500.0, 0.0), Vec3::new(0.0, -1.0, 0.0), white);
+    b.plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0), white);
+    b.plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0), diffuse(0.65, 0.05, 0.05));
+    b.plane(Vec3::new(500.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), diffuse(0.12, 0.45, 0.15));
+
+    b.box_light(Vec3::new(0.0, 180.0, 180.0), Vec3::new(25.0, 220.0, 320.0), Vec3::splat(50.0));
+
+    b.sphere(Vec3::new(250.0, 180.0, 250.0), 80.0, glass(Vec3::ONE));
+    b.sphere(Vec3::new(150.0, 180.0, 380.0), 80.0, mirror(0.95, 0.95, 0.95));
+    b.sphere(Vec3::new(350.0, 160.0, 300.0), 60.0, diffuse(0.65, 0.05, 0.05));
+    b.build()
+}
+
+fn caustics() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.0, 160.0, -280.0),
+            lookat: Vec3::new(20.0, 30.0, 120.0),
+            vfov: 38.0,
+        },
+        Background::Color(Vec3::ZERO),
+    );
+
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.8, 0.8, 0.8));
+    b.plane(Vec3::new(0.0, 0.0, 520.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.6, 0.6, 0.6));
+
+    // One downward-facing panel above and well behind the spheres. A glass
+    // sphere focuses light about 1.5 radii from its centre; with the light
+    // ~48 degrees from vertical that focus lands on the floor just in front
+    // of each sphere, facing the camera. A smaller panel gives sharper
+    // caustics but makes them slower to converge.
+    b.rect_light(
+        Vec3::new(-30.0, 250.0, 330.0),
+        Vec3::new(30.0, 250.0, 390.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(90.0),
+    );
+
+    b.sphere(Vec3::new(-75.0, 45.0, 110.0), 45.0, glass(Vec3::ONE));
+    b.sphere(Vec3::new(40.0, 38.0, 80.0), 38.0, glass(Vec3::new(1.0, 0.7, 0.3)));
+    b.sphere(Vec3::new(125.0, 30.0, 150.0), 30.0, glass(Vec3::new(0.4, 0.65, 1.0)));
+    b.sphere(Vec3::new(-10.0, 35.0, 210.0), 35.0, mirror(0.9, 0.9, 0.9));
+    b.build()
+}
+
+/// A grid of small random spheres around three large feature spheres (after
+/// "Ray Tracing in One Weekend"). A `light_fraction` of the small spheres glow.
+fn sphere_field(b: &mut SceneBuilder, light_fraction: f64) {
+    let mut rng = SmallRng::seed_from_u64(7);
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.5, 0.5, 0.5));
+
+    let big = [
+        (Vec3::new(0.0, 1.0, 0.0), glass(Vec3::ONE)),
+        (Vec3::new(-4.0, 1.0, 0.0), diffuse(0.4, 0.2, 0.1)),
+        (Vec3::new(4.0, 1.0, 0.0), mirror(0.7, 0.6, 0.5)),
+    ];
+    for (center, material) in big {
+        b.sphere(center, 1.0, material);
+    }
+
+    for a in -11..11 {
+        for c in -11..11 {
+            let center = Vec3::new(
+                a as f64 + 0.9 * rng.gen::<f64>(),
+                0.2,
+                c as f64 + 0.9 * rng.gen::<f64>(),
+            );
+            if big.iter().any(|(bc, _)| (center - Vec3::new(bc.x, 0.2, bc.z)).length() < 1.2) {
+                continue;
+            }
+            let rnd_color = |rng: &mut SmallRng| Vec3::new(rng.gen(), rng.gen(), rng.gen());
+            if rng.gen::<f64>() < light_fraction {
+                // Warm-to-white glow of varying strength.
+                let warmth: f64 = rng.gen();
+                let color = Vec3::new(1.0, 0.55 + 0.4 * warmth, 0.2 + 0.6 * warmth);
+                b.sphere_light(center, 0.2, color * (6.0 + 10.0 * rng.gen::<f64>()));
+                continue;
+            }
+            let choice: f64 = rng.gen();
+            let material = if choice < 0.75 {
+                let albedo = rnd_color(&mut rng).mul_elem(rnd_color(&mut rng));
+                Material::Diffuse { albedo }
+            } else if choice < 0.9 {
+                Material::Mirror { albedo: Vec3::splat(0.5) + 0.5 * rnd_color(&mut rng) }
+            } else {
+                glass(Vec3::ONE)
+            };
+            b.sphere(center, 0.2, material);
+        }
+    }
+}
+
+fn spheres() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(13.0, 2.0, 3.0),
+            lookat: Vec3::ZERO,
+            vfov: 20.0,
+        },
+        Background::Sky,
+    );
+    sphere_field(&mut b, 0.0);
+    b.build()
+}
+
+fn night() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(13.0, 2.0, 3.0),
+            lookat: Vec3::ZERO,
+            vfov: 20.0,
+        },
+        Background::Color(Vec3::new(0.004, 0.006, 0.015)),
+    );
+    sphere_field(&mut b, 0.12);
+    b.build()
+}
+
+fn classic() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(278.0, 278.0, -800.0),
+            lookat: Vec3::new(278.0, 278.0, 0.0),
+            vfov: 40.0,
+        },
+        Background::Color(Vec3::ZERO),
+    );
+
+    let white = diffuse(0.73, 0.73, 0.73);
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), white);
+    b.plane(Vec3::new(0.0, 555.0, 0.0), Vec3::new(0.0, -1.0, 0.0), white);
+    b.plane(Vec3::new(0.0, 0.0, 555.0), Vec3::new(0.0, 0.0, -1.0), white);
+    // The camera looks down +z with +x to its left, so red ends up on the left as usual.
+    b.plane(Vec3::new(555.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), diffuse(0.65, 0.05, 0.05));
+    b.plane(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), diffuse(0.12, 0.45, 0.15));
+
+    // A flat panel just under the ceiling, shining down only.
+    b.rect_light(
+        Vec3::new(213.0, 554.0, 227.0),
+        Vec3::new(343.0, 554.0, 332.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(15.0),
+    );
+
+    b.cuboid(Vec3::new(130.0, 0.0, 65.0), Vec3::new(295.0, 165.0, 230.0), white);
+    b.cuboid(Vec3::new(265.0, 0.0, 295.0), Vec3::new(430.0, 330.0, 460.0), white);
+    b.build()
+}

rust/src/vec3.rs

@@ -56,6 +56,11 @@ impl Vec3 {
     pub fn mean(self) -> f64 {
         (self.x + self.y + self.z) / 3.0
     }
+
+    /// Rec. 709 luminance of a linear RGB colour.
+    pub fn luminance(self) -> f64 {
+        0.2126 * self.x + 0.7152 * self.y + 0.0722 * self.z
+    }
 }
 
 impl Index<usize> for Vec3 {