Fix colorspace illuminant, LayeredBxDF transmittance, and film pixel bounds

- from_interleaved ignored its  arg, so D65 came from a separate
  hand-normalized table that was scaled to luminance 1 instead of
  CIE_Y_integral and mapped to 300nm when it starts at 360nm. sRGB
  rgb_from_xyz came out ~120x too large with the wrong white point, making
  every render ~120x too bright with a colour cast. Now built as pbrt does:
  from_interleaved(&CIE_ILLUM_D6500, true).

- LayeredBxDF::Tr:  puts the minus on the result,
  not the exponent, so transmittance was negative and grew exponentially.
  killeroo-coated-gold went from 15 Infs / min -4168 / stddev 638 to clean.

- layered.rs used Float::MIN where C++ has numeric_limits<Float>::min(). Those
  differ: the latter is f32::MIN_POSITIVE. Both guards were dead, leaving a
  zero thickness to divide by.

- FilmBase::create expanded pixel_bounds by the filter radius. pbrt never does;
  the radius widens SampleBounds() only. Output was 1372x1030, not 1368x1026.

- Gate the ENQUEUE/DEQUEUE wavefront prints behind cpu_debug; their queues
  reset every depth/batch/sample so they fired constantly (100h renders).

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<Float>::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<Float>::min()) return 1;`
+        // C++ `numeric_limits<Float>::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)]

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 lambdas: GVec<Float> = gvec_with_capacity(n + 2);
-        let mut values = gvec_with_capacity(n);
+        let mut values: GVec<Float> = 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;
+        }
     }
 }
 

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();
+        // NOTE: pbrt does NOT expand pixelBounds by the filter radius (film.cpp:97,
-        let expansion = Point2i::new(rad.x().ceil() as i32, rad.y().ceil() as i32);
+        // `pixelBounds = Bounds2i(Point2i(0, 0), fullResolution)`, then only intersected
-        pixel_bounds = pixel_bounds.expand(expansion);
+        // 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");
 

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<Float> = (0..data.len())

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::cie::{CIE_ILLUM_D6500, CIE_X, CIE_Y, CIE_Z};
-use shared::spectra::{DenselySampledSpectrum, DeviceStandardColorSpaces, RGBColorSpace};
+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<DenselySampledSpectrum> =
     LazyLock::new(|| data::create_cie(&CIE_Y));
 pub static CIE_Z_DATA: LazyLock<DenselySampledSpectrum> =
     LazyLock::new(|| data::create_cie(&CIE_Z));
-pub static CIE_D65_DATA: LazyLock<DenselySampledSpectrum> =
+/// pbrt builds D65 as `GetNamedSpectrum("stdillum-D65")`, i.e.
-    LazyLock::new(|| data::create_cie(&CIE_D65));
+/// `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<DenselySampledSpectrum> = 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))

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

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;