fog rendering

README.md

@@ -68,7 +68,73 @@ python raytracer.py 800 600 100 15 cornell_hq.png
 
 ### Workers
 
-Multiprocessing uses 4 workers by default. Change `workers=4` in the `render()` call.
+Multiprocessing uses one worker per CPU core by default. Pass `workers=N` to the `render()` call to change it.
+
+## Rust version
+
+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 --scene beam 800 600 1024 15 beam.png
+./target/release/pathtracer --list-scenes
+```
+
+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`.
+
+### 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 |
+| `fog` | The `caustics` scene in light fog: lit haze, a glow around the light and caustics under the spheres, from caustic photon maps |
+| `beam` | After the classic photon-mapping image of a volume caustic: a spotlight shines down through a glass sphere in dense fog, and the sphere focuses its beam into a cone that narrows to a focus and widens again onto the floor. Photon mapped; ~6 min at 800x600, 1024 spp |
+| `media` | Media inside objects: jade (glass with a green scattering interior), smoke, amber liquid (pure absorption) and a cloud |
+| `haze` | The `city` under a layer of atmospheric haze that fades the horizon and glows around the sun |
+| `city` | Massive scale: a procedural 400x400-lot city of ~124,000 objects (buildings with stepped tiers, mirror-glass towers, parks with trees) under a low sun. Renders 800x600 at 64 spp in ~1 s; without the BVH it is ~1000x slower |
+
+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`, `spot_light`, `fog`, `caustic_photons`, `media_photons`). Every emitter must be added with a `*_light` method so light sampling knows about it.
+
+### Volumetrics
+
+Homogeneous participating media ([`rust/src/medium.rs`](rust/src/medium.rs)), each with a coloured absorption coefficient `sigma_a`, a grey scattering coefficient `sigma_s` and a Henyey-Greenstein anisotropy `g`, can be placed:
+
+- as **fog** filling an axis-aligned region outside all objects (`SceneBuilder::fog`, e.g. a haze layer below some height, so the sky stays visible above it);
+- **inside glass**: `Material::Dielectric { medium: Some(..), .. }` (jade, milk, coloured liquids);
+- inside an **invisible boundary**: `Material::Volume { medium }` (smoke, clouds).
+
+Paths sample a scattering distance from `sigma_s` and apply absorption as an exact `exp(-sigma_a * d)` weight, so coloured absorption needs no spectral tricks. At a scattering point they sample a light like a diffuse surface does (with the phase function instead of the BRDF) and continue in a Henyey-Greenstein direction. Shadow rays are attenuated by the media they cross and pass through volume boundaries. `cargo test --release` runs white-furnace tests (a non-absorbing medium in front of a uniform background must render exactly 1), a Beer-Lambert test and phase function tests. Scenes without media render exactly as before.
+
+### Photon mapping for caustics
+
+Path tracing finds a caustic by starting at the camera: a point in fog or on the floor has to pick a direction that happens to pass through the glass and hit the light. With a small light, and especially a spotlight, that almost never happens, so volumetric caustics come out as fireflies or not at all. Scenes can opt in to photon maps ([`rust/src/photons.rs`](rust/src/photons.rs)): before rendering, photons are shot from the lights, and those that pass through glass or off mirrors are stored where they first land in fog or on a diffuse surface. Camera paths gather them along their fog segments (the beam radiance estimate of Jarosz et al.) and at diffuse surfaces, and path tracing skips exactly those caustic paths so nothing is counted twice.
+
+- `caustic_photons(emitted, volume_radius, surface_radius)`: all lights emit photons; the maps replace caustics at a path's first fog or diffuse interaction.
+- `media_photons(...)`: only spotlights emit photons, and their photons also carry their direct light into the fog (much smoother beams than sampling the spotlight from each fog point). Other lights are path traced as usual.
+
+`--photons N` overrides the photon count and `--no-photons` renders with path tracing alone. Gathering blurs by the radius, so photon-mapped images converge to a slightly blurred result; checked against pure path tracing at 4096-8192 spp, mean brightness agrees to within 0.3%. Spotlights (`spot_light`) emit within a cone with a soft edge, and light selection gives a spotlight zero probability for points outside its cone.
+
+### Beyond the Python version
+
+- **8-wide BVH with AVX-512** - bounded objects live in a bounding volume hierarchy (median split, collapsed from binary to 8 children per node); planes are tested separately. Each node stores its children's boxes as structure-of-arrays, so one AVX-512 sequence tests a ray against all 8. Leaves hold up to 8 spheres, also intersected 8 at a time, but at most one other shape (boxes and rects are still tested one by one, so tight leaves are faster for them). The CPU is checked at runtime; without AVX-512, or with `PATHTRACER_NO_AVX512=1` set, the same traversal runs with scalar kernels and renders bit-identical images. Compared with the earlier binary BVH on a Ryzen 7 9800X3D: `city` 2.6x, `spheres` 2.1x, `night` 1.8x faster; the small scenes are unchanged within ±10%.
+- **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.
+- **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
 
@@ -84,7 +150,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

@@ -1,6 +1,7 @@
 import numpy as np
 from PIL import Image
 import math
+import os
 import time
 
 EPSILON = 1e-6
@@ -25,6 +26,14 @@ def refract(incident, normal, eta):
     return eta * incident + (eta * cos_i - cos_t) * normal
 
 
+def random_unit_vector(rng):
+    while True:
+        v = rng.standard_normal(3)
+        n = np.linalg.norm(v)
+        if n > 1e-12:
+            return v / n
+
+
 def schlick(cos, ref_idx):
     r0 = (1.0 - ref_idx) / (1.0 + ref_idx)
     r0 *= r0
@@ -148,8 +157,9 @@ class Box(Hittable):
         if t_lo >= t_max or t_lo <= t_min:
             return None
         rec = HitRecord()
+        # A ray travelling +axis enters through the min face, whose outward normal is -axis.
         n = np.zeros(3)
-        n[hit_axis] = 1.0 if inv_ds[hit_axis] > 0 else -1.0
+        n[hit_axis] = -1.0 if inv_ds[hit_axis] > 0 else 1.0
         if np.dot(r.direction, n) < 0:
             rec.normal = n
             rec.front_face = True
@@ -209,7 +219,9 @@ class Diffuse(Material):
         self.albedo = np.asarray(albedo, dtype=np.float64)
 
     def scatter(self, ray_in, rec, rng):
-        target = rec.normal + np.random.randn(3) * 0.5
+        # normal + uniform unit vector is cosine-distributed, so the Lambertian
+        # BRDF * cos / pdf reduces to just the albedo.
+        target = rec.normal + random_unit_vector(rng)
         if np.linalg.norm(target) < 1e-8:
             target = rec.normal.copy()
         return True, self.albedo, Ray(rec.point, normalize(target))
@@ -241,7 +253,7 @@ class Dielectric(Material):
         cos_theta = min(-np.dot(unit_direction, rec.normal), 1.0)
         sin_theta = math.sqrt(max(0.0, 1.0 - cos_theta * cos_theta))
         cannot_refract = refraction_ratio * sin_theta > 1.0
-        if cannot_refract or schlick(cos_theta, refraction_ratio) > np.random.random():
+        if cannot_refract or schlick(cos_theta, refraction_ratio) > rng.random():
             direction = reflect(unit_direction, rec.normal)
         else:
             direction = refract(unit_direction, rec.normal, refraction_ratio)
@@ -271,15 +283,15 @@ 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
 
     def get_ray(self, s, t, rng):
         rd = np.zeros(3)
         if self.aperture > 0:
-            rd = (np.random.rand(3) - 0.5) * self.aperture
+            rd = (rng.random(3) - 0.5) * self.aperture
             rd[2] = 0.0
         direction = self.lower_left_corner + s * self.horizontal + t * self.vertical - self.origin + rd
         return Ray(self.origin, direction)
@@ -302,7 +314,7 @@ def cornell_box():
     scene.add(Plane(np.array([0, 0, 0]), np.array([1, 0, 0]), red))
     scene.add(Plane(np.array([500, 0, 0]), np.array([-1, 0, 0]), green))
 
-    scene.add(Box(np.array([0, 180, 180]), np.array([25, 220, 320]), light))
+    scene.add_light(Box(np.array([0, 180, 180]), np.array([25, 220, 320]), light))
 
     scene.add(Sphere(np.array([250, 180, 250]), 80, glass))
     scene.add(Sphere(np.array([150, 180, 380]), 80, mirror))
@@ -311,7 +323,48 @@ def cornell_box():
     return scene
 
 
-def trace(scene, ray, depth, max_depth, rng):
+SHADOW_EPSILON = 1e-3
+
+
+def sample_box_light(scene, light_obj, rec, albedo, rng):
+    """Direct light from one uniformly sampled point on the box light's surface."""
+    lo, hi = light_obj.min_pt, light_obj.max_pt
+    dx, dy, dz = hi - lo
+    face_areas = np.array([dy * dz, dx * dz, dx * dy])
+    total_area = 2.0 * face_areas.sum()
+
+    # Pick a face with probability proportional to its area, then a uniform
+    # point on it: overall pdf is 1 / total_area.
+    u = rng.random() * face_areas.sum()
+    axis = 0 if u < face_areas[0] else (1 if u < face_areas[0] + face_areas[1] else 2)
+    positive_side = rng.random() < 0.5
+    lp = lo + rng.random(3) * (hi - lo)
+    lp[axis] = hi[axis] if positive_side else lo[axis]
+    light_normal = np.zeros(3)
+    light_normal[axis] = 1.0 if positive_side else -1.0
+
+    to_light = lp - rec.point
+    dist_sq = np.dot(to_light, to_light)
+    if dist_sq < 1e-8:
+        return np.zeros(3)
+    dist = math.sqrt(dist_sq)
+    d_to_light = to_light / dist
+    cos_at_rec = np.dot(rec.normal, d_to_light)
+    cos_at_light = -np.dot(light_normal, d_to_light)
+    if cos_at_rec <= 0 or cos_at_light <= 0:
+        return np.zeros(3)
+
+    if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist - SHADOW_EPSILON):
+        return np.zeros(3)
+
+    brdf = albedo / PI
+    geometry = cos_at_rec * cos_at_light / dist_sq
+    return light_obj.material.get_emission() * brdf * geometry * total_area
+
+
+def trace(scene, ray, depth, max_depth, rng, count_emission=True):
+    """count_emission is False after a diffuse bounce: that light was already
+    gathered by next event estimation, so hitting it again would double count."""
     if depth >= max_depth:
         return np.zeros(3)
 
@@ -320,67 +373,22 @@ def trace(scene, ray, depth, max_depth, rng):
         t_val = 0.5 * (ray.direction[1] + 1.0)
         return np.array([1.0, 1.0, 1.0]) * (1.0 - t_val) + np.array([0.5, 0.7, 1.0]) * t_val
 
-    emitted = rec.material.get_emission()
-    if isinstance(rec.material, LightSource):
-        return emitted
-
     mat = rec.material
+    if isinstance(mat, LightSource):
+        return mat.get_emission() if count_emission else np.zeros(3)
+
     scattered, attenuation, scattered_ray = mat.scatter(ray, rec, rng)
     if not scattered or scattered_ray is None:
-        return emitted
-
-    indirect = trace(scene, scattered_ray, depth + 1, max_depth, rng)
-    color = emitted + attenuation * indirect
+        return np.zeros(3)
 
     if isinstance(mat, Diffuse):
         direct = np.zeros(3)
         for light_obj in scene.lights:
-            emission = light_obj.material.get_emission()
-            for _ in range(2):
-                face_idx = int(rng.random() * 6)
-                x = rng.uniform(light_obj.min_pt[0], light_obj.max_pt[0])
-                y = rng.uniform(light_obj.min_pt[1], light_obj.max_pt[1])
-                z = rng.uniform(light_obj.min_pt[2], light_obj.max_pt[2])
-                if face_idx == 0:
-                    lp = np.array([light_obj.min_pt[0], y, z])
-                elif face_idx == 1:
-                    lp = np.array([light_obj.max_pt[0], y, z])
-                elif face_idx == 2:
-                    lp = np.array([x, light_obj.min_pt[1], z])
-                elif face_idx == 3:
-                    lp = np.array([x, light_obj.max_pt[1], z])
-                elif face_idx == 4:
-                    lp = np.array([x, y, light_obj.min_pt[2]])
-                else:
-                    lp = np.array([x, y, light_obj.max_pt[2]])
+            direct += sample_box_light(scene, light_obj, rec, mat.albedo, rng)
+        indirect = trace(scene, scattered_ray, depth + 1, max_depth, rng, count_emission=False)
+        return direct + attenuation * indirect
 
-                to_light = lp - rec.point
-                dist_sq = np.dot(to_light, to_light)
-                if dist_sq < 1e-8:
-                    continue
-                dist = math.sqrt(dist_sq)
-                d_to_light = to_light / dist
-                cos_at_rec = np.dot(rec.normal, d_to_light)
-                if cos_at_rec <= 0:
-                    continue
-
-                if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist - EPSILON):
-                    continue
-
-                light_area = 6 * (light_obj.max_pt[0] - light_obj.min_pt[0]) * \
-                             (light_obj.max_pt[1] - light_obj.min_pt[1]) * \
-                             (light_obj.max_pt[2] - light_obj.min_pt[2])
-                pdf = light_area / (dist_sq * cos_at_rec)
-                if pdf < 1e-12:
-                    continue
-                brdf = mat.albedo / PI
-                direct += emission * cos_at_rec * brdf / pdf
-                break
-
-        if np.dot(direct, direct) > 0:
-            color = emitted + attenuation * (indirect + direct)
-
-    return color
+    return attenuation * trace(scene, scattered_ray, depth + 1, max_depth, rng, count_emission=True)
 
 
 def render_chunk(args):
@@ -394,10 +402,10 @@ def render_chunk(args):
     for j_local, j in enumerate(range(chunk_start, chunk_end)):
         for i in range(width):
             color = np.zeros(3)
-            rng = np.random.default_rng(seed_off + i * 1000 + j * 7)
+            rng = np.random.default_rng([seed_off, j, i])
             for _ in range(spp):
-                u = (i + rng.random()) / (width - 1)
-                v = (height - 1 - j + rng.random()) / (height - 1)
+                u = (i + rng.random()) / width
+                v = (height - 1 - j + rng.random()) / height
                 ray = cam.get_ray(u, v, rng)
                 color += trace(scene, ray, 0, max_d, rng)
             color /= spp
@@ -409,22 +417,26 @@ def render_chunk(args):
 
 def render(scene, width=800, height=600, samples_per_pixel=50, max_depth=15,
            lookfrom=None, lookat=None, vfov=40, aperture=0.0, focus_dist=1078.0,
-           output="cornell_box.png", workers=16):
+           output="cornell_box.png", workers=None):
     if lookfrom is None:
         lookfrom = np.array([278.0, 278.0, -800.0])
     if lookat is None:
         lookat = np.array([278.0, 278.0, 278.0])
+    if workers is None:
+        workers = os.cpu_count() or 4
+    workers = max(1, min(workers, height))
 
     image = np.zeros((height, width, 3), dtype=np.float64)
     start = time.time()
 
-    chunk = max(1, height // workers)
+    # Many small row chunks (rather than one per worker) keep the progress bar
+    # moving and balance load between the cheap and expensive parts of the image.
+    rows_per_chunk = max(1, height // (workers * 8))
     ranges = []
-    for w in range(workers):
-        s = w * chunk
-        e = height if w == workers - 1 else (w + 1) * chunk
+    for s in range(0, height, rows_per_chunk):
+        e = min(s + rows_per_chunk, height)
         ranges.append((scene, width, height, samples_per_pixel,
-                       max_depth, lookfrom, lookat, vfov, s, e, w * 100000 + 42))
+                       max_depth, lookfrom, lookat, vfov, s, e, 42))
 
     from concurrent.futures import ProcessPoolExecutor, as_completed
 
@@ -468,6 +480,7 @@ def render(scene, width=800, height=600, samples_per_pixel=50, max_depth=15,
     timestamps = []
     completed_rows = 0
     total_rows = height
+    last_preview = 0.0
 
     with progress:
         render_task = progress.add_task("Rendering", total=total_rows)
@@ -484,10 +497,12 @@ def render(scene, width=800, height=600, samples_per_pixel=50, max_depth=15,
                 progress.update(render_task, completed=completed_rows)
                 progress.refresh()
 
-                img = (np.clip(image, 0.0, 1.0) * 255).astype(np.uint8)
-                Image.fromarray(img, 'RGB').save(preview_path)
-
                 elapsed = time.time() - start
+                if elapsed - last_preview >= 1.0 or completed_rows == total_rows:
+                    last_preview = elapsed
+                    img = (np.clip(image, 0.0, 1.0) * 255).astype(np.uint8)
+                    Image.fromarray(img, 'RGB').save(preview_path)
+
                 brightness = float(image.mean())
                 brightness_history.append(brightness)
                 timestamps.append(elapsed)
@@ -509,12 +524,13 @@ def render(scene, width=800, height=600, samples_per_pixel=50, max_depth=15,
     region_table.add_column("Brightness", justify="right", style="white")
     region_table.add_column("Visual", justify="center")
     h_start, h_end = height // 3, 2 * height // 3
+    w_start, w_end = width // 3, 2 * width // 3
     regions = [
         ("Left Wall", image[:, :width // 4]),
         ("Right Wall", image[:, 3 * width // 4:]),
-        ("Floor", image[2 * height // 3:, h_start:h_end]),
-        ("Ceiling", image[:height // 3, h_start:h_end]),
-        ("Center", image[h_start:h_end, h_start:h_end]),
+        ("Floor", image[h_end:, w_start:w_end]),
+        ("Ceiling", image[:h_start, w_start:w_end]),
+        ("Center", image[h_start:h_end, w_start:w_end]),
     ]
     for name, region in regions:
         brightness = float(region.mean()) * 255

rust/.gitignore

@@ -0,0 +1 @@
+/target

rust/Cargo.lock

@@ -0,0 +1,583 @@
+# This file is automatically @generated by Cargo.
+# It is not intended for manual editing.
+version = 4
+
+[[package]]
+name = "adler2"
+version = "2.0.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa"
+
+[[package]]
+name = "autocfg"
+version = "1.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53"
+
+[[package]]
+name = "bitflags"
+version = "2.13.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3ded4057c258ba199e2d26386d3af3780957ecaee6c4ef4041c6b4b8b97c0b06"
+
+[[package]]
+name = "bumpalo"
+version = "3.20.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "72f5acc6cb2ba439de613abc23857ec3d78374d8ed5ac84e9d11336e87da8649"
+
+[[package]]
+name = "bytemuck"
+version = "1.25.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "95832e849adfb21180ccb6826a99da14e5d266ae5c2e668e1602cf234f153797"
+
+[[package]]
+name = "byteorder-lite"
+version = "0.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8f1fe948ff07f4bd06c30984e69f5b4899c516a3ef74f34df92a2df2ab535495"
+
+[[package]]
+name = "cfg-if"
+version = "1.0.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4e7648175b45a9a48536d676f68d918270699102aa8dab5496df06904c914600"
+
+[[package]]
+name = "console"
+version = "0.15.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "054ccb5b10f9f2cbf51eb355ca1d05c2d279ce1804688d0db74b4733a5aeafd8"
+dependencies = [
+ "encode_unicode",
+ "libc",
+ "once_cell",
+ "unicode-width",
+ "windows-sys",
+]
+
+[[package]]
+name = "crc32fast"
+version = "1.5.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "01a7799fd6b852db0e61728dde9a204c423b44d689dbd432522543614b490e78"
+dependencies = [
+ "cfg-if",
+]
+
+[[package]]
+name = "crossbeam-deque"
+version = "0.8.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "622f3fc73690be383c7214310406f28a90e6edeadc3cea882f9d71e495b9711a"
+dependencies = [
+ "crossbeam-epoch",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.9.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dc74980687109a3b14c72fd458107bf0baa1da1a1a805e178d15501ba9b86d9d"
+dependencies = [
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-utils"
+version = "0.8.23"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a31eee39dddec8330830986fcd7625edb5a24ec90ea038215273bbc3adb08ac6"
+
+[[package]]
+name = "either"
+version = "1.18.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "252afb9ae5eaa683babdc6a068b3f5726eb19e05070c731f9b2a23a7c3e8ed34"
+
+[[package]]
+name = "encode_unicode"
+version = "1.0.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "34aa73646ffb006b8f5147f3dc182bd4bcb190227ce861fc4a4844bf8e3cb2c0"
+
+[[package]]
+name = "fdeflate"
+version = "0.3.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1e6853b52649d4ac5c0bd02320cddc5ba956bdb407c4b75a2c6b75bf51500f8c"
+dependencies = [
+ "simd-adler32",
+]
+
+[[package]]
+name = "flate2"
+version = "1.1.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6e634e2e0ebac1ee034020da1ca582e17ffe4e0f5e985823721e168928136dcb"
+dependencies = [
+ "crc32fast",
+ "miniz_oxide 0.9.1",
+ "zlib-rs",
+]
+
+[[package]]
+name = "futures-core"
+version = "0.3.34"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "92d699e522242e69e3003b94ecc1f960f3a5e015aa7c5d7486e65ad01dd94f5e"
+
+[[package]]
+name = "futures-task"
+version = "0.3.34"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cd417de3d1d015fc3bfd2b1ea46dfc7bab72ef86f1cc7cc9c78e728b34a6d1fd"
+
+[[package]]
+name = "futures-util"
+version = "0.3.34"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0d50a92467f8ba5dd6e3ee5d4bd04d73ab2e4e1c44474a0674821dfce14b79bc"
+dependencies = [
+ "futures-core",
+ "futures-task",
+ "pin-project-lite",
+ "slab",
+]
+
+[[package]]
+name = "getrandom"
+version = "0.2.17"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ff2abc00be7fca6ebc474524697ae276ad847ad0a6b3faa4bcb027e9a4614ad0"
+dependencies = [
+ "cfg-if",
+ "libc",
+ "wasi",
+]
+
+[[package]]
+name = "image"
+version = "0.25.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "85ab80394333c02fe689eaf900ab500fbd0c2213da414687ebf995a65d5a6104"
+dependencies = [
+ "bytemuck",
+ "byteorder-lite",
+ "moxcms",
+ "num-traits",
+ "png",
+]
+
+[[package]]
+name = "indicatif"
+version = "0.17.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "183b3088984b400f4cfac3620d5e076c84da5364016b4f49473de574b2586235"
+dependencies = [
+ "console",
+ "number_prefix",
+ "portable-atomic",
+ "unicode-width",
+ "web-time",
+]
+
+[[package]]
+name = "js-sys"
+version = "0.3.106"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7883d941dae510fb2d978fc3fe018c71c9e2892fd38854de3e8b92c2e5ad9cc5"
+dependencies = [
+ "cfg-if",
+ "futures-util",
+ "wasm-bindgen",
+]
+
+[[package]]
+name = "libc"
+version = "0.2.189"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3eaf3ede3fee6db1a4c2ee091bf8a8b4dccdc6d17f656fb07896ee72867612f2"
+
+[[package]]
+name = "miniz_oxide"
+version = "0.8.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316"
+dependencies = [
+ "adler2",
+ "simd-adler32",
+]
+
+[[package]]
+name = "miniz_oxide"
+version = "0.9.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b63fbc4a50860e98e7b2aa7804ded1db5cbc3aff9193adaff57a6931bf7c4b4c"
+dependencies = [
+ "adler2",
+ "simd-adler32",
+]
+
+[[package]]
+name = "moxcms"
+version = "0.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "bb85c154ba489f01b25c0d36ae69a87e4a1c73a72631fc6c0eb6dde34a73e44b"
+dependencies = [
+ "num-traits",
+ "pxfm",
+]
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+]
+
+[[package]]
+name = "number_prefix"
+version = "0.4.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "830b246a0e5f20af87141b25c173cd1b609bd7779a4617d6ec582abaf90870f3"
+
+[[package]]
+name = "once_cell"
+version = "1.21.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50"
+
+[[package]]
+name = "pathtracer"
+version = "0.1.0"
+dependencies = [
+ "image",
+ "indicatif",
+ "rand",
+ "rayon",
+]
+
+[[package]]
+name = "pin-project-lite"
+version = "0.2.17"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "a89322df9ebe1c1578d689c92318e070967d1042b512afbe49518723f4e6d5cd"
+
+[[package]]
+name = "png"
+version = "0.18.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60769b8b31b2a9f263dae2776c37b1b28ae246943cf719eb6946a1db05128a61"
+dependencies = [
+ "bitflags",
+ "crc32fast",
+ "fdeflate",
+ "flate2",
+ "miniz_oxide 0.8.9",
+]
+
+[[package]]
+name = "portable-atomic"
+version = "1.15.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "05c8b63e8d9609db387f0324918f81d68fe27748f084ef092fb35954d0539a85"
+
+[[package]]
+name = "ppv-lite86"
+version = "0.2.21"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "85eae3c4ed2f50dcfe72643da4befc30deadb458a9b590d720cde2f2b1e97da9"
+dependencies = [
+ "zerocopy",
+]
+
+[[package]]
+name = "proc-macro2"
+version = "1.0.107"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "985e7ec9bb745e6ce6535b544d84d6cd6f7ad8bd711c398938ae983b91a766d9"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "pxfm"
+version = "0.1.30"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d55d956fa96f5ec02be2e13af0e20391a5aa83d6a074e3ad368959d0fab299ea"
+
+[[package]]
+name = "quote"
+version = "1.0.47"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1fbf4db142a473a8d80c26bbf18454ed458bf8d26c8219c331daecfdbd079001"
+dependencies = [
+ "proc-macro2",
+]
+
+[[package]]
+name = "rand"
+version = "0.8.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e058c7de0b26af77780c769414d6257830bb240f3c38477dbc2c16e5f54d6d4c"
+dependencies = [
+ "libc",
+ "rand_chacha",
+ "rand_core",
+]
+
+[[package]]
+name = "rand_chacha"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88"
+dependencies = [
+ "ppv-lite86",
+ "rand_core",
+]
+
+[[package]]
+name = "rand_core"
+version = "0.6.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
+dependencies = [
+ "getrandom",
+]
+
+[[package]]
+name = "rayon"
+version = "1.12.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fb39b166781f92d482534ef4b4b1b2568f42613b53e5b6c160e24cfbfa30926d"
+dependencies = [
+ "either",
+ "rayon-core",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.13.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "22e18b0f0062d30d4230b2e85ff77fdfe4326feb054b9783a3460d8435c8ab91"
+dependencies = [
+ "crossbeam-deque",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "rustversion"
+version = "1.0.23"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cf54715a573b99ac80df0bc206da022bcd442c974952c7b9720069370852e21f"
+
+[[package]]
+name = "simd-adler32"
+version = "0.3.10"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3a219298ac11a56ea9a6d2120044824d6f01aeb034955e7af7bc16858527deea"
+
+[[package]]
+name = "slab"
+version = "0.4.12"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0c790de23124f9ab44544d7ac05d60440adc586479ce501c1d6d7da3cd8c9cf5"
+
+[[package]]
+name = "syn"
+version = "2.0.119"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "872831b642d1a07999a962a351ed35b955ea2cfc8f3862091e2a240a84f17297"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "syn"
+version = "3.0.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8593e8e72159ed2257d083c7a454a85cbf854f37a0966d8d483aff8c8a3ebcee"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "unicode-ident",
+]
+
+[[package]]
+name = "unicode-ident"
+version = "1.0.26"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d245f478577f809a851594d02313b640fb437e0bb33866753cff937863096954"
+
+[[package]]
+name = "unicode-width"
+version = "0.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b4ac048d71ede7ee76d585517add45da530660ef4390e49b098733c6e897f254"
+
+[[package]]
+name = "wasi"
+version = "0.11.1+wasi-snapshot-preview1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b"
+
+[[package]]
+name = "wasm-bindgen"
+version = "0.2.129"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9bb54f33acc68fd454578d9820b0bde1a1a3d17aa17bb7b6595806d02886d409"
+dependencies = [
+ "cfg-if",
+ "once_cell",
+ "rustversion",
+ "wasm-bindgen-macro",
+ "wasm-bindgen-shared",
+]
+
+[[package]]
+name = "wasm-bindgen-macro"
+version = "0.2.129"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2e29d0c35b16e224a7eeb5cd2d25e3e1968fbd65604117b44d3b789d00ee8535"
+dependencies = [
+ "quote",
+ "wasm-bindgen-macro-support",
+]
+
+[[package]]
+name = "wasm-bindgen-macro-support"
+version = "0.2.129"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6f501a8bc3719dba86ef8ae4728879c08001bea749eb1333ac5b91e040e2a6b7"
+dependencies = [
+ "bumpalo",
+ "proc-macro2",
+ "quote",
+ "syn 3.0.6",
+ "wasm-bindgen-shared",
+]
+
+[[package]]
+name = "wasm-bindgen-shared"
+version = "0.2.129"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "23f0c9c52aa7cd7d77769a4cfe2a9adb1b331f489a41d912ce14513d5ab995c6"
+dependencies = [
+ "unicode-ident",
+]
+
+[[package]]
+name = "web-time"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5a6580f308b1fad9207618087a65c04e7a10bc77e02c8e84e9b00dd4b12fa0bb"
+dependencies = [
+ "js-sys",
+ "wasm-bindgen",
+]
+
+[[package]]
+name = "windows-sys"
+version = "0.59.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b"
+dependencies = [
+ "windows-targets",
+]
+
+[[package]]
+name = "windows-targets"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973"
+dependencies = [
+ "windows_aarch64_gnullvm",
+ "windows_aarch64_msvc",
+ "windows_i686_gnu",
+ "windows_i686_gnullvm",
+ "windows_i686_msvc",
+ "windows_x86_64_gnu",
+ "windows_x86_64_gnullvm",
+ "windows_x86_64_msvc",
+]
+
+[[package]]
+name = "windows_aarch64_gnullvm"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3"
+
+[[package]]
+name = "windows_aarch64_msvc"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469"
+
+[[package]]
+name = "windows_i686_gnu"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b"
+
+[[package]]
+name = "windows_i686_gnullvm"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66"
+
+[[package]]
+name = "windows_i686_msvc"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66"
+
+[[package]]
+name = "windows_x86_64_gnu"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78"
+
+[[package]]
+name = "windows_x86_64_gnullvm"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d"
+
+[[package]]
+name = "windows_x86_64_msvc"
+version = "0.52.6"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec"
+
+[[package]]
+name = "zerocopy"
+version = "0.8.59"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6df92bf3d9227be3d53173901ddbffac2babc27ae50f397776ffd6dc33f800cb"
+dependencies = [
+ "zerocopy-derive",
+]
+
+[[package]]
+name = "zerocopy-derive"
+version = "0.8.59"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ac4f328cf2f05d084e496c3e9c3f33ed0a183656a16e1fcec4d464d8373aec82"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 2.0.119",
+]
+
+[[package]]
+name = "zlib-rs"
+version = "0.6.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b268e58e7c693d7c271f93ffc4ba3b380412554231c85bf61ca7af91042a4112"

rust/Cargo.toml

@@ -0,0 +1,16 @@
+[package]
+name = "pathtracer"
+version = "0.1.0"
+edition = "2021"
+description = "Cornell Box path tracer (Rust port of raytracer.py)"
+
+[dependencies]
+rand = { version = "0.8", features = ["small_rng"] }
+rayon = "1"
+indicatif = "0.17"
+image = { version = "0.25", default-features = false, features = ["png"] }
+
+[profile.release]
+opt-level = 3
+lto = true
+codegen-units = 1

rust/src/bvh.rs

@@ -0,0 +1,687 @@
+//! 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.
+//!
+//! The tree is built as a binary BVH (median split) and then collapsed into
+//! an 8-wide BVH: each node holds the bounds of up to 8 children in
+//! structure-of-arrays form, so one AVX-512 instruction sequence tests a ray
+//! against all 8 boxes. Leaves hold up to 8 objects; spheres in a leaf are
+//! intersected 8 at a time the same way, other shapes one by one through a
+//! callback. Without AVX-512 (or with PATHTRACER_NO_AVX512 set) the same
+//! traversal runs with scalar kernels.
+
+use std::mem::MaybeUninit;
+
+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
+    }
+
+    fn surface_area(&self) -> f64 {
+        let d = self.max - self.min;
+        2.0 * (d.x * d.y + d.y * d.z + d.x * d.z)
+    }
+}
+
+const WIDTH: usize = 8;
+/// Leaves hold up to WIDTH objects, since spheres are tested WIDTH at a time,
+/// but at most this many other shapes, which are tested one by one: for
+/// box-heavy scenes, small leaves (more SIMD node tests, fewer scalar box
+/// tests) measured fastest.
+const MAX_OTHERS_PER_LEAF: usize = 1;
+
+// ---------------------------------------------------------------------------
+// Build: binary BVH first.
+
+struct BinNode {
+    bounds: Aabb,
+    start: usize,
+    count: usize, // > 0 for leaves
+    left: usize,
+    right: usize,
+}
+
+/// Per-object inputs to the build.
+struct BuildInput<'a> {
+    bounds: &'a [Aabb],
+    centroids: Vec<Vec3>,
+    is_sphere: Vec<bool>,
+}
+
+fn build_binary(nodes: &mut Vec<BinNode>, prims: &mut [u32], offset: usize, input: &BuildInput) -> usize {
+    let (bounds, centroids) = (input.bounds, &input.centroids);
+    let node_bounds = prims.iter().fold(Aabb::EMPTY, |b, &p| b.union(bounds[p as usize]));
+    let idx = nodes.len();
+    nodes.push(BinNode { bounds: node_bounds, start: offset, count: prims.len(), left: 0, right: 0 });
+    let others = prims.iter().filter(|&&p| !input.is_sphere[p as usize]).count();
+    if prims.len() <= WIDTH && others <= MAX_OTHERS_PER_LEAF {
+        return idx;
+    }
+
+    // Median split along the axis where the object centroids spread widest.
+    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 half = prims.len() / 2;
+    // If all centroids coincide the order doesn't matter; just halve.
+    if extent[axis] > 0.0 {
+        prims.select_nth_unstable_by(half, |&a, &b| centroids[a as usize][axis].total_cmp(&centroids[b as usize][axis]));
+    }
+    let (lo, hi) = prims.split_at_mut(half);
+    let left = build_binary(nodes, lo, offset, input);
+    let right = build_binary(nodes, hi, offset + half, input);
+    nodes[idx].count = 0;
+    nodes[idx].left = left;
+    nodes[idx].right = right;
+    idx
+}
+
+// ---------------------------------------------------------------------------
+// 8-wide nodes.
+
+/// Child references: an interior child is its node index; a leaf child has
+/// LEAF_BIT set, flags saying whether it holds spheres and/or other shapes
+/// (so traversal can skip the pass it doesn't need), (count - 1) in bits
+/// 26..29 and its first slot in bits 0..26.
+const LEAF_BIT: u32 = 1 << 31;
+const HAS_OTHERS: u32 = 1 << 30;
+const HAS_SPHERES: u32 = 1 << 29;
+const START_BITS: u32 = 26;
+const START_MASK: u32 = (1 << START_BITS) - 1;
+
+fn leaf_ref(start: usize, count: usize, has_spheres: bool, has_others: bool) -> u32 {
+    assert!(start < (1 << START_BITS), "too many objects for the BVH leaf encoding");
+    debug_assert!((1..=WIDTH).contains(&count));
+    let mut r = LEAF_BIT | (((count - 1) as u32) << START_BITS) | start as u32;
+    if has_spheres {
+        r |= HAS_SPHERES;
+    }
+    if has_others {
+        r |= HAS_OTHERS;
+    }
+    r
+}
+
+fn decode_leaf(r: u32) -> (usize, usize) {
+    ((r & START_MASK) as usize, (((r >> START_BITS) & 7) + 1) as usize)
+}
+
+/// Child bounds in structure-of-arrays layout. Unused lanes have all bounds
+/// at +inf: every slab distance is then ±inf, so the lane can never be hit
+/// and needs no separate mask, provided t_max is finite (with t_max = inf a
+/// ray with all-positive direction would see inf <= inf and "hit" it, which
+/// is why the traversal clamps t_max to f64::MAX).
+#[repr(C, align(64))]
+struct Node8 {
+    min_x: [f64; WIDTH],
+    min_y: [f64; WIDTH],
+    min_z: [f64; WIDTH],
+    max_x: [f64; WIDTH],
+    max_y: [f64; WIDTH],
+    max_z: [f64; WIDTH],
+    child: [u32; WIDTH],
+}
+
+impl Node8 {
+    fn empty() -> Self {
+        let inf = [f64::INFINITY; WIDTH];
+        Node8 { min_x: inf, min_y: inf, min_z: inf, max_x: inf, max_y: inf, max_z: inf, child: [0; WIDTH] }
+    }
+}
+
+/// Sphere data per slot, in BVH order, padded by WIDTH entries so a leaf's
+/// 8 lanes can always be loaded. Non-sphere slots have r2 = -inf, which
+/// makes the discriminant -inf so they never register a hit.
+struct Spheres {
+    cx: Vec<f64>,
+    cy: Vec<f64>,
+    cz: Vec<f64>,
+    r2: Vec<f64>,
+}
+
+/// Per-ray values shared by all box and sphere tests.
+#[derive(Clone, Copy)]
+struct RayData {
+    o: Vec3,
+    d: Vec3,
+    inv: Vec3,
+    /// d·d, as the scalar sphere test computes it.
+    a: f64,
+}
+
+impl RayData {
+    fn new(r: &Ray) -> Self {
+        // A finite stand-in for 1/0 keeps (bound - origin) * inv free of the
+        // NaN that 0 * inf would produce.
+        let inv = |d: f64| {
+            let i = 1.0 / d;
+            if i.is_finite() {
+                i
+            } else {
+                1e300_f64.copysign(i)
+            }
+        };
+        RayData {
+            o: r.origin,
+            d: r.direction,
+            inv: Vec3::new(inv(r.direction.x), inv(r.direction.y), inv(r.direction.z)),
+            a: r.direction.dot(r.direction),
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Scalar kernels.
+
+/// Returns a bitmask of the child boxes the ray hits within [t_min, t_max],
+/// writing each child's entry distance to `tnear`.
+#[inline(always)]
+fn node_scalar(n: &Node8, r: &RayData, t_min: f64, t_max: f64, tnear: &mut [f64; WIDTH]) -> u8 {
+    let mut mask = 0u8;
+    for i in 0..WIDTH {
+        let (x0, x1) = ((n.min_x[i] - r.o.x) * r.inv.x, (n.max_x[i] - r.o.x) * r.inv.x);
+        let (y0, y1) = ((n.min_y[i] - r.o.y) * r.inv.y, (n.max_y[i] - r.o.y) * r.inv.y);
+        let (z0, z1) = ((n.min_z[i] - r.o.z) * r.inv.z, (n.max_z[i] - r.o.z) * r.inv.z);
+        let near = x0.min(x1).max(y0.min(y1)).max(z0.min(z1).max(t_min));
+        let far = x0.max(x1).min(y0.max(y1)).min(z0.max(z1).min(t_max));
+        tnear[i] = near;
+        if near <= far {
+            mask |= 1 << i;
+        }
+    }
+    mask
+}
+
+/// Intersects the ray with the spheres in slots start..start + count.
+/// Returns a bitmask of lanes hit within [t_min, t_max] and each lane's
+/// hit distance. Mirrors Shape::hit for spheres operation by operation.
+#[inline(always)]
+fn spheres_scalar(s: &Spheres, start: usize, count: usize, r: &RayData, t_min: f64, t_max: f64) -> (u8, [f64; WIDTH]) {
+    let mut mask = 0u8;
+    let mut roots = [0.0; WIDTH];
+    for i in 0..count {
+        let k = start + i;
+        let oc = Vec3::new(r.o.x - s.cx[k], r.o.y - s.cy[k], r.o.z - s.cz[k]);
+        let b = oc.dot(r.d);
+        let c = oc.dot(oc) - s.r2[k];
+        let disc = b * b - r.a * c;
+        if disc < 0.0 || disc.is_nan() {
+            continue;
+        }
+        let sq = disc.sqrt();
+        let mut root = (-b - sq) / r.a;
+        if root < t_min || root > t_max {
+            root = (-b + sq) / r.a;
+            if root < t_min || root > t_max {
+                continue;
+            }
+        }
+        roots[i] = root;
+        mask |= 1 << i;
+    }
+    (mask, roots)
+}
+
+// ---------------------------------------------------------------------------
+// AVX-512 kernels: the same arithmetic on 8 lanes of f64 at once.
+
+#[cfg(target_arch = "x86_64")]
+mod avx512 {
+    use super::{Node8, RayData, Spheres, WIDTH};
+    use std::arch::x86_64::*;
+
+    #[target_feature(enable = "avx512f")]
+    #[inline]
+    pub unsafe fn node(n: &Node8, r: &RayData, t_min: f64, t_max: f64, tnear: &mut [f64; WIDTH]) -> u8 {
+        let slab = |lo: &[f64; WIDTH], hi: &[f64; WIDTH], o: f64, inv: f64| {
+            let (o, inv) = (_mm512_set1_pd(o), _mm512_set1_pd(inv));
+            let t0 = _mm512_mul_pd(_mm512_sub_pd(_mm512_loadu_pd(lo.as_ptr()), o), inv);
+            let t1 = _mm512_mul_pd(_mm512_sub_pd(_mm512_loadu_pd(hi.as_ptr()), o), inv);
+            (_mm512_min_pd(t0, t1), _mm512_max_pd(t0, t1))
+        };
+        let (nx, fx) = slab(&n.min_x, &n.max_x, r.o.x, r.inv.x);
+        let (ny, fy) = slab(&n.min_y, &n.max_y, r.o.y, r.inv.y);
+        let (nz, fz) = slab(&n.min_z, &n.max_z, r.o.z, r.inv.z);
+        let near = _mm512_max_pd(_mm512_max_pd(nx, ny), _mm512_max_pd(nz, _mm512_set1_pd(t_min)));
+        let far = _mm512_min_pd(_mm512_min_pd(fx, fy), _mm512_min_pd(fz, _mm512_set1_pd(t_max)));
+        _mm512_storeu_pd(tnear.as_mut_ptr(), near);
+        _mm512_cmp_pd_mask::<_CMP_LE_OQ>(near, far)
+    }
+
+    #[target_feature(enable = "avx512f")]
+    #[inline]
+    pub unsafe fn spheres(
+        s: &Spheres,
+        start: usize,
+        count: usize,
+        r: &RayData,
+        t_min: f64,
+        t_max: f64,
+    ) -> (u8, [f64; WIDTH]) {
+        let mut roots = [0.0; WIDTH];
+        let lanes = ((1u32 << count) - 1) as u8;
+        let load = |v: &Vec<f64>| _mm512_loadu_pd(v.as_ptr().add(start));
+        let (add, sub, mul) = (
+            |a, b| _mm512_add_pd(a, b),
+            |a, b| _mm512_sub_pd(a, b),
+            |a, b| _mm512_mul_pd(a, b),
+        );
+
+        let ocx = sub(_mm512_set1_pd(r.o.x), load(&s.cx));
+        let ocy = sub(_mm512_set1_pd(r.o.y), load(&s.cy));
+        let ocz = sub(_mm512_set1_pd(r.o.z), load(&s.cz));
+        let (dx, dy, dz) = (_mm512_set1_pd(r.d.x), _mm512_set1_pd(r.d.y), _mm512_set1_pd(r.d.z));
+        let b = add(add(mul(ocx, dx), mul(ocy, dy)), mul(ocz, dz));
+        let c = sub(add(add(mul(ocx, ocx), mul(ocy, ocy)), mul(ocz, ocz)), load(&s.r2));
+        let a = _mm512_set1_pd(r.a);
+        let disc = sub(mul(b, b), mul(a, c));
+        let zero = _mm512_setzero_pd();
+        let mut mask = lanes & _mm512_cmp_pd_mask::<_CMP_GE_OQ>(disc, zero);
+        if mask == 0 {
+            return (0, roots);
+        }
+
+        let sq = _mm512_sqrt_pd(disc);
+        let neg_b = sub(zero, b);
+        let r1 = _mm512_div_pd(sub(neg_b, sq), a);
+        let r2 = _mm512_div_pd(add(neg_b, sq), a);
+        let (lo, hi) = (_mm512_set1_pd(t_min), _mm512_set1_pd(t_max));
+        let in_range = |t| _mm512_cmp_pd_mask::<_CMP_GE_OQ>(t, lo) & _mm512_cmp_pd_mask::<_CMP_LE_OQ>(t, hi);
+        let ok1 = in_range(r1);
+        let ok2 = in_range(r2);
+        // Nearer root when it is in range, as in the scalar version.
+        let root = _mm512_mask_blend_pd(ok1, r2, r1);
+        mask &= ok1 | ok2;
+        _mm512_storeu_pd(roots.as_mut_ptr(), root);
+        (mask, roots)
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Traversal, written once and instantiated for each set of kernels (a macro
+// rather than generics so the AVX-512 kernels inline into a function that is
+// itself compiled with AVX-512 enabled).
+
+/// Traversal pushes at most WIDTH - 1 entries per level beyond the one it pops,
+/// so this covers trees up to (STACK_SIZE - 1) / 7 = 36 levels deep (checked in build).
+const STACK_SIZE: usize = 256;
+
+macro_rules! traversal {
+    ($(#[$attr:meta])* $closest:ident, $any:ident, $node:path, $spheres:path) => {
+        /// Closest hit. `other(object, t_max)` intersects non-sphere objects
+        /// and returns the distance if hit closer than t_max. Returns
+        /// (object, t, whether it was a sphere).
+        $(#[$attr])*
+        unsafe fn $closest(
+            &self,
+            root: u32,
+            r: &RayData,
+            t_min: f64,
+            mut t_max: f64,
+            mut other: impl FnMut(usize, f64) -> Option<f64>,
+        ) -> Option<(usize, f64, bool)> {
+            let mut best = None;
+            // Uninitialised: zero-filling these on every ray is measurable in
+            // small scenes. Only entries below sp are ever read.
+            let mut stack_ref = [MaybeUninit::<u32>::uninit(); STACK_SIZE];
+            let mut stack_t = [MaybeUninit::<f64>::uninit(); STACK_SIZE];
+            stack_ref[0].write(root);
+            stack_t[0].write(t_min);
+            let mut sp = 1;
+            let mut tnear = [0.0; WIDTH];
+
+            while sp > 0 {
+                sp -= 1;
+                if stack_t[sp].assume_init() > t_max {
+                    continue; // Already found something closer than this whole subtree.
+                }
+                let node_ref = stack_ref[sp].assume_init();
+                if node_ref & LEAF_BIT != 0 {
+                    let (start, count) = decode_leaf(node_ref);
+                    if node_ref & HAS_SPHERES != 0 {
+                        let (mut mask, roots) = $spheres(&self.spheres, start, count, r, t_min, t_max);
+                        while mask != 0 {
+                            let i = mask.trailing_zeros() as usize;
+                            mask &= mask - 1;
+                            if roots[i] < t_max {
+                                t_max = roots[i];
+                                best = Some((self.slots[start + i] as usize, roots[i], true));
+                            }
+                        }
+                    }
+                    if node_ref & HAS_OTHERS != 0 {
+                        for k in start..start + count {
+                            if self.spheres.r2[k] == f64::NEG_INFINITY {
+                                let obj = self.slots[k] as usize;
+                                if let Some(t) = other(obj, t_max) {
+                                    t_max = t;
+                                    best = Some((obj, t, false));
+                                }
+                            }
+                        }
+                    }
+                    continue;
+                }
+
+                let node = &self.nodes[node_ref as usize];
+                let mut mask = $node(node, r, t_min, t_max, &mut tnear);
+                // Push hit children farthest first, so the nearest is popped
+                // next and shrinks t_max soonest.
+                let mut hits = [(0u32, 0.0f64); WIDTH];
+                let mut n = 0;
+                while mask != 0 {
+                    let i = mask.trailing_zeros() as usize;
+                    mask &= mask - 1;
+                    let mut j = n;
+                    while j > 0 && hits[j - 1].1 < tnear[i] {
+                        hits[j] = hits[j - 1];
+                        j -= 1;
+                    }
+                    hits[j] = (node.child[i], tnear[i]);
+                    n += 1;
+                }
+                debug_assert!(sp + n <= STACK_SIZE);
+                for &(child, t) in &hits[..n] {
+                    stack_ref[sp].write(child);
+                    stack_t[sp].write(t);
+                    sp += 1;
+                }
+            }
+            best
+        }
+
+        /// True as soon as anything is hit within (t_min, t_max);
+        /// `other(object)` tests non-sphere objects.
+        $(#[$attr])*
+        unsafe fn $any(
+            &self,
+            root: u32,
+            r: &RayData,
+            t_min: f64,
+            t_max: f64,
+            mut other: impl FnMut(usize) -> bool,
+        ) -> bool {
+            let mut stack = [MaybeUninit::<u32>::uninit(); STACK_SIZE];
+            stack[0].write(root);
+            let mut sp = 1;
+            let mut tnear = [0.0; WIDTH];
+
+            while sp > 0 {
+                sp -= 1;
+                let node_ref = stack[sp].assume_init();
+                if node_ref & LEAF_BIT != 0 {
+                    let (start, count) = decode_leaf(node_ref);
+                    if node_ref & HAS_SPHERES != 0 && $spheres(&self.spheres, start, count, r, t_min, t_max).0 != 0 {
+                        return true;
+                    }
+                    if node_ref & HAS_OTHERS != 0 {
+                        for k in start..start + count {
+                            if self.spheres.r2[k] == f64::NEG_INFINITY && other(self.slots[k] as usize) {
+                                return true;
+                            }
+                        }
+                    }
+                    continue;
+                }
+
+                let node = &self.nodes[node_ref as usize];
+                let mut mask = $node(node, r, t_min, t_max, &mut tnear);
+                while mask != 0 {
+                    let i = mask.trailing_zeros() as usize;
+                    mask &= mask - 1;
+                    debug_assert!(sp < STACK_SIZE);
+                    stack[sp].write(node.child[i]);
+                    sp += 1;
+                }
+            }
+            false
+        }
+    };
+}
+
+pub struct Bvh {
+    /// Reference to the root: node 0, or a single leaf when the whole scene
+    /// fits in one (then there is no node test at all). None if empty.
+    root: Option<u32>,
+    /// Bounds of everything; tested first when the root is a single leaf,
+    /// which otherwise has no box test to reject rays that miss it.
+    root_bounds: Aabb,
+    nodes: Vec<Node8>,
+    /// Object index for each slot; each leaf covers a contiguous slot range.
+    slots: Vec<u32>,
+    spheres: Spheres,
+    use_avx512: bool,
+}
+
+impl Bvh {
+    /// `spheres[i]` is Some((center, radius)) if object i is a sphere.
+    pub fn build(bounds: &[Aabb], spheres: &[Option<(Vec3, f64)>]) -> Bvh {
+        let mut slots: Vec<u32> = (0..bounds.len() as u32).collect();
+        let mut bin = Vec::with_capacity(2 * bounds.len());
+        if !bounds.is_empty() {
+            let input = BuildInput {
+                bounds,
+                centroids: bounds.iter().map(Aabb::centroid).collect(),
+                is_sphere: spheres.iter().map(Option::is_some).collect(),
+            };
+            build_binary(&mut bin, &mut slots, 0, &input);
+        }
+
+        let mut soa = Spheres {
+            cx: vec![0.0; slots.len() + WIDTH],
+            cy: vec![0.0; slots.len() + WIDTH],
+            cz: vec![0.0; slots.len() + WIDTH],
+            r2: vec![f64::NEG_INFINITY; slots.len() + WIDTH],
+        };
+        for (k, &obj) in slots.iter().enumerate() {
+            if let Some((c, radius)) = spheres[obj as usize] {
+                (soa.cx[k], soa.cy[k], soa.cz[k], soa.r2[k]) = (c.x, c.y, c.z, radius * radius);
+            }
+        }
+
+        let mut bvh = Bvh {
+            root: None,
+            root_bounds: bin.first().map_or(Aabb::EMPTY, |b| b.bounds),
+            nodes: Vec::new(),
+            slots,
+            spheres: soa,
+            use_avx512: detect_avx512(),
+        };
+        if let Some(b) = bin.first() {
+            bvh.root = Some(if b.count > 0 { bvh.leaf(b) } else { bvh.collapse(&bin, 0) as u32 });
+        }
+        assert!(
+            (WIDTH - 1) * bvh.depth() + 1 <= STACK_SIZE,
+            "BVH too deep for the traversal stack"
+        );
+        bvh
+    }
+
+    fn leaf(&self, b: &BinNode) -> u32 {
+        let r2 = &self.spheres.r2[b.start..b.start + b.count];
+        let has_spheres = r2.iter().any(|&r| r != f64::NEG_INFINITY);
+        let has_others = r2.iter().any(|&r| r == f64::NEG_INFINITY);
+        leaf_ref(b.start, b.count, has_spheres, has_others)
+    }
+
+    /// Converts binary node `b` (and its subtree) into 8-wide nodes, returning
+    /// the new node's index. Starts from b's two children and repeatedly
+    /// opens up the interior child with the largest surface area until there
+    /// are 8 children or only leaves remain.
+    fn collapse(&mut self, bin: &[BinNode], b: usize) -> usize {
+        let mut children = vec![bin[b].left, bin[b].right];
+        while children.len() < WIDTH {
+            let Some(pos) = children
+                .iter()
+                .enumerate()
+                .filter(|(_, &c)| bin[c].count == 0)
+                .max_by(|(_, &a), (_, &c)| bin[a].bounds.surface_area().total_cmp(&bin[c].bounds.surface_area()))
+                .map(|(pos, _)| pos)
+            else {
+                break;
+            };
+            let opened = children.swap_remove(pos);
+            children.push(bin[opened].left);
+            children.push(bin[opened].right);
+        }
+
+        let idx = self.nodes.len();
+        self.nodes.push(Node8::empty());
+        for (lane, &c) in children.iter().enumerate() {
+            let child_ref = if bin[c].count > 0 {
+                self.leaf(&bin[c])
+            } else {
+                self.collapse(bin, c) as u32
+            };
+            let (lo, hi) = (bin[c].bounds.min, bin[c].bounds.max);
+            let n = &mut self.nodes[idx];
+            (n.min_x[lane], n.min_y[lane], n.min_z[lane]) = (lo.x, lo.y, lo.z);
+            (n.max_x[lane], n.max_y[lane], n.max_z[lane]) = (hi.x, hi.y, hi.z);
+            n.child[lane] = child_ref;
+        }
+        idx
+    }
+
+    /// Number of 8-wide node levels (0 when the root is a single leaf).
+    pub fn depth(&self) -> usize {
+        fn go(b: &Bvh, n: usize) -> usize {
+            1 + b.nodes[n]
+                .child
+                .iter()
+                .zip(&b.nodes[n].min_x)
+                .filter(|(&c, &mx)| mx != f64::INFINITY && c & LEAF_BIT == 0)
+                .map(|(&c, _)| go(b, c as usize))
+                .max()
+                .unwrap_or(0)
+        }
+        if self.nodes.is_empty() {
+            0
+        } else {
+            go(self, 0)
+        }
+    }
+
+    /// False if the root is a single leaf whose bounds the ray misses. (A
+    /// node root does this test for its children itself.)
+    #[inline(always)]
+    fn root_hit(&self, root: u32, r: &RayData, t_min: f64, t_max: f64) -> bool {
+        if root & LEAF_BIT == 0 {
+            return true;
+        }
+        let (lo, hi) = (self.root_bounds.min, self.root_bounds.max);
+        let (mut near, mut far) = (t_min, t_max);
+        for axis in 0..3 {
+            let t0 = (lo[axis] - r.o[axis]) * r.inv[axis];
+            let t1 = (hi[axis] - r.o[axis]) * r.inv[axis];
+            near = near.max(t0.min(t1));
+            far = far.min(t0.max(t1));
+        }
+        near <= far
+    }
+
+    /// Whether traversal runs the AVX-512 kernels. When the whole scene is a
+    /// single leaf without spheres there is nothing for them to do, and code
+    /// compiled for AVX-512 measured slightly slower on such scenes.
+    fn simd_pays_off(&self, root: u32) -> bool {
+        self.use_avx512 && (root & LEAF_BIT == 0 || root & HAS_SPHERES != 0)
+    }
+
+    pub fn uses_avx512(&self) -> bool {
+        self.root.is_some_and(|root| self.simd_pays_off(root))
+    }
+
+    traversal!(closest_scalar, any_scalar, node_scalar, spheres_scalar);
+
+    #[cfg(target_arch = "x86_64")]
+    traversal!(
+        #[target_feature(enable = "avx512f")]
+        closest_avx512,
+        any_avx512,
+        avx512::node,
+        avx512::spheres
+    );
+
+    /// See the traversal macro. Returns (object, t, is_sphere) for the closest hit.
+    pub fn closest(
+        &self,
+        r: &Ray,
+        t_min: f64,
+        t_max: f64,
+        other: impl FnMut(usize, f64) -> Option<f64>,
+    ) -> Option<(usize, f64, bool)> {
+        let Some(root) = self.root else {
+            return None;
+        };
+        let rd = RayData::new(r);
+        let t_max = t_max.min(f64::MAX); // See Node8.
+        if !self.root_hit(root, &rd, t_min, t_max) {
+            return None;
+        }
+        #[cfg(target_arch = "x86_64")]
+        if self.simd_pays_off(root) {
+            // SAFETY: simd_pays_off() requires use_avx512, which is only set when the CPU supports AVX-512F.
+            return unsafe { self.closest_avx512(root, &rd, t_min, t_max, other) };
+        }
+        unsafe { self.closest_scalar(root, &rd, t_min, t_max, other) }
+    }
+
+    pub fn any(&self, r: &Ray, t_min: f64, t_max: f64, other: impl FnMut(usize) -> bool) -> bool {
+        let Some(root) = self.root else {
+            return false;
+        };
+        let rd = RayData::new(r);
+        let t_max = t_max.min(f64::MAX); // See Node8.
+        if !self.root_hit(root, &rd, t_min, t_max) {
+            return false;
+        }
+        #[cfg(target_arch = "x86_64")]
+        if self.simd_pays_off(root) {
+            // SAFETY: as in closest().
+            return unsafe { self.any_avx512(root, &rd, t_min, t_max, other) };
+        }
+        unsafe { self.any_scalar(root, &rd, t_min, t_max, other) }
+    }
+}
+
+fn detect_avx512() -> bool {
+    if std::env::var_os("PATHTRACER_NO_AVX512").is_some() {
+        return false;
+    }
+    #[cfg(target_arch = "x86_64")]
+    {
+        std::arch::is_x86_feature_detected!("avx512f")
+    }
+    #[cfg(not(target_arch = "x86_64"))]
+    {
+        false
+    }
+}

rust/src/camera.rs

@@ -0,0 +1,33 @@
+use crate::scene::Ray;
+use crate::vec3::Vec3;
+
+pub struct Camera {
+    origin: Vec3,
+    lower_left_corner: Vec3,
+    horizontal: Vec3,
+    vertical: Vec3,
+}
+
+impl Camera {
+    pub fn new(lookfrom: Vec3, lookat: Vec3, vup: Vec3, vfov: f64, aspect: f64, focus_dist: f64) -> Self {
+        let half_height = (vfov.to_radians() / 2.0).tan();
+        let half_width = aspect * half_height;
+        let w = (lookfrom - lookat).normalize();
+        let u = vup.cross(w).normalize();
+        let v = w.cross(u);
+        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,
+            lower_left_corner,
+            horizontal,
+            vertical,
+        }
+    }
+
+    pub fn get_ray(&self, s: f64, t: f64) -> Ray {
+        let direction = self.lower_left_corner + s * self.horizontal + t * self.vertical - self.origin;
+        Ray::new(self.origin, direction)
+    }
+}

rust/src/main.rs

@@ -0,0 +1,691 @@
+mod bvh;
+mod camera;
+mod medium;
+mod photons;
+mod scene;
+mod scenes;
+mod vec3;
+
+use std::f64::consts::PI;
+use std::time::Instant;
+
+use indicatif::{ProgressBar, ProgressStyle};
+use rand::rngs::SmallRng;
+use rand::{Rng, SeedableRng};
+use rayon::prelude::*;
+
+use camera::Camera;
+use medium::{hg_phase, sample_hg, MediumState};
+use photons::CausticMaps;
+use scene::{spot_emission, Light, Material, Ray, Scene, EPSILON};
+use vec3::Vec3;
+
+/// 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;
+
+/// Where light sampled by next event estimation arrives.
+enum Receiver {
+    /// A diffuse surface: Lambertian BRDF times the cosine at the surface.
+    Surface { normal: Vec3, albedo: Vec3 },
+    /// A scattering point in a medium: the phase function for turning from
+    /// the path's direction `dir` towards the light.
+    Medium { dir: Vec3, g: f64 },
+}
+
+/// Direct light at `point` from one uniformly sampled point on `light`'s
+/// surface, attenuated by whatever media the shadow ray passes through.
+fn sample_direct<R: Rng>(
+    scene: &Scene,
+    state: MediumState,
+    light: &Light,
+    point: Vec3,
+    receiver: Receiver,
+    rng: &mut R,
+) -> Vec3 {
+    let (lp, light_normal) = light.sample(rng);
+
+    let to_light = lp - point;
+    let dist_sq = to_light.dot(to_light);
+    if dist_sq < 1e-12 {
+        return Vec3::ZERO;
+    }
+    let dist = dist_sq.sqrt();
+    let d_to_light = to_light / dist;
+    let cos_at_light = -light_normal.dot(d_to_light);
+    let (response, cos_at_rec, origin) = match receiver {
+        Receiver::Surface { normal, albedo } => (albedo / PI, normal.dot(d_to_light), point + EPSILON * normal),
+        Receiver::Medium { dir, g } => (Vec3::splat(hg_phase(dir.dot(d_to_light), g)), 1.0, point),
+    };
+    if cos_at_rec <= 0.0 || cos_at_light <= 0.0 {
+        return Vec3::ZERO;
+    }
+
+    let tr = scene.transmittance(state, origin, d_to_light, dist * (1.0 - SHADOW_RELATIVE_EPSILON));
+    if tr == Vec3::ZERO {
+        return Vec3::ZERO;
+    }
+
+    // Uniform area sampling has pdf 1 / area.
+    let geometry = cos_at_rec * cos_at_light / dist_sq;
+    let emission = spot_emission(light.emission, light.spot, cos_at_light);
+    (emission.mul_elem(response) * (geometry * light.area())).mul_elem(tr)
+}
+
+/// Bounces before Russian roulette may terminate a path.
+const RR_MIN_DEPTH: u32 = 3;
+/// Guard against paths that keep crossing invisible volume boundaries.
+const MAX_VOLUME_CROSSINGS: u32 = 256;
+
+/// Russian roulette: after RR_MIN_DEPTH bounces, end dim paths early,
+/// boosting survivors so the estimate stays unbiased. Glass/mirror paths
+/// (throughput ~1) survive. Returns false if the path ends.
+fn russian_roulette<R: Rng>(throughput: &mut Vec3, bounces: u32, rng: &mut R) -> bool {
+    if bounces >= RR_MIN_DEPTH {
+        let survive = throughput.x.max(throughput.y).max(throughput.z).min(0.95);
+        if rng.gen::<f64>() >= survive {
+            return false;
+        }
+        *throughput = *throughput / survive;
+    }
+    true
+}
+
+/// With `caustics`, caustic light (light that went through glass or off
+/// mirrors before reaching fog or a diffuse surface) comes from the photon
+/// maps, and path tracing skips the matching paths (from a diffuse or medium
+/// interaction through mirrors/glass to a light):
+/// - normally (all lights emit photons), only at the path's first such
+///   interaction: the maps are gathered along fog segments and at the first
+///   diffuse surface until then, and later caustics are path traced;
+/// - with `direct_in_media` (only spotlights emit photons), for spotlights
+///   at every interaction: the volume map also holds their direct light, so
+///   points in fog don't sample spotlights, and every fog segment and
+///   diffuse surface of the path gathers the maps. Other lights are path
+///   traced as usual.
+fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R, caustics: Option<&CausticMaps>) -> 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 or medium scattering: 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;
+    let mut state = MediumState::Outside;
+    let mut bounces = 0;
+    let mut crossings = 0;
+    // Diffuse and medium interactions so far.
+    let mut non_specular = 0;
+    let direct_in_media = caustics.is_some_and(|m| m.direct_in_media);
+    // Whether this interaction gathers from the photon maps.
+    let maps_here = |non_specular: u32| caustics.filter(|_| direct_in_media || non_specular == 0);
+
+    while bounces < max_depth {
+        let hit = scene.hit(&ray, EPSILON, f64::INFINITY);
+
+        // Travel through the medium (if any) up to the surface: either
+        // scatter somewhere along the way, or arrive with the absorption
+        // along the way applied.
+        let t_hit = hit.as_ref().map_or(f64::INFINITY, |h| h.t);
+        if let Some((medium, t0, t1)) = scene.medium_segment(state, &ray, t_hit) {
+            if let Some(maps) = maps_here(non_specular) {
+                radiance += throughput.mul_elem(maps.beam(&medium, &ray, t0, t1));
+            }
+            let t = t0 + medium.sample_distance(rng);
+            if t < t1 {
+                throughput = throughput.mul_elem(medium.absorption(t - t0));
+                let p = ray.at(t);
+                // Spotlights are left out when the photon map carries their
+                // direct light into the fog.
+                if let Some((light, pick_prob)) = scene.pick_light(p, direct_in_media, rng) {
+                    let receiver = Receiver::Medium { dir: ray.direction, g: medium.g };
+                    let direct = sample_direct(scene, state, light, p, receiver, rng) / pick_prob;
+                    radiance += throughput.mul_elem(direct);
+                }
+                count_emission = false;
+                non_specular += 1;
+                ray = Ray::from_unit(p, sample_hg(ray.direction, medium.g, rng));
+                bounces += 1;
+                if !russian_roulette(&mut throughput, bounces, rng) {
+                    break;
+                }
+                continue;
+            }
+            throughput = throughput.mul_elem(medium.absorption(t1 - t0));
+        }
+
+        let Some(rec) = hit else {
+            radiance += throughput.mul_elem(scene.background.radiance(ray.direction));
+            break;
+        };
+
+        match rec.material {
+            Material::Light { emission, spot } => {
+                // Lights emit from their front side only (matters for one-sided rects).
+                // Reached from a diffuse/medium interaction via mirrors or
+                // glass: a caustic, which the photon maps may already cover.
+                let from_maps = match caustics {
+                    Some(_) if non_specular == 0 => false,
+                    Some(m) if m.direct_in_media => spot.is_some(),
+                    Some(_) => non_specular == 1,
+                    None => false,
+                };
+                if count_emission && rec.front_face && !from_maps {
+                    let emission = spot_emission(emission, spot, -ray.direction.dot(rec.normal));
+                    radiance += throughput.mul_elem(emission);
+                }
+                break;
+            }
+            Material::Volume { medium } => {
+                // Invisible boundary: carry on in the same direction, now
+                // inside (or back outside) the volume. Not a bounce.
+                state = if rec.front_face { MediumState::Inside(Some(medium)) } else { MediumState::Outside };
+                ray = Ray::from_unit(rec.point, ray.direction);
+                crossings += 1;
+                if crossings > MAX_VOLUME_CROSSINGS {
+                    break;
+                }
+                continue;
+            }
+            _ => {}
+        }
+
+        let Some((attenuation, scattered)) = rec.material.scatter(&ray, &rec, rng) else {
+            break;
+        };
+
+        match rec.material {
+            Material::Diffuse { albedo } => {
+                if let Some(maps) = maps_here(non_specular) {
+                    radiance += throughput.mul_elem(maps.surface(rec.point, rec.normal, albedo));
+                }
+                if let Some((light, pick_prob)) = scene.pick_light(rec.point, false, rng) {
+                    let receiver = Receiver::Surface { normal: rec.normal, albedo };
+                    let direct = sample_direct(scene, state, light, rec.point, receiver, rng) / pick_prob;
+                    radiance += throughput.mul_elem(direct);
+                }
+                count_emission = false;
+                non_specular += 1;
+            }
+            Material::Dielectric { medium, .. } => {
+                // Refracted through the surface (rather than reflected):
+                // entering or leaving the object's interior medium.
+                if scattered.direction.dot(rec.normal) < 0.0 {
+                    state = if rec.front_face { MediumState::Inside(medium) } else { MediumState::Outside };
+                }
+                count_emission = true;
+            }
+            _ => count_emission = true,
+        }
+
+        throughput = throughput.mul_elem(attenuation);
+        ray = scattered;
+        bounces += 1;
+        if !russian_roulette(&mut throughput, bounces, rng) {
+            break;
+        }
+    }
+
+    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,
+}
+
+/// 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;
+
+/// 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,
+    caustics: Option<&'a CausticMaps>,
+}
+
+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, self.caustics);
+            px.sum += sample;
+            px.stats.add(sample.luminance());
+        }
+    }
+
+    /// 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]);
+                    }
+                }
+                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, caustics: Option<&CausticMaps>) -> (Vec<Vec3>, u64) {
+    let (width, height) = (s.width, s.height);
+    let view = scene.view;
+    let r = Renderer {
+        scene,
+        cam: Camera::new(
+            view.lookfrom,
+            view.lookat,
+            Vec3::new(0.0, 1.0, 0.0),
+            view.vfov,
+            width as f64 / height as f64,
+            (view.lookfrom - view.lookat).length(),
+        ),
+        s,
+        caustics,
+    };
+
+    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();
+
+    // 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<()> {
+    let bytes: Vec<u8> = image
+        .iter()
+        .flat_map(|c| [c.x, c.y, c.z])
+        .map(|c| (c.clamp(0.0, 1.0) * 255.0) as u8)
+        .collect();
+    image::save_buffer(path, &bytes, width as u32, height as u32, image::ColorType::Rgb8)
+}
+
+fn region_mean(image: &[Vec3], width: usize, rows: std::ops::Range<usize>, cols: std::ops::Range<usize>) -> f64 {
+    let n = rows.len() * cols.len();
+    if n == 0 {
+        return 0.0;
+    }
+    let sum: f64 = rows
+        .flat_map(|j| cols.clone().map(move |i| image[j * width + i].mean()))
+        .sum();
+    sum / n as f64
+}
+
+fn print_region_table(image: &[Vec3], width: usize, height: usize) {
+    let (h_start, h_end) = (height / 3, 2 * height / 3);
+    let (w_start, w_end) = (width / 3, 2 * width / 3);
+    let regions = [
+        ("Left Wall", 0..height, 0..width / 4),
+        ("Right Wall", 0..height, 3 * width / 4..width),
+        ("Floor", h_end..height, w_start..w_end),
+        ("Ceiling", 0..h_start, w_start..w_end),
+        ("Center", h_start..h_end, w_start..w_end),
+    ];
+
+    println!("Region Brightness Analysis");
+    println!("{:<12} {:>10}  Visual", "Region", "Brightness");
+    for (name, rows, cols) in regions {
+        let brightness = region_mean(image, width, rows, cols) * 255.0;
+        let filled = ((brightness / 255.0 * 20.0) as usize).min(20);
+        let bar = format!("{}{}", "#".repeat(filled), "-".repeat(20 - filled));
+        let color = if brightness > 80.0 { 32 } else if brightness > 40.0 { 33 } else { 31 };
+        println!("{name:<12} {brightness:>10.1}  \x1b[{color}m{bar}\x1b[0m");
+    }
+}
+
+const USAGE: &str = "usage: pathtracer [--scene NAME] [--list-scenes] [--photons N | --no-photons] \
+[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() {
+    // Pull out the --scene / --list-scenes options; the rest are positional.
+    let mut scene_name = String::from("cornell");
+    // Override for the number of caustic photons (0: path trace only).
+    let mut photons_override: Option<usize> = None;
+    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 == "--photons" {
+            photons_override = raw.next().and_then(|n| n.parse().ok());
+            if photons_override.is_none() {
+                eprintln!("{USAGE}");
+                std::process::exit(2);
+            }
+        } else if a == "--no-photons" {
+            photons_override = Some(0);
+        } 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}");
+            std::process::exit(2);
+        })
+    }
+
+    let settings = Settings {
+        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);
+    }
+
+    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);
+    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!(
+        "BVH:               8-wide, depth {}, {}",
+        scene.bvh_depth(),
+        if scene.uses_avx512() { "AVX-512" } else { "scalar" }
+    );
+    println!("Output:            {}", settings.output);
+    println!();
+
+    let start = Instant::now();
+    let photon_settings = scene.caustic_photons.map(|mut p| {
+        if let Some(n) = photons_override {
+            p.emitted = n;
+        }
+        p
+    });
+    let caustics = match photon_settings {
+        Some(p) if p.emitted > 0 => {
+            let t = Instant::now();
+            let maps = CausticMaps::build(&scene, p);
+            let (volume, surface) = maps.counts();
+            println!(
+                "Caustic photons:   {} emitted, {} stored in media, {} on surfaces ({:.2}s)",
+                p.emitted,
+                volume,
+                surface,
+                t.elapsed().as_secs_f64()
+            );
+            Some(maps)
+        }
+        _ => None,
+    };
+    let (image, total_samples) = render(&scene, &settings, caustics.as_ref());
+    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);
+        std::process::exit(1);
+    }
+
+    println!();
+    println!("Rendered {}", settings.output);
+    println!(
+        "Total time: {:.2}s | Speed: {:.2} MPix*samp/s | Avg samples/pixel: {:.1}",
+        elapsed,
+        total_samples as f64 / elapsed / 1e6,
+        total_samples as f64 / pixels
+    );
+    // 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);
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use medium::Medium;
+    use scene::{Background, SceneBuilder, View};
+
+    fn view() -> View {
+        View { lookfrom: Vec3::new(0.0, 0.0, -10.0), lookat: Vec3::ZERO, vfov: 40.0 }
+    }
+
+    /// Average of `n` path samples along one ray.
+    fn mean_radiance(scene: &Scene, ray: Ray, n: usize) -> Vec3 {
+        let mut rng = SmallRng::seed_from_u64(1);
+        let mut sum = Vec3::ZERO;
+        for _ in 0..n {
+            sum += trace(scene, &ray, 10_000, &mut rng, None);
+        }
+        sum / n as f64
+    }
+
+    fn assert_close(v: Vec3, expected: Vec3, tol: f64) {
+        for a in 0..3 {
+            assert!((v[a] - expected[a]).abs() < tol, "got {v:?}, expected {expected:?}");
+        }
+    }
+
+    const INTO_CENTER: fn() -> Ray = || Ray::from_unit(Vec3::new(0.0, 0.0, -10.0), Vec3::new(0.0, 0.0, 1.0));
+
+    /// White furnace: with no absorption and a uniform background of 1,
+    /// every path eventually leaves and sees 1, so the mean must be 1 however
+    /// much the medium scatters.
+    #[test]
+    fn furnace_fog() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        b.fog(Vec3::splat(-1.0), Vec3::splat(1.0), Medium::new(Vec3::ZERO, 3.0, 0.6));
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 200_000), Vec3::ONE, 0.01);
+    }
+
+    #[test]
+    fn furnace_volume_sphere() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        let smoke = Medium::new(Vec3::ZERO, 4.0, -0.3);
+        b.sphere(Vec3::ZERO, 1.0, Material::Volume { medium: smoke });
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 200_000), Vec3::ONE, 0.01);
+    }
+
+    #[test]
+    fn furnace_glass_with_scattering_interior() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        let milk = Medium::new(Vec3::ZERO, 5.0, 0.2);
+        b.sphere(Vec3::ZERO, 1.0, Material::Dielectric { ref_idx: 1.5, albedo: Vec3::ONE, medium: Some(milk) });
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 200_000), Vec3::ONE, 0.01);
+    }
+
+    /// Henyey-Greenstein: sampled directions have mean cosine g, and the
+    /// phase function (used by next event estimation) integrates to 1.
+    #[test]
+    fn henyey_greenstein() {
+        let mut rng = SmallRng::seed_from_u64(2);
+        let dir = Vec3::new(0.3, -0.5, 0.8).normalize();
+        for g in [-0.7, 0.0, 0.4, 0.9] {
+            let n = 400_000;
+            let mean_cos: f64 = (0..n).map(|_| sample_hg(dir, g, &mut rng).dot(dir)).sum::<f64>() / n as f64;
+            assert!((mean_cos - g).abs() < 0.005, "g = {g}: mean cos {mean_cos}");
+
+            // Midpoint rule over cos_theta (the phase function doesn't depend on phi).
+            let steps = 200_000;
+            let integral: f64 = (0..steps)
+                .map(|i| {
+                    let c = -1.0 + (i as f64 + 0.5) * 2.0 / steps as f64;
+                    hg_phase(c, g) * 2.0 * PI * (2.0 / steps as f64)
+                })
+                .sum();
+            assert!((integral - 1.0).abs() < 1e-3, "g = {g}: integral {integral}");
+        }
+    }
+
+    /// A purely absorbing slab of thickness 2: Beer's law, exactly, per channel.
+    #[test]
+    fn beer_lambert_absorption() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        let sigma_a = Vec3::new(0.1, 0.5, 1.5);
+        b.fog(Vec3::splat(-1.0), Vec3::splat(1.0), Medium::new(sigma_a, 0.0, 0.0));
+        let expected = (sigma_a * -2.0).map(f64::exp);
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 10), expected, 1e-12);
+    }
+}

rust/src/medium.rs

@@ -0,0 +1,109 @@
+//! Participating media: homogeneous volumes that absorb and scatter light
+//! along a ray, not just at surfaces.
+
+use rand::Rng;
+
+use crate::bvh::Aabb;
+use crate::vec3::Vec3;
+
+/// A homogeneous medium. Coefficients are per scene unit of distance.
+///
+/// Scattering is grey so that distances can be sampled with one scalar
+/// rate; absorption may be coloured and is applied as an exact per-channel
+/// weight exp(-sigma_a * d) rather than sampled. (Sampling scattering at
+/// rate sigma_s and weighting by exp(-sigma_a d) gives exactly the
+/// transmittance / pdf ratio of the full extinction, so this stays unbiased.)
+#[derive(Clone, Copy, Debug)]
+pub struct Medium {
+    pub sigma_a: Vec3,
+    pub sigma_s: f64,
+    /// Henyey-Greenstein anisotropy: 0 scatters evenly in all directions,
+    /// towards +1 mostly forwards, towards -1 mostly backwards.
+    pub g: f64,
+}
+
+impl Medium {
+    pub fn new(sigma_a: Vec3, sigma_s: f64, g: f64) -> Self {
+        Medium { sigma_a, sigma_s, g }
+    }
+
+    /// Fraction of light that survives absorption over distance `d`.
+    pub fn absorption(&self, d: f64) -> Vec3 {
+        (self.sigma_a * -d).map(f64::exp)
+    }
+
+    /// Fraction of light that survives both absorption and out-scattering
+    /// over distance `d` (what a shadow ray sees).
+    pub fn transmittance(&self, d: f64) -> Vec3 {
+        self.absorption(d) * (-self.sigma_s * d).exp()
+    }
+
+    /// Distance to the next scattering event, sampled with pdf
+    /// sigma_s * exp(-sigma_s * t). Infinite if the medium does not scatter.
+    pub fn sample_distance<R: Rng>(&self, rng: &mut R) -> f64 {
+        if self.sigma_s > 0.0 {
+            -(1.0 - rng.gen::<f64>()).ln() / self.sigma_s
+        } else {
+            f64::INFINITY
+        }
+    }
+}
+
+/// The Henyey-Greenstein phase function for the angle between the incoming
+/// propagation direction and the outgoing one (cos_theta = 1: straight on).
+pub fn hg_phase(cos_theta: f64, g: f64) -> f64 {
+    let denom = 1.0 + g * g - 2.0 * g * cos_theta;
+    (1.0 - g * g) / (4.0 * std::f64::consts::PI * denom * denom.sqrt())
+}
+
+/// A direction distributed by the Henyey-Greenstein phase function around
+/// the propagation direction `dir` (unit length). The weight phase / pdf is
+/// exactly 1, so the caller only multiplies by the medium's albedo.
+pub fn sample_hg<R: Rng>(dir: Vec3, g: f64, rng: &mut R) -> Vec3 {
+    let (u1, u2): (f64, f64) = (rng.gen(), rng.gen());
+    let cos_theta = if g.abs() < 1e-3 {
+        1.0 - 2.0 * u1
+    } else {
+        let s = (1.0 - g * g) / (1.0 - g + 2.0 * g * u1);
+        ((1.0 + g * g - s * s) / (2.0 * g)).clamp(-1.0, 1.0)
+    };
+    let sin_theta = (1.0 - cos_theta * cos_theta).max(0.0).sqrt();
+    let phi = 2.0 * std::f64::consts::PI * u2;
+    let (t, b) = crate::scene::orthonormal_basis(dir);
+    t * (sin_theta * phi.cos()) + b * (sin_theta * phi.sin()) + dir * cos_theta
+}
+
+/// A medium filling an axis-aligned region of the scene outside all objects,
+/// e.g. a layer of haze below some height. Bounds may be infinite.
+#[derive(Clone, Copy, Debug)]
+pub struct Fog {
+    pub region: Aabb,
+    pub medium: Medium,
+}
+
+impl Fog {
+    /// The part of the ray segment [t_min, t_max] that lies inside the region.
+    pub fn clip(&self, origin: Vec3, dir: Vec3, t_min: f64, t_max: f64) -> Option<(f64, f64)> {
+        let (mut t0, mut t1) = (t_min, t_max);
+        for axis in 0..3 {
+            // A finite stand-in for 1/0 keeps (bound - origin) * inv free of NaNs.
+            let inv = 1.0 / dir[axis];
+            let inv = if inv.is_finite() { inv } else { 1e300_f64.copysign(inv) };
+            let a = (self.region.min[axis] - origin[axis]) * inv;
+            let b = (self.region.max[axis] - origin[axis]) * inv;
+            t0 = t0.max(a.min(b));
+            t1 = t1.min(a.max(b));
+        }
+        (t0 < t1).then_some((t0, t1))
+    }
+}
+
+/// Which medium a path is currently travelling through.
+#[derive(Clone, Copy, Debug)]
+pub enum MediumState {
+    /// Outside all objects: in the scene's fog, if any, where the fog region reaches.
+    Outside,
+    /// Inside an object, with that object's interior medium (None: clear,
+    /// like plain glass).
+    Inside(Option<Medium>),
+}

rust/src/photons.rs

@@ -0,0 +1,357 @@
+//! Caustic photon maps (after Jensen's photon mapping, with the beam
+//! radiance estimate of Jarosz et al. 2008 for media).
+//!
+//! Path tracing finds caustics by starting at the camera: a point in fog or on
+//! a surface has to pick a direction that happens to go through the glass and
+//! land on the light. With a small light or a spotlight that almost never
+//! happens. Photons go the other way: they are shot from the lights, and
+//! those that pass through glass or off mirrors ("caustic photons") are
+//! stored where they first hit fog or a diffuse surface. The renderer then
+//! gathers nearby photons (a) along camera rays through fog, (b) at the first
+//! diffuse surface a camera path reaches, and path tracing skips exactly
+//! those caustic paths so nothing is counted twice (see trace() in main.rs).
+//!
+//! Gathering blurs by the radius (bias), which shrinks as the radius does.
+
+use rand::rngs::SmallRng;
+use rand::SeedableRng;
+use rayon::prelude::*;
+
+use crate::medium::{hg_phase, Medium, MediumState};
+use crate::scene::{CausticPhotons, Material, Ray, Scene, EPSILON};
+use crate::vec3::Vec3;
+
+#[derive(Clone, Copy, Debug)]
+pub struct Photon {
+    pub pos: Vec3,
+    /// Direction of travel when stored.
+    pub dir: Vec3,
+    pub power: Vec3,
+}
+
+/// Longest specular chain a photon may follow.
+const MAX_PHOTON_BOUNCES: u32 = 32;
+
+/// Traces `emitted` photons and returns the photons stored in media and on
+/// surfaces: caustic photons, plus (with `direct_in_media`) every spotlight
+/// photon's first scattering in a medium.
+fn trace_photons(scene: &Scene, emitted: usize, direct_in_media: bool) -> (Vec<Photon>, Vec<Photon>) {
+    const CHUNK: usize = 1 << 16;
+    let chunks = emitted.div_ceil(CHUNK);
+    let per_photon = 1.0 / emitted as f64;
+    let results: Vec<(Vec<Photon>, Vec<Photon>)> = (0..chunks)
+        .into_par_iter()
+        .map(|c| {
+            let mut rng = SmallRng::seed_from_u64(0x5EED_0000 + c as u64);
+            let (mut volume, mut surface) = (Vec::new(), Vec::new());
+            let count = CHUNK.min(emitted - c * CHUNK);
+            for _ in 0..count {
+                trace_photon(scene, per_photon, direct_in_media, &mut rng, &mut volume, &mut surface);
+            }
+            (volume, surface)
+        })
+        .collect();
+    let mut volume = Vec::new();
+    let mut surface = Vec::new();
+    for (v, s) in results {
+        volume.extend(v);
+        surface.extend(s);
+    }
+    (volume, surface)
+}
+
+fn trace_photon(
+    scene: &Scene,
+    per_photon: f64,
+    direct_in_media: bool,
+    rng: &mut SmallRng,
+    volume: &mut Vec<Photon>,
+    surface: &mut Vec<Photon>,
+) {
+    let Some((light, pick_prob)) = scene.pick_emitting_light(direct_in_media, rng) else {
+        return;
+    };
+    let (origin, dir, power) = light.emit(rng);
+    let store_direct = direct_in_media; // Only spotlights emit in this mode.
+    let mut power = power * (per_photon / pick_prob);
+    let mut ray = Ray::from_unit(origin, dir);
+    let mut state = MediumState::Outside;
+    // Whether the photon has been through glass or off a mirror yet.
+    let mut specular = false;
+
+    for _ in 0..MAX_PHOTON_BOUNCES {
+        let hit = scene.hit(&ray, EPSILON, f64::INFINITY);
+        let t_hit = hit.as_ref().map_or(f64::INFINITY, |h| h.t);
+        if let Some((medium, t0, t1)) = scene.medium_segment(state, &ray, t_hit) {
+            let t = t0 + medium.sample_distance(rng);
+            if t < t1 {
+                // First non-specular event: stored if it is a caustic photon
+                // (or a spotlight's, with direct_in_media), and either way
+                // this photon's job is done (multiply scattered light is path
+                // traced).
+                if specular || store_direct {
+                    let power = power.mul_elem(medium.absorption(t - t0));
+                    volume.push(Photon { pos: ray.at(t), dir: ray.direction, power });
+                }
+                return;
+            }
+            power = power.mul_elem(medium.absorption(t1 - t0));
+        }
+
+        let Some(rec) = hit else {
+            return;
+        };
+        match rec.material {
+            Material::Light { .. } => return,
+            Material::Diffuse { .. } => {
+                if specular {
+                    surface.push(Photon { pos: rec.point, dir: ray.direction, power });
+                }
+                return;
+            }
+            Material::Volume { medium } => {
+                state = if rec.front_face { MediumState::Inside(Some(medium)) } else { MediumState::Outside };
+                ray = Ray::from_unit(rec.point, ray.direction);
+            }
+            Material::Mirror { .. } | Material::Dielectric { .. } => {
+                let Some((attenuation, scattered)) = rec.material.scatter(&ray, &rec, rng) else {
+                    return;
+                };
+                if let Material::Dielectric { medium, .. } = rec.material {
+                    if scattered.direction.dot(rec.normal) < 0.0 {
+                        state = if rec.front_face { MediumState::Inside(medium) } else { MediumState::Outside };
+                    }
+                }
+                power = power.mul_elem(attenuation);
+                ray = scattered;
+                specular = true;
+            }
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// A binary tree over photon positions (median split, leaves of up to 8),
+// with every node's bounds grown by the gathering radius so a query only
+// has to test whether its ray or point is inside a node's box.
+
+struct Node {
+    lo: Vec3,
+    hi: Vec3,
+    start: u32,
+    count: u32, // > 0 for leaves
+    right: u32, // left child is the next node
+}
+
+pub struct PhotonTree {
+    photons: Vec<Photon>,
+    nodes: Vec<Node>,
+    radius: f64,
+}
+
+impl PhotonTree {
+    fn new(mut photons: Vec<Photon>, radius: f64) -> Self {
+        let mut nodes = Vec::new();
+        if !photons.is_empty() {
+            build(&mut nodes, &mut photons, 0, radius);
+        }
+        PhotonTree { photons, nodes, radius }
+    }
+
+    pub fn len(&self) -> usize {
+        self.photons.len()
+    }
+
+    /// Calls `f(photon, t, d2)` for each photon within the radius of the ray
+    /// segment o + t d, t in [t0, t1], where t is the photon's projection on
+    /// the ray and d2 its squared distance from it.
+    fn near_segment(&self, o: Vec3, d: Vec3, t0: f64, t1: f64, mut f: impl FnMut(&Photon, f64, f64)) {
+        if self.nodes.is_empty() {
+            return;
+        }
+        let inv = |x: f64| {
+            let i = 1.0 / x;
+            if i.is_finite() {
+                i
+            } else {
+                1e300_f64.copysign(i)
+            }
+        };
+        let inv = Vec3::new(inv(d.x), inv(d.y), inv(d.z));
+        let r2 = self.radius * self.radius;
+        let mut stack = [0u32; 64];
+        let mut sp = 1;
+        while sp > 0 {
+            sp -= 1;
+            let idx = stack[sp];
+            let node = &self.nodes[idx as usize];
+            let (mut near, mut far) = (t0, t1);
+            for a in 0..3 {
+                let ta = (node.lo[a] - o[a]) * inv[a];
+                let tb = (node.hi[a] - o[a]) * inv[a];
+                near = near.max(ta.min(tb));
+                far = far.min(ta.max(tb));
+            }
+            if near > far {
+                continue;
+            }
+            if node.count > 0 {
+                for p in &self.photons[node.start as usize..(node.start + node.count) as usize] {
+                    let v = p.pos - o;
+                    let t = v.dot(d);
+                    if t < t0 || t > t1 {
+                        continue;
+                    }
+                    let d2 = v.dot(v) - t * t;
+                    if d2 < r2 {
+                        f(p, t, d2);
+                    }
+                }
+            } else {
+                stack[sp] = idx + 1;
+                stack[sp + 1] = node.right;
+                sp += 2;
+            }
+        }
+    }
+
+    /// Calls `f(photon, d2)` for each photon within the radius of `x`.
+    fn near_point(&self, x: Vec3, mut f: impl FnMut(&Photon, f64)) {
+        if self.nodes.is_empty() {
+            return;
+        }
+        let r2 = self.radius * self.radius;
+        let mut stack = [0u32; 64];
+        let mut sp = 1;
+        while sp > 0 {
+            sp -= 1;
+            let idx = stack[sp];
+            let node = &self.nodes[idx as usize];
+            if (0..3).any(|a| x[a] < node.lo[a] || x[a] > node.hi[a]) {
+                continue;
+            }
+            if node.count > 0 {
+                for p in &self.photons[node.start as usize..(node.start + node.count) as usize] {
+                    let v = p.pos - x;
+                    let d2 = v.dot(v);
+                    if d2 < r2 {
+                        f(p, d2);
+                    }
+                }
+            } else {
+                stack[sp] = idx + 1;
+                stack[sp + 1] = node.right;
+                sp += 2;
+            }
+        }
+    }
+}
+
+fn build(nodes: &mut Vec<Node>, photons: &mut [Photon], offset: usize, radius: f64) -> usize {
+    let (mut lo, mut hi) = (Vec3::splat(f64::INFINITY), Vec3::splat(f64::NEG_INFINITY));
+    for p in photons.iter() {
+        lo = Vec3::new(lo.x.min(p.pos.x), lo.y.min(p.pos.y), lo.z.min(p.pos.z));
+        hi = Vec3::new(hi.x.max(p.pos.x), hi.y.max(p.pos.y), hi.z.max(p.pos.z));
+    }
+    let idx = nodes.len();
+    nodes.push(Node {
+        lo: lo - Vec3::splat(radius),
+        hi: hi + Vec3::splat(radius),
+        start: offset as u32,
+        count: photons.len() as u32,
+        right: 0,
+    });
+    if photons.len() <= 8 {
+        return idx;
+    }
+    let extent = hi - lo;
+    let axis = if extent.x >= extent.y && extent.x >= extent.z {
+        0
+    } else if extent.y >= extent.z {
+        1
+    } else {
+        2
+    };
+    let half = photons.len() / 2;
+    photons.select_nth_unstable_by(half, |a, b| a.pos[axis].total_cmp(&b.pos[axis]));
+    let (left, right) = photons.split_at_mut(half);
+    build(nodes, left, offset, radius);
+    let right_idx = build(nodes, right, offset + half, radius);
+    nodes[idx].count = 0;
+    nodes[idx].right = right_idx as u32;
+    idx
+}
+
+/// Silverman's 2D biweight kernel: smooth, integrates to 1 over the disc.
+fn kernel(d2: f64, r2: f64) -> f64 {
+    let x = 1.0 - d2 / r2;
+    3.0 / (std::f64::consts::PI * r2) * x * x
+}
+
+pub struct CausticMaps {
+    volume: PhotonTree,
+    surface: PhotonTree,
+    /// The volume map includes the spotlights' direct light, so the path
+    /// tracer must not also sample spotlights from points in media.
+    pub direct_in_media: bool,
+}
+
+impl CausticMaps {
+    pub fn build(scene: &Scene, settings: CausticPhotons) -> Self {
+        let (volume, surface) = trace_photons(scene, settings.emitted, settings.direct_in_media);
+        CausticMaps {
+            volume: PhotonTree::new(volume, settings.volume_radius),
+            surface: PhotonTree::new(surface, settings.surface_radius),
+            direct_in_media: settings.direct_in_media,
+        }
+    }
+
+    pub fn counts(&self) -> (usize, usize) {
+        (self.volume.len(), self.surface.len())
+    }
+
+    /// Caustic light scattered towards the camera along the part [t0, t1] of
+    /// `ray` inside `medium`: the beam radiance estimate, summing every
+    /// photon near the segment, weighted by the kernel, the phase function
+    /// (photon direction into -ray direction) and the transmittance back
+    /// along the ray to where the segment enters the medium.
+    pub fn beam(&self, medium: &Medium, ray: &Ray, t0: f64, t1: f64) -> Vec3 {
+        let r2 = self.volume.radius * self.volume.radius;
+        let mut sum = Vec3::ZERO;
+        // A ray crossing a dense beam can pass hundreds of photons, so the
+        // common cases skip work: isotropic scattering has a constant phase
+        // function, and grey extinction needs one exp instead of three.
+        let isotropic = medium.g == 0.0;
+        let grey = medium.sigma_a.x == medium.sigma_a.y && medium.sigma_a.y == medium.sigma_a.z;
+        let sigma_t = medium.sigma_a.x + medium.sigma_s;
+        self.volume.near_segment(ray.origin, ray.direction, t0, t1, |p, t, d2| {
+            let phase = if isotropic {
+                1.0 / (4.0 * std::f64::consts::PI)
+            } else {
+                hg_phase(-p.dir.dot(ray.direction), medium.g)
+            };
+            let weight = kernel(d2, r2) * phase;
+            if grey {
+                sum += p.power * (weight * (-sigma_t * (t - t0)).exp());
+            } else {
+                sum += p.power.mul_elem(medium.transmittance(t - t0)) * weight;
+            }
+        });
+        sum
+    }
+
+    /// Caustic light reflected by a diffuse surface at `x` (normal `n`
+    /// facing the viewer): albedo / pi times the photon irradiance estimate.
+    /// Photons must lie on the surface and arrive from the viewer's side.
+    pub fn surface(&self, x: Vec3, n: Vec3, albedo: Vec3) -> Vec3 {
+        let r = self.surface.radius;
+        let r2 = r * r;
+        let mut sum = Vec3::ZERO;
+        self.surface.near_point(x, |p, d2| {
+            if p.dir.dot(n) < 0.0 && (p.pos - x).dot(n).abs() < 0.05 * r {
+                sum += p.power * kernel(d2, r2);
+            }
+        });
+        sum.mul_elem(albedo / std::f64::consts::PI)
+    }
+}

rust/src/scene.rs

@@ -0,0 +1,762 @@
+use rand::Rng;
+
+use crate::bvh::{Aabb, Bvh};
+use crate::medium::{Fog, Medium, MediumState};
+use crate::vec3::{reflect, refract, schlick, Vec3};
+
+pub const EPSILON: f64 = 1e-6;
+
+/// Two unit vectors that form an orthonormal basis with the unit vector `n`,
+/// using the branchless construction of Duff et al. 2017 ("Building an
+/// Orthonormal Basis, Revisited").
+pub fn orthonormal_basis(n: Vec3) -> (Vec3, Vec3) {
+    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, bt)
+}
+
+/// Cosine-weighted direction in the hemisphere around the unit vector `n`;
+/// 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 (t, bt) = orthonormal_basis(n);
+    t * x + bt * y + n * z
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct Ray {
+    pub origin: Vec3,
+    pub direction: Vec3,
+}
+
+impl Ray {
+    pub fn new(origin: Vec3, direction: Vec3) -> Self {
+        Ray {
+            origin,
+            direction: direction.normalize(),
+        }
+    }
+
+    /// 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
+    }
+}
+
+#[derive(Clone, Copy, Debug)]
+pub enum Material {
+    Diffuse { albedo: Vec3 },
+    Mirror { albedo: Vec3 },
+    /// Glass-like surface; `medium` fills its interior (None: clear).
+    Dielectric { ref_idx: f64, albedo: Vec3, medium: Option<Medium> },
+    /// `spot` restricts emission to a cone around the surface normal.
+    Light { emission: Vec3, spot: Option<Spot> },
+    /// Invisible boundary of a volume (smoke, cloud): rays pass straight
+    /// through, but inside it they travel through `medium`.
+    Volume { medium: Medium },
+}
+
+pub struct HitRecord {
+    pub t: f64,
+    pub point: Vec3,
+    pub normal: Vec3,
+    pub front_face: bool,
+    pub material: Material,
+}
+
+impl HitRecord {
+    /// Orients the normal against the incoming ray, as in the Python version.
+    fn new(r: &Ray, t: f64, outward: Vec3, material: Material) -> Self {
+        let front_face = r.direction.dot(outward) < 0.0;
+        HitRecord {
+            t,
+            point: r.at(t),
+            normal: if front_face { outward } else { -outward },
+            front_face,
+            material,
+        }
+    }
+}
+
+impl Material {
+    /// Returns (attenuation, scattered ray), or None if the ray is absorbed.
+    pub fn scatter<R: Rng>(&self, ray_in: &Ray, rec: &HitRecord, rng: &mut R) -> Option<(Vec3, Ray)> {
+        match *self {
+            Material::Diffuse { albedo } => {
+                // 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);
+                if reflected.dot(rec.normal) <= 0.0 {
+                    return None;
+                }
+                Some((albedo, Ray::new(rec.point, reflected)))
+            }
+            Material::Dielectric { ref_idx, albedo, .. } => {
+                let ratio = if rec.front_face { 1.0 / ref_idx } else { ref_idx };
+                let d = ray_in.direction;
+                let cos_theta = (-d.dot(rec.normal)).min(1.0);
+                let sin_theta = (1.0 - cos_theta * cos_theta).max(0.0).sqrt();
+                let cannot_refract = ratio * sin_theta > 1.0;
+                let direction = if cannot_refract || schlick(cos_theta, ratio) > rng.gen::<f64>() {
+                    reflect(d, rec.normal)
+                } else {
+                    refract(d, rec.normal, ratio).unwrap_or_else(|| reflect(d, rec.normal))
+                };
+                Some((albedo, Ray::new(rec.point, direction)))
+            }
+            Material::Light { .. } => None,
+            Material::Volume { .. } => Some((Vec3::ONE, Ray::from_unit(rec.point, ray_in.direction))),
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug)]
+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 {
+    fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<(f64, Vec3)> {
+        match *self {
+            Shape::Sphere { center, radius } => {
+                let oc = r.origin - center;
+                let a = r.direction.dot(r.direction);
+                let b = oc.dot(r.direction);
+                let c = oc.dot(oc) - radius * radius;
+                let disc = b * b - a * c;
+                if disc < 0.0 {
+                    return None;
+                }
+                let sq = disc.sqrt();
+                let mut root = (-b - sq) / a;
+                if root < t_min || root > t_max {
+                    root = (-b + sq) / a;
+                    if root < t_min || root > t_max {
+                        return None;
+                    }
+                }
+                let outward = (r.at(root) - center).normalize();
+                Some((root, outward))
+            }
+            Shape::Plane { point, normal } => {
+                let denom = r.direction.dot(normal);
+                if denom.abs() < 1e-12 {
+                    return None;
+                }
+                let t = (point - r.origin).dot(normal) / denom;
+                if t < t_min || t > t_max {
+                    return None;
+                }
+                Some((t, normal))
+            }
+            Shape::Box { min, max } => {
+                let mut t_lo = t_min;
+                let mut t_hi = t_max;
+                // Python indexes with -1 when no slab raises t_lo, which means the z axis.
+                let mut hit_axis = 2;
+                let mut inv_ds = [0.0; 3];
+                for axis in 0..3 {
+                    let d = r.direction[axis];
+                    let inv_d = if d.abs() > 1e-12 { 1.0 / d } else { 1e12 };
+                    inv_ds[axis] = inv_d;
+                    let mut t0 = (min[axis] - r.origin[axis]) * inv_d;
+                    let mut t1 = (max[axis] - r.origin[axis]) * inv_d;
+                    if inv_d < 0.0 {
+                        std::mem::swap(&mut t0, &mut t1);
+                    }
+                    if t0 > t_lo {
+                        t_lo = t0;
+                        hit_axis = axis;
+                    }
+                    t_hi = t_hi.min(t1);
+                    if t_hi < t_lo {
+                        return None;
+                    }
+                }
+                if t_lo >= t_max || t_lo <= t_min {
+                    return None;
+                }
+                // A ray travelling +axis enters through the min face, whose outward normal is -axis.
+                let sign = if inv_ds[hit_axis] > 0.0 { -1.0 } else { 1.0 };
+                let mut n = [0.0; 3];
+                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,
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct Object {
+    pub shape: Shape,
+    pub material: Material,
+}
+
+/// Emission limited to a cone around a light's surface normal, fading out
+/// smoothly between the inner and outer angle.
+#[derive(Clone, Copy, Debug)]
+pub struct Spot {
+    outer: f64, // radians
+    cos_outer: f64,
+    cos_inner: f64,
+}
+
+impl Spot {
+    pub fn new(outer_deg: f64, inner_deg: f64) -> Self {
+        Spot {
+            outer: outer_deg.to_radians(),
+            cos_outer: outer_deg.to_radians().cos(),
+            cos_inner: inner_deg.to_radians().cos(),
+        }
+    }
+
+    /// Fraction of the full emission leaving at an angle with cosine `cos`
+    /// to the normal.
+    pub fn profile(&self, cos: f64) -> f64 {
+        let t = ((cos - self.cos_outer) / (self.cos_inner - self.cos_outer)).clamp(0.0, 1.0);
+        t * t * (3.0 - 2.0 * t)
+    }
+}
+
+/// Emitted radiance leaving at an angle with cosine `cos` to the normal.
+pub fn spot_emission(emission: Vec3, spot: Option<Spot>, cos: f64) -> Vec3 {
+    match spot {
+        Some(s) => emission * s.profile(cos),
+        None => emission,
+    }
+}
+
+/// An emitter that next event estimation can sample: a box, sphere or rectangle.
+#[derive(Clone, Copy, Debug)]
+pub struct Light {
+    shape: Shape,
+    pub emission: Vec3,
+    pub spot: Option<Spot>,
+    /// 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, spot: Option<Spot>) -> Self {
+        let bounds = shape.bounds().expect("planes cannot be lights");
+        let mut light = Light {
+            shape,
+            emission,
+            spot,
+            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"),
+        }
+    }
+
+    /// A photon leaving the light: its origin, direction and power (the
+    /// light's whole emitted power; the caller divides by the photon count).
+    pub fn emit<R: Rng>(&self, rng: &mut R) -> (Vec3, Vec3, Vec3) {
+        let (p, n) = self.sample(rng);
+        match self.spot {
+            None => {
+                // Cosine-weighted directions: power = L * area * pi.
+                let dir = cosine_hemisphere(n, rng);
+                (p, dir, self.emission * (self.area() * std::f64::consts::PI))
+            }
+            Some(spot) => {
+                // Uniform directions within the outer cone, weighted by the
+                // profile and the cosine over the direction pdf.
+                let cos = 1.0 - rng.gen::<f64>() * (1.0 - spot.cos_outer);
+                let sin = (1.0 - cos * cos).max(0.0).sqrt();
+                let phi = 2.0 * std::f64::consts::PI * rng.gen::<f64>();
+                let (t, b) = orthonormal_basis(n);
+                let dir = t * (sin * phi.cos()) + b * (sin * phi.sin()) + n * cos;
+                let pdf = 1.0 / (2.0 * std::f64::consts::PI * (1.0 - spot.cos_outer));
+                (p, dir, self.emission * (spot.profile(cos) * cos * self.area() / pdf))
+            }
+        }
+    }
+
+    /// 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.
+    ///
+    /// A spotlight gets exactly zero outside its cone: widened by the
+    /// light's angular radius seen from `point`, so that no point on the
+    /// light could reach it (otherwise light selection keeps picking a
+    /// powerful spotlight for points it can't light, and every other light
+    /// gets rare, heavily weighted samples: speckle noise).
+    fn importance(&self, point: Vec3) -> f64 {
+        let d = point - self.center;
+        if let (Some(spot), Shape::Rect { normal, .. }) = (self.spot, self.shape) {
+            let dist = d.length();
+            if dist > self.extent {
+                let angle = (normal.dot(d) / dist).clamp(-1.0, 1.0).acos();
+                if angle > spot.outer + (self.extent / dist).asin() {
+                    return 0.0;
+                }
+            }
+        }
+        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,
+    fog: Option<Fog>,
+    caustic_photons: Option<CausticPhotons>,
+}
+
+/// Photon mapping settings: photons emitted, the gathering radius in fog and
+/// on surfaces (larger: smoother but blurrier), and whether the fog's map
+/// also carries the spotlights' direct light.
+#[derive(Clone, Copy, Debug)]
+pub struct CausticPhotons {
+    pub emitted: usize,
+    pub volume_radius: f64,
+    pub surface_radius: f64,
+    pub direct_in_media: bool,
+}
+
+impl SceneBuilder {
+    pub fn new(view: View, background: Background) -> Self {
+        SceneBuilder {
+            objects: Vec::new(),
+            planes: Vec::new(),
+            lights: Vec::new(),
+            background,
+            view,
+            fog: None,
+            caustic_photons: None,
+        }
+    }
+
+    /// Fills the region min..max (bounds may be infinite) with `medium`,
+    /// everywhere outside objects.
+    pub fn fog(&mut self, min: Vec3, max: Vec3, medium: Medium) {
+        self.fog = Some(Fog { region: Aabb { min, max }, medium });
+    }
+
+    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, None);
+    }
+
+    pub fn sphere_light(&mut self, center: Vec3, radius: f64, emission: Vec3) {
+        self.add_light(Shape::Sphere { center, radius }, emission, None);
+    }
+
+    /// 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, None);
+    }
+
+    /// A rectangular spotlight: like rect_light, but emitting only within
+    /// `outer_deg` of `normal`, at full strength within `inner_deg`.
+    pub fn spot_light(&mut self, min: Vec3, max: Vec3, normal: Vec3, emission: Vec3, outer_deg: f64, inner_deg: f64) {
+        self.add_light(Shape::Rect { min, max, normal }, emission, Some(Spot::new(outer_deg, inner_deg)));
+    }
+
+    /// Shoot photons to render caustics (light through glass or off mirrors
+    /// onto fog and diffuse surfaces) from photon maps instead of by path
+    /// tracing alone. See photons.rs.
+    pub fn caustic_photons(&mut self, emitted: usize, volume_radius: f64, surface_radius: f64) {
+        self.caustic_photons = Some(CausticPhotons { emitted, volume_radius, surface_radius, direct_in_media: false });
+    }
+
+    /// Photons for the spotlights only, carrying both their caustics and
+    /// their direct light into fog (replacing light sampling of spotlights
+    /// at points in fog): much smoother beams, at the cost of storing (and
+    /// gathering) many more photons. Other lights are fully path traced.
+    pub fn media_photons(&mut self, emitted: usize, volume_radius: f64, surface_radius: f64) {
+        self.caustic_photons = Some(CausticPhotons { emitted, volume_radius, surface_radius, direct_in_media: true });
+    }
+
+    fn add_light(&mut self, shape: Shape, emission: Vec3, spot: Option<Spot>) {
+        self.lights.push(Light::new(shape, emission, spot));
+        self.add(shape, Material::Light { emission, spot });
+    }
+
+    pub fn build(self) -> Scene {
+        let bounds: Vec<Aabb> = self.objects.iter().map(|o| o.shape.bounds().unwrap()).collect();
+        let spheres: Vec<Option<(Vec3, f64)>> = self
+            .objects
+            .iter()
+            .map(|o| match o.shape {
+                Shape::Sphere { center, radius } => Some((center, radius)),
+                _ => None,
+            })
+            .collect();
+        let has_volume_boundaries = self
+            .objects
+            .iter()
+            .chain(&self.planes)
+            .any(|o| matches!(o.material, Material::Volume { .. }));
+        Scene {
+            bvh: Bvh::build(&bounds, &spheres),
+            objects: self.objects,
+            planes: self.planes,
+            lights: self.lights,
+            background: self.background,
+            view: self.view,
+            fog: self.fog,
+            has_volume_boundaries,
+            caustic_photons: self.caustic_photons,
+        }
+    }
+}
+
+pub struct Scene {
+    objects: Vec<Object>,
+    planes: Vec<Object>,
+    bvh: Bvh,
+    lights: Vec<Light>,
+    pub background: Background,
+    pub view: View,
+    fog: Option<Fog>,
+    /// Whether shadow rays may need to pass through invisible volume
+    /// boundaries (otherwise any hit blocks them).
+    has_volume_boundaries: bool,
+    pub caustic_photons: Option<CausticPhotons>,
+}
+
+impl Scene {
+    pub fn object_count(&self) -> usize {
+        self.objects.len() + self.planes.len()
+    }
+
+    pub fn bvh_depth(&self) -> usize {
+        self.bvh.depth()
+    }
+
+    pub fn uses_avx512(&self) -> bool {
+        self.bvh.uses_avx512()
+    }
+
+    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, zero outside a
+    /// spotlight's cone), so a scene with many lights costs one shadow ray
+    /// per shading point and mostly samples lights that matter. Spotlights
+    /// are left out with `exclude_spots` (when photons carry their light).
+    /// 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, exclude_spots: bool, rng: &mut R) -> Option<(&Light, f64)> {
+        let weight = |l: &Light| if exclude_spots && l.spot.is_some() { 0.0 } else { l.importance(point) };
+        if self.lights.len() == 1 {
+            let only = &self.lights[0];
+            return (weight(only) > 0.0).then_some((only, 1.0));
+        }
+        let total: f64 = self.lights.iter().map(weight).sum();
+        if total <= 0.0 {
+            return None;
+        }
+        let u = rng.gen::<f64>() * total;
+        let mut acc = 0.0;
+        let mut last = None;
+        for light in &self.lights {
+            let w = weight(light);
+            if w <= 0.0 {
+                continue;
+            }
+            acc += w;
+            if u < acc {
+                return Some((light, w / total));
+            }
+            last = Some((light, w / total));
+        }
+        // Rounding left u just past the last bucket.
+        last
+    }
+
+    /// Picks a light to emit a photon from (only spotlights if
+    /// `spots_only`), proportionally to its power. Returns the light and the
+    /// probability of picking it.
+    pub fn pick_emitting_light<R: Rng>(&self, spots_only: bool, rng: &mut R) -> Option<(&Light, f64)> {
+        let eligible = |l: &&Light| !spots_only || l.spot.is_some();
+        let total: f64 = self.lights.iter().filter(eligible).map(|l| l.power).sum();
+        if total <= 0.0 {
+            return None;
+        }
+        let u = rng.gen::<f64>() * total;
+        let mut acc = 0.0;
+        let mut last = None;
+        for light in self.lights.iter().filter(eligible) {
+            acc += light.power;
+            if u < acc {
+                return Some((light, light.power / total));
+            }
+            last = Some(light);
+        }
+        last.map(|l| (l, l.power / 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.planes {
+            if let Some((t, outward)) = obj.shape.hit(r, t_min, min_t) {
+                min_t = t;
+                closest = Some((t, outward, obj.material));
+            }
+        }
+        // Spheres are intersected inside the BVH (8 at a time); other shapes
+        // come back through this callback, which remembers their normal.
+        let mut other_normal = Vec3::ZERO;
+        let bvh_hit = self.bvh.closest(r, t_min, min_t, |i, t_max| {
+            let (t, outward) = self.objects[i].shape.hit(r, t_min, t_max)?;
+            other_normal = outward;
+            Some(t)
+        });
+        if let Some((i, t, is_sphere)) = bvh_hit {
+            let obj = &self.objects[i];
+            let outward = match obj.shape {
+                Shape::Sphere { center, .. } if is_sphere => (r.at(t) - center).normalize(),
+                _ => other_normal,
+            };
+            closest = Some((t, outward, obj.material));
+        }
+        closest.map(|(t, outward, material)| HitRecord::new(r, t, outward, material))
+    }
+
+    /// The medium along the ray segment [0, t_end) for a path in `state`, and
+    /// the part of the segment it covers.
+    pub fn medium_segment(&self, state: MediumState, r: &Ray, t_end: f64) -> Option<(Medium, f64, f64)> {
+        match state {
+            MediumState::Outside => {
+                let fog = self.fog.as_ref()?;
+                let (t0, t1) = fog.clip(r.origin, r.direction, 0.0, t_end)?;
+                Some((fog.medium, t0, t1))
+            }
+            MediumState::Inside(medium) => medium.map(|m| (m, 0.0, t_end)),
+        }
+    }
+
+    fn segment_transmittance(&self, state: MediumState, r: &Ray, t_end: f64) -> Vec3 {
+        match self.medium_segment(state, r, t_end) {
+            Some((m, t0, t1)) => m.transmittance(t1 - t0),
+            None => Vec3::ONE,
+        }
+    }
+
+    /// Fraction of light that travels from `origin` a distance `dist` along
+    /// `direction` for a path in `state`: zero if a surface is in the way,
+    /// otherwise the transmittance of the media passed through, including
+    /// any invisible volume boundaries crossed.
+    pub fn transmittance(&self, state: MediumState, origin: Vec3, direction: Vec3, dist: f64) -> Vec3 {
+        if !self.has_volume_boundaries {
+            // Any hit blocks the ray, so it never leaves its starting medium.
+            if self.blocked(origin, direction, dist) {
+                return Vec3::ZERO;
+            }
+            return self.segment_transmittance(state, &Ray::from_unit(origin, direction), dist);
+        }
+
+        let (mut origin, mut remaining, mut state) = (origin, dist, state);
+        let mut tr = Vec3::ONE;
+        // Walk from boundary to boundary; the bound only guards against
+        // pathological scenes.
+        for _ in 0..64 {
+            let r = Ray::from_unit(origin, direction);
+            let hit = self.hit(&r, EPSILON, remaining);
+            tr = tr.mul_elem(self.segment_transmittance(state, &r, hit.as_ref().map_or(remaining, |h| h.t)));
+            let Some(rec) = hit else {
+                return tr;
+            };
+            let Material::Volume { medium } = rec.material else {
+                return Vec3::ZERO;
+            };
+            state = if rec.front_face { MediumState::Inside(Some(medium)) } else { MediumState::Outside };
+            origin = rec.point;
+            remaining -= rec.t;
+        }
+        Vec3::ZERO
+    }
+
+    pub fn blocked(&self, origin: Vec3, direction: Vec3, t_max: f64) -> bool {
+        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())
+    }
+}

rust/src/scenes.rs

@@ -0,0 +1,487 @@
+//! The built-in scenes, selectable with `--scene NAME`.
+
+use rand::rngs::SmallRng;
+use rand::{Rng, SeedableRng};
+
+use crate::medium::Medium;
+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,
+    },
+    SceneInfo {
+        name: "city",
+        description: "Massive scale: a procedural city of ~125,000 objects under a low sun",
+        build: city,
+    },
+    SceneInfo {
+        name: "fog",
+        description: "The caustics scene in light fog: beams of focused light and shadow shafts",
+        build: fog,
+    },
+    SceneInfo {
+        name: "media",
+        description: "Participating media inside objects: jade, smoke, amber liquid and cloud",
+        build: media,
+    },
+    SceneInfo {
+        name: "beam",
+        description: "A spotlight through a glass sphere in dense fog: a volumetric caustic (photon mapped)",
+        build: beam,
+    },
+    SceneInfo {
+        name: "haze",
+        description: "The city under a layer of atmospheric haze, with a glow around the sun",
+        build: haze,
+    },
+];
+
+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, medium: None }
+}
+
+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 {
+    caustics_room(30.0).build()
+}
+
+/// The caustics scene with a square light of the given half-size (same
+/// total power whatever the size).
+fn caustics_room(light_half_size: f64) -> SceneBuilder {
+    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.
+    let (h, z) = (light_half_size, 360.0);
+    b.rect_light(
+        Vec3::new(-h, 250.0, z - h),
+        Vec3::new(h, 250.0, z + h),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(90.0 * (30.0 * 30.0) / (h * h)),
+    );
+
+    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
+}
+
+fn fog() -> Scene {
+    // A light 4x the area (same power): light that reaches it through the
+    // glass spheres after scattering in the fog (volumetric caustics, the
+    // beams) is found 4x as often at a quarter of the brightness, which cuts
+    // the fireflies those rare paths cause.
+    let mut b = caustics_room(60.0);
+    // Light fog (mean free path ~650 units, about the size of the room).
+    // Forward scattering, so the beams stand out when looking towards the
+    // light, as the camera does. It stops well below the light (at y = 250):
+    // sampling the light from fog right next to it divides by a tiny squared
+    // distance and turns into fireflies.
+    b.fog(
+        Vec3::new(f64::NEG_INFINITY, 0.0, f64::NEG_INFINITY),
+        Vec3::new(f64::INFINITY, 190.0, f64::INFINITY),
+        Medium::new(Vec3::splat(0.0001), 0.0015, 0.5),
+    );
+    // The caustics in the fog and on the floor come from photon maps.
+    b.caustic_photons(4_000_000, 4.0, 3.0);
+    b.build()
+}
+
+/// After the classic photon-mapping image of a volume caustic: a spotlight
+/// shines down through a glass sphere in a foggy room, and the sphere
+/// focuses its beam into a cone of light that narrows to a focus and widens
+/// again towards the floor.
+fn beam() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.0, 200.0, -420.0),
+            lookat: Vec3::new(0.0, 190.0, 150.0),
+            vfov: 50.0,
+        },
+        Background::Color(Vec3::ZERO),
+    );
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.6, 0.6, 0.6));
+    b.plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.5, 0.5, 0.5));
+
+    // A small spotlight pointing straight down; its 12 degree cone is
+    // narrower than the sphere, so the whole beam goes through the glass.
+    b.spot_light(
+        Vec3::new(-8.0, 400.0, 142.0),
+        Vec3::new(8.0, 400.0, 158.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::new(1.0, 0.95, 0.85) * 10000.0,
+        12.0,
+        8.0,
+    );
+    // A large, dim, cool fill light in the ceiling, so the room and the
+    // glass show.
+    b.rect_light(
+        Vec3::new(-400.0, 420.0, -300.0),
+        Vec3::new(400.0, 420.0, 500.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::new(0.15, 0.165, 0.21),
+    );
+    // A glass sphere focuses light about 1.5 radii from its centre: here
+    // at y = 130, halfway between the sphere and the floor.
+    b.sphere(Vec3::new(0.0, 220.0, 150.0), 60.0, glass(Vec3::ONE));
+
+    // Dense fog (mean free path ~330, twice as dense as the fog scene)
+    // below the light, scattering evenly in all directions so the beam
+    // shows from the side.
+    b.fog(
+        Vec3::new(f64::NEG_INFINITY, 0.0, f64::NEG_INFINITY),
+        Vec3::new(f64::INFINITY, 390.0, f64::INFINITY),
+        Medium::new(Vec3::splat(0.0003), 0.003, 0.0),
+    );
+    // Photons carry both the direct beam and the focused one through the
+    // fog: path tracing alone leaves a small spotlight's beam grainy.
+    b.media_photons(600_000, 3.5, 3.0);
+    b.build()
+}
+
+fn media() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.0, 150.0, -330.0),
+            lookat: Vec3::new(0.0, 40.0, 100.0),
+            vfov: 36.0,
+        },
+        Background::Color(Vec3::new(0.02, 0.02, 0.025)),
+    );
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.75, 0.75, 0.75));
+    b.plane(Vec3::new(0.0, 0.0, 400.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.55, 0.55, 0.6));
+    // A large soft key light above and slightly behind, plus a dim fill from the left.
+    b.rect_light(
+        Vec3::new(-120.0, 320.0, 60.0),
+        Vec3::new(120.0, 320.0, 220.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(7.0),
+    );
+    b.rect_light(
+        Vec3::new(400.0, 40.0, -60.0),
+        Vec3::new(400.0, 200.0, 140.0),
+        Vec3::new(-1.0, 0.0, 0.0),
+        Vec3::new(1.2, 1.3, 1.6),
+    );
+
+    let r = 42.0;
+    // Jade: glass filled with a dense, strongly scattering medium that
+    // absorbs red and blue, so light diffuses through it and comes out green.
+    let jade = Medium::new(Vec3::new(0.030, 0.003, 0.018), 0.2, 0.3);
+    b.sphere(Vec3::new(-165.0, r, 120.0), r, Material::Dielectric { ref_idx: 1.5, albedo: Vec3::ONE, medium: Some(jade) });
+    // Smoke: no surface at all, a grey absorbing and scattering volume.
+    let smoke = Medium::new(Vec3::splat(0.008), 0.018, 0.0);
+    b.sphere(Vec3::new(-55.0, r, 120.0), r, Material::Volume { medium: smoke });
+    // Amber liquid: clear glass holding a purely absorbing medium that
+    // absorbs blue most (Beer's law), so thicker parts look deeper amber.
+    let amber = Medium::new(Vec3::new(0.002, 0.011, 0.035), 0.0, 0.0);
+    b.sphere(Vec3::new(55.0, r, 120.0), r, Material::Dielectric { ref_idx: 1.33, albedo: Vec3::ONE, medium: Some(amber) });
+    // Cloud: dense, non-absorbing and strongly forward scattering.
+    let cloud = Medium::new(Vec3::ZERO, 0.15, 0.8);
+    b.sphere(Vec3::new(165.0, r, 120.0), r, Material::Volume { medium: cloud });
+    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()
+}
+
+/// A procedural city on a GRID x GRID lot layout: buildings (with stepped
+/// tiers on the tall ones) that rise towards downtown, some mirror-glass
+/// towers, and parks with sphere trees. Lit by a low sun and the sky.
+fn city() -> Scene {
+    city_builder().build()
+}
+
+fn haze() -> Scene {
+    let mut b = city_builder();
+    // Haze in the lowest 400 units of air, thin enough to see across a few
+    // kilometres (mean free path 2500) but enough to fade the horizon.
+    // Forward scattering puts a glow around the sun.
+    b.fog(
+        Vec3::new(f64::NEG_INFINITY, 0.0, f64::NEG_INFINITY),
+        Vec3::new(f64::INFINITY, 400.0, f64::INFINITY),
+        Medium::new(Vec3::splat(0.00004), 0.0004, 0.7),
+    );
+    b.build()
+}
+
+fn city_builder() -> SceneBuilder {
+    const GRID: usize = 400;
+    const LOT: f64 = 10.0;
+    // Every BLOCK-th row/column of lots is left empty as an avenue.
+    const BLOCK: usize = 5;
+    let extent = GRID as f64 * LOT;
+    let downtown = Vec3::new(0.5 * extent, 0.0, 0.5 * extent);
+
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.35 * extent, 120.0, 0.29 * extent),
+            lookat: Vec3::new(0.5 * extent, 0.0, 0.5 * extent),
+            vfov: 42.0,
+        },
+        Background::Sky,
+    );
+    let mut rng = SmallRng::seed_from_u64(2024);
+
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.22, 0.22, 0.23));
+
+    // The sun: a distant sphere light 1.5 degrees in radius, low and off to
+    // the side of the view so the buildings cast long shadows across it.
+    let sun_dir = Vec3::new(0.85, 0.26, -0.35).normalize();
+    let sun_dist = 50_000.0;
+    b.sphere_light(
+        downtown + sun_dir * sun_dist,
+        sun_dist * 1.5_f64.to_radians().tan(),
+        Vec3::new(1.0, 0.88, 0.72) * 1600.0,
+    );
+
+    let concrete = [
+        Vec3::new(0.72, 0.70, 0.66),
+        Vec3::new(0.56, 0.56, 0.60),
+        Vec3::new(0.76, 0.70, 0.60),
+        Vec3::new(0.46, 0.43, 0.41),
+        Vec3::new(0.82, 0.80, 0.77),
+        Vec3::new(0.62, 0.52, 0.46),
+    ];
+    let grass = diffuse(0.18, 0.36, 0.12);
+
+    for gx in 0..GRID {
+        for gz in 0..GRID {
+            if gx % BLOCK == BLOCK - 1 || gz % BLOCK == BLOCK - 1 {
+                continue;
+            }
+            let x0 = gx as f64 * LOT;
+            let z0 = gz as f64 * LOT;
+            let lot_center = Vec3::new(x0 + 0.5 * LOT, 0.0, z0 + 0.5 * LOT);
+            let r = (lot_center - downtown).length() / extent;
+
+            if rng.gen::<f64>() < 0.05 {
+                // Park: a lawn with a few round trees.
+                b.cuboid(Vec3::new(x0 + 0.5, 0.0, z0 + 0.5), Vec3::new(x0 + LOT - 0.5, 0.15, z0 + LOT - 0.5), grass);
+                for _ in 0..rng.gen_range(2..6) {
+                    let radius = 1.0 + 1.2 * rng.gen::<f64>();
+                    let c = Vec3::new(
+                        x0 + radius + 0.5 + (LOT - 2.0 * radius - 1.0) * rng.gen::<f64>(),
+                        radius + 0.8,
+                        z0 + radius + 0.5 + (LOT - 2.0 * radius - 1.0) * rng.gen::<f64>(),
+                    );
+                    let shade = 0.7 + 0.3 * rng.gen::<f64>();
+                    b.sphere(c, radius, diffuse(0.12 * shade, 0.32 * shade, 0.08 * shade));
+                }
+                continue;
+            }
+
+            // Height falls off from downtown, with a long tail of towers.
+            let skyline = 7.0 + 95.0 * (-(r * r) / 0.0068).exp();
+            let height = (skyline * (0.35 + 0.9 * rng.gen::<f64>().powi(2))).max(3.0);
+            // Tall buildings fill more of their lot.
+            let inset = if height > 40.0 { 0.4 + 0.6 * rng.gen::<f64>() } else { 0.8 + 1.4 * rng.gen::<f64>() };
+            let (mut lo, mut hi) = (
+                Vec3::new(x0 + inset, 0.0, z0 + inset),
+                Vec3::new(x0 + LOT - inset, height, z0 + LOT - inset),
+            );
+            let material = if height > 45.0 && rng.gen::<f64>() < 0.15 {
+                mirror(0.35, 0.45, 0.55)
+            } else {
+                let c = concrete[rng.gen_range(0..concrete.len())];
+                Material::Diffuse { albedo: c * (0.85 + 0.15 * rng.gen::<f64>()) }
+            };
+            b.cuboid(lo, hi, material);
+
+            // Tall buildings step in with up to three narrower tiers.
+            let mut tiers = if height > 30.0 { rng.gen_range(1..4) } else { 0 };
+            while tiers > 0 && hi.x - lo.x > 2.5 {
+                let step = 0.6 + 0.6 * rng.gen::<f64>();
+                let tier_height = (0.15 + 0.25 * rng.gen::<f64>()) * height;
+                lo = Vec3::new(lo.x + step, hi.y, lo.z + step);
+                hi = Vec3::new(hi.x - step, hi.y + tier_height, hi.z - step);
+                b.cuboid(lo, hi, material);
+                tiers -= 1;
+            }
+        }
+    }
+    b
+}

rust/src/vec3.rs

@@ -0,0 +1,149 @@
+use std::ops::{Add, AddAssign, Div, DivAssign, Index, Mul, Neg, Sub};
+
+#[derive(Clone, Copy, Debug, Default, PartialEq)]
+pub struct Vec3 {
+    pub x: f64,
+    pub y: f64,
+    pub z: f64,
+}
+
+impl Vec3 {
+    pub const ZERO: Vec3 = Vec3::new(0.0, 0.0, 0.0);
+    pub const ONE: Vec3 = Vec3::new(1.0, 1.0, 1.0);
+
+    pub const fn new(x: f64, y: f64, z: f64) -> Self {
+        Vec3 { x, y, z }
+    }
+
+    pub fn splat(v: f64) -> Self {
+        Vec3::new(v, v, v)
+    }
+
+    pub fn dot(self, o: Vec3) -> f64 {
+        self.x * o.x + self.y * o.y + self.z * o.z
+    }
+
+    pub fn cross(self, o: Vec3) -> Vec3 {
+        Vec3::new(
+            self.y * o.z - self.z * o.y,
+            self.z * o.x - self.x * o.z,
+            self.x * o.y - self.y * o.x,
+        )
+    }
+
+    pub fn length(self) -> f64 {
+        self.dot(self).sqrt()
+    }
+
+    pub fn normalize(self) -> Vec3 {
+        let n = self.length();
+        if n > 1e-12 {
+            self / n
+        } else {
+            self
+        }
+    }
+
+    /// Component-wise product.
+    pub fn mul_elem(self, o: Vec3) -> Vec3 {
+        Vec3::new(self.x * o.x, self.y * o.y, self.z * o.z)
+    }
+
+    pub fn map(self, f: impl Fn(f64) -> f64) -> Vec3 {
+        Vec3::new(f(self.x), f(self.y), f(self.z))
+    }
+
+    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 {
+    type Output = f64;
+    fn index(&self, i: usize) -> &f64 {
+        match i {
+            0 => &self.x,
+            1 => &self.y,
+            _ => &self.z,
+        }
+    }
+}
+
+impl Add for Vec3 {
+    type Output = Vec3;
+    fn add(self, o: Vec3) -> Vec3 {
+        Vec3::new(self.x + o.x, self.y + o.y, self.z + o.z)
+    }
+}
+
+impl AddAssign for Vec3 {
+    fn add_assign(&mut self, o: Vec3) {
+        *self = *self + o;
+    }
+}
+
+impl Sub for Vec3 {
+    type Output = Vec3;
+    fn sub(self, o: Vec3) -> Vec3 {
+        Vec3::new(self.x - o.x, self.y - o.y, self.z - o.z)
+    }
+}
+
+impl Neg for Vec3 {
+    type Output = Vec3;
+    fn neg(self) -> Vec3 {
+        Vec3::new(-self.x, -self.y, -self.z)
+    }
+}
+
+impl Mul<f64> for Vec3 {
+    type Output = Vec3;
+    fn mul(self, s: f64) -> Vec3 {
+        Vec3::new(self.x * s, self.y * s, self.z * s)
+    }
+}
+
+impl Mul<Vec3> for f64 {
+    type Output = Vec3;
+    fn mul(self, v: Vec3) -> Vec3 {
+        v * self
+    }
+}
+
+impl Div<f64> for Vec3 {
+    type Output = Vec3;
+    fn div(self, s: f64) -> Vec3 {
+        Vec3::new(self.x / s, self.y / s, self.z / s)
+    }
+}
+
+impl DivAssign<f64> for Vec3 {
+    fn div_assign(&mut self, s: f64) {
+        *self = *self / s;
+    }
+}
+
+pub fn reflect(d: Vec3, n: Vec3) -> Vec3 {
+    d - 2.0 * d.dot(n) * n
+}
+
+pub fn refract(incident: Vec3, normal: Vec3, eta: f64) -> Option<Vec3> {
+    let cos_i = (-incident.dot(normal)).min(1.0);
+    let sin_t_sq = eta * eta * (1.0 - cos_i * cos_i);
+    if sin_t_sq > 1.0 {
+        return None;
+    }
+    let cos_t = (1.0 - sin_t_sq).max(0.0).sqrt();
+    Some(eta * incident + (eta * cos_i - cos_t) * normal)
+}
+
+pub fn schlick(cos: f64, ref_idx: f64) -> f64 {
+    let r0 = (1.0 - ref_idx) / (1.0 + ref_idx);
+    let r0 = r0 * r0;
+    r0 + (1.0 - r0) * (1.0 - cos).powi(5)
+}