Restructured image code, refactored spectrum code, removed some less idiomatic code to make compiler happy, added mipmap methods

Cargo.toml

@@ -10,7 +10,7 @@ bumpalo = "3.19.0"
 enum_dispatch = "0.3.13"
 exr = "1.73.0"
 half = "2.7.1"
-image = "0.25.8"
+image_rs = { package = "image", version = "0.25.8" }
 num = "0.4.3"
 num-integer = "0.1.46"
 num-traits = "0.2.19"
@@ -24,3 +24,9 @@ thiserror = "2.0.17"
 [features]
 default = []
 use_f64 = []
+
+[lints.clippy]
+excessive_precision = "allow"
+approx_constant = "allow"
+upper_case_acronyms = "allow"
+wrong_self_convention = "allow"

src/camera/mod.rs

@@ -13,7 +13,7 @@ use crate::core::medium::Medium;
 use crate::core::pbrt::{Float, RenderingCoordinateSystem, lerp};
 use crate::core::sampler::CameraSample;
 use crate::geometry::{Normal3f, Point3f, Ray, RayDifferential, Vector3f, VectorLike};
-use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::transform::{AnimatedTransform, Transform};
 
 use std::sync::Arc;
@@ -108,11 +108,7 @@ pub trait CameraTrait {
         sample: CameraSample,
         lambda: &SampledWavelengths,
     ) -> Option<CameraRay> {
-        let mut central_cam_ray = match self.generate_ray(sample, lambda) {
+        let mut central_cam_ray = self.generate_ray(sample, lambda)?;
-            Some(cr) => cr,
-            None => return None,
-        };
-
         let mut rd = RayDifferential::default();
         let mut rx_found = false;
         let mut ry_found = false;

src/camera/orthographic.rs

@@ -5,8 +5,8 @@ use crate::core::sampler::CameraSample;
 use crate::geometry::{
     Bounds2f, Point2f, Point3f, Ray, RayDifferential, Vector2f, Vector3f, VectorLike,
 };
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::sampling::sample_uniform_disk_concentric;
-use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
 use crate::utils::transform::Transform;
 
 pub struct OrthographicCamera {
@@ -115,11 +115,7 @@ impl CameraTrait for OrthographicCamera {
         sample: CameraSample,
         lambda: &SampledWavelengths,
     ) -> Option<CameraRay> {
-        let mut central_cam_ray = match self.generate_ray(sample, lambda) {
+        let mut central_cam_ray = self.generate_ray(sample, lambda)?;
-            Some(cr) => cr,
-            None => return None,
-        };
-
         let mut rd = RayDifferential::default();
         if self.lens_radius > 0.0 {
             return self.generate_ray_differential(sample, lambda);

src/camera/perspective.rs

@@ -7,8 +7,8 @@ use crate::core::sampler::CameraSample;
 use crate::geometry::{
     Bounds2f, Point2f, Point3f, Ray, RayDifferential, Vector2f, Vector3f, VectorLike,
 };
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::sampling::sample_uniform_disk_concentric;
-use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
 use crate::utils::transform::Transform;
 
 pub struct PerspectiveCamera {
@@ -48,7 +48,7 @@ impl PerspectiveCamera {
         );
         let raster_from_screen = raster_from_ndc * ndc_from_screen;
         let screen_from_raster = raster_from_screen.inverse();
-        let camera_from_raster = screen_from_camera.inverse() * screen_from_raster.clone();
+        let camera_from_raster = screen_from_camera.inverse() * screen_from_raster;
         let dx_camera = camera_from_raster.apply_to_point(Point3f::new(1., 0., 0.))
             - camera_from_raster.apply_to_point(Point3f::new(0., 0., 0.));
         let dy_camera = camera_from_raster.apply_to_point(Point3f::new(0., 1., 0.))
@@ -60,7 +60,7 @@ impl PerspectiveCamera {
         let cos_total_width = w_corner_camera.z();
         Self {
             base,
-            screen_from_camera: screen_from_camera.clone(),
+            screen_from_camera: *screen_from_camera,
             camera_from_raster,
             raster_from_screen,
             screen_from_raster,

src/camera/realistic.rs

@@ -5,10 +5,10 @@ use crate::core::sampler::CameraSample;
 use crate::geometry::{
     Bounds2f, Normal3f, Point2f, Point2i, Point3f, Ray, Vector2i, Vector3f, VectorLike,
 };
-use crate::utils::image::Image;
+use crate::image::Image;
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::math::{quadratic, square};
 use crate::utils::scattering::refract;
-use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
 
 pub struct RealisticCamera {
     base: CameraBase,
@@ -61,14 +61,15 @@ impl RealisticCamera {
             element_interface.push(el_int);
         }
 
-        let mut exit_pupil_bounds: Vec<Bounds2f> = Vec::new();
-
         let n_samples = 64;
-        for i in 0..64 {
+        let half_diag = base.film.diagonal() / 2.0;
-            let r0 = i as Float / 64. * base.film.diagonal() / 2.;
+        let exit_pupil_bounds: Vec<_> = (0..n_samples)
-            let r1 = (i + 1) as Float / n_samples as Float * base.film.diagonal() / 2.;
+            .map(|i| {
-            exit_pupil_bounds[i] = self.bound_exit_pupil(r0, r1);
+                let r0 = (i as Float / n_samples as Float) * half_diag;
-        }
+                let r1 = ((i + 1) as Float / n_samples as Float) * half_diag;
+                self.bound_exit_pupil(r0, r1)
+            })
+            .collect();
 
         Self {
             base,

src/camera/spherical.rs

@@ -3,8 +3,8 @@ use crate::core::film::FilmTrait;
 use crate::core::pbrt::{Float, PI};
 use crate::core::sampler::CameraSample;
 use crate::geometry::{Bounds2f, Point2f, Point3f, Ray, Vector3f, spherical_direction};
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::math::{equal_area_square_to_sphere, wrap_equal_area_square};
-use crate::utils::spectrum::{SampledSpectrum, SampledWavelengths};
 
 #[derive(PartialEq)]
 pub struct EquiRectangularMapping;

src/core/aggregates.rs

@@ -61,6 +61,7 @@ impl BVHPrimitiveInfo {
     }
 }
 
+#[derive(Clone, Debug)]
 pub enum BVHBuildNode {
     Leaf {
         first_prim_offset: usize,
@@ -349,7 +350,12 @@ impl BVHAggregate {
             })
             .collect();
 
-        Self::build_upper_sah(treelet_roots, total_nodes)
+        let mut contiguous_nodes: Vec<BVHBuildNode> = treelet_roots
+            .into_iter()
+            .map(|node_box| *node_box)
+            .collect();
+
+        Self::build_upper_sah(&mut contiguous_nodes, total_nodes)
     }
 
     fn emit_lbvh(
@@ -384,7 +390,7 @@ impl BVHAggregate {
         }
         let mask = 1 << bit_index;
         let first_code = morton_prims[0].morton_code;
-        let last_match_index = find_interval(n_primitives as usize, |index| {
+        let last_match_index = find_interval(n_primitives, |index| {
             let current_code = morton_prims[index].morton_code;
             (current_code & mask) == (first_code & mask)
         });
@@ -430,13 +436,10 @@ impl BVHAggregate {
         Box::new(BVHBuildNode::new_interior(axis, child0, child1))
     }
 
-    fn build_upper_sah(
+    fn build_upper_sah(nodes: &mut [BVHBuildNode], total_nodes: &AtomicUsize) -> Box<BVHBuildNode> {
-        nodes: Vec<Box<BVHBuildNode>>,
-        total_nodes: &AtomicUsize,
-    ) -> Box<BVHBuildNode> {
         let n_nodes = nodes.len();
         if n_nodes == 1 {
-            return nodes.into_iter().next().unwrap();
+            return Box::new(nodes[0].clone());
         }
 
         total_nodes.fetch_add(1, AtomicOrdering::Relaxed);
@@ -450,21 +453,16 @@ impl BVHAggregate {
 
         let dim = centroid_bounds.max_dimension();
         if centroid_bounds.p_max[dim] == centroid_bounds.p_min[dim] {
-            // Fallback: Split equally
             let mid = n_nodes / 2;
-            // We have to sort to ensure deterministic split if we are just splitting by index,
+            let (left_part, right_part) = nodes.split_at_mut(mid);
-            // though without a spatial dimension difference, ordering doesn't matter much.
-            // Let's just split the vector.
-            let mut nodes = nodes;
-            let right_nodes = nodes.split_off(mid);
-            let left_nodes = nodes;
 
             return Box::new(BVHBuildNode::new_interior(
                 dim as u8,
-                Self::build_upper_sah(left_nodes, total_nodes),
+                Self::build_upper_sah(left_part, total_nodes),
-                Self::build_upper_sah(right_nodes, total_nodes),
+                Self::build_upper_sah(right_part, total_nodes),
             ));
         }
+
         const N_BUCKETS: usize = 12;
         #[derive(Copy, Clone, Default)]
         struct Bucket {
@@ -472,89 +470,105 @@ impl BVHAggregate {
             bounds: Bounds3f,
         }
         let mut buckets = [Bucket::default(); N_BUCKETS];
-        // Initialize _Bucket_ for HLBVH SAH partition buckets
+        let get_bucket_idx = |node: &BVHBuildNode| -> usize {
-        for node in &nodes {
             let offset = centroid_bounds.offset(&node.bounds().centroid())[dim];
-
             let mut b = (N_BUCKETS as Float * offset) as usize;
             if b == N_BUCKETS {
                 b = N_BUCKETS - 1;
             }
+            b
+        };
+
+        // Initialize _Bucket_ for HLBVH SAH partition buckets
+        for node in nodes.iter() {
+            let b = get_bucket_idx(node);
             buckets[b].count += 1;
             buckets[b].bounds = buckets[b].bounds.union(node.bounds());
         }
 
         // Compute costs for splitting after each bucket
         let mut cost = [0.0; N_BUCKETS - 1];
-        for i in 0..(N_BUCKETS - 1) {
-            let mut b0 = Bounds3f::default();
-            let mut b1 = Bounds3f::default();
-            let mut count0 = 0;
-            let mut count1 = 0;
 
-            for j in 0..=i {
+        // Forward Pass: Accumulate Left side (0 -> N-1)
-                b0 = b0.union(buckets[j].bounds);
+        let mut left_area = [0.0; N_BUCKETS];
-                count0 += buckets[j].count;
+        let mut left_count = [0; N_BUCKETS];
-            }
+        let mut b_left = Bounds3f::default();
-            for j in (i + 1)..N_BUCKETS {
+        let mut c_left = 0;
-                b1 = b1.union(buckets[j].bounds);
+
-                count1 += buckets[j].count;
+        for i in 0..N_BUCKETS {
-            }
+            b_left = b_left.union(buckets[i].bounds);
+            c_left += buckets[i].count;
+            left_area[i] = b_left.surface_area();
+            left_count[i] = c_left;
+        }
+
+        // Backward Pass: Accumulate Right side (N-1 -> 0) and compute cost
+        let mut b_right = Bounds3f::default();
+        let mut c_right = 0;
+        let inv_total_sa = 1.0 / bounds.surface_area();
+
+        for i in (0..N_BUCKETS - 1).rev() {
+            b_right = b_right.union(buckets[i + 1].bounds);
+            c_right += buckets[i + 1].count;
+
+            let count_left = left_count[i];
+            let sa_left = left_area[i];
+            let sa_right = b_right.surface_area();
 
             cost[i] = 0.125
-                + (count0 as Float * b0.surface_area() + count1 as Float * b1.surface_area())
+                + (count_left as Float * sa_left + c_right as Float * sa_right) * inv_total_sa;
-                    / bounds.surface_area();
         }
 
         // Find bucket to split at that minimizes SAH metric
         let mut min_cost = cost[0];
         let mut min_cost_split_bucket = 0;
-        for (i, cost) in cost[1..].iter().enumerate() {
+        for (i, &c) in cost.iter().enumerate().skip(1) {
-            if *cost < min_cost {
+            if c < min_cost {
-                min_cost = *cost;
+                min_cost = c;
                 min_cost_split_bucket = i;
             }
         }
 
         // Split nodes and create interior HLBVH SAH node
-        let (left_nodes, right_nodes): (Vec<_>, Vec<_>) = nodes.into_iter().partition(|node| {
+        let mid = {
-            let offset = centroid_bounds.offset(&node.bounds().centroid())[dim];
+            let mut left = 0;
-            let mut b = (N_BUCKETS as Float * offset) as usize;
+            for i in 0..n_nodes {
-            if b == N_BUCKETS {
+                let b = get_bucket_idx(&nodes[i]);
-                b = N_BUCKETS - 1;
+                if b <= min_cost_split_bucket {
+                    nodes.swap(left, i);
+                    left += 1;
+                }
             }
-            b <= min_cost_split_bucket
+            left
-        });
+        };
 
-        if left_nodes.is_empty() || right_nodes.is_empty() {
+        if mid == 0 || mid == n_nodes {
-            // Recombine to split by count
+            let mid = n_nodes / 2;
-            let mut all_nodes = left_nodes;
-            all_nodes.extend(right_nodes);
 
-            let mid = all_nodes.len() / 2;
+            // Partially sort so the median is in the middle and elements are partitioned around it
-
+            nodes.select_nth_unstable_by(mid, |a, b| {
-            // We must perform a partial sort or select to make the split meaningful spatially
-            all_nodes.select_nth_unstable_by(mid, |a, b| {
                 a.bounds().centroid()[dim]
                     .partial_cmp(&b.bounds().centroid()[dim])
                     .unwrap_or(std::cmp::Ordering::Equal)
             });
 
-            let right = all_nodes.split_off(mid);
+            let (left_part, right_part) = nodes.split_at_mut(mid);
-            let left = all_nodes;
 
-            return Box::new(BVHBuildNode::new_interior(
+            Box::new(BVHBuildNode::new_interior(
                 dim as u8,
-                Self::build_upper_sah(left, total_nodes),
+                Self::build_upper_sah(left_part, total_nodes),
-                Self::build_upper_sah(right, total_nodes),
+                Self::build_upper_sah(right_part, total_nodes),
-            ));
+            ))
-        }
+        } else {
+            // Standard SAH Split
+            let (left_part, right_part) = nodes.split_at_mut(mid);
 
-        Box::new(BVHBuildNode::new_interior(
+            Box::new(BVHBuildNode::new_interior(
-            dim as u8,
+                dim as u8,
-            Self::build_upper_sah(left_nodes, total_nodes),
+                Self::build_upper_sah(left_part, total_nodes),
-            Self::build_upper_sah(right_nodes, total_nodes),
+                Self::build_upper_sah(right_part, total_nodes),
-        ))
+            ))
+        }
     }
 
     fn build_recursive(
@@ -596,7 +610,6 @@ impl BVHAggregate {
             ));
         }
 
-        // let mut mid = n_primitives / 2;
         let mut mid: usize;
         match split_method {
             SplitMethod::Middle => {
@@ -660,9 +673,9 @@ impl BVHAggregate {
                     // Find bucket to split at that minimizes SAH metric>
                     let mut min_cost = Float::INFINITY;
                     let mut min_cost_split_bucket = 0;
-                    for i in 0..N_SPLITS {
+                    for (i, &cost) in costs.iter().enumerate().take(N_SPLITS) {
-                        if costs[i] < min_cost {
+                        if cost < min_cost {
-                            min_cost = costs[i];
+                            min_cost = cost;
                             min_cost_split_bucket = i;
                         }
                     }
@@ -774,7 +787,7 @@ impl BVHAggregate {
                 if node.n_primitives > 0 {
                     // Intersect ray with all primitives in this leaf
                     for i in 0..node.n_primitives {
-                        let prim_idx = node.primitives_offset as usize + i as usize;
+                        let prim_idx = node.primitives_offset + i as usize;
                         let prim = &self.primitives[prim_idx];
 
                         if let Some(si) = prim.intersect(r, Some(hit_t)) {
@@ -800,14 +813,14 @@ impl BVHAggregate {
                         to_visit_offset += 1;
 
                         // Visit Near immediately
-                        current_node_index = node.primitives_offset as usize;
+                        current_node_index = node.primitives_offset;
                     } else {
                         // Ray is positive (Left -> Right).
                         // Push Far
-                        nodes_to_visit[to_visit_offset] = node.primitives_offset as usize;
+                        nodes_to_visit[to_visit_offset] = node.primitives_offset;
                         to_visit_offset += 1;
 
-                        current_node_index = current_node_index + 1;
+                        current_node_index += 1;
                     }
                 }
             } else {
@@ -848,7 +861,7 @@ impl BVHAggregate {
             if node.bounds.intersect_p(r.o, t_max, inv_dir, &dir_is_neg) {
                 if node.n_primitives > 0 {
                     for i in 0..node.n_primitives {
-                        let prim_idx = node.primitives_offset as usize + i as usize;
+                        let prim_idx = node.primitives_offset + i as usize;
                         let prim = &self.primitives[prim_idx];
 
                         if prim.intersect_p(r, Some(t_max)) {
@@ -867,17 +880,13 @@ impl BVHAggregate {
                     // closer to the origin faster, potentially saving work.
 
                     if dir_is_neg[node.axis as usize] == 1 {
-                        // Ray is negative (Right -> Left)
+                        nodes_to_visit[to_visit_offset] = current_node_index + 1;
-                        // Visit Right (Second) first
-                        nodes_to_visit[to_visit_offset] = current_node_index + 1; // Push Left (Far)
                         to_visit_offset += 1;
-                        current_node_index = node.primitives_offset as usize; // Go Right (Near)
+                        current_node_index = node.primitives_offset;
                     } else {
-                        // Ray is positive (Left -> Right)
+                        nodes_to_visit[to_visit_offset] = node.primitives_offset;
-                        // Visit Left (First) first
-                        nodes_to_visit[to_visit_offset] = node.primitives_offset as usize; // Push Right (Far)
                         to_visit_offset += 1;
-                        current_node_index = current_node_index + 1; // Go Left (Near)
+                        current_node_index += 1;
                     }
                 }
             } else {

src/core/bxdf.rs

@@ -7,13 +7,16 @@ use std::sync::{Arc, RwLock};
 
 use enum_dispatch::enum_dispatch;
 
-use crate::core::pbrt::{Float, INV_PI, PI, PI_OVER_2};
+use crate::core::pbrt::{Float, INV_2_PI, INV_PI, PI, PI_OVER_2};
 use crate::geometry::{
     Frame, Normal3f, Point2f, Vector3f, VectorLike, abs_cos_theta, cos_theta, same_hemisphere,
     spherical_direction, spherical_theta,
 };
+use crate::spectra::{
+    N_SPECTRUM_SAMPLES, RGBUnboundedSpectrum, SampledSpectrum, SampledWavelengths,
+};
 use crate::utils::color::RGB;
-use crate::utils::colorspace::RGBColorspace;
+use crate::utils::colorspace::RGBColorSpace;
 use crate::utils::math::{
     fast_exp, i0, log_i0, radians, safe_acos, safe_asin, safe_sqrt, sample_discrete, square,
     trimmed_logistic,
@@ -26,9 +29,6 @@ use crate::utils::scattering::{
     TrowbridgeReitzDistribution, fr_complex, fr_complex_from_spectrum, fr_dielectric, reflect,
     refract,
 };
-use crate::utils::spectrum::{
-    N_SPECTRUM_SAMPLES, RGBUnboundedSpectrum, SampledSpectrum, SampledWavelengths,
-};
 
 bitflags! {
     #[derive(Debug, Clone, Copy, PartialEq, Eq)]
@@ -114,6 +114,19 @@ pub struct BSDFSample {
     pub pdf_is_proportional: bool,
 }
 
+impl Default for BSDFSample {
+    fn default() -> Self {
+        Self {
+            f: SampledSpectrum::default(),
+            wi: Vector3f::default(),
+            pdf: 0.0,
+            flags: BxDFFlags::empty(),
+            eta: 1.0,
+            pdf_is_proportional: false,
+        }
+    }
+}
+
 impl BSDFSample {
     pub fn new(
         f: SampledSpectrum,
@@ -133,17 +146,6 @@ impl BSDFSample {
         }
     }
 
-    pub fn default() -> Self {
-        Self {
-            f: SampledSpectrum::default(),
-            wi: Vector3f::default(),
-            pdf: 0.0,
-            flags: BxDFFlags::empty(),
-            eta: 1.0,
-            pdf_is_proportional: false,
-        }
-    }
-
     #[inline]
     pub fn is_reflective(&self) -> bool {
         self.flags.is_reflective()
@@ -320,7 +322,7 @@ impl HairBxDF {
         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(RGBColorspace::srgb(), sigma_a)
+        RGBUnboundedSpectrum::new(RGBColorSpace::srgb(), sigma_a)
     }
 }
 
@@ -554,11 +556,12 @@ impl BSDF {
         if wo.z() == 0. || !(guard.flags().contains(sampling_flags)) {
             return None;
         }
-        if let Some(mut bs) = guard.sample_f(&wo, u, u2, f_args) {
+        if let Some(mut bs) = guard.sample_f(&wo, u, u2, f_args)
-            if bs.pdf > 0.0 && bs.wi.z() != 0.0 {
+            && bs.pdf > 0.0
-                bs.wi = self.local_to_render(bs.wi);
+            && bs.wi.z() != 0.0
-                return Some(bs);
+        {
-            }
+            bs.wi = self.local_to_render(bs.wi);
+            return Some(bs);
         }
         None
     }
@@ -607,15 +610,16 @@ impl BSDF {
         if wo.z() == 0. || !guard.flags().contains(sampling_flags) {
             return None;
         }
-        if let Some(specific_bxdf) = guard.as_any().downcast_ref::<B>() {
+
-            if let Some(mut bs) = specific_bxdf.sample_f(&wo, u, u2, f_args) {
+        guard
-                if bs.pdf > 0.0 && bs.wi.z() != 0.0 {
+            .as_any()
-                    bs.wi = self.local_to_render(bs.wi);
+            .downcast_ref::<B>()
-                    return Some(bs);
+            .and_then(|specific_bxdf| specific_bxdf.sample_f(&wo, u, u2, f_args))
-                }
+            .filter(|bs| bs.pdf > 0.0 && bs.wi.z() != 0.0)
-            }
+            .map(|mut bs| {
-        }
+                bs.wi = self.local_to_render(bs.wi);
-        None
+                bs
+            })
     }
 
     pub fn pdf_specific<B: BxDFTrait>(
@@ -682,11 +686,13 @@ impl BxDFTrait for DiffuseBxDF {
             wi[2] *= -1.;
         }
         let pdf = cosine_hemisphere_pdf(abs_cos_theta(wi));
-        let mut bsdf = BSDFSample::default();
+        let bsdf = BSDFSample {
-        bsdf.f = self.r * INV_PI;
+            f: self.r * INV_PI,
-        bsdf.wi = wi;
+            wi,
-        bsdf.pdf = pdf;
+            pdf,
-        bsdf.flags = BxDFFlags::DIFFUSE_REFLECTION;
+            flags: BxDFFlags::DIFFUSE_REFLECTION,
+            ..Default::default()
+        };
         Some(bsdf)
     }
 
@@ -707,9 +713,9 @@ impl BxDFTrait for DiffuseBxDF {
 impl BxDFTrait for ConductorBxDF {
     fn flags(&self) -> BxDFFlags {
         if self.mf_distrib.effectively_smooth() {
-            return BxDFFlags::SPECULAR_REFLECTION;
+            BxDFFlags::SPECULAR_REFLECTION
         } else {
-            return BxDFFlags::GLOSSY_REFLECTION;
+            BxDFFlags::GLOSSY_REFLECTION
         }
     }
     fn sample_f(&self, wo: &Vector3f, _uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
@@ -723,12 +729,16 @@ impl BxDFTrait for ConductorBxDF {
             let wi = Vector3f::new(-wo.x(), -wo.y(), wo.z());
             let f =
                 fr_complex_from_spectrum(abs_cos_theta(wi), self.eta, self.k) / abs_cos_theta(wi);
-            let mut bsd = BSDFSample::default();
+
-            bsd.f = f;
+            let bsdf = BSDFSample {
-            bsd.wi = wi;
+                f,
-            bsd.pdf = 1.;
+                wi,
-            bsd.flags = BxDFFlags::SPECULAR_REFLECTION;
+                pdf: 1.,
-            return Some(bsd);
+                flags: BxDFFlags::SPECULAR_REFLECTION,
+                ..Default::default()
+            };
+
+            return Some(bsdf);
         }
 
         if wo.z() == 0. {
@@ -752,11 +762,13 @@ impl BxDFTrait for ConductorBxDF {
         let f = self.mf_distrib.d(wm) * f_spectrum * self.mf_distrib.g(*wo, wi)
             / (4. * cos_theta_i * cos_theta_o);
 
-        let mut bsdf = BSDFSample::default();
+        let bsdf = BSDFSample {
-        bsdf.f = f;
+            f,
-        bsdf.wi = wi;
+            wi,
-        bsdf.pdf = pdf;
+            pdf,
-        bsdf.flags = BxDFFlags::GLOSSY_REFLECTION;
+            flags: BxDFFlags::GLOSSY_REFLECTION,
+            ..Default::default()
+        };
 
         Some(bsdf)
     }
@@ -862,10 +874,10 @@ impl BxDFTrait for DielectricBxDF {
         let fr = fr_dielectric((*wo).dot(wm.into()), self.eta);
 
         if reflect {
-            return SampledSpectrum::new(
+            SampledSpectrum::new(
                 self.mf_distrib.d(wm.into()) * self.mf_distrib.g(*wo, *wi) * fr
                     / (4. * cos_theta_i * cos_theta_o).abs(),
-            );
+            )
         } else {
             let denom = square((*wi).dot(wm.into()) + (*wo).dot(wm.into()) / etap)
                 * cos_theta_i
@@ -878,7 +890,7 @@ impl BxDFTrait for DielectricBxDF {
                 ft /= square(etap)
             }
 
-            return SampledSpectrum::new(ft);
+            SampledSpectrum::new(ft)
         }
     }
 
@@ -930,18 +942,16 @@ impl BxDFTrait for DielectricBxDF {
             return 0.;
         }
 
-        let pdf: Float;
         if reflect {
-            pdf = self.mf_distrib.pdf(
+            self.mf_distrib.pdf(
                 *wo,
                 Vector3f::from(wm) / (4. * (*wo).dot(wm.into()).abs()) * pr / (pt + pr),
-            );
+            )
         } else {
             let denom = square((*wi).dot(wm.into()) + (*wo).dot(wm.into()) / etap);
             let dwm_dwi = (*wi).dot(wm.into()).abs() / denom;
-            pdf = self.mf_distrib.pdf(*wo, wm.into()) * dwm_dwi * pr / (pr + pt);
+            self.mf_distrib.pdf(*wo, wm.into()) * dwm_dwi * pr / (pr + pt)
         }
-        pdf
     }
 
     fn sample_f(&self, wo: &Vector3f, uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
@@ -968,12 +978,14 @@ impl BxDFTrait for DielectricBxDF {
             if uc < pr / (pr + pt) {
                 let wi = Vector3f::new(-wo.x(), -wo.y(), wo.z());
                 let fr = SampledSpectrum::new(r / abs_cos_theta(wi));
-                let mut bsd = BSDFSample::default();
+                let bsdf = BSDFSample {
-                bsd.f = fr;
+                    f: fr,
-                bsd.wi = wi;
+                    wi,
-                bsd.pdf = pr / (pr + pt);
+                    pdf: pr / (pr + pt),
-                bsd.flags = BxDFFlags::SPECULAR_REFLECTION;
+                    flags: BxDFFlags::SPECULAR_REFLECTION,
-                return Some(bsd);
+                    ..Default::default()
+                };
+                Some(bsdf)
             } else {
                 // Compute ray direction for specular transmission
                 if let Some((wi, etap)) = refract(*wo, Normal3f::new(0., 0., 1.), self.eta) {
@@ -981,15 +993,17 @@ impl BxDFTrait for DielectricBxDF {
                     if f_args.mode == TransportMode::Radiance {
                         ft /= square(etap);
                     }
-                    let mut bsd = BSDFSample::default();
+                    let bsdf = BSDFSample {
-                    bsd.f = ft;
+                        f: ft,
-                    bsd.wi = wi;
+                        wi,
-                    bsd.pdf = pt / (pr + pt);
+                        pdf: pt / (pr + pt),
-                    bsd.flags = BxDFFlags::SPECULAR_TRANSMISSION;
+                        flags: BxDFFlags::SPECULAR_TRANSMISSION,
-                    bsd.eta = etap;
+                        eta: etap,
-                    return Some(bsd);
+                        ..Default::default()
+                    };
+                    Some(bsdf)
                 } else {
-                    return None;
+                    None
                 }
             }
         } else {
@@ -1025,17 +1039,19 @@ impl BxDFTrait for DielectricBxDF {
                     self.mf_distrib.d(wm) * self.mf_distrib.g(*wo, wi) * r
                         / (4. * cos_theta(wi) * cos_theta(*wo)),
                 );
-                let mut bsd = BSDFSample::default();
+                let bsdf = BSDFSample {
-                bsd.f = f;
+                    f,
-                bsd.wi = wi;
+                    wi,
-                bsd.pdf = pdf;
+                    pdf,
-                bsd.flags = BxDFFlags::GLOSSY_REFLECTION;
+                    flags: BxDFFlags::GLOSSY_REFLECTION,
-                return Some(bsd);
+                    ..Default::default()
+                };
+                Some(bsdf)
             } else {
                 // Sample transmission at rough dielectric interface
                 if let Some((wi, etap)) = refract(*wo, wm.into(), self.eta) {
                     if same_hemisphere(*wo, wi) || wi.z() == 0. {
-                        return None;
+                        None
                     } else {
                         let denom = square(wi.dot(wm) + wo.dot(wm) / etap);
                         let dwm_mi = wi.dot(wm).abs() / denom;
@@ -1049,16 +1065,18 @@ impl BxDFTrait for DielectricBxDF {
                         if f_args.mode == TransportMode::Radiance {
                             ft /= square(etap);
                         }
-                        let mut bsd = BSDFSample::default();
+                        let bsdf = BSDFSample {
-                        bsd.f = ft;
+                            f: ft,
-                        bsd.wi = wi;
+                            wi,
-                        bsd.pdf = pdf;
+                            pdf,
-                        bsd.flags = BxDFFlags::GLOSSY_TRANSMISSION;
+                            flags: BxDFFlags::GLOSSY_TRANSMISSION,
-                        bsd.eta = etap;
+                            eta: etap,
-                        return Some(bsd);
+                            ..Default::default()
+                        };
+                        Some(bsdf)
                     }
                 } else {
-                    return None;
+                    None
                 }
             }
         }
@@ -1111,23 +1129,27 @@ impl BxDFTrait for ThinDielectricBxDF {
 
         if uc < pr / (pr + pt) {
             let wi = Vector3f::new(-wo.x(), -wo.y(), wo.z());
-            let fr = SampledSpectrum::new(r / abs_cos_theta(wi));
+            let f = SampledSpectrum::new(r / abs_cos_theta(wi));
-            let mut bsdf = BSDFSample::default();
+            let bsdf = BSDFSample {
-            bsdf.f = fr;
+                f,
-            bsdf.wi = wi;
+                wi,
-            bsdf.pdf = pr / (pr + pt);
+                pdf: pr / (pr + pt),
-            bsdf.flags = BxDFFlags::SPECULAR_REFLECTION;
+                flags: BxDFFlags::SPECULAR_REFLECTION,
-            return Some(bsdf);
+                ..Default::default()
+            };
+            Some(bsdf)
         } else {
             // Perfect specular transmission
             let wi = -(*wo);
-            let ft = SampledSpectrum::new(t / abs_cos_theta(wi));
+            let f = SampledSpectrum::new(t / abs_cos_theta(wi));
-            let mut bsdf = BSDFSample::default();
+            let bsdf = BSDFSample {
-            bsdf.f = ft;
+                f,
-            bsdf.wi = wi;
+                wi,
-            bsdf.pdf = pr / (pr + pt);
+                pdf: pr / (pr + pt),
-            bsdf.flags = BxDFFlags::SPECULAR_TRANSMISSION;
+                flags: BxDFFlags::SPECULAR_TRANSMISSION,
-            return Some(bsdf);
+                ..Default::default()
+            };
+            Some(bsdf)
         }
     }
 
@@ -1163,9 +1185,9 @@ impl BxDFTrait for MeasuredBxDF {
         let wm = wm_curr.normalize();
 
         // Map vectors to unit square
-        let theta_o = spherical_theta::<Float>(wo_curr);
+        let theta_o = spherical_theta(wo_curr);
         let phi_o = wo_curr.y().atan2(wo_curr.x());
-        let theta_m = spherical_theta::<Float>(wm);
+        let theta_m = spherical_theta(wm);
         let phi_m = wm.y().atan2(wm.x());
         let u_wo = Point2f::new(MeasuredBxDF::theta2u(theta_o), MeasuredBxDF::phi2u(phi_o));
         let u_wm_phi = if self.brdf.isotropic {
@@ -1181,14 +1203,13 @@ impl BxDFTrait for MeasuredBxDF {
 
         // Inverse parametrization
         let ui = self.brdf.vndf.invert(u_wm, [phi_o, theta_o]);
-        let mut fr = SampledSpectrum::default();
+        let fr = SampledSpectrum::from_fn(|i| {
-        for i in 0..N_SPECTRUM_SAMPLES {
+            self.brdf
-            fr[i] = 0_f32.max(
+                .spectra
-                self.brdf
+                .evaluate(ui.p, [phi_o, theta_o, self.lambda[i]])
-                    .spectra
+                .max(0.0)
-                    .evaluate(ui.p, [phi_o, theta_o, self.lambda[i]]),
+        });
-            );
+
-        }
         fr * self.brdf.ndf.evaluate(u_wm, [])
             / (4. * self.brdf.sigma.evaluate(u_wo, []) * cos_theta(*wi))
     }
@@ -1215,15 +1236,17 @@ impl BxDFTrait for MeasuredBxDF {
             return 0.;
         }
         let wm = wm_curr.normalize();
-        let theta_o = spherical_theta::<Float>(wo_curr);
+        let theta_o = spherical_theta(wo_curr);
         let phi_o = wo_curr.y().atan2(wo_curr.x());
-        let theta_m = spherical_theta::<Float>(wm);
+        let theta_m = spherical_theta(wm);
         let phi_m = wm.y().atan2(wm.x());
+
         let u_wm_phi = if self.brdf.isotropic {
             phi_m - phi_o
         } else {
             phi_m
         };
+
         let mut u_wm = Point2f::new(
             MeasuredBxDF::theta2u(theta_m),
             MeasuredBxDF::phi2u(u_wm_phi),
@@ -1254,7 +1277,7 @@ impl BxDFTrait for MeasuredBxDF {
             flip_w = true;
         }
 
-        let theta_o = spherical_theta::<Float>(wo_curr);
+        let theta_o = spherical_theta(wo_curr);
         let phi_o = wo_curr.y().atan2(wo_curr.x());
         // Warp sample using luminance distribution
         let mut s = self.brdf.luminance.sample(u, [phi_o, theta_o]);
@@ -1281,28 +1304,30 @@ impl BxDFTrait for MeasuredBxDF {
         }
 
         // Interpolate spectral BRDF
-        let mut fr = SampledSpectrum::new(0.);
+        let mut f = SampledSpectrum::from_fn(|i| {
-        for i in 0..N_SPECTRUM_SAMPLES {
+            self.brdf
-            fr[i] = 0_f32.max(
+                .spectra
-                self.brdf
+                .evaluate(u, [phi_o, theta_o, self.lambda[i]])
-                    .spectra
+                .max(0.0)
-                    .evaluate(u, [phi_o, theta_o, self.lambda[i]]),
+        });
-            );
-        }
 
         let u_wo = Point2f::new(MeasuredBxDF::theta2u(theta_o), MeasuredBxDF::phi2u(phi_o));
-        fr *= self.brdf.ndf.evaluate(u_wm, [])
+        f *= self.brdf.ndf.evaluate(u_wm, [])
             / (4. * self.brdf.sigma.evaluate(u_wo, []) * abs_cos_theta(wi));
         pdf /= 4. * wm.dot(wo_curr) * f32::max(2. * square(PI) * u_wm.x(), 1e-6);
 
         if flip_w {
             wi = -wi;
         }
-        let mut bsdf = BSDFSample::default();
+
-        bsdf.f = fr;
+        let bsdf = BSDFSample {
-        bsdf.wi = wi;
+            f,
-        bsdf.pdf = pdf * lum_pdf;
+            wi,
-        bsdf.flags = BxDFFlags::GLOSSY_REFLECTION;
+            pdf: pdf * lum_pdf,
+            flags: BxDFFlags::GLOSSY_REFLECTION,
+            ..Default::default()
+        };
+
         Some(bsdf)
     }
 
@@ -1339,30 +1364,24 @@ impl BxDFTrait for HairBxDF {
         let sampled_value = -self.sigma_a * (2. * cos_gamma_t / cos_theta_t);
         let t = sampled_value.exp();
         let phi = phi_i - phi_o;
-        let ap = Self::ap(cos_theta_o, self.eta, self.h, t);
+        let ap_pdf = Self::ap(cos_theta_o, self.eta, self.h, t);
         let mut f_sum = SampledSpectrum::new(0.);
-        for p in 0..P_MAX {
+        for (p, &ap) in ap_pdf.iter().enumerate().take(P_MAX) {
-            let sin_thetap_o: Float;
+            let (sin_thetap_o, cos_thetap_o) = match p {
-            let cos_thetap_o: Float;
+                0 => (
-            if p == 0 {
+                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
-                sin_thetap_o =
+                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
-                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1];
+                ),
-                cos_thetap_o =
+                1 => (
-                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1];
+                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
-            } else if p == 1 {
+                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
-                sin_thetap_o =
+                ),
-                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0];
+                2 => (
-                cos_thetap_o =
+                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
-                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0];
+                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
-            } else if p == 2 {
+                ),
-                sin_thetap_o =
+                _ => (sin_theta_o, cos_theta_o),
-                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2];
+            };
-                cos_thetap_o =
-                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2];
-            } else {
-                sin_thetap_o = sin_theta_o;
-                cos_thetap_o = cos_theta_o;
-            }
 
             f_sum += Self::mp(
                 cos_theta_i,
@@ -1370,7 +1389,7 @@ impl BxDFTrait for HairBxDF {
                 sin_theta_i,
                 sin_thetap_o,
                 self.v[p],
-            ) * ap[p]
+            ) * ap
                 * Self::np(phi, p as i32, self.s, gamma_o, gamma_t);
         }
         if abs_cos_theta(*wi) > 0. {
@@ -1394,21 +1413,21 @@ impl BxDFTrait for HairBxDF {
         let ap_pdf = self.ap_pdf(cos_theta_o);
         let p =
             sample_discrete(&ap_pdf, uc, None, Some(&mut uc)).expect("Could not sample from AP");
-        let mut sin_thetap_o: Float;
+        let (sin_thetap_o, mut cos_thetap_o) = match p {
-        let mut cos_thetap_o: Float;
+            0 => (
-        if p == 0 {
+                sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
-            sin_thetap_o = sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1];
+                cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
-            cos_thetap_o = cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1];
+            ),
-        } else if p == 1 {
+            1 => (
-            sin_thetap_o = sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0];
+                sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
-            cos_thetap_o = cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0];
+                cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
-        } else if p == 2 {
+            ),
-            sin_thetap_o = sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2];
+            2 => (
-            cos_thetap_o = cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2];
+                sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
-        } else {
+                cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
-            sin_thetap_o = sin_theta_o;
+            ),
-            cos_thetap_o = cos_theta_o;
+            _ => (sin_theta_o, cos_theta_o),
-        }
+        };
 
         cos_thetap_o = cos_thetap_o.abs();
         let cos_theta =
@@ -1434,34 +1453,30 @@ impl BxDFTrait for HairBxDF {
         );
 
         let mut pdf = 0.;
-        for p in 0..P_MAX {
+        for (p, &ap) in ap_pdf.iter().enumerate().take(P_MAX) {
-            if p == 0 {
+            let (sin_thetap_o, cos_thetap_o_raw) = match p {
-                sin_thetap_o =
+                0 => (
-                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1];
+                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
-                cos_thetap_o =
+                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
-                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1];
+                ),
-            } else if p == 1 {
+                1 => (
-                sin_thetap_o =
+                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
-                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0];
+                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
-                cos_thetap_o =
+                ),
-                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0];
+                2 => (
-            } else if p == 2 {
+                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
-                sin_thetap_o =
+                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
-                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2];
+                ),
-                cos_thetap_o =
+                _ => (sin_theta_o, cos_theta_o),
-                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2];
+            };
-            } else {
+            let cos_thetap_o = cos_thetap_o_raw.abs();
-                sin_thetap_o = sin_theta_o;
-                cos_thetap_o = cos_theta_o;
-            }
-            cos_thetap_o = cos_thetap_o.abs();
             pdf += Self::mp(
                 cos_theta_i,
                 cos_thetap_o,
                 sin_theta_i,
                 sin_thetap_o,
                 self.v[p],
-            ) * ap_pdf[p]
+            ) * ap
                 * Self::np(dphi, p as i32, self.s, gamma_o, gamma_t);
         }
         pdf += Self::mp(
@@ -1471,13 +1486,16 @@ impl BxDFTrait for HairBxDF {
             sin_theta_o,
             self.v[P_MAX],
         ) * ap_pdf[P_MAX]
-            * (1. / (2. * PI));
+            * INV_2_PI;
+
+        let bsd = BSDFSample {
+            f: self.f(wo, &wi, f_args.mode),
+            wi,
+            pdf,
+            flags: self.flags(),
+            ..Default::default()
+        };
 
-        let mut bsd = BSDFSample::default();
-        bsd.f = self.f(wo, &wi, f_args.mode);
-        bsd.wi = wi;
-        bsd.pdf = pdf;
-        bsd.flags = self.flags();
         Some(bsd)
     }
 
@@ -1497,38 +1515,34 @@ impl BxDFTrait for HairBxDF {
         // Compute PDF for $A_p$ terms
         let ap_pdf = self.ap_pdf(cos_theta_o);
         let phi = phi_i - phi_o;
-        let mut pdf = 0.;
-        let mut sin_thetap_o: Float;
-        let mut cos_thetap_o: Float;
 
-        for p in 0..P_MAX {
+        let mut pdf = 0.;
-            if p == 0 {
+        for (p, &ap) in ap_pdf.iter().enumerate().take(P_MAX) {
-                sin_thetap_o =
+            let (sin_thetap_o, raw_cos_thetap_o) = match p {
-                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1];
+                0 => (
-                cos_thetap_o =
+                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
-                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1];
+                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
-            } else if p == 1 {
+                ),
-                sin_thetap_o =
+                1 => (
-                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0];
+                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
-                cos_thetap_o =
+                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
-                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0];
+                ),
-            } else if p == 2 {
+                2 => (
-                sin_thetap_o =
+                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
-                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2];
+                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
-                cos_thetap_o =
+                ),
-                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2];
+                _ => (sin_theta_o, cos_theta_o),
-            } else {
+            };
-                sin_thetap_o = sin_theta_o;
+
-                cos_thetap_o = cos_theta_o;
+            let cos_thetap_o = raw_cos_thetap_o.abs();
-            }
+
-            cos_thetap_o = cos_thetap_o.abs();
             pdf += Self::mp(
                 cos_theta_i,
                 cos_thetap_o,
                 sin_theta_i,
                 sin_thetap_o,
                 self.v[p],
-            ) * ap_pdf[p]
+            ) * ap
                 * Self::np(phi, p as i32, self.s, gamma_o, gamma_t);
         }
 
@@ -1539,7 +1553,7 @@ impl BxDFTrait for HairBxDF {
             sin_theta_o,
             self.v[P_MAX],
         ) * ap_pdf[P_MAX]
-            * (1. / (2. * PI));
+            * INV_2_PI;
         pdf
     }
 

src/core/film.rs

@@ -5,19 +5,18 @@ use crate::geometry::{
     Bounds2f, Bounds2fi, Bounds2i, Normal3f, Point2f, Point2i, Point3f, Vector2f, Vector2fi,
     Vector2i, Vector3f,
 };
-use crate::utils::color::{RGB, SRGBEncoding, Triplet, XYZ, white_balance};
+use crate::image::{Image, ImageChannelDesc, ImageChannelValues, ImageMetadata, PixelFormat};
-use crate::utils::colorspace::RGBColorspace;
+use crate::spectra::sampled::{LAMBDA_MAX, LAMBDA_MIN};
-use crate::utils::containers::Array2D;
+use crate::spectra::{
-use crate::utils::image::{
+    ConstantSpectrum, DenselySampledSpectrum, N_SPECTRUM_SAMPLES, SampledSpectrum,
-    Image, ImageChannelDesc, ImageChannelValues, ImageMetadata, PixelFormat,
+    SampledWavelengths, Spectrum, SpectrumTrait, cie_x, cie_y, cie_z,
 };
+use crate::utils::color::{MatrixMulColor, RGB, SRGB, Triplet, XYZ, white_balance};
+use crate::utils::colorspace::RGBColorSpace;
+use crate::utils::containers::Array2D;
 use crate::utils::math::SquareMatrix;
 use crate::utils::math::linear_least_squares;
 use crate::utils::sampling::VarianceEstimator;
-use crate::utils::spectrum::{
-    ConstantSpectrum, DenselySampledSpectrum, LAMBDA_MAX, LAMBDA_MIN, N_SPECTRUM_SAMPLES,
-    SampledSpectrum, SampledWavelengths, Spectrum, inner_product, spectra,
-};
 use crate::utils::transform::AnimatedTransform;
 use rayon::prelude::*;
 use std::sync::{Arc, atomic::AtomicUsize, atomic::Ordering};
@@ -54,7 +53,7 @@ impl Default for RGBPixel {
 impl RGBFilm {
     pub fn new(
         base: FilmBase,
-        colorspace: RGBColorspace,
+        colorspace: RGBColorSpace,
         max_component_value: Float,
         write_fp16: bool,
     ) -> Self {
@@ -82,7 +81,7 @@ pub struct GBufferBFilm {
     output_from_render: AnimatedTransform,
     apply_inverse: bool,
     pixels: Array2D<GBufferPixel>,
-    colorspace: RGBColorspace,
+    colorspace: RGBColorSpace,
     max_component_value: Float,
     write_fp16: bool,
     filter_integral: Float,
@@ -131,13 +130,13 @@ impl PixelSensor {
         r: &Spectrum,
         g: &Spectrum,
         b: &Spectrum,
-        output_colorspace: RGBColorspace,
+        output_colorspace: RGBColorSpace,
         sensor_illum: &Spectrum,
         imaging_ratio: Float,
     ) -> Result<Self, Box<dyn Error>> {
-        let r_bar = DenselySampledSpectrum::from_spectrum(&r, LAMBDA_MIN, LAMBDA_MAX);
+        let r_bar = DenselySampledSpectrum::from_spectrum(r, LAMBDA_MIN, LAMBDA_MAX);
-        let g_bar = DenselySampledSpectrum::from_spectrum(&g, LAMBDA_MIN, LAMBDA_MAX);
+        let g_bar = DenselySampledSpectrum::from_spectrum(g, LAMBDA_MIN, LAMBDA_MAX);
-        let b_bar = DenselySampledSpectrum::from_spectrum(&b, LAMBDA_MIN, LAMBDA_MAX);
+        let b_bar = DenselySampledSpectrum::from_spectrum(b, LAMBDA_MIN, LAMBDA_MAX);
         let mut rgb_camera = [[0.; 3]; N_SWATCH_REFLECTANCES];
 
         for i in 0..N_SWATCH_REFLECTANCES {
@@ -154,16 +153,16 @@ impl PixelSensor {
         }
 
         let mut xyz_output = [[0.; 3]; N_SWATCH_REFLECTANCES];
-        let sensor_white_g = inner_product(sensor_illum, &Spectrum::DenselySampled(g_bar.clone()));
+        let sensor_white_g = sensor_illum.inner_product(&Spectrum::DenselySampled(g_bar.clone()));
-        let sensor_white_y = inner_product(sensor_illum, &spectra::Y);
+        let sensor_white_y = sensor_illum.inner_product(cie_y());
         for i in 0..N_SWATCH_REFLECTANCES {
             let s = SWATCH_REFLECTANCES[i].clone();
             let xyz = Self::project_reflectance::<XYZ>(
                 &s,
                 &output_colorspace.illuminant,
-                &spectra::X,
+                cie_x(),
-                &spectra::Y,
+                cie_y(),
-                &spectra::Z,
+                cie_z(),
             ) * (sensor_white_y / sensor_white_g);
             for c in 0..3 {
                 xyz_output[i][c] = xyz[c];
@@ -182,13 +181,13 @@ impl PixelSensor {
     }
 
     pub fn new_with_white_balance(
-        output_colorspace: RGBColorspace,
+        output_colorspace: RGBColorSpace,
         sensor_illum: Option<Spectrum>,
         imaging_ratio: Float,
     ) -> Self {
-        let r_bar = DenselySampledSpectrum::from_spectrum(&spectra::X, LAMBDA_MIN, LAMBDA_MAX);
+        let r_bar = DenselySampledSpectrum::from_spectrum(cie_x(), LAMBDA_MIN, LAMBDA_MAX);
-        let g_bar = DenselySampledSpectrum::from_spectrum(&spectra::Y, LAMBDA_MIN, LAMBDA_MAX);
+        let g_bar = DenselySampledSpectrum::from_spectrum(cie_y(), LAMBDA_MIN, LAMBDA_MAX);
-        let b_bar = DenselySampledSpectrum::from_spectrum(&spectra::Z, LAMBDA_MIN, LAMBDA_MAX);
+        let b_bar = DenselySampledSpectrum::from_spectrum(cie_z(), LAMBDA_MIN, LAMBDA_MAX);
         let xyz_from_sensor_rgb: SquareMatrix<Float, 3>;
 
         if let Some(illum) = sensor_illum {
@@ -220,13 +219,13 @@ impl PixelSensor {
 
         for lambda_ind in LAMBDA_MIN..=LAMBDA_MAX {
             let lambda = lambda_ind as Float;
-            let illum_val = illum.sample_at(lambda);
+            let illum_val = illum.evaluate(lambda);
 
-            g_integral += b2.sample_at(lambda) * illum_val;
+            g_integral += b2.evaluate(lambda) * illum_val;
-            let refl_illum = refl.sample_at(lambda) * illum_val;
+            let refl_illum = refl.evaluate(lambda) * illum_val;
-            result[0] += b1.sample_at(lambda) * refl_illum;
+            result[0] += b1.evaluate(lambda) * refl_illum;
-            result[1] += b2.sample_at(lambda) * refl_illum;
+            result[1] += b2.evaluate(lambda) * refl_illum;
-            result[2] += b3.sample_at(lambda) * refl_illum;
+            result[2] += b3.evaluate(lambda) * refl_illum;
         }
 
         if g_integral > 0. {
@@ -268,9 +267,10 @@ impl VisibleSurface {
         albedo: SampledSpectrum,
         _lambda: SampledWavelengths,
     ) -> Self {
-        let mut vs = VisibleSurface::default();
+        VisibleSurface {
-        vs.albedo = albedo;
+            albedo,
-        vs
+            ..Default::default()
+        }
     }
 }
 
@@ -392,7 +392,7 @@ pub trait FilmTrait: Sync {
             })
             .collect();
 
-        let mut image = Image::new(format, resolution, channel_names, Arc::new(SRGBEncoding));
+        let mut image = Image::new(format, resolution, channel_names, SRGB);
         let rgb_desc = ImageChannelDesc::new(&[0, 1, 2]);
 
         for (iy, row_data) in processed_rows.into_iter().enumerate() {
@@ -556,7 +556,7 @@ impl FilmTrait for RGBFilm {
                 rgb[c] += rgb_splat[c] as Float / self.filter_integral;
             }
         }
-        self.output_rgbf_from_sensor_rgb * rgb
+        self.output_rgbf_from_sensor_rgb.mul_rgb(rgb)
     }
 
     fn to_output_rgb(&self, l: SampledSpectrum, lambda: &SampledWavelengths) -> RGB {
@@ -564,7 +564,7 @@ impl FilmTrait for RGBFilm {
             .get_pixel_sensor()
             .expect("Sensor must exist")
             .to_sensor_rgb(l, lambda);
-        self.output_rgbf_from_sensor_rgb * sensor_rgb
+        self.output_rgbf_from_sensor_rgb.mul_rgb(sensor_rgb)
     }
 
     fn uses_visible_surface(&self) -> bool {
@@ -627,7 +627,7 @@ impl FilmTrait for GBufferBFilm {
             .get_pixel_sensor()
             .expect("Sensor must exist")
             .to_sensor_rgb(l, lambda);
-        self.output_rgbf_from_sensor_rgb * sensor_rgb
+        self.output_rgbf_from_sensor_rgb.mul_rgb(sensor_rgb)
     }
 
     fn get_pixel_rgb(&self, p: Point2i, splat_scale: Option<Float>) -> RGB {
@@ -652,7 +652,7 @@ impl FilmTrait for GBufferBFilm {
                 rgb[c] += rgb_splat[c] as Float / self.filter_integral;
             }
         }
-        self.output_rgbf_from_sensor_rgb * rgb
+        self.output_rgbf_from_sensor_rgb.mul_rgb(rgb)
     }
 
     fn uses_visible_surface(&self) -> bool {

src/core/interaction.rs

@@ -485,22 +485,23 @@ impl SurfaceInteraction {
         }
         self.shading.dpdu = dpdus;
         self.shading.dpdv = dpdvs;
-        self.shading.dndu = dndus.into();
+        self.shading.dndu = dndus;
-        self.shading.dndv = dndvs.into();
+        self.shading.dndv = dndvs;
     }
 
     pub fn new_simple(pi: Point3fi, n: Normal3f, uv: Point2f) -> Self {
-        let mut si = Self::default();
+        Self {
-        si.common = InteractionData {
+            common: InteractionData {
-            pi,
+                pi,
-            n,
+                n,
-            time: 0.,
+                time: 0.,
-            wo: Vector3f::zero(),
+                wo: Vector3f::zero(),
-            medium_interface: None,
+                medium_interface: None,
-            medium: None,
+                medium: None,
-        };
+            },
-        si.uv = uv;
+            uv,
-        si
+            ..Default::default()
+        }
     }
 
     pub fn set_intersection_properties(
@@ -514,7 +515,7 @@ impl SurfaceInteraction {
         self.area_light = Some(area);
         if prim_medium_interface
             .as_ref()
-            .map_or(false, |mi| mi.is_medium_transition())
+            .is_some_and(|mi| mi.is_medium_transition())
         {
             self.common.medium_interface = prim_medium_interface;
         } else {

src/core/pbrt.rs

@@ -29,22 +29,22 @@ impl FloatBitOps for f32 {
         self.to_bits()
     }
 
+    // Shift 23, Mask 8 bits (0xFF), Bias 127
     #[inline(always)]
     fn exponent_val(self) -> i32 {
         let bits = self.to_bits();
-        // Shift 23, Mask 8 bits (0xFF), Bias 127
         ((bits >> 23) & 0xFF) as i32 - 127
     }
 
+    // Mask bottom 23 bits
     #[inline(always)]
     fn significand_val(self) -> u32 {
-        // Mask bottom 23 bits
         self.to_bits() & ((1 << 23) - 1)
     }
 
+    // Mask top bit (31)
     #[inline(always)]
     fn sign_bit_val(self) -> u32 {
-        // Mask top bit (31)
         self.to_bits() & 0x8000_0000
     }
 
@@ -62,22 +62,22 @@ impl FloatBitOps for f64 {
         self.to_bits()
     }
 
+    // Shift 52, Mask 11 bits (0x7FF), Bias 1023
     #[inline(always)]
     fn exponent_val(self) -> i32 {
         let bits = self.to_bits();
-        // Shift 52, Mask 11 bits (0x7FF), Bias 1023
         ((bits >> 52) & 0x7FF) as i32 - 1023
     }
 
+    // Mask bottom 52 bits
     #[inline(always)]
     fn significand_val(self) -> u64 {
-        // Mask bottom 52 bits
         self.to_bits() & ((1u64 << 52) - 1)
     }
 
+    // Mask top bit (63)
     #[inline(always)]
     fn sign_bit_val(self) -> u64 {
-        // Mask top bit (63)
         self.to_bits() & 0x8000_0000_0000_0000
     }
 
@@ -87,7 +87,7 @@ impl FloatBitOps for f64 {
     }
 }
 
-pub const MACHINE_EPSILON: Float = std::f32::EPSILON * 0.5;
+pub const MACHINE_EPSILON: Float = Float::EPSILON * 0.5;
 pub const SHADOW_EPSILON: Float = 0.0001;
 pub const ONE_MINUS_EPSILON: Float = 0.99999994;
 pub const PI: Float = std::f32::consts::PI;
@@ -99,12 +99,12 @@ pub const PI_OVER_4: Float = 0.785_398_163_397_448_309_61;
 pub const SQRT_2: Float = 1.414_213_562_373_095_048_80;
 
 #[inline]
-pub fn lerp<Value, Factor>(t: Factor, a: Value, b: Value) -> Value
+pub fn lerp<T, F>(t: F, a: T, b: T) -> T
 where
-    Factor: Copy + Num,
+    T: Lerp<F>,
-    Value: Copy + Lerp<Factor>,
+    F: Copy,
 {
-    Value::lerp(t, a, b)
+    T::lerp(t, a, b)
 }
 
 pub fn linear_pdf<T>(x: T, a: T, b: T) -> T
@@ -190,7 +190,7 @@ where
 
 #[inline]
 pub fn gamma(n: i32) -> Float {
-    return (n as Float * MACHINE_EPSILON) / (1. - n as Float * MACHINE_EPSILON);
+    n as Float * MACHINE_EPSILON / (1. - n as Float * MACHINE_EPSILON)
 }
 
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]

src/core/primitive.rs

@@ -187,7 +187,7 @@ pub struct KdTreeAggregate;
 pub enum Primitive {
     Geometric(GeometricPrimitive),
     Transformed(TransformedPrimitive),
-    Animated(AnimatedPrimitive),
+    Animated(Box<AnimatedPrimitive>),
     BVH(BVHAggregatePrimitive),
     KdTree(KdTreeAggregate),
 }

src/core/sampler.rs

@@ -1,5 +1,6 @@
 use std::ops::RangeFull;
 
+use enum_dispatch::enum_dispatch;
 use rand::seq::index::sample;
 
 use crate::core::filter::FilterTrait;
@@ -44,12 +45,12 @@ where
     F: FilterTrait,
 {
     let fs = filter.sample(sampler.get_pixel2d());
-    let mut cs = CameraSample::default();
+    CameraSample {
-    cs.p_film = Point2f::from(p_pixel) + Vector2f::from(fs.p) + Vector2f::new(0.5, 0.5);
+        p_film: Point2f::from(p_pixel) + Vector2f::from(fs.p) + Vector2f::new(0.5, 0.5),
-    cs.time = sampler.get1d();
+        time: sampler.get1d(),
-    cs.p_lens = sampler.get2d();
+        p_lens: sampler.get2d(),
-    cs.filter_weight = fs.weight;
+        filter_weight: fs.weight,
-    cs
+    }
 }
 
 #[derive(Default, Debug, Clone)]
@@ -67,11 +68,13 @@ impl IndependentSampler {
             rng: Rng::default(),
         }
     }
+}
 
-    pub fn samples_per_pixel(&self) -> usize {
+impl SamplerTrait for IndependentSampler {
+    fn samples_per_pixel(&self) -> usize {
         self.samples_per_pixel
     }
-    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
         let hash_input = [p.x() as u64, p.y() as u64, self.seed];
         let sequence_index = hash_buffer(&hash_input, 0);
         self.rng.set_sequence(sequence_index);
@@ -79,13 +82,13 @@ impl IndependentSampler {
             .advance((sample_index as u64) * 65536 + (dim.unwrap_or(0) as u64));
     }
 
-    pub fn get1d(&mut self) -> Float {
+    fn get1d(&mut self) -> Float {
         self.rng.uniform::<Float>()
     }
-    pub fn get2d(&mut self) -> Point2f {
+    fn get2d(&mut self) -> Point2f {
         Point2f::new(self.rng.uniform::<Float>(), self.rng.uniform::<Float>())
     }
-    pub fn get_pixel2d(&mut self) -> Point2f {
+    fn get_pixel2d(&mut self) -> Point2f {
         self.get2d()
     }
 }
@@ -160,82 +163,22 @@ impl HaltonSampler {
             dim: 0,
         }
     }
-    pub fn samples_per_pixel(&self) -> usize {
-        self.samples_per_pixel
-    }
-
-    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
-        self.halton_index = 0;
-
-        let sample_stride = self.base_scales[0] * self.base_scales[1];
-        if sample_stride > 1 {
-            let pm_x = p.x().rem_euclid(MAX_HALTON_RESOLUTION) as u64;
-            let pm_y = p.y().rem_euclid(MAX_HALTON_RESOLUTION) as u64;
-
-            let dim_offset_x = inverse_radical_inverse(pm_x, 2, self.base_exponents[0] as u64);
-
-            self.halton_index = self.halton_index.wrapping_add(
-                dim_offset_x
-                    .wrapping_mul(sample_stride / self.base_scales[0])
-                    .wrapping_mul(self.mult_inverse[0]),
-            );
-
-            let dim_offset_y = inverse_radical_inverse(pm_y, 3, self.base_exponents[1] as u64);
-
-            self.halton_index = self.halton_index.wrapping_add(
-                dim_offset_y
-                    .wrapping_mul(sample_stride / self.base_scales[1])
-                    .wrapping_mul(self.mult_inverse[1]),
-            );
-
-            self.halton_index %= sample_stride;
-        }
-
-        self.halton_index = self
-            .halton_index
-            .wrapping_add((sample_index as u64).wrapping_mul(sample_stride));
-
-        self.dim = 2.max(dim.unwrap_or(0));
-    }
-
-    pub fn get1d(&mut self) -> Float {
-        if self.dim > PRIME_TABLE_SIZE {
-            self.dim = 2;
-        }
-        self.sample_dimension(self.dim)
-    }
-
-    pub fn get2d(&mut self) -> Point2f {
-        if self.dim > PRIME_TABLE_SIZE {
-            self.dim = 2;
-        }
-        let dim = self.dim;
-        self.dim += 2;
-        Point2f::new(self.sample_dimension(dim), self.sample_dimension(dim + 1))
-    }
-
-    pub fn get_pixel2d(&mut self) -> Point2f {
-        Point2f::new(
-            radical_inverse(0, self.halton_index >> self.base_exponents[0]),
-            radical_inverse(1, self.halton_index >> self.base_exponents[1]),
-        )
-    }
 
     fn sample_dimension(&self, dimension: usize) -> Float {
         if self.randomize == RandomizeStrategy::None {
-            return radical_inverse(dimension, self.halton_index);
+            radical_inverse(dimension, self.halton_index)
         } else if self.randomize == RandomizeStrategy::PermuteDigits {
-            return scrambled_radical_inverse(
+            scrambled_radical_inverse(
                 dimension,
                 self.halton_index,
                 &self.digit_permutations[dimension],
-            );
+            )
         } else {
-            return owen_scrambled_radical_inverse(
+            owen_scrambled_radical_inverse(
                 dimension,
                 self.halton_index,
                 mix_bits(1 + ((dimension as u64) << 4)) as u32,
-            );
+            )
         }
     }
 
@@ -257,6 +200,69 @@ impl HaltonSampler {
     }
 }
 
+impl SamplerTrait for HaltonSampler {
+    fn samples_per_pixel(&self) -> usize {
+        self.samples_per_pixel
+    }
+
+    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+        self.halton_index = 0;
+
+        let sample_stride = self.base_scales[0] * self.base_scales[1];
+        if sample_stride > 1 {
+            let pm_x = p.x().rem_euclid(MAX_HALTON_RESOLUTION) as u64;
+            let pm_y = p.y().rem_euclid(MAX_HALTON_RESOLUTION) as u64;
+
+            let dim_offset_x = inverse_radical_inverse(pm_x, 2, self.base_exponents[0]);
+
+            self.halton_index = self.halton_index.wrapping_add(
+                dim_offset_x
+                    .wrapping_mul(sample_stride / self.base_scales[0])
+                    .wrapping_mul(self.mult_inverse[0]),
+            );
+
+            let dim_offset_y = inverse_radical_inverse(pm_y, 3, self.base_exponents[1]);
+
+            self.halton_index = self.halton_index.wrapping_add(
+                dim_offset_y
+                    .wrapping_mul(sample_stride / self.base_scales[1])
+                    .wrapping_mul(self.mult_inverse[1]),
+            );
+
+            self.halton_index %= sample_stride;
+        }
+
+        self.halton_index = self
+            .halton_index
+            .wrapping_add((sample_index as u64).wrapping_mul(sample_stride));
+
+        self.dim = 2.max(dim.unwrap_or(0));
+    }
+
+    fn get1d(&mut self) -> Float {
+        if self.dim > PRIME_TABLE_SIZE {
+            self.dim = 2;
+        }
+        self.sample_dimension(self.dim)
+    }
+
+    fn get2d(&mut self) -> Point2f {
+        if self.dim > PRIME_TABLE_SIZE {
+            self.dim = 2;
+        }
+        let dim = self.dim;
+        self.dim += 2;
+        Point2f::new(self.sample_dimension(dim), self.sample_dimension(dim + 1))
+    }
+
+    fn get_pixel2d(&mut self) -> Point2f {
+        Point2f::new(
+            radical_inverse(0, self.halton_index >> self.base_exponents[0]),
+            radical_inverse(1, self.halton_index >> self.base_exponents[1]),
+        )
+    }
+}
+
 #[derive(Default, Debug, Clone)]
 pub struct StratifiedSampler {
     x_pixel_samples: usize,
@@ -287,12 +293,14 @@ impl StratifiedSampler {
             dim: 0,
         }
     }
+}
 
-    pub fn samples_per_pixel(&self) -> usize {
+impl SamplerTrait for StratifiedSampler {
+    fn samples_per_pixel(&self) -> usize {
         self.x_pixel_samples * self.y_pixel_samples
     }
 
-    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
         self.pixel = p;
         self.sample_index = sample_index;
         let hash_input = [p.x() as u64, p.y() as u64, self.seed];
@@ -302,7 +310,7 @@ impl StratifiedSampler {
             .advance((sample_index as u64) * 65536 + (dim.unwrap_or(0) as u64));
     }
 
-    pub fn get1d(&mut self) -> Float {
+    fn get1d(&mut self) -> Float {
         let hash_input = [
             self.pixel.x() as u64,
             self.pixel.y() as u64,
@@ -325,7 +333,7 @@ impl StratifiedSampler {
         (stratum as Float + delta) / (self.samples_per_pixel() as Float)
     }
 
-    pub fn get2d(&mut self) -> Point2f {
+    fn get2d(&mut self) -> Point2f {
         let hash_input = [
             self.pixel.x() as u64,
             self.pixel.y() as u64,
@@ -357,7 +365,7 @@ impl StratifiedSampler {
         )
     }
 
-    pub fn get_pixel2d(&mut self) -> Point2f {
+    fn get_pixel2d(&mut self) -> Point2f {
         self.get2d()
     }
 }
@@ -383,16 +391,35 @@ impl PaddedSobolSampler {
             dim: 0,
         }
     }
-    pub fn samples_per_pixel(&self) -> usize {
+
+    fn sample_dimension(&self, dimension: usize, a: u32, hash: u32) -> Float {
+        if self.randomize == RandomizeStrategy::None {
+            return sobol_sample(a as u64, dimension, NoRandomizer);
+        }
+        match self.randomize {
+            RandomizeStrategy::PermuteDigits => {
+                sobol_sample(a as u64, dimension, BinaryPermuteScrambler::new(hash))
+            }
+            RandomizeStrategy::FastOwen => {
+                sobol_sample(a as u64, dimension, FastOwenScrambler::new(hash))
+            }
+            RandomizeStrategy::Owen => sobol_sample(a as u64, dimension, OwenScrambler::new(hash)),
+            RandomizeStrategy::None => unreachable!(),
+        }
+    }
+}
+
+impl SamplerTrait for PaddedSobolSampler {
+    fn samples_per_pixel(&self) -> usize {
         self.samples_per_pixel
     }
-    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
         self.pixel = p;
         self.sample_index = sample_index;
         self.dim = dim.unwrap_or(0);
     }
 
-    pub fn get1d(&mut self) -> Float {
+    fn get1d(&mut self) -> Float {
         let hash_input = [
             self.pixel.x() as u64,
             self.pixel.y() as u64,
@@ -407,7 +434,7 @@ impl PaddedSobolSampler {
         );
         self.sample_dimension(0, index, hash >> 32)
     }
-    pub fn get2d(&mut self) -> Point2f {
+    fn get2d(&mut self) -> Point2f {
         let hash_input = [
             self.pixel.x() as u64,
             self.pixel.y() as u64,
@@ -427,25 +454,9 @@ impl PaddedSobolSampler {
         )
     }
 
-    pub fn get_pixel2d(&mut self) -> Point2f {
+    fn get_pixel2d(&mut self) -> Point2f {
         self.get2d()
     }
-
-    fn sample_dimension(&self, dimension: usize, a: u32, hash: u32) -> Float {
-        if self.randomize == RandomizeStrategy::None {
-            return sobol_sample(a as u64, dimension, NoRandomizer);
-        }
-        match self.randomize {
-            RandomizeStrategy::PermuteDigits => {
-                sobol_sample(a as u64, dimension, BinaryPermuteScrambler::new(hash))
-            }
-            RandomizeStrategy::FastOwen => {
-                sobol_sample(a as u64, dimension, FastOwenScrambler::new(hash))
-            }
-            RandomizeStrategy::Owen => sobol_sample(a as u64, dimension, OwenScrambler::new(hash)),
-            RandomizeStrategy::None => unreachable!(),
-        }
-    }
 }
 
 #[derive(Default, Debug, Clone)]
@@ -477,55 +488,6 @@ impl SobolSampler {
             sobol_index: 0,
         }
     }
-    pub fn samples_per_pixel(&self) -> usize {
-        self.samples_per_pixel
-    }
-    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
-        self.pixel = p;
-        self.dim = 2.max(dim.unwrap_or(0));
-        self.sobol_index =
-            sobol_interval_to_index(log2_int(self.scale as Float) as u32, sample_index as u64, p)
-    }
-
-    pub fn get1d(&mut self) -> Float {
-        if self.dim >= N_SOBOL_DIMENSIONS {
-            self.dim = 2;
-        }
-
-        let dim = self.dim;
-        self.dim += 1;
-        self.sample_dimension(dim)
-    }
-
-    pub fn get2d(&mut self) -> Point2f {
-        if self.dim >= N_SOBOL_DIMENSIONS {
-            self.dim = 2;
-        }
-        let u = Point2f::new(
-            self.sample_dimension(self.dim),
-            self.sample_dimension(self.dim + 1),
-        );
-        self.dim += 2;
-        u
-    }
-
-    pub fn get_pixel2d(&mut self) -> Point2f {
-        let mut u = Point2f::new(
-            sobol_sample(self.sobol_index, 0, NoRandomizer),
-            sobol_sample(self.sobol_index, 1, NoRandomizer),
-        );
-        u[0] = clamp_t(
-            u[0] * self.scale as Float - self.pixel[0] as Float,
-            0.,
-            ONE_MINUS_EPSILON,
-        ) as Float;
-        u[1] = clamp_t(
-            u[1] * self.scale as Float - self.pixel[1] as Float,
-            1.,
-            ONE_MINUS_EPSILON,
-        ) as Float;
-        u
-    }
 
     fn sample_dimension(&self, dimension: usize) -> Float {
         if self.randomize == RandomizeStrategy::None {
@@ -550,6 +512,58 @@ impl SobolSampler {
     }
 }
 
+impl SamplerTrait for SobolSampler {
+    fn samples_per_pixel(&self) -> usize {
+        self.samples_per_pixel
+    }
+    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+        self.pixel = p;
+        self.dim = 2.max(dim.unwrap_or(0));
+        self.sobol_index =
+            sobol_interval_to_index(log2_int(self.scale as Float) as u32, sample_index as u64, p)
+    }
+
+    fn get1d(&mut self) -> Float {
+        if self.dim >= N_SOBOL_DIMENSIONS {
+            self.dim = 2;
+        }
+
+        let dim = self.dim;
+        self.dim += 1;
+        self.sample_dimension(dim)
+    }
+
+    fn get2d(&mut self) -> Point2f {
+        if self.dim >= N_SOBOL_DIMENSIONS {
+            self.dim = 2;
+        }
+        let u = Point2f::new(
+            self.sample_dimension(self.dim),
+            self.sample_dimension(self.dim + 1),
+        );
+        self.dim += 2;
+        u
+    }
+
+    fn get_pixel2d(&mut self) -> Point2f {
+        let mut u = Point2f::new(
+            sobol_sample(self.sobol_index, 0, NoRandomizer),
+            sobol_sample(self.sobol_index, 1, NoRandomizer),
+        );
+        u[0] = clamp_t(
+            u[0] * self.scale as Float - self.pixel[0] as Float,
+            0.,
+            ONE_MINUS_EPSILON,
+        ) as Float;
+        u[1] = clamp_t(
+            u[1] * self.scale as Float - self.pixel[1] as Float,
+            1.,
+            ONE_MINUS_EPSILON,
+        ) as Float;
+        u
+    }
+}
+
 #[derive(Default, Debug, Clone)]
 pub struct ZSobolSampler {
     randomize: RandomizeStrategy,
@@ -569,7 +583,7 @@ impl ZSobolSampler {
     ) -> Self {
         let log2_samples_per_pixel = log2_int(samples_per_pixel as Float) as u32;
         let res = round_up_pow2(full_resolution.x().max(full_resolution.y()));
-        let log4_samples_per_pixel = (log2_samples_per_pixel + 1) / 2;
+        let log4_samples_per_pixel = log2_samples_per_pixel.div_ceil(2);
         let n_base4_digits = log2_int(res as Float) as u32 + log4_samples_per_pixel;
         Self {
             randomize,
@@ -581,71 +595,6 @@ impl ZSobolSampler {
         }
     }
 
-    pub fn samples_per_pixel(&self) -> usize {
-        todo!()
-    }
-    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
-        self.dim = dim.unwrap_or(0);
-        self.morton_index = (encode_morton_2(p.x() as u32, p.y() as u32)
-            << self.log2_samples_per_pixel)
-            | sample_index as u64
-    }
-
-    pub fn get1d(&mut self) -> Float {
-        let sample_index = self.get_sample_index();
-        let hash_input = [self.dim as u64, self.seed];
-        let hash = hash_buffer(&hash_input, 0) as u32;
-        self.dim += 1;
-        if self.randomize == RandomizeStrategy::None {
-            return sobol_sample(sample_index, self.dim, NoRandomizer);
-        }
-        match self.randomize {
-            RandomizeStrategy::PermuteDigits => {
-                sobol_sample(sample_index, self.dim, BinaryPermuteScrambler::new(hash))
-            }
-            RandomizeStrategy::FastOwen => {
-                sobol_sample(sample_index, self.dim, FastOwenScrambler::new(hash))
-            }
-            RandomizeStrategy::Owen => {
-                sobol_sample(sample_index, self.dim, OwenScrambler::new(hash))
-            }
-            RandomizeStrategy::None => unreachable!(),
-        }
-    }
-
-    pub fn get2d(&mut self) -> Point2f {
-        let sample_index = self.get_sample_index();
-        self.dim += 2;
-        let hash_input = [self.dim as u64, self.seed];
-        let hash = hash_buffer(&hash_input, 0);
-        let sample_hash = [hash as u32, (hash >> 32) as u32];
-        if self.randomize == RandomizeStrategy::None {
-            return Point2f::new(
-                sobol_sample(sample_index, 0, NoRandomizer),
-                sobol_sample(sample_index, 1, NoRandomizer),
-            );
-        }
-        match self.randomize {
-            RandomizeStrategy::PermuteDigits => Point2f::new(
-                sobol_sample(sample_index, 0, BinaryPermuteScrambler::new(sample_hash[0])),
-                sobol_sample(sample_index, 1, BinaryPermuteScrambler::new(sample_hash[1])),
-            ),
-            RandomizeStrategy::FastOwen => Point2f::new(
-                sobol_sample(sample_index, 0, FastOwenScrambler::new(sample_hash[0])),
-                sobol_sample(sample_index, 1, FastOwenScrambler::new(sample_hash[1])),
-            ),
-            RandomizeStrategy::Owen => Point2f::new(
-                sobol_sample(sample_index, 0, OwenScrambler::new(sample_hash[0])),
-                sobol_sample(sample_index, 1, OwenScrambler::new(sample_hash[1])),
-            ),
-            RandomizeStrategy::None => unreachable!(),
-        }
-    }
-
-    pub fn get_pixel2d(&mut self) -> Point2f {
-        self.get2d()
-    }
-
     fn get_sample_index(&self) -> u64 {
         const PERMUTATIONS: [[u8; 4]; 24] = [
             [0, 1, 2, 3],
@@ -701,35 +650,107 @@ impl ZSobolSampler {
         sample_index
     }
 }
+
+impl SamplerTrait for ZSobolSampler {
+    fn samples_per_pixel(&self) -> usize {
+        todo!()
+    }
+    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+        self.dim = dim.unwrap_or(0);
+        self.morton_index = (encode_morton_2(p.x() as u32, p.y() as u32)
+            << self.log2_samples_per_pixel)
+            | sample_index as u64
+    }
+
+    fn get1d(&mut self) -> Float {
+        let sample_index = self.get_sample_index();
+        let hash_input = [self.dim as u64, self.seed];
+        let hash = hash_buffer(&hash_input, 0) as u32;
+        self.dim += 1;
+        if self.randomize == RandomizeStrategy::None {
+            return sobol_sample(sample_index, self.dim, NoRandomizer);
+        }
+        match self.randomize {
+            RandomizeStrategy::PermuteDigits => {
+                sobol_sample(sample_index, self.dim, BinaryPermuteScrambler::new(hash))
+            }
+            RandomizeStrategy::FastOwen => {
+                sobol_sample(sample_index, self.dim, FastOwenScrambler::new(hash))
+            }
+            RandomizeStrategy::Owen => {
+                sobol_sample(sample_index, self.dim, OwenScrambler::new(hash))
+            }
+            RandomizeStrategy::None => unreachable!(),
+        }
+    }
+
+    fn get2d(&mut self) -> Point2f {
+        let sample_index = self.get_sample_index();
+        self.dim += 2;
+        let hash_input = [self.dim as u64, self.seed];
+        let hash = hash_buffer(&hash_input, 0);
+        let sample_hash = [hash as u32, (hash >> 32) as u32];
+        if self.randomize == RandomizeStrategy::None {
+            return Point2f::new(
+                sobol_sample(sample_index, 0, NoRandomizer),
+                sobol_sample(sample_index, 1, NoRandomizer),
+            );
+        }
+        match self.randomize {
+            RandomizeStrategy::PermuteDigits => Point2f::new(
+                sobol_sample(sample_index, 0, BinaryPermuteScrambler::new(sample_hash[0])),
+                sobol_sample(sample_index, 1, BinaryPermuteScrambler::new(sample_hash[1])),
+            ),
+            RandomizeStrategy::FastOwen => Point2f::new(
+                sobol_sample(sample_index, 0, FastOwenScrambler::new(sample_hash[0])),
+                sobol_sample(sample_index, 1, FastOwenScrambler::new(sample_hash[1])),
+            ),
+            RandomizeStrategy::Owen => Point2f::new(
+                sobol_sample(sample_index, 0, OwenScrambler::new(sample_hash[0])),
+                sobol_sample(sample_index, 1, OwenScrambler::new(sample_hash[1])),
+            ),
+            RandomizeStrategy::None => unreachable!(),
+        }
+    }
+
+    fn get_pixel2d(&mut self) -> Point2f {
+        self.get2d()
+    }
+}
+
 #[derive(Default, Debug, Clone)]
 pub struct MLTSampler;
-impl MLTSampler {
+impl SamplerTrait for MLTSampler {
-    pub fn samples_per_pixel(&self) -> usize {
+    fn samples_per_pixel(&self) -> usize {
         todo!()
     }
-    pub fn start_pixel_sample(&mut self, _p: Point2i, _sample_index: usize, _dim: Option<usize>) {
+    fn start_pixel_sample(&mut self, _p: Point2i, _sample_index: usize, _dim: Option<usize>) {
         todo!()
     }
-    pub fn get1d(&mut self) -> Float {
+    fn get1d(&mut self) -> Float {
         todo!()
     }
-    pub fn get2d(&mut self) -> Point2f {
+    fn get2d(&mut self) -> Point2f {
         todo!()
     }
-    pub fn get_pixel2d(&mut self) -> Point2f {
+    fn get_pixel2d(&mut self) -> Point2f {
         todo!()
     }
 }
 
+#[enum_dispatch]
 pub trait SamplerTrait {
     fn samples_per_pixel(&self) -> usize;
     fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>);
     fn get1d(&mut self) -> Float;
     fn get2d(&mut self) -> Point2f;
     fn get_pixel2d(&mut self) -> Point2f;
-    fn clone_sampler(&self) -> Box<Sampler>;
+    // fn clone_sampler(&self) -> Box<Sampler> {
+    //     Box::new(self.clone())
+    // }
 }
 
+#[enum_dispatch(SamplerTrait)]
 #[derive(Debug, Clone)]
 pub enum Sampler {
     Independent(IndependentSampler),
@@ -740,68 +761,3 @@ pub enum Sampler {
     ZSobol(ZSobolSampler),
     MLT(MLTSampler),
 }
-
-impl SamplerTrait for Sampler {
-    fn samples_per_pixel(&self) -> usize {
-        match self {
-            Self::Independent(inner) => inner.samples_per_pixel(),
-            Self::Stratified(inner) => inner.samples_per_pixel(),
-            Self::Halton(inner) => inner.samples_per_pixel(),
-            Self::PaddedSobol(inner) => inner.samples_per_pixel(),
-            Self::Sobol(inner) => inner.samples_per_pixel(),
-            Self::ZSobol(inner) => inner.samples_per_pixel(),
-            Self::MLT(inner) => inner.samples_per_pixel(),
-        }
-    }
-
-    fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
-        match self {
-            Self::Independent(inner) => inner.start_pixel_sample(p, sample_index, dim),
-            Self::Stratified(inner) => inner.start_pixel_sample(p, sample_index, dim),
-            Self::Halton(inner) => inner.start_pixel_sample(p, sample_index, dim),
-            Self::PaddedSobol(inner) => inner.start_pixel_sample(p, sample_index, dim),
-            Self::Sobol(inner) => inner.start_pixel_sample(p, sample_index, dim),
-            Self::ZSobol(inner) => inner.start_pixel_sample(p, sample_index, dim),
-            Self::MLT(inner) => inner.start_pixel_sample(p, sample_index, dim),
-        }
-    }
-
-    fn get1d(&mut self) -> Float {
-        match self {
-            Self::Independent(inner) => inner.get1d(),
-            Self::Stratified(inner) => inner.get1d(),
-            Self::Halton(inner) => inner.get1d(),
-            Self::PaddedSobol(inner) => inner.get1d(),
-            Self::Sobol(inner) => inner.get1d(),
-            Self::ZSobol(inner) => inner.get1d(),
-            Self::MLT(inner) => inner.get1d(),
-        }
-    }
-
-    fn get2d(&mut self) -> Point2f {
-        match self {
-            Self::Independent(inner) => inner.get2d(),
-            Self::Stratified(inner) => inner.get2d(),
-            Self::Halton(inner) => inner.get2d(),
-            Self::PaddedSobol(inner) => inner.get2d(),
-            Self::Sobol(inner) => inner.get2d(),
-            Self::ZSobol(inner) => inner.get2d(),
-            Self::MLT(inner) => inner.get2d(),
-        }
-    }
-    fn get_pixel2d(&mut self) -> Point2f {
-        match self {
-            Self::Independent(inner) => inner.get_pixel2d(),
-            Self::Stratified(inner) => inner.get_pixel2d(),
-            Self::Halton(inner) => inner.get_pixel2d(),
-            Self::PaddedSobol(inner) => inner.get_pixel2d(),
-            Self::Sobol(inner) => inner.get_pixel2d(),
-            Self::ZSobol(inner) => inner.get_pixel2d(),
-            Self::MLT(inner) => inner.get_pixel2d(),
-        }
-    }
-
-    fn clone_sampler(&self) -> Box<Sampler> {
-        Box::new(self.clone())
-    }
-}

src/core/texture.rs

@@ -1,34 +1,279 @@
-use crate::core::pbrt::Float;
-use crate::geometry::{Point3f, Vector3f, Normal3f, Point2f};
 use crate::core::interaction::SurfaceInteraction;
+use crate::core::pbrt::Float;
+use crate::core::pbrt::{INV_2_PI, INV_PI, PI};
+use crate::geometry::{Normal3f, Point2f, Point3f, Vector3f, VectorLike};
+use crate::geometry::{spherical_phi, spherical_theta};
+use crate::image::WrapMode;
+use crate::spectra::{SampledSpectrum, Spectrum};
+use crate::spectra::{SampledWavelengths, SpectrumTrait};
+use crate::utils::color::ColorEncoding;
+use crate::utils::math::square;
+use crate::utils::mipmap::MIPMap;
+use crate::utils::mipmap::{MIPMapFilterOptions, MIPMapSample};
+use crate::utils::transform::Transform;
+
+use std::collections::HashMap;
+use std::sync::{Arc, Mutex, OnceLock};
+
+use enum_dispatch::enum_dispatch;
+
+struct TexCoord2D {
+    st: Point2f,
+    dsdx: Float,
+    dsdy: Float,
+    dtdx: Float,
+    dtdy: Float,
+}
+
+#[enum_dispatch]
+trait TextureMapping2DTrait {
+    fn map(&self, ctx: TextureEvalContext) -> TexCoord2D;
+}
+
+#[enum_dispatch(TextureMapping2DTrait)]
+pub enum TextureMapping2D {
+    UV(UVMapping),
+    Spherical(SphericalMapping),
+    Cylindrical(CylindricalMapping),
+    Planar(PlanarMapping),
+}
+
+pub struct UVMapping {
+    su: Float,
+    sv: Float,
+    du: Float,
+    dv: Float,
+}
+
+impl Default for UVMapping {
+    fn default() -> Self {
+        Self {
+            su: 1.0,
+            sv: 1.0,
+            du: 0.0,
+            dv: 0.0,
+        }
+    }
+}
+
+impl TextureMapping2DTrait for UVMapping {
+    fn map(&self, ctx: TextureEvalContext) -> TexCoord2D {
+        let dsdx = self.su * ctx.dudx;
+        let dsdy = self.su * ctx.dudy;
+        let dtdx = self.sv * ctx.dvdx;
+        let dtdy = self.sv * ctx.dvdy;
+        let st = Point2f::new(self.su * ctx.uv[0] + self.du, self.sv * ctx.uv[1] * self.dv);
+        TexCoord2D {
+            st,
+            dsdx,
+            dsdy,
+            dtdx,
+            dtdy,
+        }
+    }
+}
+pub struct SphericalMapping {
+    texture_from_render: Transform<Float>,
+}
+
+impl SphericalMapping {
+    pub fn new(texture_from_render: &Transform<Float>) -> Self {
+        Self {
+            texture_from_render: *texture_from_render,
+        }
+    }
+}
+
+impl TextureMapping2DTrait for SphericalMapping {
+    fn map(&self, ctx: TextureEvalContext) -> TexCoord2D {
+        let pt = self.texture_from_render.apply_to_point(ctx.p);
+        let x2y2 = square(pt.x()) + square(pt.y());
+        let sqrtx2y2 = x2y2.sqrt();
+        let dsdp = Vector3f::new(-pt.y(), pt.x(), 0.) / (2. * PI * x2y2);
+        let dtdp = 1. / (PI * (x2y2 * square(pt.z())))
+            * Vector3f::new(
+                pt.x() * pt.z() / sqrtx2y2,
+                pt.y() * pt.z() / sqrtx2y2,
+                -sqrtx2y2,
+            );
+        let dpdx = self.texture_from_render.apply_to_vector(ctx.dpdx);
+        let dpdy = self.texture_from_render.apply_to_vector(ctx.dpdy);
+        let dsdx = dsdp.dot(dpdx);
+        let dsdy = dsdp.dot(dpdy);
+        let dtdx = dtdp.dot(dpdx);
+        let dtdy = dtdp.dot(dpdy);
+        let vec = (pt - Point3f::default()).normalize();
+        let st = Point2f::new(spherical_theta(vec) * INV_PI, spherical_phi(vec) * INV_2_PI);
+        TexCoord2D {
+            st,
+            dsdx,
+            dsdy,
+            dtdx,
+            dtdy,
+        }
+    }
+}
+
+pub struct CylindricalMapping {
+    texture_from_render: Transform<Float>,
+}
+
+impl CylindricalMapping {
+    pub fn new(texture_from_render: &Transform<Float>) -> Self {
+        Self {
+            texture_from_render: *texture_from_render,
+        }
+    }
+}
+
+impl TextureMapping2DTrait for CylindricalMapping {
+    fn map(&self, ctx: TextureEvalContext) -> TexCoord2D {
+        let pt = self.texture_from_render.apply_to_point(ctx.p);
+        let x2y2 = square(pt.x()) + square(pt.y());
+        let dsdp = Vector3f::new(-pt.y(), pt.x(), 0.) / (2. * PI * x2y2);
+        let dtdp = Vector3f::new(1., 0., 0.);
+        let dpdx = self.texture_from_render.apply_to_vector(ctx.dpdx);
+        let dpdy = self.texture_from_render.apply_to_vector(ctx.dpdy);
+        let dsdx = dsdp.dot(dpdx);
+        let dsdy = dsdp.dot(dpdy);
+        let dtdx = dtdp.dot(dpdx);
+        let dtdy = dtdp.dot(dpdy);
+        let st = Point2f::new((PI * pt.y().atan2(pt.x())) * INV_2_PI, pt.z());
+        TexCoord2D {
+            st,
+            dsdx,
+            dsdy,
+            dtdx,
+            dtdy,
+        }
+    }
+}
+
+pub struct PlanarMapping {
+    texture_from_render: Transform<Float>,
+    vs: Vector3f,
+    vt: Vector3f,
+    ds: Float,
+    dt: Float,
+}
+
+impl PlanarMapping {
+    pub fn new(
+        texture_from_render: &Transform<Float>,
+        vs: Vector3f,
+        vt: Vector3f,
+        ds: Float,
+        dt: Float,
+    ) -> Self {
+        Self {
+            texture_from_render: *texture_from_render,
+            vs,
+            vt,
+            ds,
+            dt,
+        }
+    }
+}
+
+impl TextureMapping2DTrait for PlanarMapping {
+    fn map(&self, ctx: TextureEvalContext) -> TexCoord2D {
+        let vec: Vector3f = self.texture_from_render.apply_to_point(ctx.p).into();
+        let dpdx = self.texture_from_render.apply_to_vector(ctx.dpdx);
+        let dpdy = self.texture_from_render.apply_to_vector(ctx.dpdy);
+        let dsdx = self.vs.dot(dpdx);
+        let dsdy = self.vs.dot(dpdy);
+        let dtdx = self.vt.dot(dpdx);
+        let dtdy = self.vt.dot(dpdy);
+        let st = Point2f::new(self.ds + vec.dot(self.vs), self.dt + vec.dot(self.vt));
+        TexCoord2D {
+            st,
+            dsdx,
+            dsdy,
+            dtdx,
+            dtdy,
+        }
+    }
+}
+
+pub struct TexCoord3D {
+    p: Point3f,
+    dpdx: Vector3f,
+    dpdy: Vector3f,
+}
+
+pub trait TextureMapping3DTrait {
+    fn map(&self, ctx: &TextureEvalContext) -> TexCoord3D;
+}
+
+#[derive(Clone, Debug)]
+#[enum_dispatch(TextureMapping3DTrait)]
+pub enum TextureMapping3D {
+    PointTransform(PointTransformMapping),
+}
+
+#[derive(Clone, Debug)]
+pub struct PointTransformMapping {
+    texture_from_render: Transform<Float>,
+}
+
+impl PointTransformMapping {
+    pub fn new(texture_from_render: &Transform<Float>) -> Self {
+        Self {
+            texture_from_render: *texture_from_render,
+        }
+    }
+}
+
+impl TextureMapping3DTrait for PointTransformMapping {
+    fn map(&self, ctx: &TextureEvalContext) -> TexCoord3D {
+        TexCoord3D {
+            p: self.texture_from_render.apply_to_point(ctx.p),
+            dpdx: self.texture_from_render.apply_to_vector(ctx.dpdx),
+            dpdy: self.texture_from_render.apply_to_vector(ctx.dpdy),
+        }
+    }
+}
 
 pub struct TextureEvalContext {
     p: Point3f,
-    dpdx: Vector3f, 
+    dpdx: Vector3f,
     dpdy: Vector3f,
     n: Normal3f,
     uv: Point2f,
     // All 0
     dudx: Float,
-    dudy: Float, 
+    dudy: Float,
-    dvdx: Float, 
+    dvdx: Float,
     dvdy: Float,
     face_index: usize,
 }
 
 impl TextureEvalContext {
-    pub fn new(p: Point3f,
+    #[allow(clippy::too_many_arguments)]
-    dpdx: Vector3f, 
+    pub fn new(
-    dpdy: Vector3f,
+        p: Point3f,
-    n: Normal3f,
+        dpdx: Vector3f,
-    uv: Point2f,
+        dpdy: Vector3f,
-    dudx: Float,
+        n: Normal3f,
-    dudy: Float, 
+        uv: Point2f,
-    dvdx: Float, 
+        dudx: Float,
-    dvdy: Float,
+        dudy: Float,
-    face_index: usize,
+        dvdx: Float,
+        dvdy: Float,
+        face_index: usize,
     ) -> Self {
-        Self {p, dpdx, dpdy, n, uv, dudx, dudy, dvdx, dvdy , face_index }
+        Self {
+            p,
+            dpdx,
+            dpdy,
+            n,
+            uv,
+            dudx,
+            dudy,
+            dvdx,
+            dvdy,
+            face_index,
+        }
     }
 }
 
@@ -49,18 +294,45 @@ impl From<&SurfaceInteraction> for TextureEvalContext {
     }
 }
 
+#[enum_dispatch]
 pub trait FloatTextureTrait: Send + Sync + std::fmt::Debug {
     fn evaluate(&self, _ctx: &TextureEvalContext) -> Float {
         todo!()
     }
 }
 
+#[enum_dispatch(FloatTextureTrait)]
 #[derive(Debug, Clone)]
-pub struct FloatImageTexture;
+pub enum FloatTexture {
+    FloatConstant(FloatConstantTexture),
+    GPUFloatImage(GPUFloatImageTexture),
+    FloatMix(FloatMixTexture),
+    FloatDirectionMix(FloatDirectionMixTexture),
+    FloatScaled(FloatScaledTexture),
+    FloatBilerp(FloatBilerpTexture),
+    FloatCheckerboard(FloatCheckerboardTexture),
+    FloatDots(FloatDotsTexture),
+    FBm(FBmTexture),
+    FloatPtex(FloatPtexTexture),
+    GPUFLoatPtex(GPUFloatPtexTexture),
+    Windy(WindyTexture),
+    Wrinkled(WrinkledTexture),
+}
 
-impl FloatTextureTrait for FloatImageTexture {
+#[derive(Debug, Clone)]
+pub struct FloatConstantTexture {
+    value: Float,
+}
+
+impl FloatConstantTexture {
+    pub fn new(value: Float) -> Self {
+        Self { value }
+    }
+}
+
+impl FloatTextureTrait for FloatConstantTexture {
     fn evaluate(&self, _ctx: &TextureEvalContext) -> Float {
-        todo!();
+        self.value
     }
 }
 
@@ -69,20 +341,56 @@ pub struct GPUFloatImageTexture;
 impl FloatTextureTrait for GPUFloatImageTexture {}
 
 #[derive(Debug, Clone)]
-pub struct FloatMixTexture;
+pub struct FloatMixTexture {
-impl FloatTextureTrait for FloatMixTexture {}
+    tex1: Box<FloatTexture>,
+    tex2: Box<FloatTexture>,
+    amount: Box<FloatTexture>,
+}
+
+impl FloatMixTexture {
+    pub fn new(
+        tex1: Box<FloatTexture>,
+        tex2: Box<FloatTexture>,
+        amount: Box<FloatTexture>,
+    ) -> Self {
+        Self { tex1, tex2, amount }
+    }
+}
+
+impl FloatTextureTrait for FloatMixTexture {
+    fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
+        let amt = self.amount.evaluate(ctx);
+        let mut t1 = 0.;
+        let mut t2 = 0.;
+        if amt != 1. {
+            t1 = self.tex1.evaluate(ctx);
+        }
+        if amt != 0. {
+            t2 = self.tex2.evaluate(ctx);
+        }
+        (1. - amt) * t1 + amt * t2
+    }
+}
 
 #[derive(Debug, Clone)]
 pub struct FloatDirectionMixTexture;
 impl FloatTextureTrait for FloatDirectionMixTexture {}
 
 #[derive(Debug, Clone)]
-pub struct FloatScaledTexture;
+pub struct FloatScaledTexture {
-impl FloatTextureTrait for FloatScaledTexture {}
+    tex: Box<FloatTexture>,
+    scale: Box<FloatTexture>,
+}
 
-#[derive(Debug, Clone)]
+impl FloatTextureTrait for FloatScaledTexture {
-pub struct FloatConstantTexture;
+    fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
-impl FloatTextureTrait for FloatConstantTexture {}
+        let sc = self.scale.evaluate(ctx);
+        if sc == 0. {
+            return 0.;
+        }
+        self.tex.evaluate(ctx)
+    }
+}
 
 #[derive(Debug, Clone)]
 pub struct FloatBilerpTexture;
@@ -105,8 +413,8 @@ pub struct FloatPtexTexture;
 impl FloatTextureTrait for FloatPtexTexture {}
 
 #[derive(Debug, Clone)]
-pub struct GPUFloatPtex;
+pub struct GPUFloatPtexTexture;
-impl FloatTextureTrait for GPUFloatPtex {}
+impl FloatTextureTrait for GPUFloatPtexTexture {}
 
 #[derive(Debug, Clone)]
 pub struct WindyTexture;
@@ -116,60 +424,15 @@ impl FloatTextureTrait for WindyTexture {}
 pub struct WrinkledTexture;
 impl FloatTextureTrait for WrinkledTexture {}
 
-#[derive(Debug, Clone)]
+#[enum_dispatch]
-pub enum FloatTexture {
+pub trait SpectrumTextureTrait: Send + Sync + std::fmt::Debug {
-    FloatImage(FloatImageTexture),
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
-    GPUFloatImage(GPUFloatImageTexture),
+        todo!()
-    FloatMix(FloatMixTexture),
-    FloatDirectionMix(FloatDirectionMixTexture),
-    FloatScaled(FloatScaledTexture),
-    FloatConstant(FloatConstantTexture),
-    FloatBilerp(FloatBilerpTexture),
-    FloatCheckerboard(FloatCheckerboardTexture),
-    FloatDots(FloatDotsTexture),
-    FBm(FBmTexture),
-    FloatPtex(FloatPtexTexture),
-    GPUFloatPtex(GPUFloatPtex),
-    Windy(WindyTexture),
-    Wrinkled(WrinkledTexture),
-}
-
-impl FloatTextureTrait for FloatTexture {
-    fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
-        match self {
-            FloatTexture::FloatImage(texture) => texture.evaluate(ctx),
-            FloatTexture::GPUFloatImage(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatMix(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatDirectionMix(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatScaled(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatConstant(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatBilerp(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatCheckerboard(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatDots(texture) => texture.evaluate(ctx),
-            FloatTexture::FBm(texture) => texture.evaluate(ctx),
-            FloatTexture::FloatPtex(texture) => texture.evaluate(ctx),
-            FloatTexture::GPUFloatPtex(texture) => texture.evaluate(ctx),
-            FloatTexture::Windy(texture) => texture.evaluate(ctx),
-            FloatTexture::Wrinkled(texture) => texture.evaluate(ctx),
-        }
     }
 }
 
-pub struct RGBConstantTexture;
+#[derive(Clone, Debug)]
-pub struct RGBReflectanceConstantTexture;
+#[enum_dispatch(SpectrumTextureTrait)]
-pub struct SpectrumConstantTexture;
-pub struct SpectrumBilerpTexture;
-pub struct SpectrumCheckerboardTexture;
-pub struct SpectrumImageTexture;
-pub struct GPUSpectrumImageTexture;
-pub struct MarbleTexture;
-pub struct SpectrumMixTexture;
-pub struct SpectrumDirectionMixTexture;
-pub struct SpectrumDotsTexture;
-pub struct SpectrumPtexTexture;
-pub struct GPUSpectrumPtexTexture;
-pub struct SpectrumScaledTexture;
-
 pub enum SpectrumTexture {
     RGBConstant(RGBConstantTexture),
     RGBReflectanceConstant(RGBReflectanceConstantTexture),
@@ -187,23 +450,149 @@ pub enum SpectrumTexture {
     SpectrumScaled(SpectrumScaledTexture),
 }
 
-impl SpectrumTexture {
+#[derive(Clone, Debug)]
-    //     pub fn evaluate(&self, ctx: TextureEvalContext) -> f32 {
+pub struct RGBConstantTexture;
-    //         match self {
+impl SpectrumTextureTrait for RGBConstantTexture {
-    //             SpectrumTexture::FloatImage(texture) => texture.evaluate(ctx),
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
-    //             SpectrumTexture::GPUFloatImage(texture) => texture.evaluate(ctx),
+        todo!()
-    //             SpectrumTexture::FloatMix(texture) => texture.evaluate(ctx),
+    }
-    //             SpectrumTexture::FloatDirectionMix(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FloatScaled(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FloatConstant(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FloatBilerp(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FloatCheckerboard(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FloatDots(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FBm(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::FloatPtex(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::GPUFloatPtex(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::Windy(texture) => texture.evaluate(ctx),
-    //             SpectrumTexture::Wrinkled(texture) => texture.evaluate(ctx),
-    //         }
-    //     }
 }
+
+#[derive(Clone, Debug)]
+pub struct RGBReflectanceConstantTexture;
+impl SpectrumTextureTrait for RGBReflectanceConstantTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumConstantTexture {
+    value: Spectrum,
+}
+impl SpectrumTextureTrait for SpectrumConstantTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, lambda: &SampledWavelengths) -> SampledSpectrum {
+        self.value.sample(lambda)
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumBilerpTexture;
+impl SpectrumTextureTrait for SpectrumBilerpTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumCheckerboardTexture;
+impl SpectrumTextureTrait for SpectrumCheckerboardTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumImageTexture;
+impl SpectrumTextureTrait for SpectrumImageTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct GPUSpectrumImageTexture;
+impl SpectrumTextureTrait for GPUSpectrumImageTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct MarbleTexture;
+impl SpectrumTextureTrait for MarbleTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumMixTexture;
+impl SpectrumTextureTrait for SpectrumMixTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumDirectionMixTexture {
+    tex1: Box<SpectrumTexture>,
+    tex2: Box<SpectrumTexture>,
+    dir: Vector3f,
+}
+impl SpectrumTextureTrait for SpectrumDirectionMixTexture {
+    fn evaluate(&self, ctx: &TextureEvalContext, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let amt = ctx.n.abs_dot(self.dir.into());
+        let mut t1 = SampledSpectrum::default();
+        let mut t2 = SampledSpectrum::default();
+        if amt != 0. {
+            t1 = self.tex1.evaluate(ctx, lambda);
+        }
+        if amt != 1. {
+            t2 = self.tex2.evaluate(ctx, lambda);
+        }
+        amt * t1 + (1. - amt) * t2
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumDotsTexture;
+impl SpectrumTextureTrait for SpectrumDotsTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumPtexTexture;
+impl SpectrumTextureTrait for SpectrumPtexTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct GPUSpectrumPtexTexture;
+impl SpectrumTextureTrait for GPUSpectrumPtexTexture {
+    fn evaluate(&self, _ctx: &TextureEvalContext, _lambda: &SampledWavelengths) -> SampledSpectrum {
+        todo!()
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct SpectrumScaledTexture {
+    tex: Box<SpectrumTexture>,
+    scale: Box<FloatTexture>,
+}
+
+impl SpectrumTextureTrait for SpectrumScaledTexture {
+    fn evaluate(&self, ctx: &TextureEvalContext, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let sc = self.scale.evaluate(ctx);
+        if sc == 0. {
+            return SampledSpectrum::new(0.);
+        }
+        self.tex.evaluate(ctx, lambda) * sc
+    }
+}
+
+#[derive(Debug, Hash, PartialEq, Eq, Clone)]
+struct TexInfo {
+    filename: String,
+    filter_options: MIPMapFilterOptions,
+    wrap_mode: WrapMode,
+    encoding: ColorEncoding,
+}
+
+static TEXTURE_CACHE: OnceLock<Mutex<HashMap<TexInfo, Arc<MIPMap>>>> = OnceLock::new();
+
+pub struct ImageTextureBase {}

src/geometry/bounds.rs

@@ -137,15 +137,13 @@ where
     }
 
     pub fn corner(&self, corner_index: usize) -> Point<T, N> {
-        let mut p_arr = [self.p_min[0]; N];
+        Point(std::array::from_fn(|i| {
-        for i in 0..N {
+            if (corner_index >> i) & 1 == 1 {
-            p_arr[i] = if ((corner_index >> i) & 1) == 1 {
                 self.p_max[i]
             } else {
                 self.p_min[i]
             }
-        }
+        }))
-        Point(p_arr)
     }
 
     pub fn overlaps(&self, rhs: &Self) -> bool {

src/geometry/mod.rs

@@ -58,23 +58,21 @@ pub fn tan2_theta(w: Vector3f) -> Float {
 #[inline]
 pub fn cos_phi(w: Vector3f) -> Float {
     let sin_theta = sin_theta(w);
-    let result = if sin_theta == 0. {
+    if sin_theta == 0. {
         1.
     } else {
         clamp_t(w.x() / sin_theta, -1., 1.)
-    };
+    }
-    result
 }
 
 #[inline]
 pub fn sin_phi(w: Vector3f) -> Float {
     let sin_theta = sin_theta(w);
-    let result = if sin_theta == 0. {
+    if sin_theta == 0. {
         0.
     } else {
         clamp_t(w.y() / sin_theta, -1., 1.)
-    };
+    }
-    result
 }
 
 pub fn same_hemisphere(w: Vector3f, wp: Vector3f) -> bool {
@@ -85,7 +83,7 @@ pub fn spherical_direction(sin_theta: Float, cos_theta: Float, phi: Float) -> Ve
     Vector3f::new(sin_theta * phi.cos(), sin_theta * phi.sin(), cos_theta)
 }
 
-pub fn spherical_triangle_area<T: NumFloat>(a: Vector3f, b: Vector3f, c: Vector3f) -> Float {
+pub fn spherical_triangle_area(a: Vector3f, b: Vector3f, c: Vector3f) -> Float {
     (2.0 * (a.dot(b.cross(c))).atan2(1.0 + a.dot(b) + a.dot(c) + b.dot(c))).abs()
 }
 
@@ -114,11 +112,11 @@ pub fn spherical_quad_area(a: Vector3f, b: Vector3f, c: Vector3f, d: Vector3f) -
     (alpha + beta + gamma + delta - 2. * PI).abs()
 }
 
-pub fn spherical_theta<T: NumFloat>(v: Vector3f) -> Float {
+pub fn spherical_theta(v: Vector3f) -> Float {
     clamp_t(v.z(), -1.0, 1.0).acos()
 }
 
-pub fn spherical_phi<T: NumFloat>(v: Vector3f) -> Float {
+pub fn spherical_phi(v: Vector3f) -> Float {
     let p = v.y().atan2(v.x());
     if p < 0.0 { p + 2.0 * PI } else { p }
 }

src/geometry/primitives.rs

@@ -1,4 +1,4 @@
-use super::traits::{Lerp, Sqrt, Tuple, VectorLike};
+use super::traits::{Sqrt, Tuple, VectorLike};
 use super::{Float, NumFloat, PI};
 use crate::utils::interval::Interval;
 use crate::utils::math::{difference_of_products, quadratic, safe_asin};
@@ -39,18 +39,6 @@ macro_rules! impl_tuple_core {
             }
         }
 
-        impl<T, Factor, const N: usize> Lerp<Factor> for $Struct<T, N>
-        where
-            Factor: Copy + Num,
-            T: Copy + Mul<Factor, Output = T> + Add<Output = T>,
-        {
-            #[inline]
-            fn lerp(t: Factor, a: Self, b: Self) -> Self {
-                let result = std::array::from_fn(|i| a[i] * (Factor::one() - t) + b[i] * t);
-                Self::from_array(result)
-            }
-        }
-
         impl<T: Default + Copy, const N: usize> Default for $Struct<T, N> {
             fn default() -> Self {
                 Self([T::default(); N])
@@ -389,12 +377,12 @@ impl Point2f {
 
         if let Some((v0, v1)) = quadratic(k2, k1, k0) {
             let u = (h.x() - f.x() * v0) / (e.x() + g.x() * v0);
-            if u < 0. || u > 1. || v0 < 0. || v0 > 1. {
+            if !(0.0..=1.).contains(&u) || !(0.0..=1.0).contains(&v0) {
                 return Point2f::new((h.x() - f.x() * v1) / (e.x() + g.x() * v1), v1);
             }
-            return Point2f::new(u, v0);
+            Point2f::new(u, v0)
         } else {
-            return Point2f::zero();
+            Point2f::zero()
         }
     }
 }
@@ -699,7 +687,7 @@ where
     T: Num + PartialOrd + Copy + Neg<Output = T> + Sqrt,
 {
     pub fn face_forward(self, v: Vector3<T>) -> Self {
-        if Vector3::<T>::from(self).dot(v.into()) < T::zero() {
+        if Vector3::<T>::from(self).dot(v) < T::zero() {
             -self
         } else {
             self

src/geometry/ray.rs

@@ -42,7 +42,7 @@ impl Ray {
     }
 
     pub fn offset_origin(p: &Point3fi, n: &Normal3f, w: &Vector3f) -> Point3f {
-        let d: Float = Vector3f::from(n.abs()).dot(p.error().into());
+        let d: Float = Vector3f::from(n.abs()).dot(p.error());
         let normal: Vector3f = Vector3f::from(*n);
 
         let mut offset = p.midpoint();

src/geometry/traits.rs

@@ -166,16 +166,18 @@ impl Sqrt for Interval {
     }
 }
 
-pub trait Lerp<Factor: Copy + Num>: Sized + Copy {
+pub trait Lerp<Factor = Float>: Sized + Copy {
     fn lerp(t: Factor, a: Self, b: Self) -> Self;
 }
 
-impl<T> Lerp<T> for T
+impl<T, F, Diff> Lerp<F> for T
 where
-    T: Num + Copy + Mul<T, Output = T> + Add<T, Output = T>,
+    T: Copy + Sub<Output = Diff> + Add<Diff, Output = T>,
+    Diff: Mul<F, Output = Diff>,
+    F: Copy,
 {
-    #[inline]
+    #[inline(always)]
-    fn lerp(t: T, a: Self, b: Self) -> Self {
+    fn lerp(t: F, a: Self, b: Self) -> Self {
-        a * (T::one() - t) + b * t
+        a + (b - a) * t
     }
 }

src/image/io.rs

@@ -0,0 +1,345 @@
+use crate::core::pbrt::Float;
+use crate::image::{Image, ImageMetadata, PixelData, PixelFormat, Point2i, WrapMode};
+use crate::utils::color::{ColorEncoding, LINEAR, SRGB};
+use crate::utils::error::ImageError;
+use anyhow::{Context, Result, bail};
+use exr::prelude::{read_first_rgba_layer_from_file, write_rgba_file};
+use image_rs::ImageReader;
+use image_rs::{DynamicImage, ImageBuffer, Rgb, Rgba};
+use std::fs::File;
+use std::io::{BufRead, BufReader, BufWriter, Read, Write};
+use std::path::Path;
+
+impl Image {
+    pub fn read(path: &Path, encoding: Option<ColorEncoding>) -> Result<Self> {
+        let ext = path
+            .extension()
+            .and_then(|s| s.to_str())
+            .unwrap_or("")
+            .to_lowercase();
+
+        match ext.as_str() {
+            "exr" => read_exr(path),
+            "pfm" => read_pfm(path),
+            _ => read_generic(path, encoding),
+        }
+    }
+
+    pub fn write(&self, filename: &str, metadata: &ImageMetadata) -> Result<(), ImageError> {
+        let path = Path::new(filename);
+        let ext = path.extension().and_then(|s| s.to_str()).unwrap_or("");
+        let res = match ext.to_lowercase().as_str() {
+            "exr" => self.write_exr(path, metadata),
+            "png" => self.write_png(path),
+            "pfm" => self.write_pfm(path),
+            "qoi" => self.write_qoi(path),
+            _ => Err(anyhow::anyhow!("Unsupported write format: {}", ext)),
+        };
+        res.map_err(|e| ImageError::Io(std::io::Error::other(e)))
+    }
+
+    fn write_png(&self, path: &Path) -> Result<()> {
+        let w = self.resolution.x() as u32;
+        let h = self.resolution.y() as u32;
+
+        // Convert whatever we have to u8 [0..255]
+        let data = self.to_u8_buffer();
+        let channels = self.n_channels();
+
+        match channels {
+            1 => {
+                // Luma
+                image_rs::save_buffer_with_format(
+                    path,
+                    &data,
+                    w,
+                    h,
+                    image_rs::ColorType::L8,
+                    image_rs::ImageFormat::Png,
+                )?;
+            }
+            3 => {
+                // RGB
+                image_rs::save_buffer_with_format(
+                    path,
+                    &data,
+                    w,
+                    h,
+                    image_rs::ColorType::Rgb8,
+                    image_rs::ImageFormat::Png,
+                )?;
+            }
+            4 => {
+                // RGBA
+                image_rs::save_buffer_with_format(
+                    path,
+                    &data,
+                    w,
+                    h,
+                    image_rs::ColorType::Rgba8,
+                    image_rs::ImageFormat::Png,
+                )?;
+            }
+            _ => bail!("PNG writer only supports 1, 3, or 4 channels"),
+        }
+        Ok(())
+    }
+
+    fn write_qoi(&self, path: &Path) -> Result<()> {
+        let w = self.resolution.x() as u32;
+        let h = self.resolution.y() as u32;
+        let data = self.to_u8_buffer();
+
+        let color_type = match self.n_channels() {
+            3 => image_rs::ColorType::Rgb8,
+            4 => image_rs::ColorType::Rgba8,
+            _ => bail!("QOI only supports 3 or 4 channels"),
+        };
+
+        image_rs::save_buffer_with_format(
+            path,
+            &data,
+            w,
+            h,
+            color_type,
+            image_rs::ImageFormat::Qoi,
+        )?;
+        Ok(())
+    }
+
+    fn write_exr(&self, path: &Path, _metadata: &ImageMetadata) -> Result<()> {
+        // EXR requires F32. If we have others, we rely on the float accessors or convert.
+        let w = self.resolution.x() as usize;
+        let h = self.resolution.y() as usize;
+        let c = self.n_channels();
+
+        write_rgba_file(path, w, h, |x, y| {
+            // Helper to get float value regardless of internal storage
+            let get =
+                |ch| self.get_channel(Point2i::new(x as i32, y as i32), ch, WrapMode::Clamp.into());
+
+            if c == 1 {
+                let v = get(0);
+                (v, v, v, 1.0)
+            } else if c == 3 {
+                (get(0), get(1), get(2), 1.0)
+            } else {
+                (get(0), get(1), get(2), get(3))
+            }
+        })
+        .context("Failed to write EXR")?;
+
+        Ok(())
+    }
+
+    fn write_pfm(&self, path: &Path) -> Result<()> {
+        let file = File::create(path)?;
+        let mut writer = BufWriter::new(file);
+
+        if self.n_channels() != 3 {
+            bail!("PFM writing currently only supports 3 channels (RGB)");
+        }
+
+        // Header
+        writeln!(writer, "PF")?;
+        writeln!(writer, "{} {}", self.resolution.x(), self.resolution.y())?;
+        // Scale: Negative means Little Endian
+        let scale = if cfg!(target_endian = "little") {
+            -1.0
+        } else {
+            1.0
+        };
+        writeln!(writer, "{}", scale)?;
+
+        // PFM stores bottom-to-top. PBRT stores top-to-bottom.
+        for y in (0..self.resolution.y()).rev() {
+            for x in 0..self.resolution.x() {
+                for c in 0..3 {
+                    let val = self.get_channel(Point2i::new(x, y), c, WrapMode::Clamp.into());
+                    writer.write_all(&val.to_le_bytes())?;
+                }
+            }
+        }
+
+        Ok(())
+    }
+
+    fn to_u8_buffer(&self) -> Vec<u8> {
+        match &self.pixels {
+            PixelData::U8(data) => data.clone(),
+            PixelData::F16(data) => data
+                .iter()
+                .map(|v| (v.to_f32().clamp(0.0, 1.0) * 255.0 + 0.5) as u8)
+                .collect(),
+            PixelData::F32(data) => data
+                .iter()
+                .map(|v| (v.clamp(0.0, 1.0) * 255.0 + 0.5) as u8)
+                .collect(),
+        }
+    }
+}
+
+fn read_generic(path: &Path, encoding: Option<ColorEncoding>) -> Result<Image> {
+    // 1. Load via 'image' crate
+    let dyn_img = ImageReader::open(path)
+        .with_context(|| format!("Failed to open image: {:?}", path))?
+        .decode()?;
+
+    let w = dyn_img.width() as i32;
+    let h = dyn_img.height() as i32;
+    let res = Point2i::new(w, h);
+
+    // 2. Convert to PBRT Image
+    // Check if it was loaded as high precision (HDR/EXR via image crate) or standard
+    match dyn_img {
+        DynamicImage::ImageRgb32F(buf) => Ok(Image {
+            format: PixelFormat::F32,
+            resolution: res,
+            channel_names: vec!["R".into(), "G".into(), "B".into()],
+            encoding: LINEAR,
+            pixels: PixelData::F32(buf.into_raw()),
+        }),
+        DynamicImage::ImageRgba32F(buf) => Ok(Image {
+            format: PixelFormat::F32,
+            resolution: res,
+            channel_names: vec!["R".into(), "G".into(), "B".into(), "A".into()],
+            encoding: LINEAR,
+            pixels: PixelData::F32(buf.into_raw()),
+        }),
+        _ => {
+            // Default to RGB8 for everything else (PNG, JPG)
+            // Note: PBRT handles alpha, so we check if the source had alpha
+            if dyn_img.color().has_alpha() {
+                let buf = dyn_img.to_rgba8();
+                Ok(Image {
+                    format: PixelFormat::U8,
+                    resolution: res,
+                    channel_names: vec!["R".into(), "G".into(), "B".into(), "A".into()],
+                    encoding: encoding.unwrap_or(SRGB),
+                    pixels: PixelData::U8(buf.into_raw()),
+                })
+            } else {
+                let buf = dyn_img.to_rgb8();
+                Ok(Image {
+                    format: PixelFormat::U8,
+                    resolution: res,
+                    channel_names: vec!["R".into(), "G".into(), "B".into()],
+                    encoding: encoding.unwrap_or(SRGB),
+                    pixels: PixelData::U8(buf.into_raw()),
+                })
+            }
+        }
+    }
+}
+
+fn read_exr(path: &Path) -> Result<Image> {
+    let image = read_first_rgba_layer_from_file(
+        path,
+        |resolution, _| {
+            let size = resolution.width() * resolution.height() * 4;
+            vec![0.0 as Float; size]
+        },
+        |buffer, position, pixel| {
+            let width = position.width();
+            let idx = (position.y() * width + position.x()) * 4;
+            // Map exr pixel struct to our buffer
+            buffer[idx] = pixel.0;
+            buffer[idx + 1] = pixel.1;
+            buffer[idx + 2] = pixel.2;
+            buffer[idx + 3] = pixel.3;
+        },
+    )
+    .with_context(|| format!("Failed to read EXR: {:?}", path))?;
+
+    let w = image.layer_data.size.width() as i32;
+    let h = image.layer_data.size.height() as i32;
+
+    Ok(Image {
+        format: PixelFormat::F32,
+        resolution: Point2i::new(w, h),
+        channel_names: vec!["R".into(), "G".into(), "B".into(), "A".into()],
+        encoding: LINEAR,
+        pixels: PixelData::F32(image.layer_data.channel_data.pixels),
+    })
+}
+
+fn read_pfm(path: &Path) -> Result<Image> {
+    let file = File::open(path)?;
+    let mut reader = BufReader::new(file);
+
+    // PFM Headers are: "PF\nwidth height\nscale\n" (or Pf for grayscale)
+    let mut header_word = String::new();
+    reader.read_line(&mut header_word)?;
+    let header_word = header_word.trim();
+
+    let channels = match header_word {
+        "PF" => 3,
+        "Pf" => 1,
+        _ => bail!("Invalid PFM header: {}", header_word),
+    };
+
+    let mut dims_line = String::new();
+    reader.read_line(&mut dims_line)?;
+    let dims: Vec<i32> = dims_line
+        .split_whitespace()
+        .map(|s| s.parse().unwrap_or(0))
+        .collect();
+
+    if dims.len() < 2 {
+        bail!("Invalid PFM dimensions");
+    }
+    let w = dims[0];
+    let h = dims[1];
+
+    let mut scale_line = String::new();
+    reader.read_line(&mut scale_line)?;
+    let scale: f32 = scale_line.trim().parse().context("Invalid PFM scale")?;
+
+    let file_is_little_endian = scale < 0.0;
+    let abs_scale = scale.abs();
+
+    let mut buffer = Vec::new();
+    reader.read_to_end(&mut buffer)?;
+
+    let expected_bytes = (w * h * channels) as usize * 4;
+    if buffer.len() < expected_bytes {
+        bail!("PFM file too short");
+    }
+
+    let mut pixels = vec![0.0 as Float; (w * h * channels) as usize];
+
+    // PFM is Bottom-to-Top. We must flip Y while reading.
+    for y in 0..h {
+        let src_y = h - 1 - y; // Flip logic
+        for x in 0..w {
+            for c in 0..channels {
+                let src_idx = ((src_y * w + x) * channels + c) as usize * 4;
+                let dst_idx = ((y * w + x) * channels + c) as usize;
+
+                let bytes: [u8; 4] = buffer[src_idx..src_idx + 4].try_into()?;
+
+                let val = if file_is_little_endian {
+                    f32::from_le_bytes(bytes)
+                } else {
+                    f32::from_be_bytes(bytes)
+                };
+
+                pixels[dst_idx] = val * abs_scale;
+            }
+        }
+    }
+
+    let names = if channels == 1 {
+        vec!["Y".into()]
+    } else {
+        vec!["R".into(), "G".into(), "B".into()]
+    };
+
+    Ok(Image {
+        format: PixelFormat::F32,
+        resolution: Point2i::new(w, h),
+        channel_names: names,
+        encoding: LINEAR,
+        pixels: PixelData::F32(pixels),
+    })
+}

src/image/metadata.rs

@@ -0,0 +1,60 @@
+use crate::core::pbrt::Float;
+use crate::geometry::{Bounds2i, Point2i};
+use crate::utils::colorspace::RGBColorSpace;
+use crate::utils::math::SquareMatrix;
+use smallvec::SmallVec;
+use std::collections::HashMap;
+
+pub type ImageChannelValues = SmallVec<[Float; 4]>;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+pub enum WrapMode {
+    Black,
+    Clamp,
+    Repeat,
+    OctahedralSphere,
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+pub struct WrapMode2D {
+    pub uv: [WrapMode; 2],
+}
+
+impl From<WrapMode> for WrapMode2D {
+    fn from(w: WrapMode) -> Self {
+        Self { uv: [w, w] }
+    }
+}
+
+#[derive(Debug, Clone, Default)]
+pub struct ImageChannelDesc {
+    pub offset: Vec<usize>,
+}
+
+impl ImageChannelDesc {
+    pub fn new(offset: &[usize]) -> Self {
+        Self {
+            offset: offset.into(),
+        }
+    }
+    pub fn len(&self) -> usize {
+        self.offset.len()
+    }
+    pub fn is_empty(&self) -> bool {
+        self.offset.is_empty()
+    }
+}
+
+#[derive(Debug, Default)]
+pub struct ImageMetadata {
+    pub render_time_seconds: Option<Float>,
+    pub camera_from_world: Option<SquareMatrix<Float, 4>>,
+    pub ndc_from_world: Option<SquareMatrix<Float, 4>>,
+    pub pixel_bounds: Option<Bounds2i>,
+    pub full_resolution: Option<Point2i>,
+    pub samples_per_pixel: Option<i32>,
+    pub mse: Option<Float>,
+    pub color_space: Option<&'static RGBColorSpace>,
+    pub strings: HashMap<String, String>,
+    pub string_vectors: HashMap<String, Vec<String>>,
+}

src/image/mod.rs

@@ -0,0 +1,197 @@
+pub mod io;
+pub mod metadata;
+pub mod ops;
+pub mod pixel;
+
+use crate::core::pbrt::{Float, lerp};
+use crate::geometry::{Point2f, Point2i};
+use crate::utils::color::{ColorEncoding, ColorEncodingTrait};
+use half::f16;
+use pixel::PixelStorage;
+
+pub use metadata::{ImageChannelDesc, ImageChannelValues, ImageMetadata, WrapMode, WrapMode2D};
+
+#[derive(Debug, Copy, Clone, PartialEq, Eq)]
+pub enum PixelFormat {
+    U8,
+    F16,
+    F32,
+}
+
+impl PixelFormat {
+    pub fn is_8bit(&self) -> bool {
+        matches!(self, PixelFormat::U8)
+    }
+
+    pub fn is_16bit(&self) -> bool {
+        matches!(self, PixelFormat::F16)
+    }
+
+    pub fn is_32bit(&self) -> bool {
+        matches!(self, PixelFormat::F32)
+    }
+
+    pub fn texel_bytes(&self) -> usize {
+        match self {
+            PixelFormat::U8 => 1,
+            PixelFormat::F16 => 2,
+            PixelFormat::F32 => 4,
+        }
+    }
+}
+
+impl std::fmt::Display for PixelFormat {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            PixelFormat::U8 => write!(f, "U256"),
+            PixelFormat::F16 => write!(f, "Half"),
+            PixelFormat::F32 => write!(f, "Float"),
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub enum PixelData {
+    U8(Vec<u8>),
+    F16(Vec<f16>),
+    F32(Vec<f32>),
+}
+
+#[derive(Debug, Clone)]
+pub struct Image {
+    pub format: PixelFormat,
+    pub resolution: Point2i,
+    pub channel_names: Vec<String>,
+    pub encoding: ColorEncoding,
+    pub pixels: PixelData,
+}
+
+#[derive(Debug)]
+pub struct ImageAndMetadata {
+    pub image: Image,
+    pub metadata: ImageMetadata,
+}
+
+impl Image {
+    pub fn new(
+        format: PixelFormat,
+        resolution: Point2i,
+        channel_names: Vec<String>,
+        encoding: ColorEncoding,
+    ) -> Self {
+        let size = (resolution.x() * resolution.y()) as usize * channel_names.len();
+        let pixels = match format {
+            PixelFormat::U8 => PixelData::U8(vec![0; size]),
+            PixelFormat::F16 => PixelData::F16(vec![f16::ZERO; size]),
+            PixelFormat::F32 => PixelData::F32(vec![0.0; size]),
+        };
+        Self {
+            resolution,
+            channel_names,
+            format,
+            pixels,
+            encoding,
+        }
+    }
+
+    pub fn format(&self) -> PixelFormat {
+        self.format
+    }
+    pub fn resolution(&self) -> Point2i {
+        self.resolution
+    }
+    pub fn n_channels(&self) -> usize {
+        self.channel_names.len()
+    }
+    pub fn channel_names(&self) -> &[String] {
+        &self.channel_names
+    }
+    pub fn encoding(&self) -> ColorEncoding {
+        self.encoding
+    }
+
+    pub fn pixel_offset(&self, p: Point2i) -> usize {
+        (p.y() as usize * self.resolution.x() as usize + p.x() as usize) * self.n_channels()
+    }
+
+    pub fn get_channel(&self, p: Point2i, c: usize, wrap: WrapMode2D) -> Float {
+        let mut pp = p;
+        if !self.remap_pixel_coords(&mut pp, wrap) {
+            return 0.0;
+        }
+
+        let idx = self.pixel_offset(pp) + c;
+        match &self.pixels {
+            PixelData::U8(d) => u8::to_linear(d[idx], self.encoding),
+            PixelData::F16(d) => f16::to_linear(d[idx], self.encoding),
+            PixelData::F32(d) => f32::to_linear(d[idx], self.encoding),
+        }
+    }
+
+    pub fn all_channels_desc(&self) -> ImageChannelDesc {
+        ImageChannelDesc {
+            offset: (0..self.n_channels()).collect(),
+        }
+    }
+
+    pub fn set_channel(&mut self, p: Point2i, c: usize, value: Float) {
+        let val_no_nan = if value.is_nan() { 0.0 } else { value };
+        let offset = self.pixel_offset(p) + c;
+        match &mut self.pixels {
+            PixelData::U8(data) => {
+                let linear = [val_no_nan];
+                self.encoding
+                    .from_linear_slice(&linear, &mut data[offset..offset + 1]);
+            }
+            PixelData::F16(data) => data[offset] = f16::from_f32(val_no_nan),
+            PixelData::F32(data) => data[offset] = val_no_nan,
+        }
+    }
+
+    pub fn set_channels(
+        &mut self,
+        p: Point2i,
+        desc: &ImageChannelDesc,
+        values: &ImageChannelValues,
+    ) {
+        assert_eq!(desc.len(), values.len());
+        for i in 0..desc.len() {
+            self.set_channel(p, desc.offset[i], values[i]);
+        }
+    }
+
+    pub fn set_channels_all(&mut self, p: Point2i, values: &ImageChannelValues) {
+        self.set_channels(p, &self.all_channels_desc(), values)
+    }
+
+    fn remap_pixel_coords(&self, p: &mut Point2i, wrap_mode: WrapMode2D) -> bool {
+        for i in 0..2 {
+            if p[i] >= 0 && p[i] < self.resolution[i] {
+                continue;
+            }
+            match wrap_mode.uv[i] {
+                WrapMode::Black => return false,
+                WrapMode::Clamp => p[i] = p[i].clamp(0, self.resolution[i] - 1),
+                WrapMode::Repeat => p[i] = p[i].rem_euclid(self.resolution[i]),
+                WrapMode::OctahedralSphere => {
+                    p[i] = p[i].clamp(0, self.resolution[i] - 1);
+                }
+            }
+        }
+        true
+    }
+
+    pub fn bilerp_channel(&self, p: Point2f, c: usize, wrap_mode: WrapMode2D) -> Float {
+        let x = p.x() * self.resolution.x() as Float - 0.5;
+        let y = p.y() * self.resolution.y() as Float - 0.5;
+        let xi = x.floor() as i32;
+        let yi = y.floor() as i32;
+        let dx = x - xi as Float;
+        let dy = y - yi as Float;
+        let v00 = self.get_channel(Point2i::new(xi, yi), c, wrap_mode);
+        let v10 = self.get_channel(Point2i::new(xi + 1, yi), c, wrap_mode);
+        let v01 = self.get_channel(Point2i::new(xi, yi + 1), c, wrap_mode);
+        let v11 = self.get_channel(Point2i::new(xi + 1, yi + 1), c, wrap_mode);
+        lerp(dy, lerp(dx, v00, v10), lerp(dx, v01, v11))
+    }
+}

src/image/ops.rs

@@ -0,0 +1,418 @@
+// use rayon::prelude::*;
+use crate::core::pbrt::Float;
+use crate::geometry::{Bounds2i, Point2i};
+use crate::image::pixel::PixelStorage;
+use crate::image::{Image, PixelData, PixelFormat, WrapMode, WrapMode2D};
+use crate::utils::math::windowed_sinc;
+use rayon::prelude::*;
+use std::sync::{Arc, Mutex};
+
+#[derive(Debug, Clone, Copy)]
+pub struct ResampleWeight {
+    pub first_pixel: i32,
+    pub weight: [Float; 4],
+}
+
+impl Image {
+    pub fn flip_y(&mut self) {
+        let res = self.resolution;
+        let nc = self.n_channels();
+
+        match &mut self.pixels {
+            PixelData::U8(d) => flip_y_kernel(d, res, nc),
+            PixelData::F16(d) => flip_y_kernel(d, res, nc),
+            PixelData::F32(d) => flip_y_kernel(d, res, nc),
+        }
+    }
+
+    pub fn crop(&self, bounds: Bounds2i) -> Image {
+        let new_res = Point2i::new(
+            bounds.p_max.x() - bounds.p_min.x(),
+            bounds.p_max.y() - bounds.p_min.y(),
+        );
+
+        let mut new_image = Image::new(
+            self.format,
+            new_res,
+            self.channel_names.clone(),
+            self.encoding,
+        );
+
+        match (&self.pixels, &mut new_image.pixels) {
+            (PixelData::U8(src), PixelData::U8(dst)) => {
+                crop_kernel(src, dst, self.resolution, bounds, self.n_channels())
+            }
+            (PixelData::F16(src), PixelData::F16(dst)) => {
+                crop_kernel(src, dst, self.resolution, bounds, self.n_channels())
+            }
+            (PixelData::F32(src), PixelData::F32(dst)) => {
+                crop_kernel(src, dst, self.resolution, bounds, self.n_channels())
+            }
+            _ => panic!("Format mismatch in crop"),
+        }
+
+        new_image
+    }
+
+    pub fn copy_rect_out(&self, extent: Bounds2i, buf: &mut [Float], wrap: WrapMode2D) {
+        match &self.pixels {
+            PixelData::U8(d) => copy_rect_out_kernel(d, self, extent, buf, wrap),
+            PixelData::F16(d) => copy_rect_out_kernel(d, self, extent, buf, wrap),
+            PixelData::F32(d) => copy_rect_out_kernel(d, self, extent, buf, wrap),
+        }
+    }
+
+    pub fn copy_rect_in(&mut self, extent: Bounds2i, buf: &[Float]) {
+        let resolution = self.resolution;
+        let n_channels = self.n_channels();
+        let encoding = self.encoding;
+
+        match &mut self.pixels {
+            PixelData::U8(d) => {
+                copy_rect_in_kernel(d, resolution, n_channels, encoding, extent, buf)
+            }
+            PixelData::F16(d) => {
+                copy_rect_in_kernel(d, resolution, n_channels, encoding, extent, buf)
+            }
+            PixelData::F32(d) => {
+                copy_rect_in_kernel(d, resolution, n_channels, encoding, extent, buf)
+            }
+        }
+    }
+
+    pub fn float_resize_up(&self, new_res: Point2i, wrap_mode: WrapMode2D) -> Image {
+        assert!(new_res.x() >= self.resolution.x() && new_res.y() >= self.resolution.y());
+        assert!(
+            matches!(self.format, PixelFormat::F32),
+            "ResizeUp requires Float format"
+        );
+
+        // Create destination image wrapped in Mutex for tile-based write access
+        let resampled_image = Arc::new(Mutex::new(Image::new(
+            PixelFormat::F32, // Force float output
+            new_res,
+            self.channel_names.clone(),
+            self.encoding,
+        )));
+
+        // Precompute weights
+        let x_weights = resample_weights(self.resolution.x() as usize, new_res.x() as usize);
+        let y_weights = resample_weights(self.resolution.y() as usize, new_res.y() as usize);
+        let n_channels = self.n_channels();
+
+        // Generate Tiles
+        let tile_size = 16;
+        let tiles = generate_tiles(new_res, tile_size);
+
+        // Parallel Execution
+        tiles.par_iter().for_each(|out_extent| {
+            let x_start = x_weights[out_extent.p_min.x() as usize].first_pixel;
+            let x_end = x_weights[(out_extent.p_max.x() - 1) as usize].first_pixel + 4;
+            let y_start = y_weights[out_extent.p_min.y() as usize].first_pixel;
+            let y_end = y_weights[(out_extent.p_max.y() - 1) as usize].first_pixel + 4;
+
+            let in_extent =
+                Bounds2i::from_points(Point2i::new(x_start, y_start), Point2i::new(x_end, y_end));
+
+            let mut in_buf = vec![0.0; in_extent.area() as usize * n_channels];
+            self.copy_rect_out(in_extent, &mut in_buf, wrap_mode);
+
+            let out_buf = compute_resize_tile(
+                &in_buf,
+                in_extent,
+                *out_extent,
+                n_channels,
+                &x_weights,
+                &y_weights,
+            );
+
+            let mut guard = resampled_image.lock().unwrap();
+            guard.copy_rect_in(*out_extent, &out_buf);
+        });
+
+        Arc::try_unwrap(resampled_image)
+            .unwrap()
+            .into_inner()
+            .unwrap()
+    }
+
+    pub fn generate_pyramid(base: Image, _wrap: WrapMode) -> Vec<Image> {
+        let mut levels = vec![base];
+        let internal_wrap = WrapMode2D {
+            uv: [WrapMode::Clamp; 2],
+        };
+
+        loop {
+            let prev = levels.last().unwrap();
+            let old = prev.resolution;
+            if old.x() == 1 && old.y() == 1 {
+                break;
+            }
+
+            let new_res = Point2i::new((old.x() / 2).max(1), (old.y() / 2).max(1));
+            let mut next = Image::new(
+                prev.format,
+                new_res,
+                prev.channel_names.clone(),
+                prev.encoding,
+            );
+
+            match &mut next.pixels {
+                PixelData::U8(d) => downsample_kernel(d, new_res, prev, internal_wrap),
+                PixelData::F16(d) => downsample_kernel(d, new_res, prev, internal_wrap),
+                PixelData::F32(d) => downsample_kernel(d, new_res, prev, internal_wrap),
+            }
+            levels.push(next);
+        }
+        levels
+    }
+}
+
+fn flip_y_kernel<T: PixelStorage>(pixels: &mut [T], res: Point2i, channels: usize) {
+    let w = res.x() as usize;
+    let h = res.y() as usize;
+    let stride = w * channels;
+
+    for y in 0..(h / 2) {
+        let bot = h - 1 - y;
+        for i in 0..stride {
+            pixels.swap(y * stride + i, bot * stride + i);
+        }
+    }
+}
+
+fn crop_kernel<T: PixelStorage>(
+    src: &[T],
+    dst: &mut [T],
+    src_res: Point2i,
+    bounds: Bounds2i,
+    channels: usize,
+) {
+    let dst_w = (bounds.p_max.x() - bounds.p_min.x()) as usize;
+    // let dst_h = (bounds.p_max.y() - bounds.p_min.y()) as usize;
+
+    dst.par_chunks_mut(dst_w * channels)
+        .enumerate()
+        .for_each(|(dy, dst_row)| {
+            let sy = bounds.p_min.y() as usize + dy;
+            let sx_start = bounds.p_min.x() as usize;
+
+            let src_offset = (sy * src_res.x() as usize + sx_start) * channels;
+            let count = dst_w * channels;
+
+            dst_row.copy_from_slice(&src[src_offset..src_offset + count]);
+        });
+}
+
+fn copy_rect_out_kernel<T: PixelStorage>(
+    src: &[T],
+    image: &Image,
+    extent: Bounds2i,
+    buf: &mut [Float],
+    wrap: WrapMode2D,
+) {
+    let w = (extent.p_max.x() - extent.p_min.x()) as usize;
+    let channels = image.n_channels();
+    let enc = image.encoding;
+    let res = image.resolution;
+
+    buf.par_chunks_mut(w * channels)
+        .enumerate()
+        .for_each(|(y_rel, row_buf)| {
+            let y = extent.p_min.y() + y_rel as i32;
+            for x_rel in 0..w {
+                let x = extent.p_min.x() + x_rel as i32;
+
+                // This allows us to use 'src' directly (Fast Path).
+                if x >= 0 && x < res.x() && y >= 0 && y < res.y() {
+                    let offset = (y as usize * res.x() as usize + x as usize) * channels;
+
+                    for c in 0..channels {
+                        row_buf[x_rel * channels + c] = T::to_linear(src[offset + c], enc);
+                    }
+                } else {
+                    // Slow path: Out of bounds, requires Wrap Mode logic.
+                    // We fall back to get_channel which handles the wrapping math.
+                    let p = Point2i::new(x, y);
+                    for c in 0..channels {
+                        row_buf[x_rel * channels + c] = image.get_channel(p, c, wrap);
+                    }
+                }
+            }
+        });
+}
+
+fn copy_rect_in_kernel<T: PixelStorage>(
+    dst: &mut [T],
+    res: Point2i,
+    channels: usize,
+    enc: crate::utils::color::ColorEncoding,
+    extent: Bounds2i,
+    buf: &[Float],
+) {
+    let w = (extent.p_max.x() - extent.p_min.x()) as usize;
+    let res_x = res.x() as usize;
+
+    let rows = buf.chunks(w * channels);
+    for (y_rel, row) in rows.enumerate() {
+        let y = extent.p_min.y() + y_rel as i32;
+        if y < 0 || y >= res.y() {
+            continue;
+        }
+
+        let dst_row_start = (y as usize * res_x) * channels;
+
+        for (x_rel, &val) in row.iter().enumerate() {
+            let c = x_rel % channels;
+            let x_pixel = x_rel / channels;
+            let x = extent.p_min.x() + x_pixel as i32;
+
+            if x >= 0 && x < res.x() {
+                let idx = dst_row_start + (x as usize * channels) + c;
+                dst[idx] = T::from_linear(val, enc);
+            }
+        }
+    }
+}
+
+fn downsample_kernel<T: PixelStorage>(
+    dst: &mut [T],
+    dst_res: Point2i,
+    prev: &Image,
+    wrap: WrapMode2D,
+) {
+    let w = dst_res.x() as usize;
+    let channels = prev.n_channels();
+    let enc = prev.encoding;
+    let old_res = prev.resolution;
+
+    dst.par_chunks_mut(w * channels)
+        .enumerate()
+        .for_each(|(y, row)| {
+            let src_y = y * 2;
+            for x in 0..w {
+                let src_x = x * 2;
+                for c in 0..channels {
+                    let mut sum = 0.0;
+                    let mut count = 0.0;
+
+                    for dy in 0..2 {
+                        for dx in 0..2 {
+                            let sx = src_x as i32 + dx;
+                            let sy = src_y as i32 + dy;
+                            if sx < old_res.x() && sy < old_res.y() {
+                                sum += prev.get_channel(Point2i::new(sx, sy), c, wrap);
+                                count += 1.0;
+                            }
+                        }
+                    }
+
+                    let avg = if count > 0.0 { sum / count } else { 0.0 };
+                    row[x * channels + c] = T::from_linear(avg, enc);
+                }
+            }
+        });
+}
+
+fn resample_weights(old_res: usize, new_res: usize) -> Vec<ResampleWeight> {
+    let filter_radius = 2.0;
+    let tau = 2.0;
+    (0..new_res)
+        .map(|i| {
+            let center = (i as Float + 0.5) * old_res as Float / new_res as Float;
+            let first_pixel = ((center - filter_radius) + 0.5).floor() as i32;
+            let mut weights = [0.0; 4];
+            let mut sum = 0.0;
+            for j in 0..4 {
+                let pos = (first_pixel + j) as Float + 0.5;
+                weights[j as usize] = windowed_sinc(pos - center, filter_radius, tau);
+                sum += weights[j as usize];
+            }
+            let inv_sum = 1.0 / sum;
+            for w in &mut weights {
+                *w *= inv_sum;
+            }
+            ResampleWeight {
+                first_pixel,
+                weight: weights,
+            }
+        })
+        .collect()
+}
+
+fn generate_tiles(res: Point2i, tile_size: i32) -> Vec<Bounds2i> {
+    let nx = (res.x() + tile_size - 1) / tile_size;
+    let ny = (res.y() + tile_size - 1) / tile_size;
+    let mut tiles = Vec::new();
+    for y in 0..ny {
+        for x in 0..nx {
+            let p_min = Point2i::new(x * tile_size, y * tile_size);
+            let p_max = Point2i::new(
+                (x * tile_size + tile_size).min(res.x()),
+                (y * tile_size + tile_size).min(res.y()),
+            );
+            tiles.push(Bounds2i::from_points(p_min, p_max));
+        }
+    }
+    tiles
+}
+
+fn compute_resize_tile(
+    in_buf: &[Float],
+    in_extent: Bounds2i,
+    out_extent: Bounds2i,
+    n_channels: usize,
+    x_w: &[ResampleWeight],
+    y_w: &[ResampleWeight],
+) -> Vec<Float> {
+    let nx_out = (out_extent.p_max.x() - out_extent.p_min.x()) as usize;
+    let ny_out = (out_extent.p_max.y() - out_extent.p_min.y()) as usize;
+    let nx_in = (in_extent.p_max.x() - in_extent.p_min.x()) as usize;
+    let ny_in = (in_extent.p_max.y() - in_extent.p_min.y()) as usize;
+
+    let mut x_buf = vec![0.0; n_channels * ny_in * nx_out];
+
+    // Resize X
+    for y in 0..ny_in {
+        for x in 0..nx_out {
+            let x_global = out_extent.p_min.x() + x as i32;
+            let w = &x_w[x_global as usize];
+            let x_in_start = (w.first_pixel - in_extent.p_min.x()) as usize;
+
+            let in_idx_base = (y * nx_in + x_in_start) * n_channels;
+            let out_idx_base = (y * nx_out + x) * n_channels;
+
+            for c in 0..n_channels {
+                let mut val = 0.0;
+                for k in 0..4 {
+                    val += w.weight[k] * in_buf[in_idx_base + k * n_channels + c];
+                }
+                x_buf[out_idx_base + c] = val;
+            }
+        }
+    }
+
+    let mut out_buf = vec![0.0; n_channels * nx_out * ny_out];
+
+    // Resize Y
+    for x in 0..nx_out {
+        for y in 0..ny_out {
+            let y_global = out_extent.p_min.y() + y as i32;
+            let w = &y_w[y_global as usize];
+            let y_in_start = (w.first_pixel - in_extent.p_min.y()) as usize;
+
+            let in_idx_base = (y_in_start * nx_out + x) * n_channels;
+            let out_idx_base = (y * nx_out + x) * n_channels;
+            let stride = nx_out * n_channels;
+
+            for c in 0..n_channels {
+                let mut val = 0.0;
+                for k in 0..4 {
+                    val += w.weight[k] * x_buf[in_idx_base + k * stride + c];
+                }
+                out_buf[out_idx_base + c] = val.max(0.0);
+            }
+        }
+    }
+    out_buf
+}

src/image/pixel.rs

@@ -0,0 +1,47 @@
+use crate::core::pbrt::Float;
+use crate::utils::color::{ColorEncoding, ColorEncodingTrait};
+use half::f16;
+
+/// Allows writing generic algorithms that work on any image format.
+pub trait PixelStorage: Copy + Send + Sync + 'static + PartialEq {
+    fn from_linear(val: Float, encoding: ColorEncoding) -> Self;
+    fn to_linear(self, encoding: ColorEncoding) -> Float;
+}
+
+impl PixelStorage for f32 {
+    #[inline(always)]
+    fn from_linear(val: Float, _enc: ColorEncoding) -> Self {
+        val
+    }
+    #[inline(always)]
+    fn to_linear(self, _enc: ColorEncoding) -> Float {
+        self
+    }
+}
+
+impl PixelStorage for f16 {
+    #[inline(always)]
+    fn from_linear(val: Float, _enc: ColorEncoding) -> Self {
+        f16::from_f32(val)
+    }
+    #[inline(always)]
+    fn to_linear(self, _enc: ColorEncoding) -> Float {
+        self.to_f32()
+    }
+}
+
+impl PixelStorage for u8 {
+    #[inline(always)]
+    fn from_linear(val: Float, enc: ColorEncoding) -> Self {
+        let mut out = [0u8];
+        enc.from_linear_slice(&[val], &mut out);
+        out[0]
+    }
+
+    #[inline(always)]
+    fn to_linear(self, enc: ColorEncoding) -> Float {
+        let mut out = [0.0];
+        enc.to_linear_slice(&[self], &mut out);
+        out[0]
+    }
+}

src/lib.rs

@@ -5,6 +5,8 @@
 mod camera;
 mod core;
 mod geometry;
+mod image;
 mod lights;
 mod shapes;
+mod spectra;
 mod utils;

src/lights/diffuse.rs

@@ -1,7 +1,7 @@
 use crate::core::medium::MediumInterface;
 use crate::core::pbrt::Float;
 use crate::shapes::Shape;
-use crate::utils::spectrum::Spectrum;
+use crate::spectra::Spectrum;
 use crate::utils::transform::Transform;
 
 use std::sync::Arc;

src/shapes/bilinear.rs

@@ -61,6 +61,7 @@ impl BilinearPatchShape {
         } else {
             const NA: usize = 3;
             let mut p = [[Point3f::default(); NA + 1]; NA + 1];
+            #[allow(clippy::needless_range_loop)]
             for i in 0..=NA {
                 let u = i as Float / NA as Float;
                 for j in 0..=NA {
@@ -94,10 +95,10 @@ impl BilinearPatchShape {
         let mesh = self.mesh();
         let start_index = 4 * self.blp_index;
         let v = &mesh.vertex_indices[start_index..start_index + 4];
-        let p00: Point3f = mesh.p[v[0] as usize];
+        let p00: Point3f = mesh.p[v[0]];
-        let p10: Point3f = mesh.p[v[1] as usize];
+        let p10: Point3f = mesh.p[v[1]];
-        let p01: Point3f = mesh.p[v[2] as usize];
+        let p01: Point3f = mesh.p[v[2]];
-        let p11: Point3f = mesh.p[v[3] as usize];
+        let p11: Point3f = mesh.p[v[3]];
         let n = mesh
             .n
             .as_ref()
@@ -141,7 +142,7 @@ impl BilinearPatchShape {
                 return false;
             }
         }
-        return true;
+        true
     }
 
     fn intersect_bilinear_patch(
@@ -628,13 +629,13 @@ impl BilinearPatchShape {
         }
     }
 
-    pub fn pdf(&self, intr: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, intr: &dyn Interaction) -> Float {
         let Some(si) = intr.as_any().downcast_ref::<SurfaceInteraction>() else {
             return 0.;
         };
         let data = self.get_data();
         let uv = if let Some(uvs) = &data.mesh.uv {
-            Point2f::invert_bilinear(si.uv, &uvs)
+            Point2f::invert_bilinear(si.uv, uvs)
         } else {
             si.uv
         };
@@ -675,18 +676,14 @@ impl BilinearPatchShape {
             || data.mesh.image_distribution.is_some()
             || spherical_quad_area(v00, v10, v01, v11) <= Self::MIN_SPHERICAL_SAMPLE_AREA;
         if use_area_sampling {
-            let isect_pdf = self.pdf(Arc::new(isect.intr.as_ref()));
+            let isect_pdf = self.pdf(isect.intr.as_ref());
             let distsq = ctx.p().distance_squared(isect.intr.p());
             let absdot = Vector3f::from(isect.intr.n()).abs_dot(-wi);
             if absdot == 0. {
                 return 0.;
             }
             let pdf = isect_pdf * distsq / absdot;
-            if pdf.is_infinite() {
+            if pdf.is_infinite() { 0. } else { pdf }
-                return 0.;
-            } else {
-                return pdf;
-            }
         } else {
             let mut pdf = 1. / spherical_quad_area(v00, v10, v01, v11);
             if ctx.ns != Normal3f::zero() {

src/shapes/curves.rs

@@ -76,7 +76,7 @@ impl CurveShape {
         };
 
         let context = IntersectionContext {
-            ray: ray,
+            ray,
             object_from_ray: Arc::new(ray_from_object.inverse()),
             common: self.common.clone(),
         };
@@ -208,6 +208,7 @@ impl CurveShape {
         true
     }
 
+    #[allow(clippy::too_many_arguments)]
     fn intersection_result(
         &self,
         context: &IntersectionContext,
@@ -257,7 +258,10 @@ impl CurveShape {
             flip_normal,
         );
 
-        Some(ShapeIntersection { intr: Box::new(intr), t_hit })
+        Some(ShapeIntersection {
+            intr: Box::new(intr),
+            t_hit,
+        })
     }
 
     pub fn bounds(&self) -> Bounds3f {
@@ -296,7 +300,7 @@ impl CurveShape {
         self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY))
     }
 
-    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: &dyn Interaction) -> Float {
         todo!()
     }
 
@@ -308,7 +312,11 @@ impl CurveShape {
         todo!()
     }
 
-    pub fn sample_from_context(&self, _ctx: &ShapeSampleContext, _u: Point2f) -> Option<ShapeSample> {
+    pub fn sample_from_context(
+        &self,
+        _ctx: &ShapeSampleContext,
+        _u: Point2f,
+    ) -> Option<ShapeSample> {
         todo!()
     }
 }

src/shapes/cylinder.rs

@@ -13,7 +13,7 @@ use std::sync::Arc;
 impl CylinderShape {
     pub fn new(
         render_from_object: Arc<Transform<Float>>,
-        object_from_render:  Arc<Transform<Float>>,
+        object_from_render: Arc<Transform<Float>>,
         reverse_orientation: bool,
         radius: Float,
         z_min: Float,
@@ -55,12 +55,11 @@ impl CylinderShape {
             return None;
         }
         let root_discrim = discrim.sqrt();
-        let q: Interval;
+        let q = if Float::from(b) < 0. {
-        if Float::from(b) < 0. {
+            -0.5 * (b - root_discrim)
-            q = -0.5 * (b - root_discrim);
         } else {
-            q = -0.5 * (b + root_discrim);
+            -0.5 * (b + root_discrim)
-        }
+        };
         let mut t0 = q / a;
         let mut t1 = c / q;
         if t0.low > t1.low {
@@ -155,7 +154,7 @@ impl CylinderShape {
         let flip_normal = self.reverse_orientation ^ self.transform_swap_handedness;
         let wo_object = self.object_from_render.apply_to_vector(wo);
         // (*renderFromObject)
-        let surf_point = SurfaceInteraction::new(
+        SurfaceInteraction::new(
             Point3fi::new_with_error(p_hit, p_error),
             Point2f::new(u, v),
             wo_object,
@@ -165,8 +164,7 @@ impl CylinderShape {
             dndv,
             time,
             flip_normal,
-        );
+        )
-        surf_point
     }
 
     pub fn area(&self) -> Float {
@@ -189,7 +187,7 @@ impl CylinderShape {
         let t = t_max.unwrap_or(Float::INFINITY);
         if let Some(isect) = self.basic_intersect(ray, t) {
             let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time);
-            return Some(ShapeIntersection::new(intr, isect.t_hit));
+            Some(ShapeIntersection::new(intr, isect.t_hit))
         } else {
             None
         }
@@ -203,7 +201,7 @@ impl CylinderShape {
         }
     }
 
-    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: &dyn Interaction) -> Float {
         1. / self.area()
     }
 
@@ -216,9 +214,9 @@ impl CylinderShape {
             if pdf.is_infinite() {
                 return 0.;
             }
-            return pdf;
+            pdf
         } else {
-            return 0.;
+            0.
         }
     }
 
@@ -264,6 +262,6 @@ impl CylinderShape {
         if ss.pdf.is_infinite() {
             return None;
         }
-        return Some(ss);
+        Some(ss)
     }
 }

src/shapes/disk.rs

@@ -88,7 +88,7 @@ impl DiskShape {
         let p_error = Vector3f::zero();
         let flip_normal = self.reverse_orientation ^ self.transform_swap_handedness;
         let wo_object = self.object_from_render.apply_to_vector(wo);
-        let surf_point = SurfaceInteraction::new(
+        SurfaceInteraction::new(
             Point3fi::new_with_error(p_hit, p_error),
             Point2f::new(u, v),
             wo_object,
@@ -98,8 +98,7 @@ impl DiskShape {
             dndv,
             time,
             flip_normal,
-        );
+        )
-        surf_point
     }
 
     pub fn area(&self) -> Float {
@@ -127,7 +126,7 @@ impl DiskShape {
         let t = t_max.unwrap_or(Float::INFINITY);
         if let Some(isect) = self.basic_intersect(ray, t) {
             let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time);
-            return Some(ShapeIntersection::new(intr, isect.t_hit));
+            Some(ShapeIntersection::new(intr, isect.t_hit))
         } else {
             None
         }
@@ -182,10 +181,10 @@ impl DiskShape {
         if ss.pdf.is_infinite() {
             return None;
         }
-        return Some(ss);
+        Some(ss)
     }
 
-    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: &dyn Interaction) -> Float {
         1. / self.area()
     }
 
@@ -198,9 +197,9 @@ impl DiskShape {
             if pdf.is_infinite() {
                 return 0.;
             }
-            return pdf;
+            pdf
         } else {
-            return 0.;
+            0.
         }
     }
 }

src/shapes/mod.rs

@@ -92,6 +92,7 @@ pub struct CurveCommon {
 }
 
 impl CurveCommon {
+    #[allow(clippy::too_many_arguments)]
     pub fn new(
         c: &[Point3f],
         w0: Float,
@@ -105,10 +106,7 @@ impl CurveCommon {
         let transform_swap_handedness = render_from_object.swaps_handedness();
         let width = [w0, w1];
         assert_eq!(c.len(), 4);
-        let mut cp_obj = [Point3f::default(); 4];
+        let cp_obj: [Point3f; 4] = c[..4].try_into().unwrap();
-        for i in 0..4 {
-            cp_obj[i] = c[i];
-        }
 
         let mut n = [Normal3f::default(); 2];
         let mut normal_angle: Float = 0.;
@@ -166,7 +164,7 @@ impl ShapeIntersection {
     }
 
     pub fn intr_mut(&mut self) -> &mut SurfaceInteraction {
-        &mut *self.intr
+        &mut self.intr
     }
 
     pub fn t_hit(&self) -> Float {
@@ -361,7 +359,7 @@ impl Shape {
         }
     }
 
-    pub fn pdf(&self, interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, interaction: &dyn Interaction) -> Float {
         match self {
             Shape::None => 0.,
             Shape::Sphere(s) => s.pdf(interaction),

src/shapes/sphere.rs

@@ -1,7 +1,7 @@
 use super::{
     Bounds3f, DirectionCone, Float, Interaction, Normal3f, PI, Point2f, Point3f, Point3fi,
-    QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext,
+    QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, SphereShape,
-    SphereShape, SurfaceInteraction, Transform, Vector3f, Vector3fi,
+    SurfaceInteraction, Transform, Vector3f, Vector3fi,
 };
 use crate::core::pbrt::{clamp_t, gamma};
 use crate::geometry::{Frame, Sqrt, VectorLike, spherical_direction};
@@ -61,13 +61,12 @@ impl SphereShape {
         }
 
         let root_discrim = discrim.sqrt();
-        let q: Interval;
 
-        if Float::from(b) < 0. {
+        let q = if Float::from(b) < 0. {
-            q = -0.5 * (b - root_discrim);
+            -0.5 * (b - root_discrim)
         } else {
-            q = -0.5 * (b + root_discrim);
+            -0.5 * (b + root_discrim)
-        }
+        };
 
         let mut t0 = q / a;
         let mut t1 = c / q;
@@ -179,7 +178,7 @@ impl SphereShape {
         let p_error = gamma(5) * Vector3f::from(p_hit).abs();
         let flip_normal = self.reverse_orientation ^ self.transform_swap_handedness;
         let wo_object = self.object_from_render.apply_to_vector(wo);
-        let surf_point = SurfaceInteraction::new(
+        SurfaceInteraction::new(
             Point3fi::new_with_error(p_hit, p_error),
             Point2f::new(u, v),
             wo_object,
@@ -189,9 +188,7 @@ impl SphereShape {
             dndv,
             time,
             flip_normal,
-        );
+        )
-        // self.render_from_object.apply_to_point(surf_point)
-        surf_point
     }
 
     pub fn bounds(&self) -> Bounds3f {
@@ -210,7 +207,7 @@ impl SphereShape {
         self.phi_max * self.radius * (self.z_max - self.z_min)
     }
 
-    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: &dyn Interaction) -> Float {
         1. / self.area()
     }
 
@@ -218,7 +215,7 @@ impl SphereShape {
         let t = t_max.unwrap_or(Float::INFINITY);
         if let Some(isect) = self.basic_intersect(ray, t) {
             let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time);
-            return Some(ShapeIntersection::new(intr, isect.t_hit));
+            Some(ShapeIntersection::new(intr, isect.t_hit))
         } else {
             None
         }

src/shapes/triangle.rs

@@ -1,7 +1,7 @@
 use super::{
     Bounds3f, DirectionCone, Float, Interaction, Normal3f, Point2f, Point3f, Point3fi, Ray,
-    ShapeIntersection, ShapeSample, ShapeSampleContext, SurfaceInteraction,
+    ShapeIntersection, ShapeSample, ShapeSampleContext, SurfaceInteraction, TriangleIntersection,
-    TriangleIntersection, TriangleShape, Vector2f, Vector3f,
+    TriangleShape, Vector2f, Vector3f,
 };
 use crate::core::pbrt::gamma;
 use crate::geometry::{Sqrt, Tuple, VectorLike, spherical_triangle_area};
@@ -74,7 +74,7 @@ impl TriangleShape {
     fn solid_angle(&self, p: Point3f) -> Float {
         let data = self.get_data();
         let [p0, p1, p2] = data.vertices;
-        spherical_triangle_area::<Float>(
+        spherical_triangle_area(
             (p0 - p).normalize(),
             (p1 - p).normalize(),
             (p2 - p).normalize(),
@@ -296,7 +296,7 @@ impl TriangleShape {
 
         let mut ts = Vector3f::from(ns).cross(ss);
         if ts.norm_squared() > 0. {
-            ss = ts.cross(Vector3f::from(ns)).into();
+            ss = ts.cross(Vector3f::from(ns));
         } else {
             (ss, ts) = Vector3f::from(ns).coordinate_system();
         }
@@ -319,7 +319,7 @@ impl TriangleShape {
 
         isect.shading.n = ns;
         isect.shading.dpdu = ss;
-        isect.shading.dpdv = ts.into();
+        isect.shading.dpdv = ts;
         isect.dndu = dndu;
         isect.dndv = dndv;
     }
@@ -344,14 +344,13 @@ impl TriangleShape {
         self.get_data().area
     }
 
-    pub fn pdf(&self, _interaction: Arc<&dyn Interaction>) -> Float {
+    pub fn pdf(&self, _interaction: &dyn Interaction) -> Float {
         1. / self.area()
     }
 
     pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float {
         let solid_angle = self.solid_angle(ctx.p());
-        if solid_angle < Self::MIN_SPHERICAL_SAMPLE_AREA
+        if (Self::MIN_SPHERICAL_SAMPLE_AREA..Self::MAX_SPHERICAL_SAMPLE_AREA).contains(&solid_angle)
-            || solid_angle > Self::MAX_SPHERICAL_SAMPLE_AREA
         {
             let ray = ctx.spawn_ray(wi);
             return self.intersect(&ray, None).map_or(0., |isect| {
@@ -389,27 +388,23 @@ impl TriangleShape {
         let data = self.get_data();
         let [p0, p1, p2] = data.vertices;
         let solid_angle = self.solid_angle(ctx.p());
-        if solid_angle < Self::MIN_SPHERICAL_SAMPLE_AREA
+        if (Self::MIN_SPHERICAL_SAMPLE_AREA..Self::MAX_SPHERICAL_SAMPLE_AREA).contains(&solid_angle)
-            || solid_angle > Self::MAX_SPHERICAL_SAMPLE_AREA
         {
             // Sample shape by area and compute incident direction wi
-            return self
+            return self.sample(u).and_then(|mut ss| {
-                .sample(u)
+                let mut intr_clone = (*ss.intr).clone();
-                .map(|mut ss| {
+                intr_clone.common.time = ctx.time;
-                    let mut intr_clone = (*ss.intr).clone();
+                ss.intr = Arc::new(intr_clone);
-                    intr_clone.common.time = ctx.time;
-                    ss.intr = Arc::new(intr_clone);
 
-                    let wi = (ss.intr.p() - ctx.p()).normalize();
+                let wi = (ss.intr.p() - ctx.p()).normalize();
-                    if wi.norm_squared() == 0. {
+                if wi.norm_squared() == 0. {
-                        return None;
+                    return None;
-                    }
+                }
-                    let absdot = Vector3f::from(ss.intr.n()).abs_dot(-wi);
+                let absdot = Vector3f::from(ss.intr.n()).abs_dot(-wi);
-                    let d2 = ctx.p().distance_squared(ss.intr.p());
+                let d2 = ctx.p().distance_squared(ss.intr.p());
-                    ss.pdf /= absdot / d2;
+                ss.pdf /= absdot / d2;
-                    if ss.pdf.is_infinite() { None } else { Some(ss) }
+                if ss.pdf.is_infinite() { None } else { Some(ss) }
-                })
+            });
-                .flatten();
         }
 
         // Sample spherical triangle from reference point
@@ -432,9 +427,7 @@ impl TriangleShape {
             pdf = bilinear_pdf(u, &w);
         }
 
-        let Some((b, tri_pdf)) = sample_spherical_triangle(&[p0, p1, p2], ctx.p(), u) else {
+        let (b, tri_pdf) = sample_spherical_triangle(&[p0, p1, p2], ctx.p(), u)?;
-            return None;
-        };
         if tri_pdf == 0. {
             return None;
         }
@@ -505,7 +498,10 @@ impl TriangleShape {
         self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY))
             .map(|ti| {
                 let intr = self.interaction_from_intersection(ti, ray.time, -ray.d);
-                ShapeIntersection { intr: Box::new(intr), t_hit: ti.t }
+                ShapeIntersection {
+                    intr: Box::new(intr),
+                    t_hit: ti.t,
+                }
             })
     }
 

src/spectra/data.rs

@@ -0,0 +1,33 @@
+use super::Spectrum;
+use super::sampled::{LAMBDA_MAX, LAMBDA_MIN};
+use super::simple::{DenselySampledSpectrum, PiecewiseLinearSpectrum};
+use crate::core::cie;
+use crate::core::pbrt::Float;
+use once_cell::sync::Lazy;
+
+fn create_cie_spectrum(data: &[Float]) -> Spectrum {
+    let pls = PiecewiseLinearSpectrum::from_interleaved(data);
+
+    let dss = DenselySampledSpectrum::from_spectrum(
+        &Spectrum::PiecewiseLinear(pls),
+        LAMBDA_MIN,
+        LAMBDA_MAX,
+    );
+
+    Spectrum::DenselySampled(dss)
+}
+
+pub(crate) fn cie_x() -> &'static Spectrum {
+    static X: Lazy<Spectrum> = Lazy::new(|| create_cie_spectrum(&cie::CIE_X));
+    &X
+}
+
+pub(crate) fn cie_y() -> &'static Spectrum {
+    static Y: Lazy<Spectrum> = Lazy::new(|| create_cie_spectrum(&cie::CIE_Y));
+    &Y
+}
+
+pub(crate) fn cie_z() -> &'static Spectrum {
+    static Z: Lazy<Spectrum> = Lazy::new(|| create_cie_spectrum(&cie::CIE_Z));
+    &Z
+}

src/spectra/mod.rs

@@ -0,0 +1,68 @@
+pub mod data;
+pub mod rgb;
+pub mod sampled;
+pub mod simple;
+
+use crate::core::pbrt::Float;
+use crate::utils::color::{RGB, XYZ};
+use crate::utils::colorspace::RGBColorSpace;
+use enum_dispatch::enum_dispatch;
+
+pub use data::*;
+pub use rgb::*;
+use sampled::{CIE_Y_INTEGRAL, LAMBDA_MAX, LAMBDA_MIN};
+pub use sampled::{N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths};
+pub use simple::*; // CIE_X, etc
+//
+#[enum_dispatch]
+pub trait SpectrumTrait {
+    fn evaluate(&self, lambda: Float) -> Float;
+    fn max_value(&self) -> Float;
+}
+
+#[enum_dispatch(SpectrumTrait)]
+#[derive(Debug, Clone)]
+pub enum Spectrum {
+    Constant(ConstantSpectrum),
+    DenselySampled(DenselySampledSpectrum),
+    PiecewiseLinear(PiecewiseLinearSpectrum),
+    Blackbody(BlackbodySpectrum),
+    RGBAlbedo(RGBAlbedoSpectrum),
+}
+
+impl Spectrum {
+    pub fn std_illuminant_d65() -> Self {
+        todo!()
+    }
+
+    pub fn to_xyz(&self) -> XYZ {
+        let x = self.inner_product(data::cie_x());
+        let y = self.inner_product(data::cie_y());
+        let z = self.inner_product(data::cie_z());
+
+        XYZ::new(x, y, z) / CIE_Y_INTEGRAL
+    }
+
+    fn to_rgb(&self, cs: &RGBColorSpace) -> RGB {
+        let xyz = self.to_xyz();
+        cs.to_rgb(xyz)
+    }
+
+    pub fn sample(&self, wavelengths: &SampledWavelengths) -> SampledSpectrum {
+        let mut s = SampledSpectrum::default();
+        for i in 0..N_SPECTRUM_SAMPLES {
+            s[i] = self.evaluate(wavelengths[i]);
+        }
+        s
+    }
+
+    pub fn inner_product(&self, other: &Spectrum) -> Float {
+        let mut integral = 0.0;
+        // Iterate integer wavelengths.
+        for lambda in LAMBDA_MIN..=LAMBDA_MAX {
+            let l = lambda as Float;
+            integral += self.evaluate(l) * other.evaluate(l);
+        }
+        integral
+    }
+}

src/spectra/rgb.rs

@@ -0,0 +1,172 @@
+use super::sampled::{
+    LAMBDA_MAX, LAMBDA_MIN, N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths,
+};
+use crate::core::pbrt::Float;
+use crate::spectra::{DenselySampledSpectrum, SpectrumTrait};
+use crate::utils::color::{RGB, RGBSigmoidPolynomial, XYZ};
+use crate::utils::colorspace::RGBColorSpace;
+use std::sync::Arc;
+
+#[derive(Debug, Clone, Copy)]
+pub struct RGBAlbedoSpectrum {
+    rsp: RGBSigmoidPolynomial,
+}
+
+impl RGBAlbedoSpectrum {
+    pub fn new(cs: &RGBColorSpace, rgb: RGB) -> Self {
+        Self {
+            rsp: cs.to_rgb_coeffs(rgb),
+        }
+    }
+
+    fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let mut s = SampledSpectrum::default();
+        for i in 0..N_SPECTRUM_SAMPLES {
+            s[i] = self.rsp.evaluate(lambda[i]);
+        }
+        s
+    }
+}
+
+impl SpectrumTrait for RGBAlbedoSpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        self.rsp.evaluate(lambda)
+    }
+
+    fn max_value(&self) -> Float {
+        self.rsp.max_value()
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct UnboundedRGBSpectrum {
+    scale: Float,
+    rsp: RGBSigmoidPolynomial,
+}
+
+impl UnboundedRGBSpectrum {
+    pub fn new(cs: RGBColorSpace, rgb: RGB) -> Self {
+        let m = rgb.r.max(rgb.g).max(rgb.b);
+        let scale = 2.0 * m;
+        let scaled_rgb = if scale != 0.0 {
+            rgb / scale
+        } else {
+            RGB::new(0.0, 0.0, 0.0)
+        };
+        Self {
+            scale,
+            rsp: cs.to_rgb_coeffs(scaled_rgb),
+        }
+    }
+}
+
+impl SpectrumTrait for UnboundedRGBSpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        self.scale * self.rsp.evaluate(lambda)
+    }
+
+    fn max_value(&self) -> Float {
+        self.scale * self.rsp.max_value()
+    }
+}
+
+#[derive(Debug, Clone, Default)]
+pub struct RGBIlluminantSpectrum {
+    scale: Float,
+    rsp: RGBSigmoidPolynomial,
+    illuminant: Option<Arc<DenselySampledSpectrum>>,
+}
+
+impl RGBIlluminantSpectrum {
+    pub fn new(cs: RGBColorSpace, rgb: RGB) -> Self {
+        let illuminant = &cs.illuminant;
+        let densely_sampled =
+            DenselySampledSpectrum::from_spectrum(illuminant, LAMBDA_MIN, LAMBDA_MAX);
+        let m = rgb.max();
+        let scale = 2. * m;
+        let rsp = cs.to_rgb_coeffs(if scale == 1. {
+            rgb / scale
+        } else {
+            RGB::new(0., 0., 0.)
+        });
+        Self {
+            scale,
+            rsp,
+            illuminant: Some(Arc::new(densely_sampled)),
+        }
+    }
+}
+
+impl SpectrumTrait for RGBIlluminantSpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        match &self.illuminant {
+            Some(illuminant) => {
+                self.scale * self.rsp.evaluate(lambda) * illuminant.evaluate(lambda)
+            }
+            None => 0.0,
+        }
+    }
+
+    fn max_value(&self) -> Float {
+        match &self.illuminant {
+            Some(illuminant) => self.scale * self.rsp.max_value() * illuminant.max_value(),
+            None => 0.0,
+        }
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct RGBSpectrum {
+    pub c: [Float; 3],
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct RGBUnboundedSpectrum {
+    scale: Float,
+    rsp: RGBSigmoidPolynomial,
+}
+
+impl Default for RGBUnboundedSpectrum {
+    fn default() -> Self {
+        Self {
+            scale: 0.0,
+            rsp: RGBSigmoidPolynomial::default(),
+        }
+    }
+}
+
+impl RGBUnboundedSpectrum {
+    pub fn new(cs: &RGBColorSpace, rgb: RGB) -> Self {
+        let m = rgb.max();
+
+        let scale = 2.0 * m;
+
+        let rgb_norm = if scale > 0.0 {
+            rgb / scale
+        } else {
+            RGB::new(0.0, 0.0, 0.0)
+        };
+
+        let rsp = cs.to_rgb_coeffs(rgb_norm);
+
+        Self { scale, rsp }
+    }
+
+    pub fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let mut s = SampledSpectrum::default();
+        for i in 0..N_SPECTRUM_SAMPLES {
+            s[i] = self.scale * self.rsp.evaluate(lambda[i]);
+        }
+        s
+    }
+}
+
+impl SpectrumTrait for RGBUnboundedSpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        self.scale * self.rsp.evaluate(lambda)
+    }
+
+    fn max_value(&self) -> Float {
+        self.scale * self.rsp.max_value()
+    }
+}

src/spectra/sampled.rs

@@ -0,0 +1,359 @@
+use crate::core::pbrt::Float;
+use std::ops::{
+    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
+};
+
+pub const CIE_Y_INTEGRAL: Float = 106.856895;
+
+pub const N_SPECTRUM_SAMPLES: usize = 1200;
+pub const LAMBDA_MIN: i32 = 360;
+pub const LAMBDA_MAX: i32 = 830;
+
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct SampledSpectrum {
+    pub values: [Float; N_SPECTRUM_SAMPLES],
+}
+
+impl Default for SampledSpectrum {
+    fn default() -> Self {
+        Self {
+            values: [0.0; N_SPECTRUM_SAMPLES],
+        }
+    }
+}
+
+impl SampledSpectrum {
+    pub fn new(c: Float) -> Self {
+        Self {
+            values: [c; N_SPECTRUM_SAMPLES],
+        }
+    }
+
+    #[inline]
+    pub fn from_fn<F>(cb: F) -> Self
+    where
+        F: FnMut(usize) -> Float,
+    {
+        Self {
+            values: std::array::from_fn(cb),
+        }
+    }
+
+    pub fn from_vector(v: Vec<Float>) -> Self {
+        let mut values = [0.0; N_SPECTRUM_SAMPLES];
+        let count = v.len().min(N_SPECTRUM_SAMPLES);
+        values[..count].copy_from_slice(&v[..count]);
+        Self { values }
+    }
+
+    pub fn is_black(&self) -> bool {
+        self.values.iter().all(|&sample| sample == 0.0)
+    }
+
+    pub fn has_nans(&self) -> bool {
+        self.values.iter().any(|&v| v.is_nan())
+    }
+
+    pub fn min_component_value(&self) -> Float {
+        self.values.iter().fold(Float::INFINITY, |a, &b| a.min(b))
+    }
+
+    pub fn max_component_value(&self) -> Float {
+        self.values
+            .iter()
+            .fold(Float::NEG_INFINITY, |a, &b| a.max(b))
+    }
+
+    pub fn average(&self) -> Float {
+        self.values.iter().sum::<Float>() / (N_SPECTRUM_SAMPLES as Float)
+    }
+
+    pub fn safe_div(&self, rhs: SampledSpectrum) -> Self {
+        let mut r = SampledSpectrum::new(0.0);
+        for i in 0..N_SPECTRUM_SAMPLES {
+            r.values[i] = if rhs[i] != 0.0 {
+                self.values[i] / rhs.values[i]
+            } else {
+                0.0
+            }
+        }
+        r
+    }
+
+    pub fn exp(&self) -> SampledSpectrum {
+        let values = self.values.map(|v| v.exp());
+        let ret = Self { values };
+        debug_assert!(!ret.has_nans());
+        ret
+    }
+
+    pub fn pow_int(&self, mut power: usize) -> Self {
+        let mut result = Self::new(1.0);
+        let mut base = *self;
+        while power > 0 {
+            if power % 2 == 1 {
+                result *= base;
+            }
+            base *= base;
+            power /= 2;
+        }
+        result
+    }
+}
+
+impl<'a> IntoIterator for &'a SampledSpectrum {
+    type Item = &'a Float;
+    type IntoIter = std::slice::Iter<'a, Float>;
+
+    fn into_iter(self) -> Self::IntoIter {
+        self.values.iter()
+    }
+}
+
+impl<'a> IntoIterator for &'a mut SampledSpectrum {
+    type Item = &'a mut Float;
+    type IntoIter = std::slice::IterMut<'a, Float>;
+
+    fn into_iter(self) -> Self::IntoIter {
+        self.values.iter_mut()
+    }
+}
+
+impl Index<usize> for SampledSpectrum {
+    type Output = Float;
+    fn index(&self, i: usize) -> &Self::Output {
+        &self.values[i]
+    }
+}
+
+impl IndexMut<usize> for SampledSpectrum {
+    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
+        &mut self.values[i]
+    }
+}
+
+impl Add for SampledSpectrum {
+    type Output = Self;
+    fn add(self, rhs: Self) -> Self {
+        let mut ret = self;
+        for i in 0..N_SPECTRUM_SAMPLES {
+            ret.values[i] += rhs.values[i];
+        }
+        ret
+    }
+}
+
+impl AddAssign for SampledSpectrum {
+    fn add_assign(&mut self, rhs: Self) {
+        for i in 0..N_SPECTRUM_SAMPLES {
+            self.values[i] += rhs.values[i];
+        }
+    }
+}
+
+impl Sub for SampledSpectrum {
+    type Output = Self;
+    fn sub(self, rhs: Self) -> Self {
+        let mut ret = self;
+        for i in 0..N_SPECTRUM_SAMPLES {
+            ret.values[i] -= rhs.values[i];
+        }
+        ret
+    }
+}
+
+impl SubAssign for SampledSpectrum {
+    fn sub_assign(&mut self, rhs: Self) {
+        for i in 0..N_SPECTRUM_SAMPLES {
+            self.values[i] -= rhs.values[i];
+        }
+    }
+}
+
+impl Sub<SampledSpectrum> for Float {
+    type Output = SampledSpectrum;
+    fn sub(self, rhs: SampledSpectrum) -> SampledSpectrum {
+        SampledSpectrum::from_fn(|i| self - rhs[i])
+    }
+}
+
+impl Mul for SampledSpectrum {
+    type Output = Self;
+    fn mul(self, rhs: Self) -> Self {
+        let mut ret = self;
+        for i in 0..N_SPECTRUM_SAMPLES {
+            ret.values[i] *= rhs.values[i];
+        }
+        ret
+    }
+}
+
+impl MulAssign for SampledSpectrum {
+    fn mul_assign(&mut self, rhs: Self) {
+        for i in 0..N_SPECTRUM_SAMPLES {
+            self.values[i] *= rhs.values[i];
+        }
+    }
+}
+
+impl Mul<Float> for SampledSpectrum {
+    type Output = Self;
+    fn mul(self, rhs: Float) -> Self {
+        let mut ret = self;
+        for i in 0..N_SPECTRUM_SAMPLES {
+            ret.values[i] *= rhs;
+        }
+        ret
+    }
+}
+
+impl Mul<SampledSpectrum> for Float {
+    type Output = SampledSpectrum;
+    fn mul(self, rhs: SampledSpectrum) -> SampledSpectrum {
+        rhs * self
+    }
+}
+
+impl MulAssign<Float> for SampledSpectrum {
+    fn mul_assign(&mut self, rhs: Float) {
+        for i in 0..N_SPECTRUM_SAMPLES {
+            self.values[i] *= rhs;
+        }
+    }
+}
+
+impl DivAssign for SampledSpectrum {
+    fn div_assign(&mut self, rhs: Self) {
+        for i in 0..N_SPECTRUM_SAMPLES {
+            debug_assert_ne!(0.0, rhs.values[i]);
+            self.values[i] /= rhs.values[i];
+        }
+    }
+}
+
+impl Div for SampledSpectrum {
+    type Output = Self;
+    fn div(self, rhs: Self) -> Self::Output {
+        let mut ret = self;
+        ret /= rhs;
+        ret
+    }
+}
+impl Div<Float> for SampledSpectrum {
+    type Output = Self;
+
+    fn div(self, rhs: Float) -> Self::Output {
+        debug_assert_ne!(rhs, 0.0);
+        let mut ret = self;
+        for i in 0..N_SPECTRUM_SAMPLES {
+            ret.values[i] /= rhs;
+        }
+        ret
+    }
+}
+
+impl DivAssign<Float> for SampledSpectrum {
+    fn div_assign(&mut self, rhs: Float) {
+        debug_assert_ne!(rhs, 0.0);
+        for i in 0..N_SPECTRUM_SAMPLES {
+            self.values[i] /= rhs;
+        }
+    }
+}
+
+impl Neg for SampledSpectrum {
+    type Output = Self;
+
+    fn neg(self) -> Self::Output {
+        let mut ret = SampledSpectrum::new(0.0);
+        for i in 0..N_SPECTRUM_SAMPLES {
+            ret.values[i] = -self.values[i];
+        }
+        ret
+    }
+}
+
+#[derive(Debug, Copy, Clone, PartialEq)]
+pub struct SampledWavelengths {
+    pub lambda: [Float; N_SPECTRUM_SAMPLES],
+    pub pdf: [Float; N_SPECTRUM_SAMPLES],
+}
+
+impl SampledWavelengths {
+    pub fn pdf(&self) -> SampledSpectrum {
+        SampledSpectrum::from_vector(self.pdf.to_vec())
+    }
+
+    pub fn secondary_terminated(&self) -> bool {
+        for i in 1..N_SPECTRUM_SAMPLES {
+            if self.pdf[i] != 0.0 {
+                return false;
+            }
+        }
+        true
+    }
+
+    pub fn terminate_secondary(&mut self) {
+        if !self.secondary_terminated() {
+            for i in 1..N_SPECTRUM_SAMPLES {
+                self.pdf[i] = 0.0;
+            }
+            self.pdf[0] /= N_SPECTRUM_SAMPLES as Float;
+        }
+    }
+
+    pub fn sample_uniform(u: Float, lambda_min: Float, lambda_max: Float) -> Self {
+        let mut lambda = [0.0; N_SPECTRUM_SAMPLES];
+        lambda[0] = crate::core::pbrt::lerp(u, lambda_min, lambda_min);
+        let delta = (lambda_max - lambda_min) / N_SPECTRUM_SAMPLES as Float;
+        for i in 1..N_SPECTRUM_SAMPLES {
+            lambda[i] = lambda[i - 1] + delta;
+            if lambda[i] > lambda_max {
+                lambda[i] = lambda_min + (lambda[i] - lambda_max);
+            }
+        }
+
+        let pdf = [1. / (lambda_max - lambda_min); N_SPECTRUM_SAMPLES];
+
+        Self { lambda, pdf }
+    }
+
+    pub fn sample_visible_wavelengths(u: Float) -> Float {
+        538.0 - 138.888889 * Float::atanh(0.85691062 - 1.82750197 * u)
+    }
+
+    pub fn visible_wavelengths_pdf(lambda: Float) -> Float {
+        if !(360.0..830.0).contains(&lambda) {
+            return 0.0;
+        }
+        0.0039398042 / (Float::cosh(0.0072 * (lambda - 538.0))).sqrt()
+    }
+
+    pub fn sample_visible(u: Float) -> Self {
+        let mut lambda = [0.0; N_SPECTRUM_SAMPLES];
+        let mut pdf = [0.0; N_SPECTRUM_SAMPLES];
+
+        for i in 0..N_SPECTRUM_SAMPLES {
+            let mut up = u + i as Float / N_SPECTRUM_SAMPLES as Float;
+            if up > 1.0 {
+                up -= 1.0;
+            }
+            lambda[i] = Self::sample_visible_wavelengths(up);
+            pdf[i] = Self::visible_wavelengths_pdf(lambda[i]);
+        }
+        Self { lambda, pdf }
+    }
+}
+
+impl Index<usize> for SampledWavelengths {
+    type Output = Float;
+    fn index(&self, i: usize) -> &Self::Output {
+        &self.lambda[i]
+    }
+}
+
+impl IndexMut<usize> for SampledWavelengths {
+    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
+        &mut self.lambda[i]
+    }
+}

src/spectra/simple.rs

@@ -0,0 +1,230 @@
+use super::sampled::{LAMBDA_MAX, LAMBDA_MIN};
+use crate::core::pbrt::Float;
+use crate::spectra::{
+    N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths, Spectrum, SpectrumTrait,
+};
+
+#[derive(Debug, Clone, Copy)]
+pub struct ConstantSpectrum {
+    c: Float,
+}
+
+impl ConstantSpectrum {
+    pub fn new(c: Float) -> Self {
+        Self { c }
+    }
+}
+
+impl SpectrumTrait for ConstantSpectrum {
+    fn evaluate(&self, _lambda: Float) -> Float {
+        self.c
+    }
+
+    fn max_value(&self) -> Float {
+        self.c
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct DenselySampledSpectrum {
+    lambda_min: i32,
+    lambda_max: i32,
+    values: Vec<Float>,
+}
+
+impl DenselySampledSpectrum {
+    pub fn new(lambda_min: i32, lambda_max: i32) -> Self {
+        let n_values = (lambda_max - lambda_min + 1).max(0) as usize;
+        Self {
+            lambda_min,
+            lambda_max,
+            values: vec![0.0; n_values],
+        }
+    }
+
+    pub fn from_spectrum(spec: &Spectrum, lambda_min: i32, lambda_max: i32) -> Self {
+        let mut s = Self::new(lambda_min, lambda_max);
+        if s.values.is_empty() {
+            return s;
+        }
+        for lambda in lambda_min..=lambda_max {
+            let index = (lambda - lambda_min) as usize;
+            s.values[index] = spec.evaluate(lambda as Float);
+        }
+        s
+    }
+
+    pub fn from_function<F>(f: F, lambda_min: i32, lambda_max: i32) -> Self
+    where
+        F: Fn(Float) -> Float,
+    {
+        let mut s = Self::new(lambda_min, lambda_max);
+        if s.values.is_empty() {
+            return s;
+        }
+        for lambda in lambda_min..=lambda_max {
+            let index = (lambda - lambda_min) as usize;
+            s.values[index] = f(lambda as Float);
+        }
+        s
+    }
+
+    pub fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let mut s = SampledSpectrum::default();
+
+        for i in 0..N_SPECTRUM_SAMPLES {
+            let offset = lambda[i].round() as usize - LAMBDA_MIN as usize;
+            if offset >= self.values.len() {
+                s[i] = 0.;
+            } else {
+                s[i] = self.values[offset];
+            }
+        }
+        s
+    }
+
+    pub fn min_component_value(&self) -> Float {
+        self.values.iter().fold(Float::INFINITY, |a, &b| a.min(b))
+    }
+
+    pub fn max_component_value(&self) -> Float {
+        self.values
+            .iter()
+            .fold(Float::NEG_INFINITY, |a, &b| a.max(b))
+    }
+
+    pub fn average(&self) -> Float {
+        self.values.iter().sum::<Float>() / (N_SPECTRUM_SAMPLES as Float)
+    }
+
+    pub fn safe_div(&self, rhs: SampledSpectrum) -> Self {
+        let mut r = Self::new(1, 1);
+        for i in 0..N_SPECTRUM_SAMPLES {
+            r.values[i] = if rhs[i] != 0.0 {
+                self.values[i] / rhs.values[i]
+            } else {
+                0.0
+            }
+        }
+        r
+    }
+}
+
+impl SpectrumTrait for DenselySampledSpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        let offset = (lambda.round() as i32) - self.lambda_min;
+        if offset < 0 || offset as usize >= self.values.len() {
+            0.0
+        } else {
+            self.values[offset as usize]
+        }
+    }
+
+    fn max_value(&self) -> Float {
+        self.values.iter().fold(Float::MIN, |a, b| a.max(*b))
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct PiecewiseLinearSpectrum {
+    samples: Vec<(Float, Float)>,
+}
+
+impl PiecewiseLinearSpectrum {
+    pub fn from_interleaved(data: &[Float]) -> Self {
+        assert!(
+            data.len().is_multiple_of(2),
+            "Interleaved data must have an even number of elements"
+        );
+        let mut samples = Vec::new();
+        for pair in data.chunks(2) {
+            samples.push((pair[0], pair[1]));
+        }
+        samples.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
+        Self { samples }
+    }
+}
+
+impl SpectrumTrait for PiecewiseLinearSpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        if self.samples.is_empty() {
+            return 0.0;
+        }
+        // Handle boundary conditions
+        if lambda <= self.samples[0].0 {
+            return self.samples[0].1;
+        }
+        if lambda >= self.samples.last().unwrap().0 {
+            return self.samples.last().unwrap().1;
+        }
+
+        let i = self.samples.partition_point(|s| s.0 < lambda);
+        let s1 = self.samples[i - 1];
+        let s2 = self.samples[i];
+
+        let t = (lambda - s1.0) / (s2.0 - s1.0);
+        (1.0 - t) * s1.1 + t * s2.1
+    }
+
+    fn max_value(&self) -> Float {
+        if self.samples.is_empty() {
+            return 0.0;
+        }
+        self.samples
+            .iter()
+            .map(|(_, value)| *value)
+            .fold(0.0, |a, b| a.max(b))
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct BlackbodySpectrum {
+    temperature: Float,
+    normalization_factor: Float,
+}
+
+// Planck's Law
+impl BlackbodySpectrum {
+    const C: Float = 299792458.0;
+    const H: Float = 6.62606957e-34;
+    const KB: Float = 1.3806488e-23;
+
+    pub fn new(temperature: Float) -> Self {
+        // Physical constants
+        let lambda_max = 2.8977721e-3 / temperature * 1e9;
+        let max_val = Self::planck_law(lambda_max, temperature);
+        Self {
+            temperature,
+            normalization_factor: if max_val > 0.0 { 1.0 / max_val } else { 0.0 },
+        }
+    }
+
+    fn planck_law(lambda_nm: Float, temp: Float) -> Float {
+        if temp <= 0.0 {
+            return 0.0;
+        }
+        let lambda_m = lambda_nm * 1e-9;
+        let c1 = 2.0 * Self::H * Self::C * Self::C;
+        let c2 = (Self::H * Self::C) / Self::KB;
+
+        let numerator = c1 / lambda_m.powi(5);
+        let denominator = (c2 / (lambda_m * temp)).exp() - 1.0;
+
+        if denominator.is_infinite() {
+            0.0
+        } else {
+            numerator / denominator
+        }
+    }
+}
+
+impl SpectrumTrait for BlackbodySpectrum {
+    fn evaluate(&self, lambda: Float) -> Float {
+        Self::planck_law(lambda, self.temperature) * self.normalization_factor
+    }
+
+    fn max_value(&self) -> Float {
+        let lambda_max = 2.8977721e-3 / self.temperature * 1e9;
+        Self::planck_law(lambda_max, self.temperature)
+    }
+}

src/utils/color.rs

@@ -4,12 +4,14 @@ use std::ops::{
     Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
 };
 
-use super::spectrum::Spectrum;
 use crate::core::pbrt::{Float, lerp};
 use crate::geometry::Point2f;
+use crate::spectra::Spectrum;
 use crate::utils::math::SquareMatrix;
 use once_cell::sync::Lazy;
 
+use enum_dispatch::enum_dispatch;
+
 pub trait Triplet {
     fn from_triplet(c1: Float, c2: Float, c3: Float) -> Self;
 }
@@ -250,7 +252,7 @@ impl fmt::Display for XYZ {
     }
 }
 
-#[derive(Debug, Clone)]
+#[derive(Debug, Copy, Clone)]
 pub struct RGB {
     pub r: Float,
     pub g: Float,
@@ -453,6 +455,52 @@ impl fmt::Display for RGB {
     }
 }
 
+impl Mul<XYZ> for SquareMatrix<Float, 3> {
+    type Output = RGB;
+
+    fn mul(self, v: XYZ) -> RGB {
+        let r = self[0][0] * v.x + self[0][1] * v.y + self[0][2] * v.z;
+        let g = self[1][0] * v.x + self[1][1] * v.y + self[1][2] * v.z;
+        let b = self[2][0] * v.x + self[2][1] * v.y + self[2][2] * v.z;
+        RGB::new(r, g, b)
+    }
+}
+
+impl Mul<RGB> for SquareMatrix<Float, 3> {
+    type Output = XYZ;
+    fn mul(self, v: RGB) -> XYZ {
+        let x = self[0][0] * v.r + self[0][1] * v.g + self[0][2] * v.b;
+        let y = self[1][0] * v.r + self[1][1] * v.g + self[1][2] * v.b;
+        let z = self[2][0] * v.r + self[2][1] * v.g + self[2][2] * v.b;
+        XYZ::new(x, y, z)
+    }
+}
+
+pub trait MatrixMulColor {
+    fn mul_rgb(&self, v: RGB) -> RGB;
+    fn mul_xyz(&self, v: XYZ) -> XYZ;
+}
+
+impl MatrixMulColor for SquareMatrix<Float, 3> {
+    fn mul_rgb(&self, v: RGB) -> RGB {
+        let m = self;
+        RGB::new(
+            m[0][0] * v.r + m[0][1] * v.g + m[0][2] * v.b,
+            m[1][0] * v.r + m[1][1] * v.g + m[1][2] * v.b,
+            m[2][0] * v.r + m[2][1] * v.g + m[2][2] * v.b,
+        )
+    }
+
+    fn mul_xyz(&self, v: XYZ) -> XYZ {
+        let m = self;
+        XYZ::new(
+            m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z,
+            m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z,
+            m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z,
+        )
+    }
+}
+
 pub const RES: usize = 64;
 pub type CoefficientArray = [[[[[Float; 3]; RES]; RES]; RES]; 3];
 
@@ -497,7 +545,7 @@ impl RGBSigmoidPolynomial {
     pub fn max_value(&self) -> Float {
         let lambda = -self.c1 / (2.0 * self.c0);
         let result = self.evaluate(360.0).max(self.evaluate(830.0));
-        if lambda >= 360.0 && lambda <= 830.0 {
+        if (360.0..830.0).contains(&lambda) {
             return result.max(self.evaluate(lambda));
         }
         result
@@ -531,12 +579,10 @@ impl RGBToSpectrumTable {
             } else {
                 maxc = 2;
             }
+        } else if rgb[1] > rgb[2] {
+            maxc = 1;
         } else {
-            if rgb[1] > rgb[2] {
+            maxc = 2;
-                maxc = 1;
-            } else {
-                maxc = 2;
-            }
         }
 
         let z = rgb[maxc];
@@ -550,6 +596,7 @@ impl RGBToSpectrumTable {
         let dy = (y - yi) as usize;
         let dz = (z - self.z_nodes[zi]) / (self.z_nodes[zi + 1] - self.z_nodes[zi]);
         let mut c = [0.0; 3];
+        #[allow(clippy::needless_range_loop)]
         for i in 0..3 {
             let co = |dx: usize, dy: usize, dz: usize| {
                 self.coeffs[maxc][zi as usize + dz][yi as usize + dy][xi as usize + dx][i]
@@ -604,7 +651,8 @@ pub fn white_balance(src_white: Point2f, target_white: Point2f) -> SquareMatrix<
     XYZ_FROM_LMS * lms_correct * LMS_FROM_XYZ
 }
 
-pub trait ColorEncoding: 'static + Send + Sync + fmt::Debug + fmt::Display {
+#[enum_dispatch]
+pub trait ColorEncodingTrait: 'static + Send + Sync + fmt::Debug + fmt::Display {
     fn from_linear_slice(&self, vin: &[Float], vout: &mut [u8]);
     fn to_linear_slice(&self, vin: &[u8], vout: &mut [Float]);
     fn to_float_linear(&self, v: Float) -> Float;
@@ -613,9 +661,22 @@ pub trait ColorEncoding: 'static + Send + Sync + fmt::Debug + fmt::Display {
     }
 }
 
-#[derive(Debug)]
+#[enum_dispatch(ColorEncodingTrait)]
+#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)]
+pub enum ColorEncoding {
+    Linear(LinearEncoding),
+    SRGB(SRGBEncoding),
+}
+
+impl fmt::Display for ColorEncoding {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "Encoding")
+    }
+}
+
+#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)]
 pub struct LinearEncoding;
-impl ColorEncoding for LinearEncoding {
+impl ColorEncodingTrait for LinearEncoding {
     fn from_linear_slice(&self, vin: &[Float], vout: &mut [u8]) {
         for (i, &v) in vin.iter().enumerate() {
             vout[i] = (v.clamp(0.0, 1.0) * 255.0 + 0.5) as u8;
@@ -637,9 +698,9 @@ impl fmt::Display for LinearEncoding {
     }
 }
 
-#[derive(Debug)]
+#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)]
 pub struct SRGBEncoding;
-impl ColorEncoding for SRGBEncoding {
+impl ColorEncodingTrait for SRGBEncoding {
     fn from_linear_slice(&self, vin: &[Float], vout: &mut [u8]) {
         for (i, &v_linear) in vin.iter().enumerate() {
             let v = v_linear.clamp(0.0, 1.0);
@@ -674,8 +735,8 @@ impl fmt::Display for SRGBEncoding {
     }
 }
 
-pub static LINEAR: Lazy<&'static dyn ColorEncoding> = Lazy::new(|| &LinearEncoding);
+pub const LINEAR: ColorEncoding = ColorEncoding::Linear(LinearEncoding);
-pub static SRGB: Lazy<&'static dyn ColorEncoding> = Lazy::new(|| &SRGBEncoding);
+pub const SRGB: ColorEncoding = ColorEncoding::SRGB(SRGBEncoding);
 
 const SRGB_TO_LINEAR_LUT: [Float; 256] = [
     0.0000000000,

src/utils/colorspace.rs

@@ -1,8 +1,8 @@
 use super::color::{RGB, RGBSigmoidPolynomial, RGBToSpectrumTable, XYZ};
 use super::math::SquareMatrix;
-use super::spectrum::{DenselySampledSpectrum, SampledSpectrum, Spectrum};
 use crate::core::pbrt::Float;
 use crate::geometry::Point2f;
+use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, Spectrum};
 
 use once_cell::sync::Lazy;
 
@@ -16,7 +16,7 @@ pub enum ColorEncoding {
 }
 
 #[derive(Debug, Clone)]
-pub struct RGBColorspace {
+pub struct RGBColorSpace {
     pub r: Point2f,
     pub g: Point2f,
     pub b: Point2f,
@@ -27,7 +27,7 @@ pub struct RGBColorspace {
     pub rgb_from_xyz: SquareMatrix<Float, 3>,
 }
 
-impl RGBColorspace {
+impl RGBColorSpace {
     pub fn new(
         r: Point2f,
         g: Point2f,
@@ -35,7 +35,7 @@ impl RGBColorspace {
         illuminant: Spectrum,
         rgb_to_spectrum_table: RGBToSpectrumTable,
     ) -> Result<Self, Box<dyn Error>> {
-        let w_xyz = illuminant.to_xyz();
+        let w_xyz: XYZ = illuminant.to_xyz();
         let w = w_xyz.xy();
         let r_xyz = XYZ::from_xyy(r, Some(1.0));
         let g_xyz = XYZ::from_xyy(g, Some(1.0));
@@ -46,12 +46,11 @@ impl RGBColorspace {
             [r_xyz.z(), g_xyz.z(), g_xyz.z()],
         ];
         let rgb = SquareMatrix::new(rgb_values);
-        let c = rgb.inverse()? * w_xyz;
+        let c: RGB = rgb.inverse()? * w_xyz;
-        let xyz_from_rgb_m = [[c[0], 0.0, 0.0], [0.0, c[1], 0.0], [0.0, 0.0, c[2]]];
+        let xyz_from_rgb = rgb * SquareMatrix::diag(&[c.r, c.g, c.b]);
-        let xyz_from_rgb = rgb * SquareMatrix::new(xyz_from_rgb_m);
         let rgb_from_xyz = xyz_from_rgb
             .inverse()
-            .expect("Failed to invert the XYZfromRGB matrix. Is it singular?");
+            .expect("XYZ from RGB matrix is singular");
 
         Ok(Self {
             r,
@@ -66,18 +65,18 @@ impl RGBColorspace {
     }
 
     pub fn to_xyz(&self, rgb: RGB) -> XYZ {
-        self.xyz_from_rgb.transform_to_xyz(rgb)
+        self.xyz_from_rgb * rgb
     }
 
     pub fn to_rgb(&self, xyz: XYZ) -> RGB {
-        self.rgb_from_xyz.transform_to_rgb(xyz)
+        self.rgb_from_xyz * xyz
     }
 
     pub fn to_rgb_coeffs(&self, rgb: RGB) -> RGBSigmoidPolynomial {
         self.rgb_to_spectrum_table.to_polynomial(rgb)
     }
 
-    pub fn convert_colorspace(&self, other: &RGBColorspace) -> SquareMatrix<Float, 3> {
+    pub fn convert_colorspace(&self, other: &RGBColorSpace) -> SquareMatrix<Float, 3> {
         if self == other {
             return SquareMatrix::default();
         }
@@ -86,7 +85,7 @@ impl RGBColorspace {
     }
 
     pub fn srgb() -> &'static Self {
-        static SRGB_SPACE: Lazy<RGBColorspace> = Lazy::new(|| {
+        static SRGB_SPACE: Lazy<RGBColorSpace> = Lazy::new(|| {
             let r = Point2f::new(0.64, 0.33);
             let g = Point2f::new(0.30, 0.60);
             let b = Point2f::new(0.15, 0.06);
@@ -94,7 +93,7 @@ impl RGBColorspace {
             let illuminant = Spectrum::std_illuminant_d65();
             let table = RGBToSpectrumTable::srgb();
 
-            RGBColorspace::new(r, g, b, illuminant, table)
+            RGBColorSpace::new(r, g, b, illuminant, table)
                 .expect("Failed to initialize standard sRGB color space")
         });
 
@@ -102,7 +101,7 @@ impl RGBColorspace {
     }
 }
 
-impl PartialEq for RGBColorspace {
+impl PartialEq for RGBColorSpace {
     fn eq(&self, other: &Self) -> bool {
         self.r == other.r
             && self.g == other.g

src/utils/containers.rs

@@ -76,7 +76,7 @@ impl<T> Index<Point2i> for Array2D<T> {
     type Output = T;
 
     fn index(&self, mut p: Point2i) -> &Self::Output {
-        p = p - Vector2i::from(self.extent.p_min);
+        p -= Vector2i::from(self.extent.p_min);
         let width = self.extent.p_max.x() - self.extent.p_min.x();
         let idx = (p.x() + width * p.y()) as usize;
         &self.values[idx]
@@ -85,7 +85,7 @@ impl<T> Index<Point2i> for Array2D<T> {
 
 impl<T> IndexMut<Point2i> for Array2D<T> {
     fn index_mut(&mut self, mut p: Point2i) -> &mut Self::Output {
-        p = p - Vector2i::from(self.extent.p_min);
+        p -= Vector2i::from(self.extent.p_min);
         let width = self.extent.p_max.x() - self.extent.p_min.x();
         let idx = (p.x() + width * p.y()) as usize;
         &mut self.values[idx]

src/utils/error.rs

@@ -1,8 +1,8 @@
-use image::error;
+use image_rs::{ImageError as IError, error};
 use std::fmt;
 use thiserror::Error;
 
-use crate::utils::image::PixelFormat;
+use crate::image::PixelFormat;
 
 #[derive(Error, Debug)]
 pub enum LlsError {
@@ -40,7 +40,7 @@ pub enum ImageError {
     Io(#[from] std::io::Error),
 
     #[error("Image file error: {0}")]
-    Image(#[from] image::ImageError),
+    Image(#[from] IError),
 
     #[error("EXR file error: {0}")]
     Exr(#[from] exr::error::Error),

src/utils/interval.rs

@@ -46,20 +46,20 @@ impl Interval {
     }
 
     pub fn abs(&self) -> Self {
-    if self.low >= 0.0 {
+        if self.low >= 0.0 {
-        return *self;
+            return *self;
-    }
+        }
 
-    if self.high < 0.0 {
+        if self.high < 0.0 {
-        return -(*self);
+            return -(*self);
-    }
+        }
 
-    Self {
+        Self {
-        low: 0.0,
+            low: 0.0,
-        high: next_float_up((-self.low).max(self.high)),
+            high: next_float_up((-self.low).max(self.high)),
+        }
     }
 }
-}
 
 impl Default for Interval {
     fn default() -> Self {

src/utils/math.rs

@@ -12,6 +12,7 @@ use num_traits::{Float as NumFloat, Num, One, Signed, Zero};
 use rayon::prelude::*;
 use std::error::Error;
 use std::fmt::{self, Display};
+use std::iter::{Product, Sum};
 use std::mem;
 use std::ops::{Add, Div, Index, IndexMut, Mul, Neg, Rem};
 
@@ -49,7 +50,7 @@ where
     let cd = c * d;
     let difference_of_products = fma(a, b, -cd);
     let error = fma(-c, d, cd);
-    return difference_of_products + error;
+    difference_of_products + error
 }
 
 #[inline]
@@ -60,7 +61,7 @@ where
     let cd = c * d;
     let sum_of_products = fma(a, b, cd);
     let error = fma(c, d, -cd);
-    return sum_of_products + error;
+    sum_of_products + error
 }
 
 #[inline]
@@ -82,7 +83,7 @@ pub fn safe_asin<T: NumFloat>(x: T) -> T {
 
 #[inline]
 pub fn safe_acos(x: Float) -> Float {
-    if x >= -1.0001 && x <= 1.0001 {
+    if (-1.001..1.001).contains(&x) {
         clamp_t(x, -1., 1.).asin()
     } else {
         panic!("Not valid value for acos")
@@ -94,12 +95,12 @@ pub fn sinx_over_x(x: Float) -> Float {
     if 1. - x * x == 1. {
         return 1.;
     }
-    return x.sin() / x;
+    x.sin() / x
 }
 
 #[inline]
 pub fn sinc(x: Float) -> Float {
-    return sinx_over_x(PI * x);
+    sinx_over_x(PI * x)
 }
 
 #[inline]
@@ -107,7 +108,7 @@ pub fn windowed_sinc(x: Float, radius: Float, tau: Float) -> Float {
     if x.abs() > radius {
         return 0.;
     }
-    return sinc(x) * sinc(x / tau);
+    sinc(x) * sinc(x / tau)
 }
 
 #[inline]
@@ -284,11 +285,11 @@ pub fn equal_area_square_to_sphere(p: Point2f) -> Vector3f {
     // Compute $\cos\phi$ and $\sin\phi$ for original quadrant and return vector
     let cos_phi = phi.cos().copysign(u);
     let sin_phi = phi.sin().copysign(v);
-    return Vector3f::new(
+    Vector3f::new(
         cos_phi * r * (2. - square(r)).sqrt(),
         sin_phi * r * (2. - square(r)).sqrt(),
         z,
-    );
+    )
 }
 
 pub fn gaussian(x: Float, y: Float, sigma: Float) -> Float {
@@ -420,9 +421,9 @@ pub fn sample_tent(u: Float, r: Float) -> Float {
     let offset = sample_discrete(&[0.5, 0.5], u, None, Some(&mut u_remapped))
         .expect("Discrete sampling shouldn't fail");
     if offset == 0 {
-        return -r + r * sample_linear(u, 0., 1.);
+        -r + r * sample_linear(u, 0., 1.)
     } else {
-        return r * sample_linear(u, 1., 0.);
+        r * sample_linear(u, 1., 0.)
     }
 }
 
@@ -589,14 +590,14 @@ impl DigitPermutation {
             inv_base_m *= inv_base;
         }
 
-        let mut permutations = vec![0u16; (n_digits * base) as usize];
+        let mut permutations = vec![0u16; n_digits * base];
 
         for digit_index in 0..n_digits {
-            let hash_input = [base as u64, digit_index as u64, seed as u64];
+            let hash_input = [base as u64, digit_index as u64, seed];
             let dseed = hash_buffer(&hash_input, 0);
 
             for digit_value in 0..base {
-                let index = (digit_index * base + digit_value) as usize;
+                let index = digit_index * base + digit_value;
 
                 permutations[index] =
                     permutation_element(digit_value as u32, base as u32, dseed as u32) as u16;
@@ -870,7 +871,7 @@ pub fn sobol_sample<S: Scrambler>(mut a: u64, dimension: usize, randomizer: S) -
 
     while a != 0 {
         if (a & 1) != 0 {
-            v ^= SOBOL_MATRICES_32[i] as u32;
+            v ^= SOBOL_MATRICES_32[i];
         }
         a >>= 1;
         i += 1;
@@ -1059,6 +1060,7 @@ impl<T: Display, const R: usize, const C: usize> Display for Matrix<T, R, C> {
 
         for i in 0..R {
             write!(f, "[")?;
+            #[allow(clippy::needless_range_loop)]
             for j in 0..C {
                 write!(f, "{: >width$} ", self.m[i][j], width = col_widths[j])?;
             }
@@ -1091,11 +1093,11 @@ impl<T, const N: usize> SquareMatrix<T, N> {
     where
         T: Copy + Zero + One,
     {
-        let mut m = [[T::zero(); N]; N];
+        Self {
-        for i in 0..N {
+            m: std::array::from_fn(|i| {
-            m[i][i] = T::one();
+                std::array::from_fn(|j| if i == j { T::one() } else { T::zero() })
+            }),
         }
-        Self { m }
     }
 
     pub fn diag(v: &[T]) -> Self
@@ -1139,7 +1141,10 @@ impl<T, const N: usize> SquareMatrix<T, N> {
     }
 }
 
-impl<T: NumFloat + Copy + PartialEq, const N: usize> SquareMatrix<T, N> {
+impl<T, const N: usize> SquareMatrix<T, N>
+where
+    T: NumFloat + Sum + Product + Copy,
+{
     pub fn inverse(&self) -> Result<Self, InversionError> {
         if N == 0 {
             return Err(InversionError::EmptyMatrix);
@@ -1149,12 +1154,9 @@ impl<T: NumFloat + Copy + PartialEq, const N: usize> SquareMatrix<T, N> {
         let mut inv = Self::identity();
 
         for i in 0..N {
-            let mut pivot_row = i;
+            let pivot_row = (i..N)
-            for j in (i + 1)..N {
+                .max_by(|&a, &b| mat[a][i].abs().partial_cmp(&mat[b][i].abs()).unwrap())
-                if mat[j][i].abs() > mat[pivot_row][i].abs() {
+                .unwrap_or(i);
-                    pivot_row = j;
-                }
-            }
 
             if pivot_row != i {
                 mat.swap(i, pivot_row);
@@ -1166,20 +1168,19 @@ impl<T: NumFloat + Copy + PartialEq, const N: usize> SquareMatrix<T, N> {
                 return Err(InversionError::SingularMatrix);
             }
 
-            for j in i..N {
+            let inv_pivot = T::one() / pivot;
-                mat[i][j] = mat[i][j] / pivot;
+            mat[i][i..].iter_mut().for_each(|x| *x = *x * inv_pivot);
-            }
+            inv[i].iter_mut().for_each(|x| *x = *x * inv_pivot);
-
-            for j in 0..N {
-                inv.m[i][j] = inv.m[i][j] / pivot;
-            }
 
             for j in 0..N {
                 if i != j {
                     let factor = mat[j][i];
+                    #[allow(clippy::needless_range_loop)]
                     for k in i..N {
                         mat[j][k] = mat[j][k] - factor * mat[i][k];
                     }
+
+                    #[allow(clippy::needless_range_loop)]
                     for k in 0..N {
                         inv.m[j][k] = inv.m[j][k] - factor * inv.m[i][k];
                     }
@@ -1237,12 +1238,9 @@ impl<T: NumFloat + Copy + PartialEq, const N: usize> SquareMatrix<T, N> {
                 let mut parity = 0;
 
                 for i in 0..N {
-                    let mut max_row = i;
+                    let max_row = (i..N)
-                    for k in (i + 1)..N {
+                        .max_by(|&a, &b| lum[a][i].abs().partial_cmp(&lum[b][i].abs()).unwrap())
-                        if lum[k][i].abs() > lum[max_row][i].abs() {
+                        .unwrap_or(i);
-                            max_row = k;
-                        }
-                    }
 
                     if max_row != i {
                         lum.swap(i, max_row);
@@ -1257,22 +1255,17 @@ impl<T: NumFloat + Copy + PartialEq, const N: usize> SquareMatrix<T, N> {
                     // Gaussian elimination
                     for j in (i + 1)..N {
                         let factor = lum[j][i] / lum[i][i];
+
+                        #[allow(clippy::needless_range_loop)]
                         for k in i..N {
-                            lum[j][k] = lum[j][k] - factor * lum[i][k];
+                            let val_i = lum[i][k];
+                            lum[j][k] = lum[j][k] - factor * val_i;
                         }
                     }
                 }
 
-                let mut det = T::one();
+                let det: T = (0..N).map(|i| lum[i][i]).product();
-                for i in 0..N {
+                if parity < 0 { -det } else { det }
-                    det = det * lum[i][i];
-                }
-
-                if parity % 2 != 0 {
-                    det = -det;
-                }
-
-                det
             }
         }
     }
@@ -1292,81 +1285,25 @@ impl<T, const N: usize> IndexMut<usize> for SquareMatrix<T, N> {
     }
 }
 
-impl<T: Num + Copy, const N: usize> Mul<Vector<T, N>> for SquareMatrix<T, N> {
+impl<T, const N: usize> Mul<Vector<T, N>> for SquareMatrix<T, N>
+where
+    T: Copy + Mul<Output = T> + Sum + Default,
+{
     type Output = Vector<T, N>;
     fn mul(self, rhs: Vector<T, N>) -> Self::Output {
-        let mut result_arr = [T::zero(); N];
+        let arr = std::array::from_fn(|i| self.m[i].iter().zip(&rhs.0).map(|(m, v)| *m * *v).sum());
-        for i in 0..N {
+        Vector(arr)
-            for j in 0..N {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs.0[j];
-            }
-        }
-        Vector(result_arr)
     }
 }
 
-impl<T: Num + Copy, const N: usize> Mul<Point<T, N>> for SquareMatrix<T, N> {
+impl<T, const N: usize> Mul<Point<T, N>> for SquareMatrix<T, N>
+where
+    T: Copy + Mul<Output = T> + Sum + Default,
+{
     type Output = Point<T, N>;
     fn mul(self, rhs: Point<T, N>) -> Self::Output {
-        let mut result_arr = [T::zero(); N];
+        let arr = std::array::from_fn(|i| self.m[i].iter().zip(&rhs.0).map(|(m, v)| *m * *v).sum());
-        for i in 0..N {
+        Point(arr)
-            for j in 0..N {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs.0[j];
-            }
-        }
-        Point(result_arr)
-    }
-}
-
-impl SquareMatrix<Float, 3> {
-    pub fn transform_to_xyz(&self, rgb: RGB) -> XYZ {
-        let r = rgb[0];
-        let g = rgb[1];
-        let b = rgb[2];
-
-        let x = self.m[0][0] * r + self.m[0][1] * g + self.m[0][2] * b;
-        let y = self.m[1][0] * r + self.m[1][1] * g + self.m[1][2] * b;
-        let z = self.m[2][0] * r + self.m[2][1] * g + self.m[2][2] * b;
-
-        XYZ::new(x, y, z)
-    }
-
-    pub fn transform_to_rgb(&self, xyz: XYZ) -> RGB {
-        let x = xyz[0];
-        let y = xyz[1];
-        let z = xyz[2];
-
-        let r = self.m[0][0] * x + self.m[0][1] * y + self.m[0][2] * z;
-        let g = self.m[1][0] * x + self.m[1][1] * y + self.m[1][2] * z;
-        let b = self.m[2][0] * x + self.m[2][1] * y + self.m[2][2] * z;
-
-        RGB::new(r, g, b)
-    }
-}
-
-impl Mul<XYZ> for SquareMatrix<Float, 3> {
-    type Output = XYZ;
-    fn mul(self, rhs: XYZ) -> Self::Output {
-        let mut result_arr = [0.0; 3];
-        for i in 0..3 {
-            for j in 0..3 {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs[j];
-            }
-        }
-        XYZ::new(result_arr[0], result_arr[1], result_arr[2])
-    }
-}
-
-impl Mul<RGB> for SquareMatrix<Float, 3> {
-    type Output = RGB;
-    fn mul(self, rhs: RGB) -> Self::Output {
-        let mut result_arr = [0.0; 3];
-        for i in 0..3 {
-            for j in 0..3 {
-                result_arr[i] = result_arr[i] + self.m[i][j] * rhs[j];
-            }
-        }
-        RGB::new(result_arr[0], result_arr[1], result_arr[2])
     }
 }
 

src/utils/mesh.rs

@@ -19,6 +19,7 @@ pub struct TriangleMesh {
 }
 
 impl TriangleMesh {
+    #[allow(clippy::too_many_arguments)]
     pub fn new(
         render_from_object: &Transform<Float>,
         reverse_orientation: bool,

src/utils/mipmap.rs

@@ -0,0 +1,260 @@
+use crate::core::pbrt::{Float, lerp};
+use crate::geometry::{Lerp, Point2f, Point2i, Vector2f, VectorLike};
+use crate::image::{Image, PixelData, PixelFormat, WrapMode, WrapMode2D};
+use crate::utils::color::RGB;
+use crate::utils::colorspace::RGBColorSpace;
+
+use std::hash::{Hash, Hasher};
+use std::ops::{Add, Mul, Sub};
+use std::sync::Arc;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
+pub enum FilterFunction {
+    Point,
+    Bilinear,
+    Trilinear,
+    Ewa,
+}
+
+impl std::fmt::Display for FilterFunction {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        let s = match self {
+            FilterFunction::Ewa => "EWA",
+            FilterFunction::Trilinear => "trilinear",
+            FilterFunction::Bilinear => "bilinear",
+            FilterFunction::Point => "point",
+        };
+        write!(f, "{}", s)
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct MIPMapFilterOptions {
+    pub filter: FilterFunction,
+    pub max_anisotropy: f32,
+}
+
+impl Default for MIPMapFilterOptions {
+    fn default() -> Self {
+        Self {
+            filter: FilterFunction::Ewa,
+            max_anisotropy: 8.0,
+        }
+    }
+}
+
+impl PartialEq for MIPMapFilterOptions {
+    fn eq(&self, other: &Self) -> bool {
+        self.filter == other.filter
+            && self.max_anisotropy.to_bits() == other.max_anisotropy.to_bits()
+    }
+}
+
+impl Eq for MIPMapFilterOptions {}
+
+impl Hash for MIPMapFilterOptions {
+    fn hash<H: Hasher>(&self, state: &mut H) {
+        self.filter.hash(state);
+        // Hash the bits, not the float value
+        self.max_anisotropy.to_bits().hash(state);
+    }
+}
+
+pub trait MIPMapSample:
+    Copy + Add<Output = Self> + Sub<Output = Self> + Mul<Float, Output = Self> + std::fmt::Debug
+{
+    fn zero() -> Self;
+    fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self;
+    fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self;
+}
+
+impl MIPMapSample for Float {
+    fn zero() -> Self {
+        0.
+    }
+
+    fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self {
+        image.bilerp_channel(st, 0, wrap)
+    }
+
+    fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self {
+        image.get_channel(st, 0, wrap)
+    }
+}
+
+impl MIPMapSample for RGB {
+    fn zero() -> Self {
+        RGB::new(0., 0., 0.)
+    }
+
+    fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self {
+        let nc = image.n_channels();
+        if nc >= 3 {
+            let r = image.bilerp_channel(st, 0, wrap);
+            let g = image.bilerp_channel(st, 1, wrap);
+            let b = image.bilerp_channel(st, 2, wrap);
+            RGB::new(r, g, b)
+        } else {
+            let v = image.bilerp_channel(st, 0, wrap);
+            RGB::new(v, v, v)
+        }
+    }
+
+    fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self {
+        let nc = image.n_channels();
+        if nc >= 3 {
+            let r = image.get_channel(st, 0, wrap);
+            let g = image.get_channel(st, 1, wrap);
+            let b = image.get_channel(st, 2, wrap);
+            RGB::new(r, g, b)
+        } else {
+            let v = image.get_channel(st, 0, wrap);
+            RGB::new(v, v, v)
+        }
+    }
+}
+
+#[derive(Debug)]
+pub struct MIPMap {
+    pyramid: Vec<Image>,
+    color_space: Arc<RGBColorSpace>,
+    wrap_mode: WrapMode,
+    options: MIPMapFilterOptions,
+}
+
+impl MIPMap {
+    pub fn new(
+        image: Image,
+        color_space: Arc<RGBColorSpace>,
+        wrap_mode: WrapMode,
+        options: MIPMapFilterOptions,
+    ) -> Self {
+        let pyramid = Image::generate_pyramid(image, wrap_mode);
+        Self {
+            pyramid,
+            color_space,
+            wrap_mode,
+            options,
+        }
+    }
+
+    pub fn level_resolution(&self, level: usize) -> Point2i {
+        self.pyramid[level].resolution()
+    }
+
+    pub fn levels(&self) -> usize {
+        self.pyramid.len()
+    }
+
+    pub fn get_rgb_colorspace(&self) -> Arc<RGBColorSpace> {
+        self.color_space.clone()
+    }
+
+    pub fn get_level(&self, level: usize) -> &Image {
+        &self.pyramid[level]
+    }
+
+    pub fn filter<T: MIPMapSample>(
+        &self,
+        st: Point2f,
+        mut dst0: Vector2f,
+        mut dst1: Vector2f,
+    ) -> T {
+        if self.options.filter != FilterFunction::Ewa {
+            // Compute largest change in texture coordinates
+            let width = 2.0
+                * [
+                    dst0.x().abs(),
+                    dst0.y().abs(),
+                    dst1.x().abs(),
+                    dst1.y().abs(),
+                ]
+                .into_iter()
+                .reduce(Float::max)
+                .unwrap_or(0.0);
+
+            // Compute MIP Map level
+            // n_levels - 1 + log2(width) maps width=1.0 to the top level (1x1)
+            let n_levels = self.levels() as Float;
+            let level = n_levels - 1.0 + width.max(1e-8).log2();
+
+            // Handle very wide filter (return the 1x1 average pixel)
+            if level >= n_levels - 1.0 {
+                return self.texel(self.levels() - 1, Point2i::new(0, 0));
+            }
+
+            let i_level = level.floor() as usize;
+
+            return match self.options.filter {
+                FilterFunction::Point => {
+                    // Nearest Neighbor
+                    let resolution = self.level_resolution(i_level);
+                    let sti = Point2i::new(
+                        (st.x() * resolution.x() as Float - 0.5).round() as i32,
+                        (st.y() * resolution.y() as Float - 0.5).round() as i32,
+                    );
+                    self.texel(i_level, sti)
+                }
+                FilterFunction::Bilinear => self.bilerp(i_level, st),
+                FilterFunction::Trilinear => {
+                    // Interpolate between current level and next level
+                    let v0 = self.bilerp(i_level, st);
+                    let v1 = self.bilerp(i_level + 1, st);
+                    let t = level - i_level as Float;
+                    lerp(t, v0, v1)
+                }
+                FilterFunction::Ewa => unreachable!(),
+            };
+        }
+
+        if dst0.norm_squared() < dst1.norm_squared() {
+            std::mem::swap(&mut dst0, &mut dst1);
+        }
+
+        let longer_len = dst0.norm();
+        let mut shorter_len = dst1.norm();
+
+        // If the ellipse is too thin, we fatten the minor axis to limit the number
+        // of texels we have to filter, trading accuracy for performance/cache hits.
+        if shorter_len * self.options.max_anisotropy < longer_len && shorter_len > 0.0 {
+            let scale = longer_len / (shorter_len * self.options.max_anisotropy);
+            dst1 *= scale;
+            shorter_len *= scale;
+        }
+
+        if shorter_len == 0.0 {
+            return self.bilerp(0, st);
+        }
+
+        // EWA also interpolates between two MIP levels, but does complex filtering on both.
+        let lod = (self.levels() as Float - 1.0 + shorter_len.log2()).max(0.0);
+        let ilod = lod.floor() as usize;
+
+        let v0 = self.ewa(ilod, st, dst0, dst1);
+        let v1 = self.ewa(ilod + 1, st, dst0, dst1);
+
+        lerp(lod - ilod as Float, v0, v1)
+    }
+
+    fn texel<T: MIPMapSample>(&self, _level: usize, _st: Point2i) -> T {
+        todo!()
+    }
+
+    fn bilerp<T: MIPMapSample>(&self, level: usize, st: Point2f) -> T {
+        let image = &self.pyramid[level];
+        let wrap_2d = WrapMode2D {
+            uv: [self.wrap_mode; 2],
+        };
+        T::sample_bilerp(image, st, wrap_2d)
+    }
+
+    fn ewa<T: MIPMapSample>(
+        &self,
+        _scale: usize,
+        _st: Point2f,
+        _dst0: Vector2f,
+        _dst1: Vector2f,
+    ) -> T {
+        todo!()
+    }
+}

src/utils/mod.rs

@@ -3,16 +3,15 @@ pub mod colorspace;
 pub mod containers;
 pub mod error;
 pub mod hash;
-pub mod image;
 pub mod interval;
 pub mod math;
 pub mod mesh;
+pub mod mipmap;
 pub mod quaternion;
 pub mod rng;
 pub mod sampling;
 pub mod scattering;
 pub mod sobol;
-pub mod spectrum;
 pub mod splines;
 pub mod transform;
 

src/utils/quaternion.rs

@@ -134,17 +134,17 @@ impl Quaternion {
     #[inline]
     pub fn angle_between(&self, rhs: Quaternion) -> Float {
         if self.dot(rhs) < 0.0 {
-            return PI - 2. * safe_asin((self.v + rhs.v).norm() / 2.);
+            PI - 2. * safe_asin((self.v + rhs.v).norm() / 2.)
         } else {
-            return 2. * safe_asin((rhs.v - self.v).norm() / 2.);
+            2. * safe_asin((rhs.v - self.v).norm() / 2.)
         }
     }
 
     pub fn slerp(t: Float, q1: Quaternion, q2: Quaternion) -> Quaternion {
         let theta = q1.angle_between(q2);
         let sin_theta_over_theta = sinx_over_x(theta);
-        return q1 * (1. - t) * sinx_over_x((1. - t) * theta) / sin_theta_over_theta
+        q1 * (1. - t) * sinx_over_x((1. - t) * theta) / sin_theta_over_theta
-            + q2 * t * sinx_over_x(t * theta) / sin_theta_over_theta;
+            + q2 * t * sinx_over_x(t * theta) / sin_theta_over_theta
     }
 
     pub fn length(&self) -> Float {

src/utils/sampling.rs

@@ -289,12 +289,11 @@ pub fn invert_spherical_rectangle_sample(
     let hv = [lerp(u1[0], h0, h1), lerp(u1[1], h0, h1)];
     let hvsq = [square(hv[0]), square(hv[1])];
     let yz = [(hv[0] * dd) / (1. - hvsq[0]), (hv[1] * dd) / (1. - hvsq[1])];
-    let u = if (yz[0] - yv).abs() < (yz[1] - yv).abs() {
+    if (yz[0] - yv).abs() < (yz[1] - yv).abs() {
         Point2f::new(clamp_t(u0, 0., 1.), u1[0])
     } else {
         Point2f::new(clamp_t(u0, 0., 1.), u1[1])
-    };
+    }
-    u
 }
 
 pub fn sample_spherical_triangle(
@@ -552,13 +551,13 @@ impl PiecewiseConstant1D {
 
         let func_integral = cdf[n];
         if func_integral == 0. {
-            for i in 1..=n {
+            let n_float = n as Float;
-                cdf[i] = i as Float / n as Float;
+            cdf.iter_mut()
-            }
+                .enumerate()
+                .for_each(|(i, c)| *c = i as Float / n_float);
         } else {
-            for i in 1..=n {
+            let inv_integral = 1.0 / func_integral;
-                cdf[i] /= func_integral;
+            cdf.iter_mut().for_each(|c| *c *= inv_integral);
-            }
         }
 
         Self {
@@ -629,7 +628,7 @@ impl PiecewiseConstant2D {
     }
 
     pub fn new_with_bounds(data: &Array2D<Float>, domain: Bounds2f) -> Self {
-        Self::new(data, data.x_size() as usize, data.y_size() as usize, domain)
+        Self::new(data, data.x_size(), data.y_size(), domain)
     }
 
     pub fn new_with_data(data: &Array2D<Float>, nx: usize, ny: usize) -> Self {
@@ -1047,14 +1046,16 @@ impl<const N: usize> PiecewiseLinear2D<N> {
         for mask in 0..num_corners {
             let mut offset = 0u32;
             let mut weight: Float = 1.0;
-            for d in 0..N {
+            let mut current_mask = mask;
-                let bit = (mask >> d) & 1;
+            for (weights, &stride) in param_weight.iter().zip(&self.param_strides) {
-                if bit == 1 {
+                if (current_mask & 1) == 1 {
-                    offset += self.param_strides[d] * size;
+                    offset += stride * size;
-                    weight *= param_weight[d].1;
+                    weight *= weights.1;
                 } else {
-                    weight *= param_weight[d].0;
+                    weight *= weights.0;
                 }
+
+                current_mask >>= 1;
             }
             result += weight * data[(i0 + offset) as usize];
         }

src/utils/scattering.rs

@@ -1,11 +1,11 @@
 use super::math::safe_sqrt;
 use super::sampling::sample_uniform_disk_polar;
-use super::spectrum::{N_SPECTRUM_SAMPLES, SampledSpectrum};
 use crate::core::pbrt::{Float, PI, clamp_t, lerp};
 use crate::geometry::{
     Normal3f, Point2f, Vector2f, Vector3f, VectorLike, abs_cos_theta, cos_phi, cos2_theta, sin_phi,
     tan2_theta,
 };
+use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum};
 use crate::utils::math::square;
 
 use num::complex::Complex;
@@ -151,7 +151,7 @@ pub fn fr_dielectric(cos_theta_i: Float, eta: Float) -> Float {
 pub fn fr_complex(cos_theta_i: Float, eta: Complex<Float>) -> Float {
     let cos_corr = clamp_t(cos_theta_i, 0., 1.);
     let sin2_theta_i = 1. - square(cos_corr);
-    let sin2_theta_t: Complex<Float> = (sin2_theta_i / square(eta)).into();
+    let sin2_theta_t: Complex<Float> = sin2_theta_i / square(eta);
     let cos2_theta_t: Complex<Float> = (1. - sin2_theta_t).sqrt();
 
     let r_parl = (eta * cos_corr - cos2_theta_t) / (eta * cos_corr + cos2_theta_t);

src/utils/spectrum.rs

@@ -1,789 +0,0 @@
-use once_cell::sync::Lazy;
-use std::fmt;
-use std::ops::{
-    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
-};
-use std::sync::Arc;
-
-use super::color::{RGB, RGBSigmoidPolynomial, XYZ};
-use super::colorspace::RGBColorspace;
-use crate::core::cie;
-use crate::core::pbrt::Float;
-
-pub const CIE_Y_INTEGRAL: Float = 106.856895;
-pub const N_SPECTRUM_SAMPLES: usize = 1200;
-pub const LAMBDA_MIN: i32 = 360;
-pub const LAMBDA_MAX: i32 = 830;
-
-#[derive(Debug, Copy, Clone)]
-pub struct SampledSpectrum {
-    values: [Float; N_SPECTRUM_SAMPLES],
-}
-
-impl Default for SampledSpectrum {
-    fn default() -> Self {
-        Self {
-            values: [0.0; N_SPECTRUM_SAMPLES],
-        }
-    }
-}
-
-impl SampledSpectrum {
-    pub fn new(c: Float) -> Self {
-        Self {
-            values: [c; N_SPECTRUM_SAMPLES],
-        }
-    }
-
-    pub fn from_vector(v: Vec<Float>) -> Self {
-        let mut s = Self::new(0.0);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            s.values[i] = v[i];
-        }
-        s
-    }
-
-    pub fn is_black(&self) -> bool {
-        self.values.iter().all(|&sample| sample == 0.0)
-    }
-
-    pub fn has_nans(&self) -> bool {
-        self.values.iter().any(|&v| v.is_nan())
-    }
-
-    pub fn min_component_value(&self) -> Float {
-        self.values.iter().fold(Float::INFINITY, |a, &b| a.min(b))
-    }
-
-    pub fn max_component_value(&self) -> Float {
-        self.values
-            .iter()
-            .fold(Float::NEG_INFINITY, |a, &b| a.max(b))
-    }
-
-    pub fn average(&self) -> Float {
-        self.values.iter().sum::<Float>() / (N_SPECTRUM_SAMPLES as Float)
-    }
-
-    pub fn safe_div(&self, rhs: SampledSpectrum) -> Self {
-        let mut r = SampledSpectrum::new(0.0);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            r.values[i] = if rhs[i] != 0.0 {
-                self.values[i] / rhs.values[i]
-            } else {
-                0.0
-            }
-        }
-        r
-    }
-
-    pub fn to_xyz(&self, lambda: &SampledWavelengths) -> XYZ {
-        let x = spectra::X.sample(lambda);
-        let y = spectra::Y.sample(lambda);
-        let z = spectra::Z.sample(lambda);
-        let pdf = lambda.pdf();
-        XYZ::new(
-            (*self * x).safe_div(pdf).average(),
-            (*self * y).safe_div(pdf).average(),
-            (*self * z).safe_div(pdf).average(),
-        ) / CIE_Y_INTEGRAL
-    }
-
-    pub fn to_rgb(&self, lambda: &SampledWavelengths, c: &RGBColorspace) -> RGB {
-        let xyz = self.to_xyz(lambda);
-        c.to_rgb(xyz)
-    }
-
-    pub fn exp(&self) -> SampledSpectrum {
-        let values = self.values.map(|v| v.exp());
-        let ret = Self { values };
-        debug_assert!(!ret.has_nans());
-        ret
-    }
-
-    pub fn pow_int(&self, mut power: usize) -> Self {
-        let mut result = Self::new(1.0);
-        let mut base = *self;
-        while power > 0 {
-            if power % 2 == 1 {
-                result = result * base;
-            }
-            base = base * base;
-            power /= 2;
-        }
-        result
-    }
-}
-
-impl Index<usize> for SampledSpectrum {
-    type Output = Float;
-    fn index(&self, i: usize) -> &Self::Output {
-        &self.values[i]
-    }
-}
-
-impl IndexMut<usize> for SampledSpectrum {
-    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
-        &mut self.values[i]
-    }
-}
-
-impl Add for SampledSpectrum {
-    type Output = Self;
-    fn add(self, rhs: Self) -> Self {
-        let mut ret = self;
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] += rhs.values[i];
-        }
-        ret
-    }
-}
-
-impl AddAssign for SampledSpectrum {
-    fn add_assign(&mut self, rhs: Self) {
-        for i in 0..N_SPECTRUM_SAMPLES {
-            self.values[i] += rhs.values[i];
-        }
-    }
-}
-
-impl Sub for SampledSpectrum {
-    type Output = Self;
-    fn sub(self, rhs: Self) -> Self {
-        let mut ret = self;
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] -= rhs.values[i];
-        }
-        ret
-    }
-}
-
-impl SubAssign for SampledSpectrum {
-    fn sub_assign(&mut self, rhs: Self) {
-        for i in 0..N_SPECTRUM_SAMPLES {
-            self.values[i] -= rhs.values[i];
-        }
-    }
-}
-
-impl Sub<SampledSpectrum> for Float {
-    type Output = SampledSpectrum;
-    fn sub(self, rhs: SampledSpectrum) -> SampledSpectrum {
-        let mut ret = SampledSpectrum::new(0.0);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] = self - rhs.values[i];
-        }
-        ret
-    }
-}
-
-impl Mul for SampledSpectrum {
-    type Output = Self;
-    fn mul(self, rhs: Self) -> Self {
-        let mut ret = self;
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] *= rhs.values[i];
-        }
-        ret
-    }
-}
-
-impl MulAssign for SampledSpectrum {
-    fn mul_assign(&mut self, rhs: Self) {
-        for i in 0..N_SPECTRUM_SAMPLES {
-            self.values[i] *= rhs.values[i];
-        }
-    }
-}
-
-impl Mul<Float> for SampledSpectrum {
-    type Output = Self;
-    fn mul(self, rhs: Float) -> Self {
-        let mut ret = self;
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] *= rhs;
-        }
-        ret
-    }
-}
-
-impl Mul<SampledSpectrum> for Float {
-    type Output = SampledSpectrum;
-    fn mul(self, rhs: SampledSpectrum) -> SampledSpectrum {
-        rhs * self
-    }
-}
-
-impl MulAssign<Float> for SampledSpectrum {
-    fn mul_assign(&mut self, rhs: Float) {
-        for i in 0..N_SPECTRUM_SAMPLES {
-            self.values[i] *= rhs;
-        }
-    }
-}
-
-impl DivAssign for SampledSpectrum {
-    fn div_assign(&mut self, rhs: Self) {
-        for i in 0..N_SPECTRUM_SAMPLES {
-            debug_assert_ne!(0.0, rhs.values[i]);
-            self.values[i] /= rhs.values[i];
-        }
-    }
-}
-
-impl Div for SampledSpectrum {
-    type Output = Self;
-    fn div(self, rhs: Self) -> Self::Output {
-        let mut ret = self;
-        ret /= rhs;
-        ret
-    }
-}
-impl Div<Float> for SampledSpectrum {
-    type Output = Self;
-
-    fn div(self, rhs: Float) -> Self::Output {
-        debug_assert_ne!(rhs, 0.0);
-        let mut ret = self;
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] /= rhs;
-        }
-        ret
-    }
-}
-
-impl DivAssign<Float> for SampledSpectrum {
-    fn div_assign(&mut self, rhs: Float) {
-        debug_assert_ne!(rhs, 0.0);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            self.values[i] /= rhs;
-        }
-    }
-}
-
-impl Neg for SampledSpectrum {
-    type Output = Self;
-
-    fn neg(self) -> Self::Output {
-        let mut ret = SampledSpectrum::new(0.0);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            ret.values[i] = -self.values[i];
-        }
-        ret
-    }
-}
-
-#[derive(Debug, Copy, Clone, PartialEq)]
-pub struct SampledWavelengths {
-    pub lambda: [Float; N_SPECTRUM_SAMPLES],
-    pub pdf: [Float; N_SPECTRUM_SAMPLES],
-}
-
-impl SampledWavelengths {
-    pub fn pdf(&self) -> SampledSpectrum {
-        SampledSpectrum::from_vector(self.pdf.to_vec())
-    }
-
-    pub fn secondary_terminated(&self) -> bool {
-        for i in 1..N_SPECTRUM_SAMPLES {
-            if self.pdf[i] != 0.0 {
-                return false;
-            }
-        }
-        true
-    }
-
-    pub fn terminate_secondary(&mut self) {
-        if !self.secondary_terminated() {
-            for i in 1..N_SPECTRUM_SAMPLES {
-                self.pdf[i] = 0.0;
-            }
-            self.pdf[0] /= N_SPECTRUM_SAMPLES as Float;
-        }
-    }
-
-    pub fn sample_uniform(u: Float, lambda_min: Float, lambda_max: Float) -> Self {
-        let mut lambda = [0.0; N_SPECTRUM_SAMPLES];
-        lambda[0] = crate::core::pbrt::lerp(u, lambda_min, lambda_min);
-        let delta = (lambda_max - lambda_min) / N_SPECTRUM_SAMPLES as Float;
-        for i in 1..N_SPECTRUM_SAMPLES {
-            lambda[i] = lambda[i - 1] + delta;
-            if lambda[i] > lambda_max {
-                lambda[i] = lambda_min + (lambda[i] - lambda_max);
-            }
-        }
-
-        let mut pdf = [0.0; N_SPECTRUM_SAMPLES];
-        for i in 0..N_SPECTRUM_SAMPLES {
-            pdf[i] = 1.0 / (lambda_max - lambda_min);
-        }
-
-        Self { lambda, pdf }
-    }
-
-    pub fn sample_visible_wavelengths(u: Float) -> Float {
-        538.0 - 138.888889 * Float::atanh(0.85691062 - 1.82750197 * u)
-    }
-
-    pub fn visible_wavelengths_pdf(lambda: Float) -> Float {
-        if lambda < 360.0 || lambda > 830.0 {
-            return 0.0;
-        }
-        0.0039398042 / (Float::cosh(0.0072 * (lambda - 538.0))).sqrt()
-    }
-
-    pub fn sample_visible(u: Float) -> Self {
-        let mut lambda = [0.0; N_SPECTRUM_SAMPLES];
-        let mut pdf = [0.0; N_SPECTRUM_SAMPLES];
-
-        for i in 0..N_SPECTRUM_SAMPLES {
-            let mut up = u + i as Float / N_SPECTRUM_SAMPLES as Float;
-            if up > 1.0 {
-                up -= 1.0;
-            }
-            lambda[i] = Self::sample_visible_wavelengths(up);
-            pdf[i] = Self::visible_wavelengths_pdf(lambda[i]);
-        }
-        Self { lambda, pdf }
-    }
-}
-
-impl Index<usize> for SampledWavelengths {
-    type Output = Float;
-    fn index(&self, i: usize) -> &Self::Output {
-        &self.lambda[i]
-    }
-}
-
-impl IndexMut<usize> for SampledWavelengths {
-    fn index_mut(&mut self, i: usize) -> &mut Self::Output {
-        &mut self.lambda[i]
-    }
-}
-
-#[derive(Debug, Clone)]
-pub enum Spectrum {
-    Constant(ConstantSpectrum),
-    DenselySampled(DenselySampledSpectrum),
-    PiecewiseLinear(PiecewiseLinearSpectrum),
-    Blackbody(BlackbodySpectrum),
-    RGBAlbedo(RGBAlbedoSpectrum),
-}
-
-impl Spectrum {
-    pub fn std_illuminant_d65() -> Self {
-        todo!()
-    }
-    pub fn sample_at(&self, lambda: Float) -> Float {
-        match self {
-            Spectrum::Constant(s) => s.sample_at(lambda),
-            Spectrum::DenselySampled(s) => s.sample_at(lambda),
-            Spectrum::PiecewiseLinear(s) => s.sample_at(lambda),
-            Spectrum::Blackbody(s) => s.sample_at(lambda),
-            Spectrum::RGBAlbedo(s) => s.sample_at(lambda),
-        }
-    }
-
-    pub fn max_value(&self) -> Float {
-        match self {
-            Spectrum::Constant(_s) => 0.0,
-            Spectrum::DenselySampled(_s) => 0.0,
-            Spectrum::PiecewiseLinear(_s) => 0.0,
-            Spectrum::Blackbody(_s) => 0.0,
-            Spectrum::RGBAlbedo(s) => s.max_value(),
-        }
-    }
-
-    pub fn sample(&self, wavelengths: &SampledWavelengths) -> SampledSpectrum {
-        let mut s = SampledSpectrum::default();
-        for i in 0..N_SPECTRUM_SAMPLES {
-            s[i] = self.sample_at(wavelengths[i]);
-        }
-        s
-    }
-
-    pub fn to_xyz(&self) -> XYZ {
-        XYZ::new(
-            inner_product(&spectra::X, self),
-            inner_product(&spectra::Y, self),
-            inner_product(&spectra::Z, self),
-        ) / CIE_Y_INTEGRAL
-    }
-}
-
-pub fn inner_product(f: &Spectrum, g: &Spectrum) -> Float {
-    let mut integral = 0.0;
-    for lambda in LAMBDA_MIN..=LAMBDA_MAX {
-        integral += f.sample_at(lambda as f32) * g.sample_at(lambda as f32);
-    }
-    integral
-}
-
-#[derive(Debug, Clone, Copy)]
-pub struct ConstantSpectrum {
-    c: Float,
-}
-
-impl ConstantSpectrum {
-    pub fn new(c: Float) -> Self {
-        Self { c }
-    }
-
-    pub fn sample_at(&self, _lambda: Float) -> Float {
-        self.c
-    }
-
-    fn max_value(&self) -> Float {
-        self.c
-    }
-}
-
-#[derive(Debug, Clone, Copy)]
-pub struct RGBAlbedoSpectrum {
-    rsp: RGBSigmoidPolynomial,
-}
-
-impl RGBAlbedoSpectrum {
-    pub fn new(cs: &RGBColorspace, rgb: RGB) -> Self {
-        Self {
-            rsp: cs.to_rgb_coeffs(rgb),
-        }
-    }
-
-    pub fn sample_at(&self, lambda: Float) -> Float {
-        self.rsp.evaluate(lambda)
-    }
-
-    pub fn max_value(&self) -> Float {
-        self.rsp.max_value()
-    }
-
-    fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
-        let mut s = SampledSpectrum::default();
-        for i in 0..N_SPECTRUM_SAMPLES {
-            s[i] = self.rsp.evaluate(lambda[i]);
-        }
-        s
-    }
-}
-
-#[derive(Debug, Clone, Copy)]
-pub struct UnboundedRGBSpectrum {
-    scale: Float,
-    rsp: RGBSigmoidPolynomial,
-}
-
-impl UnboundedRGBSpectrum {
-    pub fn new(cs: RGBColorspace, rgb: RGB) -> Self {
-        let m = rgb.r.max(rgb.g).max(rgb.b);
-        let scale = 2.0 * m;
-        let scaled_rgb = if scale != 0.0 {
-            rgb / scale
-        } else {
-            RGB::new(0.0, 0.0, 0.0)
-        };
-        Self {
-            scale,
-            rsp: cs.to_rgb_coeffs(scaled_rgb),
-        }
-    }
-    pub fn sample_at(&self, lambda: Float) -> Float {
-        self.scale * self.rsp.evaluate(lambda)
-    }
-
-    pub fn max_value(&self) -> Float {
-        self.scale * self.rsp.max_value()
-    }
-}
-
-#[derive(Debug, Clone, Default)]
-pub struct RGBIlluminantSpectrum {
-    scale: Float,
-    rsp: RGBSigmoidPolynomial,
-    illuminant: Option<Arc<DenselySampledSpectrum>>,
-}
-
-impl RGBIlluminantSpectrum {
-    pub fn sample_at(&self, lambda: Float) -> Float {
-        match &self.illuminant {
-            Some(illuminant) => {
-                self.scale * self.rsp.evaluate(lambda) * illuminant.sample_at(lambda)
-            }
-            None => 0.0,
-        }
-    }
-
-    pub fn max_value(&self) -> Float {
-        match &self.illuminant {
-            Some(illuminant) => self.scale * self.rsp.max_value() * illuminant.max_value(),
-            None => 0.0,
-        }
-    }
-}
-
-#[derive(Debug, Clone)]
-pub struct DenselySampledSpectrum {
-    lambda_min: i32,
-    lambda_max: i32,
-    values: Vec<Float>,
-}
-
-impl DenselySampledSpectrum {
-    pub fn new(lambda_min: i32, lambda_max: i32) -> Self {
-        let n_values = (lambda_max - lambda_min + 1).max(0) as usize;
-        Self {
-            lambda_min,
-            lambda_max,
-            values: vec![0.0; n_values],
-        }
-    }
-
-    pub fn from_spectrum(spec: &Spectrum, lambda_min: i32, lambda_max: i32) -> Self {
-        let mut s = Self::new(lambda_min, lambda_max);
-        if s.values.is_empty() {
-            return s;
-        }
-        for lambda in lambda_min..=lambda_max {
-            let index = (lambda - lambda_min) as usize;
-            s.values[index] = spec.sample_at(lambda as Float);
-        }
-        s
-    }
-
-    pub fn from_function<F>(f: F, lambda_min: i32, lambda_max: i32) -> Self
-    where
-        F: Fn(Float) -> Float,
-    {
-        let mut s = Self::new(lambda_min, lambda_max);
-        if s.values.is_empty() {
-            return s;
-        }
-        for lambda in lambda_min..=lambda_max {
-            let index = (lambda - lambda_min) as usize;
-            s.values[index] = f(lambda as Float);
-        }
-        s
-    }
-
-    pub fn sample_at(&self, lambda: Float) -> Float {
-        let offset = (lambda.round() as i32) - self.lambda_min;
-        if offset < 0 || offset as usize >= self.values.len() {
-            0.0
-        } else {
-            self.values[offset as usize]
-        }
-    }
-
-    pub fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
-        let mut s = SampledSpectrum::default();
-
-        for i in 0..N_SPECTRUM_SAMPLES {
-            let offset = lambda[i].round() as usize - LAMBDA_MIN as usize;
-            if offset >= self.values.len() {
-                s[i] = 0.;
-            } else {
-                s[i] = self.values[offset];
-            }
-        }
-        s
-    }
-
-    pub fn max_value(&self) -> Float {
-        self.values.iter().fold(Float::MIN, |a, b| a.max(*b))
-    }
-
-    pub fn min_component_value(&self) -> Float {
-        self.values.iter().fold(Float::INFINITY, |a, &b| a.min(b))
-    }
-
-    pub fn max_component_value(&self) -> Float {
-        self.values
-            .iter()
-            .fold(Float::NEG_INFINITY, |a, &b| a.max(b))
-    }
-
-    pub fn average(&self) -> Float {
-        self.values.iter().sum::<Float>() / (N_SPECTRUM_SAMPLES as Float)
-    }
-
-    pub fn safe_div(&self, rhs: SampledSpectrum) -> Self {
-        let mut r = Self::new(1, 1);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            r.values[i] = if rhs[i] != 0.0 {
-                self.values[i] / rhs.values[i]
-            } else {
-                0.0
-            }
-        }
-        r
-    }
-}
-
-#[derive(Debug, Clone)]
-pub struct PiecewiseLinearSpectrum {
-    samples: Vec<(Float, Float)>,
-}
-
-impl PiecewiseLinearSpectrum {
-    pub fn from_interleaved(data: &[Float]) -> Self {
-        assert!(
-            data.len() % 2 == 0,
-            "Interleaved data must have an even number of elements"
-        );
-        let mut samples = Vec::new();
-        for pair in data.chunks(2) {
-            samples.push((pair[0], pair[1]));
-        }
-        // PBRT requires the samples to be sorted by wavelength for interpolation.
-        samples.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
-        Self { samples }
-    }
-
-    fn sample_at(&self, lambda: Float) -> Float {
-        if self.samples.is_empty() {
-            return 0.0;
-        }
-        // Handle boundary conditions
-        if lambda <= self.samples[0].0 {
-            return self.samples[0].1;
-        }
-        if lambda >= self.samples.last().unwrap().0 {
-            return self.samples.last().unwrap().1;
-        }
-
-        let i = self.samples.partition_point(|s| s.0 < lambda);
-        let s1 = self.samples[i - 1];
-        let s2 = self.samples[i];
-
-        let t = (lambda - s1.0) / (s2.0 - s1.0);
-        (1.0 - t) * s1.1 + t * s2.1
-    }
-}
-
-#[derive(Debug, Clone, Copy)]
-pub struct BlackbodySpectrum {
-    temperature: Float,
-    normalization_factor: Float,
-}
-
-// Planck's Law
-impl BlackbodySpectrum {
-    const C: Float = 299792458.0;
-    const H: Float = 6.62606957e-34;
-    const KB: Float = 1.3806488e-23;
-
-    pub fn new(temperature: Float) -> Self {
-        let lambda_max = 2.8977721e-3 / temperature * 1e9;
-        let max_val = Self::planck_law(lambda_max, temperature);
-        Self {
-            temperature,
-            normalization_factor: if max_val > 0.0 { 1.0 / max_val } else { 0.0 },
-        }
-    }
-
-    fn planck_law(lambda_nm: Float, temp: Float) -> Float {
-        if temp <= 0.0 {
-            return 0.0;
-        }
-        let lambda_m = lambda_nm * 1e-9; // Convert nm to meters
-        let c1 = 2.0 * Self::H * Self::C * Self::C;
-        let c2 = (Self::H * Self::C) / Self::KB;
-
-        let numerator = c1 / lambda_m.powi(5);
-        let denominator = (c2 / (lambda_m * temp)).exp() - 1.0;
-
-        if denominator.is_infinite() {
-            0.0
-        } else {
-            numerator / denominator
-        }
-    }
-
-    fn sample_at(&self, lambda: Float) -> Float {
-        Self::planck_law(lambda, self.temperature) * self.normalization_factor
-    }
-}
-
-#[derive(Debug, Clone, Copy)]
-pub struct RGBSpectrum {
-    pub c: [Float; 3],
-}
-
-#[derive(Debug, Clone, Copy)]
-pub struct RGBUnboundedSpectrum {
-    scale: Float,
-    rsp: RGBSigmoidPolynomial,
-}
-
-impl Default for RGBUnboundedSpectrum {
-    fn default() -> Self {
-        Self {
-            scale: 0.0,
-            rsp: RGBSigmoidPolynomial::default(),
-        }
-    }
-}
-
-impl RGBUnboundedSpectrum {
-    pub fn new(cs: &RGBColorspace, rgb: RGB) -> Self {
-        let m = rgb.max();
-
-        let scale = 2.0 * m;
-
-        let rgb_norm = if scale > 0.0 {
-            rgb / scale
-        } else {
-            RGB::new(0.0, 0.0, 0.0)
-        };
-
-        let rsp = cs.to_rgb_coeffs(rgb_norm);
-
-        Self { scale, rsp }
-    }
-    pub fn evaluate(&self, lambda: Float) -> Float {
-        self.scale * self.rsp.evaluate(lambda)
-    }
-
-    pub fn max_value(&self) -> Float {
-        self.scale * self.rsp.max_value()
-    }
-
-    pub fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
-        let mut s = SampledSpectrum::default();
-        for i in 0..N_SPECTRUM_SAMPLES {
-            s[i] = self.scale * self.rsp.evaluate(lambda[i]);
-        }
-        s
-    }
-}
-
-pub mod spectra {
-    use super::*;
-    pub static X: Lazy<Spectrum> = Lazy::new(|| {
-        let pls = PiecewiseLinearSpectrum::from_interleaved(&cie::CIE_X);
-        let dss = DenselySampledSpectrum::from_spectrum(
-            &Spectrum::PiecewiseLinear(pls),
-            LAMBDA_MIN,
-            LAMBDA_MAX,
-        );
-        Spectrum::DenselySampled(dss)
-    });
-    pub static Y: Lazy<Spectrum> = Lazy::new(|| {
-        let pls = PiecewiseLinearSpectrum::from_interleaved(&cie::CIE_Y);
-        let dss = DenselySampledSpectrum::from_spectrum(
-            &Spectrum::PiecewiseLinear(pls),
-            LAMBDA_MIN,
-            LAMBDA_MAX,
-        );
-        Spectrum::DenselySampled(dss)
-    });
-    pub static Z: Lazy<Spectrum> = Lazy::new(|| {
-        let pls = PiecewiseLinearSpectrum::from_interleaved(&cie::CIE_Z);
-        let dss = DenselySampledSpectrum::from_spectrum(
-            &Spectrum::PiecewiseLinear(pls),
-            LAMBDA_MIN,
-            LAMBDA_MAX,
-        );
-        Spectrum::DenselySampled(dss)
-    });
-}

src/utils/splines.rs

@@ -65,12 +65,11 @@ pub fn evaluate_cubic_bezier(cp: &[Point3f], u: Float) -> (Point3f, Vector3f) {
         lerp(u, cp[2], cp[3]),
     ];
     let cp2 = [lerp(u, cp1[0], cp1[1]), lerp(u, cp1[1], cp1[2])];
-    let deriv: Vector3f;
+    let deriv = if (cp2[1] - cp2[0]).norm_squared() > 0. {
-    if (cp2[1] - cp2[0]).norm_squared() > 0. {
+        (cp2[1] - cp2[0]) * 3.
-        deriv = (cp2[1] - cp2[0]) * 3.;
     } else {
-        deriv = cp[3] - cp[0]
+        cp[3] - cp[0]
-    }
+    };
 
     (lerp(u, cp2[0], cp2[1]), deriv)
 }

src/utils/transform.rs

@@ -1,13 +1,16 @@
 use num_traits::Float as NumFloat;
 use std::error::Error;
 use std::fmt::{self, Display};
+use std::iter::{Product, Sum};
 use std::ops::{Add, Div, Index, IndexMut, Mul};
 use std::sync::Arc;
 
 use super::color::{RGB, XYZ};
 use super::math::{SquareMatrix, safe_acos};
 use super::quaternion::Quaternion;
-use crate::core::interaction::{SurfaceInteraction, MediumInteraction, Interaction, InteractionData};
+use crate::core::interaction::{
+    Interaction, InteractionData, MediumInteraction, SurfaceInteraction,
+};
 use crate::core::pbrt::{Float, gamma};
 use crate::geometry::{
     Bounds3f, Normal, Normal3f, Point, Point3f, Point3fi, Ray, Vector, Vector3f, Vector3fi,
@@ -21,7 +24,7 @@ pub struct Transform<T: NumFloat> {
     m_inv: SquareMatrix<T, 4>,
 }
 
-impl<T: NumFloat> Transform<T> {
+impl<T: NumFloat + Sum + Product> Transform<T> {
     pub fn new(m: SquareMatrix<T, 4>, m_inv: SquareMatrix<T, 4>) -> Self {
         Self { m, m_inv }
     }
@@ -33,8 +36,7 @@ impl<T: NumFloat> Transform<T> {
 
     pub fn identity() -> Self {
         let m: SquareMatrix<T, 4> = SquareMatrix::identity();
-        let m_inv = m.clone();
+        Self { m, m_inv: m }
-        Self { m, m_inv }
     }
 
     pub fn is_identity(&self) -> bool {
@@ -43,26 +45,28 @@ impl<T: NumFloat> Transform<T> {
 
     pub fn inverse(&self) -> Self {
         Self {
-            m: self.m_inv.clone(),
+            m: self.m_inv,
-            m_inv: self.m.clone(),
+            m_inv: self.m,
         }
     }
 
     pub fn apply_inverse(&self, p: Point<T, 3>) -> Point<T, 3> {
-        self.clone() * p
+        *self * p
     }
 
     pub fn apply_inverse_vector(&self, v: Vector<T, 3>) -> Vector<T, 3> {
         let x = v.x();
         let y = v.y();
         let z = v.z();
-        Vector::<T, 3>::new(self.m_inv[0][0] * x + self.m_inv[0][1] * y + self.m_inv[0][2] * z,
+        Vector::<T, 3>::new(
-                          self.m_inv[1][0] * x + self.m_inv[1][1] * y + self.m_inv[1][2] * z,
+            self.m_inv[0][0] * x + self.m_inv[0][1] * y + self.m_inv[0][2] * z,
-                          self.m_inv[2][0] * x + self.m_inv[2][1] * y + self.m_inv[2][2] * z)
+            self.m_inv[1][0] * x + self.m_inv[1][1] * y + self.m_inv[1][2] * z,
-    }    
+            self.m_inv[2][0] * x + self.m_inv[2][1] * y + self.m_inv[2][2] * z,
+        )
+    }
 
     pub fn apply_inverse_normal(&self, n: Normal<T, 3>) -> Normal<T, 3> {
-        self.clone() * n
+        *self * n
     }
 
     pub fn swaps_handedness(&self) -> bool {
@@ -71,7 +75,7 @@ impl<T: NumFloat> Transform<T> {
             [self.m[1][0], self.m[1][1], self.m[1][2]],
             [self.m[2][0], self.m[2][1], self.m[2][2]],
         ]);
-        return s.determinant() < T::zero();
+        s.determinant() < T::zero()
     }
 }
 
@@ -86,9 +90,9 @@ impl Transform<Float> {
         let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z + self.m[3][3];
 
         if wp == 1. {
-            return Point3f::new(xp, yp, zp);
+            Point3f::new(xp, yp, zp)
         } else {
-            return Point3f::new(xp / wp, yp / wp, zp / wp);
+            Point3f::new(xp / wp, yp / wp, zp / wp)
         }
     }
 
@@ -102,9 +106,9 @@ impl Transform<Float> {
         let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z;
 
         if wp == 1. {
-            return Vector3f::new(xp, yp, zp);
+            Vector3f::new(xp, yp, zp)
         } else {
-            return Vector3f::new(xp / wp, yp / wp, zp / wp);
+            Vector3f::new(xp / wp, yp / wp, zp / wp)
         }
     }
 
@@ -118,9 +122,9 @@ impl Transform<Float> {
         let wp = self.m[3][0] * x + self.m[3][1] * y + self.m[3][2] * z;
 
         if wp == 1. {
-            return Normal3f::new(xp, yp, zp);
+            Normal3f::new(xp, yp, zp)
         } else {
-            return Normal3f::new(xp / wp, yp / wp, zp / wp);
+            Normal3f::new(xp / wp, yp / wp, zp / wp)
         }
     }
 
@@ -178,9 +182,9 @@ impl Transform<Float> {
         }
 
         if wp == 1. {
-            return Point3fi::new_with_error(Point([xp, yp, zp]), p_error);
+            Point3fi::new_with_error(Point([xp, yp, zp]), p_error)
         } else {
-            return Point3fi::new_with_error(Point([xp / wp, yp / wp, zp / wp]), p_error);
+            Point3fi::new_with_error(Point([xp / wp, yp / wp, zp / wp]), p_error)
         }
     }
 
@@ -209,9 +213,9 @@ impl Transform<Float> {
         }
 
         if wp == 1. {
-            return Vector3fi::new_with_error(Vector3f::new(xp, yp, zp), v_error);
+            Vector3fi::new_with_error(Vector3f::new(xp, yp, zp), v_error)
         } else {
-            return Vector3fi::new_with_error(Vector3f::new(xp / wp, yp / wp, zp / wp), v_error);
+            Vector3fi::new_with_error(Vector3f::new(xp / wp, yp / wp, zp / wp), v_error)
         }
     }
 
@@ -240,19 +244,19 @@ impl Transform<Float> {
 
             Box::new(ret)
         } else if let Some(mi) = inter.as_any().downcast_ref::<MediumInteraction>() {
-        let ret = MediumInteraction {
+            let ret = MediumInteraction {
-            common: InteractionData {
+                common: InteractionData {
-                pi: self.apply_to_interval(&mi.common.pi),
+                    pi: self.apply_to_interval(&mi.common.pi),
-                n: self.apply_to_normal(mi.common.n).normalize(),
+                    n: self.apply_to_normal(mi.common.n).normalize(),
-                wo: self.apply_to_vector(mi.common.wo).normalize(),
+                    wo: self.apply_to_vector(mi.common.wo).normalize(),
-                time: mi.common.time,
+                    time: mi.common.time,
-                medium_interface: mi.common.medium_interface.clone(),
+                    medium_interface: mi.common.medium_interface.clone(),
-                medium: mi.common.medium.clone(),
+                    medium: mi.common.medium.clone(),
-            },
+                },
-            medium: mi.medium.clone(),
+                medium: mi.medium.clone(),
-            phase: mi.phase.clone(),
+                phase: mi.phase.clone(),
-        };
+            };
-        Box::new(ret)
+            Box::new(ret)
         } else {
             panic!("Unhandled Interaction type in Transform::apply_to_interaction");
         }
@@ -271,44 +275,50 @@ impl Transform<Float> {
         let wp = (m_inv[3][0] * x + m_inv[3][1] * y) + (m_inv[3][2] * z + m_inv[3][3]);
 
         // Compute absolute error for transformed point
-        let p_out_error: Vector3f;
         let g3 = gamma(3);
 
-        if p.is_exact() {
+        let p_out_error = if p.is_exact() {
-            p_out_error = Vector3f::new(
+            Vector3f::new(
                 g3 * (m_inv[0][0] * x).abs() + (m_inv[0][1] * y).abs() + (m_inv[0][2] * z).abs(),
                 g3 * (m_inv[1][0] * x).abs() + (m_inv[1][1] * y).abs() + (m_inv[1][2] * z).abs(),
                 g3 * (m_inv[2][0] * x).abs() + (m_inv[2][1] * y).abs() + (m_inv[2][2] * z).abs(),
-            );
+            )
         } else {
             let p_in_error = p.error();
             let g3_plus_1 = g3 + 1.0;
-            
+
-            p_out_error = Vector3f::new(
+            Vector3f::new(
-                g3_plus_1 * (m_inv[0][0].abs() * p_in_error.x() +
+                g3_plus_1
-                             m_inv[0][1].abs() * p_in_error.y() +
+                    * (m_inv[0][0].abs() * p_in_error.x()
-                             m_inv[0][2].abs() * p_in_error.z()) +
+                        + m_inv[0][1].abs() * p_in_error.y()
-                g3 * ((m_inv[0][0] * x).abs() + (m_inv[0][1] * y).abs() +
+                        + m_inv[0][2].abs() * p_in_error.z())
-                      (m_inv[0][2] * z).abs() + m_inv[0][3].abs()),
+                    + g3 * ((m_inv[0][0] * x).abs()
-                
+                        + (m_inv[0][1] * y).abs()
-                g3_plus_1 * (m_inv[1][0].abs() * p_in_error.x() +
+                        + (m_inv[0][2] * z).abs()
-                             m_inv[1][1].abs() * p_in_error.y() +
+                        + m_inv[0][3].abs()),
-                             m_inv[1][2].abs() * p_in_error.z()) +
+                g3_plus_1
-                g3 * ((m_inv[1][0] * x).abs() + (m_inv[1][1] * y).abs() +
+                    * (m_inv[1][0].abs() * p_in_error.x()
-                      (m_inv[1][2] * z).abs() + m_inv[1][3].abs()),
+                        + m_inv[1][1].abs() * p_in_error.y()
-                
+                        + m_inv[1][2].abs() * p_in_error.z())
-                g3_plus_1 * (m_inv[2][0].abs() * p_in_error.x() +
+                    + g3 * ((m_inv[1][0] * x).abs()
-                             m_inv[2][1].abs() * p_in_error.y() +
+                        + (m_inv[1][1] * y).abs()
-                             m_inv[2][2].abs() * p_in_error.z()) +
+                        + (m_inv[1][2] * z).abs()
-                g3 * ((m_inv[2][0] * x).abs() + (m_inv[2][1] * y).abs() +
+                        + m_inv[1][3].abs()),
-                      (m_inv[2][2] * z).abs() + m_inv[2][3].abs()),
+                g3_plus_1
-            );
+                    * (m_inv[2][0].abs() * p_in_error.x()
-        }
+                        + m_inv[2][1].abs() * p_in_error.y()
+                        + m_inv[2][2].abs() * p_in_error.z())
+                    + g3 * ((m_inv[2][0] * x).abs()
+                        + (m_inv[2][1] * y).abs()
+                        + (m_inv[2][2] * z).abs()
+                        + m_inv[2][3].abs()),
+            )
+        };
 
         if wp == 1.0 {
             Point3fi::new_with_error(Point3f::new(xp, yp, zp), p_out_error)
         } else {
-            Point3fi::new_with_error(Point3f::new(xp/wp, yp/wp, zp/wp), p_out_error)
+            Point3fi::new_with_error(Point3f::new(xp / wp, yp / wp, zp / wp), p_out_error)
         }
     }
 
@@ -326,10 +336,13 @@ impl Transform<Float> {
                 t = t_max - dt;
             }
         }
-        (Ray::new(Point3f::from(o), d, Some(r.time), r.medium.clone()), t)
+        (
+            Ray::new(Point3f::from(o), d, Some(r.time), r.medium.clone()),
+            t,
+        )
     }
 
-    pub fn to_quaternion(&self) -> Quaternion {
+    pub fn to_quaternion(self) -> Quaternion {
         let trace = self.m.trace();
         let mut quat = Quaternion::default();
         if trace > 0. {
@@ -430,7 +443,7 @@ impl Transform<Float> {
         m[2][1] = a.y() * a.z() * (1. - cos_theta) + a.x() * sin_theta;
         m[2][2] = a.z() * a.z() + (1. - a.z() * a.z()) * cos_theta;
         m[2][3] = 0.;
-        Transform::new(m.clone(), m.transpose())
+        Transform::new(m, m.transpose())
     }
 
     pub fn rotate_from_to(from: Vector3f, to: Vector3f) -> Self {
@@ -457,7 +470,7 @@ impl Transform<Float> {
                     + 4. * uv / (uu * vv) * v[i] * u[j];
             }
         }
-        Transform::new(r.clone(), r.transpose())
+        Transform::new(r, r.transpose())
     }
 
     pub fn rotate_x(theta: Float) -> Self {
@@ -470,7 +483,7 @@ impl Transform<Float> {
             [0., 0., 0., 1.],
         ]);
         Self {
-            m: m.clone(),
+            m,
             m_inv: m.transpose(),
         }
     }
@@ -485,7 +498,7 @@ impl Transform<Float> {
             [0., 0., 0., 1.],
         ]);
         Self {
-            m: m.clone(),
+            m,
             m_inv: m.transpose(),
         }
     }
@@ -500,14 +513,14 @@ impl Transform<Float> {
             [0., 0., 0., 1.],
         ]);
         Self {
-            m: m.clone(),
+            m,
             m_inv: m.transpose(),
         }
     }
 
     pub fn decompose(&self) -> (Vector3f, SquareMatrix<Float, 4>, SquareMatrix<Float, 4>) {
         let t = Vector3f::new(self.m[0][3], self.m[1][3], self.m[2][3]);
-        let mut m = self.m.clone();
+        let mut m = self.m;
         for i in 0..3 {
             m[i][3] = 0.;
             m[3][i] = m[i][3];
@@ -516,14 +529,14 @@ impl Transform<Float> {
         m[3][3] = 1.;
 
         let mut norm: Float;
-        let mut r = m.clone();
+        let mut r = m;
         let mut count = 0;
         loop {
             let rit = r
                 .transpose()
                 .inverse()
                 .expect("Transform is not decomposable");
-            let rnext = (r.clone() + rit.clone()) / 2.;
+            let rnext = (r + rit) / 2.;
             norm = 0.0;
             for i in 0..3 {
                 let n = (r[i][0] - rnext[i][0]).abs()
@@ -563,15 +576,15 @@ impl<T: NumFloat> Mul for Transform<T> {
     }
 }
 
-impl<'a, 'b, T> Mul<&'b Transform<T>> for &'a Transform<T>
+impl<'b, T> Mul<&'b Transform<T>> for &Transform<T>
 where
     T: NumFloat,
 {
     type Output = Transform<T>;
     fn mul(self, rhs: &'b Transform<T>) -> Self::Output {
         Transform {
-            m: self.m.clone() * rhs.m.clone(),
+            m: self.m * rhs.m,
-            m_inv: rhs.m_inv.clone() * self.m_inv.clone(),
+            m_inv: rhs.m_inv * self.m_inv,
         }
     }
 }
@@ -730,8 +743,8 @@ impl AnimatedTransform {
 
         if !actually_animated {
             return Self {
-                start_transform: start_transform.clone(),
+                start_transform: *start_transform,
-                end_transform: end_transform.clone(),
+                end_transform: *end_transform,
                 start_time,
                 end_time,
                 actually_animated: false,
@@ -1874,7 +1887,7 @@ impl AnimatedTransform {
                 std::array::from_fn(|_| DerivativeTerm::default()),
             )
         };
-        return AnimatedTransform {
+        AnimatedTransform {
             start_transform: *start_transform,
             end_transform: *end_transform,
             start_time,
@@ -1889,7 +1902,7 @@ impl AnimatedTransform {
             c3,
             c4,
             c5,
-        };
+        }
     }
 
     pub fn apply_point(&self, p: Point3f, time: Float) -> Point3f {
@@ -1924,7 +1937,7 @@ impl AnimatedTransform {
 
     pub fn apply_interaction(&self, si: &dyn Interaction) -> Box<dyn Interaction> {
         if !self.actually_animated {
-            return self.start_transform.apply_to_interaction(si)
+            return self.start_transform.apply_to_interaction(si);
         }
         let t = self.interpolate(si.time());
         t.apply_to_interaction(si)
@@ -1947,35 +1960,37 @@ impl AnimatedTransform {
         let scale_transform =
             Transform::from_matrix(scale).expect("Scale matrix is not inversible");
 
-        return Transform::translate(trans) * Transform::from(rotate) * scale_transform;
+        Transform::translate(trans) * Transform::from(rotate) * scale_transform
     }
 
     pub fn apply_inverse_point(&self, p: Point3f, time: Float) -> Point3f {
         if !self.actually_animated {
             return self.start_transform.apply_inverse(p);
         }
-        return self.interpolate(time).apply_inverse(p);
+        self.interpolate(time).apply_inverse(p)
     }
 
     pub fn apply_inverse_vector(&self, v: Vector3f, time: Float) -> Vector3f {
         if !self.actually_animated {
             return self.start_transform.apply_inverse_vector(v);
         }
-        return self.interpolate(time).apply_inverse_vector(v);
+        self.interpolate(time).apply_inverse_vector(v)
     }
 
     pub fn apply_inverse_normal(&self, n: Normal3f, time: Float) -> Normal3f {
         if !self.actually_animated {
             return self.start_transform.apply_inverse_normal(n);
         }
-        return self.interpolate(time).apply_inverse_normal(n);
+        self.interpolate(time).apply_inverse_normal(n)
     }
 
     pub fn motion_bounds(&self, b: &Bounds3f) -> Bounds3f {
         if !self.actually_animated {
-            return self.start_transform.apply_to_bounds(*b)
+            return self.start_transform.apply_to_bounds(*b);
         }
-        return self.start_transform.apply_to_bounds(*b).union(self.end_transform.apply_to_bounds(*b))
+        self.start_transform
+            .apply_to_bounds(*b)
+            .union(self.end_transform.apply_to_bounds(*b))
     }
 }