diff --git a/README.md b/README.md index 65a6f45..e262fe3 100644 --- a/README.md +++ b/README.md @@ -72,17 +72,34 @@ Multiprocessing uses one worker per CPU core by default. Pass `workers=N` to the ## Rust version -A much faster port lives in [`rust/`](rust/) (renders 800x600 at 100 spp in seconds). It takes the same arguments: +A much faster port lives in [`rust/`](rust/) (renders 800x600 at 100 spp in seconds). It takes the same arguments, plus an optional `--scene`: ```bash cd rust cargo build --release ./target/release/pathtracer 800 600 100 15 cornell_hq.png +./target/release/pathtracer --scene caustics 800 600 1024 15 caustics.png +./target/release/pathtracer --list-scenes ``` -It uses all cores via rayon, shows a progress bar and the region brightness table, but does not write `preview.png` or `render_metrics.png`. +It uses all cores via rayon, shows a progress bar and (for the `cornell` scene) the region brightness table, but does not write `preview.png` or `render_metrics.png`. Without an output name it writes `.png`. -Beyond the Python version it adds: +### Scenes + +| Scene | What it shows | +|-------|---------------| +| `cornell` (default) | The scene from `raytracer.py` | +| `caustics` | Clear, amber and blue glass spheres and a mirror under one light, placed so each sphere focuses a sharp caustic into its shadow. Caustics are the hardest light paths for a path tracer, so use ~1000+ spp | +| `spheres` | Large scene: ~470 random diffuse/mirror/glass spheres under a sky (after *Ray Tracing in One Weekend*). The BVH makes it ~11x faster than testing every object | +| `night` | The same field at night, lit by 61 glowing spheres. Each shading point picks one light, weighted by power / distance², so many lights cost no more shadow rays than one | +| `classic` | The classic Cornell box: one-sided ceiling panel light, two boxes, colour bleeding from the walls | + +Scenes are defined in [`rust/src/scenes.rs`](rust/src/scenes.rs) with a small builder (`plane`, `sphere`, `cuboid`, `box_light`, `sphere_light`, `rect_light`). Every emitter must be added with a `*_light` method so light sampling knows about it. + +### Beyond the Python version + +- **BVH** - bounded objects live in a bounding volume hierarchy (median split); planes are tested separately. +- **Many-light sampling** - box, sphere and one-sided rectangle lights, one picked per shading point in proportion to its estimated contribution. - **Russian roulette** - after 3 bounces, dim paths are terminated randomly (survivors are re-weighted, so the image stays unbiased). About 2.5x more image quality per second of rendering in this scene. - **Cheap cosine-weighted diffuse sampling** - two uniform random numbers and an orthonormal basis instead of Gaussian samples. @@ -108,7 +125,7 @@ The scene is the classic **Cornell Box**: - Position: (278, 278, -800) - Look-at: (278, 278, 278) (box center) -- FOV: 40° +- FOV: 40° vertical (earlier versions effectively used half of this, zooming the image in 2x) - Focus distance: 1078 ## Materials diff --git a/raytracer.py b/raytracer.py index 5c2c1cf..a2d5120 100644 --- a/raytracer.py +++ b/raytracer.py @@ -283,8 +283,8 @@ class Camera: w = normalize(self.origin - np.asarray(lookat, dtype=np.float64)) u = normalize(np.cross(vup, w)) v = np.cross(w, u) - self.horizontal = focus_dist * half_width * u - self.vertical = focus_dist * half_height * v + self.horizontal = 2.0 * focus_dist * half_width * u + self.vertical = 2.0 * focus_dist * half_height * v self.lower_left_corner = self.origin - self.horizontal / 2.0 - self.vertical / 2.0 - focus_dist * w self.aperture = aperture diff --git a/rust/caustics.png b/rust/caustics.png new file mode 100644 index 0000000..39a0a7d Binary files /dev/null and b/rust/caustics.png differ diff --git a/rust/classic.png b/rust/classic.png new file mode 100644 index 0000000..73722cc Binary files /dev/null and b/rust/classic.png differ diff --git a/rust/night.png b/rust/night.png new file mode 100644 index 0000000..21de2ce Binary files /dev/null and b/rust/night.png differ diff --git a/rust/spheres.png b/rust/spheres.png new file mode 100644 index 0000000..6dbd75b Binary files /dev/null and b/rust/spheres.png differ diff --git a/rust/src/bvh.rs b/rust/src/bvh.rs new file mode 100644 index 0000000..ec1d6bc --- /dev/null +++ b/rust/src/bvh.rs @@ -0,0 +1,185 @@ +//! Bounding volume hierarchy over the scene's bounded objects, so a ray only +//! tests the objects whose boxes it passes through instead of all of them. + +use crate::scene::Ray; +use crate::vec3::Vec3; + +#[derive(Clone, Copy, Debug)] +pub struct Aabb { + pub min: Vec3, + pub max: Vec3, +} + +impl Aabb { + pub const EMPTY: Aabb = Aabb { + min: Vec3::new(f64::INFINITY, f64::INFINITY, f64::INFINITY), + max: Vec3::new(f64::NEG_INFINITY, f64::NEG_INFINITY, f64::NEG_INFINITY), + }; + + pub fn union(self, o: Aabb) -> Aabb { + Aabb { + min: Vec3::new(self.min.x.min(o.min.x), self.min.y.min(o.min.y), self.min.z.min(o.min.z)), + max: Vec3::new(self.max.x.max(o.max.x), self.max.y.max(o.max.y), self.max.z.max(o.max.z)), + } + } + + pub fn centroid(&self) -> Vec3 { + (self.min + self.max) * 0.5 + } + + /// Slab test. `inv_dir` components may be infinite for axis-parallel rays; + /// the resulting NaNs are discarded by f64::min/max. + fn hit(&self, origin: Vec3, inv_dir: Vec3, mut t_min: f64, mut t_max: f64) -> bool { + for axis in 0..3 { + let mut t0 = (self.min[axis] - origin[axis]) * inv_dir[axis]; + let mut t1 = (self.max[axis] - origin[axis]) * inv_dir[axis]; + if inv_dir[axis] < 0.0 { + std::mem::swap(&mut t0, &mut t1); + } + t_min = t_min.max(t0); + t_max = t_max.min(t1); + if t_max < t_min { + return false; + } + } + true + } +} + +const LEAF_SIZE: usize = 4; + +/// Flattened node: a leaf covers `prims[start..start + count]`; an interior +/// node (count == 0) has its left child at the next index and its right child at `right`. +struct Node { + bounds: Aabb, + start: u32, + count: u32, + right: u32, + axis: u8, +} + +pub struct Bvh { + nodes: Vec, + /// Object indices, reordered so each leaf's objects are contiguous. + prims: Vec, +} + +impl Bvh { + pub fn build(bounds: &[Aabb]) -> Bvh { + let mut bvh = Bvh { + nodes: Vec::with_capacity(2 * bounds.len()), + prims: (0..bounds.len() as u32).collect(), + }; + if !bounds.is_empty() { + let centroids: Vec = bounds.iter().map(Aabb::centroid).collect(); + bvh.build_range(0, bounds.len(), bounds, ¢roids); + } + bvh + } + + /// Median split along the axis where the object centroids spread widest. + fn build_range(&mut self, start: usize, end: usize, bounds: &[Aabb], centroids: &[Vec3]) -> usize { + let prims = &self.prims[start..end]; + let node_bounds = prims.iter().fold(Aabb::EMPTY, |b, &p| b.union(bounds[p as usize])); + let centroid_bounds = prims.iter().fold(Aabb::EMPTY, |b, &p| { + let c = centroids[p as usize]; + b.union(Aabb { min: c, max: c }) + }); + let extent = centroid_bounds.max - centroid_bounds.min; + let axis = if extent.x >= extent.y && extent.x >= extent.z { + 0 + } else if extent.y >= extent.z { + 1 + } else { + 2 + }; + + let idx = self.nodes.len(); + self.nodes.push(Node { + bounds: node_bounds, + start: start as u32, + count: (end - start) as u32, + right: 0, + axis: axis as u8, + }); + if end - start <= LEAF_SIZE || extent[axis] <= 0.0 { + return idx; + } + + let half = (end - start) / 2; + self.prims[start..end].select_nth_unstable_by(half, |&a, &b| { + centroids[a as usize][axis].total_cmp(¢roids[b as usize][axis]) + }); + self.build_range(start, start + half, bounds, centroids); + let right = self.build_range(start + half, end, bounds, centroids); + self.nodes[idx].count = 0; + self.nodes[idx].right = right as u32; + idx + } + + /// Calls `visit` for each object in a leaf whose box the ray passes through + /// within (t_min, t_max). `visit` returns Some(t) to shrink t_max after a + /// hit, or Err(()) to stop traversal entirely. + fn traverse( + &self, + r: &Ray, + t_min: f64, + mut t_max: f64, + mut visit: impl FnMut(usize, f64) -> Result, ()>, + ) { + if self.nodes.is_empty() { + return; + } + let inv_dir = Vec3::new(1.0 / r.direction.x, 1.0 / r.direction.y, 1.0 / r.direction.z); + let mut stack = [0u32; 64]; + let mut sp = 1; + while sp > 0 { + sp -= 1; + let node = &self.nodes[stack[sp] as usize]; + if !node.bounds.hit(r.origin, inv_dir, t_min, t_max) { + continue; + } + if node.count > 0 { + for &p in &self.prims[node.start as usize..(node.start + node.count) as usize] { + match visit(p as usize, t_max) { + Ok(Some(t)) => t_max = t, + Ok(None) => {} + Err(()) => return, + } + } + } else { + // Push the far child first so the near one is visited first, + // which shrinks t_max sooner. + let left = stack[sp] + 1; + let (near, far) = if r.direction[node.axis as usize] < 0.0 { + (node.right, left) + } else { + (left, node.right) + }; + stack[sp] = far; + stack[sp + 1] = near; + sp += 2; + } + } + } + + /// Closest hit: `hit(object, t_max)` returns the hit distance if the object + /// is hit closer than t_max. + pub fn closest(&self, r: &Ray, t_min: f64, t_max: f64, mut hit: impl FnMut(usize, f64) -> Option) { + self.traverse(r, t_min, t_max, |p, t_max| Ok(hit(p, t_max))); + } + + /// True as soon as `hit(object)` reports any hit. + pub fn any(&self, r: &Ray, t_min: f64, t_max: f64, mut hit: impl FnMut(usize) -> bool) -> bool { + let mut found = false; + self.traverse(r, t_min, t_max, |p, _| { + if hit(p) { + found = true; + Err(()) + } else { + Ok(None) + } + }); + found + } +} diff --git a/rust/src/camera.rs b/rust/src/camera.rs index 7b00d19..7f6c6b0 100644 --- a/rust/src/camera.rs +++ b/rust/src/camera.rs @@ -15,8 +15,8 @@ impl Camera { let w = (lookfrom - lookat).normalize(); let u = vup.cross(w).normalize(); let v = w.cross(u); - let horizontal = focus_dist * half_width * u; - let vertical = focus_dist * half_height * v; + let horizontal = 2.0 * focus_dist * half_width * u; + let vertical = 2.0 * focus_dist * half_height * v; let lower_left_corner = lookfrom - horizontal / 2.0 - vertical / 2.0 - focus_dist * w; Camera { origin: lookfrom, diff --git a/rust/src/main.rs b/rust/src/main.rs index 690b9b8..bb70b6f 100644 --- a/rust/src/main.rs +++ b/rust/src/main.rs @@ -1,5 +1,7 @@ +mod bvh; mod camera; mod scene; +mod scenes; mod vec3; use std::f64::consts::PI; @@ -11,43 +13,20 @@ use rand::{Rng, SeedableRng}; use rayon::prelude::*; use camera::Camera; -use scene::{cornell_box, BoxLight, HitRecord, Material, Ray, Scene, EPSILON}; +use scene::{HitRecord, Light, Material, Ray, Scene, EPSILON}; use vec3::Vec3; -const SHADOW_EPSILON: f64 = 1e-3; +/// Shadow rays stop this fraction short of the light, so they don't hit the +/// light's own surface. Relative, because scenes range from ~1 to ~1000 units. +const SHADOW_RELATIVE_EPSILON: f64 = 1e-4; -/// Direct light from one uniformly sampled point on the box light's surface. -fn sample_box_light(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::() * 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::() < 0.5; - let mut p = [ - lo.x + size.x * rng.gen::(), - lo.y + size.y * rng.gen::(), - lo.z + size.z * rng.gen::(), - ]; - p[axis] = if positive_side { hi[axis] } else { lo[axis] }; - let mut n = [0.0; 3]; - n[axis] = if positive_side { 1.0 } else { -1.0 }; - let (lp, light_normal) = (Vec3::new(p[0], p[1], p[2]), Vec3::new(n[0], n[1], n[2])); +/// Direct light from one uniformly sampled point on `light`'s surface. +fn sample_direct(scene: &Scene, light: &Light, rec: &HitRecord, albedo: Vec3, rng: &mut R) -> Vec3 { + let (lp, light_normal) = light.sample(rng); let to_light = lp - rec.point; let dist_sq = to_light.dot(to_light); - if dist_sq < 1e-8 { + if dist_sq < 1e-12 { return Vec3::ZERO; } let dist = dist_sq.sqrt(); @@ -58,13 +37,14 @@ fn sample_box_light(scene: &Scene, light: &BoxLight, rec: &HitRecord, al return Vec3::ZERO; } - if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist - SHADOW_EPSILON) { + if scene.blocked(rec.point + EPSILON * rec.normal, d_to_light, dist * (1.0 - SHADOW_RELATIVE_EPSILON)) { return Vec3::ZERO; } + // Uniform area sampling has pdf 1 / area. let brdf = albedo / PI; let geometry = cos_at_rec * cos_at_light / dist_sq; - light.emission.mul_elem(brdf) * (geometry * total_area) + light.emission.mul_elem(brdf) * (geometry * light.area()) } /// Bounces before Russian roulette may terminate a path. @@ -82,14 +62,13 @@ fn trace(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> V for depth in 0..max_depth { let Some(rec) = scene.hit(&ray, EPSILON, f64::INFINITY) else { - let t = 0.5 * (ray.direction.y + 1.0); - let sky = Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t; - radiance += throughput.mul_elem(sky); + radiance += throughput.mul_elem(scene.background.radiance(ray.direction)); break; }; if let Material::Light { emission } = rec.material { - if count_emission { + // Lights emit from their front side only (matters for one-sided rects). + if count_emission && rec.front_face { radiance += throughput.mul_elem(emission); } break; @@ -100,8 +79,9 @@ fn trace(scene: &Scene, primary: &Ray, max_depth: u32, rng: &mut R) -> V }; if let Material::Diffuse { albedo } = rec.material { - for light in &scene.lights { - radiance += throughput.mul_elem(sample_box_light(scene, light, &rec, albedo, rng)); + if let Some((light, pick_prob)) = scene.pick_light(rec.point, rng) { + let direct = sample_direct(scene, light, &rec, albedo, rng) / pick_prob; + radiance += throughput.mul_elem(direct); } count_emission = false; } else { @@ -142,10 +122,6 @@ const ADAPTIVE_MIN_SPP: u32 = 32; /// ...then later passes grow each still-noisy pixel's count by this fraction. const ADAPTIVE_GROWTH: f64 = 0.5; -fn luminance(c: Vec3) -> f64 { - 0.2126 * c.x + 0.7152 * c.y + 0.0722 * c.z -} - /// Running mean and variance of a pixel's luminance (Welford's algorithm). #[derive(Default)] struct PixelStats { @@ -194,7 +170,7 @@ impl Renderer<'_> { let v = ((height - 1 - j) as f64 + px.rng.gen::()) / height as f64; let sample = trace(self.scene, &self.cam.get_ray(u, v), self.s.max_depth, &mut px.rng); px.sum += sample; - px.stats.add(luminance(sample)); + px.stats.add(sample.luminance()); } } @@ -234,15 +210,16 @@ fn bar(len: u64, template: &str) -> ProgressBar { /// Returns the gamma-encoded image and the total number of samples taken. fn render(scene: &Scene, s: &Settings) -> (Vec, u64) { let (width, height) = (s.width, s.height); + let view = scene.view; let r = Renderer { scene, cam: Camera::new( - Vec3::new(278.0, 278.0, -800.0), - Vec3::new(278.0, 278.0, 278.0), + view.lookfrom, + view.lookat, Vec3::new(0.0, 1.0, 0.0), - 40.0, + view.vfov, width as f64 / height as f64, - 1078.0, + (view.lookfrom - view.lookat).length(), ), s, }; @@ -356,32 +333,73 @@ fn print_region_table(image: &[Vec3], width: usize, height: usize) { } } +const USAGE: &str = "usage: pathtracer [--scene NAME] [--list-scenes] \ +[width] [height] [samples_per_pixel] [max_bounces] [output] [noise_threshold]"; + +fn list_scenes() { + println!("Scenes:"); + for s in scenes::SCENES { + println!(" {:<10} {}", s.name, s.description); + } +} + fn main() { - let args: Vec = std::env::args().collect(); + // Pull out the --scene / --list-scenes options; the rest are positional. + let mut scene_name = String::from("cornell"); + let mut args: Vec = Vec::new(); + let mut raw = std::env::args().skip(1); + while let Some(a) = raw.next() { + if a == "--list-scenes" { + list_scenes(); + return; + } else if a == "--scene" { + scene_name = raw.next().unwrap_or_else(|| { + eprintln!("{USAGE}"); + std::process::exit(2); + }); + } else if let Some(name) = a.strip_prefix("--scene=") { + scene_name = name.to_string(); + } else if a == "--help" || a == "-h" { + println!("{USAGE}"); + list_scenes(); + return; + } else { + args.push(a); + } + } + let Some(scene_info) = scenes::find(&scene_name) else { + eprintln!("unknown scene '{scene_name}'"); + list_scenes(); + std::process::exit(2); + }; + let arg = |i: usize, default: &str| args.get(i).cloned().unwrap_or_else(|| default.to_string()); fn parse(s: String) -> T { s.parse().unwrap_or_else(|_| { - eprintln!("usage: pathtracer [width] [height] [samples_per_pixel] [max_bounces] [output] [noise_threshold]"); + eprintln!("{USAGE}"); std::process::exit(2); }) } let settings = Settings { - width: parse(arg(1, "400")), - height: parse(arg(2, "300")), - spp: parse(arg(3, "50")), - max_depth: parse(arg(4, "15")), - output: arg(5, "cornell_box.png"), - noise_threshold: parse(arg(6, "0")), + width: parse(arg(0, "400")), + height: parse(arg(1, "300")), + spp: parse(arg(2, "50")), + max_depth: parse(arg(3, "15")), + output: arg(4, &format!("{scene_name}.png")), + noise_threshold: parse(arg(5, "0")), }; if settings.width == 0 || settings.height == 0 || settings.spp == 0 || settings.noise_threshold < 0.0 { eprintln!("width, height and samples_per_pixel must be >= 1, noise_threshold >= 0"); std::process::exit(2); } - println!("=== Cornell Box Path Tracer (Rust) ==="); - println!("Features: reflections, refraction/caustics, multiple bounces, NEE"); + let scene = (scene_info.build)(); + + println!("=== Path Tracer (Rust) ==="); + println!("Features: reflections, refraction/caustics, multiple bounces, NEE, BVH"); println!(); + println!("Scene: {} ({} objects, {} lights)", scene_name, scene.object_count(), scene.light_count()); println!("Resolution: {}x{}", settings.width, settings.height); println!("Samples per pixel: {} (max)", settings.spp); println!("Max bounces: {} (Russian roulette after {})", settings.max_depth, RR_MIN_DEPTH); @@ -394,7 +412,6 @@ fn main() { println!("Output: {}", settings.output); println!(); - let scene = cornell_box(); let start = Instant::now(); let (image, total_samples) = render(&scene, &settings); let elapsed = start.elapsed().as_secs_f64(); @@ -413,6 +430,9 @@ fn main() { total_samples as f64 / elapsed / 1e6, total_samples as f64 / pixels ); - println!(); - print_region_table(&image, settings.width, settings.height); + // The regions (left wall, floor, ...) only mean something for the Cornell box. + if scene_name == "cornell" { + println!(); + print_region_table(&image, settings.width, settings.height); + } } diff --git a/rust/src/scene.rs b/rust/src/scene.rs index b3e4be4..549476f 100644 --- a/rust/src/scene.rs +++ b/rust/src/scene.rs @@ -1,5 +1,6 @@ use rand::Rng; +use crate::bvh::{Aabb, Bvh}; use crate::vec3::{reflect, refract, schlick, Vec3}; pub const EPSILON: f64 = 1e-6; @@ -114,6 +115,10 @@ pub enum Shape { Sphere { center: Vec3, radius: f64 }, Plane { point: Vec3, normal: Vec3 }, Box { min: Vec3, max: Vec3 }, + /// Axis-aligned rectangle spanning min..max, flat along the axis of + /// `normal` (which must be ±x, ±y or ±z). As a light it emits only on the + /// `normal` side. + Rect { min: Vec3, max: Vec3, normal: Vec3 }, } impl Shape { @@ -183,6 +188,38 @@ impl Shape { n[hit_axis] = sign; Some((t_lo, Vec3::new(n[0], n[1], n[2]))) } + Shape::Rect { min, max, normal } => { + let denom = r.direction.dot(normal); + if denom.abs() < 1e-12 { + return None; + } + let t = (min - r.origin).dot(normal) / denom; + if t < t_min || t > t_max { + return None; + } + let p = r.at(t); + let tol = 1e-9; + let inside = (0..3).all(|a| p[a] >= min[a] - tol && p[a] <= max[a] + tol); + inside.then_some((t, normal)) + } + } + } +} + +impl Shape { + /// None for unbounded shapes (planes), which are tested outside the BVH. + fn bounds(&self) -> Option { + match *self { + Shape::Sphere { center, radius } => Some(Aabb { + min: center - Vec3::splat(radius), + max: center + Vec3::splat(radius), + }), + Shape::Box { min, max } => Some(Aabb { min, max }), + Shape::Rect { min, max, .. } => Some(Aabb { + min: min - Vec3::splat(1e-9), + max: max + Vec3::splat(1e-9), + }), + Shape::Plane { .. } => None, } } } @@ -193,79 +230,278 @@ pub struct Object { pub material: Material, } -/// An axis-aligned box light, sampled explicitly for next event estimation. +/// An emitter that next event estimation can sample: a box, sphere or rectangle. #[derive(Clone, Copy, Debug)] -pub struct BoxLight { - pub min: Vec3, - pub max: Vec3, +pub struct Light { + shape: Shape, pub emission: Vec3, + /// Total emitted power (up to a constant): luminance * area. + power: f64, + center: Vec3, + /// Half the bounding box diagonal: how close a point can be to the centre + /// while still outside the light. + extent: f64, +} + +impl Light { + fn new(shape: Shape, emission: Vec3) -> Self { + let bounds = shape.bounds().expect("planes cannot be lights"); + let mut light = Light { + shape, + emission, + power: 0.0, + center: bounds.centroid(), + extent: 0.5 * (bounds.max - bounds.min).length(), + }; + light.power = emission.luminance() * light.area(); + light + } + + pub fn area(&self) -> f64 { + match self.shape { + Shape::Box { min, max } => { + let d = max - min; + 2.0 * (d.x * d.y + d.y * d.z + d.x * d.z) + } + Shape::Sphere { radius, .. } => 4.0 * std::f64::consts::PI * radius * radius, + Shape::Rect { min, max, normal } => { + let d = max - min; + if normal.x != 0.0 { + d.y * d.z + } else if normal.y != 0.0 { + d.x * d.z + } else { + d.x * d.y + } + } + Shape::Plane { .. } => unreachable!("planes cannot be lights"), + } + } + + /// Rough estimate of how much this light contributes at `point`: its + /// power over the squared distance, clamped so points near or inside the + /// light's bounds don't blow up. + fn importance(&self, point: Vec3) -> f64 { + let d = point - self.center; + self.power / d.dot(d).max(self.extent * self.extent) + } + + /// A uniformly distributed point on the surface (pdf 1 / area) and the + /// outward normal there. + pub fn sample(&self, rng: &mut R) -> (Vec3, Vec3) { + match self.shape { + Shape::Box { min, max } => { + // Pick a face with probability proportional to its area, then + // a uniform point on it. + let size = max - min; + let face_areas = [size.y * size.z, size.x * size.z, size.x * size.y]; + let u = rng.gen::() * (face_areas[0] + face_areas[1] + face_areas[2]); + let axis = if u < face_areas[0] { + 0 + } else if u < face_areas[0] + face_areas[1] { + 1 + } else { + 2 + }; + let positive_side = rng.gen::() < 0.5; + let mut p = [ + min.x + size.x * rng.gen::(), + min.y + size.y * rng.gen::(), + min.z + size.z * rng.gen::(), + ]; + p[axis] = if positive_side { max[axis] } else { min[axis] }; + let mut n = [0.0; 3]; + n[axis] = if positive_side { 1.0 } else { -1.0 }; + (Vec3::new(p[0], p[1], p[2]), Vec3::new(n[0], n[1], n[2])) + } + Shape::Sphere { center, radius } => { + let z = 1.0 - 2.0 * rng.gen::(); + let r = (1.0 - z * z).max(0.0).sqrt(); + let phi = 2.0 * std::f64::consts::PI * rng.gen::(); + let n = Vec3::new(r * phi.cos(), r * phi.sin(), z); + (center + radius * n, n) + } + Shape::Rect { min, max, normal } => { + let size = max - min; + let p = Vec3::new( + min.x + size.x * rng.gen::(), + min.y + size.y * rng.gen::(), + min.z + size.z * rng.gen::(), + ); + (p, normal) + } + Shape::Plane { .. } => unreachable!("planes cannot be lights"), + } + } +} + +/// What a ray that leaves the scene sees. +#[derive(Clone, Copy, Debug)] +pub enum Background { + /// White-to-blue gradient, as in raytracer.py. + Sky, + Color(Vec3), +} + +impl Background { + pub fn radiance(&self, direction: Vec3) -> Vec3 { + match *self { + Background::Sky => { + let t = 0.5 * (direction.y + 1.0); + Vec3::ONE * (1.0 - t) + Vec3::new(0.5, 0.7, 1.0) * t + } + Background::Color(c) => c, + } + } +} + +#[derive(Clone, Copy, Debug)] +pub struct View { + pub lookfrom: Vec3, + pub lookat: Vec3, + pub vfov: f64, +} + +pub struct SceneBuilder { + objects: Vec, + planes: Vec, + lights: Vec, + background: Background, + view: View, +} + +impl SceneBuilder { + pub fn new(view: View, background: Background) -> Self { + SceneBuilder { + objects: Vec::new(), + planes: Vec::new(), + lights: Vec::new(), + background, + view, + } + } + + pub fn add(&mut self, shape: Shape, material: Material) { + let obj = Object { shape, material }; + match shape { + Shape::Plane { .. } => self.planes.push(obj), + _ => self.objects.push(obj), + } + } + + pub fn plane(&mut self, point: Vec3, normal: Vec3, material: Material) { + self.add(Shape::Plane { point, normal: normal.normalize() }, material); + } + + pub fn sphere(&mut self, center: Vec3, radius: f64, material: Material) { + self.add(Shape::Sphere { center, radius }, material); + } + + pub fn cuboid(&mut self, min: Vec3, max: Vec3, material: Material) { + self.add(Shape::Box { min, max }, material); + } + + /// Every emitter must be added through a *_light method: after a diffuse + /// bounce, trace() only counts light found by sampling these. + pub fn box_light(&mut self, min: Vec3, max: Vec3, emission: Vec3) { + self.add_light(Shape::Box { min, max }, emission); + } + + pub fn sphere_light(&mut self, center: Vec3, radius: f64, emission: Vec3) { + self.add_light(Shape::Sphere { center, radius }, emission); + } + + /// One-sided rectangular light emitting towards `normal` (±x, ±y or ±z). + pub fn rect_light(&mut self, min: Vec3, max: Vec3, normal: Vec3, emission: Vec3) { + self.add_light(Shape::Rect { min, max, normal }, emission); + } + + fn add_light(&mut self, shape: Shape, emission: Vec3) { + self.lights.push(Light::new(shape, emission)); + self.add(shape, Material::Light { emission }); + } + + pub fn build(self) -> Scene { + let bounds: Vec = self.objects.iter().map(|o| o.shape.bounds().unwrap()).collect(); + Scene { + bvh: Bvh::build(&bounds), + objects: self.objects, + planes: self.planes, + lights: self.lights, + background: self.background, + view: self.view, + } + } } -#[derive(Default)] pub struct Scene { - pub objects: Vec, - pub lights: Vec, + objects: Vec, + planes: Vec, + bvh: Bvh, + lights: Vec, + pub background: Background, + pub view: View, } impl Scene { - pub fn add(&mut self, shape: Shape, material: Material) { - self.objects.push(Object { shape, material }); + pub fn object_count(&self) -> usize { + self.objects.len() + self.planes.len() } - pub fn add_box_light(&mut self, min: Vec3, max: Vec3, emission: Vec3) { - self.lights.push(BoxLight { min, max, emission }); - self.add(Shape::Box { min, max }, Material::Light { emission }); + pub fn light_count(&self) -> usize { + self.lights.len() + } + + /// Picks one light for shading `point`, with probability proportional to + /// its estimated contribution there (power / distance^2), so a scene with + /// many lights costs one shadow ray per shading point and mostly samples + /// nearby lights. Returns the light and the probability of picking it. + /// This is O(lights) per call; thousands of lights would want a light tree. + pub fn pick_light(&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::() * total; + let mut acc = 0.0; + for light in &self.lights { + let w = light.importance(point); + acc += w; + if u < acc { + return Some((light, w / total)); + } + } + // Rounding left u just past the last bucket. + let last = self.lights.last().unwrap(); + Some((last, last.importance(point) / total)) + } + } } pub fn hit(&self, r: &Ray, t_min: f64, t_max: f64) -> Option { let mut closest: Option<(f64, Vec3, Material)> = None; let mut min_t = t_max; - for obj in &self.objects { + for obj in &self.planes { if let Some((t, outward)) = obj.shape.hit(r, t_min, min_t) { - if t < min_t { - min_t = t; - closest = Some((t, outward, obj.material)); - } + min_t = t; + closest = Some((t, outward, obj.material)); } } + self.bvh.closest(r, t_min, min_t, |i, t_max| { + let obj = &self.objects[i]; + let (t, outward) = obj.shape.hit(r, t_min, t_max)?; + closest = Some((t, outward, obj.material)); + Some(t) + }); closest.map(|(t, outward, material)| HitRecord::new(r, t, outward, material)) } pub fn blocked(&self, origin: Vec3, direction: Vec3, t_max: f64) -> bool { - let r = Ray::new(origin, direction); - self.objects - .iter() - .any(|obj| obj.shape.hit(&r, EPSILON, t_max).is_some()) + let r = Ray::from_unit(origin, direction); + self.planes.iter().any(|obj| obj.shape.hit(&r, EPSILON, t_max).is_some()) + || self + .bvh + .any(&r, EPSILON, t_max, |i| self.objects[i].shape.hit(&r, EPSILON, t_max).is_some()) } } - -pub fn cornell_box() -> Scene { - let mut scene = Scene::default(); - - let white = Material::Diffuse { albedo: Vec3::new(0.73, 0.73, 0.73) }; - let gray = Material::Diffuse { albedo: Vec3::new(0.5, 0.5, 0.5) }; - let red = Material::Diffuse { albedo: Vec3::new(0.65, 0.05, 0.05) }; - let green = Material::Diffuse { albedo: Vec3::new(0.12, 0.45, 0.15) }; - let mirror = Material::Mirror { albedo: Vec3::new(0.95, 0.95, 0.95) }; - let glass = Material::Dielectric { ref_idx: 1.5, albedo: Vec3::ONE }; - - let plane = |p: Vec3, n: Vec3| Shape::Plane { point: p, normal: n.normalize() }; - scene.add(plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0)), gray); - scene.add(plane(Vec3::new(0.0, 500.0, 0.0), Vec3::new(0.0, -1.0, 0.0)), white); - scene.add(plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0)), white); - scene.add(plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0)), red); - scene.add(plane(Vec3::new(500.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0)), green); - - scene.add_box_light( - Vec3::new(0.0, 180.0, 180.0), - Vec3::new(25.0, 220.0, 320.0), - Vec3::splat(50.0), - ); - - let sphere = |c: Vec3, r: f64| Shape::Sphere { center: c, radius: r }; - scene.add(sphere(Vec3::new(250.0, 180.0, 250.0), 80.0), glass); - scene.add(sphere(Vec3::new(150.0, 180.0, 380.0), 80.0), mirror); - scene.add(sphere(Vec3::new(350.0, 160.0, 300.0), 60.0), red); - - scene -} diff --git a/rust/src/scenes.rs b/rust/src/scenes.rs new file mode 100644 index 0000000..fc79d8d --- /dev/null +++ b/rust/src/scenes.rs @@ -0,0 +1,220 @@ +//! The built-in scenes, selectable with `--scene NAME`. + +use rand::rngs::SmallRng; +use rand::{Rng, SeedableRng}; + +use crate::scene::{Background, Material, Scene, SceneBuilder, View}; +use crate::vec3::Vec3; + +pub struct SceneInfo { + pub name: &'static str, + pub description: &'static str, + pub build: fn() -> Scene, +} + +pub const SCENES: &[SceneInfo] = &[ + SceneInfo { + name: "cornell", + description: "Cornell box from raytracer.py: glass, mirror and diffuse spheres, wall light", + build: cornell, + }, + SceneInfo { + name: "caustics", + description: "Clear, amber and blue glass plus a mirror under one light, focusing caustics onto the floor", + build: caustics, + }, + SceneInfo { + name: "spheres", + description: "Large scene: ~480 random spheres under a sky (exercises the BVH)", + build: spheres, + }, + SceneInfo { + name: "night", + description: "The same sphere field at night, lit by dozens of glowing spheres (many-light sampling)", + build: night, + }, + SceneInfo { + name: "classic", + description: "The classic Cornell box: ceiling light, two boxes, colour bleeding", + build: classic, + }, +]; + +pub fn find(name: &str) -> Option<&'static SceneInfo> { + SCENES.iter().find(|s| s.name == name) +} + +fn diffuse(r: f64, g: f64, b: f64) -> Material { + Material::Diffuse { albedo: Vec3::new(r, g, b) } +} + +fn mirror(r: f64, g: f64, b: f64) -> Material { + Material::Mirror { albedo: Vec3::new(r, g, b) } +} + +/// Glass; `tint` is applied at every reflection/refraction, so colour +/// deepens with each pass through the surface. +fn glass(tint: Vec3) -> Material { + Material::Dielectric { ref_idx: 1.5, albedo: tint } +} + +fn cornell() -> Scene { + let mut b = SceneBuilder::new( + View { + lookfrom: Vec3::new(278.0, 278.0, -800.0), + lookat: Vec3::new(278.0, 278.0, 278.0), + vfov: 40.0, + }, + Background::Sky, + ); + + let white = diffuse(0.73, 0.73, 0.73); + b.plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(0.0, 1.0, 0.0), diffuse(0.5, 0.5, 0.5)); + b.plane(Vec3::new(0.0, 500.0, 0.0), Vec3::new(0.0, -1.0, 0.0), white); + b.plane(Vec3::new(0.0, 0.0, 500.0), Vec3::new(0.0, 0.0, -1.0), white); + b.plane(Vec3::new(0.0, 0.0, 0.0), Vec3::new(1.0, 0.0, 0.0), diffuse(0.65, 0.05, 0.05)); + b.plane(Vec3::new(500.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), diffuse(0.12, 0.45, 0.15)); + + b.box_light(Vec3::new(0.0, 180.0, 180.0), Vec3::new(25.0, 220.0, 320.0), Vec3::splat(50.0)); + + b.sphere(Vec3::new(250.0, 180.0, 250.0), 80.0, glass(Vec3::ONE)); + b.sphere(Vec3::new(150.0, 180.0, 380.0), 80.0, mirror(0.95, 0.95, 0.95)); + b.sphere(Vec3::new(350.0, 160.0, 300.0), 60.0, diffuse(0.65, 0.05, 0.05)); + b.build() +} + +fn caustics() -> Scene { + let mut b = SceneBuilder::new( + View { + lookfrom: Vec3::new(0.0, 160.0, -280.0), + lookat: Vec3::new(20.0, 30.0, 120.0), + vfov: 38.0, + }, + Background::Color(Vec3::ZERO), + ); + + b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.8, 0.8, 0.8)); + b.plane(Vec3::new(0.0, 0.0, 520.0), Vec3::new(0.0, 0.0, -1.0), diffuse(0.6, 0.6, 0.6)); + + // One downward-facing panel above and well behind the spheres. A glass + // sphere focuses light about 1.5 radii from its centre; with the light + // ~48 degrees from vertical that focus lands on the floor just in front + // of each sphere, facing the camera. A smaller panel gives sharper + // caustics but makes them slower to converge. + b.rect_light( + Vec3::new(-30.0, 250.0, 330.0), + Vec3::new(30.0, 250.0, 390.0), + Vec3::new(0.0, -1.0, 0.0), + Vec3::splat(90.0), + ); + + b.sphere(Vec3::new(-75.0, 45.0, 110.0), 45.0, glass(Vec3::ONE)); + b.sphere(Vec3::new(40.0, 38.0, 80.0), 38.0, glass(Vec3::new(1.0, 0.7, 0.3))); + b.sphere(Vec3::new(125.0, 30.0, 150.0), 30.0, glass(Vec3::new(0.4, 0.65, 1.0))); + b.sphere(Vec3::new(-10.0, 35.0, 210.0), 35.0, mirror(0.9, 0.9, 0.9)); + b.build() +} + +/// A grid of small random spheres around three large feature spheres (after +/// "Ray Tracing in One Weekend"). A `light_fraction` of the small spheres glow. +fn sphere_field(b: &mut SceneBuilder, light_fraction: f64) { + let mut rng = SmallRng::seed_from_u64(7); + b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), diffuse(0.5, 0.5, 0.5)); + + let big = [ + (Vec3::new(0.0, 1.0, 0.0), glass(Vec3::ONE)), + (Vec3::new(-4.0, 1.0, 0.0), diffuse(0.4, 0.2, 0.1)), + (Vec3::new(4.0, 1.0, 0.0), mirror(0.7, 0.6, 0.5)), + ]; + for (center, material) in big { + b.sphere(center, 1.0, material); + } + + for a in -11..11 { + for c in -11..11 { + let center = Vec3::new( + a as f64 + 0.9 * rng.gen::(), + 0.2, + c as f64 + 0.9 * rng.gen::(), + ); + if big.iter().any(|(bc, _)| (center - Vec3::new(bc.x, 0.2, bc.z)).length() < 1.2) { + continue; + } + let rnd_color = |rng: &mut SmallRng| Vec3::new(rng.gen(), rng.gen(), rng.gen()); + if rng.gen::() < light_fraction { + // Warm-to-white glow of varying strength. + let warmth: f64 = rng.gen(); + let color = Vec3::new(1.0, 0.55 + 0.4 * warmth, 0.2 + 0.6 * warmth); + b.sphere_light(center, 0.2, color * (6.0 + 10.0 * rng.gen::())); + continue; + } + let choice: f64 = rng.gen(); + let material = if choice < 0.75 { + let albedo = rnd_color(&mut rng).mul_elem(rnd_color(&mut rng)); + Material::Diffuse { albedo } + } else if choice < 0.9 { + Material::Mirror { albedo: Vec3::splat(0.5) + 0.5 * rnd_color(&mut rng) } + } else { + glass(Vec3::ONE) + }; + b.sphere(center, 0.2, material); + } + } +} + +fn spheres() -> Scene { + let mut b = SceneBuilder::new( + View { + lookfrom: Vec3::new(13.0, 2.0, 3.0), + lookat: Vec3::ZERO, + vfov: 20.0, + }, + Background::Sky, + ); + sphere_field(&mut b, 0.0); + b.build() +} + +fn night() -> Scene { + let mut b = SceneBuilder::new( + View { + lookfrom: Vec3::new(13.0, 2.0, 3.0), + lookat: Vec3::ZERO, + vfov: 20.0, + }, + Background::Color(Vec3::new(0.004, 0.006, 0.015)), + ); + sphere_field(&mut b, 0.12); + b.build() +} + +fn classic() -> Scene { + let mut b = SceneBuilder::new( + View { + lookfrom: Vec3::new(278.0, 278.0, -800.0), + lookat: Vec3::new(278.0, 278.0, 0.0), + vfov: 40.0, + }, + Background::Color(Vec3::ZERO), + ); + + let white = diffuse(0.73, 0.73, 0.73); + b.plane(Vec3::ZERO, Vec3::new(0.0, 1.0, 0.0), white); + b.plane(Vec3::new(0.0, 555.0, 0.0), Vec3::new(0.0, -1.0, 0.0), white); + b.plane(Vec3::new(0.0, 0.0, 555.0), Vec3::new(0.0, 0.0, -1.0), white); + // The camera looks down +z with +x to its left, so red ends up on the left as usual. + b.plane(Vec3::new(555.0, 0.0, 0.0), Vec3::new(-1.0, 0.0, 0.0), diffuse(0.65, 0.05, 0.05)); + b.plane(Vec3::ZERO, Vec3::new(1.0, 0.0, 0.0), diffuse(0.12, 0.45, 0.15)); + + // A flat panel just under the ceiling, shining down only. + b.rect_light( + Vec3::new(213.0, 554.0, 227.0), + Vec3::new(343.0, 554.0, 332.0), + Vec3::new(0.0, -1.0, 0.0), + Vec3::splat(15.0), + ); + + b.cuboid(Vec3::new(130.0, 0.0, 65.0), Vec3::new(295.0, 165.0, 230.0), white); + b.cuboid(Vec3::new(265.0, 0.0, 295.0), Vec3::new(430.0, 330.0, 460.0), white); + b.build() +} diff --git a/rust/src/vec3.rs b/rust/src/vec3.rs index 7396f94..f4894d6 100644 --- a/rust/src/vec3.rs +++ b/rust/src/vec3.rs @@ -56,6 +56,11 @@ impl Vec3 { pub fn mean(self) -> f64 { (self.x + self.y + self.z) / 3.0 } + + /// Rec. 709 luminance of a linear RGB colour. + pub fn luminance(self) -> f64 { + 0.2126 * self.x + 0.7152 * self.y + 0.0722 * self.z + } } impl Index for Vec3 {