diff --git a/Cargo.toml b/Cargo.toml
index e61eccf..6758722 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -7,6 +7,7 @@ edition = "2024"
 anyhow = "1.0.100"
 bitflags = "2.10.0"
 bumpalo = "3.19.0"
+enum_dispatch = "0.3.13"
 exr = "1.73.0"
 half = "2.7.1"
 image = "0.25.8"
@@ -19,3 +20,7 @@ rand = "0.9.2"
 rayon = "1.11.0"
 smallvec = "1.15.1"
 thiserror = "2.0.17"
+
+[features]
+default = []
+use_f64 = []
diff --git a/src/core/bxdf.rs b/src/core/bxdf.rs
index f4098de..22f2e8d 100644
--- a/src/core/bxdf.rs
+++ b/src/core/bxdf.rs
@@ -5,21 +5,30 @@ use std::fmt;
 use std::ops::Not;
 use std::sync::{Arc, RwLock};
 
+use enum_dispatch::enum_dispatch;
+
 use crate::core::pbrt::{Float, 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::utils::math::square;
+use crate::utils::color::RGB;
+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,
+};
 use crate::utils::sampling::{
-    PiecewiseLinear2D, cosine_hemisphere_pdf, sample_cosine_hemisphere, sample_uniform_hemisphere,
-    uniform_hemisphere_pdf,
+    PiecewiseLinear2D, cosine_hemisphere_pdf, sample_cosine_hemisphere, sample_trimmed_logistic,
+    sample_uniform_hemisphere, uniform_hemisphere_pdf,
 };
 use crate::utils::scattering::{
     TrowbridgeReitzDistribution, fr_complex, fr_complex_from_spectrum, fr_dielectric, reflect,
     refract,
 };
-use crate::utils::spectrum::{N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths};
+use crate::utils::spectrum::{
+    N_SPECTRUM_SAMPLES, RGBUnboundedSpectrum, SampledSpectrum, SampledWavelengths,
+};
 
 bitflags! {
     #[derive(Debug, Clone, Copy, PartialEq, Eq)]
@@ -194,8 +203,126 @@ pub struct CoatedDiffuseBxDF;
 #[derive(Debug, Clone)]
 pub struct CoatedConductorBxDF;
 
+static P_MAX: usize = 3;
 #[derive(Debug, Clone)]
-pub struct HairBxDF;
+pub struct HairBxDF {
+    h: Float,
+    eta: Float,
+    sigma_a: SampledSpectrum,
+    beta_m: Float,
+    beta_n: Float,
+    v: [Float; P_MAX + 1],
+    s: Float,
+    sin_2k_alpha: [Float; P_MAX],
+    cos_2k_alpha: [Float; P_MAX],
+}
+
+impl HairBxDF {
+    pub fn new(
+        h: Float,
+        eta: Float,
+        sigma_a: SampledSpectrum,
+        beta_m: Float,
+        beta_n: Float,
+        alpha: Float,
+    ) -> Self {
+        let mut sin_2k_alpha = [0.; P_MAX];
+        let mut cos_2k_alpha = [0.; P_MAX];
+        sin_2k_alpha[0] = radians(alpha).sin();
+        cos_2k_alpha[0] = safe_sqrt(1. - square(sin_2k_alpha[0]));
+
+        for i in 0..P_MAX {
+            sin_2k_alpha[i] = 2. * cos_2k_alpha[i - 1] * sin_2k_alpha[i - 1];
+            cos_2k_alpha[i] = square(cos_2k_alpha[i - 1]) - square(sin_2k_alpha[i - 1]);
+        }
+
+        Self {
+            h,
+            eta,
+            sigma_a,
+            beta_m,
+            beta_n,
+            v: [0.; P_MAX + 1],
+            s: 0.,
+            sin_2k_alpha,
+            cos_2k_alpha,
+        }
+    }
+
+    fn ap(
+        cos_theta_o: Float,
+        eta: Float,
+        h: Float,
+        t: SampledSpectrum,
+    ) -> [SampledSpectrum; P_MAX + 1] {
+        let cos_gamma_o = safe_sqrt(1. - square(h));
+        let cos_theta = cos_theta_o * cos_gamma_o;
+        let f = fr_dielectric(cos_theta, eta);
+        let ap0 = SampledSpectrum::new(f);
+        let ap1 = t * (1.0 - f).powi(2);
+        let tf = t * f;
+        std::array::from_fn(|p| match p {
+            0 => ap0,
+            1 => ap1,
+            _ if p < P_MAX => ap1 * tf.pow_int(p - 1),
+            _ => ap1 * tf.pow_int(p - 1) / (SampledSpectrum::new(1.0) - tf),
+        })
+    }
+
+    fn mp(
+        cos_theta_i: Float,
+        cos_theta_o: Float,
+        sin_theta_i: Float,
+        sin_theta_o: Float,
+        v: Float,
+    ) -> Float {
+        let a = cos_theta_i * cos_theta_o / v;
+        let b = sin_theta_i * sin_theta_o / v;
+        if v <= 0.1 {
+            fast_exp(log_i0(a) - b - 1. / v + 0.6931 + (1. / (2. * v).ln()))
+        } else {
+            fast_exp(-b) * i0(a) / ((1. / v).sinh() * 2. * v)
+        }
+    }
+
+    fn np(phi: Float, p: i32, s: Float, gamma_o: Float, gamma_t: Float) -> Float {
+        let mut dphi = phi - Self::phi(p, gamma_o, gamma_t);
+        while dphi > PI {
+            dphi -= 2. * PI;
+        }
+        while dphi < -PI {
+            dphi += 2. * PI;
+        }
+
+        trimmed_logistic(dphi, s, -PI, PI)
+    }
+
+    fn phi(p: i32, gamma_o: Float, gamma_t: Float) -> Float {
+        2. * p as Float * gamma_t - 2. * gamma_o + p as Float * PI
+    }
+
+    fn ap_pdf(&self, cos_theta_o: Float) -> [Float; P_MAX + 1] {
+        let sin_theta_o = safe_sqrt(1. - square(cos_theta_o));
+        let sin_theta_t = sin_theta_o / self.eta;
+        let cos_theta_t = safe_sqrt(1. - square(sin_theta_t));
+        let etap = safe_sqrt(square(self.eta) - square(sin_theta_o)) / cos_theta_o;
+        let sin_gamma_t = self.h / etap;
+        let cos_gamma_t = safe_sqrt(1. - square(sin_gamma_t));
+        // let gamma_t = safe_asin(sin_gamma_t);
+        let t_value = -self.sigma_a * (2. * cos_gamma_t / cos_theta_t);
+        let t = t_value.exp();
+        let ap = Self::ap(cos_theta_o, self.eta, self.h, t);
+        let sum_y: Float = ap.iter().map(|s| s.average()).sum();
+        std::array::from_fn(|i| ap[i].average() / sum_y)
+    }
+
+    fn sigma_a_from_concentration(ce: Float, cp: Float) -> RGBUnboundedSpectrum {
+        let eumelanin_sigma_a = RGB::new(0.419, 0.697, 1.37);
+        let pheomelanin_sigma_a = RGB::new(0.187, 0.4, 1.05);
+        let sigma_a = ce * eumelanin_sigma_a + cp * pheomelanin_sigma_a;
+        RGBUnboundedSpectrum::new(RGBColorspace::srgb(), sigma_a)
+    }
+}
 
 #[derive(Debug, Clone)]
 pub struct MeasuredBxDFData {
@@ -274,6 +401,7 @@ impl Default for FArgs {
     }
 }
 
+#[enum_dispatch]
 pub trait BxDFTrait: Any + Send + Sync + std::fmt::Debug {
     fn flags(&self) -> BxDFFlags;
 
@@ -344,6 +472,7 @@ impl BxDFTrait for EmptyBxDF {
     }
 }
 
+#[enum_dispatch(BxDFTrait)]
 #[derive(Debug, Clone)]
 pub enum BxDF {
     Diffuse(DiffuseBxDF),
@@ -351,70 +480,8 @@ pub enum BxDF {
     ThinDielectric(ThinDielectricBxDF),
     Conductor(ConductorBxDF),
     Measured(MeasuredBxDF),
+    Hair(HairBxDF),
     // DiffuseTransmission(DiffuseTransmissionBxDF),
-    // Hair(HairBxDF),
-}
-
-impl BxDFTrait for BxDF {
-    fn flags(&self) -> BxDFFlags {
-        match self {
-            BxDF::Diffuse(b) => b.flags(),
-            BxDF::Dielectric(b) => b.flags(),
-            BxDF::ThinDielectric(b) => b.flags(),
-            BxDF::Conductor(b) => b.flags(),
-            BxDF::Measured(b) => b.flags(),
-        }
-    }
-
-    fn f(&self, wo: &Vector3f, wi: &Vector3f, mode: TransportMode) -> SampledSpectrum {
-        match self {
-            BxDF::Diffuse(b) => b.f(wo, wi, mode),
-            BxDF::Dielectric(b) => b.f(wo, wi, mode),
-            BxDF::ThinDielectric(b) => b.f(wo, wi, mode),
-            BxDF::Conductor(b) => b.f(wo, wi, mode),
-            BxDF::Measured(b) => b.f(wo, wi, mode),
-        }
-    }
-
-    fn sample_f(&self, wo: &Vector3f, uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
-        match self {
-            BxDF::Diffuse(b) => b.sample_f(wo, uc, u, f_args),
-            BxDF::Dielectric(b) => b.sample_f(wo, uc, u, f_args),
-            BxDF::ThinDielectric(b) => b.sample_f(wo, uc, u, f_args),
-            BxDF::Conductor(b) => b.sample_f(wo, uc, u, f_args),
-            BxDF::Measured(b) => b.sample_f(wo, uc, u, f_args),
-        }
-    }
-
-    fn pdf(&self, wo: &Vector3f, wi: &Vector3f, f_args: FArgs) -> Float {
-        match self {
-            BxDF::Diffuse(b) => b.pdf(wo, wi, f_args),
-            BxDF::Dielectric(b) => b.pdf(wo, wi, f_args),
-            BxDF::ThinDielectric(b) => b.pdf(wo, wi, f_args),
-            BxDF::Conductor(b) => b.pdf(wo, wi, f_args),
-            BxDF::Measured(b) => b.pdf(wo, wi, f_args),
-        }
-    }
-
-    fn regularize(&mut self) {
-        match self {
-            BxDF::Diffuse(b) => b.regularize(),
-            BxDF::Dielectric(b) => b.regularize(),
-            BxDF::ThinDielectric(b) => b.regularize(),
-            BxDF::Conductor(b) => b.regularize(),
-            BxDF::Measured(b) => b.regularize(),
-        }
-    }
-
-    fn as_any(&self) -> &dyn Any {
-        match self {
-            BxDF::Diffuse(b) => b.as_any(),
-            BxDF::Dielectric(b) => b.as_any(),
-            BxDF::ThinDielectric(b) => b.as_any(),
-            BxDF::Conductor(b) => b.as_any(),
-            BxDF::Measured(b) => b.as_any(),
-        }
-    }
 }
 
 #[derive(Debug, Clone)]
@@ -1246,3 +1313,237 @@ impl BxDFTrait for MeasuredBxDF {
         todo!()
     }
 }
+
+impl BxDFTrait for HairBxDF {
+    fn flags(&self) -> BxDFFlags {
+        BxDFFlags::GLOSSY_REFLECTION
+    }
+
+    fn f(&self, wo: &Vector3f, wi: &Vector3f, _mode: TransportMode) -> SampledSpectrum {
+        // Compute hair coordinate system terms related to wo
+        let sin_theta_o = wo.x();
+        let cos_theta_o = safe_sqrt(1. - square(sin_theta_o));
+        let phi_o = wo.z().atan2(wo.y());
+        let gamma_o = safe_asin(self.h);
+        // Compute hair coordinate system terms related to wi
+        let sin_theta_i = wi.x();
+        let cos_theta_i = safe_sqrt(1. - square(sin_theta_i));
+        let phi_i = wi.z().atan2(wi.y());
+
+        let sin_theta_t = sin_theta_o / self.eta;
+        let cos_theta_t = safe_sqrt(1. - square(sin_theta_t));
+        let etap = safe_sqrt(square(self.eta) - square(sin_theta_o)) / cos_theta_o;
+        let sin_gamma_t = self.h / etap;
+        let cos_gamma_t = safe_sqrt(1. - square(sin_gamma_t));
+        let gamma_t = safe_asin(sin_gamma_t);
+        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 mut f_sum = SampledSpectrum::new(0.);
+        for p in 0..P_MAX {
+            let sin_thetap_o: Float;
+            let cos_thetap_o: Float;
+            if p == 0 {
+                sin_thetap_o =
+                    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];
+            } else if p == 1 {
+                sin_thetap_o =
+                    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];
+            } else if p == 2 {
+                sin_thetap_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,
+                cos_thetap_o,
+                sin_theta_i,
+                sin_thetap_o,
+                self.v[p],
+            ) * ap[p]
+                * Self::np(phi, p as i32, self.s, gamma_o, gamma_t);
+        }
+        if abs_cos_theta(*wi) > 0. {
+            f_sum /= abs_cos_theta(*wi);
+        }
+        f_sum
+    }
+
+    fn sample_f(
+        &self,
+        wo: &Vector3f,
+        mut uc: Float,
+        u: Point2f,
+        f_args: FArgs,
+    ) -> Option<BSDFSample> {
+        let sin_theta_o = wo.x();
+        let cos_theta_o = safe_sqrt(1. - square(sin_theta_o));
+        let phi_o = wo.z().atan2(wo.y());
+        let gamma_o = safe_asin(self.h);
+        // Determine which term to sample for hair scattering
+        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 mut cos_thetap_o: Float;
+        if p == 0 {
+            sin_thetap_o = 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];
+        } else if p == 1 {
+            sin_thetap_o = 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];
+        } else if p == 2 {
+            sin_thetap_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;
+        }
+
+        cos_thetap_o = cos_thetap_o.abs();
+        let cos_theta =
+            1. + self.v[p] * (u[0].max(1e-5) + (1. - u[0]) * fast_exp(-2. / self.v[p])).ln();
+        let sin_theta = safe_sqrt(1. - square(cos_theta));
+        let cos_phi = (2. * PI * u[1]).cos();
+        let sin_theta_i = -cos_theta * sin_thetap_o + sin_theta * cos_phi * cos_thetap_o;
+        let cos_theta_i = safe_sqrt(1. - square(sin_theta_i));
+        let etap = safe_sqrt(square(self.eta) - square(sin_theta_o)) / cos_theta_o;
+        let sin_gamma_t = self.h / etap;
+        // let cos_gamma_t = safe_sqrt(1. - square(sin_gamma_t));
+        let gamma_t = safe_asin(sin_gamma_t);
+        let dphi = if p < P_MAX {
+            Self::phi(p as i32, gamma_o, gamma_t) + sample_trimmed_logistic(uc, self.s, -PI, PI)
+        } else {
+            2. * PI * uc
+        };
+        let phi_i = phi_o + dphi;
+        let wi = Vector3f::new(
+            sin_theta_i,
+            cos_theta_i * phi_i.cos(),
+            cos_theta_i * phi_i.sin(),
+        );
+
+        let mut pdf = 0.;
+        for p in 0..P_MAX {
+            if p == 0 {
+                sin_thetap_o =
+                    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];
+            } else if p == 1 {
+                sin_thetap_o =
+                    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];
+            } else if p == 2 {
+                sin_thetap_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;
+            }
+            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]
+                * Self::np(dphi, p as i32, self.s, gamma_o, gamma_t);
+        }
+        pdf += Self::mp(
+            cos_theta_i,
+            cos_theta_o,
+            sin_theta_i,
+            sin_theta_o,
+            self.v[P_MAX],
+        ) * ap_pdf[P_MAX]
+            * (1. / (2. * PI));
+
+        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)
+    }
+
+    fn pdf(&self, wo: &Vector3f, wi: &Vector3f, _f_args: FArgs) -> Float {
+        let sin_theta_o = wo.x();
+        let cos_theta_o = safe_sqrt(1. - square(sin_theta_o));
+        let phi_o = wo.z().atan2(wo.y());
+        let gamma_o = safe_asin(self.h);
+        // Determine which term to sample for hair scattering
+        let sin_theta_i = wi.x();
+        let cos_theta_i = safe_sqrt(1. - square(sin_theta_i));
+        let phi_i = wi.z().atan2(wi.y());
+        // Compute $\gammat$ for refracted ray
+        let etap = safe_sqrt(self.eta * self.eta - square(sin_theta_o)) / cos_theta_o;
+        let sin_gamma_t = self.h / etap;
+        let gamma_t = safe_asin(sin_gamma_t);
+        // 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 {
+            if p == 0 {
+                sin_thetap_o =
+                    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];
+            } else if p == 1 {
+                sin_thetap_o =
+                    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];
+            } else if p == 2 {
+                sin_thetap_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;
+            }
+            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]
+                * Self::np(phi, p as i32, self.s, gamma_o, gamma_t);
+        }
+
+        pdf += Self::mp(
+            cos_theta_i,
+            cos_theta_o,
+            sin_theta_i,
+            sin_theta_o,
+            self.v[P_MAX],
+        ) * ap_pdf[P_MAX]
+            * (1. / (2. * PI));
+        pdf
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+}
diff --git a/src/core/filter.rs b/src/core/filter.rs
index 8c839fb..e53f304 100644
--- a/src/core/filter.rs
+++ b/src/core/filter.rs
@@ -1,16 +1,18 @@
 use crate::core::pbrt::{Float, lerp};
-use crate::core::sampler::PiecewiseConstant2D;
 use crate::geometry::{Bounds2f, Bounds2i, Point2f, Point2i, Vector2f};
 use crate::utils::containers::Array2D;
 use crate::utils::math::{gaussian, gaussian_integral, sample_tent, windowed_sinc};
+use crate::utils::sampling::PiecewiseConstant2D;
+use enum_dispatch::enum_dispatch;
 
 use rand::Rng;
 use std::hash::Hash;
 
 pub struct FilterSample {
-    p: Point2f,
-    weight: Float,
+    pub p: Point2f,
+    pub weight: Float,
 }
+
 #[derive(Debug)]
 pub struct FilterSampler {
     domain: Bounds2f,
@@ -19,7 +21,7 @@ pub struct FilterSampler {
 }
 
 impl FilterSampler {
-    pub fn new<F>(radius: Vector2f, resolution: Point2i, evaluate_fn: F) -> Self
+    pub fn new<F>(radius: Vector2f, func: F) -> Self
     where
         F: Fn(Point2f) -> Float,
     {
@@ -27,15 +29,18 @@ impl FilterSampler {
             Point2f::new(-radius.x(), -radius.y()),
             Point2f::new(radius.x(), radius.y()),
         );
-        let array_bounds = Bounds2i::from_points(Point2i::new(0, 0), resolution);
-        let mut f = Array2D::new(array_bounds);
-        for j in 0..resolution.y() {
-            for i in 0..resolution.x() {
+
+        let nx = (32.0 * radius.x()) as i32;
+        let ny = (32.0 * radius.y()) as i32;
+
+        let mut f = Array2D::new_with_dims(nx, ny);
+        for y in 0..f.y_size() {
+            for x in 0..f.x_size() {
                 let p = domain.lerp(Point2f::new(
-                    (i as Float + 0.5) / resolution.x() as Float,
-                    (j as Float + 0.5) / resolution.y() as Float,
+                    (x as Float + 0.5) / f.x_size() as Float,
+                    (y as Float + 0.5) / f.y_size() as Float,
                 ));
-                f[Point2i::new(i, j)] = evaluate_fn(p);
+                f[(x as i32, y as i32)] = func(p);
             }
         }
         let distrib = PiecewiseConstant2D::new_with_bounds(&f, domain);
@@ -53,6 +58,7 @@ impl FilterSampler {
     }
 }
 
+#[enum_dispatch]
 pub trait FilterTrait {
     fn radius(&self) -> Vector2f;
     fn evaluate(&self, p: Point2f) -> Float;
@@ -60,6 +66,7 @@ pub trait FilterTrait {
     fn sample(&self, u: Point2f) -> FilterSample;
 }
 
+#[enum_dispatch(FilterTrait)]
 #[derive(Debug)]
 pub enum Filter {
     Box(BoxFilter),
@@ -69,48 +76,6 @@ pub enum Filter {
     Triangle(TriangleFilter),
 }
 
-impl FilterTrait for Filter {
-    fn radius(&self) -> Vector2f {
-        match self {
-            Filter::Box(c) => c.radius(),
-            Filter::Gaussian(c) => c.radius(),
-            Filter::Mitchell(c) => c.radius(),
-            Filter::LanczosSinc(c) => c.radius(),
-            Filter::Triangle(c) => c.radius(),
-        }
-    }
-
-    fn evaluate(&self, p: Point2f) -> Float {
-        match self {
-            Filter::Box(c) => c.evaluate(p),
-            Filter::Gaussian(c) => c.evaluate(p),
-            Filter::Mitchell(c) => c.evaluate(p),
-            Filter::LanczosSinc(c) => c.evaluate(p),
-            Filter::Triangle(c) => c.evaluate(p),
-        }
-    }
-
-    fn integral(&self) -> Float {
-        match self {
-            Filter::Box(c) => c.integral(),
-            Filter::Gaussian(c) => c.integral(),
-            Filter::Mitchell(c) => c.integral(),
-            Filter::LanczosSinc(c) => c.integral(),
-            Filter::Triangle(c) => c.integral(),
-        }
-    }
-
-    fn sample(&self, u: Point2f) -> FilterSample {
-        match self {
-            Filter::Box(c) => c.sample(u),
-            Filter::Gaussian(c) => c.sample(u),
-            Filter::Mitchell(c) => c.sample(u),
-            Filter::LanczosSinc(c) => c.sample(u),
-            Filter::Triangle(c) => c.sample(u),
-        }
-    }
-}
-
 #[derive(Debug)]
 pub struct BoxFilter {
     pub radius: Vector2f,
@@ -162,14 +127,12 @@ impl GaussianFilter {
         let exp_x = gaussian(radius.x(), 0., sigma);
         let exp_y = gaussian(radius.y(), 0., sigma);
 
-        let sampler = FilterSampler::new(
-            radius,
-            Point2i::new((32.0 * radius.x()) as i32, (32.0 * radius.y()) as i32),
-            |p: Point2f| {
-                (gaussian(p.x(), 0., sigma) - exp_x).max(0.)
-                    * (gaussian(p.y(), 0., sigma) - exp_y).max(0.)
-            },
-        );
+        let sampler = FilterSampler::new(radius, move |p: Point2f| {
+            let gx = (gaussian(p.x(), 0., sigma) - exp_x).max(0.0);
+            let gy = (gaussian(p.y(), 0., sigma) - exp_y).max(0.0);
+            gx * gy
+        });
+
         Self {
             radius,
             sigma,
@@ -212,30 +175,12 @@ pub struct MitchellFilter {
 
 impl MitchellFilter {
     pub fn new(radius: Vector2f, b: Float, c: Float) -> Self {
-        let sampler = FilterSampler::new(
-            radius,
-            Point2i::new((32.0 * radius.x()) as i32, (32.0 * radius.y()) as i32),
-            move |p: Point2f| {
-                let mitchell_1d = |x: Float| {
-                    let x = x.abs();
-                    if x <= 1.0 {
-                        ((12.0 - 9.0 * b - 6.0 * c) * x.powi(3)
-                            + (-18.0 + 12.0 * b + 6.0 * c) * x.powi(2)
-                            + (6.0 - 2.0 * b))
-                            * (1.0 / 6.0)
-                    } else if x <= 2.0 {
-                        ((-b - 6.0 * c) * x.powi(3)
-                            + (6.0 * b + 30.0 * c) * x.powi(2)
-                            + (-12.0 * b - 48.0 * c) * x
-                            + (8.0 * b + 24.0 * c))
-                            * (1.0 / 6.0)
-                    } else {
-                        0.0
-                    }
-                };
-                mitchell_1d(2.0 * p.x() / radius.x()) * mitchell_1d(2.0 * p.y() / radius.y())
-            },
-        );
+        let sampler = FilterSampler::new(radius, move |p: Point2f| {
+            let nx = 2.0 * p.x() / radius.x();
+            let ny = 2.0 * p.y() / radius.y();
+            Self::mitchell_1d_eval(b, c, nx) * Self::mitchell_1d_eval(b, c, ny)
+        });
+
         Self {
             radius,
             b,
@@ -244,23 +189,27 @@ impl MitchellFilter {
         }
     }
 
-    fn mitchell_1d(&self, x: Float) -> Float {
+    fn mitchell_1d_eval(b: Float, c: Float, x: Float) -> Float {
         let x = x.abs();
         if x <= 1.0 {
-            ((12.0 - 9.0 * self.b - 6.0 * self.c) * x.powi(3)
-                + (-18.0 + 12.0 * self.b + 6.0 * self.c) * x.powi(2)
-                + (6.0 - 2.0 * self.b))
+            ((12.0 - 9.0 * b - 6.0 * c) * x.powi(3)
+                + (-18.0 + 12.0 * b + 6.0 * c) * x.powi(2)
+                + (6.0 - 2.0 * b))
                 * (1.0 / 6.0)
         } else if x <= 2.0 {
-            ((-self.b - 6.0 * self.c) * x.powi(3)
-                + (6.0 * self.b + 30.0 * self.c) * x.powi(2)
-                + (-12.0 * self.b - 48.0 * self.c) * x
-                + (8.0 * self.b + 24.0 * self.c))
+            ((-b - 6.0 * c) * x.powi(3)
+                + (6.0 * b + 30.0 * c) * x.powi(2)
+                + (-12.0 * b - 48.0 * c) * x
+                + (8.0 * b + 24.0 * c))
                 * (1.0 / 6.0)
         } else {
             0.0
         }
     }
+
+    fn mitchell_1d(&self, x: Float) -> Float {
+        Self::mitchell_1d_eval(self.b, self.c, x)
+    }
 }
 
 impl FilterTrait for MitchellFilter {
@@ -291,13 +240,10 @@ pub struct LanczosSincFilter {
 
 impl LanczosSincFilter {
     pub fn new(radius: Vector2f, tau: Float) -> Self {
-        let sampler = FilterSampler::new(
-            radius,
-            Point2i::new((32.0 * radius.x()) as i32, (32.0 * radius.y()) as i32),
-            move |p: Point2f| {
-                windowed_sinc(p.x(), radius.x(), tau) * windowed_sinc(p.y(), radius.y(), tau)
-            },
-        );
+        let sampler = FilterSampler::new(radius, move |p: Point2f| {
+            windowed_sinc(p.x(), radius.x(), tau) * windowed_sinc(p.y(), radius.y(), tau)
+        });
+
         Self {
             radius,
             tau,
diff --git a/src/core/interaction.rs b/src/core/interaction.rs
index 8cdb239..d8a76b1 100644
--- a/src/core/interaction.rs
+++ b/src/core/interaction.rs
@@ -1,10 +1,14 @@
-use crate::core::bxdf::BSDF;
+use crate::camera::{Camera, CameraTrait};
+use crate::core::bxdf::{BSDF, BxDFFlags};
 use crate::core::material::MaterialTrait;
 use crate::core::medium::{Medium, MediumInterface};
-use crate::core::pbrt::Float;
-use crate::geometry::{Normal3f, Point2f, Point3f, Point3fi, Ray, Vector3f, VectorLike};
+use crate::core::pbrt::{Float, clamp_t};
+use crate::geometry::{
+    Normal3f, Point2f, Point3f, Point3fi, Ray, RayDifferential, Vector3f, VectorLike,
+};
 use crate::lights::Light;
 use crate::shapes::Shape;
+use crate::utils::math::{difference_of_products, square};
 
 use std::any::Any;
 use std::sync::Arc;
@@ -23,6 +27,8 @@ pub trait Interaction: Send + Sync {
     fn get_common(&self) -> &InteractionData;
     fn get_common_mut(&mut self) -> &mut InteractionData;
     fn as_any(&self) -> &dyn Any;
+    fn as_any_mut(&mut self) -> &mut dyn Any;
+    fn into_any(self: Box<Self>) -> Box<dyn Any>;
 
     fn p(&self) -> Point3f {
         self.get_common().pi.into()
@@ -111,6 +117,213 @@ pub struct SurfaceInteraction {
 }
 
 impl SurfaceInteraction {
+    pub fn compute_differentials(&mut self, r: &Ray, camera: Camera, samples_per_pixel: i32) {
+        let computed = if let Some(diff) = &r.differential {
+            let dot_rx = self.common.n.dot(diff.rx_direction.into());
+            let dot_ry = self.common.n.dot(diff.ry_direction.into());
+
+            if dot_rx != 0.0 && dot_ry != 0.0 {
+                //  Estimate screen-space change in p using ray differentials>
+                let p_as_vec = Normal3f::new(self.p().x(), self.p().y(), self.p().z());
+                let d = -self.common.n.dot(p_as_vec);
+
+                // Compute t for x-auxiliary ray
+                let rx_origin_vec =
+                    Normal3f::new(diff.rx_origin.x(), diff.rx_origin.y(), diff.rx_origin.z());
+                let tx = (-self.common.n.dot(rx_origin_vec) - d) / dot_rx;
+
+                // Compute intersection point px
+                let px = diff.rx_origin + diff.rx_direction * tx;
+
+                // Compute t for y-auxiliary ray
+                let ry_origin_vec =
+                    Normal3f::new(diff.ry_origin.x(), diff.ry_origin.y(), diff.ry_origin.z());
+                let ty = (-self.common.n.dot(ry_origin_vec) - d) / dot_ry;
+
+                let py = diff.ry_origin + diff.ry_direction * ty;
+
+                self.dpdx = px - self.p();
+                self.dpdy = py - self.p();
+
+                true
+            } else {
+                false
+            }
+        } else {
+            false
+        };
+
+        if !computed {
+            camera.approximate_dp_dxy(
+                self.p(),
+                self.n(),
+                self.time(),
+                samples_per_pixel,
+                &mut self.dpdx,
+                &mut self.dpdy,
+            );
+        }
+
+        let ata00 = self.dpdu.dot(self.dpdu);
+        let ata01 = self.dpdu.dot(self.dpdv);
+        let ata11 = self.dpdv.dot(self.dpdv);
+        let mut inv_det = 1. / difference_of_products(ata00, ata11, ata01, ata01);
+        inv_det = if inv_det.is_finite() { inv_det } else { 0. };
+        let atb0x = self.dpdu.dot(self.dpdx);
+        let atb1x = self.dpdv.dot(self.dpdx);
+        let atb0y = self.dpdu.dot(self.dpdy);
+        let atb1y = self.dpdv.dot(self.dpdy);
+        // Compute u and v derivatives in x and y
+        self.dudx = difference_of_products(ata11, atb0x, ata01, atb1x) * inv_det;
+        self.dvdx = difference_of_products(ata00, atb1x, ata01, atb0x) * inv_det;
+        self.dudy = difference_of_products(ata11, atb0y, ata01, atb1y) * inv_det;
+        self.dvdy = difference_of_products(ata00, atb1y, ata01, atb0y) * inv_det;
+        // Clamp derivatives
+        self.dudx = if self.dudx.is_finite() {
+            clamp_t(self.dudx, -1e8, 1e8)
+        } else {
+            0.
+        };
+        self.dvdx = if self.dvdx.is_finite() {
+            clamp_t(self.dvdx, -1e8, 1e8)
+        } else {
+            0.
+        };
+        self.dudy = if self.dudy.is_finite() {
+            clamp_t(self.dudy, -1e8, 1e8)
+        } else {
+            0.
+        };
+        self.dvdy = if self.dvdy.is_finite() {
+            clamp_t(self.dvdy, -1e8, 1e8)
+        } else {
+            0.
+        };
+    }
+
+    pub fn skip_intersection(&self, ray: &mut Ray, t: Float) {
+        let new_ray = Ray::spawn(&self.pi(), &self.n(), ray.time, ray.d);
+        ray.o = new_ray.o;
+        // Skipping other variables, since they should not change when passing through surface
+        if let Some(diff) = &mut ray.differential {
+            diff.rx_origin += diff.rx_direction * t;
+            diff.ry_origin += diff.ry_direction * t;
+        }
+    }
+
+    pub fn spawn_ray(&self, d: Vector3f) -> Ray {
+        let mut ray = Ray::spawn(&self.pi(), &self.n(), self.time(), d);
+        ray.medium = self.get_medium(d);
+        ray
+    }
+
+    pub fn spawn_ray_to(&self, p2: Point3f) -> Ray {
+        let mut ray = Ray::spawn_to_point(&self.pi(), &self.n(), self.time(), p2);
+        ray.medium = self.get_medium(ray.d);
+        ray
+    }
+
+    // fn get_medium(&self, w: &Vector3f) -> Option<Arc<Medium>> {
+    //     if w.dot(self.n().into()) > 0.0 {
+    //         self.medium_interface.outside.clone()
+    //     } else {
+    //         self.medium_interface.inside.clone()
+    //     }
+    // }
+
+    pub fn spawn_ray_with_differentials(
+        &self,
+        ray_i: &Ray,
+        wi: Vector3f,
+        flags: BxDFFlags,
+        eta: Float,
+    ) -> Ray {
+        let mut rd = self.spawn_ray(wi);
+
+        if let Some(diff_i) = &ray_i.differential {
+            let mut n = self.shading.n;
+
+            let mut dndx = self.shading.dndu * self.dudx + self.shading.dndv * self.dvdx;
+            let mut dndy = self.shading.dndu * self.dudy + self.shading.dndv * self.dvdy;
+
+            let dwodx = -diff_i.rx_direction - self.wo();
+            let dwody = -diff_i.ry_direction - self.wo();
+
+            let new_diff_rx_origin = self.p() + self.dpdx;
+            let new_diff_ry_origin = self.p() + self.dpdy;
+            let mut new_diff_rx_dir = Vector3f::default();
+            let mut new_diff_ry_dir = Vector3f::default();
+
+            let mut valid_differentials = false;
+
+            if flags.contains(BxDFFlags::SPECULAR_REFLECTION) {
+                valid_differentials = true;
+
+                let d_wo_dot_n_dx = dwodx.dot(n.into()) + self.wo().dot(dndx.into());
+                let d_wo_dot_n_dy = dwody.dot(n.into()) + self.wo().dot(dndy.into());
+
+                let wo_dot_n = self.wo().dot(n.into());
+
+                new_diff_rx_dir = wi - dwodx
+                    + (Vector3f::from(dndx) * wo_dot_n + Vector3f::from(n) * d_wo_dot_n_dx) * 2.0;
+
+                new_diff_ry_dir = wi - dwody
+                    + (Vector3f::from(dndy) * wo_dot_n + Vector3f::from(n) * d_wo_dot_n_dy) * 2.0;
+            } else if flags.contains(BxDFFlags::SPECULAR_TRANSMISSION) {
+                valid_differentials = true;
+
+                if self.wo().dot(n.into()) < 0.0 {
+                    n = -n;
+                    dndx = -dndx;
+                    dndy = -dndy;
+                }
+
+                // Compute partial derivatives
+                let d_wo_dot_n_dx = dwodx.dot(n.into()) + self.wo().dot(dndx.into());
+                let d_wo_dot_n_dy = dwody.dot(n.into()) + self.wo().dot(dndy.into());
+
+                let wo_dot_n = self.wo().dot(n.into());
+                let wi_dot_n = wi.dot(Vector3f::from(n));
+                let abs_wi_dot_n = wi.abs_dot(n.into());
+
+                let mu = wo_dot_n / eta - abs_wi_dot_n;
+
+                let f_eta = 1.0 / eta;
+                let f_eta2 = 1.0 / square(eta);
+
+                let term = f_eta + (f_eta2 * wo_dot_n / wi_dot_n);
+
+                let dmudx = d_wo_dot_n_dx * term;
+                let dmudy = d_wo_dot_n_dy * term;
+
+                new_diff_rx_dir =
+                    wi - dwodx * eta + (Vector3f::from(dndx) * mu + Vector3f::from(n) * dmudx);
+                new_diff_ry_dir =
+                    wi - dwody * eta + (Vector3f::from(dndy) * mu + Vector3f::from(n) * dmudy);
+            }
+
+            if valid_differentials {
+                let threshold = 1e16;
+
+                if new_diff_rx_dir.norm_squared() > threshold
+                    || new_diff_ry_dir.norm_squared() > threshold
+                    || Vector3f::from(new_diff_rx_origin).norm_squared() > threshold
+                    || Vector3f::from(new_diff_ry_origin).norm_squared() > threshold
+                {
+                    rd.differential = None;
+                } else {
+                    rd.differential = Some(RayDifferential {
+                        rx_origin: new_diff_rx_origin,
+                        ry_origin: new_diff_ry_origin,
+                        rx_direction: new_diff_rx_dir,
+                        ry_direction: new_diff_ry_dir,
+                    });
+                }
+            }
+        }
+
+        rd
+    }
 }
 
 #[derive(Default, Debug, Clone)]
@@ -123,6 +336,18 @@ pub struct MediumInteraction {
 }
 
 impl Interaction for SurfaceInteraction {
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn as_any_mut(&mut self) -> &mut dyn Any {
+        self
+    }
+
+    fn into_any(self: Box<Self>) -> Box<dyn Any> {
+        self
+    }
+
     fn get_common(&self) -> &InteractionData {
         &self.common
     }
@@ -131,10 +356,6 @@ impl Interaction for SurfaceInteraction {
         &mut self.common
     }
 
-    fn as_any(&self) -> &dyn Any {
-        self
-    }
-
     fn get_medium(&self, w: Vector3f) -> Option<Arc<Medium>> {
         self.medium_interface.as_ref().and_then(|interface| {
             if self.n().dot(w.into()) > 0.0 {
@@ -197,7 +418,14 @@ impl SurfaceInteraction {
         }
 
         Self {
-            common: InteractionData { pi, n, time, wo, medium_interface: None, medium: None },
+            common: InteractionData {
+                pi,
+                n,
+                time,
+                wo,
+                medium_interface: None,
+                medium: None,
+            },
             uv,
             dpdu,
             dpdv,
@@ -284,21 +512,30 @@ impl SurfaceInteraction {
     ) {
         self.material = Some(mtl);
         self.area_light = Some(area);
-        if prim_medium_interface.as_ref().map_or(false, |mi| mi.is_medium_transition()) {
+        if prim_medium_interface
+            .as_ref()
+            .map_or(false, |mi| mi.is_medium_transition())
+        {
             self.common.medium_interface = prim_medium_interface;
         } else {
             self.common.medium = Some(ray_medium);
         }
     }
-
 }
 
-// pub enum InteractionEnum {
-//     Surface(SurfaceInteraction),
-//     Medium(MediumInteraction),
-// }
-
 impl Interaction for MediumInteraction {
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn as_any_mut(&mut self) -> &mut dyn Any {
+        self
+    }
+
+    fn into_any(self: Box<Self>) -> Box<dyn Any> {
+        self
+    }
+
     fn get_common(&self) -> &InteractionData {
         &self.common
     }
@@ -307,10 +544,6 @@ impl Interaction for MediumInteraction {
         &mut self.common
     }
 
-    fn as_any(&self) -> &dyn Any {
-        self
-    }
-
     fn get_medium(&self, _w: Vector3f) -> Option<Arc<Medium>> {
         Some(self.medium.clone())
     }
diff --git a/src/core/medium.rs b/src/core/medium.rs
index c7bc2b2..bfc0715 100644
--- a/src/core/medium.rs
+++ b/src/core/medium.rs
@@ -32,26 +32,10 @@ impl MediumInterface {
     }
 
     pub fn is_medium_transition(&self) -> bool {
-        if let Some(ref inside) = self.inside {
-            // self.inside == Some
-            if let Some(ref outside) = self.outside {
-                // self.outside == Some
-                let pi = inside as *const _ as *const usize;
-                let po = outside as *const _ as *const usize;
-                pi != po
-            } else {
-                // self.outside == None
-                true
-            }
-        } else {
-            // self.inside == None
-            if let Some(ref _outside) = self.outside {
-                // self.outside == Some
-                true
-            } else {
-                // self.outside == None
-                false
-            }
+        match (&self.inside, &self.outside) {
+            (Some(inside), Some(outside)) => !Arc::ptr_eq(inside, outside),
+            (None, None) => false,
+            _ => true,
         }
     }
 }
diff --git a/src/core/mod.rs b/src/core/mod.rs
index 559a170..bdea598 100644
--- a/src/core/mod.rs
+++ b/src/core/mod.rs
@@ -1,3 +1,4 @@
+pub mod aggregates;
 pub mod bxdf;
 pub mod cie;
 pub mod film;
diff --git a/src/core/pbrt.rs b/src/core/pbrt.rs
index 72efc07..852ffb4 100644
--- a/src/core/pbrt.rs
+++ b/src/core/pbrt.rs
@@ -1,10 +1,92 @@
 use crate::geometry::{Lerp, Point2f, Vector2f, Vector3f};
-use num_traits::Num;
+use num_traits::{Num, PrimInt};
 use std::ops::{Add, Mul};
 use std::sync::atomic::{AtomicU64, Ordering as SyncOrdering};
 
 pub type Float = f32;
 
+#[cfg(not(feature = "use_f64"))]
+pub type FloatBits = u32;
+
+#[cfg(feature = "use_f64")]
+pub type FloatBits = u64;
+
+pub trait FloatBitOps: Copy {
+    type Bits;
+
+    fn to_bits_val(self) -> Self::Bits;
+    fn exponent_val(self) -> i32;
+    fn significand_val(self) -> Self::Bits;
+    fn sign_bit_val(self) -> Self::Bits;
+    fn from_bits_val(bits: Self::Bits) -> Self;
+}
+
+impl FloatBitOps for f32 {
+    type Bits = u32;
+
+    #[inline(always)]
+    fn to_bits_val(self) -> u32 {
+        self.to_bits()
+    }
+
+    #[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
+    }
+
+    #[inline(always)]
+    fn significand_val(self) -> u32 {
+        // Mask bottom 23 bits
+        self.to_bits() & ((1 << 23) - 1)
+    }
+
+    #[inline(always)]
+    fn sign_bit_val(self) -> u32 {
+        // Mask top bit (31)
+        self.to_bits() & 0x8000_0000
+    }
+
+    #[inline(always)]
+    fn from_bits_val(bits: u32) -> Self {
+        Self::from_bits(bits)
+    }
+}
+
+impl FloatBitOps for f64 {
+    type Bits = u64;
+
+    #[inline(always)]
+    fn to_bits_val(self) -> u64 {
+        self.to_bits()
+    }
+
+    #[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
+    }
+
+    #[inline(always)]
+    fn significand_val(self) -> u64 {
+        // Mask bottom 52 bits
+        self.to_bits() & ((1u64 << 52) - 1)
+    }
+
+    #[inline(always)]
+    fn sign_bit_val(self) -> u64 {
+        // Mask top bit (63)
+        self.to_bits() & 0x8000_0000_0000_0000
+    }
+
+    #[inline(always)]
+    fn from_bits_val(bits: u64) -> Self {
+        Self::from_bits(bits)
+    }
+}
+
 pub const MACHINE_EPSILON: Float = std::f32::EPSILON * 0.5;
 pub const SHADOW_EPSILON: Float = 0.0001;
 pub const ONE_MINUS_EPSILON: Float = 0.99999994;
@@ -75,28 +157,35 @@ pub fn evaluate_polynomial(t: Float, coeffs: &[Float]) -> Option<Float> {
 }
 
 #[inline]
-pub fn find_interval<P>(sz: usize, pred: P) -> usize
+pub fn find_interval<T, P>(sz: T, pred: P) -> T
 where
-    P: Fn(usize) -> bool,
+    T: PrimInt,
+    P: Fn(T) -> bool,
 {
-    if sz <= 2 {
-        return 0;
+    let zero = T::zero();
+    let one = T::one();
+    let two = one + one;
+
+    if sz <= two {
+        return zero;
     }
 
-    let mut low = 1;
-    let mut high = sz - 1;
+    let mut low = one;
+    let mut high = sz - one;
 
     while low < high {
-        let mid = low + (high - low) / 2;
+        // mid = low + (high - low) / 2
+        let mid = low + (high - low) / two;
         if pred(mid) {
-            low = mid + 1;
+            low = mid + one;
         } else {
             high = mid;
         }
     }
 
-    let result = low - 1;
-    clamp_t(result, 0, sz - 2)
+    let result = low - one;
+
+    num_traits::clamp(result, zero, sz - two)
 }
 
 #[inline]
diff --git a/src/core/primitive.rs b/src/core/primitive.rs
index 73f836c..e04a405 100644
--- a/src/core/primitive.rs
+++ b/src/core/primitive.rs
@@ -1,11 +1,12 @@
+use crate::core::aggregates::LinearBVHNode;
 use crate::core::interaction::{Interaction, SurfaceInteraction};
-use crate::core::pbrt::Float;
-use crate::geometry::{Bounds3f, Ray};
-use crate::shapes::{ShapeIntersection, Shape};
 use crate::core::material::MaterialTrait;
-use crate::core::medium::MediumInterface;
-use crate::lights::Light;
+use crate::core::medium::{Medium, MediumInterface};
+use crate::core::pbrt::Float;
 use crate::core::texture::{FloatTextureTrait, TextureEvalContext};
+use crate::geometry::{Bounds3f, Ray};
+use crate::lights::Light;
+use crate::shapes::{Shape, ShapeIntersection};
 use crate::utils::hash::hash_float;
 use crate::utils::transform::{AnimatedTransform, Transform};
 
@@ -24,31 +25,47 @@ pub struct GeometricPrimitive {
     area_light: Arc<Light>,
     medium_interface: Arc<MediumInterface>,
     alpha: Option<Arc<dyn FloatTextureTrait>>,
-
 }
 
 impl PrimitiveTrait for GeometricPrimitive {
     fn bounds(&self) -> Bounds3f {
         self.shape.bounds()
     }
+
     fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let mut si = self.shape.intersect(r, t_max)?;
-        let ctx: TextureEvalContext = si.intr().into();
-        let Some(ref alpha) = self.alpha else { return None };
-        let a = alpha.evaluate(&ctx);
-        if a < 1. {
-            let u = if a <= 0. { 1. } else { hash_float((r.o, r.d)) };
-            if u > a {
-                let r_next = si.intr().spawn_ray(r.d);
-                let new_t_max = t_max.map(|t| t - si.t_hit())?;
-                let mut si_next = self.intersect(&r_next, Some(new_t_max - si.t_hit()))?;
-                si_next.set_t_hit(si_next.t_hit() + si.t_hit());
-                return Some(si_next)
+        if let Some(ref alpha) = self.alpha {
+            let ctx: TextureEvalContext = (&*si.intr_mut()).into();
+            let a = alpha.evaluate(&ctx);
+            if a < 1.0 {
+                let u = if a <= 0.0 {
+                    1.0
+                } else {
+                    hash_float((r.o, r.d))
+                };
+
+                if u > a {
+                    let r_next = si.intr().spawn_ray(r.d);
+
+                    let new_t_max = t_max.map(|t| t - si.t_hit());
+
+                    if let Some(mut si_next) = self.intersect(&r_next, new_t_max) {
+                        si_next.set_t_hit(si_next.t_hit() + si.t_hit());
+                        return Some(si_next);
+                    } else {
+                        return None;
+                    }
+                }
             }
         }
-        if let Some(si) = si.intr().downcast_mut::<SurfaceInteraction>() {
-            si.as_ref().set_intersection_properties(self.material.clone(), self.area_light.clone(), Some(self.medium_interface.clone()), r.medium.clone()?);
-        }
+
+        si.intr_mut().set_intersection_properties(
+            self.material.clone(),
+            self.area_light.clone(),
+            Some(self.medium_interface.clone()),
+            r.medium.clone().expect("Medium not set"),
+        );
+
         Some(si)
     }
 
@@ -75,11 +92,25 @@ pub struct TransformedPrimitive {
 
 impl PrimitiveTrait for TransformedPrimitive {
     fn bounds(&self) -> Bounds3f {
-        self.render_from_primitive.apply_to_bounds(self.primitive.bounds())
+        self.render_from_primitive
+            .apply_to_bounds(self.primitive.bounds())
     }
 
-    fn intersect(&self, _r: &Ray, _t_max: Option<Float>) -> Option<ShapeIntersection> {
-        todo!()
+    fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+        let (ray, t_max) = self.render_from_primitive.apply_inverse_ray(r, t_max);
+        let mut si = self.primitive.intersect(&ray, Some(t_max))?;
+        let new_render: Box<dyn Interaction> = self
+            .render_from_primitive
+            .apply_to_interaction(si.intr.as_ref());
+        let new_render_surf = new_render
+            .into_any()
+            .downcast::<SurfaceInteraction>()
+            .expect(
+                "Failed to downcast Interaction to SurfaceInteraction. This should not happen.",
+            );
+
+        si.intr = new_render_surf;
+        Some(si)
     }
 
     fn intersect_p(&self, _r: &Ray, _t_max: Option<Float>) -> bool {
@@ -90,29 +121,67 @@ impl PrimitiveTrait for TransformedPrimitive {
 #[derive(Debug, Clone)]
 pub struct AnimatedPrimitive {
     primitive: Arc<dyn PrimitiveTrait>,
-    render_from_primitive: AnimatedTransform
+    render_from_primitive: AnimatedTransform,
 }
 
 impl PrimitiveTrait for AnimatedPrimitive {
     fn bounds(&self) -> Bounds3f {
-        self.render_from_primitive.motion_bounds(&self.primitive.bounds())
+        self.render_from_primitive
+            .motion_bounds(&self.primitive.bounds())
     }
 
     fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
         let interp_render_from_primitive = self.render_from_primitive.interpolate(r.time);
-        let ray = interp_render_from_primitive.apply_inverse_ray(r, t_max);
-        let mut si = self.primitive.intersect(r, t_max) else { return None };
-        let new_render = interp_render_from_primitive.apply_to_interaction(si?.intr());
-        let transform = new_render.as_any().downcast_mut::<SurfaceInteraction>().expect(
-        "Failed to downcast Interaction to SurfaceInteraction. This should not happen."
-        );
+        let (ray, t_max) = interp_render_from_primitive.apply_inverse_ray(r, t_max);
+        let mut si = self.primitive.intersect(&ray, Some(t_max))?;
+        let new_render: Box<dyn Interaction> =
+            interp_render_from_primitive.apply_to_interaction(si.intr.as_ref());
+        let new_render_surf = new_render
+            .into_any()
+            .downcast::<SurfaceInteraction>()
+            .expect(
+                "Failed to downcast Interaction to SurfaceInteraction. This should not happen.",
+            );
 
-        *si?.intr = new_render;
-        si
+        si.intr = new_render_surf;
+        Some(si)
     }
 
+    fn intersect_p(&self, _r: &Ray, _t_max: Option<Float>) -> bool {
+        todo!()
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct BVHAggregatePrimitive {
+    max_prims_in_node: usize,
+    primitives: Vec<Arc<dyn PrimitiveTrait>>,
+    nodes: Vec<LinearBVHNode>,
+}
+
+impl BVHAggregatePrimitive {
+    fn bounds(&self) -> Bounds3f {
+        if !self.nodes.is_empty() {
+            self.nodes[0].bounds
+        } else {
+            Bounds3f::default()
+        }
+    }
+
+    fn intersect(&self, r: &Ray, t_max: Option<Float>) -> Option<ShapeIntersection> {
+        if self.nodes.is_empty() {
+            return None;
+        }
+        self.intersect(r, t_max)
+    }
+
+    fn intersect_p(&self, r: &Ray, t_max: Option<Float>) -> bool {
+        if self.nodes.is_empty() {
+            return false;
+        }
+        self.intersect_p(r, t_max)
+    }
 }
-pub struct BVHAggregatePrimitive;
 pub struct KdTreeAggregate;
 
 pub enum Primitive {
@@ -122,29 +191,3 @@ pub enum Primitive {
     BVH(BVHAggregatePrimitive),
     KdTree(KdTreeAggregate),
 }
-
-
-impl Primitive {
-    // pub fn bounds(&self) -> Bounds3f {
-    //     match self {
-    //         Primitive::Geometric(primitive) => primitive.bounds(),
-    //         Primitive::Transformed(primitive) => primitive.bounds(),
-    //         Primitive::Animated(primitive) => primitive.bounds(),
-    //         Primitive::BVH(primitive) => primitive.bounds(),
-    //         Primitive::KdTree(primitive) => primitive.bounds(),
-    //     }
-    // }
-}
-
-// struct GeometricPrimitive {
-//     shape: Shape,
-//     material: Material,
-//     area_light: Light,
-//     medium_interface: MediumInterface,
-//     alpha: Texture<f64>,
-// }
-//
-//
-// impl GeometricPrimitive {
-//     fn new(shape: Shape, material: Material, medium_interface: MediumInterface, alpha: Texture<f64>)
-// }
diff --git a/src/core/sampler.rs b/src/core/sampler.rs
index 48d1145..6947734 100644
--- a/src/core/sampler.rs
+++ b/src/core/sampler.rs
@@ -1,6 +1,23 @@
-use crate::core::pbrt::{Float, PI, PI_OVER_2, PI_OVER_4, find_interval, lerp};
+use std::ops::RangeFull;
+
+use rand::seq::index::sample;
+
+use crate::core::filter::FilterTrait;
+use crate::core::pbrt::{
+    Float, ONE_MINUS_EPSILON, PI, PI_OVER_2, PI_OVER_4, clamp_t, find_interval, lerp,
+};
 use crate::geometry::{Bounds2f, Point2f, Point2i, Vector2f};
 use crate::utils::containers::Array2D;
+use crate::utils::math::{
+    BinaryPermuteScrambler, DigitPermutation, FastOwenScrambler, NoRandomizer, OwenScrambler,
+    PRIME_TABLE_SIZE, Scrambler, compute_radical_inverse_permutations, encode_morton_2,
+    inverse_radical_inverse, log2_int, owen_scrambled_radical_inverse, permutation_element,
+    radical_inverse, round_up_pow2, scrambled_radical_inverse, sobol_interval_to_index,
+    sobol_sample,
+};
+use crate::utils::rng::Rng;
+use crate::utils::sobol::N_SOBOL_DIMENSIONS;
+use crate::utils::{hash::*, sobol};
 
 #[derive(Debug, Clone, Copy)]
 pub struct CameraSample {
@@ -21,156 +38,770 @@ impl Default for CameraSample {
     }
 }
 
-#[derive(Debug, Clone)]
-pub struct PiecewiseConstant1D {
-    pub func: Vec<Float>,
-    pub cdf: Vec<Float>,
-    pub min: Float,
-    pub max: Float,
-    pub func_integral: Float,
+pub fn get_camera_sample<S, F>(sampler: &mut S, p_pixel: Point2i, filter: &F) -> CameraSample
+where
+    S: SamplerTrait,
+    F: FilterTrait,
+{
+    let fs = filter.sample(sampler.get_pixel2d());
+    let mut cs = CameraSample::default();
+    cs.p_film = Point2f::from(p_pixel) + Vector2f::from(fs.p) + Vector2f::new(0.5, 0.5);
+    cs.time = sampler.get1d();
+    cs.p_lens = sampler.get2d();
+    cs.filter_weight = fs.weight;
+    cs
 }
-impl PiecewiseConstant1D {
-    pub fn new(f: &[Float]) -> Self {
-        Self::new_with_bounds(f, 0., 1.)
+
+#[derive(Default, Debug, Clone)]
+pub struct IndependentSampler {
+    samples_per_pixel: usize,
+    seed: u64,
+    rng: Rng,
+}
+
+impl IndependentSampler {
+    pub fn new(samples_per_pixel: usize, seed: u64) -> Self {
+        Self {
+            samples_per_pixel,
+            seed,
+            rng: Rng::default(),
+        }
     }
 
-    pub fn new_with_bounds(f: &[Float], min: Float, max: Float) -> Self {
-        assert!(max > min);
-        let n = f.len();
-        let mut func = Vec::with_capacity(n);
-        for &val in f {
-            func.push(val.abs());
+    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>) {
+        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);
+        self.rng
+            .advance((sample_index as u64) * 65536 + (dim.unwrap_or(0) as u64));
+    }
+
+    pub fn get1d(&mut self) -> Float {
+        self.rng.uniform::<Float>()
+    }
+    pub fn get2d(&mut self) -> Point2f {
+        Point2f::new(self.rng.uniform::<Float>(), self.rng.uniform::<Float>())
+    }
+    pub fn get_pixel2d(&mut self) -> Point2f {
+        self.get2d()
+    }
+}
+
+const MAX_HALTON_RESOLUTION: i32 = 128;
+
+#[derive(Debug, Default, Clone, PartialEq, Eq)]
+pub enum RandomizeStrategy {
+    #[default]
+    None,
+    PermuteDigits,
+    FastOwen,
+    Owen,
+}
+
+#[derive(Default, Debug, Clone)]
+pub struct HaltonSampler {
+    samples_per_pixel: usize,
+    randomize: RandomizeStrategy,
+    digit_permutations: Vec<DigitPermutation>,
+    base_scales: [u64; 2],
+    base_exponents: [u64; 2],
+    mult_inverse: [u64; 2],
+    halton_index: u64,
+    dim: usize,
+}
+
+impl HaltonSampler {
+    pub fn new(
+        samples_per_pixel: usize,
+        full_res: Point2i,
+        randomize: RandomizeStrategy,
+        seed: u64,
+    ) -> Self {
+        let digit_permutations = compute_radical_inverse_permutations(seed);
+        let mut base_scales = [0u64; 2];
+        let mut base_exponents = [0u64; 2];
+        let bases = [2, 3];
+        let res_coords = [full_res.x(), full_res.y()];
+
+        for i in 0..2 {
+            let base = bases[i] as u64;
+            let mut scale = 1u64;
+            let mut exp = 0u64;
+
+            let limit = std::cmp::min(res_coords[i], MAX_HALTON_RESOLUTION) as u64;
+
+            while scale < limit {
+                scale *= base;
+                exp += 1;
+            }
+
+            base_scales[i] = scale;
+            base_exponents[i] = exp;
         }
 
-        let mut cdf = vec![0.; n + 1];
-        for i in 1..=n {
-            debug_assert!(func[i - 1] >= 0.);
-            cdf[i] = cdf[i - 1] + func[i - 1] * (max - min) / n as Float;
-        }
+        let mut mult_inverse = [0u64; 2];
 
-        let func_integral = cdf[n];
-        if func_integral == 0. {
-            for i in 1..=n {
-                cdf[i] = i as Float / n as Float;
-            }
-        } else {
-            for i in 1..=n {
-                cdf[i] /= func_integral;
-            }
-        }
+        mult_inverse[0] =
+            Self::multiplicative_inverse(base_scales[0] as i64, base_scales[0] as i64);
+        mult_inverse[1] =
+            Self::multiplicative_inverse(base_scales[1] as i64, base_scales[1] as i64);
 
         Self {
-            func,
-            cdf,
-            func_integral,
-            min,
-            max,
+            samples_per_pixel,
+            randomize,
+            digit_permutations,
+            base_scales,
+            base_exponents,
+            mult_inverse,
+            halton_index: 0,
+            dim: 0,
         }
     }
-
-    pub fn integral(&self) -> Float {
-        self.func_integral
+    pub fn samples_per_pixel(&self) -> usize {
+        self.samples_per_pixel
     }
 
-    pub fn size(&self) -> usize {
-        self.func.len()
-    }
+    pub fn start_pixel_sample(&mut self, p: Point2i, sample_index: usize, dim: Option<usize>) {
+        self.halton_index = 0;
 
-    pub fn sample(&self, u: Float) -> (Float, Float, usize) {
-        let o = find_interval(self.cdf.len(), |idx| self.cdf[idx] <= u);
-        let mut du = u - self.cdf[o];
-        if self.cdf[o + 1] - self.cdf[o] > 0. {
-            du /= self.cdf[o + 1] - self.cdf[o];
+        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;
         }
-        debug_assert!(!du.is_nan());
-        let value = lerp((o as Float + du) / self.size() as Float, self.min, self.max);
-        let pdf_val = if self.func_integral > 0. {
-            self.func[o] / self.func_integral
+
+        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);
+        } else if self.randomize == RandomizeStrategy::PermuteDigits {
+            return scrambled_radical_inverse(
+                dimension,
+                self.halton_index,
+                &self.digit_permutations[dimension],
+            );
         } else {
-            0.
-        };
-        (value, pdf_val, o)
+            return owen_scrambled_radical_inverse(
+                dimension,
+                self.halton_index,
+                mix_bits(1 + ((dimension as u64) << 4)) as u32,
+            );
+        }
     }
-}
-#[derive(Debug, Clone)]
-pub struct PiecewiseConstant2D {
-    pub p_conditional_v: Vec<PiecewiseConstant1D>,
-    pub p_marginal: PiecewiseConstant1D,
-    pub domain: Bounds2f,
-}
 
-impl PiecewiseConstant2D {
-    pub fn new(data: &Array2D<Float>, nu: usize, nv: usize, domain: Bounds2f) -> Self {
-        let mut p_conditional_v = Vec::with_capacity(nv);
-        for v in 0..nv {
-            let start = v * nu;
-            let end = start + nu;
-            p_conditional_v.push(PiecewiseConstant1D::new_with_bounds(
-                &data.as_slice()[start..end],
-                domain.p_min.x(),
-                domain.p_max.x(),
-            ));
+    fn multiplicative_inverse(a: i64, n: i64) -> u64 {
+        let (x, _) = Self::extended_gcd(a as u64, n as u64);
+        x.rem_euclid(n) as u64
+    }
+
+    fn extended_gcd(a: u64, b: u64) -> (i64, i64) {
+        if b == 0 {
+            return (1, 0);
         }
 
-        let marginal_func: Vec<Float> = p_conditional_v.iter().map(|p| p.integral()).collect();
-        let p_marginal = PiecewiseConstant1D::new_with_bounds(
-            &marginal_func,
-            domain.p_min.y(),
-            domain.p_max.y(),
+        let (xp, yp) = Self::extended_gcd(b, a % b);
+
+        let d = (a / b) as i64;
+
+        (yp, xp - d * yp)
+    }
+}
+
+#[derive(Default, Debug, Clone)]
+pub struct StratifiedSampler {
+    x_pixel_samples: usize,
+    y_pixel_samples: usize,
+    jitter: bool,
+    seed: u64,
+    rng: Rng,
+    pixel: Point2i,
+    sample_index: usize,
+    dim: usize,
+}
+
+impl StratifiedSampler {
+    pub fn new(
+        x_pixel_samples: usize,
+        y_pixel_samples: usize,
+        seed: Option<u64>,
+        jitter: bool,
+    ) -> Self {
+        Self {
+            x_pixel_samples,
+            y_pixel_samples,
+            jitter,
+            seed: seed.unwrap_or(0),
+            rng: Rng::default(),
+            pixel: Point2i::default(),
+            sample_index: 0,
+            dim: 0,
+        }
+    }
+
+    pub 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>) {
+        self.pixel = p;
+        self.sample_index = sample_index;
+        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);
+        self.rng
+            .advance((sample_index as u64) * 65536 + (dim.unwrap_or(0) as u64));
+    }
+
+    pub fn get1d(&mut self) -> Float {
+        let hash_input = [
+            self.pixel.x() as u64,
+            self.pixel.y() as u64,
+            self.dim as u64,
+            self.seed,
+        ];
+        let hash = hash_buffer(&hash_input, 0);
+        let stratum = permutation_element(
+            self.sample_index as u32,
+            self.samples_per_pixel() as u32,
+            hash as u32,
         );
+        self.dim += 1;
 
+        let delta = if self.jitter {
+            self.rng.uniform::<Float>()
+        } else {
+            0.5
+        };
+        (stratum as Float + delta) / (self.samples_per_pixel() as Float)
+    }
+
+    pub fn get2d(&mut self) -> Point2f {
+        let hash_input = [
+            self.pixel.x() as u64,
+            self.pixel.y() as u64,
+            self.dim as u64,
+            self.seed,
+        ];
+        let hash = hash_buffer(&hash_input, 0);
+        let stratum = permutation_element(
+            self.sample_index as u32,
+            self.samples_per_pixel() as u32,
+            hash as u32,
+        );
+        self.dim += 2;
+        let x = stratum % self.x_pixel_samples as u32;
+        let y = stratum / self.y_pixel_samples as u32;
+        let dx = if self.jitter {
+            self.rng.uniform::<Float>()
+        } else {
+            0.5
+        };
+        let dy = if self.jitter {
+            self.rng.uniform::<Float>()
+        } else {
+            0.5
+        };
+        Point2f::new(
+            (x as Float + dx) / self.x_pixel_samples as Float,
+            (y as Float + dy) / self.y_pixel_samples as Float,
+        )
+    }
+
+    pub fn get_pixel2d(&mut self) -> Point2f {
+        self.get2d()
+    }
+}
+
+#[derive(Default, Debug, Clone)]
+pub struct PaddedSobolSampler {
+    samples_per_pixel: usize,
+    seed: u64,
+    randomize: RandomizeStrategy,
+    pixel: Point2i,
+    sample_index: usize,
+    dim: usize,
+}
+
+impl PaddedSobolSampler {
+    pub fn new(samples_per_pixel: usize, randomize: RandomizeStrategy, seed: Option<u64>) -> Self {
         Self {
-            p_conditional_v,
-            p_marginal,
-            domain,
+            samples_per_pixel,
+            seed: seed.unwrap_or(0),
+            randomize,
+            pixel: Point2i::default(),
+            sample_index: 0,
+            dim: 0,
         }
     }
-
-    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)
+    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.sample_index = sample_index;
+        self.dim = dim.unwrap_or(0);
     }
 
-    pub fn new_with_data(data: &Array2D<Float>, nx: usize, ny: usize) -> Self {
-        Self::new(
-            data,
-            nx,
-            ny,
-            Bounds2f::from_points(Point2f::new(0., 0.), Point2f::new(1., 1.)),
+    pub fn get1d(&mut self) -> Float {
+        let hash_input = [
+            self.pixel.x() as u64,
+            self.pixel.y() as u64,
+            self.dim as u64,
+            self.seed,
+        ];
+        let hash = hash_buffer(&hash_input, 0) as u32;
+        let index = permutation_element(
+            self.sample_index as u32,
+            self.samples_per_pixel as u32,
+            hash,
+        );
+        self.sample_dimension(0, index, hash >> 32)
+    }
+    pub fn get2d(&mut self) -> Point2f {
+        let hash_input = [
+            self.pixel.x() as u64,
+            self.pixel.y() as u64,
+            self.dim as u64,
+            self.seed,
+        ];
+        let hash = hash_buffer(&hash_input, 0) as u32;
+        let index = permutation_element(
+            self.sample_index as u32,
+            self.samples_per_pixel as u32,
+            hash,
+        );
+        self.dim += 2;
+        Point2f::new(
+            self.sample_dimension(0, index, hash),
+            self.sample_dimension(1, index, hash >> 32),
         )
     }
 
-    pub fn resolution(&self) -> Point2i {
-        Point2i::new(
-            self.p_conditional_v[0].size() as i32,
-            self.p_conditional_v[1].size() as i32,
-        )
+    pub fn get_pixel2d(&mut self) -> Point2f {
+        self.get2d()
     }
 
-    pub fn integral(&self) -> f32 {
-        self.p_marginal.integral()
-    }
-
-    pub fn sample(&self, u: Point2f) -> (Point2f, f32, Point2i) {
-        let (d1, pdf1, off_y) = self.p_marginal.sample(u.y());
-        let (d0, pdf0, off_x) = self.p_conditional_v[off_y].sample(u.x());
-        let pdf = pdf0 * pdf1;
-        let offset = Point2i::new(off_x as i32, off_y as i32);
-        (Point2f::new(d0, d1), pdf, offset)
-    }
-
-    pub fn pdf(&self, p: Point2f) -> f32 {
-        let p_offset = self.domain.offset(&p);
-        let nu = self.p_conditional_v[0].size();
-        let nv = self.p_marginal.size();
-
-        let iu = (p_offset.x() * nu as f32).clamp(0.0, nu as f32 - 1.0) as usize;
-        let iv = (p_offset.y() * nv as f32).clamp(0.0, nv as f32 - 1.0) as usize;
-
-        let integral = self.p_marginal.integral();
-        if integral == 0.0 {
-            0.0
-        } else {
-            self.p_conditional_v[iv].func[iu] / integral
+    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)]
+pub struct SobolSampler {
+    samples_per_pixel: usize,
+    scale: i32,
+    seed: u64,
+    randomize: RandomizeStrategy,
+    pixel: Point2i,
+    dim: usize,
+    sobol_index: u64,
+}
+
+impl SobolSampler {
+    pub fn new(
+        samples_per_pixel: usize,
+        full_resolution: Point2i,
+        randomize: RandomizeStrategy,
+        seed: Option<u64>,
+    ) -> Self {
+        let scale = round_up_pow2(full_resolution.x().max(full_resolution.y()));
+        Self {
+            samples_per_pixel,
+            scale,
+            seed: seed.unwrap_or(0),
+            randomize,
+            pixel: Point2i::default(),
+            dim: 0,
+            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 {
+            return sobol_sample(self.sobol_index, dimension, NoRandomizer);
+        }
+        let hash_input = [self.pixel.x() as u64, self.pixel.y() as u64, self.seed];
+        let hash = hash_buffer(&hash_input, 0) as u32;
+        match self.randomize {
+            RandomizeStrategy::PermuteDigits => sobol_sample(
+                self.sobol_index,
+                dimension,
+                BinaryPermuteScrambler::new(hash),
+            ),
+            RandomizeStrategy::FastOwen => {
+                sobol_sample(self.sobol_index, dimension, FastOwenScrambler::new(hash))
+            }
+            RandomizeStrategy::Owen => {
+                sobol_sample(self.sobol_index, dimension, OwenScrambler::new(hash))
+            }
+            RandomizeStrategy::None => unreachable!(),
+        }
+    }
+}
+
+#[derive(Default, Debug, Clone)]
+pub struct ZSobolSampler {
+    randomize: RandomizeStrategy,
+    seed: u64,
+    log2_samples_per_pixel: u32,
+    n_base4_digits: u32,
+    morton_index: u64,
+    dim: usize,
+}
+
+impl ZSobolSampler {
+    pub fn new(
+        samples_per_pixel: u32,
+        full_resolution: Point2i,
+        randomize: RandomizeStrategy,
+        seed: Option<u64>,
+    ) -> 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 n_base4_digits = log2_int(res as Float) as u32 + log4_samples_per_pixel;
+        Self {
+            randomize,
+            seed: seed.unwrap_or(0),
+            log2_samples_per_pixel,
+            n_base4_digits,
+            morton_index: 0,
+            dim: 0,
+        }
+    }
+
+    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],
+            [0, 1, 3, 2],
+            [0, 2, 1, 3],
+            [0, 2, 3, 1],
+            [0, 3, 2, 1],
+            [0, 3, 1, 2],
+            [1, 0, 2, 3],
+            [1, 0, 3, 2],
+            [1, 2, 0, 3],
+            [1, 2, 3, 0],
+            [1, 3, 2, 0],
+            [1, 3, 0, 2],
+            [2, 1, 0, 3],
+            [2, 1, 3, 0],
+            [2, 0, 1, 3],
+            [2, 0, 3, 1],
+            [2, 3, 0, 1],
+            [2, 3, 1, 0],
+            [3, 1, 2, 0],
+            [3, 1, 0, 2],
+            [3, 2, 1, 0],
+            [3, 2, 0, 1],
+            [3, 0, 2, 1],
+            [3, 0, 1, 2],
+        ];
+
+        let mut sample_index = 0;
+        let pow2_samples = (self.log2_samples_per_pixel & 1) != 0;
+        let last_digit = if pow2_samples { 1 } else { 0 };
+
+        for i in (last_digit..self.n_base4_digits).rev() {
+            let digit_shift = (2 * i) - if pow2_samples { 1 } else { 0 };
+
+            let mut digit = (self.morton_index >> digit_shift) & 3;
+
+            let higher_digits = self.morton_index >> (digit_shift + 2);
+
+            let mix_input = higher_digits ^ (0x55555555 * self.dim as u64);
+            let p = (mix_bits(mix_input) >> 24) % 24;
+
+            digit = PERMUTATIONS[p as usize][digit as usize] as u64;
+
+            sample_index |= digit << digit_shift;
+        }
+
+        if pow2_samples {
+            let lsb = self.morton_index & 1;
+            sample_index |= lsb;
+        }
+
+        sample_index
+    }
+}
+#[derive(Default, Debug, Clone)]
+pub struct MLTSampler;
+impl MLTSampler {
+    pub fn samples_per_pixel(&self) -> usize {
+        todo!()
+    }
+    pub fn start_pixel_sample(&mut self, _p: Point2i, _sample_index: usize, _dim: Option<usize>) {
+        todo!()
+    }
+    pub fn get1d(&mut self) -> Float {
+        todo!()
+    }
+    pub fn get2d(&mut self) -> Point2f {
+        todo!()
+    }
+    pub fn get_pixel2d(&mut self) -> Point2f {
+        todo!()
+    }
+}
+
+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>;
+}
+
+#[derive(Debug, Clone)]
+pub enum Sampler {
+    Independent(IndependentSampler),
+    Stratified(StratifiedSampler),
+    Halton(HaltonSampler),
+    PaddedSobol(PaddedSobolSampler),
+    Sobol(SobolSampler),
+    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())
+    }
+}
diff --git a/src/geometry/bounds.rs b/src/geometry/bounds.rs
index 58fb452..009f9b6 100644
--- a/src/geometry/bounds.rs
+++ b/src/geometry/bounds.rs
@@ -77,6 +77,11 @@ where
         self.p_max - self.p_min
     }
 
+    pub fn centroid(&self) -> Point<T, N> {
+        let two = T::one() + T::one();
+        self.p_min + (self.diagonal() / two)
+    }
+
     pub fn volume(&self) -> T {
         let d = self.diagonal();
         d.0.iter().fold(T::one(), |acc, &val| acc * val)
@@ -248,6 +253,58 @@ where
 }
 
 impl Bounds3f {
+    #[inline(always)]
+    pub fn intersect_p(
+        &self,
+        o: Point3f,
+        ray_t_max: Float,
+        inv_dir: Vector3f,
+        dir_is_neg: &[usize; 3],
+    ) -> bool {
+        // We create a temporary array of references to access p_min/p_max by index
+        // The compiler optimizes this away entirely.
+        let bounds = [&self.p_min, &self.p_max];
+
+        // Check X
+        let mut t_min = (bounds[dir_is_neg[0]].x() - o.x()) * inv_dir.x();
+        let mut t_max = (bounds[1 - dir_is_neg[0]].x() - o.x()) * inv_dir.x();
+
+        // Check Y
+        let ty_min = (bounds[dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+        let ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y();
+
+        // Optional: PBRT Gamma correction for float precision (t_max *= 1.00000024)
+        // t_max *= 1.00000024;
+
+        if t_min > ty_max || ty_min > t_max {
+            return false;
+        }
+        if ty_min > t_min {
+            t_min = ty_min;
+        }
+        if ty_max < t_max {
+            t_max = ty_max;
+        }
+
+        // Check Z
+        let tz_min = (bounds[dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+        let tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z();
+
+        // t_max *= 1.00000024;
+
+        if t_min > tz_max || tz_min > t_max {
+            return false;
+        }
+        if tz_min > t_min {
+            t_min = tz_min;
+        }
+        if tz_max < t_max {
+            t_max = tz_max;
+        }
+
+        (t_min < ray_t_max) && (t_max > 0.0)
+    }
+
     pub fn intersect_with_inverse(&self, o: Point3f, d: Vector3f, ray_t_max: Float) -> bool {
         let inv_dir = Vector3::new(1.0 / d.x(), 1.0 / d.y(), 1.0 / d.z());
         let dir_is_neg: [usize; 3] = [
diff --git a/src/geometry/primitives.rs b/src/geometry/primitives.rs
index 350c016..a88b1cf 100644
--- a/src/geometry/primitives.rs
+++ b/src/geometry/primitives.rs
@@ -67,7 +67,6 @@ macro_rules! impl_tuple_core {
             }
         }
 
-
         impl<const N: usize> $Struct<f32, N> {
             #[inline]
             pub fn floor(&self) -> $Struct<i32, N> {
@@ -241,7 +240,6 @@ macro_rules! impl_float_vector_ops {
     };
 }
 
-
 macro_rules! impl_abs {
     ($Struct:ident) => {
         impl<T, const N: usize> $Struct<T, N>
@@ -538,7 +536,7 @@ where
 impl<const N: usize> Hash for Vector<Float, N> {
     fn hash<H: Hasher>(&self, state: &mut H) {
         for item in self.0.iter() {
-        item.to_bits().hash(state);
+            item.to_bits().hash(state);
         }
     }
 }
@@ -546,7 +544,7 @@ impl<const N: usize> Hash for Vector<Float, N> {
 impl<const N: usize> Hash for Point<Float, N> {
     fn hash<H: Hasher>(&self, state: &mut H) {
         for item in self.0.iter() {
-        item.to_bits().hash(state);
+            item.to_bits().hash(state);
         }
     }
 }
@@ -554,7 +552,7 @@ impl<const N: usize> Hash for Point<Float, N> {
 impl<const N: usize> Hash for Normal<Float, N> {
     fn hash<H: Hasher>(&self, state: &mut H) {
         for item in self.0.iter() {
-        item.to_bits().hash(state);
+            item.to_bits().hash(state);
         }
     }
 }
diff --git a/src/lib.rs b/src/lib.rs
index d13d191..59c5f0a 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,7 +1,6 @@
 #![allow(unused_imports, dead_code)]
 #![feature(float_erf)]
 #![feature(f16)]
-#![feature(generic_const_exprs)]
 
 mod camera;
 mod core;
diff --git a/src/shapes/bilinear.rs b/src/shapes/bilinear.rs
index 9ad0566..048bfd3 100644
--- a/src/shapes/bilinear.rs
+++ b/src/shapes/bilinear.rs
@@ -90,7 +90,7 @@ impl BilinearPatchShape {
         &meshes[self.mesh_index]
     }
 
-    fn get_data(&self) -> PatchData {
+    fn get_data(&self) -> PatchData<'_> {
         let mesh = self.mesh();
         let start_index = 4 * self.blp_index;
         let v = &mesh.vertex_indices[start_index..start_index + 4];
@@ -628,7 +628,7 @@ impl BilinearPatchShape {
         }
     }
 
-    pub fn pdf(&self, intr: &dyn Interaction) -> Float {
+    pub fn pdf(&self, intr: Arc<&dyn Interaction>) -> Float {
         let Some(si) = intr.as_any().downcast_ref::<SurfaceInteraction>() else {
             return 0.;
         };
@@ -675,7 +675,7 @@ 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(isect.intr.as_ref());
+            let isect_pdf = self.pdf(Arc::new(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. {
diff --git a/src/shapes/mod.rs b/src/shapes/mod.rs
index 6ab6847..f93383e 100644
--- a/src/shapes/mod.rs
+++ b/src/shapes/mod.rs
@@ -6,11 +6,14 @@ pub mod sphere;
 pub mod triangle;
 
 use crate::core::interaction::{Interaction, MediumInteraction, SurfaceInteraction};
+use crate::core::material::MaterialTrait;
+use crate::core::medium::{Medium, MediumInterface};
 use crate::core::pbrt::{Float, PI};
 use crate::geometry::{
     Bounds3f, DirectionCone, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector2f, Vector3f,
-    Vector3fi, VectorLike
+    Vector3fi, VectorLike,
 };
+use crate::lights::Light;
 use crate::utils::math::{next_float_down, next_float_up};
 use crate::utils::transform::Transform;
 use std::sync::{Arc, Mutex};
@@ -148,10 +151,15 @@ pub struct ShapeIntersection {
 
 impl ShapeIntersection {
     pub fn new(intr: SurfaceInteraction, t_hit: Float) -> Self {
-        Self { intr: Box::new(intr), t_hit }
+        Self {
+            intr: Box::new(intr),
+            t_hit,
+        }
     }
 
-    pub fn as_ref(&self) -> &dyn Interaction { &*self.intr }
+    pub fn as_ref(&self) -> &dyn Interaction {
+        &*self.intr
+    }
 
     pub fn intr(&self) -> &dyn Interaction {
         &*self.intr
@@ -168,6 +176,18 @@ impl ShapeIntersection {
     pub fn set_t_hit(&mut self, new_t: Float) {
         self.t_hit = new_t;
     }
+
+    pub fn set_intersection_properties(
+        &mut self,
+        mtl: Arc<dyn MaterialTrait>,
+        area: Arc<Light>,
+        prim_medium_interface: Option<Arc<MediumInterface>>,
+        ray_medium: Option<Arc<Medium>>,
+    ) {
+        let ray_medium = ray_medium.expect("Ray medium must be defined for intersection");
+        self.intr
+            .set_intersection_properties(mtl, area, prim_medium_interface, ray_medium);
+    }
 }
 
 #[derive(Debug, Clone)]
@@ -280,7 +300,7 @@ pub enum Shape {
     Curve(CurveShape),
 }
 
-impl<'a> Shape {
+impl Shape {
     pub fn bounds(&self) -> Bounds3f {
         match self {
             Shape::None => Bounds3f::default(),
diff --git a/src/utils/color.rs b/src/utils/color.rs
index 22571af..e4bd779 100644
--- a/src/utils/color.rs
+++ b/src/utils/color.rs
@@ -280,6 +280,10 @@ impl RGB {
     pub fn average(&self) -> Float {
         (self.r + self.g + self.b) / 3.0
     }
+
+    pub fn max(&self) -> Float {
+        self.r.max(self.g).max(self.b)
+    }
 }
 
 impl Index<usize> for RGB {
@@ -458,6 +462,14 @@ pub struct RGBToSpectrumTable {
     coeffs: &'static CoefficientArray,
 }
 
+impl RGBToSpectrumTable {
+    pub fn srgb() -> Self {
+        // use crate::core::constants::{RGB_TO_SPECTRUM_Z_NODES, RGB_TO_SPECTRUM_COEFFS};
+        // Self::new(&RGB_TO_SPECTRUM_Z_NODES, &RGB_TO_SPECTRUM_COEFFS)
+        todo!("Link the static constant arrays for sRGB coefficients here")
+    }
+}
+
 #[derive(Debug, Default, Copy, Clone)]
 pub struct RGBSigmoidPolynomial {
     c0: Float,
diff --git a/src/utils/colorspace.rs b/src/utils/colorspace.rs
index 8e42bc0..b045b4b 100644
--- a/src/utils/colorspace.rs
+++ b/src/utils/colorspace.rs
@@ -4,6 +4,8 @@ use super::spectrum::{DenselySampledSpectrum, SampledSpectrum, Spectrum};
 use crate::core::pbrt::Float;
 use crate::geometry::Point2f;
 
+use once_cell::sync::Lazy;
+
 use std::cmp::{Eq, PartialEq};
 use std::error::Error;
 use std::sync::Arc;
@@ -82,6 +84,22 @@ impl RGBColorspace {
 
         self.rgb_from_xyz * other.xyz_from_rgb
     }
+
+    pub fn srgb() -> &'static Self {
+        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);
+
+            let illuminant = Spectrum::std_illuminant_d65();
+            let table = RGBToSpectrumTable::srgb();
+
+            RGBColorspace::new(r, g, b, illuminant, table)
+                .expect("Failed to initialize standard sRGB color space")
+        });
+
+        &SRGB_SPACE
+    }
 }
 
 impl PartialEq for RGBColorspace {
diff --git a/src/utils/containers.rs b/src/utils/containers.rs
index d204d25..1f37d39 100644
--- a/src/utils/containers.rs
+++ b/src/utils/containers.rs
@@ -3,7 +3,9 @@ use std::hash::{BuildHasher, Hash, Hasher};
 use std::ops::{Index, IndexMut};
 use std::sync::RwLock;
 
-use crate::geometry::{Bounds2i, Bounds3i, Point2i, Point3f, Point3i, Vector3f, Vector3i};
+use crate::geometry::{
+    Bounds2i, Bounds3i, Point2i, Point3f, Point3i, Vector2i, Vector3f, Vector3i,
+};
 
 #[derive(Debug)]
 pub struct Array2D<T> {
@@ -22,6 +24,16 @@ impl<T> Array2D<T> {
         Self { values, extent }
     }
 
+    pub fn new_with_dims(nx: i32, ny: i32) -> Self
+    where
+        T: Default,
+    {
+        Self::new(Bounds2i::from_points(
+            Point2i::new(0, 0),
+            Point2i::new(nx, ny),
+        ))
+    }
+
     pub fn new_from_bounds(extent: Bounds2i, default_val: T) -> Self
     where
         T: Clone,
@@ -62,15 +74,33 @@ impl<T> Array2D<T> {
 
 impl<T> Index<Point2i> for Array2D<T> {
     type Output = T;
-    fn index(&self, p: Point2i) -> &Self::Output {
-        let idx = self.get_index(p);
+
+    fn index(&self, mut p: Point2i) -> &Self::Output {
+        p = 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]
     }
 }
 
 impl<T> IndexMut<Point2i> for Array2D<T> {
-    fn index_mut(&mut self, p: Point2i) -> &mut Self::Output {
-        let idx = self.get_index(p);
+    fn index_mut(&mut self, mut p: Point2i) -> &mut Self::Output {
+        p = 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]
     }
 }
+
+impl<T> Index<(i32, i32)> for Array2D<T> {
+    type Output = T;
+    fn index(&self, index: (i32, i32)) -> &Self::Output {
+        self.index(Point2i::new(index.0, index.1))
+    }
+}
+
+impl<T> IndexMut<(i32, i32)> for Array2D<T> {
+    fn index_mut(&mut self, index: (i32, i32)) -> &mut Self::Output {
+        self.index_mut(Point2i::new(index.0, index.1))
+    }
+}
diff --git a/src/utils/hash.rs b/src/utils/hash.rs
index aec39a4..3c2b7f2 100644
--- a/src/utils/hash.rs
+++ b/src/utils/hash.rs
@@ -1,14 +1,102 @@
-use std::hash::{Hash, Hasher};
 use crate::core::pbrt::Float;
 use std::collections::hash_map::DefaultHasher;
+use std::hash::{Hash, Hasher};
 
 const U32_TO_F32_SCALE: f32 = 1.0 / 4294967296.0;
 
-pub fn hash_float<T: Hash>(args: T) -> Float {
-    let mut hasher = DefaultHasher::new();
-    args.hash(&mut hasher);
-    let hash_u64 = hasher.finish();
-    let hash_u32 = hash_u64 as u32;
-    (hash_u32 as f32) * U32_TO_F32_SCALE
+#[inline(always)]
+pub fn mix_bits(mut v: u64) -> u64 {
+    v ^= v >> 31;
+    v = v.wrapping_mul(0x7fb5d329728ea185);
+    v ^= v >> 27;
+    v = v.wrapping_mul(0x81dadef4bc2dd44d);
+    v ^= v >> 33;
+    v
 }
 
+pub fn murmur_hash_64a(key: &[u8], seed: u64) -> u64 {
+    const M: u64 = 0xc6a4a7935bd1e995;
+    const R: i32 = 47;
+
+    let len = key.len();
+    let mut h = seed ^ ((len as u64).wrapping_mul(M));
+
+    // We chunk the slice into 8-byte segments
+    let chunks = key.chunks_exact(8);
+    let remainder = chunks.remainder();
+
+    for chunk in chunks {
+        // Safe conversion from [u8; 8] to u64
+        let mut k = u64::from_ne_bytes(chunk.try_into().unwrap());
+
+        k = k.wrapping_mul(M);
+        k ^= k >> R;
+        k = k.wrapping_mul(M);
+
+        h ^= k;
+        h = h.wrapping_mul(M);
+    }
+
+    // Handle the tail (remaining bytes)
+    if !remainder.is_empty() {
+        // We handle the switch-case fallthrough logic by building a u64
+        let mut k_tail = 0u64;
+        // Load bytes into the u64 based on length
+        for (i, &byte) in remainder.iter().enumerate() {
+            k_tail ^= (byte as u64) << (i * 8);
+        }
+
+        h ^= k_tail;
+        h = h.wrapping_mul(M);
+    }
+
+    h ^= h >> R;
+    h = h.wrapping_mul(M);
+    h ^= h >> R;
+
+    h
+}
+
+pub fn hash_buffer<T: ?Sized>(data: &T, seed: u64) -> u64 {
+    // In Rust, we need to treat the type as a byte slice.
+    // The safest way to do this generically without 'bytemuck' or 'unsafe' in signature
+    // is to require T to be castable, but usually in PBRT we just use unsafe for the hasher.
+    let len = std::mem::size_of_val(data);
+    let ptr = data as *const T as *const u8;
+    let bytes = unsafe { std::slice::from_raw_parts(ptr, len) };
+    murmur_hash_64a(bytes, seed)
+}
+
+#[macro_export]
+macro_rules! hash_values {
+    ( $( $x:expr ),* ) => {
+        {
+            // Create a packed buffer on stack matching C++ behavior
+            // We use a temporary tuple or array to pack bits
+            #[repr(C, packed)]
+            struct PackedData {
+                $(_field: Box<[u8]>),* // Phantom, logic below is simpler
+            }
+
+            // Calculate total size and create a byte buffer
+            let mut buffer = Vec::new(); // Or use a small stack array if size known
+            $(
+                let s = std::slice::from_raw_parts(
+                    &$x as *const _ as *const u8,
+                    std::mem::size_of_val(&$x)
+                );
+                buffer.extend_from_slice(s);
+            )*
+
+            $crate::utils::hash::murmur_hash_64a(&buffer, 0)
+        }
+    }
+}
+
+pub fn hash_float<T>(data: T) -> Float
+where
+    T: Copy, // Ensure it's plain old data
+{
+    let h = hash_buffer(&data, 0);
+    (h as u32) as Float * 2.3283064365386963e-10 // 0x1p-32f
+}
diff --git a/src/utils/image.rs b/src/utils/image.rs
index 5d15f0e..3dd4c78 100644
--- a/src/utils/image.rs
+++ b/src/utils/image.rs
@@ -1460,13 +1460,13 @@ impl Image {
     fn write_exr(&self, name: &Path, metadata: &ImageMetadata) -> Result<(), ImageError> {
         if self.format.is_8bit() {
             let conv_img = self
-                .convert_to_format(PixelFormat::F16, self.encoding)
+                .convert_to_format(PixelFormat::F16, self.encoding.clone())
                 .map_err(|_| {
                     ImageError::Unsupported(
                         "Failed to convert 8-bit image to F16 for EXR output.".to_string(),
                     )
                 })?;
-            conv_img.write_exr(name, metadata);
+            let _ = conv_img.write_exr(name, metadata);
         }
 
         // return Err(ImageError::Unsupported(
@@ -1536,7 +1536,7 @@ impl Image {
         }
 
         if let Some(cs) = metadata.get_color_space() {
-            if *cs != *RGBColorspace {
+            if cs != RGBColorspace::srgb() {
                 image_attributes.chromaticities = Some(ExrChromaticities {
                     red: Vec2(cs.r.x(), cs.r.y()),
                     green: Vec2(cs.g.x(), cs.g.y()),
diff --git a/src/utils/math.rs b/src/utils/math.rs
index 9caf422..fa29aac 100644
--- a/src/utils/math.rs
+++ b/src/utils/math.rs
@@ -1,13 +1,19 @@
 use super::color::{RGB, XYZ};
 use super::error::{InversionError, LlsError};
-use crate::core::pbrt::{Float, ONE_MINUS_EPSILON, PI, PI_OVER_4, clamp_t, lerp};
-use crate::geometry::{Point, Point2f, Vector, Vector3f};
+use crate::core::pbrt::{
+    Float, FloatBitOps, FloatBits, ONE_MINUS_EPSILON, PI, PI_OVER_4, clamp_t, evaluate_polynomial,
+    lerp,
+};
+use crate::geometry::{Point, Point2f, Point2i, Vector, Vector3f};
+use crate::utils::hash::{hash_buffer, mix_bits};
+use crate::utils::sobol::{SOBOL_MATRICES_32, VDC_SOBOL_MATRICES, VDC_SOBOL_MATRICES_INV};
 
 use num_traits::{Float as NumFloat, Num, One, Signed, Zero};
+use rayon::prelude::*;
 use std::error::Error;
 use std::fmt::{self, Display};
 use std::mem;
-use std::ops::{Add, Div, Index, IndexMut, Mul, Neg};
+use std::ops::{Add, Div, Index, IndexMut, Mul, Neg, Rem};
 
 #[inline]
 pub fn degrees(a: Float) -> Float {
@@ -56,6 +62,7 @@ where
     let error = fma(c, d, -cd);
     return sum_of_products + error;
 }
+
 #[inline]
 pub fn safe_sqrt(x: Float) -> Float {
     assert!(x > -1e-3);
@@ -103,6 +110,100 @@ pub fn windowed_sinc(x: Float, radius: Float, tau: Float) -> Float {
     return sinc(x) * sinc(x / tau);
 }
 
+#[inline]
+pub fn fast_exp(x: Float) -> Float {
+    let xp = x * 1.442695041;
+    let fxp = xp.floor();
+    let f = xp - fxp;
+    let i = fxp as i32;
+    let two_to_f = evaluate_polynomial(f, &[1., 0.695556856, 0.226173572, 0.0781455737])
+        .expect("Could not evaluate polynomial");
+    let exponent = exponent(two_to_f) + i;
+    if exponent < -126 {
+        return 0.;
+    }
+    if exponent > 127 {
+        return Float::INFINITY;
+    }
+    let mut bits = float_to_bits(two_to_f);
+    bits &= 0b10000000011111111111111111111111;
+    bits |= ((exponent + 127) as u32) << 23;
+    bits_to_float(bits)
+}
+
+#[inline]
+pub fn i0(x: Float) -> Float {
+    let mut val: Float = 0.0;
+    let mut x2i: Float = 1.0;
+    let mut ifact: i64 = 1;
+    let mut i4: i32 = 1;
+
+    for i in 0..10 {
+        if i > 1 {
+            ifact *= i as i64;
+        }
+
+        let denominator = (i4 as Float) * (square(ifact) as Float);
+        val += x2i / denominator;
+
+        x2i *= x * x;
+        i4 *= 4;
+    }
+
+    val
+}
+
+#[inline]
+pub fn logistic(x: Float, s: Float) -> Float {
+    let y = x.abs();
+    (-y / s).exp() / (s * square(1. + (-y / s).exp()))
+}
+
+#[inline]
+pub fn logistic_cdf(x: Float, s: Float) -> Float {
+    1. / (1. + (-x / s).exp())
+}
+
+#[inline]
+pub fn trimmed_logistic(x: Float, s: Float, a: Float, b: Float) -> Float {
+    logistic(x, s) / (logistic_cdf(b, s) - logistic_cdf(a, s))
+}
+
+#[inline]
+pub fn log_i0(x: Float) -> Float {
+    if x > 12.0 {
+        let inv_x = 1.0 / x;
+        let two_pi = 2.0 * PI as Float;
+
+        x + 0.5 * (-two_pi.ln() + inv_x.ln() + 1.0 / (8.0 * x))
+    } else {
+        i0(x).ln()
+    }
+}
+
+#[inline]
+pub fn log2_int(v: f32) -> i32 {
+    debug_assert!(v > 0.0, "log2_int requires positive input");
+
+    if v < 1.0 {
+        return -log2_int(1.0 / v);
+    }
+
+    // https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog
+    // midsignif = Significand(pow(2.0, 1.5))
+    const MID_SIGNIF: u32 = 0b00000000001101010000010011110011;
+
+    let bits = v.to_bits();
+    let exponent = ((bits >> 23) & 0xFF) as i32 - 127;
+    let significand = bits & 0x7FFFFF;
+
+    if significand >= MID_SIGNIF {
+        exponent + 1
+    } else {
+        exponent
+    }
+}
+
 pub fn quadratic(a: Float, b: Float, c: Float) -> Option<(Float, Float)> {
     if a == 0. {
         if b == 0. {
@@ -209,19 +310,29 @@ pub fn sample_linear(u: Float, a: Float, b: Float) -> Float {
     x.min(ONE_MINUS_EPSILON)
 }
 
-#[inline]
-pub fn float_to_bits(f: Float) -> u32 {
-    f.to_bits()
+#[inline(always)]
+pub fn bits_to_float(bits: FloatBits) -> Float {
+    Float::from_bits_val(bits)
 }
 
-#[inline]
-pub fn bits_to_float(ui: u32) -> Float {
-    f32::from_bits(ui)
+#[inline(always)]
+pub fn float_to_bits(v: Float) -> FloatBits {
+    v.to_bits_val()
 }
 
-#[inline]
-pub fn f64_to_bits(f: f64) -> u64 {
-    f.to_bits()
+#[inline(always)]
+pub fn exponent(v: Float) -> i32 {
+    v.exponent_val()
+}
+
+#[inline(always)]
+pub fn significand(v: Float) -> FloatBits {
+    v.significand_val()
+}
+
+#[inline(always)]
+pub fn sign_bit(v: Float) -> FloatBits {
+    v.sign_bit_val()
 }
 
 #[inline]
@@ -315,7 +426,36 @@ pub fn sample_tent(u: Float, r: Float) -> Float {
     }
 }
 
-fn round_up_pow2(mut n: i32) -> i32 {
+fn left_shift2(mut x: u64) -> u64 {
+    x &= 0xffffffff;
+    x = (x ^ (x << 16)) & 0x0000ffff0000ffff;
+    x = (x ^ (x << 8)) & 0x00ff00ff00ff00ff;
+    x = (x ^ (x << 4)) & 0x0f0f0f0f0f0f0f0f;
+    x = (x ^ (x << 2)) & 0x3333333333333333;
+    x = (x ^ (x << 1)) & 0x5555555555555555;
+    x
+}
+
+fn left_shift3(mut x: u32) -> u32 {
+    if x == (1 << 10) {
+        x -= 1;
+    }
+    x = (x | (x << 16)) & 0b00000011000000000000000011111111;
+    x = (x | (x << 8)) & 0b00000011000000001111000000001111;
+    x = (x | (x << 4)) & 0b00000011000011000011000011000011;
+    x = (x | (x << 2)) & 0b00001001001001001001001001001001;
+    x
+}
+
+pub fn encode_morton_2(x: u32, y: u32) -> u64 {
+    left_shift2(y as u64) << 1 | left_shift2(x as u64)
+}
+
+pub fn encode_morton_3(x: Float, y: Float, z: Float) -> u32 {
+    (left_shift3(x as u32) << 2) | (left_shift3(y as u32) << 1) | left_shift3(z as u32)
+}
+
+pub fn round_up_pow2(mut n: i32) -> i32 {
     if n <= 0 {
         return 1;
     }
@@ -328,6 +468,458 @@ fn round_up_pow2(mut n: i32) -> i32 {
     n + 1
 }
 
+pub const PRIME_TABLE_SIZE: usize = 1000;
+
+const PRIMES: [i32; PRIME_TABLE_SIZE] = [
+    2, 3, 5, 7, 11, // Subsequent prime numbers
+    13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107,
+    109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
+    223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317,
+    331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439,
+    443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569,
+    571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
+    683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821,
+    823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947,
+    953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061,
+    1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181,
+    1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291,
+    1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429,
+    1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531,
+    1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627,
+    1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759,
+    1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
+    1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017,
+    2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137,
+    2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273,
+    2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383,
+    2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
+    2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663,
+    2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753,
+    2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887,
+    2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023,
+    3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
+    3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313,
+    3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449,
+    3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557,
+    3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677,
+    3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
+    3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931,
+    3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079,
+    4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219,
+    4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349,
+    4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493,
+    4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639,
+    4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783,
+    4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931,
+    4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023,
+    5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179,
+    5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333,
+    5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471,
+    5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591,
+    5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737,
+    5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
+    5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029,
+    6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151,
+    6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287,
+    6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397,
+    6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571,
+    6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709,
+    6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857,
+    6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983,
+    6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129,
+    7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297,
+    7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477,
+    7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577,
+    7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703,
+    7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873,
+    7877, 7879, 7883, 7901, 7907, 7919,
+];
+
+#[inline]
+pub fn radical_inverse(base_index: usize, mut a: u64) -> Float {
+    let base = PRIMES[base_index] as u64;
+
+    let limit = (u64::MAX / base).saturating_sub(base);
+
+    let inv_base = (1.0 as Float) / (base as Float);
+    let mut inv_base_m = 1.0 as Float;
+    let mut reversed_digits = 0u64;
+
+    // Loop until we run out of digits or hit the overflow safety limit
+    while a > 0 && reversed_digits < limit {
+        // Rust's div_rem optimization handles / and % efficiently together
+        let next = a / base;
+        let digit = a % base;
+
+        reversed_digits = reversed_digits.wrapping_mul(base).wrapping_add(digit);
+        inv_base_m *= inv_base;
+        a = next;
+    }
+
+    // Ensure result is strictly less than 1.0
+    (reversed_digits as Float * inv_base_m).min(ONE_MINUS_EPSILON)
+}
+
+pub fn inverse_radical_inverse(mut inverse: u64, base: u64, n_digits: u64) -> u64 {
+    let mut index = 0;
+    for _ in 0..n_digits {
+        let digit = inverse % base;
+        inverse /= base;
+        index = index * base + digit;
+    }
+    index
+}
+
+// Digit scrambling
+#[derive(Default, Debug, Clone)]
+pub struct DigitPermutation {
+    base: usize,
+    n_digits: usize,
+    permutations: Vec<u16>,
+}
+
+impl DigitPermutation {
+    pub fn new(base: usize, seed: u64) -> Self {
+        let mut n_digits = 0;
+        let inv_base = 1. / base as Float;
+        let mut inv_base_m = 1.;
+
+        while 1.0 - ((base as Float - 1.0) * inv_base_m) < 1.0 {
+            n_digits += 1;
+            inv_base_m *= inv_base;
+        }
+
+        let mut permutations = vec![0u16; (n_digits * base) as usize];
+
+        for digit_index in 0..n_digits {
+            let hash_input = [base as u64, digit_index as u64, seed as u64];
+            let dseed = hash_buffer(&hash_input, 0);
+
+            for digit_value in 0..base {
+                let index = (digit_index * base + digit_value) as usize;
+
+                permutations[index] =
+                    permutation_element(digit_value as u32, base as u32, dseed as u32) as u16;
+            }
+        }
+
+        Self {
+            base,
+            n_digits,
+            permutations,
+        }
+    }
+
+    #[inline(always)]
+    pub fn permute(&self, digit_index: i32, digit_value: i32) -> i32 {
+        let idx = (digit_index * self.base as i32 + digit_value) as usize;
+        self.permutations[idx] as i32
+    }
+}
+
+pub fn compute_radical_inverse_permutations(seed: u64) -> Vec<DigitPermutation> {
+    PRIMES
+        .par_iter()
+        .map(|&base| DigitPermutation::new(base as usize, seed))
+        .collect()
+}
+
+pub fn scrambled_radical_inverse(base_index: usize, mut a: u64, perm: &DigitPermutation) -> Float {
+    // Access global PRIMES array (or pass base explicitly if preferred)
+    let base = PRIMES[base_index] as u64;
+
+    let limit = (u64::MAX / base).saturating_sub(base);
+
+    let inv_base = 1.0 / (base as Float);
+    let mut inv_base_m = 1.0;
+    let mut reversed_digits = 0u64;
+    let mut digit_index = 0;
+
+    // Note: 'a' is not part of the condition; we pad with 0s if a runs out.
+    while 1.0 - ((base as Float - 1.0) * inv_base_m) < 1.0 && reversed_digits < limit {
+        let next = a / base;
+        let digit_value = (a - next * base) as i32;
+
+        // Permute the digit
+        let permuted = perm.permute(digit_index, digit_value) as u64;
+
+        reversed_digits = reversed_digits.wrapping_mul(base).wrapping_add(permuted);
+        inv_base_m *= inv_base;
+        digit_index += 1;
+        a = next;
+    }
+
+    (inv_base_m * reversed_digits as Float).min(ONE_MINUS_EPSILON)
+}
+
+pub fn owen_scrambled_radical_inverse(base_index: usize, mut a: u64, hash: u32) -> Float {
+    let base = PRIMES[base_index] as u64;
+
+    let limit = (u64::MAX / base).saturating_sub(base);
+    let inv_base = 1.0 / (base as Float);
+    let mut inv_base_m = 1.0;
+    let mut reversed_digits = 0u64;
+
+    let mut _digit_index = 0;
+
+    while 1.0 - inv_base_m < 1.0 && reversed_digits < limit {
+        let next = a / base;
+        let digit_value = (a - next * base) as u32;
+
+        // Compute Owen-scrambled digit
+        // XOR the current seed (hash) with the accumulated reversed digits so far
+        let digit_hash = mix_bits((hash as u64) ^ reversed_digits) as u32;
+
+        let permuted = permutation_element(digit_value, base as u32, digit_hash) as u64;
+
+        reversed_digits = reversed_digits.wrapping_mul(base).wrapping_add(permuted);
+        inv_base_m *= inv_base;
+        _digit_index += 1;
+        a = next;
+    }
+
+    (inv_base_m * reversed_digits as Float).min(ONE_MINUS_EPSILON)
+}
+
+pub fn permutation_element(mut i: u32, l: u32, p: u32) -> u32 {
+    // 1. Compute the mask 'w' (next power of 2 minus 1)
+    let mut w = l - 1;
+    w |= w >> 1;
+    w |= w >> 2;
+    w |= w >> 4;
+    w |= w >> 8;
+    w |= w >> 16;
+
+    // 2. Cycle Walking loop
+    loop {
+        // We use wrapping_mul to replicate C++ unsigned overflow behavior
+        i ^= p;
+        i = i.wrapping_mul(0xe170893d);
+        i ^= p >> 16;
+        i ^= (i & w) >> 4;
+        i ^= p >> 8;
+        i = i.wrapping_mul(0x0929eb3f);
+        i ^= p >> 23;
+        i ^= (i & w) >> 1;
+        i = i.wrapping_mul(1 | (p >> 27));
+        i = i.wrapping_mul(0x6935fa69);
+        i ^= (i & w) >> 11;
+        i = i.wrapping_mul(0x74dcb303);
+        i ^= (i & w) >> 2;
+        i = i.wrapping_mul(0x9e501cc3);
+        i ^= (i & w) >> 2;
+        i = i.wrapping_mul(0xc860a3df);
+
+        i &= w;
+        i ^= i >> 5;
+
+        // If the resulting index is within range [0, l), we are done.
+        // Otherwise, we loop again (cycle walking).
+        if i < l {
+            break;
+        }
+    }
+
+    // 3. Final shift by p
+    (i.wrapping_add(p)) % l
+}
+
+pub fn multiply_generator(c: &[u32], mut a: u32) -> u32 {
+    let mut v = 0;
+    let mut i = 0;
+
+    while a != 0 && i < c.len() {
+        if (a & 1) != 0 {
+            v ^= c[i];
+        }
+        a >>= 1;
+        i += 1;
+    }
+    v
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct NoRandomizer;
+impl NoRandomizer {
+    pub fn scramble(&self, v: u32) -> u32 {
+        v
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct BinaryPermuteScrambler {
+    pub permutation: u32,
+}
+
+impl BinaryPermuteScrambler {
+    pub fn new(perm: u32) -> Self {
+        Self { permutation: perm }
+    }
+
+    #[inline(always)]
+    pub fn scramble(&self, v: u32) -> u32 {
+        self.permutation ^ v
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct FastOwenScrambler {
+    pub seed: u32,
+}
+
+impl FastOwenScrambler {
+    pub fn new(seed: u32) -> Self {
+        Self { seed }
+    }
+
+    #[inline(always)]
+    pub fn scramble(&self, mut v: u32) -> u32 {
+        v = v.reverse_bits();
+
+        // v ^= v * 0x3d20adea;
+        v ^= v.wrapping_mul(0x3d20adea);
+
+        // v += seed;
+        v = v.wrapping_add(self.seed);
+
+        // v *= (seed >> 16) | 1;
+        v = v.wrapping_mul((self.seed >> 16) | 1);
+
+        // v ^= v * 0x05526c56;
+        v ^= v.wrapping_mul(0x05526c56);
+
+        // v ^= v * 0x53a22864;
+        v ^= v.wrapping_mul(0x53a22864);
+
+        v.reverse_bits()
+    }
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct OwenScrambler {
+    pub seed: u32,
+}
+
+impl OwenScrambler {
+    pub fn new(seed: u32) -> Self {
+        Self { seed }
+    }
+
+    #[inline]
+    pub fn scramble(&self, mut v: u32) -> u32 {
+        if (self.seed & 1) != 0 {
+            v ^= 1 << 31;
+        }
+
+        for b in 1..32 {
+            let mask = (!0u32) << (32 - b);
+            let input = ((v & mask) ^ self.seed) as u64;
+            if (mix_bits(input) as u32 & (1 << b)) != 0 {
+                v ^= 1 << (31 - b);
+            }
+        }
+        v
+    }
+}
+
+pub trait Scrambler {
+    fn scramble(&self, v: u32) -> u32;
+}
+
+impl Scrambler for NoRandomizer {
+    fn scramble(&self, v: u32) -> u32 {
+        self.scramble(v)
+    }
+}
+
+impl Scrambler for BinaryPermuteScrambler {
+    fn scramble(&self, v: u32) -> u32 {
+        self.scramble(v)
+    }
+}
+impl Scrambler for FastOwenScrambler {
+    fn scramble(&self, v: u32) -> u32 {
+        self.scramble(v)
+    }
+}
+impl Scrambler for OwenScrambler {
+    fn scramble(&self, v: u32) -> u32 {
+        self.scramble(v)
+    }
+}
+
+impl<F: Fn(u32) -> u32> Scrambler for F {
+    fn scramble(&self, v: u32) -> u32 {
+        (*self)(v)
+    }
+}
+
+const N_SOBOL_DIMENSIONS: usize = 1024;
+const SOBOL_MATRIX_SIZE: usize = 52;
+#[inline]
+pub fn sobol_sample<S: Scrambler>(mut a: u64, dimension: usize, randomizer: S) -> Float {
+    debug_assert!(
+        dimension < N_SOBOL_DIMENSIONS,
+        "Sobol dimension out of bounds"
+    );
+
+    debug_assert!(a < (1u64 << SOBOL_MATRIX_SIZE), "Sobol index too large");
+
+    let mut v: u32 = 0;
+    let mut i = dimension * SOBOL_MATRIX_SIZE;
+
+    while a != 0 {
+        if (a & 1) != 0 {
+            v ^= SOBOL_MATRICES_32[i] as u32;
+        }
+        a >>= 1;
+        i += 1;
+    }
+
+    v = randomizer.scramble(v);
+
+    let float_val = (v as Float) * 2.3283064365386963e-10;
+
+    float_val.min(ONE_MINUS_EPSILON)
+}
+
+pub fn sobol_interval_to_index(m: u32, frame: u64, p: Point2i) -> u64 {
+    if m == 0 {
+        return frame;
+    }
+
+    let m2 = m << 1;
+    let mut index = frame << m2;
+
+    let mut delta = 0u64;
+    let mut current_frame = frame;
+    let mut c = 0;
+
+    while current_frame != 0 {
+        if (current_frame & 1) != 0 {
+            delta ^= VDC_SOBOL_MATRICES[(m - 1) as usize][c];
+        }
+        current_frame >>= 1;
+        c += 1;
+    }
+
+    let px = p.x() as u32 as u64;
+    let py = p.y() as u32 as u64;
+
+    let mut b = ((px << m) | py) ^ delta;
+
+    let mut c = 0;
+    while b != 0 {
+        if (b & 1) != 0 {
+            index ^= VDC_SOBOL_MATRICES_INV[(m - 1) as usize][c];
+        }
+        b >>= 1;
+        c += 1;
+    }
+
+    index
+}
+
 // MATRIX STUFF (TEST THOROUGHLY)
 #[derive(Debug, Copy, Clone)]
 pub struct Matrix<T, const R: usize, const C: usize> {
diff --git a/src/utils/mesh.rs b/src/utils/mesh.rs
index 7f1e0bb..50713a2 100644
--- a/src/utils/mesh.rs
+++ b/src/utils/mesh.rs
@@ -1,6 +1,6 @@
 use crate::core::pbrt::Float;
-use crate::core::sampler::PiecewiseConstant2D;
 use crate::geometry::{Normal3f, Point2f, Point3f, Vector3f};
+use crate::utils::sampling::PiecewiseConstant2D;
 use crate::utils::transform::Transform;
 use std::sync::Arc;
 
diff --git a/src/utils/mod.rs b/src/utils/mod.rs
index 2d9f898..c1abe3e 100644
--- a/src/utils/mod.rs
+++ b/src/utils/mod.rs
@@ -8,8 +8,25 @@ pub mod interval;
 pub mod math;
 pub mod mesh;
 pub mod quaternion;
+pub mod rng;
 pub mod sampling;
 pub mod scattering;
+pub mod sobol;
 pub mod spectrum;
 pub mod splines;
 pub mod transform;
+
+#[inline]
+pub fn partition_slice<T, F>(data: &mut [T], predicate: F) -> usize
+where
+    F: Fn(&T) -> bool,
+{
+    let mut i = 0;
+    for j in 0..data.len() {
+        if predicate(&data[j]) {
+            data.swap(i, j);
+            i += 1;
+        }
+    }
+    i
+}
diff --git a/src/utils/sampling.rs b/src/utils/sampling.rs
index 7c01086..27199fd 100644
--- a/src/utils/sampling.rs
+++ b/src/utils/sampling.rs
@@ -1,11 +1,15 @@
 use super::math::safe_sqrt;
 use crate::check_rare;
 use crate::core::pbrt::{
-    Float, INV_2_PI, INV_PI, ONE_MINUS_EPSILON, PI, PI_OVER_2, PI_OVER_4, clamp_t, lerp,
+    Float, INV_2_PI, INV_PI, ONE_MINUS_EPSILON, PI, PI_OVER_2, PI_OVER_4, clamp_t, find_interval,
+    lerp,
 };
 use crate::core::pbrt::{RARE_EVENT_CONDITION_MET, RARE_EVENT_TOTAL_CALLS};
-use crate::geometry::{Frame, Point2f, Point3f, Vector2f, Vector2i, Vector3f, VectorLike};
-use crate::utils::math::{difference_of_products, square, sum_of_products};
+use crate::geometry::{
+    Bounds2f, Frame, Point2f, Point2i, Point3f, Vector2f, Vector2i, Vector3f, VectorLike,
+};
+use crate::utils::containers::Array2D;
+use crate::utils::math::{difference_of_products, logistic, square, sum_of_products};
 use std::sync::atomic::{AtomicU64, Ordering as SyncOrdering};
 
 pub fn sample_linear(u: Float, a: Float, b: Float) -> Float {
@@ -425,6 +429,29 @@ pub fn invert_spherical_triangle_sample(
     Some(Point2f::new(clamp_t(u0, 0.0, 1.0), clamp_t(u1, 0.0, 1.0)))
 }
 
+pub fn sample_logistic(u: Float, s: Float) -> Float {
+    -s * (1. / u - 1.).ln()
+}
+
+pub fn invert_logistic_sample(x: Float, s: Float) -> Float {
+    1. / (1. + (-x / s).exp())
+}
+
+pub fn trimmed_logistic_pdf(x: Float, s: Float, a: Float, b: Float) -> Float {
+    if x < a || x > b {
+        return 0.;
+    }
+    let p = |val: Float| invert_logistic_sample(val, s);
+    logistic(x, s) / (p(b) - p(a))
+}
+
+pub fn sample_trimmed_logistic(u: Float, s: Float, a: Float, b: Float) -> Float {
+    let p = |val: Float| invert_logistic_sample(val, s);
+    let u = lerp(u, p(a), p(b));
+    let x = sample_logistic(u, s);
+    clamp_t(x, a, b)
+}
+
 pub fn uniform_hemisphere_pdf() -> Float {
     INV_2_PI
 }
@@ -496,6 +523,160 @@ pub struct PLSample {
     pub pdf: Float,
 }
 
+#[derive(Debug, Clone)]
+pub struct PiecewiseConstant1D {
+    pub func: Vec<Float>,
+    pub cdf: Vec<Float>,
+    pub min: Float,
+    pub max: Float,
+    pub func_integral: Float,
+}
+impl PiecewiseConstant1D {
+    pub fn new(f: &[Float]) -> Self {
+        Self::new_with_bounds(f, 0., 1.)
+    }
+
+    pub fn new_with_bounds(f: &[Float], min: Float, max: Float) -> Self {
+        assert!(max > min);
+        let n = f.len();
+        let mut func = Vec::with_capacity(n);
+        for &val in f {
+            func.push(val.abs());
+        }
+
+        let mut cdf = vec![0.; n + 1];
+        for i in 1..=n {
+            debug_assert!(func[i - 1] >= 0.);
+            cdf[i] = cdf[i - 1] + func[i - 1] * (max - min) / n as Float;
+        }
+
+        let func_integral = cdf[n];
+        if func_integral == 0. {
+            for i in 1..=n {
+                cdf[i] = i as Float / n as Float;
+            }
+        } else {
+            for i in 1..=n {
+                cdf[i] /= func_integral;
+            }
+        }
+
+        Self {
+            func,
+            cdf,
+            func_integral,
+            min,
+            max,
+        }
+    }
+
+    pub fn integral(&self) -> Float {
+        self.func_integral
+    }
+
+    pub fn size(&self) -> usize {
+        self.func.len()
+    }
+
+    pub fn sample(&self, u: Float) -> (Float, Float, usize) {
+        let o = find_interval(self.cdf.len(), |idx| self.cdf[idx] <= u);
+        let mut du = u - self.cdf[o];
+        if self.cdf[o + 1] - self.cdf[o] > 0. {
+            du /= self.cdf[o + 1] - self.cdf[o];
+        }
+        debug_assert!(!du.is_nan());
+        let value = lerp((o as Float + du) / self.size() as Float, self.min, self.max);
+        let pdf_val = if self.func_integral > 0. {
+            self.func[o] / self.func_integral
+        } else {
+            0.
+        };
+        (value, pdf_val, o)
+    }
+}
+#[derive(Debug, Clone)]
+pub struct PiecewiseConstant2D {
+    pub p_conditional_v: Vec<PiecewiseConstant1D>,
+    pub p_marginal: PiecewiseConstant1D,
+    pub domain: Bounds2f,
+}
+
+impl PiecewiseConstant2D {
+    pub fn new(data: &Array2D<Float>, nu: usize, nv: usize, domain: Bounds2f) -> Self {
+        let mut p_conditional_v = Vec::with_capacity(nv);
+        for v in 0..nv {
+            let start = v * nu;
+            let end = start + nu;
+            p_conditional_v.push(PiecewiseConstant1D::new_with_bounds(
+                &data.as_slice()[start..end],
+                domain.p_min.x(),
+                domain.p_max.x(),
+            ));
+        }
+
+        let marginal_func: Vec<Float> = p_conditional_v.iter().map(|p| p.integral()).collect();
+        let p_marginal = PiecewiseConstant1D::new_with_bounds(
+            &marginal_func,
+            domain.p_min.y(),
+            domain.p_max.y(),
+        );
+
+        Self {
+            p_conditional_v,
+            p_marginal,
+            domain,
+        }
+    }
+
+    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)
+    }
+
+    pub fn new_with_data(data: &Array2D<Float>, nx: usize, ny: usize) -> Self {
+        Self::new(
+            data,
+            nx,
+            ny,
+            Bounds2f::from_points(Point2f::new(0., 0.), Point2f::new(1., 1.)),
+        )
+    }
+
+    pub fn resolution(&self) -> Point2i {
+        Point2i::new(
+            self.p_conditional_v[0].size() as i32,
+            self.p_conditional_v[1].size() as i32,
+        )
+    }
+
+    pub fn integral(&self) -> f32 {
+        self.p_marginal.integral()
+    }
+
+    pub fn sample(&self, u: Point2f) -> (Point2f, f32, Point2i) {
+        let (d1, pdf1, off_y) = self.p_marginal.sample(u.y());
+        let (d0, pdf0, off_x) = self.p_conditional_v[off_y].sample(u.x());
+        let pdf = pdf0 * pdf1;
+        let offset = Point2i::new(off_x as i32, off_y as i32);
+        (Point2f::new(d0, d1), pdf, offset)
+    }
+
+    pub fn pdf(&self, p: Point2f) -> f32 {
+        let p_offset = self.domain.offset(&p);
+        let nu = self.p_conditional_v[0].size();
+        let nv = self.p_marginal.size();
+
+        let iu = (p_offset.x() * nu as f32).clamp(0.0, nu as f32 - 1.0) as usize;
+        let iv = (p_offset.y() * nv as f32).clamp(0.0, nv as f32 - 1.0) as usize;
+
+        let integral = self.p_marginal.integral();
+        if integral == 0.0 {
+            0.0
+        } else {
+            self.p_conditional_v[iv].func[iu] / integral
+        }
+    }
+}
+
 #[derive(Debug, Clone)]
 pub struct PiecewiseLinear2D<const N: usize> {
     size: Vector2i,
@@ -640,7 +821,7 @@ impl<const N: usize> PiecewiseLinear2D<N> {
                 &param_weights,
             )
         };
-        let row = Self::find_interval(marginal_size, |idx| fetch_marginal(idx) <= sample.y());
+        let row = find_interval(marginal_size, |idx| fetch_marginal(idx) <= sample.y());
         let marginal_cdf_row = fetch_marginal(row);
         sample[1] -= marginal_cdf_row;
         let r0 = self.lookup(
@@ -680,7 +861,7 @@ impl<const N: usize> PiecewiseLinear2D<N> {
             );
             (1.0 - sample.y()) * v0 + sample.y() * v1
         };
-        let col = Self::find_interval(self.size.x() as u32, |idx| {
+        let col = find_interval(self.size.x() as u32, |idx| {
             fetch_conditional(idx) <= sample.x()
         });
         sample[0] -= fetch_conditional(col);
@@ -825,22 +1006,6 @@ impl<const N: usize> PiecewiseLinear2D<N> {
         pdf * self.inv_patch_size.x() * self.inv_patch_size.y()
     }
 
-    fn find_interval(size: u32, pred: impl Fn(u32) -> bool) -> u32 {
-        let mut first = 1u32;
-        let mut size = size - 1;
-        while size > 0 {
-            let half = size >> 1;
-            let middle = first + half;
-            if pred(middle) {
-                first = middle + 1;
-                size -= half + 1;
-            } else {
-                size = half;
-            }
-        }
-        first.saturating_sub(1)
-    }
-
     fn get_slice_info(&self, params: [Float; N]) -> (u32, [(Float, Float); N]) {
         let mut param_weight = [(0.0, 0.0); N];
         let mut slice_offset = 0u32;
@@ -851,7 +1016,7 @@ impl<const N: usize> PiecewiseLinear2D<N> {
                 continue;
             }
 
-            let param_index = Self::find_interval(self.param_size[dim], |idx| {
+            let param_index = find_interval(self.param_size[dim], |idx| {
                 self.param_values[dim][idx as usize] <= params[dim]
             });
 
diff --git a/src/utils/spectrum.rs b/src/utils/spectrum.rs
index a2329e7..849a59c 100644
--- a/src/utils/spectrum.rs
+++ b/src/utils/spectrum.rs
@@ -93,6 +93,26 @@ impl SampledSpectrum {
         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 {
@@ -351,6 +371,9 @@ pub enum Spectrum {
 }
 
 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),
@@ -687,7 +710,52 @@ pub struct RGBSpectrum {
 }
 
 #[derive(Debug, Clone, Copy)]
-pub struct RGBUnboundedSpectrum(pub RGBSpectrum);
+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::*;