fixes

README.md

@@ -68,7 +68,19 @@ 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:
+
+```bash
+cd rust
+cargo build --release
+./target/release/pathtracer 800 600 100 15 cornell_hq.png
+```
+
+It uses all cores via rayon, shows a progress bar and the region brightness table, but does not write `preview.png` or `render_metrics.png`.
 
 ## Scene Description
 

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

rust/src/main.rs

@@ -11,10 +11,65 @@ use rand::{Rng, SeedableRng};
 use rayon::prelude::*;
 
 use camera::Camera;
-use scene::{cornell_box, Material, Ray, Scene, EPSILON};
+use scene::{cornell_box, BoxLight, HitRecord, Material, Ray, Scene, EPSILON};
 use vec3::Vec3;
 
-fn trace<R: Rng>(scene: &Scene, ray: &Ray, depth: u32, max_depth: u32, rng: &mut R) -> Vec3 {
+const SHADOW_EPSILON: f64 = 1e-3;
+
+/// Direct light from one uniformly sampled point on the box light's surface.
+fn sample_box_light<R: Rng>(scene: &Scene, light: &BoxLight, rec: &HitRecord, albedo: Vec3, rng: &mut R) -> Vec3 {
+    let (lo, hi) = (light.min, light.max);
+    let size = hi - lo;
+    let face_areas = [size.y * size.z, size.x * size.z, size.x * size.y];
+    let half_area: f64 = face_areas.iter().sum();
+    let total_area = 2.0 * half_area;
+
+    // Pick a face with probability proportional to its area, then a uniform
+    // point on it: overall pdf is 1 / total_area.
+    let u = rng.gen::<f64>() * half_area;
+    let axis = if u < face_areas[0] {
+        0
+    } else if u < face_areas[0] + face_areas[1] {
+        1
+    } else {
+        2
+    };
+    let positive_side = rng.gen::<f64>() < 0.5;
+    let mut p = [
+        lo.x + size.x * rng.gen::<f64>(),
+        lo.y + size.y * rng.gen::<f64>(),
+        lo.z + size.z * rng.gen::<f64>(),
+    ];
+    p[axis] = if positive_side { hi[axis] } else { lo[axis] };
+    let mut n = [0.0; 3];
+    n[axis] = if positive_side { 1.0 } else { -1.0 };
+    let (lp, light_normal) = (Vec3::new(p[0], p[1], p[2]), Vec3::new(n[0], n[1], n[2]));
+
+    let to_light = lp - rec.point;
+    let dist_sq = to_light.dot(to_light);
+    if dist_sq < 1e-8 {
+        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 - SHADOW_EPSILON) {
+        return Vec3::ZERO;
+    }
+
+    let brdf = albedo / PI;
+    let geometry = cos_at_rec * cos_at_light / dist_sq;
+    light.emission.mul_elem(brdf) * (geometry * total_area)
+}
+
+/// `count_emission` is false after a diffuse bounce: that light was already
+/// gathered by next event estimation, so hitting it again would double count.
+fn trace<R: Rng>(scene: &Scene, ray: &Ray, depth: u32, max_depth: u32, rng: &mut R, count_emission: bool) -> Vec3 {
     if depth >= max_depth {
         return Vec3::ZERO;
     }
@@ -24,69 +79,24 @@ fn trace<R: Rng>(scene: &Scene, ray: &Ray, depth: u32, max_depth: u32, rng: &mut
         return Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t;
     };
 
-    let emitted = rec.material.emission();
-    if let Material::Light { .. } = rec.material {
-        return emitted;
+    if let Material::Light { emission } = rec.material {
+        return if count_emission { emission } else { Vec3::ZERO };
     }
 
     let Some((attenuation, scattered)) = rec.material.scatter(ray, &rec, rng) else {
-        return emitted;
+        return Vec3::ZERO;
     };
 
-    let indirect = trace(scene, &scattered, depth + 1, max_depth, rng);
-    let mut color = emitted + attenuation.mul_elem(indirect);
-
     if let Material::Diffuse { albedo } = rec.material {
         let mut direct = Vec3::ZERO;
         for light in &scene.lights {
-            let (lo, hi) = (light.min, light.max);
-            for _ in 0..2 {
-                let face = (rng.gen::<f64>() * 6.0) as usize;
-                let x = lo.x + (hi.x - lo.x) * rng.gen::<f64>();
-                let y = lo.y + (hi.y - lo.y) * rng.gen::<f64>();
-                let z = lo.z + (hi.z - lo.z) * rng.gen::<f64>();
-                let lp = match face {
-                    0 => Vec3::new(lo.x, y, z),
-                    1 => Vec3::new(hi.x, y, z),
-                    2 => Vec3::new(x, lo.y, z),
-                    3 => Vec3::new(x, hi.y, z),
-                    4 => Vec3::new(x, y, lo.z),
-                    _ => Vec3::new(x, y, hi.z),
-                };
-
-                let to_light = lp - rec.point;
-                let dist_sq = to_light.dot(to_light);
-                if dist_sq < 1e-8 {
-                    continue;
-                }
-                let dist = dist_sq.sqrt();
-                let d_to_light = to_light / dist;
-                let cos_at_rec = rec.normal.dot(d_to_light);
-                if cos_at_rec <= 0.0 {
-                    continue;
-                }
-
-                if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist - EPSILON) {
-                    continue;
-                }
-
-                let light_area = 6.0 * (hi.x - lo.x) * (hi.y - lo.y) * (hi.z - lo.z);
-                let pdf = light_area / (dist_sq * cos_at_rec);
-                if pdf < 1e-12 {
-                    continue;
-                }
-                let brdf = albedo / PI;
-                direct += light.emission.mul_elem(brdf) * (cos_at_rec / pdf);
-                break;
-            }
-        }
-
-        if direct.dot(direct) > 0.0 {
-            color = emitted + attenuation.mul_elem(indirect + direct);
+            direct += sample_box_light(scene, light, &rec, albedo, rng);
         }
+        let indirect = trace(scene, &scattered, depth + 1, max_depth, rng, false);
+        return direct + attenuation.mul_elem(indirect);
     }
 
-    color
+    attenuation.mul_elem(trace(scene, &scattered, depth + 1, max_depth, rng, true))
 }
 
 struct Settings {
@@ -126,10 +136,10 @@ fn render(scene: &Scene, s: &Settings) -> Vec<Vec3> {
                 let mut rng = SmallRng::seed_from_u64((j * width + i) as u64 ^ 0x9E37_79B9_7F4A_7C15);
                 let mut color = Vec3::ZERO;
                 for _ in 0..s.spp {
-                    let u = (i as f64 + rng.gen::<f64>()) / (width - 1) as f64;
-                    let v = ((height - 1 - j) as f64 + rng.gen::<f64>()) / (height - 1) as f64;
+                    let u = (i as f64 + rng.gen::<f64>()) / width as f64;
+                    let v = ((height - 1 - j) as f64 + rng.gen::<f64>()) / height as f64;
                     let ray = cam.get_ray(u, v);
-                    color += trace(scene, &ray, 0, s.max_depth, &mut rng);
+                    color += trace(scene, &ray, 0, s.max_depth, &mut rng, true);
                 }
                 color /= s.spp as f64;
                 *px = color.map(|c| c.max(0.0).sqrt().min(1.0));
@@ -162,16 +172,14 @@ fn region_mean(image: &[Vec3], width: usize, rows: std::ops::Range<usize>, cols:
 }
 
 fn print_region_table(image: &[Vec3], width: usize, height: usize) {
-    // Same regions as the Python version (which slices the Floor/Ceiling/Center
-    // columns by height-based indices).
     let (h_start, h_end) = (height / 3, 2 * height / 3);
-    let clamp_cols = |r: std::ops::Range<usize>| r.start.min(width)..r.end.min(width);
+    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", 2 * height / 3..height, clamp_cols(h_start..h_end)),
-        ("Ceiling", 0..height / 3, clamp_cols(h_start..h_end)),
-        ("Center", h_start..h_end, clamp_cols(h_start..h_end)),
+        ("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");
@@ -202,8 +210,8 @@ fn main() {
         max_depth: parse(4, "15") as u32,
         output: arg(5, "cornell_box.png"),
     };
-    if settings.width < 2 || settings.height < 2 || settings.spp == 0 {
-        eprintln!("width and height must be >= 2 and samples_per_pixel >= 1");
+    if settings.width == 0 || settings.height == 0 || settings.spp == 0 {
+        eprintln!("width, height and samples_per_pixel must be >= 1");
         std::process::exit(2);
     }
 

rust/src/scene.rs

@@ -5,6 +5,20 @@ use crate::vec3::{reflect, refract, schlick, Vec3};
 
 pub const EPSILON: f64 = 1e-6;
 
+pub fn random_unit_vector<R: Rng>(rng: &mut R) -> Vec3 {
+    loop {
+        let v = Vec3::new(
+            rng.sample(StandardNormal),
+            rng.sample(StandardNormal),
+            rng.sample(StandardNormal),
+        );
+        let n = v.length();
+        if n > 1e-12 {
+            return v / n;
+        }
+    }
+}
+
 #[derive(Clone, Copy, Debug)]
 pub struct Ray {
     pub origin: Vec3,
@@ -53,23 +67,13 @@ impl HitRecord {
 }
 
 impl Material {
-    pub fn emission(&self) -> Vec3 {
-        match self {
-            Material::Light { emission } => *emission,
-            _ => Vec3::ZERO,
-        }
-    }
-
     /// 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 } => {
-                let jitter = Vec3::new(
-                    rng.sample(StandardNormal),
-                    rng.sample(StandardNormal),
-                    rng.sample(StandardNormal),
-                );
-                let mut target = rec.normal + jitter * 0.5;
+                // normal + uniform unit vector is cosine-distributed, so the
+                // Lambertian BRDF * cos / pdf reduces to just the albedo.
+                let mut target = rec.normal + random_unit_vector(rng);
                 if target.length() < 1e-8 {
                     target = rec.normal;
                 }
@@ -168,7 +172,8 @@ impl Shape {
                 if t_lo >= t_max || t_lo <= t_min {
                     return None;
                 }
-                let sign = if inv_ds[hit_axis] > 0.0 { 1.0 } else { -1.0 };
+                // 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])))
@@ -202,7 +207,6 @@ impl Scene {
         self.objects.push(Object { shape, material });
     }
 
-    #[allow(dead_code)] // see cornell_box()
     pub fn add_box_light(&mut self, min: Vec3, max: Vec3, emission: Vec3) {
         self.lights.push(BoxLight { min, max, emission });
         self.add(Shape::Box { min, max }, Material::Light { emission });
@@ -247,16 +251,10 @@ pub fn cornell_box() -> Scene {
     scene.add(plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0)), red);
     scene.add(plane(Vec3::new(500.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0)), green);
 
-    // Like raytracer.py, the light is added as a plain object rather than via
-    // add_box_light, so `lights` stays empty and NEE in `trace` never runs.
-    // Using add_box_light here enables the Python NEE code as-written, which
-    // over-brightens the image (it double-counts light and its pdf is off).
-    scene.add(
-        Shape::Box {
-            min: Vec3::new(0.0, 180.0, 180.0),
-            max: Vec3::new(25.0, 220.0, 320.0),
-        },
-        Material::Light { emission: Vec3::splat(50.0) },
+    scene.add_box_light(
+        Vec3::new(0.0, 180.0, 180.0),
+        Vec3::new(25.0, 220.0, 320.0),
+        Vec3::splat(50.0),
     );
 
     let sphere = |c: Vec3, r: f64| Shape::Sphere { center: c, radius: r };