fog rendering

README.md

@@ -79,6 +79,7 @@ cd rust
 cargo build --release
 ./target/release/pathtracer 800 600 100 15 cornell_hq.png
 ./target/release/pathtracer --scene caustics 800 600 1024 15 caustics.png
+./target/release/pathtracer --scene beam 800 600 1024 15 beam.png
 ./target/release/pathtracer --list-scenes
 ```
 
@@ -93,9 +94,32 @@ It uses all cores via rayon, shows a progress bar and (for the `cornell` scene)
 | `spheres` | Large scene: ~470 random diffuse/mirror/glass spheres under a sky (after *Ray Tracing in One Weekend*). The BVH makes it ~11x faster than testing every object |
 | `night` | The same field at night, lit by 61 glowing spheres. Each shading point picks one light, weighted by power / distance², so many lights cost no more shadow rays than one |
 | `classic` | The classic Cornell box: one-sided ceiling panel light, two boxes, colour bleeding from the walls |
+| `fog` | The `caustics` scene in light fog: lit haze, a glow around the light and caustics under the spheres, from caustic photon maps |
+| `beam` | After the classic photon-mapping image of a volume caustic: a spotlight shines down through a glass sphere in dense fog, and the sphere focuses its beam into a cone that narrows to a focus and widens again onto the floor. Photon mapped; ~6 min at 800x600, 1024 spp |
+| `media` | Media inside objects: jade (glass with a green scattering interior), smoke, amber liquid (pure absorption) and a cloud |
+| `haze` | The `city` under a layer of atmospheric haze that fades the horizon and glows around the sun |
 | `city` | Massive scale: a procedural 400x400-lot city of ~124,000 objects (buildings with stepped tiers, mirror-glass towers, parks with trees) under a low sun. Renders 800x600 at 64 spp in ~1 s; without the BVH it is ~1000x slower |
 
-Scenes are defined in [`rust/src/scenes.rs`](rust/src/scenes.rs) with a small builder (`plane`, `sphere`, `cuboid`, `box_light`, `sphere_light`, `rect_light`). Every emitter must be added with a `*_light` method so light sampling knows about it.
+Scenes are defined in [`rust/src/scenes.rs`](rust/src/scenes.rs) with a small builder (`plane`, `sphere`, `cuboid`, `box_light`, `sphere_light`, `rect_light`, `spot_light`, `fog`, `caustic_photons`, `media_photons`). Every emitter must be added with a `*_light` method so light sampling knows about it.
+
+### Volumetrics
+
+Homogeneous participating media ([`rust/src/medium.rs`](rust/src/medium.rs)), each with a coloured absorption coefficient `sigma_a`, a grey scattering coefficient `sigma_s` and a Henyey-Greenstein anisotropy `g`, can be placed:
+
+- as **fog** filling an axis-aligned region outside all objects (`SceneBuilder::fog`, e.g. a haze layer below some height, so the sky stays visible above it);
+- **inside glass**: `Material::Dielectric { medium: Some(..), .. }` (jade, milk, coloured liquids);
+- inside an **invisible boundary**: `Material::Volume { medium }` (smoke, clouds).
+
+Paths sample a scattering distance from `sigma_s` and apply absorption as an exact `exp(-sigma_a * d)` weight, so coloured absorption needs no spectral tricks. At a scattering point they sample a light like a diffuse surface does (with the phase function instead of the BRDF) and continue in a Henyey-Greenstein direction. Shadow rays are attenuated by the media they cross and pass through volume boundaries. `cargo test --release` runs white-furnace tests (a non-absorbing medium in front of a uniform background must render exactly 1), a Beer-Lambert test and phase function tests. Scenes without media render exactly as before.
+
+### Photon mapping for caustics
+
+Path tracing finds a caustic by starting at the camera: a point in fog or on the floor has to pick a direction that happens to pass through the glass and hit the light. With a small light, and especially a spotlight, that almost never happens, so volumetric caustics come out as fireflies or not at all. Scenes can opt in to photon maps ([`rust/src/photons.rs`](rust/src/photons.rs)): before rendering, photons are shot from the lights, and those that pass through glass or off mirrors are stored where they first land in fog or on a diffuse surface. Camera paths gather them along their fog segments (the beam radiance estimate of Jarosz et al.) and at diffuse surfaces, and path tracing skips exactly those caustic paths so nothing is counted twice.
+
+- `caustic_photons(emitted, volume_radius, surface_radius)`: all lights emit photons; the maps replace caustics at a path's first fog or diffuse interaction.
+- `media_photons(...)`: only spotlights emit photons, and their photons also carry their direct light into the fog (much smoother beams than sampling the spotlight from each fog point). Other lights are path traced as usual.
+
+`--photons N` overrides the photon count and `--no-photons` renders with path tracing alone. Gathering blurs by the radius, so photon-mapped images converge to a slightly blurred result; checked against pure path tracing at 4096-8192 spp, mean brightness agrees to within 0.3%. Spotlights (`spot_light`) emit within a cone with a soft edge, and light selection gives a spotlight zero probability for points outside its cone.
 
 ### Beyond the Python version
 

rust/src/main.rs

@@ -1,5 +1,7 @@
 mod bvh;
 mod camera;
+mod medium;
+mod photons;
 mod scene;
 mod scenes;
 mod vec3;
@@ -13,92 +15,216 @@ use rand::{Rng, SeedableRng};
 use rayon::prelude::*;
 
 use camera::Camera;
-use scene::{HitRecord, Light, Material, Ray, Scene, EPSILON};
+use medium::{hg_phase, sample_hg, MediumState};
+use photons::CausticMaps;
+use scene::{spot_emission, Light, Material, Ray, Scene, EPSILON};
 use vec3::Vec3;
 
 /// Shadow rays stop this fraction short of the light, so they don't hit the
 /// light's own surface. Relative, because scenes range from ~1 to ~1000 units.
 const SHADOW_RELATIVE_EPSILON: f64 = 1e-4;
 
-/// 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 {
+/// Where light sampled by next event estimation arrives.
+enum Receiver {
+    /// A diffuse surface: Lambertian BRDF times the cosine at the surface.
+    Surface { normal: Vec3, albedo: Vec3 },
+    /// A scattering point in a medium: the phase function for turning from
+    /// the path's direction `dir` towards the light.
+    Medium { dir: Vec3, g: f64 },
+}
+
+/// Direct light at `point` from one uniformly sampled point on `light`'s
+/// surface, attenuated by whatever media the shadow ray passes through.
+fn sample_direct<R: Rng>(
+    scene: &Scene,
+    state: MediumState,
+    light: &Light,
+    point: Vec3,
+    receiver: Receiver,
+    rng: &mut R,
+) -> Vec3 {
     let (lp, light_normal) = light.sample(rng);
 
-    let to_light = lp - rec.point;
+    let to_light = lp - point;
     let dist_sq = to_light.dot(to_light);
     if dist_sq < 1e-12 {
         return Vec3::ZERO;
     }
     let dist = dist_sq.sqrt();
     let d_to_light = to_light / dist;
-    let cos_at_rec = rec.normal.dot(d_to_light);
     let cos_at_light = -light_normal.dot(d_to_light);
+    let (response, cos_at_rec, origin) = match receiver {
+        Receiver::Surface { normal, albedo } => (albedo / PI, normal.dot(d_to_light), point + EPSILON * normal),
+        Receiver::Medium { dir, g } => (Vec3::splat(hg_phase(dir.dot(d_to_light), g)), 1.0, point),
+    };
     if cos_at_rec <= 0.0 || cos_at_light <= 0.0 {
         return Vec3::ZERO;
     }
 
-    if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist * (1.0 - SHADOW_RELATIVE_EPSILON)) {
+    let tr = scene.transmittance(state, origin, d_to_light, dist * (1.0 - SHADOW_RELATIVE_EPSILON));
+    if tr == Vec3::ZERO {
         return Vec3::ZERO;
     }
 
     // Uniform area sampling has pdf 1 / area.
-    let brdf = albedo / PI;
     let geometry = cos_at_rec * cos_at_light / dist_sq;
-    light.emission.mul_elem(brdf) * (geometry * light.area())
+    let emission = spot_emission(light.emission, light.spot, cos_at_light);
+    (emission.mul_elem(response) * (geometry * light.area())).mul_elem(tr)
 }
 
 /// Bounces before Russian roulette may terminate a path.
 const RR_MIN_DEPTH: u32 = 3;
+/// Guard against paths that keep crossing invisible volume boundaries.
+const MAX_VOLUME_CROSSINGS: u32 = 256;
 
-fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> Vec3 {
+/// Russian roulette: after RR_MIN_DEPTH bounces, end dim paths early,
+/// boosting survivors so the estimate stays unbiased. Glass/mirror paths
+/// (throughput ~1) survive. Returns false if the path ends.
+fn russian_roulette<R: Rng>(throughput: &mut Vec3, bounces: u32, rng: &mut R) -> bool {
+    if bounces >= RR_MIN_DEPTH {
+        let survive = throughput.x.max(throughput.y).max(throughput.z).min(0.95);
+        if rng.gen::<f64>() >= survive {
+            return false;
+        }
+        *throughput = *throughput / survive;
+    }
+    true
+}
+
+/// With `caustics`, caustic light (light that went through glass or off
+/// mirrors before reaching fog or a diffuse surface) comes from the photon
+/// maps, and path tracing skips the matching paths (from a diffuse or medium
+/// interaction through mirrors/glass to a light):
+/// - normally (all lights emit photons), only at the path's first such
+///   interaction: the maps are gathered along fog segments and at the first
+///   diffuse surface until then, and later caustics are path traced;
+/// - with `direct_in_media` (only spotlights emit photons), for spotlights
+///   at every interaction: the volume map also holds their direct light, so
+///   points in fog don't sample spotlights, and every fog segment and
+///   diffuse surface of the path gathers the maps. Other lights are path
+///   traced as usual.
+fn trace<R: Rng>(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R, caustics: Option<&CausticMaps>) -> Vec3 {
     let mut radiance = Vec3::ZERO;
     // Product of the attenuations along the path so far (divided by the RR
     // survival probabilities), i.e. how much of any light found here reaches the camera.
     let mut throughput = Vec3::ONE;
-    // False after a diffuse bounce: that light was already gathered by next
-    // event estimation, so hitting it again would double count.
+    // False after a diffuse bounce or medium scattering: that light was
+    // already gathered by next event estimation, so hitting it again would
+    // double count.
     let mut count_emission = true;
     let mut ray = *primary;
+    let mut state = MediumState::Outside;
+    let mut bounces = 0;
+    let mut crossings = 0;
+    // Diffuse and medium interactions so far.
+    let mut non_specular = 0;
+    let direct_in_media = caustics.is_some_and(|m| m.direct_in_media);
+    // Whether this interaction gathers from the photon maps.
+    let maps_here = |non_specular: u32| caustics.filter(|_| direct_in_media || non_specular == 0);
 
-    for depth in 0..max_depth {
-        let Some(rec) = scene.hit(&ray, EPSILON, f64::INFINITY) else {
+    while bounces < max_depth {
+        let hit = scene.hit(&ray, EPSILON, f64::INFINITY);
+
+        // Travel through the medium (if any) up to the surface: either
+        // scatter somewhere along the way, or arrive with the absorption
+        // along the way applied.
+        let t_hit = hit.as_ref().map_or(f64::INFINITY, |h| h.t);
+        if let Some((medium, t0, t1)) = scene.medium_segment(state, &ray, t_hit) {
+            if let Some(maps) = maps_here(non_specular) {
+                radiance += throughput.mul_elem(maps.beam(&medium, &ray, t0, t1));
+            }
+            let t = t0 + medium.sample_distance(rng);
+            if t < t1 {
+                throughput = throughput.mul_elem(medium.absorption(t - t0));
+                let p = ray.at(t);
+                // Spotlights are left out when the photon map carries their
+                // direct light into the fog.
+                if let Some((light, pick_prob)) = scene.pick_light(p, direct_in_media, rng) {
+                    let receiver = Receiver::Medium { dir: ray.direction, g: medium.g };
+                    let direct = sample_direct(scene, state, light, p, receiver, rng) / pick_prob;
+                    radiance += throughput.mul_elem(direct);
+                }
+                count_emission = false;
+                non_specular += 1;
+                ray = Ray::from_unit(p, sample_hg(ray.direction, medium.g, rng));
+                bounces += 1;
+                if !russian_roulette(&mut throughput, bounces, rng) {
+                    break;
+                }
+                continue;
+            }
+            throughput = throughput.mul_elem(medium.absorption(t1 - t0));
+        }
+
+        let Some(rec) = hit else {
             radiance += throughput.mul_elem(scene.background.radiance(ray.direction));
             break;
         };
 
-        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);
+        match rec.material {
+            Material::Light { emission, spot } => {
+                // Lights emit from their front side only (matters for one-sided rects).
+                // Reached from a diffuse/medium interaction via mirrors or
+                // glass: a caustic, which the photon maps may already cover.
+                let from_maps = match caustics {
+                    Some(_) if non_specular == 0 => false,
+                    Some(m) if m.direct_in_media => spot.is_some(),
+                    Some(_) => non_specular == 1,
+                    None => false,
+                };
+                if count_emission && rec.front_face && !from_maps {
+                    let emission = spot_emission(emission, spot, -ray.direction.dot(rec.normal));
+                    radiance += throughput.mul_elem(emission);
+                }
+                break;
             }
-            break;
+            Material::Volume { medium } => {
+                // Invisible boundary: carry on in the same direction, now
+                // inside (or back outside) the volume. Not a bounce.
+                state = if rec.front_face { MediumState::Inside(Some(medium)) } else { MediumState::Outside };
+                ray = Ray::from_unit(rec.point, ray.direction);
+                crossings += 1;
+                if crossings > MAX_VOLUME_CROSSINGS {
+                    break;
+                }
+                continue;
+            }
+            _ => {}
         }
 
         let Some((attenuation, scattered)) = rec.material.scatter(&ray, &rec, rng) else {
             break;
         };
 
-        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);
+        match rec.material {
+            Material::Diffuse { albedo } => {
+                if let Some(maps) = maps_here(non_specular) {
+                    radiance += throughput.mul_elem(maps.surface(rec.point, rec.normal, albedo));
+                }
+                if let Some((light, pick_prob)) = scene.pick_light(rec.point, false, rng) {
+                    let receiver = Receiver::Surface { normal: rec.normal, albedo };
+                    let direct = sample_direct(scene, state, light, rec.point, receiver, rng) / pick_prob;
+                    radiance += throughput.mul_elem(direct);
+                }
+                count_emission = false;
+                non_specular += 1;
             }
-            count_emission = false;
-        } else {
-            count_emission = true;
+            Material::Dielectric { medium, .. } => {
+                // Refracted through the surface (rather than reflected):
+                // entering or leaving the object's interior medium.
+                if scattered.direction.dot(rec.normal) < 0.0 {
+                    state = if rec.front_face { MediumState::Inside(medium) } else { MediumState::Outside };
+                }
+                count_emission = true;
+            }
+            _ => count_emission = true,
         }
 
         throughput = throughput.mul_elem(attenuation);
         ray = scattered;
-
-        // 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;
+        bounces += 1;
+        if !russian_roulette(&mut throughput, bounces, rng) {
+            break;
         }
     }
 
@@ -160,6 +286,7 @@ struct Renderer<'a> {
     scene: &'a Scene,
     cam: Camera,
     s: &'a Settings,
+    caustics: Option<&'a CausticMaps>,
 }
 
 impl Renderer<'_> {
@@ -168,7 +295,7 @@ impl Renderer<'_> {
         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);
+            let sample = trace(self.scene, &self.cam.get_ray(u, v), self.s.max_depth, &mut px.rng, self.caustics);
             px.sum += sample;
             px.stats.add(sample.luminance());
         }
@@ -208,7 +335,7 @@ fn bar(len: u64, template: &str) -> ProgressBar {
 }
 
 /// Returns the gamma-encoded image and the total number of samples taken.
-fn render(scene: &Scene, s: &Settings) -> (Vec<Vec3>, u64) {
+fn render(scene: &Scene, s: &Settings, caustics: Option<&CausticMaps>) -> (Vec<Vec3>, u64) {
     let (width, height) = (s.width, s.height);
     let view = scene.view;
     let r = Renderer {
@@ -222,6 +349,7 @@ fn render(scene: &Scene, s: &Settings) -> (Vec<Vec3>, u64) {
             (view.lookfrom - view.lookat).length(),
         ),
         s,
+        caustics,
     };
 
     let mut pixels: Vec<Pixel> = (0..width * height)
@@ -333,7 +461,7 @@ fn print_region_table(image: &[Vec3], width: usize, height: usize) {
     }
 }
 
-const USAGE: &str = "usage: pathtracer [--scene NAME] [--list-scenes] \
+const USAGE: &str = "usage: pathtracer [--scene NAME] [--list-scenes] [--photons N | --no-photons] \
 [width] [height] [samples_per_pixel] [max_bounces] [output] [noise_threshold]";
 
 fn list_scenes() {
@@ -346,6 +474,8 @@ fn list_scenes() {
 fn main() {
     // Pull out the --scene / --list-scenes options; the rest are positional.
     let mut scene_name = String::from("cornell");
+    // Override for the number of caustic photons (0: path trace only).
+    let mut photons_override: Option<usize> = None;
     let mut args: Vec<String> = Vec::new();
     let mut raw = std::env::args().skip(1);
     while let Some(a) = raw.next() {
@@ -359,6 +489,14 @@ fn main() {
             });
         } else if let Some(name) = a.strip_prefix("--scene=") {
             scene_name = name.to_string();
+        } else if a == "--photons" {
+            photons_override = raw.next().and_then(|n| n.parse().ok());
+            if photons_override.is_none() {
+                eprintln!("{USAGE}");
+                std::process::exit(2);
+            }
+        } else if a == "--no-photons" {
+            photons_override = Some(0);
         } else if a == "--help" || a == "-h" {
             println!("{USAGE}");
             list_scenes();
@@ -418,7 +556,29 @@ fn main() {
     println!();
 
     let start = Instant::now();
-    let (image, total_samples) = render(&scene, &settings);
+    let photon_settings = scene.caustic_photons.map(|mut p| {
+        if let Some(n) = photons_override {
+            p.emitted = n;
+        }
+        p
+    });
+    let caustics = match photon_settings {
+        Some(p) if p.emitted > 0 => {
+            let t = Instant::now();
+            let maps = CausticMaps::build(&scene, p);
+            let (volume, surface) = maps.counts();
+            println!(
+                "Caustic photons:   {} emitted, {} stored in media, {} on surfaces ({:.2}s)",
+                p.emitted,
+                volume,
+                surface,
+                t.elapsed().as_secs_f64()
+            );
+            Some(maps)
+        }
+        _ => None,
+    };
+    let (image, total_samples) = render(&scene, &settings, caustics.as_ref());
     let elapsed = start.elapsed().as_secs_f64();
     let pixels = (settings.width * settings.height) as f64;
 
@@ -441,3 +601,91 @@ fn main() {
         print_region_table(&image, settings.width, settings.height);
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use medium::Medium;
+    use scene::{Background, SceneBuilder, View};
+
+    fn view() -> View {
+        View { lookfrom: Vec3::new(0.0, 0.0, -10.0), lookat: Vec3::ZERO, vfov: 40.0 }
+    }
+
+    /// Average of `n` path samples along one ray.
+    fn mean_radiance(scene: &Scene, ray: Ray, n: usize) -> Vec3 {
+        let mut rng = SmallRng::seed_from_u64(1);
+        let mut sum = Vec3::ZERO;
+        for _ in 0..n {
+            sum += trace(scene, &ray, 10_000, &mut rng, None);
+        }
+        sum / n as f64
+    }
+
+    fn assert_close(v: Vec3, expected: Vec3, tol: f64) {
+        for a in 0..3 {
+            assert!((v[a] - expected[a]).abs() < tol, "got {v:?}, expected {expected:?}");
+        }
+    }
+
+    const INTO_CENTER: fn() -> Ray = || Ray::from_unit(Vec3::new(0.0, 0.0, -10.0), Vec3::new(0.0, 0.0, 1.0));
+
+    /// White furnace: with no absorption and a uniform background of 1,
+    /// every path eventually leaves and sees 1, so the mean must be 1 however
+    /// much the medium scatters.
+    #[test]
+    fn furnace_fog() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        b.fog(Vec3::splat(-1.0), Vec3::splat(1.0), Medium::new(Vec3::ZERO, 3.0, 0.6));
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 200_000), Vec3::ONE, 0.01);
+    }
+
+    #[test]
+    fn furnace_volume_sphere() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        let smoke = Medium::new(Vec3::ZERO, 4.0, -0.3);
+        b.sphere(Vec3::ZERO, 1.0, Material::Volume { medium: smoke });
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 200_000), Vec3::ONE, 0.01);
+    }
+
+    #[test]
+    fn furnace_glass_with_scattering_interior() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        let milk = Medium::new(Vec3::ZERO, 5.0, 0.2);
+        b.sphere(Vec3::ZERO, 1.0, Material::Dielectric { ref_idx: 1.5, albedo: Vec3::ONE, medium: Some(milk) });
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 200_000), Vec3::ONE, 0.01);
+    }
+
+    /// Henyey-Greenstein: sampled directions have mean cosine g, and the
+    /// phase function (used by next event estimation) integrates to 1.
+    #[test]
+    fn henyey_greenstein() {
+        let mut rng = SmallRng::seed_from_u64(2);
+        let dir = Vec3::new(0.3, -0.5, 0.8).normalize();
+        for g in [-0.7, 0.0, 0.4, 0.9] {
+            let n = 400_000;
+            let mean_cos: f64 = (0..n).map(|_| sample_hg(dir, g, &mut rng).dot(dir)).sum::<f64>() / n as f64;
+            assert!((mean_cos - g).abs() < 0.005, "g = {g}: mean cos {mean_cos}");
+
+            // Midpoint rule over cos_theta (the phase function doesn't depend on phi).
+            let steps = 200_000;
+            let integral: f64 = (0..steps)
+                .map(|i| {
+                    let c = -1.0 + (i as f64 + 0.5) * 2.0 / steps as f64;
+                    hg_phase(c, g) * 2.0 * PI * (2.0 / steps as f64)
+                })
+                .sum();
+            assert!((integral - 1.0).abs() < 1e-3, "g = {g}: integral {integral}");
+        }
+    }
+
+    /// A purely absorbing slab of thickness 2: Beer's law, exactly, per channel.
+    #[test]
+    fn beer_lambert_absorption() {
+        let mut b = SceneBuilder::new(view(), Background::Color(Vec3::ONE));
+        let sigma_a = Vec3::new(0.1, 0.5, 1.5);
+        b.fog(Vec3::splat(-1.0), Vec3::splat(1.0), Medium::new(sigma_a, 0.0, 0.0));
+        let expected = (sigma_a * -2.0).map(f64::exp);
+        assert_close(mean_radiance(&b.build(), INTO_CENTER(), 10), expected, 1e-12);
+    }
+}

rust/src/medium.rs

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

rust/src/photons.rs

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

rust/src/scene.rs

@@ -1,25 +1,32 @@
 use rand::Rng;
 
 use crate::bvh::{Aabb, Bvh};
+use crate::medium::{Fog, Medium, MediumState};
 use crate::vec3::{reflect, refract, schlick, Vec3};
 
 pub const EPSILON: f64 = 1e-6;
 
-/// 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.
+/// Two unit vectors that form an orthonormal basis with the unit vector `n`,
+/// using the branchless construction of Duff et al. 2017 ("Building an
+/// Orthonormal Basis, Revisited").
+pub fn orthonormal_basis(n: Vec3) -> (Vec3, Vec3) {
+    let sign = 1.0_f64.copysign(n.z);
+    let a = -1.0 / (sign + n.z);
+    let b = n.x * n.y * a;
+    let t = Vec3::new(1.0 + sign * n.x * n.x * a, sign * b, -sign * n.x);
+    let bt = Vec3::new(b, sign + n.y * n.y * a, -n.y);
+    (t, bt)
+}
+
+/// Cosine-weighted direction in the hemisphere around the unit vector `n`;
+/// the result is already unit length.
 pub fn cosine_hemisphere<R: Rng>(n: Vec3, rng: &mut R) -> Vec3 {
     let r1: f64 = rng.gen();
     let r2: f64 = rng.gen();
     let phi = 2.0 * std::f64::consts::PI * r1;
     let r = r2.sqrt();
     let (x, y, z) = (r * phi.cos(), r * phi.sin(), (1.0 - r2).max(0.0).sqrt());
-
-    let 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);
+    let (t, bt) = orthonormal_basis(n);
     t * x + bt * y + n * z
 }
 
@@ -51,11 +58,17 @@ impl Ray {
 pub enum Material {
     Diffuse { albedo: Vec3 },
     Mirror { albedo: Vec3 },
-    Dielectric { ref_idx: f64, albedo: Vec3 },
-    Light { emission: Vec3 },
+    /// Glass-like surface; `medium` fills its interior (None: clear).
+    Dielectric { ref_idx: f64, albedo: Vec3, medium: Option<Medium> },
+    /// `spot` restricts emission to a cone around the surface normal.
+    Light { emission: Vec3, spot: Option<Spot> },
+    /// Invisible boundary of a volume (smoke, cloud): rays pass straight
+    /// through, but inside it they travel through `medium`.
+    Volume { medium: Medium },
 }
 
 pub struct HitRecord {
+    pub t: f64,
     pub point: Vec3,
     pub normal: Vec3,
     pub front_face: bool,
@@ -67,6 +80,7 @@ impl HitRecord {
     fn new(r: &Ray, t: f64, outward: Vec3, material: Material) -> Self {
         let front_face = r.direction.dot(outward) < 0.0;
         HitRecord {
+            t,
             point: r.at(t),
             normal: if front_face { outward } else { -outward },
             front_face,
@@ -92,7 +106,7 @@ impl Material {
                 }
                 Some((albedo, Ray::new(rec.point, reflected)))
             }
-            Material::Dielectric { ref_idx, albedo } => {
+            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);
@@ -106,6 +120,7 @@ impl Material {
                 Some((albedo, Ray::new(rec.point, direction)))
             }
             Material::Light { .. } => None,
+            Material::Volume { .. } => Some((Vec3::ONE, Ray::from_unit(rec.point, ray_in.direction))),
         }
     }
 }
@@ -230,11 +245,46 @@ pub struct Object {
     pub material: Material,
 }
 
+/// Emission limited to a cone around a light's surface normal, fading out
+/// smoothly between the inner and outer angle.
+#[derive(Clone, Copy, Debug)]
+pub struct Spot {
+    outer: f64, // radians
+    cos_outer: f64,
+    cos_inner: f64,
+}
+
+impl Spot {
+    pub fn new(outer_deg: f64, inner_deg: f64) -> Self {
+        Spot {
+            outer: outer_deg.to_radians(),
+            cos_outer: outer_deg.to_radians().cos(),
+            cos_inner: inner_deg.to_radians().cos(),
+        }
+    }
+
+    /// Fraction of the full emission leaving at an angle with cosine `cos`
+    /// to the normal.
+    pub fn profile(&self, cos: f64) -> f64 {
+        let t = ((cos - self.cos_outer) / (self.cos_inner - self.cos_outer)).clamp(0.0, 1.0);
+        t * t * (3.0 - 2.0 * t)
+    }
+}
+
+/// Emitted radiance leaving at an angle with cosine `cos` to the normal.
+pub fn spot_emission(emission: Vec3, spot: Option<Spot>, cos: f64) -> Vec3 {
+    match spot {
+        Some(s) => emission * s.profile(cos),
+        None => emission,
+    }
+}
+
 /// An emitter that next event estimation can sample: a box, sphere or rectangle.
 #[derive(Clone, Copy, Debug)]
 pub struct Light {
     shape: Shape,
     pub emission: Vec3,
+    pub spot: Option<Spot>,
     /// Total emitted power (up to a constant): luminance * area.
     power: f64,
     center: Vec3,
@@ -244,11 +294,12 @@ pub struct Light {
 }
 
 impl Light {
-    fn new(shape: Shape, emission: Vec3) -> Self {
+    fn new(shape: Shape, emission: Vec3, spot: Option<Spot>) -> Self {
         let bounds = shape.bounds().expect("planes cannot be lights");
         let mut light = Light {
             shape,
             emission,
+            spot,
             power: 0.0,
             center: bounds.centroid(),
             extent: 0.5 * (bounds.max - bounds.min).length(),
@@ -278,11 +329,50 @@ impl Light {
         }
     }
 
+    /// A photon leaving the light: its origin, direction and power (the
+    /// light's whole emitted power; the caller divides by the photon count).
+    pub fn emit<R: Rng>(&self, rng: &mut R) -> (Vec3, Vec3, Vec3) {
+        let (p, n) = self.sample(rng);
+        match self.spot {
+            None => {
+                // Cosine-weighted directions: power = L * area * pi.
+                let dir = cosine_hemisphere(n, rng);
+                (p, dir, self.emission * (self.area() * std::f64::consts::PI))
+            }
+            Some(spot) => {
+                // Uniform directions within the outer cone, weighted by the
+                // profile and the cosine over the direction pdf.
+                let cos = 1.0 - rng.gen::<f64>() * (1.0 - spot.cos_outer);
+                let sin = (1.0 - cos * cos).max(0.0).sqrt();
+                let phi = 2.0 * std::f64::consts::PI * rng.gen::<f64>();
+                let (t, b) = orthonormal_basis(n);
+                let dir = t * (sin * phi.cos()) + b * (sin * phi.sin()) + n * cos;
+                let pdf = 1.0 / (2.0 * std::f64::consts::PI * (1.0 - spot.cos_outer));
+                (p, dir, self.emission * (spot.profile(cos) * cos * self.area() / pdf))
+            }
+        }
+    }
+
     /// Rough estimate of how much this light contributes at `point`: its
     /// power over the squared distance, clamped so points near or inside the
     /// light's bounds don't blow up.
+    ///
+    /// A spotlight gets exactly zero outside its cone: widened by the
+    /// light's angular radius seen from `point`, so that no point on the
+    /// light could reach it (otherwise light selection keeps picking a
+    /// powerful spotlight for points it can't light, and every other light
+    /// gets rare, heavily weighted samples: speckle noise).
     fn importance(&self, point: Vec3) -> f64 {
         let d = point - self.center;
+        if let (Some(spot), Shape::Rect { normal, .. }) = (self.spot, self.shape) {
+            let dist = d.length();
+            if dist > self.extent {
+                let angle = (normal.dot(d) / dist).clamp(-1.0, 1.0).acos();
+                if angle > spot.outer + (self.extent / dist).asin() {
+                    return 0.0;
+                }
+            }
+        }
         self.power / d.dot(d).max(self.extent * self.extent)
     }
 
@@ -368,6 +458,19 @@ pub struct SceneBuilder {
     lights: Vec<Light>,
     background: Background,
     view: View,
+    fog: Option<Fog>,
+    caustic_photons: Option<CausticPhotons>,
+}
+
+/// Photon mapping settings: photons emitted, the gathering radius in fog and
+/// on surfaces (larger: smoother but blurrier), and whether the fog's map
+/// also carries the spotlights' direct light.
+#[derive(Clone, Copy, Debug)]
+pub struct CausticPhotons {
+    pub emitted: usize,
+    pub volume_radius: f64,
+    pub surface_radius: f64,
+    pub direct_in_media: bool,
 }
 
 impl SceneBuilder {
@@ -378,9 +481,17 @@ impl SceneBuilder {
             lights: Vec::new(),
             background,
             view,
+            fog: None,
+            caustic_photons: None,
         }
     }
 
+    /// Fills the region min..max (bounds may be infinite) with `medium`,
+    /// everywhere outside objects.
+    pub fn fog(&mut self, min: Vec3, max: Vec3, medium: Medium) {
+        self.fog = Some(Fog { region: Aabb { min, max }, medium });
+    }
+
     pub fn add(&mut self, shape: Shape, material: Material) {
         let obj = Object { shape, material };
         match shape {
@@ -404,21 +515,42 @@ impl SceneBuilder {
     /// 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);
+        self.add_light(Shape::Box { min, max }, emission, None);
     }
 
     pub fn sphere_light(&mut self, center: Vec3, radius: f64, emission: Vec3) {
-        self.add_light(Shape::Sphere { center, radius }, emission);
+        self.add_light(Shape::Sphere { center, radius }, emission, None);
     }
 
     /// One-sided rectangular light emitting towards `normal` (±x, ±y or ±z).
     pub fn rect_light(&mut self, min: Vec3, max: Vec3, normal: Vec3, emission: Vec3) {
-        self.add_light(Shape::Rect { min, max, normal }, emission);
+        self.add_light(Shape::Rect { min, max, normal }, emission, None);
     }
 
-    fn add_light(&mut self, shape: Shape, emission: Vec3) {
-        self.lights.push(Light::new(shape, emission));
-        self.add(shape, Material::Light { emission });
+    /// A rectangular spotlight: like rect_light, but emitting only within
+    /// `outer_deg` of `normal`, at full strength within `inner_deg`.
+    pub fn spot_light(&mut self, min: Vec3, max: Vec3, normal: Vec3, emission: Vec3, outer_deg: f64, inner_deg: f64) {
+        self.add_light(Shape::Rect { min, max, normal }, emission, Some(Spot::new(outer_deg, inner_deg)));
+    }
+
+    /// Shoot photons to render caustics (light through glass or off mirrors
+    /// onto fog and diffuse surfaces) from photon maps instead of by path
+    /// tracing alone. See photons.rs.
+    pub fn caustic_photons(&mut self, emitted: usize, volume_radius: f64, surface_radius: f64) {
+        self.caustic_photons = Some(CausticPhotons { emitted, volume_radius, surface_radius, direct_in_media: false });
+    }
+
+    /// Photons for the spotlights only, carrying both their caustics and
+    /// their direct light into fog (replacing light sampling of spotlights
+    /// at points in fog): much smoother beams, at the cost of storing (and
+    /// gathering) many more photons. Other lights are fully path traced.
+    pub fn media_photons(&mut self, emitted: usize, volume_radius: f64, surface_radius: f64) {
+        self.caustic_photons = Some(CausticPhotons { emitted, volume_radius, surface_radius, direct_in_media: true });
+    }
+
+    fn add_light(&mut self, shape: Shape, emission: Vec3, spot: Option<Spot>) {
+        self.lights.push(Light::new(shape, emission, spot));
+        self.add(shape, Material::Light { emission, spot });
     }
 
     pub fn build(self) -> Scene {
@@ -431,6 +563,11 @@ impl SceneBuilder {
                 _ => None,
             })
             .collect();
+        let has_volume_boundaries = self
+            .objects
+            .iter()
+            .chain(&self.planes)
+            .any(|o| matches!(o.material, Material::Volume { .. }));
         Scene {
             bvh: Bvh::build(&bounds, &spheres),
             objects: self.objects,
@@ -438,6 +575,9 @@ impl SceneBuilder {
             lights: self.lights,
             background: self.background,
             view: self.view,
+            fog: self.fog,
+            has_volume_boundaries,
+            caustic_photons: self.caustic_photons,
         }
     }
 }
@@ -449,6 +589,11 @@ pub struct Scene {
     lights: Vec<Light>,
     pub background: Background,
     pub view: View,
+    fog: Option<Fog>,
+    /// Whether shadow rays may need to pass through invisible volume
+    /// boundaries (otherwise any hit blocks them).
+    has_volume_boundaries: bool,
+    pub caustic_photons: Option<CausticPhotons>,
 }
 
 impl Scene {
@@ -469,30 +614,60 @@ impl Scene {
     }
 
     /// 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))
-            }
+    /// its estimated contribution there (power / distance^2, zero outside a
+    /// spotlight's cone), so a scene with many lights costs one shadow ray
+    /// per shading point and mostly samples lights that matter. Spotlights
+    /// are left out with `exclude_spots` (when photons carry their light).
+    /// Returns the light and the probability of picking it. This is
+    /// O(lights) per call; thousands of lights would want a light tree.
+    pub fn pick_light<R: Rng>(&self, point: Vec3, exclude_spots: bool, rng: &mut R) -> Option<(&Light, f64)> {
+        let weight = |l: &Light| if exclude_spots && l.spot.is_some() { 0.0 } else { l.importance(point) };
+        if self.lights.len() == 1 {
+            let only = &self.lights[0];
+            return (weight(only) > 0.0).then_some((only, 1.0));
         }
+        let total: f64 = self.lights.iter().map(weight).sum();
+        if total <= 0.0 {
+            return None;
+        }
+        let u = rng.gen::<f64>() * total;
+        let mut acc = 0.0;
+        let mut last = None;
+        for light in &self.lights {
+            let w = weight(light);
+            if w <= 0.0 {
+                continue;
+            }
+            acc += w;
+            if u < acc {
+                return Some((light, w / total));
+            }
+            last = Some((light, w / total));
+        }
+        // Rounding left u just past the last bucket.
+        last
+    }
+
+    /// Picks a light to emit a photon from (only spotlights if
+    /// `spots_only`), proportionally to its power. Returns the light and the
+    /// probability of picking it.
+    pub fn pick_emitting_light<R: Rng>(&self, spots_only: bool, rng: &mut R) -> Option<(&Light, f64)> {
+        let eligible = |l: &&Light| !spots_only || l.spot.is_some();
+        let total: f64 = self.lights.iter().filter(eligible).map(|l| l.power).sum();
+        if total <= 0.0 {
+            return None;
+        }
+        let u = rng.gen::<f64>() * total;
+        let mut acc = 0.0;
+        let mut last = None;
+        for light in self.lights.iter().filter(eligible) {
+            acc += light.power;
+            if u < acc {
+                return Some((light, light.power / total));
+            }
+            last = Some(light);
+        }
+        last.map(|l| (l, l.power / total))
     }
 
     pub fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option<HitRecord> {
@@ -523,6 +698,60 @@ impl Scene {
         closest.map(|(t, outward, material)| HitRecord::new(r, t, outward, material))
     }
 
+    /// The medium along the ray segment [0, t_end) for a path in `state`, and
+    /// the part of the segment it covers.
+    pub fn medium_segment(&self, state: MediumState, r: &Ray, t_end: f64) -> Option<(Medium, f64, f64)> {
+        match state {
+            MediumState::Outside => {
+                let fog = self.fog.as_ref()?;
+                let (t0, t1) = fog.clip(r.origin, r.direction, 0.0, t_end)?;
+                Some((fog.medium, t0, t1))
+            }
+            MediumState::Inside(medium) => medium.map(|m| (m, 0.0, t_end)),
+        }
+    }
+
+    fn segment_transmittance(&self, state: MediumState, r: &Ray, t_end: f64) -> Vec3 {
+        match self.medium_segment(state, r, t_end) {
+            Some((m, t0, t1)) => m.transmittance(t1 - t0),
+            None => Vec3::ONE,
+        }
+    }
+
+    /// Fraction of light that travels from `origin` a distance `dist` along
+    /// `direction` for a path in `state`: zero if a surface is in the way,
+    /// otherwise the transmittance of the media passed through, including
+    /// any invisible volume boundaries crossed.
+    pub fn transmittance(&self, state: MediumState, origin: Vec3, direction: Vec3, dist: f64) -> Vec3 {
+        if !self.has_volume_boundaries {
+            // Any hit blocks the ray, so it never leaves its starting medium.
+            if self.blocked(origin, direction, dist) {
+                return Vec3::ZERO;
+            }
+            return self.segment_transmittance(state, &Ray::from_unit(origin, direction), dist);
+        }
+
+        let (mut origin, mut remaining, mut state) = (origin, dist, state);
+        let mut tr = Vec3::ONE;
+        // Walk from boundary to boundary; the bound only guards against
+        // pathological scenes.
+        for _ in 0..64 {
+            let r = Ray::from_unit(origin, direction);
+            let hit = self.hit(&r, EPSILON, remaining);
+            tr = tr.mul_elem(self.segment_transmittance(state, &r, hit.as_ref().map_or(remaining, |h| h.t)));
+            let Some(rec) = hit else {
+                return tr;
+            };
+            let Material::Volume { medium } = rec.material else {
+                return Vec3::ZERO;
+            };
+            state = if rec.front_face { MediumState::Inside(Some(medium)) } else { MediumState::Outside };
+            origin = rec.point;
+            remaining -= rec.t;
+        }
+        Vec3::ZERO
+    }
+
     pub fn blocked(&self, origin: Vec3, direction: Vec3, t_max: f64) -> bool {
         let r = Ray::from_unit(origin, direction);
         self.planes.iter().any(|obj| obj.shape.hit(&r, EPSILON, t_max).is_some())

rust/src/scenes.rs

@@ -3,6 +3,7 @@
 use rand::rngs::SmallRng;
 use rand::{Rng, SeedableRng};
 
+use crate::medium::Medium;
 use crate::scene::{Background, Material, Scene, SceneBuilder, View};
 use crate::vec3::Vec3;
 
@@ -43,6 +44,26 @@ pub const SCENES: &[SceneInfo] = &[
         description: "Massive scale: a procedural city of ~125,000 objects under a low sun",
         build: city,
     },
+    SceneInfo {
+        name: "fog",
+        description: "The caustics scene in light fog: beams of focused light and shadow shafts",
+        build: fog,
+    },
+    SceneInfo {
+        name: "media",
+        description: "Participating media inside objects: jade, smoke, amber liquid and cloud",
+        build: media,
+    },
+    SceneInfo {
+        name: "beam",
+        description: "A spotlight through a glass sphere in dense fog: a volumetric caustic (photon mapped)",
+        build: beam,
+    },
+    SceneInfo {
+        name: "haze",
+        description: "The city under a layer of atmospheric haze, with a glow around the sun",
+        build: haze,
+    },
 ];
 
 pub fn find(name: &str) -> Option<&'static SceneInfo> {
@@ -60,7 +81,7 @@ fn mirror(r: f64, g: f64, b: f64) -> Material {
 /// 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 }
+    Material::Dielectric { ref_idx: 1.5, albedo: tint, medium: None }
 }
 
 fn cornell() -> Scene {
@@ -89,6 +110,12 @@ fn cornell() -> Scene {
 }
 
 fn caustics() -> Scene {
+    caustics_room(30.0).build()
+}
+
+/// The caustics scene with a square light of the given half-size (same
+/// total power whatever the size).
+fn caustics_room(light_half_size: f64) -> SceneBuilder {
     let mut b = SceneBuilder::new(
         View {
             lookfrom: Vec3::new(0.0, 160.0, -280.0),
@@ -106,17 +133,134 @@ fn caustics() -> Scene {
     // ~48 degrees from vertical that focus lands on the floor just in front
     // of each sphere, facing the camera. A smaller panel gives sharper
     // caustics but makes them slower to converge.
+    let (h, z) = (light_half_size, 360.0);
     b.rect_light(
-        Vec3::new(-30.0, 250.0, 330.0),
-        Vec3::new(30.0, 250.0, 390.0),
+        Vec3::new(-h, 250.0, z - h),
+        Vec3::new(h, 250.0, z + h),
         Vec3::new(0.0, -1.0, 0.0),
-        Vec3::splat(90.0),
+        Vec3::splat(90.0 * (30.0 * 30.0) / (h * h)),
     );
 
     b.sphere(Vec3::new(-75.0, 45.0, 110.0), 45.0, glass(Vec3::ONE));
     b.sphere(Vec3::new(40.0, 38.0, 80.0), 38.0, glass(Vec3::new(1.0, 0.7, 0.3)));
     b.sphere(Vec3::new(125.0, 30.0, 150.0), 30.0, glass(Vec3::new(0.4, 0.65, 1.0)));
     b.sphere(Vec3::new(-10.0, 35.0, 210.0), 35.0, mirror(0.9, 0.9, 0.9));
+    b
+}
+
+fn fog() -> Scene {
+    // A light 4x the area (same power): light that reaches it through the
+    // glass spheres after scattering in the fog (volumetric caustics, the
+    // beams) is found 4x as often at a quarter of the brightness, which cuts
+    // the fireflies those rare paths cause.
+    let mut b = caustics_room(60.0);
+    // Light fog (mean free path ~650 units, about the size of the room).
+    // Forward scattering, so the beams stand out when looking towards the
+    // light, as the camera does. It stops well below the light (at y = 250):
+    // sampling the light from fog right next to it divides by a tiny squared
+    // distance and turns into fireflies.
+    b.fog(
+        Vec3::new(f64::NEG_INFINITY, 0.0, f64::NEG_INFINITY),
+        Vec3::new(f64::INFINITY, 190.0, f64::INFINITY),
+        Medium::new(Vec3::splat(0.0001), 0.0015, 0.5),
+    );
+    // The caustics in the fog and on the floor come from photon maps.
+    b.caustic_photons(4_000_000, 4.0, 3.0);
+    b.build()
+}
+
+/// After the classic photon-mapping image of a volume caustic: a spotlight
+/// shines down through a glass sphere in a foggy room, and the sphere
+/// focuses its beam into a cone of light that narrows to a focus and widens
+/// again towards the floor.
+fn beam() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.0, 200.0, -420.0),
+            lookat: Vec3::new(0.0, 190.0, 150.0),
+            vfov: 50.0,
+        },
+        Background::Color(Vec3::ZERO),
+    );
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.6, 0.6, 0.6));
+    b.plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.5, 0.5, 0.5));
+
+    // A small spotlight pointing straight down; its 12 degree cone is
+    // narrower than the sphere, so the whole beam goes through the glass.
+    b.spot_light(
+        Vec3::new(-8.0, 400.0, 142.0),
+        Vec3::new(8.0, 400.0, 158.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::new(1.0, 0.95, 0.85) * 10000.0,
+        12.0,
+        8.0,
+    );
+    // A large, dim, cool fill light in the ceiling, so the room and the
+    // glass show.
+    b.rect_light(
+        Vec3::new(-400.0, 420.0, -300.0),
+        Vec3::new(400.0, 420.0, 500.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::new(0.15, 0.165, 0.21),
+    );
+    // A glass sphere focuses light about 1.5 radii from its centre: here
+    // at y = 130, halfway between the sphere and the floor.
+    b.sphere(Vec3::new(0.0, 220.0, 150.0), 60.0, glass(Vec3::ONE));
+
+    // Dense fog (mean free path ~330, twice as dense as the fog scene)
+    // below the light, scattering evenly in all directions so the beam
+    // shows from the side.
+    b.fog(
+        Vec3::new(f64::NEG_INFINITY, 0.0, f64::NEG_INFINITY),
+        Vec3::new(f64::INFINITY, 390.0, f64::INFINITY),
+        Medium::new(Vec3::splat(0.0003), 0.003, 0.0),
+    );
+    // Photons carry both the direct beam and the focused one through the
+    // fog: path tracing alone leaves a small spotlight's beam grainy.
+    b.media_photons(600_000, 3.5, 3.0);
+    b.build()
+}
+
+fn media() -> Scene {
+    let mut b = SceneBuilder::new(
+        View {
+            lookfrom: Vec3::new(0.0, 150.0, -330.0),
+            lookat: Vec3::new(0.0, 40.0, 100.0),
+            vfov: 36.0,
+        },
+        Background::Color(Vec3::new(0.02, 0.02, 0.025)),
+    );
+    b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.75, 0.75, 0.75));
+    b.plane(Vec3::new(0.0, 0.0, 400.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.55, 0.55, 0.6));
+    // A large soft key light above and slightly behind, plus a dim fill from the left.
+    b.rect_light(
+        Vec3::new(-120.0, 320.0, 60.0),
+        Vec3::new(120.0, 320.0, 220.0),
+        Vec3::new(0.0, -1.0, 0.0),
+        Vec3::splat(7.0),
+    );
+    b.rect_light(
+        Vec3::new(400.0, 40.0, -60.0),
+        Vec3::new(400.0, 200.0, 140.0),
+        Vec3::new(-1.0, 0.0, 0.0),
+        Vec3::new(1.2, 1.3, 1.6),
+    );
+
+    let r = 42.0;
+    // Jade: glass filled with a dense, strongly scattering medium that
+    // absorbs red and blue, so light diffuses through it and comes out green.
+    let jade = Medium::new(Vec3::new(0.030, 0.003, 0.018), 0.2, 0.3);
+    b.sphere(Vec3::new(-165.0, r, 120.0), r, Material::Dielectric { ref_idx: 1.5, albedo: Vec3::ONE, medium: Some(jade) });
+    // Smoke: no surface at all, a grey absorbing and scattering volume.
+    let smoke = Medium::new(Vec3::splat(0.008), 0.018, 0.0);
+    b.sphere(Vec3::new(-55.0, r, 120.0), r, Material::Volume { medium: smoke });
+    // Amber liquid: clear glass holding a purely absorbing medium that
+    // absorbs blue most (Beer's law), so thicker parts look deeper amber.
+    let amber = Medium::new(Vec3::new(0.002, 0.011, 0.035), 0.0, 0.0);
+    b.sphere(Vec3::new(55.0, r, 120.0), r, Material::Dielectric { ref_idx: 1.33, albedo: Vec3::ONE, medium: Some(amber) });
+    // Cloud: dense, non-absorbing and strongly forward scattering.
+    let cloud = Medium::new(Vec3::ZERO, 0.15, 0.8);
+    b.sphere(Vec3::new(165.0, r, 120.0), r, Material::Volume { medium: cloud });
     b.build()
 }
 
@@ -228,6 +372,23 @@ fn classic() -> Scene {
 /// tiers on the tall ones) that rise towards downtown, some mirror-glass
 /// towers, and parks with sphere trees. Lit by a low sun and the sky.
 fn city() -> Scene {
+    city_builder().build()
+}
+
+fn haze() -> Scene {
+    let mut b = city_builder();
+    // Haze in the lowest 400 units of air, thin enough to see across a few
+    // kilometres (mean free path 2500) but enough to fade the horizon.
+    // Forward scattering puts a glow around the sun.
+    b.fog(
+        Vec3::new(f64::NEG_INFINITY, 0.0, f64::NEG_INFINITY),
+        Vec3::new(f64::INFINITY, 400.0, f64::INFINITY),
+        Medium::new(Vec3::splat(0.00004), 0.0004, 0.7),
+    );
+    b.build()
+}
+
+fn city_builder() -> SceneBuilder {
     const GRID: usize = 400;
     const LOT: f64 = 10.0;
     // Every BLOCK-th row/column of lots is left empty as an avenue.
@@ -322,5 +483,5 @@ fn city() -> Scene {
             }
         }
     }
-    b.build()
+    b
 }