diff --git a/shared/src/bxdfs/layered.rs b/shared/src/bxdfs/layered.rs index ecbd9d4..fa9423e 100644 --- a/shared/src/bxdfs/layered.rs +++ b/shared/src/bxdfs/layered.rs @@ -140,7 +140,10 @@ where Self { top, bottom, - thickness: thickness.max(Float::MIN), + // pbrt: `std::max(thickness, std::numeric_limits::min())` -- clamp to the + // smallest positive normal so the `dz / thickness` divisions stay finite. + // `Float::MIN` is the most negative finite value, so it never clamped. + thickness: thickness.max(Float::MIN_POSITIVE), g, albedo, max_depth, @@ -150,10 +153,17 @@ where } fn tr(&self, dz: Float, w: Vector3f) -> Float { - if dz.abs() <= Float::MIN { + // pbrt: `if (std::abs(dz) <= std::numeric_limits::min()) return 1;` + // C++ `numeric_limits::min()` is the smallest positive NORMAL value, which + // is `f32::MIN_POSITIVE` -- `Float::MIN` is the most negative finite value, so the + // guard could never fire. + if dz.abs() <= Float::MIN_POSITIVE { return 1.; } - -(dz / w.z()).abs().exp() + // pbrt: `FastExp(-std::abs(dz / w.z))`. The minus sign belongs on the EXPONENT; + // `-(x).abs().exp()` negates the result and leaves a growing `exp(+|x|)`, which + // made transmittance negative and unbounded. + fast_exp(-(dz / w.z()).abs()) } #[allow(clippy::too_many_arguments)] diff --git a/shared/src/spectra/simple.rs b/shared/src/spectra/simple.rs index 16ff6b7..9ac383e 100644 --- a/shared/src/spectra/simple.rs +++ b/shared/src/spectra/simple.rs @@ -237,7 +237,10 @@ impl PiecewiseLinearSpectrum { } } - pub fn from_interleaved(data: &[Float], _normalize: bool) -> Self { + /// pbrt `PiecewiseLinearSpectrum::FromInterleaved` (`util/spectrum.cpp`): `(lambda, value)` + /// pairs, extended flat to cover the full visible range, and -- when `normalize` is set -- + /// scaled so that `InnerProduct(spec, Y) == CIE_Y_integral` ("normalize to luminance 1"). + pub fn from_interleaved(data: &[Float], normalize: bool) -> Self { assert!( data.len() % 2 == 0, "Interleaved data must have even length" @@ -250,13 +253,45 @@ impl PiecewiseLinearSpectrum { } pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(core::cmp::Ordering::Equal)); - let mut lambdas = gvec_with_capacity(n); - let mut values = gvec_with_capacity(n); + let mut lambdas: GVec = gvec_with_capacity(n + 2); + let mut values: GVec = gvec_with_capacity(n + 2); + + // Extend samples to cover the range of visible wavelengths if needed. + if pairs[0].0 > LAMBDA_MIN as Float { + lambdas.push(LAMBDA_MIN as Float - 1.0); + values.push(pairs[0].1); + } for (l, v) in pairs.iter() { lambdas.push(*l); values.push(*v); } - Self::new(lambdas, values) + if *lambdas.last().unwrap() < LAMBDA_MAX as Float { + lambdas.push(LAMBDA_MAX as Float + 1.0); + values.push(*values.last().unwrap()); + } + + let mut spec = Self::new(lambdas, values); + if normalize { + // Normalize to have luminance of 1. + spec.scale(CIE_Y_INTEGRAL / spec.inner_product_with_cie_y()); + } + spec + } + + /// `InnerProduct(self, Spectra::Y())` -- pbrt sums over integer wavelengths across the + /// visible range, which is exactly the sampling of the tabulated `CIE_Y` curve. + pub fn inner_product_with_cie_y(&self) -> Float { + let mut integral = 0.0; + for (i, y) in CIE_Y.iter().enumerate() { + integral += *y * self.evaluate(LAMBDA_MIN as Float + i as Float); + } + integral + } + + pub fn scale(&mut self, s: Float) { + for v in self.values.iter_mut() { + *v *= s; + } } } diff --git a/src/core/film.rs b/src/core/film.rs index abac594..ae6d99c 100644 --- a/src/core/film.rs +++ b/src/core/film.rs @@ -282,16 +282,17 @@ impl CreateFilmBase for FilmBase { (full_resolution.y() as Float * crop.p_max.y()).ceil() as i32, ); - let mut pixel_bounds = Bounds2i::from_points(p_min, p_max); + let pixel_bounds = Bounds2i::from_points(p_min, p_max); if pixel_bounds.is_empty() { eprintln!("{}: Film crop window results in empty pixel bounds.", loc); } - let rad = filter.radius(); - let expansion = Point2i::new(rad.x().ceil() as i32, rad.y().ceil() as i32); - pixel_bounds = pixel_bounds.expand(expansion); - + // NOTE: pbrt does NOT expand pixelBounds by the filter radius (film.cpp:97, + // `pixelBounds = Bounds2i(Point2i(0, 0), fullResolution)`, then only intersected + // with "pixelbounds"/"cropwindow"). The filter radius widens SampleBounds(), + // never the film's stored pixel array. Expanding here made the film 1372x1030 + // instead of 1368x1026 and put the pixel origin at (-2,-2). let diagonal_mm = params.get_one_float("diagonal", 35.0)?; // let filename = params.get_one_string("filename", "pbrt.exr"); diff --git a/src/spectra/data.rs b/src/spectra/data.rs index 9f16eff..e2809fd 100644 --- a/src/spectra/data.rs +++ b/src/spectra/data.rs @@ -6,9 +6,11 @@ use std::collections::HashMap; use std::sync::LazyLock; pub fn create_cie(data: &[Float]) -> DenselySampledSpectrum { + // The CIE X/Y/Z curves are tabulated at 1nm over [360, 830]. (A 95-entry arm used to map + // the hand-normalized CIE_D65 table onto 300nm/5nm, but that table actually starts at + // 360nm; D65 now goes through PiecewiseLinearSpectrum::from_interleaved like pbrt.) let (start_lambda, step) = match data.len() { 471 => (360.0, 1.0), - 95 => (300.0, 5.0), n => panic!("Unexpected CIE data length: {}", n), }; let lambdas: Vec = (0..data.len()) diff --git a/src/spectra/mod.rs b/src/spectra/mod.rs index 00183d3..afb4de4 100644 --- a/src/spectra/mod.rs +++ b/src/spectra/mod.rs @@ -3,8 +3,10 @@ use crate::spectra::colorspace::CreateRGBColorSpace; use anyhow::{anyhow, Result}; use shared::core::geometry::Point2f; use shared::core::spectrum::{Spectrum, StandardSpectra}; -use shared::spectra::cie::{CIE_D65, CIE_X, CIE_Y, CIE_Z}; -use shared::spectra::{DenselySampledSpectrum, DeviceStandardColorSpaces, RGBColorSpace}; +use shared::spectra::cie::{CIE_ILLUM_D6500, CIE_X, CIE_Y, CIE_Z}; +use shared::spectra::{ + DenselySampledSpectrum, DeviceStandardColorSpaces, PiecewiseLinearSpectrum, RGBColorSpace, +}; use shared::Ptr; use std::sync::{Arc, LazyLock, OnceLock}; @@ -18,8 +20,14 @@ pub static CIE_Y_DATA: LazyLock = LazyLock::new(|| data::create_cie(&CIE_Y)); pub static CIE_Z_DATA: LazyLock = LazyLock::new(|| data::create_cie(&CIE_Z)); -pub static CIE_D65_DATA: LazyLock = - LazyLock::new(|| data::create_cie(&CIE_D65)); +/// pbrt builds D65 as `GetNamedSpectrum("stdillum-D65")`, i.e. +/// `PiecewiseLinearSpectrum::FromInterleaved(CIE_Illum_D6500, /*normalize=*/true)`, which +/// scales it so `InnerProduct(spec, Y) == CIE_Y_integral`. Do the same rather than carrying a +/// pre-normalized copy of the table. +pub static CIE_D65_DATA: LazyLock = LazyLock::new(|| { + let pls = PiecewiseLinearSpectrum::from_interleaved(&CIE_ILLUM_D6500, true); + DenselySampledSpectrum::from_spectrum(&Spectrum::Piecewise(shared::leak(pls))) +}); pub fn cie_x() -> Spectrum { Spectrum::Dense(Ptr::from(&*CIE_X_DATA)) diff --git a/src/wavefront/aggregate.rs b/src/wavefront/aggregate.rs index c1c80b9..31366d7 100644 --- a/src/wavefront/aggregate.rs +++ b/src/wavefront/aggregate.rs @@ -142,6 +142,9 @@ impl WavefrontAggregate for CpuAggregate { dndvs: intr.shading.dndv, }; if let Some(slot) = eval_q.push(item) { + // The queue is reset every depth/batch/sample, so `slot < 10` fires on + // every pass -- this print dominated render time. Gated behind cpu_debug. + #[cfg(feature = "cpu_debug")] if slot < 10 { eprintln!( "ENQUEUE[{slot}] pixel={:?} depth={} \ @@ -156,6 +159,7 @@ impl WavefrontAggregate for CpuAggregate { item.uv, item.material, item.area_light, item.face_index, ); } + let _ = slot; } }); } diff --git a/src/wavefront/integrator.rs b/src/wavefront/integrator.rs index 92b1228..099c40e 100644 --- a/src/wavefront/integrator.rs +++ b/src/wavefront/integrator.rs @@ -485,6 +485,9 @@ impl CpuWavefrontRenderer { (0..n as usize).into_par_iter().for_each(|i| { let w = unsafe { queue.storage.get(i) }; + // Fires on every material-queue pass (reset each depth/batch/sample), so it + // ran continuously and dominated render time. Gated behind cpu_debug. + #[cfg(feature = "cpu_debug")] if i < 10 { eprintln!( "DEQUEUE[{i}] pixel={:?} depth={} \ @@ -728,7 +731,6 @@ impl CpuWavefrontRenderer { if !pixel_bounds.contains_exclusive(p_pixel) { return; } - let l = self.pixel_sample_state.l.get(pixel_index); let camera_weight = self.pixel_sample_state.camera_ray_weight.get(pixel_index); let weighted_l = l * camera_weight;