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/build.rs

@@ -2,4 +2,13 @@ fn main() {
     // This allows "spirv" to be used in #[cfg(target_arch = "...")]
     // without triggering a warning.
     println!("cargo:rustc-check-cfg=cfg(target_arch, values(\"spirv\"))");
+
+    // `gpu` is set for every device backend, so host-only code can be gated once
+    // as #[cfg(not(gpu))] instead of naming each target. Adding a backend means
+    // editing this line, not 30-odd cfg attributes.
+    println!("cargo::rustc-check-cfg=cfg(gpu)");
+    let target = std::env::var("TARGET").unwrap_or_default();
+    if target.contains("spirv") || target.contains("cuda") {
+        println!("cargo::rustc-cfg=gpu");
+    }
 }

shared/src/bxdfs/complex.rs

@@ -8,7 +8,7 @@ use crate::core::scattering::{
     TrowbridgeReitzDistribution, fr_complex_from_spectrum, fr_dielectric, fresnel_moment1, reflect,
     refract,
 };
-use crate::spectra::{DeviceStandardColorSpaces, RGBUnboundedSpectrum, SampledSpectrum};
+use crate::spectra::{N_SPECTRUM_SAMPLES, RGBColorSpace, RGBUnboundedSpectrum, SampledSpectrum};
 use crate::utils::gpu_array_from_fn;
 use crate::utils::math::{
     clamp, fast_exp, i0, lerp, log_i0, radians, safe_acos, safe_asin, safe_sqrt, sample_discrete,
@@ -17,7 +17,7 @@ use crate::utils::math::{
 use crate::utils::sampling::{
     cosine_hemisphere_pdf, sample_cosine_hemisphere, sample_trimmed_logistic,
 };
-use crate::{Float, INV_2_PI, INV_PI, PI};
+use crate::{Float, INV_2_PI, INV_PI, PI, Ptr};
 use core::any::Any;
 use num_traits::Float as NumFloat;
 
@@ -34,7 +34,6 @@ pub struct HairBxDF {
     pub s: Float,
     pub sin_2k_alpha: [Float; P_MAX],
     pub cos_2k_alpha: [Float; P_MAX],
-    pub colorspaces: DeviceStandardColorSpaces,
 }
 
 impl HairBxDF {
@@ -45,7 +44,6 @@ impl HairBxDF {
         beta_m: Float,
         beta_n: Float,
         alpha: Float,
-        colorspaces: DeviceStandardColorSpaces,
     ) -> Self {
         let mut sin_2k_alpha = [0.; P_MAX];
         let mut cos_2k_alpha = [0.; P_MAX];
@@ -67,7 +65,6 @@ impl HairBxDF {
             s: 0.,
             sin_2k_alpha,
             cos_2k_alpha,
-            colorspaces,
         }
     }
 
@@ -141,12 +138,25 @@ impl HairBxDF {
     pub fn sigma_a_from_concentration(
         ce: Float,
         cp: Float,
-        stdcs: DeviceStandardColorSpaces,
+        srgb: Ptr<RGBColorSpace>,
     ) -> RGBUnboundedSpectrum {
         let eumelanin_sigma_a = RGB::new(0.419, 0.697, 1.37);
         let pheomelanin_sigma_a = RGB::new(0.187, 0.4, 1.05);
         let sigma_a = ce * eumelanin_sigma_a + cp * pheomelanin_sigma_a;
-        RGBUnboundedSpectrum::new(&stdcs.srgb, sigma_a)
+        RGBUnboundedSpectrum::new(&srgb, sigma_a)
+    }
+
+    pub fn sigma_a_from_reflectance(c: SampledSpectrum, beta_n: Float) -> SampledSpectrum {
+        let mut sigma_a = SampledSpectrum::zero();
+        for i in 0..N_SPECTRUM_SAMPLES {
+            sigma_a[i] = square(
+                c[i].ln()
+                    / (5.969 - 0.215 * beta_n + 2.532 * square(beta_n) - 10.73 * beta_n.powf(3.)
+                        + 5.574 * beta_n.powf(4.)
+                        + 0.245 * beta_n.powf(5.)),
+            );
+        }
+        sigma_a
     }
 }
 

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/cameras/orthographic.rs

@@ -24,7 +24,7 @@ pub struct OrthographicCamera {
     pub dy_camera: Vector3f,
 }
 
-#[cfg(not(target_os = "cuda"))]
+#[cfg(not(gpu))]
 impl OrthographicCamera {
     pub fn new(
         base: CameraBase,

shared/src/cameras/perspective.rs

@@ -26,7 +26,7 @@ pub struct PerspectiveCamera {
     pub cos_total_width: Float,
 }
 
-#[cfg(not(target_os = "cuda"))]
+#[cfg(not(gpu))]
 impl PerspectiveCamera {
     pub fn new(
         base: CameraBase,

shared/src/core/bssrdf.rs

@@ -3,11 +3,11 @@ use crate::core::bsdf::BSDF;
 use crate::core::geometry::{Frame, Normal3f, Point2f, Point3f, Point3fi, Vector3f};
 use crate::core::interaction::{InteractionBase, ShadingGeom, SurfaceInteraction};
 use crate::core::shape::Shape;
-use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES};
-use crate::utils::math::{catmull_rom_weights, square};
-use crate::utils::sampling::sample_catmull_rom_2d;
 use crate::core::{LightIdx, MaterialIdx};
-use crate::{gvec_with_capacity, Float, GVec, PI, Ptr};
+use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum};
+use crate::utils::math::{catmull_rom_weights, invert_catmull_rom, square};
+use crate::utils::sampling::sample_catmull_rom_2d;
+use crate::{Float, GVec, PI, Ptr, gvec_with_capacity};
 use enum_dispatch::enum_dispatch;
 use num_traits::Float as NumFloat;
 
@@ -105,20 +105,19 @@ pub struct BSSRDFTable {
 
 impl BSSRDFTable {
     pub fn new(n_rho: usize, n_radius: usize) -> Self {
-        let rho_samples: GVec<Float> = gvec_with_capacity(n_rho);
-        let radius_samples: GVec<Float> = gvec_with_capacity(n_radius);
-        let profile: GVec<Float> = gvec_with_capacity(n_radius * n_rho);
-        let rho_eff: GVec<Float> = gvec_with_capacity(n_rho);
-        let profile_cdf: GVec<Float> = gvec_with_capacity(n_radius * n_rho);
+        let filled = |n: usize| {
+            let mut v: GVec<Float> = gvec_with_capacity(n);
+            v.resize(n, 0.);
+            v
+        };
         Self {
-            n_rho: n_rho.try_into().unwrap(),
-            n_radius: n_radius.try_into().unwrap(),
-            rho_samples,
-            radius_samples,
-            profile,
-            rho_eff,
-            profile_cdf,
-
+            n_rho: n_rho as u32,
+            n_radius: n_radius as u32,
+            rho_samples: filled(n_rho),
+            radius_samples: filled(n_radius),
+            profile: filled(n_rho * n_radius),
+            rho_eff: filled(n_rho),
+            profile_cdf: filled(n_rho * n_radius),
         }
     }
 
@@ -148,6 +147,22 @@ impl BSSRDFTable {
     }
 }
 
+pub fn subsurface_from_diffuse(
+    t: &BSSRDFTable,
+    rho_eff: &SampledSpectrum,
+    mfp: &SampledSpectrum,
+) -> (SampledSpectrum, SampledSpectrum) {
+    // (sigma_a, sigma_s)
+    let mut sigma_a = SampledSpectrum::zero();
+    let mut sigma_s = SampledSpectrum::zero();
+    for c in 0..N_SPECTRUM_SAMPLES {
+        let rho = invert_catmull_rom(&t.rho_samples, &t.rho_eff, rho_eff[c]);
+        sigma_s[c] = rho / mfp[c];
+        sigma_a[c] = (1. - rho) / mfp[c];
+    }
+    (sigma_a, sigma_s)
+}
+
 #[repr(C)]
 #[derive(Copy, Clone, Default, Debug)]
 pub struct BSSRDFProbeSegment {
@@ -236,7 +251,7 @@ impl TabulatedBSSRDF {
                         sr += weight
                             * self
                                 .table
-                                .eval_profile(rho_offset + j as u32, radius_offset + k as u32);
+                                .eval_profile((rho_offset + j as i32) as u32, (radius_offset + k as i32) as u32);
                     }
                 }
             }
@@ -287,14 +302,14 @@ impl TabulatedBSSRDF {
             for (j, rho_weight) in rho_weights.iter().enumerate() {
                 if *rho_weight != 0. {
                     // Update _rhoEff_ and _sr_ for wavelength
-                    rho_eff += rhoeff_samples[rho_offset as usize + j] * rho_weight;
+                    rho_eff += rhoeff_samples[(rho_offset + j as i32) as usize] * rho_weight;
 
                     // Fix: Use .iter().enumerate() for 'k'
                     for (k, radius_weight) in radius_weights.iter().enumerate() {
                         if *radius_weight != 0. {
                             sr += self
                                 .table
-                                .eval_profile(rho_offset + j as u32, radius_offset + k as u32)
+                                .eval_profile((rho_offset + j as i32) as u32, (radius_offset + k as i32) as u32)
                                 * rho_weight
                                 * radius_weight;
                         }

shared/src/core/camera.rs

@@ -137,7 +137,7 @@ pub trait CameraTrait {
     fn generate_ray(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option<CameraRay>;
 
     fn get_film(&self) -> &Film {
-        #[cfg(not(target_os = "cuda"))]
+        #[cfg(not(gpu))]
         {
             if self.base().film.is_null() {
                 panic!(

shared/src/core/color.rs

@@ -625,7 +625,8 @@ impl RGBSigmoidPolynomial {
     }
 
     pub fn evaluate(&self, lambda: Float) -> Float {
-        let eval = evaluate_polynomial(lambda, &[self.c0, self.c1, self.c2]);
+        // pbrt: `s(EvaluatePolynomial(lambda, c2, c1, c0))` -- c2 is the constant term.
+        let eval = evaluate_polynomial(lambda, &[self.c2, self.c1, self.c0]);
         Self::s(eval)
     }
 

shared/src/core/film.rs

@@ -96,7 +96,7 @@ impl RGBFilm {
     }
 
     pub fn get_sensor(&self) -> &PixelSensor {
-        #[cfg(not(target_os = "cuda"))]
+        #[cfg(not(gpu))]
         {
             if self.base.sensor.is_null() {
                 panic!(
@@ -203,7 +203,7 @@ impl RGBFilm {
 
 #[repr(C)]
 #[derive(Debug, Clone)]
-#[cfg_attr(target_os = "cuda", derive(Copy))]
+#[cfg_attr(gpu, derive(Copy))]
 pub struct GBufferPixel {
     pub rgb_sum: [AtomicFloat; 3],
     pub weight_sum: AtomicFloat,
@@ -240,7 +240,7 @@ impl Default for GBufferPixel {
 
 #[repr(C)]
 #[derive(Debug, Clone)]
-#[cfg_attr(target_os = "cuda", derive(Copy))]
+#[cfg_attr(gpu, derive(Copy))]
 pub struct GBufferFilm {
     pub base: FilmBase,
     pub output_from_render: AnimatedTransform,
@@ -294,7 +294,7 @@ impl GBufferFilm {
     }
 
     pub fn get_sensor(&self) -> &PixelSensor {
-        #[cfg(not(target_os = "cuda"))]
+        #[cfg(not(gpu))]
         {
             if self.base.sensor.is_null() {
                 panic!(
@@ -387,7 +387,7 @@ impl GBufferFilm {
 
 #[repr(C)]
 #[derive(Debug)]
-#[cfg_attr(target_os = "cuda", derive(Copy))]
+#[cfg_attr(gpu, derive(Copy))]
 pub struct SpectralPixel {
     pub rgb_sum: [AtomicFloat; 3],
     pub rgb_weight_sum: AtomicFloat,
@@ -419,7 +419,7 @@ impl Default for SpectralPixel {
 
 #[repr(C)]
 #[derive(Debug)]
-#[cfg_attr(target_os = "cuda", derive(Copy, Clone))]
+#[cfg_attr(gpu, derive(Copy, Clone))]
 pub struct SpectralFilm {
     pub base: FilmBase,
     pub lambda_min: Float,
@@ -609,7 +609,7 @@ pub struct FilmBase {
 
 #[repr(C)]
 #[derive(Debug)]
-#[cfg_attr(target_os = "cuda", derive(Copy, Clone))]
+#[cfg_attr(gpu, derive(Copy, Clone))]
 pub enum Film {
     RGB(RGBFilm),
     GBuffer(GBufferFilm),

shared/src/core/interaction.rs

@@ -590,7 +590,7 @@ impl SurfaceInteraction {
         }
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     pub fn set_intersection_properties(
         &mut self,
         mtl: MaterialIdx,

shared/src/core/light.rs

@@ -61,7 +61,7 @@ pub struct LightLiSample {
     pub p_light: Interaction,
 }
 
-#[cfg(not(target_os = "cuda"))]
+#[cfg(not(gpu))]
 impl LightLiSample {
     pub fn new(l: SampledSpectrum, wi: Vector3f, pdf: Float, p_light: Interaction) -> Self {
         Self {
@@ -188,7 +188,7 @@ pub struct LightBounds {
     pub two_sided: bool,
 }
 
-#[cfg(not(target_os = "cuda"))]
+#[cfg(not(gpu))]
 impl LightBounds {
     pub fn new(
         bounds: &Bounds3f,
@@ -327,13 +327,13 @@ pub trait LightTrait {
         self.base().light_type
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn bounds(&self) -> Option<LightBounds>;
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn preprocess(&mut self, scene_bounds: &Bounds3f);
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum;
 }
 

shared/src/core/medium.rs

@@ -99,7 +99,7 @@ pub struct MajorantGrid {
 
 
 impl MajorantGrid {
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     pub fn new(bounds: Bounds3f, res: Point3i) -> Self {
         let n_voxels = (res.x() * res.y() * res.z()) as usize;
         let voxels = gvec_with_capacity(n_voxels);

shared/src/core/sampler.rs

@@ -192,15 +192,18 @@ impl SamplerTrait for HaltonSampler {
     }
 
     fn get1d(&mut self) -> Float {
-        if self.dim + 1 >= PRIME_TABLE_SIZE as u32 {
+        // pbrt: `SampleDimension(dimension++)` -- POST-increment. Pre-incrementing makes
+        // the next Get2D() reuse the dimension this call just consumed.
+        if self.dim >= PRIME_TABLE_SIZE as u32 {
             self.dim = 2;
         }
+        let dim = self.dim;
         self.dim += 1;
-        self.sample_dimension(self.dim)
+        self.sample_dimension(dim)
     }
 
     fn get2d(&mut self) -> Point2f {
-        if self.dim > PRIME_TABLE_SIZE as u32 {
+        if self.dim + 1 >= PRIME_TABLE_SIZE as u32 {
             self.dim = 2;
         }
         let dim = self.dim;
@@ -301,8 +304,9 @@ impl SamplerTrait for StratifiedSampler {
             hash as u32,
         );
         self.dim += 2;
+        // pbrt: both the modulus and the divisor are xPixelSamples.
         let x = stratum % self.x_pixel_samples as u32;
-        let y = stratum / self.y_pixel_samples as u32;
+        let y = stratum / self.x_pixel_samples as u32;
         let dx = if self.jitter {
             self.rng.uniform::<Float>()
         } else {
@@ -387,6 +391,7 @@ impl SamplerTrait for PaddedSobolSampler {
             self.samples_per_pixel as u32,
             hash as u32,
         );
+        self.dim += 1;
         self.sample_dimension(0, index, (hash >> 32) as u32)
     }
     fn get2d(&mut self) -> Point2f {
@@ -512,7 +517,7 @@ impl SamplerTrait for SobolSampler {
         ) as Float;
         u[1] = clamp(
             u[1] * self.scale as Float - self.pixel[1] as Float,
-            1.,
+            0.,
             ONE_MINUS_EPSILON,
         ) as Float;
         u

shared/src/core/scattering.rs

@@ -1,9 +1,9 @@
 use crate::core::geometry::{
-    abs_cos_theta, cos2_theta, cos_phi, sin_phi, tan2_theta, Normal3f, Point2f, Vector2f, Vector3f,
-    VectorLike,
+    Normal3f, Point2f, Vector2f, Vector3f, VectorLike, abs_cos_theta, cos_phi, cos2_theta, sin_phi,
+    tan2_theta,
 };
-use crate::core::pbrt::{Float, PI};
-use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES};
+use crate::core::pbrt::{Float, INV_4_PI, PI};
+use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum};
 use crate::utils::math::{clamp, lerp, safe_sqrt, square};
 use crate::utils::sampling::sample_uniform_disk_polar;
 use num_traits::Float as NumFloat;
@@ -188,11 +188,9 @@ pub fn fresnel_moment1(eta: Float) -> Float {
     let eta4 = eta3 * eta;
     let eta5 = eta4 * eta;
     if eta < 1. {
-        return 0.45966 - 1.73965 * eta + 3.37668 * eta2 - 3.904945 * eta3 + 2.49277 * eta4
-            - 0.68441 * eta5;
+        0.45966 - 1.73965 * eta + 3.37668 * eta2 - 3.904945 * eta3 + 2.49277 * eta4 - 0.68441 * eta5
     } else {
-        return -4.61686 + 11.1136 * eta - 10.4646 * eta2 + 5.11455 * eta3 - 1.27198 * eta4
-            + 0.12746 * eta5;
+        -4.61686 + 11.1136 * eta - 10.4646 * eta2 + 5.11455 * eta3 - 1.27198 * eta4 + 0.12746 * eta5
     }
 }
 
@@ -203,18 +201,28 @@ pub fn fresnel_moment2(eta: Float) -> Float {
     let eta5 = eta4 * eta;
 
     if eta < 1. {
-        return 0.27614 - 0.87350 * eta + 1.12077 * eta2 - 0.65095 * eta3
-            + 0.07883 * eta4
-            + 0.04860 * eta5;
+        0.27614 - 0.87350 * eta + 1.12077 * eta2 - 0.65095 * eta3 + 0.07883 * eta4 + 0.04860 * eta5
     } else {
         let r_eta = 1. / eta;
         let r_eta2 = r_eta * r_eta;
         let r_eta3 = r_eta2 * r_eta;
 
-        return -547.033 + 45.3087 * r_eta3 - 218.725 * r_eta2 + 458.843 * r_eta + 404.557 * eta
+        -547.033 + 45.3087 * r_eta3 - 218.725 * r_eta2 + 458.843 * r_eta + 404.557 * eta
             - 189.519 * eta2
             + 54.9327 * eta3
             - 9.00603 * eta4
-            + 0.63942 * eta5;
+            + 0.63942 * eta5
     }
 }
+
+#[inline]
+pub fn henyey_greenstein(cos_theta: Float, mut g: Float) -> Float {
+    // The Henyey-Greenstein phase function isn't suitable for |g| \approx
+    // 1 so we clamp it before it becomes numerically instable. (It's an
+    // analogous situation to BSDFs: if the BSDF is perfectly specular, one
+    // should use one based on a Dirac delta distribution rather than a
+    // very smooth microfacet distribution...)
+    g = g.clamp(-0.99, 0.99);
+    let denom = 1. + square(g) + 2. * g * cos_theta;
+    INV_4_PI * (1. - square(g)) / (denom * safe_sqrt(denom))
+}

shared/src/core/spectrum.rs

@@ -36,6 +36,15 @@ pub enum Spectrum {
     RGBUnbounded(RGBUnboundedSpectrum),
 }
 
+/// `enum_dispatch` already generates `From<Variant> for Spectrum`, so wrapping a
+/// `ConstantSpectrum` etc. is `.into()`. Only the plain-`Float` hop is missing,
+/// and it is the one written most often at default-value sites.
+impl From<Float> for Spectrum {
+    fn from(c: Float) -> Self {
+        Spectrum::Constant(ConstantSpectrum::new(c))
+    }
+}
+
 impl<T: SpectrumTrait> SpectrumTrait for Ptr<T> {
     fn evaluate(&self, lambda: Float) -> Float {
         self.get().unwrap().evaluate(lambda)

shared/src/core/texture.rs

@@ -245,7 +245,7 @@ pub struct PointTransformMapping {
 }
 
 impl PointTransformMapping {
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     pub fn new(texture_from_render: Transform) -> Self {
         Self {
             texture_from_render,

shared/src/lights/diffuse.rs

@@ -130,7 +130,7 @@ impl LightTrait for DiffuseAreaLight {
         }
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         let mut l = SampledSpectrum::new(0.);
         if !self.image.is_null() {
@@ -153,10 +153,10 @@ impl LightTrait for DiffuseAreaLight {
         PI * two_side * self.area * l
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn preprocess(&mut self, _scene_bounds: &Bounds3f) {}
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn bounds(&self) -> Option<LightBounds> {
         let mut phi = 0.;
         if !self.image.is_null() {

shared/src/lights/goniometric.rs

@@ -1,5 +1,8 @@
-use crate::core::geometry::{Bounds3f, Normal3f, Point2f, Point2i, Point3f, Ray, Vector3f};
+use crate::core::geometry::{
+    Bounds3f, Normal3f, Point2f, Point2i, Point3f, Ray, Vector3f, VectorLike,
+};
 use crate::core::image::Image;
+use crate::core::interaction::{Interaction, InteractionBase, SimpleInteraction};
 use crate::core::light::{
     LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType,
 };
@@ -34,12 +37,20 @@ impl LightTrait for GoniometricLight {
 
     fn sample_li(
         &self,
-        _ctx: &LightSampleContext,
+        ctx: &LightSampleContext,
         _u: Point2f,
-        _lambda: &SampledWavelengths,
+        lambda: &SampledWavelengths,
         _allow_incomplete_pdf: bool,
     ) -> Option<LightLiSample> {
-        todo!()
+        let render_from_light = self.base().render_from_light;
+        let p = render_from_light.apply_to_point(Point3f::new(0., 0., 0.));
+        let wi = (p - ctx.p()).normalize();
+        let wl = render_from_light.apply_inverse_vector(-wi);
+        let li = self.i(wl, lambda) / p.distance_squared(ctx.p());
+        let base = InteractionBase::new_boundary(p, 0., self.base.medium_interface);
+        let intr = SimpleInteraction::new(base);
+
+        Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
     }
 
     fn pdf_li(
@@ -76,7 +87,7 @@ impl LightTrait for GoniometricLight {
         ))
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         let resolution = self.image.resolution();
         let mut sum_y = 0.;

shared/src/lights/infinite.rs

@@ -109,7 +109,7 @@ impl LightTrait for UniformInfiniteLight {
         None
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         4. * PI * PI * square(self.scene_radius) * self.scale * self.lemit.sample(&lambda)
     }
@@ -227,21 +227,17 @@ impl LightTrait for ImageInfiniteLight {
         self.image_le(uv, lambda)
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         let mut sum_l = SampledSpectrum::new(0.);
         let width = self.image.resolution().x();
         let height = self.image.resolution().y();
         for v in 0..height {
             for u in 0..width {
-                let mut rgb = RGB::default();
-                for c in 0..3 {
-                    rgb[c] = self.image.get_channel_with_wrap(
-                        Point2i::new(u, v),
-                        c,
-                        WrapMode::OctahedralSphere.into(),
-                    );
-                }
+                let rgb = RGB::from(self.image.get_channels_with_wrap::<3>(
+                    Point2i::new(u, v),
+                    WrapMode::OctahedralSphere.into(),
+                ));
                 sum_l += RGBIlluminantSpectrum::new(&self.image_color_space, rgb.clamp_zero())
                     .sample(&lambda);
             }
@@ -341,7 +337,7 @@ impl PortalInfiniteLight {
         (self.portal[1] - self.portal[0]).norm() * (self.portal[3] - self.portal[0]).norm()
     }
 
-    pub fn render_from_image(portal_frame: Frame, uv: Point2f) -> (Vector3f, Float) {
+    pub fn render_from_image_with(portal_frame: Frame, uv: Point2f) -> (Vector3f, Float) {
         let alpha = -PI / 2.0 + uv.x() * PI;
         let beta = -PI / 2.0 + uv.y() * PI;
 
@@ -354,6 +350,11 @@ impl PortalInfiniteLight {
 
         (portal_frame.from_local(w), duv_dw)
     }
+
+    #[inline]
+    pub fn render_from_image(&self, uv: Point2f) -> (Vector3f, Float) {
+        Self::render_from_image_with(self.portal_frame, uv)
+    }
 }
 
 impl LightTrait for PortalInfiniteLight {
@@ -370,7 +371,7 @@ impl LightTrait for PortalInfiniteLight {
     ) -> Option<LightLiSample> {
         let b = self.image_bounds(ctx.p())?;
         let (uv, map_pdf) = self.distribution.sample(u, b)?;
-        let (wi, duv_dw) = Self::render_from_image(self.portal_frame, uv);
+        let (wi, duv_dw) = self.render_from_image(uv);
         if duv_dw == 0. {
             return None;
         }
@@ -403,17 +404,34 @@ impl LightTrait for PortalInfiniteLight {
         }
     }
 
-    #[cfg(not(target_os = "cuda"))]
-    fn phi(&self, _lambda: SampledWavelengths) -> SampledSpectrum {
-        todo!()
+    #[cfg(not(gpu))]
+    fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
+        let mut sum_l = SampledSpectrum::new(0.);
+        let width = self.image.resolution().x();
+        let height = self.image.resolution().y();
+        for y in 0..height {
+            for x in 0..width {
+                let rgb = RGB::from(self.image.get_channels::<3>(Point2i::new(x, y)));
+                let st = Point2f::new(
+                    (x as Float + 0.5) / width as Float,
+                    (y as Float + 0.5) / height as Float,
+                );
+                let (_, duv_dw) = self.render_from_image(st);
+                sum_l += RGBIlluminantSpectrum::new(&self.image_color_space, rgb.clamp_zero())
+                    .sample(&lambda)
+                    / duv_dw;
+            }
+        }
+
+        self.scale * self.area() * sum_l / (width * height) as Float
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn preprocess(&mut self, scene_bounds: &Bounds3f) {
         (self.scene_center, self.scene_radius) = scene_bounds.bounding_sphere();
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn bounds(&self) -> Option<LightBounds> {
         None
     }

shared/src/lights/point.rs

@@ -51,14 +51,14 @@ impl LightTrait for PointLight {
         0.
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         4. * PI * self.scale * self.i.sample(&lambda)
     }
 
     fn preprocess(&mut self, _scene_bounds: &Bounds3f) {}
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn bounds(&self) -> Option<LightBounds> {
         let p = self
             .base

shared/src/lights/projection.rs

@@ -4,6 +4,7 @@ use crate::core::geometry::{
     Bounds2f, Bounds3f, Normal3f, Point2f, Point2i, Point3f, Ray, Vector3f, VectorLike, cos_theta,
 };
 use crate::core::image::Image;
+use crate::core::interaction::{Interaction, InteractionBase, SimpleInteraction};
 use crate::core::light::{
     LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType,
 };
@@ -33,7 +34,7 @@ pub struct ProjectionLight {
 }
 
 impl ProjectionLight {
-    pub fn i(&self, w: Vector3f, lambda: SampledWavelengths) -> SampledSpectrum {
+    pub fn i(&self, w: Vector3f, lambda: &SampledWavelengths) -> SampledSpectrum {
         if w.z() < self.hither {
             return SampledSpectrum::new(0.);
         }
@@ -44,10 +45,10 @@ impl ProjectionLight {
         let uv = Point2f::from(self.screen_bounds.offset(&Point2f::new(ps.x(), ps.y())));
         let mut rgb = RGB::default();
         for c in 0..3 {
-            rgb[c] = self.image.lookup_nearest_channel(uv, c as i32);
+            rgb[c] = self.image.lookup_nearest_channel(uv, c);
         }
-        let s = RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero());
-        self.scale * s.sample(&lambda)
+        let s = RGBIlluminantSpectrum::new(&self.image_color_space, rgb.clamp_zero());
+        self.scale * s.sample(lambda)
     }
 }
 
@@ -58,12 +59,23 @@ impl LightTrait for ProjectionLight {
 
     fn sample_li(
         &self,
-        _ctx: &LightSampleContext,
+        ctx: &LightSampleContext,
         _u: Point2f,
-        _lambda: &SampledWavelengths,
+        lambda: &SampledWavelengths,
         _allow_incomplete_pdf: bool,
     ) -> Option<LightLiSample> {
-        todo!()
+        let render_from_light = self.base().render_from_light;
+        let p = render_from_light.apply_to_point(Point3f::new(0., 0., 0.));
+        let wi = (p - ctx.p()).normalize();
+        let wl = render_from_light.apply_inverse_vector(-wi);
+        let li = self.i(wl, lambda) / p.distance_squared(ctx.p());
+        if li.is_black() {
+            return None;
+        }
+        let base = InteractionBase::new_boundary(p, 0., self.base.medium_interface);
+        let intr = SimpleInteraction::new(base);
+
+        Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
     }
 
     fn pdf_li(
@@ -72,7 +84,7 @@ impl LightTrait for ProjectionLight {
         _wi: Vector3f,
         _allow_incomplete_pdf: bool,
     ) -> Float {
-        todo!()
+        0.
     }
 
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
@@ -93,10 +105,10 @@ impl LightTrait for ProjectionLight {
                 let dwda = cos_theta(w).powi(3);
                 let mut rgb = RGB::default();
                 for c in 0..3 {
-                    rgb[c] = self.image.get_channel(Point2i::new(x, y), c as i32);
+                    rgb[c] = self.image.get_channel(Point2i::new(x, y), c);
                 }
 
-                let s = RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero());
+                let s = RGBIlluminantSpectrum::new(&self.image_color_space, rgb.clamp_zero());
                 sum += s.sample(&lambda) * dwda;
             }
         }

shared/src/lights/spot.rs

@@ -65,7 +65,7 @@ impl LightTrait for SpotLight {
         0.
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         self.scale
             * self.iemit.sample(&lambda)

shared/src/materials/complex.rs

@@ -4,53 +4,57 @@ use crate::bxdfs::{
     MeasuredBxDF, MeasuredBxDFData,
 };
 use crate::core::bsdf::BSDF;
-use crate::core::bssrdf::{BSSRDF, BSSRDFTable};
+use crate::core::bssrdf::{BSSRDF, BSSRDFTable, TabulatedBSSRDF, subsurface_from_diffuse};
 use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{FloatTexture, SpectrumTexture, TextureEvaluator};
-use crate::spectra::{SampledSpectrum, SampledWavelengths};
+use crate::spectra::{RGBColorSpace, SampledSpectrum, SampledWavelengths};
 use crate::textures::SpectrumMixTexture;
 use crate::utils::Ptr;
 use crate::utils::math::clamp;
 
+#[repr(C)]
+#[derive(Copy, Clone, Debug)]
+pub enum HairAbsorption {
+    SigmaA(Ptr<SpectrumTexture>),
+    Color(Ptr<SpectrumTexture>),
+    Melanin {
+        eumelanin: Ptr<FloatTexture>,
+        pheomelanin: Ptr<FloatTexture>,
+    },
+}
+
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct HairMaterial {
-    pub sigma_a: Ptr<SpectrumTexture>,
-    pub color: Ptr<SpectrumTexture>,
-    pub eumelanin: Ptr<FloatTexture>,
-    pub pheomelanin: Ptr<FloatTexture>,
+    pub hair_absorption: HairAbsorption,
     pub eta: Ptr<FloatTexture>,
     pub beta_m: Ptr<FloatTexture>,
     pub beta_n: Ptr<FloatTexture>,
     pub alpha: Ptr<FloatTexture>,
+    pub colorspace: Ptr<RGBColorSpace>,
 }
 
 impl HairMaterial {
-    #[cfg(not(target_os = "cuda"))]
-    #[allow(clippy::too_many_arguments)]
+    #[cfg(not(gpu))]
     pub fn new(
-        sigma_a: Ptr<SpectrumTexture>,
-        color: Ptr<SpectrumTexture>,
-        eumelanin: Ptr<FloatTexture>,
-        pheomelanin: Ptr<FloatTexture>,
+        hair_absorption: HairAbsorption,
         eta: Ptr<FloatTexture>,
         beta_m: Ptr<FloatTexture>,
         beta_n: Ptr<FloatTexture>,
         alpha: Ptr<FloatTexture>,
+        colorspace: Ptr<RGBColorSpace>,
     ) -> Self {
         Self {
-            sigma_a,
-            color,
-            eumelanin,
-            pheomelanin,
+            hair_absorption,
             eta,
             beta_m,
             beta_n,
             alpha,
+            colorspace,
         }
     }
 }
@@ -58,27 +62,86 @@ impl HairMaterial {
 impl MaterialTrait for HairMaterial {
     fn get_bsdf<T: TextureEvaluator>(
         &self,
-        _tex_eval: &T,
-        _ctx: &MaterialEvalContext,
-        _lambda: &mut SampledWavelengths,
+        tex_eval: &T,
+        ctx: &MaterialEvalContext,
+        lambda: &mut SampledWavelengths,
     ) -> BSDF {
-        todo!()
+        let bm = tex_eval.evaluate_float(&self.beta_m, ctx).clamp(1e-2, 1.0);
+        let bn = tex_eval.evaluate_float(&self.beta_n, ctx).clamp(1e-2, 1.0);
+        let a = tex_eval.evaluate_float(&self.alpha, ctx);
+        let e = tex_eval.evaluate_float(&self.eta, ctx);
+        let sig_a = match self.hair_absorption {
+            // Absorption coefficient, not reflectance, can be larger than 1
+            HairAbsorption::SigmaA(sigma_a) => {
+                SampledSpectrum::clamp_zero(&tex_eval.evaluate_spectrum(&sigma_a, ctx, lambda))
+            }
+            HairAbsorption::Color(color) => {
+                let c = SampledSpectrum::clamp(
+                    &tex_eval.evaluate_spectrum(&color, ctx, lambda),
+                    0.,
+                    1.,
+                );
+                HairBxDF::sigma_a_from_reflectance(c, bn)
+            }
+            HairAbsorption::Melanin {
+                eumelanin,
+                pheomelanin,
+            } => {
+                debug_assert!(!eumelanin.is_null() || !pheomelanin.is_null());
+                let eu = if !eumelanin.is_null() {
+                    tex_eval.evaluate_float(&eumelanin, ctx)
+                } else {
+                    0.
+                };
+                let pheo = if !pheomelanin.is_null() {
+                    tex_eval.evaluate_float(&pheomelanin, ctx)
+                } else {
+                    0.
+                };
+
+                HairBxDF::sigma_a_from_concentration(eu.max(0.0), pheo.max(0.0), self.colorspace)
+                    .sample(lambda)
+            }
+        };
+
+        let h = -1. + 2. * ctx.uv[1];
+        let bxdf = BxDF::Hair(HairBxDF::new(h, e, sig_a, bm, bn, a));
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
+
     fn get_bssrdf<T>(
         &self,
         _tex_eval: &T,
         _ctx: &MaterialEvalContext,
         _lambda: &SampledWavelengths,
     ) -> Option<BSSRDF> {
-        todo!()
+        None
     }
 
-    fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool {
-        todo!()
+    fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool {
+        match self.hair_absorption {
+            HairAbsorption::SigmaA(t) | HairAbsorption::Color(t) => {
+                tex_eval.can_evaluate(&[self.eta, self.beta_m, self.beta_n, self.alpha], &[t])
+            }
+            HairAbsorption::Melanin {
+                eumelanin,
+                pheomelanin,
+            } => tex_eval.can_evaluate(
+                &[
+                    self.eta,
+                    self.beta_m,
+                    self.beta_n,
+                    self.alpha,
+                    eumelanin,
+                    pheomelanin,
+                ],
+                &[],
+            ),
+        }
     }
 
     fn get_normal_map(&self) -> Option<&Image> {
-        todo!()
+        None
     }
 
     fn get_displacement(&self) -> Ptr<FloatTexture> {
@@ -102,11 +165,11 @@ impl MaterialTrait for MeasuredMaterial {
     fn get_bsdf<T: TextureEvaluator>(
         &self,
         _tex_eval: &T,
-        _ctx: &MaterialEvalContext,
-        _lambda: &mut SampledWavelengths,
+        ctx: &MaterialEvalContext,
+        lambda: &mut SampledWavelengths,
     ) -> BSDF {
-        // MeasuredBxDF::new(&self.brdf, lambda)
-        todo!()
+        let bxdf = BxDF::Measured(MeasuredBxDF::new(&self.brdf, lambda));
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
     fn get_bssrdf<T>(
@@ -135,15 +198,25 @@ impl MaterialTrait for MeasuredMaterial {
     }
 }
 
+#[repr(C)]
+#[derive(Copy, Clone, Debug)]
+pub enum SubsurfaceScattering {
+    Coefficients {
+        sigma_a: Ptr<SpectrumTexture>,
+        sigma_s: Ptr<SpectrumTexture>,
+    },
+    Reflectance {
+        reflectance: Ptr<SpectrumTexture>,
+        mfp: Ptr<SpectrumTexture>,
+    },
+}
+
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct SubsurfaceMaterial {
     pub normal_map: Ptr<Image>,
     pub displacement: Ptr<FloatTexture>,
-    pub sigma_a: Ptr<SpectrumTexture>,
-    pub sigma_s: Ptr<SpectrumMixTexture>,
-    pub reflectance: Ptr<SpectrumMixTexture>,
-    pub mfp: Ptr<SpectrumMixTexture>,
+    pub scattering: SubsurfaceScattering,
     pub eta: Float,
     pub scale: Float,
     pub u_roughness: Ptr<FloatTexture>,
@@ -155,31 +228,78 @@ pub struct SubsurfaceMaterial {
 impl MaterialTrait for SubsurfaceMaterial {
     fn get_bsdf<T: TextureEvaluator>(
         &self,
-        _tex_eval: &T,
-        _ctx: &MaterialEvalContext,
+        tex_eval: &T,
+        ctx: &MaterialEvalContext,
         _lambda: &mut SampledWavelengths,
     ) -> BSDF {
-        todo!()
-    }
-    fn get_bssrdf<T>(
-        &self,
-        _tex_eval: &T,
-        _ctx: &MaterialEvalContext,
-        _lambda: &SampledWavelengths,
-    ) -> Option<BSSRDF> {
-        todo!()
+        let mut u_rough = tex_eval.evaluate_float(&self.u_roughness, ctx);
+        let mut v_rough = tex_eval.evaluate_float(&self.v_roughness, ctx);
+        if self.remap_roughness {
+            u_rough = TrowbridgeReitzDistribution::roughness_to_alpha(u_rough);
+            v_rough = TrowbridgeReitzDistribution::roughness_to_alpha(v_rough);
+        }
+
+        let distrib = TrowbridgeReitzDistribution::new(u_rough, v_rough);
+        let bxdf = BxDF::Dielectric(DielectricBxDF::new(self.eta, distrib));
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
-    fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool {
-        todo!()
+    fn get_bssrdf<T: TextureEvaluator>(
+        &self,
+        tex_eval: &T,
+        ctx: &MaterialEvalContext,
+        lambda: &SampledWavelengths,
+    ) -> Option<BSSRDF> {
+        let (sig_a, sig_s) = match self.scattering {
+            SubsurfaceScattering::Coefficients { sigma_a, sigma_s } => {
+                let s_a = SampledSpectrum::clamp_zero(
+                    &(self.scale * tex_eval.evaluate_spectrum(&sigma_a, ctx, lambda)),
+                );
+                let s_s = SampledSpectrum::clamp_zero(
+                    &(self.scale * tex_eval.evaluate_spectrum(&sigma_s, ctx, lambda)),
+                );
+                (s_a, s_s)
+            }
+            SubsurfaceScattering::Reflectance { reflectance, mfp } => {
+                debug_assert!(!reflectance.is_null() && !mfp.is_null());
+                let mfree =
+                    SampledSpectrum::clamp_zero(&tex_eval.evaluate_spectrum(&mfp, ctx, lambda));
+                let r = SampledSpectrum::clamp_zero(&tex_eval.evaluate_spectrum(
+                    &reflectance,
+                    ctx,
+                    lambda,
+                ));
+                subsurface_from_diffuse(&self.table, &r, &mfree)
+            }
+        };
+
+        Some(BSSRDF::Tabulated(TabulatedBSSRDF::new(
+            ctx.p,
+            ctx.wo,
+            ctx.ns,
+            self.eta,
+            &sig_a,
+            &sig_s,
+            &self.table,
+        )))
+    }
+
+    fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool {
+        // Slight divergence from PBRT, we check against reflectance and mfp as well in reflectance
+        // mode. Test thoroughly, keep as is for now (20260902)
+        let spectra = match self.scattering {
+            SubsurfaceScattering::Coefficients { sigma_a, sigma_s } => [sigma_a, sigma_s],
+            SubsurfaceScattering::Reflectance { reflectance, mfp } => [reflectance, mfp],
+        };
+        tex_eval.can_evaluate(&[self.u_roughness, self.v_roughness], &spectra)
     }
 
     fn get_normal_map(&self) -> Option<&Image> {
-        todo!()
+        Some(&*self.normal_map)
     }
 
     fn get_displacement(&self) -> Ptr<FloatTexture> {
-        todo!()
+        self.displacement
     }
 
     fn has_subsurface_scattering(&self) -> bool {

shared/src/materials/conductor.rs

@@ -27,6 +27,7 @@ pub struct ConductorMaterial {
 }
 
 impl ConductorMaterial {
+    #[allow(clippy::too_many_arguments)]
     pub fn new(
         normal_map: Ptr<Image>,
         reflectance: Ptr<SpectrumTexture>,
@@ -92,7 +93,7 @@ impl MaterialTrait for ConductorMaterial {
         _ctx: &MaterialEvalContext,
         _lambda: &SampledWavelengths,
     ) -> Option<BSSRDF> {
-        todo!()
+        None
     }
     fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool {
         tex_eval.can_evaluate(

shared/src/shapes/bilinear.rs

@@ -98,7 +98,7 @@ impl BilinearPatchShape {
         Some([mesh.n[v0], mesh.n[v1], mesh.n[v2], mesh.n[v3]])
     }
 
-    #[cfg(not(target_os = "cuda"))]
+    #[cfg(not(gpu))]
     pub fn new(mesh: Ptr<BilinearPatchMesh>, blp_index: i32) -> Self {
         let mut bp = BilinearPatchShape {
             mesh,

shared/src/shapes/triangle.rs

@@ -1,8 +1,8 @@
-use crate::core::geometry::{spherical_triangle_area, SqrtExt, Tuple, VectorLike};
 use crate::core::geometry::{
     Bounds3f, DirectionCone, Normal, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector2f, Vector3,
     Vector3f,
 };
+use crate::core::geometry::{SqrtExt, Tuple, VectorLike, spherical_triangle_area};
 use crate::core::interaction::{
     Interaction, InteractionBase, InteractionTrait, SimpleInteraction, SurfaceInteraction,
 };
@@ -13,7 +13,7 @@ use crate::utils::sampling::{
     bilinear_pdf, invert_spherical_triangle_sample, sample_bilinear, sample_spherical_triangle,
     sample_uniform_triangle,
 };
-use crate::{gamma, Float, GVec, Ptr};
+use crate::{Float, GVec, Ptr, gamma};
 
 #[repr(C)]
 #[derive(Debug, Clone, Copy)]
@@ -198,34 +198,34 @@ impl TriangleShape {
 
         // Ensure that computed triangle $t$ is conservatively greater than zero
         // Compute $\delta_z$ term for triangle $t$ error bounds
-        let maxZt = Vector3f::new(p0t.z(), p1t.z(), p2t.z())
+        let max_zt = Vector3f::new(p0t.z(), p1t.z(), p2t.z())
             .abs()
             .max_component_value();
-        let deltaZ = gamma(3) * maxZt;
+        let delta_z = gamma(3) * max_zt;
 
         // Compute $\delta_x$ and $\delta_y$ terms for triangle $t$ error bounds
-        let maxXt = Vector3f::new(p0t.x(), p1t.x(), p2t.x())
+        let max_xt = Vector3f::new(p0t.x(), p1t.x(), p2t.x())
             .abs()
             .max_component_value();
-        let maxYt = Vector3f::new(p0t.y(), p1t.y(), p2t.y())
+        let max_yt = Vector3f::new(p0t.y(), p1t.y(), p2t.y())
             .abs()
             .max_component_value();
-        let deltaX = gamma(5) * (maxXt + maxZt);
-        let deltaY = gamma(5) * (maxYt + maxZt);
+        let delta_x = gamma(5) * (max_xt + max_zt);
+        let delta_y = gamma(5) * (max_yt + max_zt);
 
         // Compute $\delta_e$ term for triangle $t$ error bounds
-        let deltaE = 2. * (gamma(2) * maxXt * maxYt + deltaY * maxXt + deltaX * maxYt);
+        let delta_e = 2. * (gamma(2) * max_xt * max_yt + delta_y * max_xt + delta_x * max_yt);
 
         // Compute $\delta_t$ term for triangle $t$ error bounds and check _t_
-        let maxE = Vector3f::new(e0, e1, e2).abs().max_component_value();
-        let deltaT =
-            3. * (gamma(3) * maxE * maxZt + deltaE * maxZt + deltaZ * maxE) * inv_det.abs();
-        if t <= deltaT {
+        let max_e = Vector3f::new(e0, e1, e2).abs().max_component_value();
+        let delta_t =
+            3. * (gamma(3) * max_e * max_zt + delta_e * max_zt + delta_z * max_e) * inv_det.abs();
+        if t <= delta_t {
             return None;
         }
 
         // Return _TriangleIntersection_ for intersection
-        return Some(TriangleIntersection { b0, b1, b2, t });
+        Some(TriangleIntersection { b0, b1, b2, t })
     }
 
     fn interaction_from_intersection(
@@ -345,11 +345,11 @@ impl TriangleShape {
 
         let mut ts = ns.cross(ss.into());
         if ts.norm_squared() > 0.0 {
-            ss = ts.cross(ns.into()).into();
+            ss = ts.cross(ns).into();
         } else {
             let (s, t) = ns.coordinate_system();
             ss = s.into();
-            ts = t.into();
+            ts = t;
         }
 
         let (dndu, dndv) = if let Some(normals) = self.get_shading_normals() {
@@ -539,7 +539,8 @@ impl ShapeTrait for TriangleShape {
 
     fn intersect(&self, ray: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let [p0, p1, p2] = self.get_points();
-        let tri_isect = self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY), p0, p1, p2)?;
+        let tri_isect =
+            self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY), p0, p1, p2)?;
         let intr = self.interaction_from_intersection(tri_isect, ray.time, -ray.d);
         Some(ShapeIntersection::new(intr, tri_isect.t))
     }

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<Float> = gvec_with_capacity(n + 2);
+        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;
+        }
     }
 }
 

shared/src/textures/bilerp.rs

@@ -39,8 +39,12 @@ impl FloatBilerpTexture {
         }
     }
 
-    pub fn evaluate(&self, _ctx: &TextureEvalContext) -> Float {
-        todo!()
+    pub fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
+        let c = self.mapping.map(ctx);
+        (1. - c.st[0]) * (1. - c.st[1]) * self.v00
+            + c.st[0] * (1. - c.st[1]) * self.v10
+            + (1. - c.st[0]) * c.st[1] * self.v01
+            + c.st[0] * c.st[1] * self.v11
     }
 }
 

shared/src/textures/marble.rs

@@ -22,8 +22,25 @@ pub struct MarbleTexture {
     pub colorspace: Ptr<RGBColorSpace>,
 }
 
-
 impl MarbleTexture {
+    pub fn new(
+        mapping: TextureMapping3D,
+        octaves: i32,
+        omega: Float,
+        scale: Float,
+        variation: Float,
+        colorspace: Ptr<RGBColorSpace>,
+    ) -> Self {
+        Self {
+            mapping,
+            octaves: octaves.try_into().unwrap(),
+            omega,
+            scale,
+            variation,
+            colorspace,
+        }
+    }
+
     pub fn evaluate(
         &self,
         ctx: &TextureEvalContext,
@@ -61,6 +78,6 @@ impl MarbleTexture {
         let (rgb_vec, _) = evaluate_cubic_bezier(&colors[first_idx..first_idx + 4], t_segment);
         let rgb = RGB::new(rgb_vec.x() * 1.5, rgb_vec.y() * 1.5, rgb_vec.z() * 1.5);
 
-        RGBAlbedoSpectrum::new(&*self.colorspace, rgb).sample(lambda)
+        RGBAlbedoSpectrum::new(&self.colorspace, rgb).sample(lambda)
     }
 }

shared/src/textures/windy.rs

@@ -10,6 +10,9 @@ pub struct WindyTexture {
 }
 
 impl WindyTexture {
+    pub fn new(mapping: TextureMapping3D) -> Self {
+        Self { mapping }
+    }
     pub fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
         let c = self.mapping.map(ctx);
         let wind_strength = fbm(

shared/src/textures/wrinkled.rs

@@ -11,6 +11,13 @@ pub struct WrinkledTexture {
 }
 
 impl WrinkledTexture {
+    pub fn new(mapping: TextureMapping3D, octaves: u32, omega: Float) -> Self {
+        Self {
+            mapping,
+            octaves,
+            omega,
+        }
+    }
     pub fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
         let c = self.mapping.map(ctx);
         turbulence(c.p, c.dpdx, c.dpdy, self.omega, self.octaves)

shared/src/utils/math.rs

@@ -1,9 +1,9 @@
 use crate::core::color::{RGB, XYZ};
 use crate::core::geometry::{Lerp, MulAdd, Point, Point2f, Point2i, Vector, Vector3f, VectorLike};
 use crate::core::pbrt::{Float, FloatBitOps, FloatBits, ONE_MINUS_EPSILON, PI, PI_OVER_4};
-use crate::utils::gpu_array_from_fn;
 use crate::utils::hash::{hash_buffer, mix_bits};
 use crate::utils::sobol::{SOBOL_MATRICES_32, VDC_SOBOL_MATRICES, VDC_SOBOL_MATRICES_INV};
+use crate::utils::{find_interval, gpu_array_from_fn};
 use crate::{GVec, Ptr, gvec, gvec_with_capacity};
 use core::fmt::{self, Display, Write};
 use core::iter::{Product, Sum};
@@ -61,11 +61,14 @@ where
     T::lerp(t, a, b)
 }
 
+/// pbrt `EvaluatePolynomial(t, c, cRemaining...) = FMA(t, EvaluatePolynomial(t, cRemaining...), c)`.
+/// `coeffs[0]` is the *constant* term and `coeffs[n-1]` the highest-degree one, so
+/// evaluation folds from the back.
 #[inline]
 pub fn evaluate_polynomial(t: Float, coeffs: &[Float]) -> Float {
     assert!(!coeffs.is_empty());
-    let mut result = coeffs[0];
-    for &c in &coeffs[1..] {
+    let mut result = coeffs[coeffs.len() - 1];
+    for &c in coeffs[..coeffs.len() - 1].iter().rev() {
         result = num_traits::Float::mul_add(t, result, c);
     }
     result
@@ -96,7 +99,7 @@ where
 
 #[inline]
 pub fn safe_sqrt(x: Float) -> Float {
-    assert!(x > -1e-3);
+    debug_assert!(x > -1e-3);
     0.0_f32.max(x).sqrt()
 }
 
@@ -104,20 +107,14 @@ pub fn safe_sqrt(x: Float) -> Float {
 pub fn safe_asin<T: NumFloat>(x: T) -> T {
     let epsilon = T::from(0.0001).unwrap();
     let one = T::one();
-    if x >= -(one + epsilon) && x <= one + epsilon {
-        clamp(x, -one, one).asin()
-    } else {
-        panic!("Not valid value for asin")
-    }
+    debug_assert!(x >= -(one + epsilon) && x <= one + epsilon);
+    clamp(x, -one, one).asin()
 }
 
 #[inline]
 pub fn safe_acos(x: Float) -> Float {
-    if (-1.001..1.001).contains(&x) {
-        clamp(x, -1., 1.).acos()
-    } else {
-        panic!("Not valid value for acos")
-    }
+    debug_assert!((-1.001..1.001).contains(&x));
+    clamp(x, -1., 1.).acos()
 }
 
 #[inline]
@@ -349,7 +346,93 @@ pub fn wrap_equal_area_square(uv: &mut Point2f) -> Point2f {
     *uv
 }
 
-pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(u32, [Float; 4])> {
+pub fn integrate_catmull_rom(nodes: &[Float], f: &[Float], cdf: &mut [Float]) -> Float {
+    debug_assert_eq!(nodes.len(), f.len());
+    debug_assert_eq!(cdf.len(), nodes.len());
+    let mut sum = 0.;
+    cdf[0] = 0.;
+    for i in 0..nodes.len() - 1 {
+        let x0 = nodes[i];
+        let x1 = nodes[i + 1];
+        let f0 = f[i];
+        let f1 = f[i + 1];
+        let width = x1 - x0;
+
+        // Approximate derivatives using finite differences
+        let d0 = if i > 0 {
+            width * (f1 - f[i - 1]) / (x1 - nodes[i - 1])
+        } else {
+            f1 - f0
+        };
+        let d1 = if i + 2 < nodes.len() {
+            width * (f[i + 2] - f0) / (nodes[i + 2] - x0)
+        } else {
+            f1 - f0
+        };
+
+        // Keep a running sum and build a cumulative distribution function
+        sum += width * ((f0 + f1) / 2. + (d0 - d1) / 12.);
+        cdf[i + 1] = sum;
+    }
+    sum
+}
+
+pub fn invert_catmull_rom(nodes: &[Float], f: &[Float], u: Float) -> Float {
+    // Stop when _u_ is out of bounds
+    if !(u > f[0]) {
+        return nodes[0];
+    } else if !(u < f[f.len() - 1]) {
+        return nodes[nodes.len() - 1];
+    }
+
+    // Map _u_ to a spline interval by inverting _f_
+    let i = find_interval(f.len() as u32, |j| f[j as usize] <= u) as usize;
+
+    // Look up $x_i$ and function values of spline segment _i_
+    let x0 = nodes[i];
+    let x1 = nodes[i + 1];
+    let f0 = f[i];
+    let f1 = f[i + 1];
+    let width = x1 - x0;
+
+    // Approximate derivatives using finite differences
+    let d0 = if i > 0 {
+        width * (f1 - f[i - 1]) / (x1 - nodes[i - 1])
+    } else {
+        f1 - f0
+    };
+    let d1 = if i + 2 < nodes.len() {
+        width * (f[i + 2] - f0) / (nodes[i + 2] - x0)
+    } else {
+        f1 - f0
+    };
+
+    // Invert the spline interpolant using Newton-Bisection
+    let eval = |t: Float| -> (Float, Float) {
+        // Compute powers of _t_
+        let t2 = t * t;
+        let t3 = t2 * t;
+
+        // Set _Fhat_ using Equation (\ref{eq:cubicspline-as-basisfunctions})
+        let f_cap_hat = (2. * t3 - 3. * t2 + 1.) * f0
+            + (-2. * t3 + 3. * t2) * f1
+            + (t3 - 2. * t2 + t) * d0
+            + (t3 - t2) * d1;
+
+        // Set _fhat_ using Equation (\ref{eq:cubicspline-derivative})
+        let f_hat = (6. * t2 - 6. * t) * f0
+            + (-6. * t2 + 6. * t) * f1
+            + (3. * t2 - 4. * t + 1.) * d0
+            + (3. * t2 - 2. * t) * d1;
+
+        return (f_cap_hat - u, f_hat);
+    };
+
+    let t = newton_bisection(0., 1., eval);
+    return x0 + t * width;
+}
+
+pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(i32, [Float; 4])> {
     if nodes.len() < 4 {
         return None;
     }
@@ -367,7 +450,10 @@ pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(u32, [Float; 4]
     // Safety clamp (though bounds check above handles most cases)
     let idx = idx.min(nodes.len() - 2);
 
-    let offset = idx.saturating_sub(1); // The C++ code uses idx - 1 for the offset
+    // pbrt: `*offset = idx - 1`, which is legitimately -1 when idx == 0. weights[0] is
+    // then 0, so callers never dereference the out-of-range slot -- but the remaining
+    // three weights must still line up with nodes[idx-1+i].
+    let offset = idx as i32 - 1;
     let x0 = nodes[idx];
     let x1 = nodes[idx + 1];
 
@@ -406,7 +492,7 @@ pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(u32, [Float; 4]
         weights[3] = 0.0;
     }
 
-    Some((offset as u32, weights))
+    Some((offset, weights))
 }
 
 pub fn equal_area_sphere_to_square(d: Vector3f) -> Point2f {
@@ -461,7 +547,7 @@ pub fn equal_area_square_to_sphere(p: Point2f) -> Vector3f {
 
     // Compute angle \phi for square to sphere mapping
     let mut phi = if r == 0. { 1. } else { (vp - up) / r + 1. };
-    phi /= PI_OVER_4;
+    phi *= PI_OVER_4;
 
     // Find z for spherical direction
     let z = (1. - square(r)).copysign(signed_distance);
@@ -500,12 +586,13 @@ pub fn gaussian_integral(x0: Float, x1: Float, mu: Float, sigma: Float) -> Float
     0.5 * (erf((mu - x0) / sigma_root2) - erf((mu - x1) / sigma_root2))
 }
 
+#[inline(always)]
 pub fn sample_linear(u: Float, a: Float, b: Float) -> Float {
     assert!(a >= 0. && b >= 0.);
     if u == 0. && a == 0. {
         return 0.;
     }
-    let x = u * (a + b) / (a + (lerp(u, square(a), square(b))));
+    let x = u * (a + b) / (a + lerp(u, square(a), square(b)).sqrt());
     x.min(ONE_MINUS_EPSILON)
 }
 
@@ -615,12 +702,13 @@ pub fn sample_discrete(
 }
 
 pub fn sample_tent(u: Float, r: Float) -> Float {
+    // pbrt rewrites `u` in place with the remapped sample and feeds *that* to SampleLinear.
     let mut u_remapped = 0.0;
     let offset = sample_discrete(&[0.5, 0.5], u, None, Some(&mut u_remapped));
     if offset == 0 {
-        -r + r * sample_linear(u, 0., 1.)
+        -r + r * sample_linear(u_remapped, 0., 1.)
     } else {
-        r * sample_linear(u, 1., 0.)
+        r * sample_linear(u_remapped, 1., 0.)
     }
 }
 
@@ -650,7 +738,7 @@ pub fn encode_morton_2(x: u32, y: u32) -> u64 {
 }
 
 pub fn encode_morton_3(x: Float, y: Float, z: Float) -> u32 {
-    (left_shift3(x as u32) << 2) | (left_shift3(y as u32) << 1) | left_shift3(z as u32)
+    (left_shift3(z as u32) << 2) | (left_shift3(y as u32) << 1) | left_shift3(x as u32)
 }
 
 pub fn round_up_pow2(mut n: i32) -> i32 {
@@ -1505,7 +1593,7 @@ where
                 }
 
                 let det: T = (0..N).map(|i| lum[i][i]).product();
-                if parity < 0 { -det } else { det }
+                if parity % 2 == 1 { -det } else { det }
             }
         }
     }
@@ -1638,7 +1726,9 @@ pub fn f16_to_f32_software(h: u16) -> f32 {
                 m <<= 1;
                 e += 1;
             }
-            (sign << 31) | ((112 - e) << 23) | ((m & 0x3FF) << 13)
+            // half subnormals have exponent -14; after `e` normalising shifts the
+            // f32 biased exponent is (-14 - e) + 127 = 113 - e.
+            (sign << 31) | ((113 - e) << 23) | ((m & 0x3FF) << 13)
         }
     } else if exp == 0x1F {
         (sign << 31) | (0xFF << 23) | (mant << 13)

shared/src/utils/sampling.rs

@@ -52,7 +52,7 @@ pub fn invert_linear_sample(x: Float, a: Float, b: Float) -> Float {
 
 pub fn sample_bilinear(u: Point2f, w: &[Float]) -> Point2f {
     let y = sample_linear(u[1], w[0] + w[1], w[2] + w[3]);
-    let x = sample_linear(u[1], lerp(y, w[0], w[2]), lerp(y, w[1], w[3]));
+    let x = sample_linear(u[0], lerp(y, w[0], w[2]), lerp(y, w[1], w[3]));
     Point2f::new(x, y)
 }
 
@@ -333,7 +333,8 @@ pub fn invert_spherical_rectangle_sample(
 }
 
 pub fn sample_exponential(u: Float, a: Float) -> Float {
-    (1. - u).ln() / a
+    debug_assert!(a > 0.);
+    -(1. - u).ln() / a
 }
 
 pub fn sample_spherical_triangle(
@@ -550,7 +551,9 @@ pub fn sample_catmull_rom_2d(
         let mut v = 0.;
         for i in 0..4 {
             if weights[i] != 0. {
-                let ind = (offset + i as u32) * n2 + idx;
+                // `offset` is pbrt's `idx - 1` and may be -1; the matching weight is
+                // then zero, so the slot is never actually read.
+                let ind = (offset + i as i32) * n2 as i32 + idx as i32;
                 v += array[ind as usize] * weights[i];
             }
         }
@@ -658,10 +661,10 @@ pub struct VarianceEstimator {
 impl VarianceEstimator {
     pub fn add(&mut self, x: Float) {
         self.n += 1;
-        let delta = x / self.mean;
+        let delta = x - self.mean;
         self.mean += delta / self.n as Float;
         let delta2 = x - self.mean;
-        self.s *= delta * delta2;
+        self.s += delta * delta2;
     }
 
     pub fn mean(&self) -> Float {
@@ -688,7 +691,7 @@ impl VarianceEstimator {
         if ve.n != 0 {
             self.s += ve.s
                 + square(ve.mean - self.mean) * (self.n * ve.n) as Float / (self.n + ve.n) as Float;
-            self.mean +=
+            self.mean =
                 (self.n as Float * self.mean + ve.n as Float * ve.mean) / (self.n + ve.n) as Float;
             self.n += ve.n;
         }
@@ -719,24 +722,37 @@ impl PiecewiseConstant1D {
     }
 
     pub fn new_with_bounds(f: &[Float], min: Float, max: Float) -> Self {
+        debug_assert!(max > min);
         let n = f.len();
+
+        // pbrt takes the absolute value of the function before integrating, so filters
+        // with negative lobes (Mitchell, LanczosSinc) still give a usable CDF.
+        let mut func: GVec<Float> = gvec_with_capacity(n);
+        for &v in f {
+            func.push(v.abs());
+        }
+
         let mut cdf = gvec_with_capacity(n + 1);
         cdf.push(0.0);
 
         let delta = (max - min) / n as Float;
         for i in 0..n {
-            cdf.push(cdf[i] + f[i] * delta);
+            cdf.push(cdf[i] + func[i] * delta);
         }
 
         let func_integral = cdf[n];
-        if func_integral > 0.0 {
-            for c in &mut cdf {
+        if func_integral == 0.0 {
+            for (i, c) in cdf.iter_mut().enumerate().skip(1) {
+                *c = i as Float / n as Float;
+            }
+        } else {
+            for c in cdf.iter_mut().skip(1) {
                 *c /= func_integral;
             }
         }
 
         Self {
-            func: gvec_from_slice(f),
+            func,
             cdf,
             min,
             max,
@@ -748,7 +764,6 @@ impl PiecewiseConstant1D {
     where
         F: Fn(Float) -> Float,
     {
-        let delta = (max - min) / n as Float;
         let delta = (max - min) / n as Float;
         let mut values = gvec_with_capacity(n);
         for i in 0..n {
@@ -816,6 +831,7 @@ impl PiecewiseConstant1D {
 pub struct PiecewiseConstant2D {
     pub conditionals: GVec<PiecewiseConstant1D>,
     pub marginal: PiecewiseConstant1D,
+    pub domain: Bounds2f,
     pub n_u: u32,
     pub n_v: u32,
 }
@@ -852,6 +868,7 @@ impl PiecewiseConstant2D {
         Self {
             conditionals,
             marginal,
+            domain,
             n_u: n_u.try_into().unwrap(),
             n_v: n_v.try_into().unwrap(),
         }
@@ -876,7 +893,10 @@ impl PiecewiseConstant2D {
         self.marginal.integral()
     }
 
-    pub fn pdf(&self, p: Point2f) -> Float {
+    pub fn pdf(&self, pr: Point2f) -> Float {
+        // pbrt maps the query point into the unit domain first; the distribution's
+        // domain is only [0,1]^2 for image distributions, not e.g. filter samplers.
+        let p = Point2f::from(self.domain.offset(&pr));
         let u_offset = ((p.x() * self.n_u as Float) as usize).min(self.n_u as usize - 1);
         let v_offset = ((p.y() * self.n_v as Float) as usize).min(self.n_v as usize - 1);
         let conditional = unsafe { &*self.conditionals.as_ptr().add(v_offset) };

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

src/core/material.rs

@@ -81,16 +81,15 @@ impl MaterialFactory for Material {
                 SubsurfaceMaterial::create(parameters, normal_map, named_materials, &loc, arena)
             }
             "mix" => MixMaterial::create(parameters, normal_map, named_materials, &loc, arena),
-
             _ => Err(anyhow!("Material type '{}' unknown at {}", name, &loc)),
         }
     }
 }
 
 pub fn default_diffuse_material(arena: &Arena) -> Material {
-    use shared::core::texture::SpectrumTexture;
-    use shared::core::texture::SpectrumConstantTexture;
     use shared::core::spectrum::{ConstantSpectrum, Spectrum};
+    use shared::core::texture::SpectrumConstantTexture;
+    use shared::core::texture::SpectrumTexture;
     use shared::materials::DiffuseMaterial;
     use shared::utils::Ptr;
 

src/core/medium.rs

@@ -1,7 +1,7 @@
 use crate::core::spectrum::spectrum_to_photometric;
 use crate::spectra::SRGB;
 use crate::{Arena, FileLoc, ParameterDictionary};
-use anyhow::{bail, Result};
+use anyhow::{Result, bail};
 use shared::core::geometry::{Bounds3f, Point3f, Point3i};
 use shared::core::medium::{
     GridMedium, HGPhaseFunction, HomogeneousMedium, MajorantGrid, Medium, RGBGridMedium,
@@ -262,7 +262,7 @@ static SUBSURFACE_TABLE: &[MeasuredSS] = &[
     },
 ];
 
-fn get_medium_scattering_properties(name: &str) -> Option<(Spectrum, Spectrum)> {
+pub fn get_medium_scattering_properties(name: &str) -> Option<(Spectrum, Spectrum)> {
     SUBSURFACE_TABLE.iter().find(|m| m.name == name).map(|m| {
         let sigma_a = Spectrum::RGBUnbounded(RGBUnboundedSpectrum::new(&SRGB, m.sigma_a.into()));
         let sigma_s =

src/core/texture.rs

@@ -7,10 +7,11 @@ use shared::Float;
 use shared::core::color::ColorEncoding;
 use shared::core::geometry::Vector3f;
 use shared::core::image::WrapMode;
+use shared::core::spectrum::Spectrum;
 use shared::core::texture::SpectrumType;
 use shared::core::texture::{
-    CylindricalMapping, PlanarMapping, PointTransformMapping, SphericalMapping,
-    TextureEvalContext, TextureMapping2D, TextureMapping3D, UVMapping,
+    CylindricalMapping, PlanarMapping, PointTransformMapping, SphericalMapping, TextureEvalContext,
+    TextureMapping2D, TextureMapping3D, UVMapping,
 };
 use shared::spectra::{SampledSpectrum, SampledWavelengths};
 use shared::textures::{
@@ -95,6 +96,21 @@ pub enum SpectrumTexture {
     DirectionMix(SpectrumDirectionMixTexture),
 }
 
+/// A bare `Spectrum` used as a texture is always a constant texture. Saves
+/// writing `SpectrumTexture::Constant(SpectrumConstantTexture::new(s))` at every
+/// default-value site.
+impl From<Spectrum> for SpectrumTexture {
+    fn from(s: Spectrum) -> Self {
+        SpectrumTexture::Constant(SpectrumConstantTexture::new(s))
+    }
+}
+
+impl From<Float> for FloatTexture {
+    fn from(v: Float) -> Self {
+        FloatTexture::Constant(FloatConstantTexture::new(v))
+    }
+}
+
 pub trait CreateSpectrumTexture {
     fn create(
         render_from_texture: Transform,
@@ -115,26 +131,53 @@ impl SpectrumTexture {
         arena: &Arena,
     ) -> Result<Self> {
         match name {
-            "constant" => {
-                SpectrumConstantTexture::create(render_from_texture, params, spectrum_type, loc, arena)
-            }
-            "scale" => {
-                SpectrumScaledTexture::create(render_from_texture, params, spectrum_type, loc, arena)
-            }
-            "mix" => SpectrumMixTexture::create(render_from_texture, params, spectrum_type, loc, arena),
-            "directionmix" => {
-                SpectrumDirectionMixTexture::create(render_from_texture, params, spectrum_type, loc, arena)
-            }
-            "bilerp" => {
-                SpectrumBilerpTexture::create(render_from_texture, params, spectrum_type, loc, arena)
+            "constant" => SpectrumConstantTexture::create(
+                render_from_texture,
+                params,
+                spectrum_type,
+                loc,
+                arena,
+            ),
+            "scale" => SpectrumScaledTexture::create(
+                render_from_texture,
+                params,
+                spectrum_type,
+                loc,
+                arena,
+            ),
+            "mix" => {
+                SpectrumMixTexture::create(render_from_texture, params, spectrum_type, loc, arena)
             }
+            "directionmix" => SpectrumDirectionMixTexture::create(
+                render_from_texture,
+                params,
+                spectrum_type,
+                loc,
+                arena,
+            ),
+            "bilerp" => SpectrumBilerpTexture::create(
+                render_from_texture,
+                params,
+                spectrum_type,
+                loc,
+                arena,
+            ),
             "imagemap" => {
                 SpectrumImageTexture::create(render_from_texture, params, spectrum_type, loc, arena)
             }
-            "checkerboard" => {
-                SpectrumCheckerboardTexture::create(render_from_texture, params, spectrum_type, loc, arena)
+            "checkerboard" => SpectrumCheckerboardTexture::create(
+                render_from_texture,
+                params,
+                spectrum_type,
+                loc,
+                arena,
+            ),
+            "dots" => {
+                SpectrumDotsTexture::create(render_from_texture, params, spectrum_type, loc, arena)
+            }
+            "marble" => {
+                MarbleTexture::create(render_from_texture, params, spectrum_type, loc, arena)
             }
-            "dots" => SpectrumDotsTexture::create(render_from_texture, params, spectrum_type, loc, arena),
             _ => Err(anyhow!(
                 "Spectrum texture type '{}' unknown at {}",
                 name,
@@ -209,7 +252,6 @@ impl CreateTextureMapping for TextureMapping3D {
         let mapping = PointTransformMapping::new(render_from_texture.inverse());
         Ok(TextureMapping3D::PointTransform(mapping))
     }
-
 }
 
 pub static TEXTURE_CACHE: OnceLock<Mutex<HashMap<TexInfo, Arc<MIPMap>>>> = OnceLock::new();

src/lights/infinite.rs

@@ -184,7 +184,7 @@ fn create_portal_light(
 
     // Build distribution
     let duv_dw = |p: Point2f| -> Float {
-        let (_, jacobian) = PortalInfiniteLight::render_from_image(portal_frame, p);
+        let (_, jacobian) = PortalInfiniteLight::render_from_image_with(portal_frame, p);
         jacobian
     };
     let d = remapped.get_sampling_distribution(
@@ -247,7 +247,7 @@ fn remap_image_through_portal(
                     (y as Float + 0.5) / height as Float,
                 );
 
-                let (w_world, _) = PortalInfiniteLight::render_from_image(*portal_frame, uv);
+                let (w_world, _) = PortalInfiniteLight::render_from_image_with(*portal_frame, uv);
                 let w_local = render_from_light.apply_inverse_vector(w_world).normalize();
                 let uv_equi = equal_area_sphere_to_square(w_local);
 

src/materials/coated.rs

@@ -4,8 +4,8 @@ use crate::core::texture::SpectrumTexture;
 use crate::globals::get_options;
 use crate::spectra::data::get_named_spectrum;
 use crate::utils::TextureParameterDictionary;
-use crate::{Arena, FileLoc, ArenaUpload};
-use anyhow::{bail, Result};
+use crate::{Arena, ArenaUpload, FileLoc};
+use anyhow::{Result, bail};
 use shared::core::material::Material;
 use shared::core::spectrum::Spectrum;
 use shared::core::texture::SpectrumType;
@@ -24,10 +24,7 @@ impl CreateMaterial for CoatedDiffuseMaterial {
         arena: &Arena,
     ) -> Result<Material> {
         let reflectance = parameters
-            .get_spectrum_texture("reflectance", None, SpectrumType::Albedo)
-            .unwrap_or(Arc::new(SpectrumTexture::Constant(
-                SpectrumConstantTexture::new(Spectrum::Constant(ConstantSpectrum::new(0.5))),
-            )));
+            .get_spectrum_texture("reflectance", 0.5, SpectrumType::Albedo);
 
         let u_roughness =
             parameters.get_float_texture_with_fallback("uroughness", "roughness", 0.5)?;
@@ -36,22 +33,19 @@ impl CreateMaterial for CoatedDiffuseMaterial {
             parameters.get_float_texture_with_fallback("vroughness", "roughness", 0.5)?;
 
         let thickness = parameters.get_float_texture("thickness", 0.01)?;
-        let eta = parameters
-            .get_float_array("eta")?
-            .first()
-            .map(|&v| Spectrum::Constant(ConstantSpectrum::new(v)))
-            .or_else(|| parameters.get_one_spectrum("eta", None, SpectrumType::Unbounded))
-            .unwrap_or_else(|| Spectrum::Constant(ConstantSpectrum::new(1.5)));
+        let eta = if let Some(&v) = parameters.get_float_array("eta")?.first() {
+            Spectrum::from(v)
+        } else {
+            parameters
+                .get_one_spectrum("eta", None, SpectrumType::Unbounded)
+                .unwrap_or_else(|| Spectrum::from(1.5))
+        };
 
         let max_depth = parameters.get_one_int("maxdepth", 10)?;
         let n_samples = parameters.get_one_int("nsamples", 1)?;
         let g = parameters.get_float_texture("g", 0.)?;
         let albedo = parameters
-            .get_spectrum_texture("albedo", None, SpectrumType::Albedo)
-            .unwrap_or_else(|| {
-                let default_spectrum = Spectrum::Constant(ConstantSpectrum::new(0.));
-                SpectrumTexture::Constant(SpectrumConstantTexture::new(default_spectrum)).into()
-            });
+            .get_spectrum_texture("albedo", 0., SpectrumType::Albedo);
         let displacement = parameters.get_float_texture("displacement", 0.)?;
         let remap_roughness = parameters.get_one_bool("remaproughness", true)?;
 
@@ -145,11 +139,7 @@ impl CreateMaterial for CoatedConductorMaterial {
         let n_samples = parameters.get_one_int("nsamples", 1)?;
         let g = parameters.get_float_texture("g", 0.)?;
         let albedo = parameters
-            .get_spectrum_texture("albedo", None, SpectrumType::Albedo)
-            .unwrap_or_else(|| {
-                let spectrum = Spectrum::Constant(ConstantSpectrum::new(0.));
-                SpectrumTexture::Constant(SpectrumConstantTexture::new(spectrum)).into()
-            });
+            .get_spectrum_texture("albedo", 0., SpectrumType::Albedo);
         let displacement = parameters.get_float_texture_or_null("displacement")?;
         let remap_roughness = parameters.get_one_bool("remaproughness", true)?;
 

src/materials/complex.rs

@@ -1,16 +1,29 @@
 use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
+use crate::core::medium::get_medium_scattering_properties;
 use crate::core::texture::SpectrumTexture;
 use crate::spectra::get_colorspace_device;
-use crate::utils::TextureParameterDictionary;
+use crate::utils::{TextureParameterDictionary, resolve_filename};
 use crate::{Arena, ArenaUpload, FileLoc};
-use anyhow::Result;
+use anyhow::{Result, anyhow, bail};
+use rayon::iter::{IndexedParallelIterator, IntoParallelRefMutIterator, ParallelIterator};
+use rayon::prelude::ParallelSliceMut;
 use shared::bxdfs::HairBxDF;
+use shared::core::bssrdf::BSSRDFTable;
+use shared::core::color::RGB;
 use shared::core::material::Material;
+use shared::core::scattering::{
+    fr_dielectric, fresnel_moment1, fresnel_moment2, henyey_greenstein,
+};
 use shared::core::spectrum::Spectrum;
 use shared::core::texture::SpectrumType;
 use shared::materials::complex::*;
+use shared::spectra::{ConstantSpectrum, RGBUnboundedSpectrum};
 use shared::textures::SpectrumConstantTexture;
+use shared::utils::math::{fast_exp, integrate_catmull_rom, safe_sqrt, square};
+use shared::utils::sampling::sample_exponential;
+use shared::{Float, INV_4_PI, PI};
+
 use std::collections::HashMap;
 use std::sync::Arc;
 
@@ -29,57 +42,284 @@ impl CreateMaterial for HairMaterial {
         let eumelanin = parameters.get_float_texture_or_null("eumelanin")?;
         let pheomelanin = parameters.get_float_texture_or_null("pheomelanin")?;
         let has_melanin = eumelanin.is_some() || pheomelanin.is_some();
-
-        let sigma_a = if sigma_a.is_none() && reflectance.is_none() && !has_melanin {
-            let stdcs = get_colorspace_device();
-            let default_rgb = HairBxDF::sigma_a_from_concentration(1.3, 0.0, stdcs);
-            let spectrum = Spectrum::RGBUnbounded(default_rgb);
-            let texture = SpectrumTexture::Constant(SpectrumConstantTexture::new(spectrum));
-            Some(Arc::new(texture))
-        } else {
-            sigma_a
-        };
-
         let eta = parameters.get_float_texture("eta", 1.55)?;
         let beta_m = parameters.get_float_texture("beta_m", 0.3)?;
         let beta_n = parameters.get_float_texture("beta_n", 0.3)?;
         let alpha = parameters.get_float_texture("alpha", 2.)?;
 
+        let stdcs = get_colorspace_device();
+
+        let absorption = if let Some(s) = sigma_a {
+            HairAbsorption::SigmaA(arena.upload(s))
+        } else if let Some(r) = reflectance {
+            HairAbsorption::Color(arena.upload(r))
+        } else if has_melanin {
+            HairAbsorption::Melanin {
+                eumelanin: arena.upload(eumelanin),
+                pheomelanin: arena.upload(pheomelanin),
+            }
+        } else {
+            let default_rgb = HairBxDF::sigma_a_from_concentration(1.3, 0.0, stdcs.srgb);
+            let spectrum = Spectrum::RGBUnbounded(default_rgb);
+            let texture = SpectrumTexture::Constant(SpectrumConstantTexture::new(spectrum));
+            HairAbsorption::SigmaA(arena.upload(Arc::new(texture)))
+        };
+
         let material = HairMaterial::new(
-            arena.upload(sigma_a),
-            arena.upload(reflectance),
-            arena.upload(eumelanin),
-            arena.upload(pheomelanin),
+            absorption,
             arena.upload(eta),
             arena.upload(beta_m),
             arena.upload(beta_n),
             arena.upload(alpha),
+            stdcs.srgb,
         );
 
         Ok(Material::Hair(material))
     }
 }
 
+fn cube(x: Float) -> Float {
+    x * x * x
+}
+
+fn beam_diffusion_ms(sigma_s: Float, sigma_a: Float, g: Float, eta: Float, r: Float) -> Float {
+    const N_SAMPLES: usize = 100;
+    let mut e_d = 0.;
+    // Precompute information for dipole integrand
+    // Compute reduced scattering coefficients $\sigmaps, \sigmapt$ and albedo $\rhop$
+    let sigmap_s = sigma_s * (1. - g);
+    let sigmap_t = sigma_a + sigmap_s;
+    let rhop = sigmap_s / sigmap_t;
+
+    // Compute non-classical diffusion coefficient $D_\roman{G}$ using Equation
+    // $(\ref{eq:diffusion-coefficient-grosjean})$
+    let d_g = (2. * sigma_a + sigmap_s) / (3. * sigmap_t * sigmap_t);
+
+    // Compute effective transport coefficient $\sigmatr$ based on $D_\roman{G}$
+    let sigma_tr = safe_sqrt(sigma_a / d_g);
+
+    // Determine linear extrapolation distance $\depthextrapolation$ using Equation
+    // $(\ref{eq:dipole-boundary-condition})$
+    let fm1 = fresnel_moment1(eta);
+    let fm2 = fresnel_moment2(eta);
+    let ze = -2. * d_g * (1. + 3. * fm2) / (1. - 2. * fm1);
+
+    // Determine exitance scale factors using Equations $(\ref{eq:kp-exitance-phi})$ and
+    // $(\ref{eq:kp-exitance-e})$
+    let c_phi = 0.25 * (1. - 2. * fm1);
+    let c_e = 0.5 * (1. - 3. * fm2);
+
+    for i in 0..N_SAMPLES {
+        // Sample real point source depth $\depthreal$
+        let zr = sample_exponential((i as Float + 0.5) / N_SAMPLES as Float, sigmap_t);
+
+        // Evaluate dipole integrand $E_{\roman{d}}$ at $\depthreal$ and add to _Ed_
+        let zv = -zr + 2. * ze;
+        let dr = (square(r) + square(zr)).sqrt();
+        let dv = (square(r) + square(zv)).sqrt();
+        // Compute dipole fluence rate $\dipole(r)$ using Equation
+        // $(\ref{eq:diffusion-dipole})$
+        let phi_d =
+            INV_4_PI / d_g * (fast_exp(-sigma_tr * dr) / dr - fast_exp(-sigma_tr * dv) / dv);
+
+        // Compute dipole vector irradiance $-\N{}\cdot\dipoleE(r)$ using Equation
+        // $(\ref{eq:diffusion-dipole-vector-irradiance-normal})$
+        let e_dn = INV_4_PI
+            * (zr * (1. + sigma_tr * dr) * fast_exp(-sigma_tr * dr) / cube(dr)
+                - zv * (1. + sigma_tr * dv) * fast_exp(-sigma_tr * dv) / cube(dv));
+
+        // Add contribution from dipole for depth $\depthreal$ to _Ed_
+        let e = phi_d * c_phi + e_dn * c_e;
+        let kappa = 1. - fast_exp(-2. * sigmap_t * (dr + zr));
+        e_d += kappa * rhop * rhop * e;
+    }
+    return e_d / N_SAMPLES as Float;
+}
+
+fn beam_diffusion_ss(sigma_s: Float, sigma_a: Float, g: Float, eta: Float, r: Float) -> Float {
+    // Compute material parameters and minimum $t$ below the critical angle
+    let sigma_t = sigma_a + sigma_s;
+    let rho = sigma_s / sigma_t;
+    let t_crit = r * safe_sqrt(square(eta) - 1.);
+
+    let mut ess = 0.0;
+    const N_SAMPLES: usize = 100;
+    for i in 0..N_SAMPLES {
+        // Evaluate single-scattering integrand and add to _Ess_
+        let ti = t_crit + sample_exponential((i as Float + 0.5) / N_SAMPLES as Float, sigma_t);
+        // Determine length $d$ of connecting segment and $\cos\theta_\roman{o}$
+        let d = (square(r) + square(ti)).sqrt();
+        let cos_theta_o = ti / d;
+
+        // Add contribution of single scattering at depth $t$
+        ess += rho * fast_exp(-sigma_t * (d + t_crit)) / square(d)
+            * henyey_greenstein(cos_theta_o, g)
+            * (1. - fr_dielectric(-cos_theta_o, eta))
+            * cos_theta_o.abs();
+    }
+    return ess / N_SAMPLES as Float;
+}
+
+fn compute_beam_diffusion_bssrdf(g: Float, eta: Float, t: &mut BSSRDFTable) {
+    let n_rho = t.rho_samples.len();
+    let n_radius = t.radius_samples.len();
+
+    t.radius_samples[0] = 0.;
+    t.radius_samples[1] = 2.5e-3;
+    for i in 2..n_radius {
+        t.radius_samples[i] = t.radius_samples[i - 1] * 1.2;
+    }
+
+    for i in 0..n_rho {
+        t.rho_samples[i] =
+            (1. - fast_exp(-8. * i as Float / (n_rho - 1) as Float)) / (1. - fast_exp(-8.));
+    }
+
+    let rho_samples = &t.rho_samples;
+    let radius_samples = &t.radius_samples;
+    t.profile
+        .par_chunks_mut(n_radius)
+        .zip(t.profile_cdf.par_chunks_mut(n_radius))
+        .zip(t.rho_eff.par_iter_mut())
+        .enumerate()
+        .for_each(|(i, ((profile, cdf), rho_eff))| {
+            // Compute the diffusion profile for the _i_th albedo sample
+            // Compute scattering profile for chosen albedo $\rho$
+            let rho = rho_samples[i];
+            for j in 0..n_radius {
+                let r = radius_samples[j];
+                profile[j] = 2.
+                    * PI
+                    * r
+                    * (beam_diffusion_ss(rho, 1. - rho, g, eta, r)
+                        + beam_diffusion_ms(rho, 1. - rho, g, eta, r));
+            }
+            *rho_eff = integrate_catmull_rom(radius_samples, profile, cdf);
+        });
+}
+
 impl CreateMaterial for SubsurfaceMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        parameters: &TextureParameterDictionary,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
-        _loc: &FileLoc,
-        _arena: &Arena,
+        loc: &FileLoc,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let mut g = parameters.get_one_float("g", 0.0)?;
+        let name = parameters.get_one_string("name", "")?;
+        let scattering = if !name.is_empty() {
+            let (sig_a, sig_s) = get_medium_scattering_properties(&name)
+                .ok_or_else(|| anyhow!("{loc}: named medium {name} not found"))?;
+            if g != 0. {
+                log::warn!("{loc}: non-zero \"g\" ignored with named scattering coefficients");
+            }
+            g = 0.;
+            let sigma_a = SpectrumTexture::Constant(SpectrumConstantTexture::new(sig_a));
+            let sigma_s = SpectrumTexture::Constant(SpectrumConstantTexture::new(sig_s));
+
+            SubsurfaceScattering::Coefficients {
+                sigma_a: arena.upload(&sigma_a),
+                sigma_s: arena.upload(&sigma_s),
+            }
+        } else {
+            let sigma_a =
+                parameters.get_spectrum_texture_or_null("sigma_a", SpectrumType::Unbounded);
+            let sigma_s =
+                parameters.get_spectrum_texture_or_null("sigma_s", SpectrumType::Unbounded);
+            match (sigma_a, sigma_s) {
+                (Some(a), Some(b)) => SubsurfaceScattering::Coefficients {
+                    sigma_a: arena.upload(a),
+                    sigma_s: arena.upload(b),
+                },
+                (Some(_), None) => bail!("{loc}: provided \"sigma_a\" without \"sigma_s\""),
+                (None, Some(_)) => bail!("{loc}: provided \"sigma_s\" without \"sigma_a\""),
+                (None, None) => match parameters
+                    .get_spectrum_texture_or_null("reflectance", SpectrumType::Albedo)
+                {
+                    Some(r) => {
+                        let one = Spectrum::Constant(ConstantSpectrum::new(1.));
+                        let mfp = parameters
+                            .get_spectrum_texture("mfp", one, SpectrumType::Unbounded);
+                        SubsurfaceScattering::Reflectance {
+                            reflectance: arena.upload(r),
+                            mfp: arena.upload(mfp),
+                        }
+                    }
+                    None => {
+                        let stdcs = get_colorspace_device();
+                        let default_sigma_a =
+                            RGBUnboundedSpectrum::new(&stdcs.srgb, RGB::new(0.0011, 0.0024, 0.014));
+                        let default_sigma_s =
+                            RGBUnboundedSpectrum::new(&stdcs.srgb, RGB::new(2.55, 3.21, 3.77));
+                        let sigma_a = SpectrumTexture::Constant(SpectrumConstantTexture::new(
+                            Spectrum::RGBUnbounded(default_sigma_a),
+                        ));
+                        let sigma_s = SpectrumTexture::Constant(SpectrumConstantTexture::new(
+                            Spectrum::RGBUnbounded(default_sigma_s),
+                        ));
+
+                        SubsurfaceScattering::Coefficients {
+                            sigma_a: arena.upload(&sigma_a),
+                            sigma_s: arena.upload(&sigma_s),
+                        }
+                    }
+                },
+            }
+        };
+
+        let scale = parameters.get_one_float("scale", 1.)?;
+        let eta = parameters.get_one_float("eta", 1.33)?;
+
+        let u_roughness =
+            parameters.get_float_texture_with_fallback("uroughness", "roughness", 0.)?;
+        let v_roughness =
+            parameters.get_float_texture_with_fallback("vroughness", "roughness", 0.)?;
+        let displacement = parameters.get_float_texture_or_null("displacement")?;
+        let remap_roughness = parameters.get_one_bool("remaproughness", true)?;
+        let mut table = BSSRDFTable::new(100, 64);
+        compute_beam_diffusion_bssrdf(g, eta, &mut table);
+
+        let mut ss_material = SubsurfaceMaterial {
+            scattering,
+            displacement: arena.upload(displacement),
+            normal_map: arena.upload(normal_map),
+            scale,
+            u_roughness: arena.upload(u_roughness),
+            v_roughness: arena.upload(v_roughness),
+            eta,
+            remap_roughness,
+            table: arena.alloc(table),
+        };
+
+        Ok(Material::Subsurface(ss_material))
     }
 }
 
+// fn brdf_data_from_filename(filename: String) -> Ptr<MeasuredBxDFData> {
+//     static std::map<std::string, MeasuredBxDFData *> loadedData;
+//     if (loadedData.find(filename) == loadedData.end())
+//         loadedData[filename] = MeasuredBxDFData::Create(filename, alloc);
+//     return loadedData[filename];
+//
+// }
+
 impl CreateMaterial for MeasuredMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        parameters: &TextureParameterDictionary,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
-        _loc: &FileLoc,
-        _arena: &Arena,
+        loc: &FileLoc,
+        arena: &Arena,
     ) -> Result<Material> {
+        // let filename = resolve_filename(parameters.get_one_string("filename", "")?);
+        // let displacement = parameters.get_float_texture_or_null("displacement")?;
+        // let brdf = MeasuredBxDF::brdf_data_from_file(filename);
+        // let mat = MeasuredMaterial {
+        //     displacement: arena.upload(displacement),
+        //     normal_map: arena.upload(normal_map)
+        //     brdf
+        // }
         todo!()
     }
 }

src/materials/dielectric.rs

@@ -1,33 +1,75 @@
-use crate::Arena;
 use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 use crate::utils::{FileLoc, TextureParameterDictionary};
+use crate::{Arena, ArenaUpload};
 use anyhow::Result;
 use shared::core::material::Material;
+use shared::core::spectrum::Spectrum;
+use shared::core::texture::SpectrumType;
 use shared::materials::{DielectricMaterial, ThinDielectricMaterial};
+use shared::spectra::ConstantSpectrum;
 use std::collections::HashMap;
 use std::sync::Arc;
 
 impl CreateMaterial for DielectricMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        parameters: &TextureParameterDictionary,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
-        _loc: &FileLoc,
-        _arena: &Arena,
+        loc: &FileLoc,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let eta = if let Some(&v) = parameters.get_float_array("eta")?.first() {
+            Spectrum::from(v)
+        } else {
+            parameters
+                .get_one_spectrum("eta", None, SpectrumType::Unbounded)
+                .unwrap_or_else(|| Spectrum::from(1.5))
+        };
+
+        let u_roughness =
+            parameters.get_float_texture_with_fallback("uroughness", "roughness", 0.)?;
+        let v_roughness =
+            parameters.get_float_texture_with_fallback("vroughness", "roughness", 0.)?;
+        let displacement = parameters.get_float_texture_or_null("displacement")?;
+        let remap_roughness = parameters.get_one_bool("remaproughness", true)?;
+
+        let mat = DielectricMaterial {
+            normal_map: arena.upload(normal_map),
+            displacement: arena.upload(displacement),
+            u_roughness: arena.upload(u_roughness),
+            v_roughness: arena.upload(v_roughness),
+            eta: arena.alloc(eta),
+            remap_roughness,
+        };
+
+        Ok(Material::Dielectric(mat))
     }
 }
 
 impl CreateMaterial for ThinDielectricMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        parameters: &TextureParameterDictionary,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let eta = if let Some(&v) = parameters.get_float_array("eta")?.first() {
+            Spectrum::from(v)
+        } else {
+            parameters
+                .get_one_spectrum("eta", None, SpectrumType::Unbounded)
+                .unwrap_or_else(|| Spectrum::from(1.5))
+        };
+
+        let displacement = parameters.get_float_texture_or_null("displacement")?;
+        let mat = ThinDielectricMaterial {
+            displacement: arena.upload(displacement),
+            normal_map: arena.upload(normal_map),
+            eta: arena.alloc(eta),
+        };
+
+        Ok(Material::ThinDielectric(mat))
     }
 }

src/materials/diffuse.rs

@@ -1,9 +1,9 @@
+use crate::Arena;
 use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
 use crate::core::texture::SpectrumTexture;
 use crate::utils::upload::ArenaUpload;
 use crate::utils::{FileLoc, TextureParameterDictionary};
-use crate::Arena;
 use anyhow::Result;
 use shared::core::material::Material;
 use shared::core::spectrum::Spectrum;
@@ -22,13 +22,7 @@ impl CreateMaterial for DiffuseMaterial {
         _loc: &FileLoc,
         arena: &Arena,
     ) -> Result<Material> {
-        let reflectance = parameters
-            .get_spectrum_texture("reflectance", None, SpectrumType::Albedo)
-            .unwrap_or_else(|| {
-                Arc::new(SpectrumTexture::Constant(
-                    SpectrumConstantTexture::new(Spectrum::Constant(ConstantSpectrum::new(0.5))),
-                ))
-            });
+        let reflectance = parameters.get_spectrum_texture("reflectance", 0.5, SpectrumType::Albedo);
         let displacement = parameters.get_float_texture_or_null("displacement")?;
 
         let specific = DiffuseMaterial {
@@ -42,12 +36,26 @@ impl CreateMaterial for DiffuseMaterial {
 
 impl CreateMaterial for DiffuseTransmissionMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        parameters: &TextureParameterDictionary,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
         _loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let reflectance =
+            parameters.get_spectrum_texture("reflectance", 0.25, SpectrumType::Albedo);
+        let transmittance =
+            parameters.get_spectrum_texture("transmittance", 0.25, SpectrumType::Albedo);
+
+        let displacement = parameters.get_float_texture_or_null("displacement")?;
+        let scale = parameters.get_one_float("scale", 1.)?;
+        let mat = DiffuseTransmissionMaterial {
+            normal_map: arena.upload(normal_map),
+            displacement: arena.upload(displacement),
+            reflectance: arena.upload(reflectance),
+            transmittance: arena.upload(transmittance),
+            scale,
+        };
+        Ok(Material::DiffuseTransmission(mat))
     }
 }

src/materials/mix.rs

@@ -1,7 +1,9 @@
 use crate::core::image::HostImage;
 use crate::core::material::CreateMaterial;
-use crate::utils::{Arena, FileLoc, TextureParameterDictionary};
+use crate::{Arena, ArenaUpload, FileLoc, utils::TextureParameterDictionary};
 use anyhow::Result;
+use anyhow::anyhow;
+use shared::Ptr;
 use shared::core::material::Material;
 use shared::materials::MixMaterial;
 use std::collections::HashMap;
@@ -9,12 +11,31 @@ use std::sync::Arc;
 
 impl CreateMaterial for MixMaterial {
     fn create(
-        _parameters: &TextureParameterDictionary,
+        parameters: &TextureParameterDictionary,
         _normal_map: Option<Arc<HostImage>>,
-        _named_materials: &HashMap<String, Material>,
-        _loc: &FileLoc,
-        _arena: &Arena,
+        named_materials: &HashMap<String, Material>,
+        loc: &FileLoc,
+        arena: &Arena,
     ) -> Result<Material> {
-        todo!()
+        let names = parameters.get_string_array("materials")?;
+        let [a, b]: [String; 2] = names
+            .try_into()
+            .map_err(|_| anyhow!("{loc}: must provide two values for \"string materials\""))?;
+
+        let mut resolve = |n: &str| -> Result<Ptr<Material>> {
+            let m = named_materials
+                .get(n)
+                .ok_or_else(|| anyhow!("{loc}: {n}: named material not found"))?;
+            Ok(arena.alloc(*m))
+        };
+
+        let materials = [resolve(&a)?, resolve(&b)?];
+        let amount = parameters.get_float_texture("amount", 0.5)?;
+
+        let mix_mat = MixMaterial {
+            materials,
+            amount: arena.upload(amount),
+        };
+        Ok(Material::Mix(mix_mat))
     }
 }

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::{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<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> =
-    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<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/textures/checkerboard.rs

@@ -64,8 +64,7 @@ impl CreateSpectrumTexture for SpectrumCheckerboardTexture {
         let tex = |name: &str, def: Spectrum| {
             arena.upload(
                 parameters
-                    .get_spectrum_texture(name, Some(def), spectrum_type)
-                    .expect("default supplied"),
+                    .get_spectrum_texture(name, def, spectrum_type),
             )
         };
 

src/textures/dots.rs

@@ -47,8 +47,7 @@ impl CreateSpectrumTexture for SpectrumDotsTexture {
 
         let get = |name: &str, def: Spectrum| {
             let t = parameters
-                .get_spectrum_texture(name, Some(def), spectrum_type)
-                .expect("default supplied");
+                .get_spectrum_texture(name, def, spectrum_type);
             arena.upload(t)
         };
         let tex = SpectrumDotsTexture::new(map, get("inside", one), get("outside", zero));

src/textures/marble.rs

@@ -1,19 +1,29 @@
-use crate::core::texture::{CreateSpectrumTexture, SpectrumTexture};
-use crate::utils::{FileLoc, TextureParameterDictionary};
 use crate::Arena;
+use crate::core::texture::{CreateSpectrumTexture, CreateTextureMapping, SpectrumTexture};
+use crate::spectra::get_colorspace_device;
+use crate::utils::{FileLoc, TextureParameterDictionary};
 use anyhow::Result;
 use shared::Transform;
-use shared::core::texture::SpectrumType;
+use shared::core::texture::{SpectrumType, TextureMapping3D};
 use shared::textures::MarbleTexture;
 
 impl CreateSpectrumTexture for MarbleTexture {
     fn create(
-        _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
-        _spectrum_type: SpectrumType,
-        _loc: FileLoc,
+        render_from_texture: Transform,
+        parameters: TextureParameterDictionary,
+        spectrum_type: SpectrumType,
+        loc: FileLoc,
         _arena: &Arena,
     ) -> Result<SpectrumTexture> {
-        todo!()
+        let map = TextureMapping3D::create(&parameters, &render_from_texture, &loc)?;
+        let tex = MarbleTexture::new(
+            map,
+            parameters.get_one_int("octaves", 8)?,
+            parameters.get_one_float("roughness", 0.5)?,
+            parameters.get_one_float("scale", 1.)?,
+            parameters.get_one_float("variation", 0.2)?,
+            get_colorspace_device().srgb,
+        );
+        Ok(SpectrumTexture::Marble(tex))
     }
 }

src/textures/mix.rs

@@ -1,20 +1,20 @@
-use crate::core::texture::{
-    CreateSpectrumTexture, FloatTexture, SpectrumTexture };
-use crate::utils::{FileLoc, TextureParameterDictionary};
 use crate::Arena;
+use crate::core::texture::{CreateSpectrumTexture, FloatTexture, SpectrumTexture};
+use crate::utils::{FileLoc, TextureParameterDictionary};
 use anyhow::Result;
-use shared::core::geometry::{Vector3f, VectorLike};
-use shared::core::texture::{SpectrumType, TextureEvalContext};
-use shared::spectra::{SampledSpectrum, SampledWavelengths};
-use shared::utils::Transform;
 use shared::Float;
+use shared::core::geometry::{Vector3f, VectorLike};
+use shared::core::spectrum::Spectrum;
+use shared::core::texture::{SpectrumType, TextureEvalContext};
+use shared::spectra::{ConstantSpectrum, SampledSpectrum, SampledWavelengths};
+use shared::utils::Transform;
 use std::sync::Arc;
 
 #[derive(Clone, Debug)]
 pub struct FloatMixTexture {
     pub tex1: Arc<FloatTexture>,
     pub tex2: Arc<FloatTexture>,
-    pub amount: Arc<FloatTexture>
+    pub amount: Arc<FloatTexture>,
 }
 
 impl FloatMixTexture {
@@ -44,7 +44,7 @@ impl FloatMixTexture {
 pub struct FloatDirectionMixTexture {
     pub tex1: Arc<FloatTexture>,
     pub tex2: Arc<FloatTexture>,
-    pub dir: Vector3f
+    pub dir: Vector3f,
 }
 
 impl FloatDirectionMixTexture {
@@ -71,18 +71,31 @@ impl FloatDirectionMixTexture {
 pub struct SpectrumMixTexture {
     pub tex1: Arc<SpectrumTexture>,
     pub tex2: Arc<SpectrumTexture>,
-    pub amount: Arc<FloatTexture>
+    pub amount: Arc<FloatTexture>,
 }
 
 impl CreateSpectrumTexture for SpectrumMixTexture {
     fn create(
         _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
-        _spectrum_type: SpectrumType,
-        _loc: FileLoc,
-        _arena: &Arena,
+        parameters: TextureParameterDictionary,
+        spectrum_type: SpectrumType,
+        loc: FileLoc,
+        arena: &Arena,
     ) -> Result<SpectrumTexture> {
-        todo!()
+        let zero = Spectrum::Constant(ConstantSpectrum::new(0.));
+        let one = Spectrum::Constant(ConstantSpectrum::new(1.));
+
+        let tex = |name: &str, def: Spectrum| {
+            parameters
+                .get_spectrum_texture(name, def, spectrum_type)
+        };
+
+        let tex = SpectrumMixTexture {
+            tex1: tex("tex1", zero),
+            tex2: tex("tex2", one),
+            amount: parameters.get_float_texture("amount", 0.5)?,
+        };
+        Ok(SpectrumTexture::Mix(tex))
     }
 }
 
@@ -90,18 +103,31 @@ impl CreateSpectrumTexture for SpectrumMixTexture {
 pub struct SpectrumDirectionMixTexture {
     pub tex1: Arc<SpectrumTexture>,
     pub tex2: Arc<SpectrumTexture>,
-    pub dir: Vector3f
+    pub dir: Vector3f,
 }
 
 impl CreateSpectrumTexture for SpectrumDirectionMixTexture {
     fn create(
-        _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
-        _spectrum_type: SpectrumType,
-        _loc: FileLoc,
+        render_from_texture: Transform,
+        parameters: TextureParameterDictionary,
+        spectrum_type: SpectrumType,
+        loc: FileLoc,
         _arena: &Arena,
     ) -> Result<SpectrumTexture> {
-        todo!()
+        let dir_raw = parameters.get_one_vector3f("dir", Vector3f::new(0., 1., 0.))?;
+        let dir = render_from_texture.apply_to_vector(dir_raw).normalize();
+        let tex = |name: &str, def: Spectrum| {
+            parameters
+                .get_spectrum_texture(name, def, spectrum_type)
+        };
+
+        let zero = Spectrum::Constant(ConstantSpectrum::new(0.));
+        let one = Spectrum::Constant(ConstantSpectrum::new(1.));
+        let tex = SpectrumDirectionMixTexture {
+            tex1: tex("tex1", zero),
+            tex2: tex("tex2", one),
+            dir,
+        };
+        Ok(SpectrumTexture::DirectionMix(tex))
     }
 }
-

src/textures/scaled.rs

@@ -66,9 +66,8 @@ impl CreateSpectrumTexture for SpectrumScaledTexture {
         _arena: &Arena,
     ) -> Result<SpectrumTexture> {
         let one = Spectrum::Constant(ConstantSpectrum::new(1.0));
-        let tex = parameters
-            .get_spectrum_texture("tex", Some(one), spectrum_type)
-            .ok_or_else(|| anyhow::anyhow!("{}: missing \"tex\" parameter", _loc))?;
+        // pbrt defaults "tex" to 1.0 rather than erroring (textures.cpp).
+        let tex = parameters.get_spectrum_texture("tex", one, spectrum_type);
         let scale = parameters.get_float_texture("scale", 1.0)?;
 
         if let FloatTexture::Constant(ref cscale) = *scale {

src/textures/windy.rs

@@ -1,20 +1,21 @@
 use crate::Arena;
 use anyhow::Result;
-use shared::{textures::WindyTexture, utils::Transform};
+use shared::{core::texture::TextureMapping3D, textures::WindyTexture, utils::Transform};
 
 use crate::{
-    core::texture::{CreateFloatTexture, FloatTexture },
-    utils::{FileLoc, TextureParameterDictionary}
+    core::texture::{CreateFloatTexture, CreateTextureMapping, FloatTexture},
+    utils::{FileLoc, TextureParameterDictionary},
 };
 
 impl CreateFloatTexture for WindyTexture {
     fn create(
-        _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
-        _loc: FileLoc,
+        render_from_texture: Transform,
+        parameters: TextureParameterDictionary,
+        loc: FileLoc,
         _arena: &Arena,
     ) -> Result<FloatTexture> {
-        todo!()
+        let map = TextureMapping3D::create(&parameters, &render_from_texture, &loc)?;
+        let tex = WindyTexture::new(map);
+        Ok(FloatTexture::Windy(tex))
     }
 }
-

src/textures/wrinkled.rs

@@ -1,20 +1,25 @@
 use crate::Arena;
 use anyhow::Result;
-use shared::{textures::WrinkledTexture, utils::Transform};
+use shared::{core::texture::TextureMapping3D, textures::WrinkledTexture, utils::Transform};
 
 use crate::{
-    core::texture::{CreateFloatTexture, FloatTexture },
-    utils::{FileLoc, TextureParameterDictionary}
+    core::texture::{CreateFloatTexture, CreateTextureMapping, FloatTexture},
+    utils::{FileLoc, TextureParameterDictionary},
 };
 
 impl CreateFloatTexture for WrinkledTexture {
     fn create(
-        _render_from_texture: Transform,
-        _parameters: TextureParameterDictionary,
-        _loc: FileLoc,
-        _arena: &Arena,
+        render_from_texture: Transform,
+        parameters: TextureParameterDictionary,
+        loc: FileLoc,
+        arena: &Arena,
     ) -> Result<FloatTexture> {
-        todo!()
+        let map = TextureMapping3D::create(&parameters, &render_from_texture, &loc)?;
+        let tex = WrinkledTexture::new(
+            map,
+            parameters.get_one_int("octaves", 8)?.try_into().unwrap(),
+            parameters.get_one_float("roughness", 0.5)?,
+        );
+        Ok(FloatTexture::Wrinkled(tex))
     }
 }
-

src/utils/parameters.rs

@@ -835,22 +835,18 @@ impl TextureParameterDictionary {
         self.dict.report_unused()
     }
 
+    /// Fetch a spectrum texture, falling back to a constant built from `def`.
+    /// A default is always supplied, so this cannot fail -- mirrors
+    /// `get_float_texture`. Use `get_spectrum_texture_or_null` when absence is
+    /// meaningful.
     pub fn get_spectrum_texture(
         &self,
         name: &str,
-        val: Option<Spectrum>,
+        def: impl Into<Spectrum>,
         stype: SpectrumType,
-    ) -> Option<Arc<SpectrumTexture>> {
-        let tex = self.get_spectrum_texture_or_null(name, stype);
-        if tex.is_some() {
-            tex
-        } else if val.is_some() {
-            Some(Arc::new(SpectrumTexture::Constant(
-                SpectrumConstantTexture::new(val.unwrap()),
-            )))
-        } else {
-            None
-        }
+    ) -> Arc<SpectrumTexture> {
+        self.get_spectrum_texture_or_null(name, stype)
+            .unwrap_or_else(|| Arc::new(SpectrumTexture::from(def.into())))
     }
 
     pub fn get_float_texture(&self, name: &str, val: Float) -> Result<Arc<FloatTexture>> {

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;