simd AVX512

README.md

@@ -68,7 +68,49 @@ 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 --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 |
+| `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`). Every emitter must be added with a `*_light` method so light sampling knows about it.
+
+### 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 +126,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.horizontal = 2.0 * focus_dist * half_width * u
-        self.vertical = focus_dist * half_height * v
+        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
+        return np.zeros(3)
-
-    indirect = trace(scene, scattered_ray, depth + 1, max_depth, rng)
-    color = emitted + attenuation * indirect
 
     if isinstance(mat, Diffuse):
         direct = np.zeros(3)
         for light_obj in scene.lights:
-            emission = light_obj.material.get_emission()
+            direct += sample_box_light(scene, light_obj, rec, mat.albedo, rng)
-            for _ in range(2):
+        indirect = trace(scene, scattered_ray, depth + 1, max_depth, rng, count_emission=False)
-                face_idx = int(rng.random() * 6)
+        return direct + attenuation * indirect
-                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]])
 
-                to_light = lp - rec.point
+    return attenuation * trace(scene, scattered_ray, depth + 1, max_depth, rng, count_emission=True)
-                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
 
 
 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)
+                u = (i + rng.random()) / width
-                v = (height - 1 - j + rng.random()) / (height - 1)
+                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):
+    for s in range(0, height, rows_per_chunk):
-        s = w * chunk
+        e = min(s + rows_per_chunk, height)
-        e = height if w == workers - 1 else (w + 1) * chunk
         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]),
+        ("Floor", image[h_end:, w_start:w_end]),
-        ("Ceiling", image[:height // 3, h_start:h_end]),
+        ("Ceiling", image[:h_start, w_start:w_end]),
-        ("Center", image[h_start:h_end, h_start:h_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,443 @@
+mod bvh;
+mod camera;
+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 scene::{HitRecord, 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;
+
+/// Direct light from one uniformly sampled point on `light`'s surface.
+fn sample_direct<R: Rng>(scene: &Scene, light: &Light, rec: &HitRecord, albedo: Vec3, rng: &mut R) -> Vec3 {
+    let (lp, light_normal) = light.sample(rng);
+
+    let to_light = lp - rec.point;
+    let dist_sq = to_light.dot(to_light);
+    if dist_sq < 1e-12 {
+        return Vec3::ZERO;
+    }
+    let dist = dist_sq.sqrt();
+    let d_to_light = to_light / dist;
+    let cos_at_rec = rec.normal.dot(d_to_light);
+    let cos_at_light = -light_normal.dot(d_to_light);
+    if cos_at_rec <= 0.0 || cos_at_light <= 0.0 {
+        return Vec3::ZERO;
+    }
+
+    if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist * (1.0 - SHADOW_RELATIVE_EPSILON)) {
+        return Vec3::ZERO;
+    }
+
+    // Uniform area sampling has pdf 1 / area.
+    let brdf = albedo / PI;
+    let geometry = cos_at_rec * cos_at_light / dist_sq;
+    light.emission.mul_elem(brdf) * (geometry * light.area())
+}
+
+/// Bounces before Russian roulette may terminate a path.
+const RR_MIN_DEPTH: u32 = 3;
+
+fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> Vec3 {
+    let mut radiance = Vec3::ZERO;
+    // Product of the attenuations along the path so far (divided by the RR
+    // survival probabilities), i.e. how much of any light found here reaches the camera.
+    let mut throughput = Vec3::ONE;
+    // False after a diffuse bounce: that light was already gathered by next
+    // event estimation, so hitting it again would double count.
+    let mut count_emission = true;
+    let mut ray = *primary;
+
+    for depth in 0..max_depth {
+        let Some(rec) = scene.hit(&ray, EPSILON, f64::INFINITY) else {
+            radiance += throughput.mul_elem(scene.background.radiance(ray.direction));
+            break;
+        };
+
+        if let Material::Light { emission } = rec.material {
+            // Lights emit from their front side only (matters for one-sided rects).
+            if count_emission && rec.front_face {
+                radiance += throughput.mul_elem(emission);
+            }
+            break;
+        }
+
+        let Some((attenuation, scattered)) = rec.material.scatter(&ray, &rec, rng) else {
+            break;
+        };
+
+        if let Material::Diffuse { albedo } = rec.material {
+            if let Some((light, pick_prob)) = scene.pick_light(rec.point, rng) {
+                let direct = sample_direct(scene, light, &rec, albedo, rng) / pick_prob;
+                radiance += throughput.mul_elem(direct);
+            }
+            count_emission = false;
+        } else {
+            count_emission = true;
+        }
+
+        throughput = throughput.mul_elem(attenuation);
+        ray = scattered;
+
+        // Russian roulette: end dim paths early, boosting survivors so the
+        // estimate stays unbiased. Glass/mirror paths (throughput ~1) survive.
+        if depth + 1 >= RR_MIN_DEPTH {
+            let survive = throughput.x.max(throughput.y).max(throughput.z).min(0.95);
+            if rng.gen::<f64>() >= survive {
+                break;
+            }
+            throughput = throughput / survive;
+        }
+    }
+
+    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,
+}
+
+impl Renderer<'_> {
+    fn sample(&self, px: &mut Pixel, i: usize, j: usize, count: u32) {
+        let (width, height) = (self.s.width, self.s.height);
+        for _ in 0..count {
+            let u = (i as f64 + px.rng.gen::<f64>()) / width as f64;
+            let v = ((height - 1 - j) as f64 + px.rng.gen::<f64>()) / height as f64;
+            let sample = trace(self.scene, &self.cam.get_ray(u, v), self.s.max_depth, &mut px.rng);
+            px.sum += sample;
+            px.stats.add(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) -> (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,
+    };
+
+    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] \
+[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");
+    let mut args: Vec<String> = Vec::new();
+    let mut raw = std::env::args().skip(1);
+    while let Some(a) = raw.next() {
+        if a == "--list-scenes" {
+            list_scenes();
+            return;
+        } else if a == "--scene" {
+            scene_name = raw.next().unwrap_or_else(|| {
+                eprintln!("{USAGE}");
+                std::process::exit(2);
+            });
+        } else if let Some(name) = a.strip_prefix("--scene=") {
+            scene_name = name.to_string();
+        } else if a == "--help" || a == "-h" {
+            println!("{USAGE}");
+            list_scenes();
+            return;
+        } else {
+            args.push(a);
+        }
+    }
+    let Some(scene_info) = scenes::find(&scene_name) else {
+        eprintln!("unknown scene '{scene_name}'");
+        list_scenes();
+        std::process::exit(2);
+    };
+
+    let arg = |i: usize, default: &str| args.get(i).cloned().unwrap_or_else(|| default.to_string());
+    fn parse<T: std::str::FromStr>(s: String) -> T {
+        s.parse().unwrap_or_else(|_| {
+            eprintln!("{USAGE}");
+            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 (image, total_samples) = render(&scene, &settings);
+    let elapsed = start.elapsed().as_secs_f64();
+    let pixels = (settings.width * settings.height) as f64;
+
+    if let Err(e) = save_png(&image, settings.width, settings.height, &settings.output) {
+        eprintln!("failed to write {}: {e}", settings.output);
+        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);
+    }
+}

rust/src/scene.rs

@@ -0,0 +1,533 @@
+use rand::Rng;
+
+use crate::bvh::{Aabb, Bvh};
+use crate::vec3::{reflect, refract, schlick, Vec3};
+
+pub const EPSILON: f64 = 1e-6;
+
+/// Cosine-weighted direction in the hemisphere around the unit vector `n`.
+/// Uses the branchless orthonormal basis of Duff et al. 2017 ("Building an
+/// Orthonormal Basis, Revisited"); the result is already unit length.
+pub fn cosine_hemisphere<R: Rng>(n: Vec3, rng: &mut R) -> Vec3 {
+    let r1: f64 = rng.gen();
+    let r2: f64 = rng.gen();
+    let phi = 2.0 * std::f64::consts::PI * r1;
+    let r = r2.sqrt();
+    let (x, y, z) = (r * phi.cos(), r * phi.sin(), (1.0 - r2).max(0.0).sqrt());
+
+    let sign = 1.0_f64.copysign(n.z);
+    let a = -1.0 / (sign + n.z);
+    let b = n.x * n.y * a;
+    let t = Vec3::new(1.0 + sign * n.x * n.x * a, sign * b, -sign * n.x);
+    let bt = Vec3::new(b, sign + n.y * n.y * a, -n.y);
+    t * x + bt * y + n * z
+}
+
+#[derive(Clone, Copy, Debug)]
+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 },
+    Dielectric { ref_idx: f64, albedo: Vec3 },
+    Light { emission: Vec3 },
+}
+
+pub struct HitRecord {
+    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 {
+            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,
+        }
+    }
+}
+
+#[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,
+}
+
+/// 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,
+    /// Total emitted power (up to a constant): luminance * area.
+    power: f64,
+    center: Vec3,
+    /// Half the bounding box diagonal: how close a point can be to the centre
+    /// while still outside the light.
+    extent: f64,
+}
+
+impl Light {
+    fn new(shape: Shape, emission: Vec3) -> Self {
+        let bounds = shape.bounds().expect("planes cannot be lights");
+        let mut light = Light {
+            shape,
+            emission,
+            power: 0.0,
+            center: bounds.centroid(),
+            extent: 0.5 * (bounds.max - bounds.min).length(),
+        };
+        light.power = emission.luminance() * light.area();
+        light
+    }
+
+    pub fn area(&self) -> f64 {
+        match self.shape {
+            Shape::Box { min, max } => {
+                let d = max - min;
+                2.0 * (d.x * d.y + d.y * d.z + d.x * d.z)
+            }
+            Shape::Sphere { radius, .. } => 4.0 * std::f64::consts::PI * radius * radius,
+            Shape::Rect { min, max, normal } => {
+                let d = max - min;
+                if normal.x != 0.0 {
+                    d.y * d.z
+                } else if normal.y != 0.0 {
+                    d.x * d.z
+                } else {
+                    d.x * d.y
+                }
+            }
+            Shape::Plane { .. } => unreachable!("planes cannot be lights"),
+        }
+    }
+
+    /// Rough estimate of how much this light contributes at `point`: its
+    /// power over the squared distance, clamped so points near or inside the
+    /// light's bounds don't blow up.
+    fn importance(&self, point: Vec3) -> f64 {
+        let d = point - self.center;
+        self.power / d.dot(d).max(self.extent * self.extent)
+    }
+
+    /// A uniformly distributed point on the surface (pdf 1 / area) and the
+    /// outward normal there.
+    pub fn sample<R: Rng>(&self, rng: &mut R) -> (Vec3, Vec3) {
+        match self.shape {
+            Shape::Box { min, max } => {
+                // Pick a face with probability proportional to its area, then
+                // a uniform point on it.
+                let size = max - min;
+                let face_areas = [size.y * size.z, size.x * size.z, size.x * size.y];
+                let u = rng.gen::<f64>() * (face_areas[0] + face_areas[1] + face_areas[2]);
+                let axis = if u < face_areas[0] {
+                    0
+                } else if u < face_areas[0] + face_areas[1] {
+                    1
+                } else {
+                    2
+                };
+                let positive_side = rng.gen::<f64>() < 0.5;
+                let mut p = [
+                    min.x + size.x * rng.gen::<f64>(),
+                    min.y + size.y * rng.gen::<f64>(),
+                    min.z + size.z * rng.gen::<f64>(),
+                ];
+                p[axis] = if positive_side { max[axis] } else { min[axis] };
+                let mut n = [0.0; 3];
+                n[axis] = if positive_side { 1.0 } else { -1.0 };
+                (Vec3::new(p[0], p[1], p[2]), Vec3::new(n[0], n[1], n[2]))
+            }
+            Shape::Sphere { center, radius } => {
+                let z = 1.0 - 2.0 * rng.gen::<f64>();
+                let r = (1.0 - z * z).max(0.0).sqrt();
+                let phi = 2.0 * std::f64::consts::PI * rng.gen::<f64>();
+                let n = Vec3::new(r * phi.cos(), r * phi.sin(), z);
+                (center + radius * n, n)
+            }
+            Shape::Rect { min, max, normal } => {
+                let size = max - min;
+                let p = Vec3::new(
+                    min.x + size.x * rng.gen::<f64>(),
+                    min.y + size.y * rng.gen::<f64>(),
+                    min.z + size.z * rng.gen::<f64>(),
+                );
+                (p, normal)
+            }
+            Shape::Plane { .. } => unreachable!("planes cannot be lights"),
+        }
+    }
+}
+
+/// What a ray that leaves the scene sees.
+#[derive(Clone, Copy, Debug)]
+pub enum Background {
+    /// White-to-blue gradient, as in raytracer.py.
+    Sky,
+    Color(Vec3),
+}
+
+impl Background {
+    pub fn radiance(&self, direction: Vec3) -> Vec3 {
+        match *self {
+            Background::Sky => {
+                let t = 0.5 * (direction.y + 1.0);
+                Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t
+            }
+            Background::Color(c) => c,
+        }
+    }
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct View {
+    pub lookfrom: Vec3,
+    pub lookat: Vec3,
+    pub vfov: f64,
+}
+
+pub struct SceneBuilder {
+    objects: Vec<Object>,
+    planes: Vec<Object>,
+    lights: Vec<Light>,
+    background: Background,
+    view: View,
+}
+
+impl SceneBuilder {
+    pub fn new(view: View, background: Background) -> Self {
+        SceneBuilder {
+            objects: Vec::new(),
+            planes: Vec::new(),
+            lights: Vec::new(),
+            background,
+            view,
+        }
+    }
+
+    pub fn add(&mut self, shape: Shape, material: Material) {
+        let obj = Object { shape, material };
+        match shape {
+            Shape::Plane { .. } => self.planes.push(obj),
+            _ => self.objects.push(obj),
+        }
+    }
+
+    pub fn plane(&mut self, point: Vec3, normal: Vec3, material: Material) {
+        self.add(Shape::Plane { point, normal: normal.normalize() }, material);
+    }
+
+    pub fn sphere(&mut self, center: Vec3, radius: f64, material: Material) {
+        self.add(Shape::Sphere { center, radius }, material);
+    }
+
+    pub fn cuboid(&mut self, min: Vec3, max: Vec3, material: Material) {
+        self.add(Shape::Box { min, max }, material);
+    }
+
+    /// Every emitter must be added through a *_light method: after a diffuse
+    /// bounce, trace() only counts light found by sampling these.
+    pub fn box_light(&mut self, min: Vec3, max: Vec3, emission: Vec3) {
+        self.add_light(Shape::Box { min, max }, emission);
+    }
+
+    pub fn sphere_light(&mut self, center: Vec3, radius: f64, emission: Vec3) {
+        self.add_light(Shape::Sphere { center, radius }, emission);
+    }
+
+    /// One-sided rectangular light emitting towards `normal` (±x, ±y or ±z).
+    pub fn rect_light(&mut self, min: Vec3, max: Vec3, normal: Vec3, emission: Vec3) {
+        self.add_light(Shape::Rect { min, max, normal }, emission);
+    }
+
+    fn add_light(&mut self, shape: Shape, emission: Vec3) {
+        self.lights.push(Light::new(shape, emission));
+        self.add(shape, Material::Light { emission });
+    }
+
+    pub fn build(self) -> Scene {
+        let bounds: Vec<Aabb> = self.objects.iter().map(|o| o.shape.bounds().unwrap()).collect();
+        let spheres: Vec<Option<(Vec3, f64)>> = self
+            .objects
+            .iter()
+            .map(|o| match o.shape {
+                Shape::Sphere { center, radius } => Some((center, radius)),
+                _ => None,
+            })
+            .collect();
+        Scene {
+            bvh: Bvh::build(&bounds, &spheres),
+            objects: self.objects,
+            planes: self.planes,
+            lights: self.lights,
+            background: self.background,
+            view: self.view,
+        }
+    }
+}
+
+pub struct Scene {
+    objects: Vec<Object>,
+    planes: Vec<Object>,
+    bvh: Bvh,
+    lights: Vec<Light>,
+    pub background: Background,
+    pub view: View,
+}
+
+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), so a scene with
+    /// many lights costs one shadow ray per shading point and mostly samples
+    /// nearby lights. Returns the light and the probability of picking it.
+    /// This is O(lights) per call; thousands of lights would want a light tree.
+    pub fn pick_light<R: Rng>(&self, point: Vec3, rng: &mut R) -> Option<(&Light, f64)> {
+        match self.lights.len() {
+            0 => None,
+            1 => Some((&self.lights[0], 1.0)),
+            _ => {
+                let total: f64 = self.lights.iter().map(|l| l.importance(point)).sum();
+                let u = rng.gen::<f64>() * total;
+                let mut acc = 0.0;
+                for light in &self.lights {
+                    let w = light.importance(point);
+                    acc += w;
+                    if u < acc {
+                        return Some((light, w / total));
+                    }
+                }
+                // Rounding left u just past the last bucket.
+                let last = self.lights.last().unwrap();
+                Some((last, last.importance(point) / total))
+            }
+        }
+    }
+
+    pub fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
+        let mut closest: Option<(f64, Vec3, Material)> = None;
+        let mut min_t = t_max;
+        for obj in &self.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))
+    }
+
+    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,326 @@
+//! The built-in scenes, selectable with `--scene NAME`.
+
+use rand::rngs::SmallRng;
+use rand::{Rng, SeedableRng};
+
+use crate::scene::{Background, Material, Scene, SceneBuilder, View};
+use crate::vec3::Vec3;
+
+pub struct SceneInfo {
+    pub name: &'static str,
+    pub description: &'static str,
+    pub build: fn() -> Scene,
+}
+
+pub const SCENES: &[SceneInfo] = &[
+    SceneInfo {
+        name: "cornell",
+        description: "Cornell box from raytracer.py: glass, mirror and diffuse spheres, wall light",
+        build: cornell,
+    },
+    SceneInfo {
+        name: "caustics",
+        description: "Clear, amber and blue glass plus a mirror under one light, focusing caustics onto the floor",
+        build: caustics,
+    },
+    SceneInfo {
+        name: "spheres",
+        description: "Large scene: ~480 random spheres under a sky (exercises the BVH)",
+        build: spheres,
+    },
+    SceneInfo {
+        name: "night",
+        description: "The same sphere field at night, lit by dozens of glowing spheres (many-light sampling)",
+        build: night,
+    },
+    SceneInfo {
+        name: "classic",
+        description: "The classic Cornell box: ceiling light, two boxes, colour bleeding",
+        build: classic,
+    },
+    SceneInfo {
+        name: "city",
+        description: "Massive scale: a procedural city of ~125,000 objects under a low sun",
+        build: city,
+    },
+];
+
+pub fn find(name: &str) -> Option<&'static SceneInfo> {
+    SCENES.iter().find(|s| s.name == name)
+}
+
+fn diffuse(r: f64, g: f64, b: f64) -> Material {
+    Material::Diffuse { albedo: Vec3::new(r, g, b) }
+}
+
+fn mirror(r: f64, g: f64, b: f64) -> Material {
+    Material::Mirror { albedo: Vec3::new(r, g, b) }
+}
+
+/// Glass; `tint` is applied at every reflection/refraction, so colour
+/// deepens with each pass through the surface.
+fn glass(tint: Vec3) -> Material {
+    Material::Dielectric { ref_idx: 1.5, albedo: tint }
+}
+
+fn cornell() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(278.0, 278.0, -800.0),
+            lookat: Vec3::new(278.0, 278.0, 278.0),
+            vfov: 40.0,
+        },
+        Background::Sky,
+    );
+
+    let white = diffuse(0.73, 0.73, 0.73);
+    b.plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), diffuse(0.5, 0.5, 0.5));
+    b.plane(Vec3::new(0.0, 500.0, 0.0), Vec3::new(0.0, -1.0, 0.0), white);
+    b.plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0), white);
+    b.plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0), diffuse(0.65, 0.05, 0.05));
+    b.plane(Vec3::new(500.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), diffuse(0.12, 0.45, 0.15));
+
+    b.box_light(Vec3::new(0.0, 180.0, 180.0), Vec3::new(25.0, 220.0, 320.0), Vec3::splat(50.0));
+
+    b.sphere(Vec3::new(250.0, 180.0, 250.0), 80.0, glass(Vec3::ONE));
+    b.sphere(Vec3::new(150.0, 180.0, 380.0), 80.0, mirror(0.95, 0.95, 0.95));
+    b.sphere(Vec3::new(350.0, 160.0, 300.0), 60.0, diffuse(0.65, 0.05, 0.05));
+    b.build()
+}
+
+fn caustics() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.0, 160.0, -280.0),
+            lookat: Vec3::new(20.0, 30.0, 120.0),
+            vfov: 38.0,
+        },
+        Background::Color(Vec3::ZERO),
+    );
+
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.8, 0.8, 0.8));
+    b.plane(Vec3::new(0.0, 0.0, 520.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.6, 0.6, 0.6));
+
+    // One downward-facing panel above and well behind the spheres. A glass
+    // sphere focuses light about 1.5 radii from its centre; with the light
+    // ~48 degrees from vertical that focus lands on the floor just in front
+    // of each sphere, facing the camera. A smaller panel gives sharper
+    // caustics but makes them slower to converge.
+    b.rect_light(
+        Vec3::new(-30.0, 250.0, 330.0),
+        Vec3::new(30.0, 250.0, 390.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(90.0),
+    );
+
+    b.sphere(Vec3::new(-75.0, 45.0, 110.0), 45.0, glass(Vec3::ONE));
+    b.sphere(Vec3::new(40.0, 38.0, 80.0), 38.0, glass(Vec3::new(1.0, 0.7, 0.3)));
+    b.sphere(Vec3::new(125.0, 30.0, 150.0), 30.0, glass(Vec3::new(0.4, 0.65, 1.0)));
+    b.sphere(Vec3::new(-10.0, 35.0, 210.0), 35.0, mirror(0.9, 0.9, 0.9));
+    b.build()
+}
+
+/// A grid of small random spheres around three large feature spheres (after
+/// "Ray Tracing in One Weekend"). A `light_fraction` of the small spheres glow.
+fn sphere_field(b: &mut SceneBuilder, light_fraction: f64) {
+    let mut rng = SmallRng::seed_from_u64(7);
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.5, 0.5, 0.5));
+
+    let big = [
+        (Vec3::new(0.0, 1.0, 0.0), glass(Vec3::ONE)),
+        (Vec3::new(-4.0, 1.0, 0.0), diffuse(0.4, 0.2, 0.1)),
+        (Vec3::new(4.0, 1.0, 0.0), mirror(0.7, 0.6, 0.5)),
+    ];
+    for (center, material) in big {
+        b.sphere(center, 1.0, material);
+    }
+
+    for a in -11..11 {
+        for c in -11..11 {
+            let center = Vec3::new(
+                a as f64 + 0.9 * rng.gen::<f64>(),
+                0.2,
+                c as f64 + 0.9 * rng.gen::<f64>(),
+            );
+            if big.iter().any(|(bc, _)| (center - Vec3::new(bc.x, 0.2, bc.z)).length() < 1.2) {
+                continue;
+            }
+            let rnd_color = |rng: &mut SmallRng| Vec3::new(rng.gen(), rng.gen(), rng.gen());
+            if rng.gen::<f64>() < light_fraction {
+                // Warm-to-white glow of varying strength.
+                let warmth: f64 = rng.gen();
+                let color = Vec3::new(1.0, 0.55 + 0.4 * warmth, 0.2 + 0.6 * warmth);
+                b.sphere_light(center, 0.2, color * (6.0 + 10.0 * rng.gen::<f64>()));
+                continue;
+            }
+            let choice: f64 = rng.gen();
+            let material = if choice < 0.75 {
+                let albedo = rnd_color(&mut rng).mul_elem(rnd_color(&mut rng));
+                Material::Diffuse { albedo }
+            } else if choice < 0.9 {
+                Material::Mirror { albedo: Vec3::splat(0.5) + 0.5 * rnd_color(&mut rng) }
+            } else {
+                glass(Vec3::ONE)
+            };
+            b.sphere(center, 0.2, material);
+        }
+    }
+}
+
+fn spheres() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(13.0, 2.0, 3.0),
+            lookat: Vec3::ZERO,
+            vfov: 20.0,
+        },
+        Background::Sky,
+    );
+    sphere_field(&mut b, 0.0);
+    b.build()
+}
+
+fn night() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(13.0, 2.0, 3.0),
+            lookat: Vec3::ZERO,
+            vfov: 20.0,
+        },
+        Background::Color(Vec3::new(0.004, 0.006, 0.015)),
+    );
+    sphere_field(&mut b, 0.12);
+    b.build()
+}
+
+fn classic() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(278.0, 278.0, -800.0),
+            lookat: Vec3::new(278.0, 278.0, 0.0),
+            vfov: 40.0,
+        },
+        Background::Color(Vec3::ZERO),
+    );
+
+    let white = diffuse(0.73, 0.73, 0.73);
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), white);
+    b.plane(Vec3::new(0.0, 555.0, 0.0), Vec3::new(0.0, -1.0, 0.0), white);
+    b.plane(Vec3::new(0.0, 0.0, 555.0), Vec3::new(0.0, 0.0, -1.0), white);
+    // The camera looks down +z with +x to its left, so red ends up on the left as usual.
+    b.plane(Vec3::new(555.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), diffuse(0.65, 0.05, 0.05));
+    b.plane(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), diffuse(0.12, 0.45, 0.15));
+
+    // A flat panel just under the ceiling, shining down only.
+    b.rect_light(
+        Vec3::new(213.0, 554.0, 227.0),
+        Vec3::new(343.0, 554.0, 332.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(15.0),
+    );
+
+    b.cuboid(Vec3::new(130.0, 0.0, 65.0), Vec3::new(295.0, 165.0, 230.0), white);
+    b.cuboid(Vec3::new(265.0, 0.0, 295.0), Vec3::new(430.0, 330.0, 460.0), white);
+    b.build()
+}
+
+/// 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 {
+    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.build()
+}

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)
+}