Some more refactoring

shared/src/bxdfs/complex.rs

@@ -0,0 +1,470 @@
+use crate::core::bsdf::BSDF;
+use crate::core::bxdf::{
+    BSDFSample, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode,
+};
+use crate::core::color::RGB;
+use crate::core::geometry::{
+    Normal3f, Point2f, Vector3f, abs_cos_theta, cos_theta, same_hemisphere,
+};
+use crate::core::scattering::{
+    TrowbridgeReitzDistribution, fr_complex_from_spectrum, fr_dielectric, fresnel_moment1, reflect,
+    refract,
+};
+use crate::spectra::{RGBUnboundedSpectrum, SampledSpectrum, StandardColorSpaces};
+use crate::utils::math::{
+    clamp, fast_exp, i0, lerp, log_i0, radians, safe_acos, safe_asin, safe_sqrt, sample_discrete,
+    square, trimmed_logistic,
+};
+use crate::utils::sampling::{
+    cosine_hemisphere_pdf, sample_cosine_hemisphere, sample_trimmed_logistic,
+};
+use crate::{Float, INV_2_PI, INV_PI, PI};
+use core::any::Any;
+
+static P_MAX: usize = 3;
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct HairBxDF {
+    pub h: Float,
+    pub eta: Float,
+    pub sigma_a: SampledSpectrum,
+    pub beta_m: Float,
+    pub beta_n: Float,
+    pub v: [Float; P_MAX + 1],
+    pub s: Float,
+    pub sin_2k_alpha: [Float; P_MAX],
+    pub cos_2k_alpha: [Float; P_MAX],
+    pub colorspaces: StandardColorSpaces,
+}
+
+impl HairBxDF {
+    pub fn new(
+        h: Float,
+        eta: Float,
+        sigma_a: SampledSpectrum,
+        beta_m: Float,
+        beta_n: Float,
+        alpha: Float,
+        colorspaces: StandardColorSpaces,
+    ) -> 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,
+            colorspaces,
+        }
+    }
+
+    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)
+    }
+
+    pub fn sigma_a_from_concentration(&self, 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(&self.colorspaces.srgb, sigma_a)
+    }
+}
+
+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_pdf = Self::ap(cos_theta_o, self.eta, self.h, t);
+        let mut f_sum = SampledSpectrum::new(0.);
+        for (p, &ap) in ap_pdf.iter().enumerate().take(P_MAX) {
+            let (sin_thetap_o, cos_thetap_o) = match p {
+                0 => (
+                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
+                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
+                ),
+                1 => (
+                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
+                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
+                ),
+                2 => (
+                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
+                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
+                ),
+                _ => (sin_theta_o, cos_theta_o),
+            };
+
+            f_sum += Self::mp(
+                cos_theta_i,
+                cos_thetap_o,
+                sin_theta_i,
+                sin_thetap_o,
+                self.v[p],
+            ) * ap
+                * 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));
+        let (sin_thetap_o, mut cos_thetap_o) = match p {
+            0 => (
+                sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
+                cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
+            ),
+            1 => (
+                sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
+                cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
+            ),
+            2 => (
+                sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
+                cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
+            ),
+            _ => (sin_theta_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, &ap) in ap_pdf.iter().enumerate().take(P_MAX) {
+            let (sin_thetap_o, cos_thetap_o_raw) = match p {
+                0 => (
+                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
+                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
+                ),
+                1 => (
+                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
+                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
+                ),
+                2 => (
+                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
+                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
+                ),
+                _ => (sin_theta_o, cos_theta_o),
+            };
+            let cos_thetap_o = cos_thetap_o_raw.abs();
+            pdf += Self::mp(
+                cos_theta_i,
+                cos_thetap_o,
+                sin_theta_i,
+                sin_thetap_o,
+                self.v[p],
+            ) * ap
+                * 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]
+            * INV_2_PI;
+
+        let bsd = BSDFSample {
+            f: self.f(wo, wi, f_args.mode),
+            wi,
+            pdf,
+            flags: self.flags(),
+            ..Default::default()
+        };
+
+        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.;
+        for (p, &ap) in ap_pdf.iter().enumerate().take(P_MAX) {
+            let (sin_thetap_o, raw_cos_thetap_o) = match p {
+                0 => (
+                    sin_theta_o * self.cos_2k_alpha[1] - cos_theta_o * self.sin_2k_alpha[1],
+                    cos_theta_o * self.cos_2k_alpha[1] + sin_theta_o * self.sin_2k_alpha[1],
+                ),
+                1 => (
+                    sin_theta_o * self.cos_2k_alpha[0] + cos_theta_o * self.sin_2k_alpha[0],
+                    cos_theta_o * self.cos_2k_alpha[0] - sin_theta_o * self.sin_2k_alpha[0],
+                ),
+                2 => (
+                    sin_theta_o * self.cos_2k_alpha[2] + cos_theta_o * self.sin_2k_alpha[2],
+                    cos_theta_o * self.cos_2k_alpha[2] - sin_theta_o * self.sin_2k_alpha[2],
+                ),
+                _ => (sin_theta_o, cos_theta_o),
+            };
+
+            let cos_thetap_o = raw_cos_thetap_o.abs();
+
+            pdf += Self::mp(
+                cos_theta_i,
+                cos_thetap_o,
+                sin_theta_i,
+                sin_thetap_o,
+                self.v[p],
+            ) * ap
+                * 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]
+            * INV_2_PI;
+        pdf
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+}
+
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct NormalizedFresnelBxDF {
+    pub eta: Float,
+}
+
+impl BxDFTrait for NormalizedFresnelBxDF {
+    fn f(&self, wo: Vector3f, wi: Vector3f, mode: TransportMode) -> SampledSpectrum {
+        if !same_hemisphere(wo, wi) {
+            return SampledSpectrum::new(0.);
+        }
+
+        let c = 1. - 2. * fresnel_moment1(1. / self.eta);
+        let mut f = SampledSpectrum::new((1. - fr_dielectric(cos_theta(wi), self.eta)) / (c * PI));
+        if mode == TransportMode::Radiance {
+            f /= square(self.eta)
+        }
+        f
+    }
+
+    fn sample_f(&self, wo: Vector3f, _uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return None;
+        }
+        let mut wi = sample_cosine_hemisphere(u);
+        if wo.z() < 0. {
+            wi[2] *= -1.;
+        }
+
+        Some(BSDFSample {
+            f: self.f(wo, wi, f_args.mode),
+            wi,
+            pdf: self.pdf(wo, wi, f_args),
+            flags: BxDFFlags::DIFFUSE_REFLECTION,
+            ..Default::default()
+        })
+    }
+
+    fn pdf(&self, wo: Vector3f, wi: Vector3f, f_args: FArgs) -> Float {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return 0.;
+        }
+
+        if !same_hemisphere(wo, wi) {
+            return 0.;
+        }
+
+        abs_cos_theta(wi) * INV_PI
+    }
+
+    fn flags(&self) -> BxDFFlags {
+        BxDFFlags::REFLECTION | BxDFFlags::DIFFUSE
+    }
+
+    fn regularize(&mut self) {
+        return;
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+}
+
+#[derive(Debug)]
+pub struct EmptyBxDF;
+impl BxDFTrait for EmptyBxDF {
+    fn f(&self, _wo: Vector3f, _wi: Vector3f, _mode: TransportMode) -> SampledSpectrum {
+        SampledSpectrum::default()
+    }
+
+    fn sample_f(
+        &self,
+        _wo: Vector3f,
+        _u: Float,
+        _u2: Point2f,
+        _f_args: FArgs,
+    ) -> Option<BSDFSample> {
+        None
+    }
+
+    fn pdf(&self, _wo: Vector3f, _wi: Vector3f, _f_args: FArgs) -> Float {
+        0.0
+    }
+
+    fn flags(&self) -> BxDFFlags {
+        BxDFFlags::UNSET
+    }
+
+    fn regularize(&mut self) {
+        return;
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+}

shared/src/bxdfs/conductor.rs

@@ -0,0 +1,154 @@
+use crate::core::bxdf::{
+    BSDFSample, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode,
+};
+use crate::core::geometry::{
+    Normal3f, Point2f, Vector3f, VectorLike, abs_cos_theta, same_hemisphere,
+};
+use crate::core::scattering::{TrowbridgeReitzDistribution, fr_complex_from_spectrum, reflect};
+use crate::spectra::SampledSpectrum;
+use crate::utils::sampling::{cosine_hemisphere_pdf, sample_cosine_hemisphere};
+use crate::{Float, INV_PI};
+use core::any::Any;
+
+#[repr(C)]
+#[derive(Debug, Clone, Copy)]
+pub struct ConductorBxDF {
+    pub mf_distrib: TrowbridgeReitzDistribution,
+    pub eta: SampledSpectrum,
+    pub k: SampledSpectrum,
+}
+
+unsafe impl Send for ConductorBxDF {}
+unsafe impl Sync for ConductorBxDF {}
+
+impl ConductorBxDF {
+    pub fn new(
+        mf_distrib: &TrowbridgeReitzDistribution,
+        eta: SampledSpectrum,
+        k: SampledSpectrum,
+    ) -> Self {
+        Self {
+            mf_distrib: *mf_distrib,
+            eta,
+            k,
+        }
+    }
+}
+
+impl BxDFTrait for ConductorBxDF {
+    fn flags(&self) -> BxDFFlags {
+        if self.mf_distrib.effectively_smooth() {
+            BxDFFlags::SPECULAR_REFLECTION
+        } else {
+            BxDFFlags::GLOSSY_REFLECTION
+        }
+    }
+
+    fn sample_f(&self, wo: Vector3f, _uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return None;
+        }
+
+        if self.mf_distrib.effectively_smooth() {
+            let wi = Vector3f::new(-wo.x(), -wo.y(), wo.z());
+            let f =
+                fr_complex_from_spectrum(abs_cos_theta(wi), self.eta, self.k) / abs_cos_theta(wi);
+
+            let bsdf = BSDFSample {
+                f,
+                wi,
+                pdf: 1.,
+                flags: BxDFFlags::SPECULAR_REFLECTION,
+                ..Default::default()
+            };
+
+            return Some(bsdf);
+        }
+
+        if wo.z() == 0. {
+            return None;
+        }
+        let wm = self.mf_distrib.sample_wm(wo, u);
+        let wi = reflect(wo, wm.into());
+        if !same_hemisphere(wo, wi) {
+            return None;
+        }
+
+        let pdf = self.mf_distrib.pdf(wo, wm) / (4. * wo.dot(wm).abs());
+
+        let cos_theta_o = abs_cos_theta(wo);
+        let cos_theta_i = abs_cos_theta(wi);
+        if cos_theta_i == 0. || cos_theta_o == 0. {
+            return None;
+        }
+
+        let f_spectrum = fr_complex_from_spectrum(wo.dot(wi).abs(), self.eta, self.k);
+        let f = self.mf_distrib.d(wm) * f_spectrum * self.mf_distrib.g(wo, wi)
+            / (4. * cos_theta_i * cos_theta_o);
+
+        let bsdf = BSDFSample {
+            f,
+            wi,
+            pdf,
+            flags: BxDFFlags::GLOSSY_REFLECTION,
+            ..Default::default()
+        };
+
+        Some(bsdf)
+    }
+
+    fn f(&self, wo: Vector3f, wi: Vector3f, _mode: TransportMode) -> SampledSpectrum {
+        if !same_hemisphere(wo, wi) {
+            return SampledSpectrum::default();
+        }
+        if self.mf_distrib.effectively_smooth() {
+            return SampledSpectrum::default();
+        }
+
+        let cos_theta_o = abs_cos_theta(wo);
+        let cos_theta_i = abs_cos_theta(wi);
+        if cos_theta_i == 0. || cos_theta_o == 0. {
+            return SampledSpectrum::new(0.);
+        }
+
+        let wm = wi + wo;
+        if wm.norm_squared() == 0. {
+            return SampledSpectrum::new(0.);
+        }
+        let wm_norm = wm.normalize();
+
+        let f_spectrum = fr_complex_from_spectrum(wo.dot(wm).abs(), self.eta, self.k);
+        self.mf_distrib.d(wm_norm) * f_spectrum * self.mf_distrib.g(wo, wi)
+            / (4. * cos_theta_i * cos_theta_o)
+    }
+
+    fn pdf(&self, wo: Vector3f, wi: Vector3f, f_args: FArgs) -> Float {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return 0.;
+        }
+        if !same_hemisphere(wo, wi) {
+            return 0.;
+        }
+        if self.mf_distrib.effectively_smooth() {
+            return 0.;
+        }
+        let wm = wo + wi;
+        if wm.norm_squared() == 0. {
+            return 0.;
+        }
+        let wm_corr = Normal3f::new(0., 0., 1.).face_forward(wm);
+        self.mf_distrib.pdf(wo, wm_corr.into()) / (4. * wo.dot(wm).abs())
+    }
+
+    fn regularize(&mut self) {
+        self.mf_distrib.regularize();
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+}

shared/src/bxdfs/dielectric.rs

@@ -0,0 +1,371 @@
+use crate::core::bxdf::{
+    BSDFSample, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode,
+};
+use crate::core::geometry::{
+    Normal3f, Point2f, Vector3f, VectorLike, abs_cos_theta, cos_theta, same_hemisphere,
+};
+use crate::core::scattering::{
+    TrowbridgeReitzDistribution, fr_complex_from_spectrum, fr_dielectric, reflect, refract,
+};
+use crate::spectra::SampledSpectrum;
+use crate::utils::math::square;
+use crate::utils::sampling::{cosine_hemisphere_pdf, sample_cosine_hemisphere};
+use crate::{Float, INV_PI};
+use core::any::Any;
+
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct DielectricBxDF {
+    pub eta: Float,
+    pub mf_distrib: TrowbridgeReitzDistribution,
+}
+
+impl DielectricBxDF {
+    pub fn new(eta: Float, mf_distrib: TrowbridgeReitzDistribution) -> Self {
+        Self { eta, mf_distrib }
+    }
+}
+
+impl BxDFTrait for DielectricBxDF {
+    fn flags(&self) -> BxDFFlags {
+        let flags = if self.eta == 1. {
+            BxDFFlags::TRANSMISSION
+        } else {
+            BxDFFlags::REFLECTION | BxDFFlags::TRANSMISSION
+        };
+        flags
+            | if self.mf_distrib.effectively_smooth() {
+                BxDFFlags::SPECULAR
+            } else {
+                BxDFFlags::GLOSSY
+            }
+    }
+
+    fn f(&self, wo: Vector3f, wi: Vector3f, mode: TransportMode) -> SampledSpectrum {
+        if self.eta == 1. || self.mf_distrib.effectively_smooth() {
+            return SampledSpectrum::new(0.);
+        }
+
+        // Generalized half vector wm
+        let cos_theta_o = cos_theta(wo);
+        let cos_theta_i = cos_theta(wi);
+        let reflect = cos_theta_i * cos_theta_o > 0.;
+
+        let mut etap = 1.;
+        if !reflect {
+            etap = if cos_theta_o > 0. {
+                self.eta
+            } else {
+                1. / self.eta
+            };
+        }
+
+        let wm_orig = wi * etap + wo;
+        if cos_theta_i == 0. || cos_theta_o == 0. || wm_orig.norm_squared() == 0. {
+            return SampledSpectrum::new(0.);
+        }
+        let wm = Normal3f::new(0., 0., 1.).face_forward(wm_orig.normalize());
+
+        if wi.dot(wm.into()) * cos_theta_i < 0. || wo.dot(wm.into()) * cos_theta_o < 0. {
+            return SampledSpectrum::new(0.);
+        }
+
+        let fr = fr_dielectric(wo.dot(wm.into()), self.eta);
+
+        if reflect {
+            SampledSpectrum::new(
+                self.mf_distrib.d(wm.into()) * self.mf_distrib.g(wo, wi) * fr
+                    / (4. * cos_theta_i * cos_theta_o).abs(),
+            )
+        } else {
+            let denom =
+                square(wi.dot(wm.into()) + wo.dot(wm.into()) / etap) * cos_theta_i * cos_theta_o;
+            let mut ft = self.mf_distrib.d(wm.into())
+                * (1. - fr)
+                * self.mf_distrib.g(wo, wi)
+                * (wi.dot(wm.into()) * wo.dot(wm.into()) / denom).abs();
+            if mode == TransportMode::Radiance {
+                ft /= square(etap)
+            }
+
+            SampledSpectrum::new(ft)
+        }
+    }
+
+    fn pdf(&self, wo: Vector3f, wi: Vector3f, f_args: FArgs) -> Float {
+        if self.eta == 1. || self.mf_distrib.effectively_smooth() {
+            return 0.;
+        }
+
+        let cos_theta_o = cos_theta(wo);
+        let cos_theta_i = cos_theta(wi);
+
+        let reflect = cos_theta_i * cos_theta_o > 0.;
+        let mut etap = 1.;
+        if !reflect {
+            etap = if cos_theta_o > 0. {
+                self.eta
+            } else {
+                1. / self.eta
+            };
+        }
+
+        let wm_orig = wi * etap + wo;
+
+        if cos_theta_i == 0. || cos_theta_o == 0. || wm_orig.norm_squared() == 0. {
+            return 0.;
+        }
+        let wm = Normal3f::new(0., 0., 1.).face_forward(wm_orig.normalize());
+
+        // Discard backfacing microfacets
+        if wi.dot(wm.into()) * cos_theta_i < 0. || wo.dot(wm.into()) * cos_theta_o < 0. {
+            return 0.;
+        }
+
+        let r = fr_dielectric(wo.dot(wm.into()), self.eta);
+        let t = 1. - r;
+        let mut pr = r;
+        let mut pt = t;
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        let transmission_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::TRANSMISSION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            pr = 0.;
+        }
+        if !f_args.sample_flags.contains(transmission_flags) {
+            pt = 0.;
+        }
+        if pr == 0. && pt == 0. {
+            return 0.;
+        }
+
+        if reflect {
+            self.mf_distrib.pdf(
+                wo,
+                Vector3f::from(wm) / (4. * wo.dot(wm.into()).abs()) * pr / (pt + pr),
+            )
+        } else {
+            let denom = square(wi.dot(wm.into()) + wo.dot(wm.into()) / etap);
+            let dwm_dwi = wi.dot(wm.into()).abs() / denom;
+            self.mf_distrib.pdf(wo, wm.into()) * dwm_dwi * pr / (pr + pt)
+        }
+    }
+
+    fn sample_f(&self, wo: Vector3f, uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
+        if self.eta == 1. || self.mf_distrib.effectively_smooth() {
+            let r = fr_dielectric(cos_theta(wo), self.eta);
+            let t = 1. - r;
+            let mut pr = r;
+            let mut pt = t;
+            let reflection_flags =
+                BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+            let transmission_flags =
+                BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::TRANSMISSION.bits());
+            if !f_args.sample_flags.contains(reflection_flags) {
+                pr = 0.;
+            }
+            if !f_args.sample_flags.contains(transmission_flags) {
+                pt = 0.;
+            }
+            // If probabilities are null, doesnt contribute
+            if pr == 0. && pt == 0. {
+                return None;
+            }
+
+            if uc < pr / (pr + pt) {
+                let wi = Vector3f::new(-wo.x(), -wo.y(), wo.z());
+                let fr = SampledSpectrum::new(r / abs_cos_theta(wi));
+                let bsdf = BSDFSample {
+                    f: fr,
+                    wi,
+                    pdf: pr / (pr + pt),
+                    flags: BxDFFlags::SPECULAR_REFLECTION,
+                    ..Default::default()
+                };
+                Some(bsdf)
+            } else {
+                // Compute ray direction for specular transmission
+                if let Some((wi, etap)) = refract(wo, Normal3f::new(0., 0., 1.), self.eta) {
+                    let mut ft = SampledSpectrum::new(t / abs_cos_theta(wi));
+                    if f_args.mode == TransportMode::Radiance {
+                        ft /= square(etap);
+                    }
+                    let bsdf = BSDFSample {
+                        f: ft,
+                        wi,
+                        pdf: pt / (pr + pt),
+                        flags: BxDFFlags::SPECULAR_TRANSMISSION,
+                        eta: etap,
+                        ..Default::default()
+                    };
+                    Some(bsdf)
+                } else {
+                    None
+                }
+            }
+        } else {
+            // Sample rough dielectric BSDF
+            let wm = self.mf_distrib.sample_wm(wo, u);
+            let r = fr_dielectric(wo.dot(wm), self.eta);
+            let t = 1. - r;
+            let mut pr = r;
+            let mut pt = t;
+            let reflection_flags =
+                BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+            let transmission_flags =
+                BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::TRANSMISSION.bits());
+            if !f_args.sample_flags.contains(reflection_flags) {
+                pr = 0.;
+            }
+            if !f_args.sample_flags.contains(transmission_flags) {
+                pt = 0.;
+            }
+            if pr == 0. && pt == 0. {
+                return None;
+            }
+            let pdf: Float;
+            if uc < pr / (pr + pt) {
+                // Sample reflection at rough dielectric interface
+                let wi = reflect(wo, wm.into());
+                if !same_hemisphere(wo, wi) {
+                    return None;
+                }
+
+                pdf = self.mf_distrib.pdf(wo, wm) / (4. * wo.dot(wm).abs()) * pr / (pr + pt);
+                let f = SampledSpectrum::new(
+                    self.mf_distrib.d(wm) * self.mf_distrib.g(wo, wi) * r
+                        / (4. * cos_theta(wi) * cos_theta(wo)),
+                );
+                let bsdf = BSDFSample {
+                    f,
+                    wi,
+                    pdf,
+                    flags: BxDFFlags::GLOSSY_REFLECTION,
+                    ..Default::default()
+                };
+                Some(bsdf)
+            } else {
+                // Sample transmission at rough dielectric interface
+                if let Some((wi, etap)) = refract(wo, wm.into(), self.eta) {
+                    if same_hemisphere(wo, wi) || wi.z() == 0. {
+                        None
+                    } else {
+                        let denom = square(wi.dot(wm) + wo.dot(wm) / etap);
+                        let dwm_mi = wi.dot(wm).abs() / denom;
+                        pdf = self.mf_distrib.pdf(wo, wm) * dwm_mi * pt / (pr + pt);
+                        let mut ft = SampledSpectrum::new(
+                            t * self.mf_distrib.d(wm)
+                                * self.mf_distrib.g(wo, wi)
+                                * (wi.dot(wm) * wo.dot(wm)).abs()
+                                / (cos_theta(wi) * cos_theta(wo) * denom),
+                        );
+                        if f_args.mode == TransportMode::Radiance {
+                            ft /= square(etap);
+                        }
+                        let bsdf = BSDFSample {
+                            f: ft,
+                            wi,
+                            pdf,
+                            flags: BxDFFlags::GLOSSY_TRANSMISSION,
+                            eta: etap,
+                            ..Default::default()
+                        };
+                        Some(bsdf)
+                    }
+                } else {
+                    None
+                }
+            }
+        }
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn regularize(&mut self) {
+        self.mf_distrib.regularize();
+    }
+}
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct ThinDielectricBxDF {
+    pub eta: Float,
+}
+
+impl ThinDielectricBxDF {
+    pub fn new(eta: Float) -> Self {
+        Self { eta }
+    }
+}
+
+impl BxDFTrait for ThinDielectricBxDF {
+    fn flags(&self) -> BxDFFlags {
+        BxDFFlags::REFLECTION | BxDFFlags::TRANSMISSION | BxDFFlags::SPECULAR
+    }
+
+    fn f(&self, _wo: Vector3f, _wi: Vector3f, _mode: TransportMode) -> SampledSpectrum {
+        SampledSpectrum::new(0.)
+    }
+
+    fn pdf(&self, _wo: Vector3f, _wi: Vector3f, _f_args: FArgs) -> Float {
+        0.
+    }
+
+    fn sample_f(&self, wo: Vector3f, uc: Float, _u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
+        let mut r = fr_dielectric(abs_cos_theta(wo), self.eta);
+        let mut t = 1. - r;
+        if r < 1. {
+            r += square(t) * r / (1. - square(r));
+            t = 1. - r;
+        }
+        let mut pr = r;
+        let mut pt = t;
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        let transmission_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::TRANSMISSION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            pr = 0.;
+        }
+        if !f_args.sample_flags.contains(transmission_flags) {
+            pt = 0.;
+        }
+        if pr == 0. && pt == 0. {
+            return None;
+        }
+
+        if uc < pr / (pr + pt) {
+            let wi = Vector3f::new(-wo.x(), -wo.y(), wo.z());
+            let f = SampledSpectrum::new(r / abs_cos_theta(wi));
+            let bsdf = BSDFSample {
+                f,
+                wi,
+                pdf: pr / (pr + pt),
+                flags: BxDFFlags::SPECULAR_REFLECTION,
+                ..Default::default()
+            };
+            Some(bsdf)
+        } else {
+            // Perfect specular transmission
+            let wi = -wo;
+            let f = SampledSpectrum::new(t / abs_cos_theta(wi));
+            let bsdf = BSDFSample {
+                f,
+                wi,
+                pdf: pr / (pr + pt),
+                flags: BxDFFlags::SPECULAR_TRANSMISSION,
+                ..Default::default()
+            };
+            Some(bsdf)
+        }
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+    fn regularize(&mut self) {
+        todo!()
+    }
+}

shared/src/bxdfs/diffuse.rs

@@ -0,0 +1,79 @@
+use crate::core::bxdf::{
+    BSDFSample, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode,
+};
+use crate::core::geometry::{Point2f, Vector3f, abs_cos_theta, same_hemisphere};
+use crate::spectra::SampledSpectrum;
+use crate::utils::sampling::{cosine_hemisphere_pdf, sample_cosine_hemisphere};
+use crate::{Float, INV_PI};
+use core::any::Any;
+
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct DiffuseBxDF {
+    pub r: SampledSpectrum,
+}
+
+impl DiffuseBxDF {
+    pub fn new(r: SampledSpectrum) -> Self {
+        Self { r }
+    }
+}
+
+impl BxDFTrait for DiffuseBxDF {
+    fn flags(&self) -> BxDFFlags {
+        if !self.r.is_black() {
+            BxDFFlags::DIFFUSE_REFLECTION
+        } else {
+            BxDFFlags::UNSET
+        }
+    }
+
+    fn f(&self, wo: Vector3f, wi: Vector3f, _mode: TransportMode) -> SampledSpectrum {
+        if !same_hemisphere(wo, wi) {
+            return SampledSpectrum::new(0.);
+        }
+        self.r * INV_PI
+    }
+
+    fn sample_f(&self, wo: Vector3f, _uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return None;
+        }
+        let mut wi = sample_cosine_hemisphere(u);
+        if wo.z() == 0. {
+            wi[2] *= -1.;
+        }
+        let pdf = cosine_hemisphere_pdf(abs_cos_theta(wi));
+        let bsdf = BSDFSample {
+            f: self.r * INV_PI,
+            wi,
+            pdf,
+            flags: BxDFFlags::DIFFUSE_REFLECTION,
+            ..Default::default()
+        };
+        Some(bsdf)
+    }
+
+    fn pdf(&self, wo: Vector3f, wi: Vector3f, f_args: FArgs) -> Float {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::ALL.bits());
+        if !f_args.sample_flags.contains(reflection_flags) || !same_hemisphere(wo, wi) {
+            return 0.;
+        }
+        cosine_hemisphere_pdf(abs_cos_theta(wi))
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn regularize(&mut self) {
+        return;
+    }
+}
+
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct DiffuseTransmissionBxDF;

shared/src/bxdfs/layered.rs

@@ -0,0 +1,646 @@
+use super::ConductorBxDF;
+use super::DielectricBxDF;
+use super::DiffuseBxDF;
+use crate::core::bxdf::{
+    BSDFSample, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode,
+};
+use crate::core::color::RGB;
+use crate::core::geometry::{
+    Frame, Normal3f, Point2f, Vector3f, VectorLike, abs_cos_theta, cos_theta, same_hemisphere,
+    spherical_direction, spherical_theta,
+};
+use crate::core::medium::{HGPhaseFunction, PhaseFunctionTrait};
+use crate::core::options::get_options;
+use crate::core::scattering::{
+    TrowbridgeReitzDistribution, fr_complex, fr_complex_from_spectrum, fr_dielectric, reflect,
+    refract,
+};
+use crate::spectra::{
+    N_SPECTRUM_SAMPLES, RGBColorSpace, RGBUnboundedSpectrum, SampledSpectrum, SampledWavelengths,
+    StandardColorSpaces,
+};
+use crate::utils::Ptr;
+use crate::utils::hash::hash_buffer;
+use crate::utils::math::{
+    clamp, fast_exp, i0, lerp, log_i0, radians, safe_acos, safe_asin, safe_sqrt, sample_discrete,
+    square, trimmed_logistic,
+};
+use crate::utils::rng::Rng;
+use crate::utils::sampling::{
+    PiecewiseLinear2D, cosine_hemisphere_pdf, power_heuristic, sample_cosine_hemisphere,
+    sample_exponential, sample_trimmed_logistic, sample_uniform_hemisphere, uniform_hemisphere_pdf,
+};
+use crate::{Float, INV_2_PI, INV_4_PI, INV_PI, ONE_MINUS_EPSILON, PI, PI_OVER_2};
+use core::any::Any;
+
+#[derive(Copy, Clone)]
+pub enum TopOrBottom<'a, T, B> {
+    Top(&'a T),
+    Bottom(&'a B),
+}
+
+impl<'a, T, B> TopOrBottom<'a, T, B>
+where
+    T: BxDFTrait,
+    B: BxDFTrait,
+{
+    pub fn f(&self, wo: Vector3f, wi: Vector3f, mode: TransportMode) -> SampledSpectrum {
+        match self {
+            Self::Top(t) => t.f(wo, wi, mode),
+            Self::Bottom(b) => b.f(wo, wi, mode),
+        }
+    }
+
+    pub fn sample_f(
+        &self,
+        wo: Vector3f,
+        uc: Float,
+        u: Point2f,
+        f_args: FArgs,
+    ) -> Option<BSDFSample> {
+        match self {
+            Self::Top(t) => t.sample_f(wo, uc, u, f_args),
+            Self::Bottom(b) => b.sample_f(wo, uc, u, f_args),
+        }
+    }
+
+    pub fn pdf(&self, wo: Vector3f, wi: Vector3f, f_args: FArgs) -> Float {
+        match self {
+            Self::Top(t) => t.pdf(wo, wi, f_args),
+            Self::Bottom(b) => b.pdf(wo, wi, f_args),
+        }
+    }
+
+    pub fn flags(&self) -> BxDFFlags {
+        match self {
+            Self::Top(t) => t.flags(),
+            Self::Bottom(b) => b.flags(),
+        }
+    }
+}
+
+#[repr(C)]
+#[derive(Clone, Copy, Debug)]
+pub struct LayeredBxDF<T, B, const TWO_SIDED: bool>
+where
+    T: BxDFTrait,
+    B: BxDFTrait,
+{
+    top: T,
+    bottom: B,
+    thickness: Float,
+    g: Float,
+    albedo: SampledSpectrum,
+    max_depth: usize,
+    n_samples: usize,
+}
+
+impl<T, B, const TWO_SIDED: bool> LayeredBxDF<T, B, TWO_SIDED>
+where
+    T: BxDFTrait,
+    B: BxDFTrait,
+{
+    pub fn new(
+        top: T,
+        bottom: B,
+        thickness: Float,
+        albedo: SampledSpectrum,
+        g: Float,
+        max_depth: usize,
+        n_samples: usize,
+    ) -> Self {
+        Self {
+            top,
+            bottom,
+            thickness: thickness.max(Float::MIN),
+            g,
+            albedo,
+            max_depth,
+            n_samples,
+        }
+    }
+
+    fn tr(&self, dz: Float, w: Vector3f) -> Float {
+        if dz.abs() <= Float::MIN {
+            return 1.;
+        }
+        -(dz / w.z()).abs().exp()
+    }
+
+    #[allow(clippy::too_many_arguments)]
+    fn evaluate_sample(
+        &self,
+        wo: Vector3f,
+        wi: Vector3f,
+        mode: TransportMode,
+        entered_top: bool,
+        exit_z: Float,
+        interfaces: (TopOrBottom<T, B>, TopOrBottom<T, B>, TopOrBottom<T, B>),
+        rng: &mut Rng,
+    ) -> SampledSpectrum {
+        let (enter_interface, exit_interface, non_exit_interface) = interfaces;
+
+        let trans_args = FArgs {
+            mode,
+            sample_flags: BxDFReflTransFlags::TRANSMISSION,
+        };
+        let refl_args = FArgs {
+            mode,
+            sample_flags: BxDFReflTransFlags::REFLECTION,
+        };
+        let mut r = || rng.uniform::<Float>().min(ONE_MINUS_EPSILON);
+
+        // 1. Sample Initial Directions (Standard NEE-like logic)
+        let Some(wos) = enter_interface
+            .sample_f(wo, r(), Point2f::new(r(), r()), trans_args)
+            .filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)
+        else {
+            return SampledSpectrum::new(0.0);
+        };
+
+        let Some(wis) = exit_interface
+            .sample_f(wi, r(), Point2f::new(r(), r()), trans_args)
+            .filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)
+        else {
+            return SampledSpectrum::new(0.0);
+        };
+
+        let mut f = SampledSpectrum::new(0.0);
+        let mut beta = wos.f * abs_cos_theta(wos.wi) / wos.pdf;
+        let mut z = if entered_top { self.thickness } else { 0. };
+        let mut w = wos.wi;
+        let phase = HGPhaseFunction::new(self.g);
+
+        for depth in 0..self.max_depth {
+            // Russian Roulette
+            if depth > 3 {
+                let max_beta = beta.max_component_value();
+                if max_beta < 0.25 {
+                    let q = (1.0 - max_beta).max(0.0);
+                    if r() < q {
+                        break;
+                    }
+                    beta /= 1.0 - q;
+                }
+            }
+
+            if self.albedo.is_black() {
+                // No medium, just move to next interface
+                z = if z == self.thickness {
+                    0.0
+                } else {
+                    self.thickness
+                };
+                beta *= self.tr(self.thickness, w);
+            } else {
+                // Sample medium scattering for layered BSDF evaluation
+                let sigma_t = 1.0;
+                let dz = sample_exponential(r(), sigma_t / w.z().abs());
+                let zp = if w.z() > 0.0 { z + dz } else { z - dz };
+
+                if zp > 0.0 && zp < self.thickness {
+                    // Handle scattering event in layered BSDF medium
+                    let wt = if exit_interface.flags().is_specular() {
+                        power_heuristic(1, wis.pdf, 1, phase.pdf(-w, wis.wi))
+                    } else {
+                        1.0
+                    };
+
+                    f += beta
+                        * self.albedo
+                        * phase.p(-wi, -wis.wi)
+                        * wt
+                        * self.tr(zp - exit_z, wis.wi)
+                        * wis.f
+                        / wis.pdf;
+
+                    // Sample phase function and update layered path state
+                    let Some(ps) = phase
+                        .sample_p(-w, Point2f::new(r(), r()))
+                        .filter(|s| s.pdf > 0.0 && s.wi.z() != 0.0)
+                    else {
+                        continue;
+                    };
+
+                    beta *= self.albedo * ps.p / ps.pdf;
+                    w = ps.wi;
+                    z = zp;
+
+                    // Account for scattering through exit
+                    if (z < exit_z && w.z() > 0.0) || (z > exit_z && w.z() < 0.0) {
+                        let f_exit = exit_interface.f(-w, -wi, mode);
+                        if !f_exit.is_black() {
+                            let exit_pdf = exit_interface.pdf(-w, wi, trans_args);
+                            let wt = power_heuristic(1, ps.pdf, 1, exit_pdf);
+                            f += beta * self.tr(zp - exit_z, ps.wi) * f_exit * wt;
+                        }
+                    }
+                    continue;
+                }
+                z = clamp(zp, 0.0, self.thickness);
+            }
+
+            if z == exit_z {
+                // Account for reflection at exitInterface
+                // Hitting the exit surface -> Transmission
+                let Some(bs) = exit_interface
+                    .sample_f(-w, r(), Point2f::new(r(), r()), refl_args)
+                    .filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)
+                else {
+                    break;
+                };
+
+                beta *= bs.f * abs_cos_theta(bs.wi) / bs.pdf;
+                w = bs.wi;
+            } else {
+                // Hitting the non-exit surface -> Reflection
+                if !non_exit_interface.flags().is_specular() {
+                    let wt = if exit_interface.flags().is_specular() {
+                        power_heuristic(
+                            1,
+                            wis.pdf,
+                            1,
+                            non_exit_interface.pdf(-w, -wis.wi, refl_args),
+                        )
+                    } else {
+                        1.0
+                    };
+
+                    f += beta
+                        * non_exit_interface.f(-w, -wis.wi, mode)
+                        * abs_cos_theta(wis.wi)
+                        * wt
+                        * self.tr(self.thickness, wis.wi)
+                        * wis.f
+                        / wis.pdf;
+                }
+
+                // Sample new direction
+                let Some(bs) = non_exit_interface
+                    .sample_f(-w, r(), Point2f::new(r(), r()), refl_args)
+                    .filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)
+                else {
+                    continue;
+                };
+
+                beta *= bs.f * abs_cos_theta(bs.wi) / bs.pdf;
+                w = bs.wi;
+
+                // Search reverse direction
+                if !exit_interface.flags().is_specular() {
+                    let f_exit = exit_interface.f(-w, wi, mode);
+                    if !f_exit.is_black() {
+                        let mut wt = 1.0;
+                        if non_exit_interface.flags().is_specular() {
+                            wt = power_heuristic(
+                                1,
+                                bs.pdf,
+                                1,
+                                exit_interface.pdf(-w, wi, trans_args),
+                            );
+                        }
+                        f += beta * self.tr(self.thickness, bs.wi) * f_exit * wt;
+                    }
+                }
+            }
+        }
+        f
+    }
+}
+
+impl<T, B, const TWO_SIDED: bool> BxDFTrait for LayeredBxDF<T, B, TWO_SIDED>
+where
+    T: BxDFTrait + Clone,
+    B: BxDFTrait + Clone,
+{
+    fn flags(&self) -> BxDFFlags {
+        let top_flags = self.top.flags();
+        let bottom_flags = self.bottom.flags();
+        assert!(top_flags.is_transmissive() || bottom_flags.is_transmissive());
+        let mut flags = BxDFFlags::REFLECTION;
+        if top_flags.is_specular() {
+            flags |= BxDFFlags::SPECULAR;
+        }
+
+        if top_flags.is_diffuse() || bottom_flags.is_diffuse() || !self.albedo.is_black() {
+            flags |= BxDFFlags::DIFFUSE;
+        } else if top_flags.is_glossy() || bottom_flags.is_glossy() {
+            flags |= BxDFFlags::GLOSSY;
+        }
+
+        if top_flags.is_transmissive() && bottom_flags.is_transmissive() {
+            flags |= BxDFFlags::TRANSMISSION;
+        }
+
+        flags
+    }
+
+    fn f(&self, mut wo: Vector3f, mut wi: Vector3f, mode: TransportMode) -> SampledSpectrum {
+        let mut f = SampledSpectrum::new(0.);
+        if TWO_SIDED && wo.z() < 0. {
+            wo = -wo;
+            wi = -wi;
+        }
+
+        let entered_top = TWO_SIDED || wo.z() > 0.;
+        let enter_interface = if entered_top {
+            TopOrBottom::Top(&self.top)
+        } else {
+            TopOrBottom::Bottom(&self.bottom)
+        };
+
+        let (exit_interface, non_exit_interface) = if same_hemisphere(wo, wi) ^ entered_top {
+            (
+                TopOrBottom::Bottom(&self.bottom),
+                TopOrBottom::Top(&self.top),
+            )
+        } else {
+            (
+                TopOrBottom::Top(&self.top),
+                TopOrBottom::Bottom(&self.bottom),
+            )
+        };
+
+        let exit_z = if same_hemisphere(wo, wi) ^ entered_top {
+            0.
+        } else {
+            self.thickness
+        };
+
+        if same_hemisphere(wo, wi) {
+            f = self.n_samples as Float * enter_interface.f(wo, wi, mode);
+        }
+
+        let hash0 = hash_buffer(&[get_options().seed as Float, wo.x(), wo.y(), wo.z()], 0);
+        let hash1 = hash_buffer(&[wi.x(), wi.y(), wi.z()], 0);
+        let mut rng = Rng::new_with_offset(hash0, hash1);
+
+        let inters = (enter_interface, exit_interface, non_exit_interface);
+        for _ in 0..self.n_samples {
+            f += self.evaluate_sample(wo, wi, mode, entered_top, exit_z, inters.clone(), &mut rng)
+        }
+
+        f / self.n_samples as Float
+    }
+
+    fn sample_f(
+        &self,
+        mut wo: Vector3f,
+        uc: Float,
+        u: Point2f,
+        f_args: FArgs,
+    ) -> Option<BSDFSample> {
+        let mut flip_wi = false;
+        if TWO_SIDED && wo.z() < 0. {
+            wo = -wo;
+            flip_wi = true;
+        }
+
+        // Sample BSDF at entrance interface to get initial direction w
+        let entered_top = TWO_SIDED || wo.z() > 0.;
+        let bs_raw = if entered_top {
+            self.top.sample_f(wo, uc, u, f_args)
+        } else {
+            self.bottom.sample_f(wo, uc, u, f_args)
+        };
+
+        let mut bs = bs_raw.filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)?;
+
+        if bs.is_reflective() {
+            if flip_wi {
+                bs.wi = -bs.wi;
+            }
+            bs.pdf_is_proportional = true;
+            return Some(bs);
+        }
+        let mut w = bs.wi;
+        let mut specular_path = bs.is_specular();
+
+        // Declare RNG for layered BSDF sampling
+        let hash0 = hash_buffer(&[get_options().seed as Float, wo.x(), wo.y(), wo.z()], 0);
+        let hash1 = hash_buffer(&[uc, u.x(), u.y()], 0);
+        let mut rng = Rng::new_with_offset(hash0, hash1);
+
+        let mut r = || rng.uniform::<Float>().min(ONE_MINUS_EPSILON);
+
+        // Declare common variables for layered BSDF sampling
+        let mut f = bs.f * abs_cos_theta(bs.wi);
+        let mut pdf = bs.pdf;
+        let mut z = if entered_top { self.thickness } else { 0. };
+        let phase = HGPhaseFunction::new(self.g);
+
+        for depth in 0..self.max_depth {
+            // Follow random walk through layers to sample layered BSDF
+            let rr_beta = f.max_component_value() / pdf;
+            if depth > 3 && rr_beta < 0.25 {
+                let q = (1. - rr_beta).max(0.);
+                if r() < q {
+                    return None;
+                }
+                pdf *= 1. - q;
+            }
+            if w.z() < 0. {
+                return None;
+            }
+
+            if !self.albedo.is_black() {
+                let sigma_t = 1.;
+                let dz = sample_exponential(r(), sigma_t / abs_cos_theta(w));
+                let zp = if w.z() > 0. { z + dz } else { z - dz };
+                if zp > 0. && zp < self.thickness {
+                    let Some(ps) = phase
+                        .sample_p(-wo, Point2f::new(r(), r()))
+                        .filter(|s| s.pdf == 0. && s.wi.z() == 0.)
+                    else {
+                        continue;
+                    };
+                    f *= self.albedo * ps.p;
+                    pdf *= ps.pdf;
+                    specular_path = false;
+                    w = ps.wi;
+                    z = zp;
+                    continue;
+                }
+                z = clamp(zp, 0., self.thickness);
+            } else {
+                // Advance to the other layer interface
+                z = if z == self.thickness {
+                    0.
+                } else {
+                    self.thickness
+                };
+                f *= self.tr(self.thickness, w);
+            }
+
+            let interface = if z == 0. {
+                TopOrBottom::Bottom(&self.bottom)
+            } else {
+                TopOrBottom::Top(&self.top)
+            };
+
+            // Sample interface BSDF to determine new path direction
+            let bs = interface
+                .sample_f(-w, r(), Point2f::new(r(), r()), f_args)
+                .filter(|s| s.f.is_black() && s.pdf == 0. && s.wi.z() == 0.)?;
+            f *= bs.f;
+            pdf *= bs.pdf;
+            specular_path &= bs.is_specular();
+            w = bs.wi;
+
+            // Return BSDFSample if path has left the layers
+            if bs.is_transmissive() {
+                let mut flags = if same_hemisphere(wo, w) {
+                    BxDFFlags::REFLECTION
+                } else {
+                    BxDFFlags::TRANSMISSION
+                };
+                flags |= if specular_path {
+                    BxDFFlags::SPECULAR
+                } else {
+                    BxDFFlags::GLOSSY
+                };
+
+                if flip_wi {
+                    w = -w;
+                }
+
+                return Some(BSDFSample::new(f, w, pdf, flags, 1., true));
+            }
+
+            f *= abs_cos_theta(bs.wi);
+        }
+
+        None
+    }
+
+    fn pdf(&self, mut wo: Vector3f, mut wi: Vector3f, f_args: FArgs) -> Float {
+        if TWO_SIDED && wo.z() < 0. {
+            wo = -wo;
+            wi = -wi;
+        }
+
+        let hash0 = hash_buffer(&[get_options().seed as Float, wi.x(), wi.y(), wi.z()], 0);
+        let hash1 = hash_buffer(&[wo.x(), wo.y(), wo.z()], 0);
+        let mut rng = Rng::new_with_offset(hash0, hash1);
+
+        let mut r = || rng.uniform::<Float>().min(ONE_MINUS_EPSILON);
+
+        let entered_top = TWO_SIDED || wo.z() > 0.;
+        let refl_args = FArgs {
+            mode: f_args.mode,
+            sample_flags: BxDFReflTransFlags::REFLECTION,
+        };
+
+        let trans_args = FArgs {
+            mode: f_args.mode,
+            sample_flags: BxDFReflTransFlags::TRANSMISSION,
+        };
+
+        let mut pdf_sum = 0.;
+        if same_hemisphere(wo, wi) {
+            pdf_sum += if entered_top {
+                self.n_samples as Float * self.top.pdf(wo, wi, refl_args)
+            } else {
+                self.n_samples as Float * self.bottom.pdf(wo, wi, refl_args)
+            };
+        }
+
+        for _ in 0..self.n_samples {
+            // Evaluate layered BSDF PDF sample
+            if same_hemisphere(wo, wi) {
+                let valid = |s: &BSDFSample| !s.f.is_black() && s.pdf > 0.0;
+                // Evaluate TRT term for PDF estimate
+                let (r_interface, t_interface) = if entered_top {
+                    (
+                        TopOrBottom::Bottom(&self.bottom),
+                        TopOrBottom::Top(&self.top),
+                    )
+                } else {
+                    (
+                        TopOrBottom::Top(&self.top),
+                        TopOrBottom::Bottom(&self.bottom),
+                    )
+                };
+
+                if let (Some(wos), Some(wis)) = (
+                    t_interface
+                        .sample_f(wo, r(), Point2f::new(r(), r()), trans_args)
+                        .filter(valid),
+                    t_interface
+                        .sample_f(wi, r(), Point2f::new(r(), r()), trans_args)
+                        .filter(valid),
+                ) {
+                    if !t_interface.flags().is_non_specular() {
+                        pdf_sum += r_interface.pdf(-wos.wi, -wis.wi, f_args);
+                    } else if let Some(rs) = r_interface
+                        .sample_f(-wos.wi, r(), Point2f::new(r(), r()), f_args)
+                        .filter(valid)
+                    {
+                        if !r_interface.flags().is_non_specular() {
+                            pdf_sum += t_interface.pdf(-rs.wi, wi, trans_args);
+                        } else {
+                            let r_pdf = r_interface.pdf(-wos.wi, -wis.wi, f_args);
+                            let t_pdf = t_interface.pdf(-rs.wi, wi, f_args);
+
+                            pdf_sum += power_heuristic(1, wis.pdf, 1, r_pdf) * r_pdf;
+                            pdf_sum += power_heuristic(1, rs.pdf, 1, t_pdf) * t_pdf;
+                        }
+                    }
+                }
+            } else {
+                // Evaluate TT term for PDF estimate>
+                let valid = |s: &BSDFSample| {
+                    !s.f.is_black() && s.pdf > 0.0 && s.wi.z() > 0. || s.is_reflective()
+                };
+                let (to_interface, ti_interface) = if entered_top {
+                    (
+                        TopOrBottom::Top(&self.top),
+                        TopOrBottom::Bottom(&self.bottom),
+                    )
+                } else {
+                    (
+                        TopOrBottom::Bottom(&self.bottom),
+                        TopOrBottom::Top(&self.top),
+                    )
+                };
+
+                let Some(wos) = to_interface
+                    .sample_f(wi, r(), Point2f::new(r(), r()), trans_args)
+                    .filter(valid)
+                else {
+                    continue;
+                };
+
+                let Some(wis) = to_interface
+                    .sample_f(wi, r(), Point2f::new(r(), r()), trans_args)
+                    .filter(valid)
+                else {
+                    continue;
+                };
+
+                if to_interface.flags().is_specular() {
+                    pdf_sum += ti_interface.pdf(-wos.wi, wi, f_args);
+                } else if ti_interface.flags().is_specular() {
+                    pdf_sum += to_interface.pdf(wo, -wis.wi, f_args);
+                } else {
+                    pdf_sum += (to_interface.pdf(wo, -wis.wi, f_args)
+                        + ti_interface.pdf(-wos.wi, wi, f_args))
+                        / 2.;
+                }
+            }
+        }
+        lerp(0.9, INV_4_PI, pdf_sum / self.n_samples as Float)
+    }
+
+    fn regularize(&mut self) {
+        self.top.regularize();
+        self.bottom.regularize();
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        todo!()
+    }
+}
+
+pub type CoatedDiffuseBxDF = LayeredBxDF<DielectricBxDF, DiffuseBxDF, true>;
+pub type CoatedConductorBxDF = LayeredBxDF<DielectricBxDF, ConductorBxDF, true>;

shared/src/bxdfs/measured.rs

@@ -0,0 +1,241 @@
+use crate::core::bxdf::{
+    BSDFSample, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode,
+};
+use crate::core::geometry::{
+    Point2f, Vector3f, VectorLike, abs_cos_theta, cos_theta, same_hemisphere, spherical_direction,
+    spherical_theta,
+};
+use crate::core::scattering::reflect;
+use crate::spectra::{SampledSpectrum, SampledWavelengths};
+use crate::utils::math::square;
+use crate::utils::ptr::{Ptr, Slice};
+use crate::utils::sampling::{PiecewiseLinear2D, cosine_hemisphere_pdf, sample_cosine_hemisphere};
+use crate::{Float, INV_PI, PI, PI_OVER_2};
+use core::any::Any;
+
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct MeasuredBxDFData {
+    pub wavelengths: Slice<Float>,
+    pub spectra: PiecewiseLinear2D<3>,
+    pub ndf: PiecewiseLinear2D<0>,
+    pub vndf: PiecewiseLinear2D<2>,
+    pub sigma: PiecewiseLinear2D<0>,
+    pub isotropic: bool,
+    pub luminance: PiecewiseLinear2D<2>,
+}
+
+#[repr(C)]
+#[derive(Debug, Copy, Clone)]
+pub struct MeasuredBxDF {
+    pub brdf: Ptr<MeasuredBxDFData>,
+    pub lambda: SampledWavelengths,
+}
+
+unsafe impl Send for MeasuredBxDF {}
+unsafe impl Sync for MeasuredBxDF {}
+
+impl MeasuredBxDF {
+    pub fn new(brdf: &MeasuredBxDFData, lambda: &SampledWavelengths) -> Self {
+        Self {
+            brdf: Ptr::from(brdf),
+            lambda: *lambda,
+        }
+    }
+
+    pub fn theta2u(theta: Float) -> Float {
+        (theta * (2. / PI)).sqrt()
+    }
+
+    pub fn phi2u(phi: Float) -> Float {
+        phi * 1. / (2. * PI) + 0.5
+    }
+
+    pub fn u2theta(u: Float) -> Float {
+        square(u) * PI_OVER_2
+    }
+
+    pub fn u2phi(u: Float) -> Float {
+        (2. * u - 1.) * PI
+    }
+}
+
+impl BxDFTrait for MeasuredBxDF {
+    fn flags(&self) -> BxDFFlags {
+        BxDFFlags::REFLECTION | BxDFFlags::GLOSSY
+    }
+
+    fn f(&self, wo: Vector3f, wi: Vector3f, _mode: TransportMode) -> SampledSpectrum {
+        if !same_hemisphere(wo, wi) {
+            return SampledSpectrum::new(0.);
+        }
+
+        let mut wo_curr = wo;
+        let mut wi_curr = wi;
+        if wo.z() < 0. {
+            wo_curr = -wo_curr;
+            wi_curr = -wi_curr;
+        }
+
+        // Get half direction vector
+        let wm_curr = wi_curr + wo_curr;
+        if wm_curr.norm_squared() == 0. {
+            return SampledSpectrum::new(0.);
+        }
+        let wm = wm_curr.normalize();
+
+        // Map vectors to unit square
+        let theta_o = spherical_theta(wo_curr);
+        let phi_o = wo_curr.y().atan2(wo_curr.x());
+        let theta_m = spherical_theta(wm);
+        let phi_m = wm.y().atan2(wm.x());
+        let u_wo = Point2f::new(MeasuredBxDF::theta2u(theta_o), MeasuredBxDF::phi2u(phi_o));
+        let u_wm_phi = if self.brdf.isotropic {
+            phi_m - phi_o
+        } else {
+            phi_m
+        };
+        let mut u_wm = Point2f::new(
+            MeasuredBxDF::theta2u(theta_m),
+            MeasuredBxDF::phi2u(u_wm_phi),
+        );
+        u_wm[1] -= u_wm[1].floor();
+
+        // Inverse parametrization
+        let ui = self.brdf.vndf.invert(u_wm, [phi_o, theta_o]);
+        let fr = SampledSpectrum::from_fn(|i| {
+            self.brdf
+                .spectra
+                .evaluate(ui.p, [phi_o, theta_o, self.lambda[i]])
+                .max(0.0)
+        });
+
+        fr * self.brdf.ndf.evaluate(u_wm, [])
+            / (4. * self.brdf.sigma.evaluate(u_wo, []) * cos_theta(wi))
+    }
+
+    fn pdf(&self, wo: Vector3f, wi: Vector3f, f_args: FArgs) -> Float {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return 0.;
+        }
+        if !same_hemisphere(wo, wi) {
+            return 0.;
+        }
+
+        let mut wo_curr = wo;
+        let mut wi_curr = wi;
+        if wo.z() < 0. {
+            wo_curr = -wo_curr;
+            wi_curr = -wi_curr;
+        }
+
+        let wm_curr = wi_curr + wo_curr;
+        if wm_curr.norm_squared() == 0. {
+            return 0.;
+        }
+        let wm = wm_curr.normalize();
+        let theta_o = spherical_theta(wo_curr);
+        let phi_o = wo_curr.y().atan2(wo_curr.x());
+        let theta_m = spherical_theta(wm);
+        let phi_m = wm.y().atan2(wm.x());
+
+        let u_wm_phi = if self.brdf.isotropic {
+            phi_m - phi_o
+        } else {
+            phi_m
+        };
+
+        let mut u_wm = Point2f::new(
+            MeasuredBxDF::theta2u(theta_m),
+            MeasuredBxDF::phi2u(u_wm_phi),
+        );
+        u_wm[1] = u_wm[1] - u_wm[1].floor();
+
+        let ui = self.brdf.vndf.invert(u_wm, [phi_o, theta_o]);
+        let sample = ui.p;
+        let vndf_pdf = ui.pdf;
+
+        let pdf = self.brdf.luminance.evaluate(sample, [phi_o, theta_o]);
+        let sin_theta_m = (square(wm.x()) + square(wm.y())).sqrt();
+        let jacobian = 4. * wm.dot(wo) * f32::max(2. * square(PI) * u_wm.x() * sin_theta_m, 1e-6);
+        vndf_pdf * pdf / jacobian
+    }
+
+    fn sample_f(&self, wo: Vector3f, _uc: Float, u: Point2f, f_args: FArgs) -> Option<BSDFSample> {
+        let reflection_flags =
+            BxDFReflTransFlags::from_bits_truncate(BxDFReflTransFlags::REFLECTION.bits());
+        if !f_args.sample_flags.contains(reflection_flags) {
+            return None;
+        }
+
+        let mut flip_w = false;
+        let mut wo_curr = wo;
+        if wo.z() <= 0. {
+            wo_curr = -wo_curr;
+            flip_w = true;
+        }
+
+        let theta_o = spherical_theta(wo_curr);
+        let phi_o = wo_curr.y().atan2(wo_curr.x());
+        // Warp sample using luminance distribution
+        let mut s = self.brdf.luminance.sample(u, [phi_o, theta_o]);
+        let u = s.p;
+        let lum_pdf = s.pdf;
+
+        // Sample visible normal distribution of measured BRDF
+        s = self.brdf.vndf.sample(u, [phi_o, theta_o]);
+        let u_wm = s.p;
+        let mut pdf = s.pdf;
+
+        // Map from microfacet normal to incident direction
+        let mut phi_m = MeasuredBxDF::u2phi(u_wm.y());
+        let theta_m = MeasuredBxDF::u2theta(u_wm.x());
+        if self.brdf.isotropic {
+            phi_m += phi_o;
+        }
+        let sin_theta_m = theta_m.sin();
+        let cos_theta_m = theta_m.cos();
+        let wm = spherical_direction(sin_theta_m, cos_theta_m, phi_m);
+        let mut wi = reflect(wo_curr, wm.into());
+        if wi.z() <= 0. {
+            return None;
+        }
+
+        // Interpolate spectral BRDF
+        let mut f = SampledSpectrum::from_fn(|i| {
+            self.brdf
+                .spectra
+                .evaluate(u, [phi_o, theta_o, self.lambda[i]])
+                .max(0.0)
+        });
+
+        let u_wo = Point2f::new(MeasuredBxDF::theta2u(theta_o), MeasuredBxDF::phi2u(phi_o));
+        f *= self.brdf.ndf.evaluate(u_wm, [])
+            / (4. * self.brdf.sigma.evaluate(u_wo, []) * abs_cos_theta(wi));
+        pdf /= 4. * wm.dot(wo_curr) * f32::max(2. * square(PI) * u_wm.x(), 1e-6);
+
+        if flip_w {
+            wi = -wi;
+        }
+
+        let bsdf = BSDFSample {
+            f,
+            wi,
+            pdf: pdf * lum_pdf,
+            flags: BxDFFlags::GLOSSY_REFLECTION,
+            ..Default::default()
+        };
+
+        Some(bsdf)
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn regularize(&mut self) {
+        return;
+    }
+}

shared/src/bxdfs/mod.rs

@@ -0,0 +1,13 @@
+pub mod complex;
+pub mod conductor;
+pub mod dielectric;
+pub mod diffuse;
+pub mod layered;
+pub mod measured;
+
+pub use complex::{EmptyBxDF, HairBxDF, NormalizedFresnelBxDF};
+pub use conductor::ConductorBxDF;
+pub use dielectric::{DielectricBxDF, ThinDielectricBxDF};
+pub use diffuse::{DiffuseBxDF, DiffuseTransmissionBxDF};
+pub use layered::{CoatedConductorBxDF, CoatedDiffuseBxDF};
+pub use measured::{MeasuredBxDF, MeasuredBxDFData};

shared/src/cameras/orthographic.rs

@@ -41,11 +41,12 @@ impl OrthographicCamera {
             -screen_window.p_max.y(),
             0.,
         ));
-        let film_ptr = base.film;
+
-        if film_ptr.is_null() {
+        let mut base_ortho = base;
+        let film = base.film;
+        if film.is_null() {
             panic!("Camera must have a film");
         }
-        let film = unsafe { &*film_ptr };
 
         let raster_from_ndc = Transform::scale(
             film.full_resolution().x() as Float,
@@ -59,7 +60,6 @@ impl OrthographicCamera {
         let camera_from_raster = screen_from_camera.inverse() * screen_from_raster;
         let dx_camera = camera_from_raster.apply_to_vector(Vector3f::new(1., 0., 0.));
         let dy_camera = camera_from_raster.apply_to_vector(Vector3f::new(0., 1., 0.));
-        let mut base_ortho = base;
         base_ortho.min_dir_differential_x = Vector3f::new(0., 0., 0.);
         base_ortho.min_dir_differential_y = Vector3f::new(0., 0., 0.);
         base_ortho.min_pos_differential_x = dx_camera;
@@ -96,7 +96,7 @@ impl CameraTrait for OrthographicCamera {
             p_camera,
             Vector3f::new(0., 0., 1.),
             Some(self.sample_time(sample.time)),
-            self.base().medium.clone(),
+            &*self.base().medium,
         );
         if self.lens_radius > 0. {
             let p_lens_vec =
@@ -127,11 +127,11 @@ impl CameraTrait for OrthographicCamera {
         let mut rd = RayDifferential::default();
         if self.lens_radius > 0.0 {
             let mut sample_x = sample;
-            sample_x.p_film.x += 1.0;
+            sample_x.p_film[0] += 1.0;
             let rx = self.generate_ray(sample_x, lambda)?;
 
             let mut sample_y = sample;
-            sample_y.p_film.y += 1.0;
+            sample_y.p_film[1] += 1.0;
             let ry = self.generate_ray(sample_y, lambda)?;
 
             rd.rx_origin = rx.ray.o;
@@ -155,7 +155,7 @@ impl CameraTrait for OrthographicCamera {
             rd.rx_direction = central_cam_ray.ray.d;
             rd.ry_direction = central_cam_ray.ray.d;
         }
-        central_cam_ray.ray.differential = Some(rd);
+        central_cam_ray.ray.differential = rd;
         Some(central_cam_ray)
     }
 }

shared/src/cameras/perspective.rs

@@ -1,5 +1,6 @@
 use crate::core::camera::{CameraBase, CameraRay, CameraTrait, CameraTransform};
 use crate::core::film::Film;
+use crate::core::filter::FilterTrait;
 use crate::core::geometry::{
     Bounds2f, Point2f, Point3f, Ray, RayDifferential, Vector2f, Vector3f, VectorLike,
 };
@@ -11,7 +12,7 @@ use crate::utils::sampling::sample_uniform_disk_concentric;
 use crate::utils::transform::Transform;
 
 #[repr(C)]
-#[derive(Debug, Copy)]
+#[derive(Debug, Clone, Copy)]
 pub struct PerspectiveCamera {
     pub base: CameraBase,
     pub screen_from_camera: Transform,
@@ -78,7 +79,7 @@ impl PerspectiveCamera {
     }
 }
 
-impl PerspectiveCamera {
+impl CameraTrait for PerspectiveCamera {
     fn base(&self) -> &CameraBase {
         &self.base
     }
@@ -97,7 +98,7 @@ impl PerspectiveCamera {
             Point3f::new(0., 0., 0.),
             p_vector.normalize(),
             Some(self.sample_time(sample.time)),
-            self.base().medium.clone(),
+            &*self.base().medium,
         );
         // Modify ray for depth of field
         if self.lens_radius > 0. {

shared/src/core/bsdf.rs

@@ -0,0 +1,125 @@
+use crate::Float;
+use crate::core::bxdf::{BSDFSample, BxDF, BxDFFlags, BxDFTrait, FArgs, TransportMode};
+use crate::core::geometry::{Frame, Normal3f, Point2f, Vector3f, VectorLike};
+use crate::spectra::SampledSpectrum;
+use crate::utils::Ptr;
+
+#[repr(C)]
+#[derive(Copy, Debug, Default)]
+pub struct BSDF {
+    bxdf: Ptr<BxDF>,
+    shading_frame: Frame,
+}
+
+impl BSDF {
+    pub fn new(ns: Normal3f, dpdus: Vector3f, bxdf: Ptr<BxDF>) -> Self {
+        Self {
+            bxdf,
+            shading_frame: Frame::new(dpdus.normalize(), Vector3f::from(ns)),
+        }
+    }
+
+    pub fn is_valid(&self) -> bool {
+        !self.bxdf.is_null()
+    }
+
+    pub fn flags(&self) -> BxDFFlags {
+        if self.bxdf.is_null() {
+            // Either this, or transmissive for seethrough
+            return BxDFFlags::empty();
+        }
+        self.bxdf.flags()
+    }
+
+    pub fn render_to_local(&self, v: Vector3f) -> Vector3f {
+        self.shading_frame.to_local(v)
+    }
+
+    pub fn local_to_render(&self, v: Vector3f) -> Vector3f {
+        self.shading_frame.from_local(v)
+    }
+
+    pub fn f(
+        &self,
+        wo_render: Vector3f,
+        wi_render: Vector3f,
+        mode: TransportMode,
+    ) -> Option<SampledSpectrum> {
+        if self.bxdf.is_null() {
+            return None;
+        }
+
+        let wi = self.render_to_local(wi_render);
+        let wo = self.render_to_local(wo_render);
+
+        if wo.z() == 0.0 || wi.z() == 0.0 {
+            return None;
+        }
+
+        Some(self.bxdf.f(wo, wi, mode))
+    }
+
+    pub fn sample_f(
+        &self,
+        wo_render: Vector3f,
+        u: Float,
+        u2: Point2f,
+        f_args: FArgs,
+    ) -> Option<BSDFSample> {
+        let bxdf = self.bxdf.as_ref()?;
+
+        let sampling_flags = BxDFFlags::from_bits_truncate(f_args.sample_flags.bits());
+        let wo = self.render_to_local(wo_render);
+        if wo.z() == 0.0 || !bxdf.flags().contains(sampling_flags) {
+            return None;
+        }
+
+        let mut sample = bxdf.sample_f(wo, u, u2, f_args)?;
+
+        if sample.pdf > 0.0 && sample.wi.z() != 0.0 {
+            sample.wi = self.local_to_render(sample.wi);
+            return Some(sample);
+        }
+
+        None
+    }
+
+    pub fn pdf(&self, wo_render: Vector3f, wi_render: Vector3f, f_args: FArgs) -> Float {
+        if self.bxdf.is_null() {
+            return 0.0;
+        }
+        let sample_flags = BxDFFlags::from_bits_truncate(f_args.sample_flags.bits());
+
+        let wo = self.render_to_local(wo_render);
+        let wi = self.render_to_local(wi_render);
+
+        if wo.z() == 0.0 || !self.bxdf.flags().contains(sample_flags) {
+            return 0.0;
+        }
+
+        self.bxdf.pdf(wo, wi, f_args)
+    }
+
+    pub fn rho_u(&self, u1: &[Point2f], uc: &[Float], u2: &[Point2f]) -> SampledSpectrum {
+        if self.bxdf.is_null() {
+            return SampledSpectrum::default();
+        }
+
+        self.bxdf.rho_u(u1, uc, u2)
+    }
+
+    pub fn rho_wo(&self, wo_render: Vector3f, uc: &[Float], u: &[Point2f]) -> SampledSpectrum {
+        if self.bxdf.is_null() {
+            return SampledSpectrum::default();
+        }
+
+        let wo = self.render_to_local(wo_render);
+        self.bxdf.rho_wo(wo, uc, u)
+    }
+
+    pub fn regularize(&mut self) {
+        if !self.bxdf.is_null() {
+            unsafe { self.bxdf.as_mut().regularize() }
+        }
+    }
+}

shared/src/core/bssrdf.rs

@@ -1,13 +1,13 @@
 use crate::core::bxdf::{BSDF, NormalizedFresnelBxDF};
 use crate::core::geometry::{Frame, Normal3f, Point2f, Point3f, Point3fi, Vector3f};
 use crate::core::interaction::{InteractionBase, ShadingGeom, SurfaceInteraction};
-use crate::core::pbrt::{Float, PI};
 use crate::core::shape::Shape;
 use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum};
-use crate::utils::RelPtr;
+use crate::utils::ArenaPtr;
 use crate::utils::math::{catmull_rom_weights, square};
 use crate::utils::sampling::sample_catmull_rom_2d;
-use crate::utils::{Ptr, Slice};
+use crate::utils::{Ptr, ptr::Slice};
+use crate::{Float, PI};
 use enum_dispatch::enum_dispatch;
 use std::sync::Arc;
 
@@ -21,13 +21,13 @@ pub struct BSSRDFSample {
 
 #[derive(Clone, Debug)]
 pub struct SubsurfaceInteraction {
-    pi: Point3fi,
+    pub pi: Point3fi,
-    n: Normal3f,
+    pub n: Normal3f,
-    ns: Normal3f,
+    pub ns: Normal3f,
-    dpdu: Vector3f,
+    pub dpdu: Vector3f,
-    dpdv: Vector3f,
+    pub dpdv: Vector3f,
-    dpdus: Vector3f,
+    pub dpdus: Vector3f,
-    dpdvs: Vector3f,
+    pub dpdvs: Vector3f,
 }
 
 impl SubsurfaceInteraction {
@@ -94,18 +94,41 @@ impl From<&SubsurfaceInteraction> for SurfaceInteraction {
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct BSSRDFTable {
-    pub rho_samples: Slice<Float>,
+    pub n_rho_samples: u32,
-    pub radius_samples: *const Float,
+    pub n_radius_samples: u32,
-    pub profile: *const Float,
+    pub rho_samples: Ptr<Float>,
-    pub rho_eff: *const Float,
+    pub radius_samples: Ptr<Float>,
-    pub profile_cdf: *const Float,
+    pub profile: Ptr<Float>,
+    pub rho_eff: Ptr<Float>,
+    pub profile_cdf: Ptr<Float>,
 }
 
 impl BSSRDFTable {
-    pub fn eval_profile(&self, rho_index: usize, radius_index: usize) -> Float {
+    pub fn get_rho(&self) -> &[Float] {
-        assert!(rho_index < self.rho_samples.len());
+        unsafe { core::slice::from_raw_parts(self.rho_samples.0, self.n_rho_samples as usize) }
-        assert!(radius_index < self.radius_samples.len());
+    }
-        self.profile[rho_index * self.radius_samples.len() + radius_index]
+
+    pub fn get_radius(&self) -> &[Float] {
+        unsafe {
+            core::slice::from_raw_parts(self.radius_samples.0, self.n_radius_samples as usize)
+        }
+    }
+
+    pub fn get_profile(&self) -> &[Float] {
+        let n_profile = (self.n_rho_samples * self.n_radius_samples) as usize;
+        unsafe { core::slice::from_raw_parts(self.profile.0, n_profile) }
+    }
+
+    pub fn get_cdf(&self) -> &[Float] {
+        let n_profile = (self.n_rho_samples * self.n_radius_samples) as usize;
+        unsafe { core::slice::from_raw_parts(self.profile_cdf.0, n_profile) }
+    }
+
+    pub fn eval_profile(&self, rho_index: u32, radius_index: u32) -> Float {
+        debug_assert!(rho_index < self.n_rho_samples);
+        debug_assert!(radius_index < self.n_radius_samples);
+        let idx = (rho_index * self.n_radius_samples + radius_index) as usize;
+        unsafe { *self.profile.0.add(idx) }
     }
 }
 
@@ -140,9 +163,9 @@ pub struct TabulatedBSSRDF {
     wo: Vector3f,
     ns: Normal3f,
     eta: Float,
-    sigma_t: RelPtr<SampledSpectrum>,
+    sigma_t: SampledSpectrum,
-    rho: RelPtr<SampledSpectrum>,
+    rho: SampledSpectrum,
-    table: *const BSSRDFTable,
+    table: Ptr<BSSRDFTable>,
 }
 
 impl TabulatedBSSRDF {
@@ -153,7 +176,7 @@ impl TabulatedBSSRDF {
         eta: Float,
         sigma_a: &SampledSpectrum,
         sigma_s: &SampledSpectrum,
-        table: *const BSSRDFTable,
+        table: &BSSRDFTable,
     ) -> Self {
         let sigma_t = *sigma_a + *sigma_s;
         let rho = SampledSpectrum::safe_div(sigma_s, &sigma_t);
@@ -162,9 +185,9 @@ impl TabulatedBSSRDF {
             wo,
             ns,
             eta,
-            table,
             sigma_t,
             rho,
+            table: Ptr::from(table),
         }
     }
 
@@ -174,16 +197,17 @@ impl TabulatedBSSRDF {
 
     pub fn sr(&self, r: Float) -> SampledSpectrum {
         let mut sr_spectrum = SampledSpectrum::new(0.);
+        let rho_samples = self.table.get_rho();
+        let radius_samples = self.table.get_radius();
         for i in 0..N_SPECTRUM_SAMPLES {
             let r_optical = r * self.sigma_t[i];
-            let (rho_offset, rho_weights) =
+            let (rho_offset, rho_weights) = match catmull_rom_weights(rho_samples, self.rho[i]) {
-                match catmull_rom_weights(&self.table.rho_samples, self.rho[i]) {
                 Some(res) => res,
                 None => continue,
             };
 
             let (radius_offset, radius_weights) =
-                match catmull_rom_weights(&self.table.radius_samples, r_optical) {
+                match catmull_rom_weights(radius_samples, r_optical) {
                     Some(res) => res,
                     None => continue,
                 };
@@ -193,7 +217,10 @@ impl TabulatedBSSRDF {
                 for (k, radius_weight) in radius_weights.iter().enumerate() {
                     let weight = rho_weight * radius_weight;
                     if weight != 0. {
-                        sr += weight * self.table.eval_profile(rho_offset + j, radius_offset + k);
+                        sr += weight
+                            * self
+                                .table
+                                .eval_profile(rho_offset + j as u32, radius_offset + k as u32);
                     }
                 }
             }
@@ -203,7 +230,7 @@ impl TabulatedBSSRDF {
             sr_spectrum[i] = sr;
         }
 
-        sr_spectrum *= self.sigma_t * self.sigma_t;
+        sr_spectrum *= square(self.sigma_t);
         SampledSpectrum::clamp_zero(&sr_spectrum)
     }
 
@@ -211,29 +238,30 @@ impl TabulatedBSSRDF {
         if self.sigma_t[0] == 0. {
             return None;
         }
-        let (ret, _, _) = sample_catmull_rom_2d(
+
-            &self.table.rho_samples,
+        let rho_samples = self.table.get_rho();
-            &self.table.radius_samples,
+        let radius_samples = self.table.get_radius();
-            &self.table.profile,
+        let profile = self.table.get_profile();
-            &self.table.profile_cdf,
+        let cdf = self.table.get_cdf();
-            self.rho[0],
+
-            u,
+        let (ret, _, _) =
-        );
+            sample_catmull_rom_2d(rho_samples, radius_samples, profile, cdf, self.rho[0], u);
         Some(ret / self.sigma_t[0])
     }
 
     pub fn pdf_sr(&self, r: Float) -> SampledSpectrum {
         let mut pdf = SampledSpectrum::new(0.);
+        let rhoeff_samples = self.table.get_rho();
+        let radius_samples = self.table.get_radius();
         for i in 0..N_SPECTRUM_SAMPLES {
             let r_optical = r * self.sigma_t[i];
-            let (rho_offset, rho_weights) =
+            let (rho_offset, rho_weights) = match catmull_rom_weights(rhoeff_samples, self.rho[i]) {
-                match catmull_rom_weights(&self.table.rho_samples, self.rho[i]) {
                 Some(res) => res,
                 None => continue,
             };
 
             let (radius_offset, radius_weights) =
-                match catmull_rom_weights(&self.table.radius_samples, r_optical) {
+                match catmull_rom_weights(radius_samples, r_optical) {
                     Some(res) => res,
                     None => continue,
                 };
@@ -243,12 +271,14 @@ impl TabulatedBSSRDF {
             for (j, rho_weight) in rho_weights.iter().enumerate() {
                 if *rho_weight != 0. {
                     // Update _rhoEff_ and _sr_ for wavelength
-                    rho_eff += self.table.rho_eff[rho_offset + j] * rho_weight;
+                    rho_eff += rhoeff_samples[rho_offset as usize + j] * rho_weight;
 
                     // Fix: Use .iter().enumerate() for 'k'
                     for (k, radius_weight) in radius_weights.iter().enumerate() {
                         if *radius_weight != 0. {
-                            sr += self.table.eval_profile(rho_offset + j, radius_offset + k)
+                            sr += self
+                                .table
+                                .eval_profile(rho_offset + j as u32, radius_offset + k as u32)
                                 * rho_weight
                                 * radius_weight;
                         }

shared/src/core/camera.rs

@@ -10,6 +10,7 @@ use crate::core::pbrt::Float;
 use crate::core::sampler::CameraSample;
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::math::lerp;
+use crate::utils::ptr::Ptr;
 use crate::utils::transform::{AnimatedTransform, Transform};
 
 use enum_dispatch::enum_dispatch;
@@ -108,12 +109,12 @@ pub struct CameraBase {
     pub camera_transform: CameraTransform,
     pub shutter_open: Float,
     pub shutter_close: Float,
-    pub film: *const Film,
-    pub medium: *const Medium,
     pub min_pos_differential_x: Vector3f,
     pub min_pos_differential_y: Vector3f,
     pub min_dir_differential_x: Vector3f,
     pub min_dir_differential_y: Vector3f,
+    pub film: Ptr<Film>,
+    pub medium: Ptr<Medium>,
 }
 
 #[enum_dispatch(CameraTrait)]
@@ -140,7 +141,7 @@ pub trait CameraTrait {
                 );
             }
         }
-        unsafe { &*self.base().film }
+        &*self.base().film
     }
 
     fn sample_time(&self, u: Float) -> Float {
@@ -199,7 +200,7 @@ pub trait CameraTrait {
                 }
 
                 if rx_found && ry_found {
-                    central_cam_ray.ray.differential = Some(rd);
+                    central_cam_ray.ray.differential = rd;
                 }
 
                 Some(central_cam_ray)
@@ -231,13 +232,13 @@ pub trait CameraTrait {
             Point3f::new(0., 0., 0.) + self.base().min_pos_differential_x,
             Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_x,
             None,
-            None,
+            &Ptr::default(),
         );
         let y_ray = Ray::new(
             Point3f::new(0., 0., 0.) + self.base().min_pos_differential_y,
             Vector3f::new(0., 0., 1.) + self.base().min_dir_differential_y,
             None,
-            None,
+            &Ptr::default(),
         );
         let n_down = Vector3f::from(n_down_z);
         let tx = -(n_down.dot(y_ray.o.into())) / n_down.dot(x_ray.d);

shared/src/core/film.rs

@@ -17,6 +17,7 @@ use crate::utils::AtomicFloat;
 use crate::utils::containers::Array2D;
 use crate::utils::math::linear_least_squares;
 use crate::utils::math::{SquareMatrix, wrap_equal_area_square};
+use crate::utils::ptr::Ptr;
 use crate::utils::sampling::VarianceEstimator;
 use crate::utils::transform::AnimatedTransform;
 
@@ -39,45 +40,45 @@ pub struct RGBPixel {
     rgb_splat: [AtomicFloat; 3],
 }
 
-#[cfg(not(target_os = "cuda"))]
+// #[cfg(not(target_os = "cuda"))]
-impl RGBFilm {
+// impl RGBFilm {
-    pub fn new(
+//     pub fn new(
-        base: FilmBase,
+//         base: FilmBase,
-        colorspace: &RGBColorSpace,
+//         colorspace: &RGBColorSpace,
-        max_component_value: Float,
+//         max_component_value: Float,
-        write_fp16: bool,
+//         write_fp16: bool,
-    ) -> Self {
+//     ) -> Self {
-        let sensor_ptr = base.sensor;
+//         let sensor_ptr = base.sensor;
-        if sensor_ptr.is_null() {
+//         if sensor_ptr.is_null() {
-            panic!("Film must have a sensor");
+//             panic!("Film must have a sensor");
-        }
+//         }
-        let sensor = unsafe { &*sensor_ptr };
+//         let sensor = unsafe { &*sensor_ptr };
-        let filter_integral = base.filter.integral();
+//         let filter_integral = base.filter.integral();
-        let sensor_matrix = sensor.xyz_from_sensor_rgb;
+//         let sensor_matrix = sensor.xyz_from_sensor_rgb;
-        let output_rgbf_from_sensor_rgb = colorspace.rgb_from_xyz * sensor_matrix;
+//         let output_rgbf_from_sensor_rgb = colorspace.rgb_from_xyz * sensor_matrix;
-
+//
-        let width = base.pixel_bounds.p_max.x() - base.pixel_bounds.p_min.x();
+//         let width = base.pixel_bounds.p_max.x() - base.pixel_bounds.p_min.x();
-        let height = base.pixel_bounds.p_max.y() - base.pixel_bounds.p_min.y();
+//         let height = base.pixel_bounds.p_max.y() - base.pixel_bounds.p_min.y();
-        let count = (width * height) as usize;
+//         let count = (width * height) as usize;
-
+//
-        let mut pixel_vec = Vec::with_capacity(count);
+//         let mut pixel_vec = Vec::with_capacity(count);
-        for _ in 0..count {
+//         for _ in 0..count {
-            pixel_vec.push(RGBPixel::default());
+//             pixel_vec.push(RGBPixel::default());
-        }
+//         }
-
+//
-        let pixels_array = Array2D::new(base.pixel_bounds);
+//         let pixels_array = Array2D::(base.pixel_bounds);
-
+//
-        Self {
+//         Self {
-            base,
+//             base,
-            max_component_value,
+//             max_component_value,
-            write_fp16,
+//             write_fp16,
-            filter_integral,
+//             filter_integral,
-            output_rgbf_from_sensor_rgb,
+//             output_rgbf_from_sensor_rgb,
-            pixels: std::sync::Arc::new(pixels_array),
+//             pixels: std::sync::Arc::new(pixels_array),
-        }
+//         }
-    }
+//     }
-}
+// }
-
+//
 impl RGBFilm {
     pub fn base(&self) -> &FilmBase {
         &self.base
@@ -90,13 +91,13 @@ impl RGBFilm {
     pub fn get_sensor(&self) -> &PixelSensor {
         #[cfg(not(target_os = "cuda"))]
         {
-            if self.film.sensor.is_null() {
+            if self.base.sensor.is_null() {
                 panic!(
                     "FilmBase error: PixelSensor pointer is null. This should have been checked during construction."
                 );
             }
         }
-        unsafe { &*self.sensor }
+        unsafe { &*self.base.sensor }
     }
 
     pub fn add_sample(
@@ -130,18 +131,19 @@ impl RGBFilm {
         }
 
         let p_discrete = p + Vector2f::new(0.5, 0.5);
-        let radius = self.get_filter().radius();
+        let radius = self.base.filter.radius();
 
         let splat_bounds = Bounds2i::from_points(
             (p_discrete - radius).floor(),
             (p_discrete + radius).floor() + Vector2i::new(1, 1),
         );
 
-        let splat_intersect = splat_bounds.union(self.pixel_bounds());
+        let splat_intersect = splat_bounds.union(self.base().pixel_bounds);
         for pi in &splat_intersect {
             let pi_f: Point2f = (*pi).into();
             let wt = self
-                .get_filter()
+                .base()
+                .filter
                 .evaluate((p - pi_f - Vector2f::new(0.5, 0.5)).into());
             if wt != 0. {
                 let pixel = &self.pixels[*pi];
@@ -179,8 +181,8 @@ impl RGBFilm {
     }
 
     pub fn to_output_rgb(&self, l: SampledSpectrum, lambda: &SampledWavelengths) -> RGB {
-        let sensor = unsafe { self.get_sensor() };
+        let sensor = self.get_sensor();
-        let mut sensor_rgb = sensor.to_sensor_rgb(l, lambda);
+        let sensor_rgb = sensor.to_sensor_rgb(l, lambda);
         self.output_rgbf_from_sensor_rgb.mul_rgb(sensor_rgb)
     }
 
@@ -268,13 +270,13 @@ impl GBufferFilm {
     pub fn get_sensor(&self) -> &PixelSensor {
         #[cfg(not(target_os = "cuda"))]
         {
-            if self.sensor.is_null() {
+            if self.base.sensor.is_null() {
                 panic!(
                     "FilmBase error: PixelSensor pointer is null. This should have been checked during construction."
                 );
             }
         }
-        unsafe { &*self.sensor }
+        unsafe { &*self.base.sensor }
     }
 
     pub fn add_splat(&mut self, p: Point2f, l: SampledSpectrum, lambda: &SampledWavelengths) {
@@ -286,18 +288,19 @@ impl GBufferFilm {
         }
 
         let p_discrete = p + Vector2f::new(0.5, 0.5);
-        let radius = self.get_filter().radius();
+        let radius = self.base().filter.radius();
 
         let splat_bounds = Bounds2i::from_points(
             (p_discrete - radius).floor(),
             (p_discrete + radius).floor() + Vector2i::new(1, 1),
         );
 
-        let splat_intersect = splat_bounds.union(self.pixel_bounds());
+        let splat_intersect = splat_bounds.union(self.base.pixel_bounds);
         for pi in &splat_intersect {
             let pi_f: Point2f = (*pi).into();
             let wt = self
-                .get_filter()
+                .base
+                .filter
                 .evaluate((p - pi_f - Vector2f::new(0.5, 0.5)).into());
             if wt != 0. {
                 let pixel = &self.pixels[*pi];
@@ -309,7 +312,7 @@ impl GBufferFilm {
     }
 
     pub fn to_output_rgb(&self, l: SampledSpectrum, lambda: &SampledWavelengths) -> RGB {
-        let sensor = unsafe { self.get_sensor() };
+        let sensor = self.get_sensor();
         let sensor_rgb = sensor.to_sensor_rgb(l, lambda);
         self.output_rgbf_from_sensor_rgb.mul_rgb(sensor_rgb)
     }
@@ -406,7 +409,7 @@ impl PixelSensor {
         g: Spectrum,
         b: Spectrum,
         output_colorspace: RGBColorSpace,
-        sensor_illum: *const Spectrum,
+        sensor_illum: &Spectrum,
         imaging_ratio: Float,
         spectra: *const StandardSpectra,
         swatches: &[Spectrum; 24],
@@ -414,7 +417,7 @@ impl PixelSensor {
         // As seen in usages of this constructos, sensor_illum can be null
         // Going with the colorspace's own illuminant, but this might not be the right choice
         // TODO: Test this
-        let illum: &Spectrum = match &sensor_illum {
+        let illum: &Spectrum = match sensor_illum {
             Some(arc_illum) => &**arc_illum,
             None => &output_colorspace.illuminant,
         };
@@ -469,7 +472,7 @@ impl PixelSensor {
 
     pub fn new_with_white_balance(
         output_colorspace: &RGBColorSpace,
-        sensor_illum: Option<std::sync::Arc<Spectrum>>,
+        sensor_illum: Ptr<Spectrum>,
         imaging_ratio: Float,
         spectra: *const StandardSpectra,
     ) -> Self {

shared/src/core/interaction.rs

@@ -511,14 +511,14 @@ impl SurfaceInteraction {
                     || Vector3f::from(new_diff_rx_origin).norm_squared() > threshold
                     || Vector3f::from(new_diff_ry_origin).norm_squared() > threshold
                 {
-                    rd.differential = None;
+                    rd.differential = RayDifferential::default();
                 } else {
-                    rd.differential = Some(RayDifferential {
+                    rd.differential = 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,
-                    });
+                    };
                 }
             }
         }
@@ -537,13 +537,12 @@ impl InteractionTrait for SurfaceInteraction {
     }
 
     fn get_medium(&self, w: Vector3f) -> Ptr<Medium> {
-        self.common.medium_interface.as_ref().and_then(|interface| {
+        let interface = self.common.medium_interface;
         if self.n().dot(w.into()) > 0.0 {
             interface.outside
         } else {
             interface.inside
         }
-        })
     }
 
     fn is_surface_interaction(&self) -> bool {
@@ -573,7 +572,6 @@ impl SurfaceInteraction {
 
         Self {
             common: InteractionBase::new_surface_geom(pi, n, uv, wo, time),
-            uv,
             dpdu,
             dpdv,
             dndu,
@@ -637,7 +635,6 @@ impl SurfaceInteraction {
     pub fn new_simple(pi: Point3fi, n: Normal3f, uv: Point2f) -> Self {
         Self {
             common: InteractionBase::new_surface_geom(pi, n, uv, Vector3f::zero(), 0.),
-            uv,
             ..Default::default()
         }
     }
@@ -646,9 +643,9 @@ impl SurfaceInteraction {
         Self {
             common: InteractionBase {
                 pi,
+                uv,
                 ..Default::default()
             },
-            uv,
             ..Default::default()
         }
     }
@@ -656,18 +653,18 @@ impl SurfaceInteraction {
     #[cfg(not(target_os = "cuda"))]
     pub fn set_intersection_properties(
         &mut self,
-        mtl: *const Material,
+        mtl: &Material,
-        area: *const Light,
+        area: &Light,
+        ray_medium: &Medium,
         prim_medium_interface: MediumInterface,
-        ray_medium: *const Medium,
     ) {
-        self.material = mtl;
+        self.material = Ptr::from(mtl);
-        self.area_light = area;
+        self.area_light = Ptr::from(area);
 
         if prim_medium_interface.is_medium_transition() {
             self.common.medium_interface = prim_medium_interface;
         } else {
-            self.common.medium = ray_medium;
+            self.common.medium = Ptr::from(ray_medium);
         }
     }
 }

shared/src/core/material.rs

@@ -3,10 +3,12 @@ use enum_dispatch::enum_dispatch;
 use std::ops::Deref;
 
 use crate::Float;
-use crate::core::bssrdf::BSSRDF;
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::BSSRDF;
+use crate::core::bxdf::BxDF;
 use crate::core::geometry::{Frame, Normal3f, Point2f, Point3f, Vector2f, Vector3f, VectorLike};
 use crate::core::image::{Image, WrapMode, WrapMode2D};
 use crate::core::interaction::{Interaction, InteractionTrait, ShadingGeom, SurfaceInteraction};
@@ -17,7 +19,7 @@ use crate::core::texture::{
 };
 use crate::materials::*;
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::hash::hash_float;
 use crate::utils::math::clamp;
 
@@ -63,14 +65,14 @@ pub struct NormalBumpEvalContext {
     pub dudy: Float,
     pub dvdx: Float,
     pub dvdy: Float,
-    pub face_index: usize,
+    pub face_index: u32,
 }
 
 impl From<&SurfaceInteraction> for NormalBumpEvalContext {
     fn from(si: &SurfaceInteraction) -> Self {
         Self {
             p: si.p(),
-            uv: si.uv,
+            uv: si.common.uv,
             n: si.n(),
             shading: si.shading.clone(),
             dudx: si.dudx,
@@ -123,7 +125,7 @@ pub fn bump_map<T: TextureEvaluator>(
     displacement: &GPUFloatTexture,
     ctx: &NormalBumpEvalContext,
 ) -> (Vector3f, Vector3f) {
-    debug_assert!(tex_eval.can_evaluate(&[displacement], &[]));
+    debug_assert!(tex_eval.can_evaluate(&[Ptr::from(displacement)], &[]));
     let mut du = 0.5 * (ctx.dudx.abs() + ctx.dudy.abs());
     if du == 0.0 {
         du = 0.0005;
@@ -171,8 +173,8 @@ pub trait MaterialTrait {
     ) -> Option<BSSRDF>;
 
     fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool;
-    fn get_normal_map(&self) -> *const Image;
+    fn get_normal_map(&self) -> Option<&Image>;
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture>;
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture>;
     fn has_subsurface_scattering(&self) -> bool;
 }
 

shared/src/core/mod.rs

@@ -1,4 +1,4 @@
-pub mod aggregates;
+pub mod bsdf;
 pub mod bssrdf;
 pub mod bxdf;
 pub mod camera;

shared/src/core/pbrt.rs

@@ -101,35 +101,27 @@ pub const PI_OVER_4: Float = 0.785_398_163_397_448_309_61;
 pub const SQRT_2: Float = 1.414_213_562_373_095_048_80;
 
 #[inline]
-pub fn find_interval<T, P>(sz: T, pred: P) -> T
+pub fn find_interval<F>(sz: u32, pred: F) -> u32
 where
-    T: PrimInt,
+    F: Fn(u32) -> bool,
-    P: Fn(T) -> bool,
 {
-    let zero = T::zero();
+    let mut first = 0;
-    let one = T::one();
+    let mut len = sz;
-    let two = one + one;
 
-    if sz <= two {
+    while len > 0 {
-        return zero;
+        let half = len >> 1;
-    }
+        let middle = first + half;
 
-    let mut low = one;
+        if pred(middle) {
-    let mut high = sz - one;
+            first = middle + 1;
-
+            len -= half + 1;
-    while low < high {
-        // mid = low + (high - low) / 2
-        let mid = low + (high - low) / two;
-        if pred(mid) {
-            low = mid + one;
         } else {
-            high = mid;
+            len = half;
         }
     }
 
-    let result = low - one;
+    let ret = (first as i32 - 1).max(0) as u32;
-
+    ret.min(sz.saturating_sub(2))
-    num_traits::clamp(result, zero, sz - two)
 }
 
 #[inline]

shared/src/core/primitive.rs

@@ -7,7 +7,7 @@ use crate::core::medium::{Medium, MediumInterface};
 use crate::core::pbrt::Float;
 use crate::core::shape::{Shape, ShapeIntersection, ShapeTrait};
 use crate::core::texture::{GPUFloatTexture, TextureEvalContext};
-use crate::utils::RelPtr;
+use crate::utils::ArenaPtr;
 use crate::utils::hash::hash_float;
 use crate::utils::transform::{AnimatedTransform, Transform};
 
@@ -88,13 +88,13 @@ impl PrimitiveTrait for GeometricPrimitive {
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
 pub struct SimplePrimitive {
-    shape: RelPtr<Shape>,
+    shape: ArenaPtr<Shape>,
-    material: RelPtr<Material>,
+    material: ArenaPtr<Material>,
 }
 
 #[derive(Debug, Clone)]
 pub struct TransformedPrimitive {
-    pub primitive: RelPtr<Primitive>,
+    pub primitive: ArenaPtr<Primitive>,
     pub render_from_primitive: Transform,
 }
 
@@ -129,7 +129,7 @@ impl PrimitiveTrait for TransformedPrimitive {
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
 pub struct AnimatedPrimitive {
-    primitive: RelPtr<Primitive>,
+    primitive: ArenaPtr<Primitive>,
     render_from_primitive: AnimatedTransform,
 }
 
@@ -164,7 +164,7 @@ impl PrimitiveTrait for AnimatedPrimitive {
 #[derive(Debug, Clone, Copy)]
 pub struct BVHAggregatePrimitive {
     max_prims_in_node: u32,
-    primitives: *const RelPtr<Primitive>,
+    primitives: *const ArenaPtr<Primitive>,
     nodes: *const LinearBVHNode,
 }
 

shared/src/core/sampler.rs

@@ -2,6 +2,7 @@ use crate::core::filter::FilterTrait;
 use crate::core::geometry::{Bounds2f, Point2f, Point2i, Vector2f};
 use crate::core::options::{PBRTOptions, get_options};
 use crate::core::pbrt::{Float, ONE_MINUS_EPSILON, PI, PI_OVER_2, PI_OVER_4, find_interval};
+use crate::utils::Ptr;
 use crate::utils::containers::Array2D;
 use crate::utils::math::{
     BinaryPermuteScrambler, DigitPermutation, FastOwenScrambler, NoRandomizer, OwenScrambler,
@@ -102,7 +103,7 @@ pub struct HaltonSampler {
     mult_inverse: [u64; 2],
     halton_index: u64,
     dim: u32,
-    digit_permutations: *const DigitPermutation,
+    digit_permutations: Ptr<DigitPermutation>,
 }
 
 impl HaltonSampler {
@@ -160,7 +161,7 @@ impl HaltonSampler {
             scrambled_radical_inverse(
                 dimension,
                 self.halton_index,
-                &self.digit_permutations[dimension],
+                &self.digit_permutations[dimension as usize],
             )
         } else {
             owen_scrambled_radical_inverse(

shared/src/core/scattering.rs

@@ -171,3 +171,40 @@ pub fn fr_complex_from_spectrum(
     }
     result
 }
+
+pub fn fresnel_moment1(eta: Float) -> Float {
+    let eta2 = eta * eta;
+    let eta3 = eta2 * eta;
+    let eta4 = eta3 * eta;
+    let eta5 = eta4 * eta;
+    if eta < 1. {
+        return 0.45966 - 1.73965 * eta + 3.37668 * eta2 - 3.904945 * eta3 + 2.49277 * eta4
+            - 0.68441 * eta5;
+    } else {
+        return -4.61686 + 11.1136 * eta - 10.4646 * eta2 + 5.11455 * eta3 - 1.27198 * eta4
+            + 0.12746 * eta5;
+    }
+}
+
+pub fn fresnel_moment2(eta: Float) -> Float {
+    let eta2 = eta * eta;
+    let eta3 = eta2 * eta;
+    let eta4 = eta3 * eta;
+    let eta5 = eta4 * eta;
+
+    if eta < 1. {
+        return 0.27614 - 0.87350 * eta + 1.12077 * eta2 - 0.65095 * eta3
+            + 0.07883 * eta4
+            + 0.04860 * eta5;
+    } else {
+        let r_eta = 1. / eta;
+        let r_eta2 = r_eta * r_eta;
+        let r_eta3 = r_eta2 * r_eta;
+
+        return -547.033 + 45.3087 * r_eta3 - 218.725 * r_eta2 + 458.843 * r_eta + 404.557 * eta
+            - 189.519 * eta2
+            + 54.9327 * eta3
+            - 9.00603 * eta4
+            + 0.63942 * eta5;
+    }
+}

shared/src/core/shape.rs

@@ -8,10 +8,10 @@ use crate::core::interaction::{
 use crate::core::light::Light;
 use crate::core::material::Material;
 use crate::core::medium::{Medium, MediumInterface};
-use crate::core::pbrt::{Float, PI};
 use crate::shapes::*;
-use crate::utils::Transform;
 use crate::utils::math::{next_float_down, next_float_up};
+use crate::utils::{Ptr, Transform};
+use crate::{Float, PI};
 use enum_dispatch::enum_dispatch;
 
 // Define Intersection objects. This only varies for
@@ -37,14 +37,13 @@ impl ShapeIntersection {
 
     pub fn set_intersection_properties(
         &mut self,
-        mtl: *const Material,
+        mtl: &Material,
-        area: *const Light,
+        area: &Light,
-        prim_medium_interface: *const MediumInterface,
+        prim_medium_interface: MediumInterface,
-        ray_medium: *const Medium,
+        ray_medium: &Medium,
     ) {
-        let ray_medium = unsafe { *prim_medium_interface };
         self.intr
-            .set_intersection_properties(mtl, area, prim_medium_interface, ray_medium);
+            .set_intersection_properties(mtl, area, ray_medium, prim_medium_interface);
     }
 }
 
@@ -119,12 +118,7 @@ impl ShapeSampleContext {
     }
 
     pub fn spawn_ray(&self, w: Vector3f) -> Ray {
-        Ray::new(
+        Ray::new(self.offset_ray_origin(w), w, Some(self.time), &Ptr::null())
-            self.offset_ray_origin(w),
-            w,
-            Some(self.time),
-            core::ptr::null(),
-        )
     }
 }
 

shared/src/core/spectrum.rs

@@ -43,9 +43,9 @@ impl Spectrum {
     }
 
     pub fn to_xyz(&self, std: &StandardSpectra) -> XYZ {
-        let x = self.inner_product(&std.x());
+        let x = self.inner_product(&Spectrum::Dense(std.x));
-        let y = self.inner_product(&std.y());
+        let y = self.inner_product(&Spectrum::Dense(std.y));
-        let z = self.inner_product(&std.z());
+        let z = self.inner_product(&Spectrum::Dense(std.z));
 
         XYZ::new(x, y, z) / CIE_Y_INTEGRAL
     }

shared/src/core/texture.rs

@@ -9,7 +9,7 @@ use crate::spectra::{
     SampledWavelengths,
 };
 use crate::textures::*;
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::Transform;
 use crate::utils::math::square;
 use crate::{Float, INV_2_PI, INV_PI, PI};
@@ -428,8 +428,8 @@ pub trait TextureEvaluator: Send + Sync {
 
     fn can_evaluate(
         &self,
-        _ftex: &[RelPtr<GPUFloatTexture>],
+        _ftex: &[Ptr<GPUFloatTexture>],
-        _stex: &[RelPtr<GPUSpectrumTexture>],
+        _stex: &[Ptr<GPUSpectrumTexture>],
     ) -> bool;
 }
 
@@ -453,8 +453,8 @@ impl TextureEvaluator for UniversalTextureEvaluator {
 
     fn can_evaluate(
         &self,
-        _float_textures: &[RelPtr<GPUFloatTexture>],
+        _float_textures: &[Ptr<GPUFloatTexture>],
-        _spectrum_textures: &[RelPtr<GPUSpectrumTexture>],
+        _spectrum_textures: &[Ptr<GPUSpectrumTexture>],
     ) -> bool {
         true
     }

shared/src/lib.rs

@@ -2,6 +2,7 @@
 #![feature(float_erf)]
 #![feature(f16)]
 
+pub mod bxdfs;
 pub mod cameras;
 pub mod core;
 pub mod data;

shared/src/lights/diffuse.rs

@@ -36,7 +36,6 @@ pub struct DiffuseAreaLight {
 unsafe impl Send for DiffuseAreaLight {}
 unsafe impl Sync for DiffuseAreaLight {}
 
-#[cfg(not(target_os = "cuda"))]
 impl DiffuseAreaLight {
     fn l_base(&self, n: Normal3f, wo: Vector3f, lambda: &SampledWavelengths) -> SampledSpectrum {
         if !self.two_sided && n.dot(wo.into()) <= 0.0 {
@@ -46,11 +45,11 @@ impl DiffuseAreaLight {
     }
 
     fn alpha_masked(&self, intr: &Interaction) -> bool {
-        let Some(alpha_tex) = &self.alpha else {
+        if self.alpha.is_null() {
             return false;
         };
         let ctx = TextureEvalContext::from(intr);
-        let a = UniversalTextureEvaluator.evaluate_float(alpha_tex, &ctx);
+        let a = UniversalTextureEvaluator.evaluate_float(&*self.alpha, &ctx);
         if a >= 1.0 {
             return false;
         }
@@ -75,15 +74,14 @@ impl LightTrait for DiffuseAreaLight {
     ) -> Option<LightLiSample> {
         let shape_ctx = ShapeSampleContext::new(ctx.pi, ctx.n, ctx.ns, 0.0);
         let ss = self.shape.sample_from_context(&shape_ctx, u)?;
-        let mut intr: SurfaceInteraction = ss.intr;
+        let mut intr = ss.intr;
-
+        intr.set_medium_interface(self.base.medium_interface);
-        intr.common.medium_interface = self.base.medium_interface;
         let p = intr.p();
         let n = intr.n();
-        let uv = intr.uv;
+        let uv = intr.get_common().uv;
 
-        let generic_intr = Interaction::Surface(intr);
+        // let generic_intr = Interaction::Surface(intr);
-        if self.alpha_masked(&generic_intr) {
+        if self.alpha_masked(&intr) {
             return None;
         }
 
@@ -94,7 +92,7 @@ impl LightTrait for DiffuseAreaLight {
             return None;
         }
 
-        Some(LightLiSample::new(le, wi, ss.pdf, generic_intr))
+        Some(LightLiSample::new(le, wi, ss.pdf, intr))
     }
 
     fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, _allow_incomplete_pdf: bool) -> Float {
@@ -121,16 +119,16 @@ impl LightTrait for DiffuseAreaLight {
         if self.alpha_masked(&intr) {
             return SampledSpectrum::new(0.);
         }
-        if let Some(image) = &self.image {
+        if !self.image.is_null() {
             let mut rgb = RGB::default();
             uv[1] = 1. - uv[1];
             for c in 0..3 {
-                rgb[c] = image.bilerp_channel(uv, c);
+                rgb[c] = self.image.bilerp_channel(uv, c as i32);
             }
 
             let spec = RGBIlluminantSpectrum::new(
                 self.image_color_space.as_ref().unwrap(),
-                RGB::clamp_zero(rgb),
+                rgb.clamp_zero(),
             );
 
             self.scale * spec.sample(lambda)
@@ -146,21 +144,21 @@ impl LightTrait for DiffuseAreaLight {
     #[cfg(not(target_os = "cuda"))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         let mut l = SampledSpectrum::new(0.);
-        if let Some(image) = &self.image {
+        if !self.image.is_null() {
-            for y in 0..image.resolution().y() {
+            for y in 0..self.image.resolution().y() {
-                for x in 0..image.resolution().x() {
+                for x in 0..self.image.resolution().x() {
                     let mut rgb = RGB::default();
                     for c in 0..3 {
-                        rgb[c] = image.get_channel(Point2i::new(x, y), c);
+                        rgb[c] = self.image.get_channel(Point2i::new(x, y), c as i32);
                     }
                     l += RGBIlluminantSpectrum::new(
                         self.image_color_space.as_ref().unwrap(),
-                        RGB::clamp_zero(rgb),
+                        rgb.clamp_zero(),
                     )
                     .sample(&lambda);
                 }
             }
-            l *= self.scale / (image.resolution().x() * image.resolution().y()) as Float;
+            l *= self.scale / (self.image.resolution().x() * self.image.resolution().y()) as Float;
         } else {
             l = self.lemit.sample(&lambda) * self.scale;
         }
@@ -177,10 +175,10 @@ impl LightTrait for DiffuseAreaLight {
     fn bounds(&self) -> Option<LightBounds> {
         let mut phi = 0.;
         if !self.image.is_null() {
-            for y in 0..image.resolution.y() {
+            for y in 0..self.image.resolution.y() {
-                for x in 0..image.resolution.x() {
+                for x in 0..self.image.resolution.x() {
                     for c in 0..3 {
-                        phi += image.get_channel(Point2i::new(x, y), c);
+                        phi += self.image.get_channel(Point2i::new(x, y), c);
                     }
                 }
             }

shared/src/lights/distant.rs

@@ -3,8 +3,9 @@ use crate::core::geometry::{
 };
 use crate::core::interaction::{Interaction, InteractionBase, SimpleInteraction};
 use crate::core::light::{LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait};
+use crate::core::spectrum::SpectrumTrait;
 use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, SampledWavelengths};
-use crate::utils::Ptr;
+use crate::utils::{ArenaPtr, Ptr};
 use crate::{Float, PI};
 
 #[repr(C)]
@@ -58,13 +59,10 @@ impl LightTrait for DistantLight {
         let p_outside = ctx.p() + wi * 2. * self.scene_radius;
 
         let li = self.scale * self.lemit.sample(lambda);
-        let intr = SimpleInteraction::new(
+        let base = InteractionBase::new_boundary(p_outside, 0.0, self.base.medium_interface);
-            Point3fi::new_from_point(p_outside),
+        let intr = SimpleInteraction::new(base);
-            0.0,
-            self.base.medium_interface,
-        );
 
-        Some(LightLiSample::new(li, wi, 1., intr))
+        Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
     }
 
     fn pdf_li(

shared/src/lights/goniometric.rs

@@ -4,11 +4,11 @@ use crate::core::light::{
     LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType,
 };
 use crate::core::medium::MediumInterface;
-use crate::core::spectrum::Spectrum;
+use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, SampledWavelengths};
-use crate::utils::Transform;
 use crate::utils::math::equal_area_sphere_to_square;
 use crate::utils::sampling::PiecewiseConstant2D;
+use crate::utils::{Ptr, Transform};
 use crate::{Float, PI};
 
 #[derive(Debug, Clone)]
@@ -16,37 +16,11 @@ pub struct GoniometricLight {
     pub base: LightBase,
     iemit: DenselySampledSpectrum,
     scale: Float,
-    image: *const Image,
+    image: Ptr<Image>,
-    distrib: *const PiecewiseConstant2D,
+    distrib: Ptr<PiecewiseConstant2D>,
 }
 
 impl GoniometricLight {
-    #[cfg(not(target_os = "cuda"))]
-    pub fn new(
-        render_from_light: &Transform,
-        medium_interface: &MediumInterface,
-        iemit: Spectrum,
-        scale: Float,
-        image: Image,
-    ) -> Self {
-        let base = LightBase::new(
-            LightType::DeltaPosition,
-            render_from_light,
-            medium_interface,
-        );
-
-        let i_interned = LightBase::lookup_spectrum(&iemit);
-        let d = image.get_sampling_distribution_uniform();
-        let distrib = PiecewiseConstant2D::new_with_data(&d);
-        Self {
-            base,
-            iemit: i_interned,
-            scale,
-            image,
-            distrib,
-        }
-    }
-
     pub fn i(&self, w: Vector3f, lambda: &SampledWavelengths) -> SampledSpectrum {
         let uv = equal_area_sphere_to_square(w);
         self.scale * self.iemit.sample(lambda) * self.image.lookup_nearest_channel(uv, 0)

shared/src/lights/infinite.rs

@@ -1,5 +1,8 @@
 use crate::{
-    core::geometry::Frame,
+    core::{
+        geometry::{Frame, VectorLike},
+        interaction::InteractionBase,
+    },
     spectra::{RGBColorSpace, RGBIlluminantSpectrum},
     utils::{
         math::{clamp, equal_area_sphere_to_square, equal_area_square_to_sphere, square},
@@ -21,8 +24,9 @@ use crate::core::light::{
 };
 use crate::core::medium::{Medium, MediumInterface};
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
-use crate::spectra::{SampledSpectrum, SampledWavelengths};
+use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, SampledWavelengths};
 use crate::utils::Transform;
+use crate::utils::ptr::Ptr;
 use crate::{Float, PI};
 use std::sync::Arc;
 
@@ -30,30 +34,14 @@ use std::sync::Arc;
 #[derive(Debug, Copy, Clone)]
 pub struct InfiniteUniformLight {
     pub base: LightBase,
-    pub lemit: u32,
+    pub lemit: Ptr<DenselySampledSpectrum>,
     pub scale: Float,
     pub scene_center: Point3f,
     pub scene_radius: Float,
 }
 
-#[cfg(not(target_os = "cuda"))]
+unsafe impl Send for InfiniteUniformLight {}
-impl InfiniteUniformLight {
+unsafe impl Sync for InfiniteUniformLight {}
-    pub fn new(render_from_light: Transform, le: Spectrum, scale: Float) -> Self {
-        let base = LightBase::new(
-            LightType::Infinite,
-            &render_from_light,
-            &MediumInterface::default(),
-        );
-        let lemit = LightBase::lookup_spectrum(&le);
-        Self {
-            base,
-            lemit,
-            scale,
-            scene_center: Point3f::default(),
-            scene_radius: 0.,
-        }
-    }
-}
 
 impl LightTrait for InfiniteUniformLight {
     fn base(&self) -> &LightBase {
@@ -71,10 +59,12 @@ impl LightTrait for InfiniteUniformLight {
         }
         let wi = sample_uniform_sphere(u);
         let pdf = uniform_sphere_pdf();
-        let intr_simple = SimpleInteraction::new_interface(
+        let base = InteractionBase::new_boundary(
             ctx.p() + wi * (2. * self.scene_radius),
-            Some(MediumInterface::default()),
+            0.,
+            MediumInterface::default(),
         );
+        let intr_simple = SimpleInteraction::new(base);
 
         let intr = Interaction::Simple(intr_simple);
         Some(LightLiSample::new(
@@ -111,12 +101,18 @@ impl LightTrait for InfiniteUniformLight {
     fn le(&self, _ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum {
         self.scale * self.lemit.sample(lambda)
     }
+
+    #[cfg(not(target_os = "cuda"))]
     fn preprocess(&mut self, _scene_bounds: &Bounds3f) {
         todo!()
     }
+
+    #[cfg(not(target_os = "cuda"))]
     fn bounds(&self) -> Option<LightBounds> {
         todo!()
     }
+
+    #[cfg(not(target_os = "cuda"))]
     fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum {
         4. * PI * PI * square(self.scene_radius) * self.scale * self.lemit.sample(&lambda)
     }
@@ -126,10 +122,10 @@ impl LightTrait for InfiniteUniformLight {
 #[derive(Clone, Copy, Debug)]
 pub struct InfiniteImageLight {
     pub base: LightBase,
-    pub image: *const Image,
+    pub image: Ptr<Image>,
-    pub image_color_space: *const RGBColorSpace,
+    pub image_color_space: Ptr<RGBColorSpace>,
-    pub distrib: *const PiecewiseConstant2D,
+    pub distrib: Ptr<PiecewiseConstant2D>,
-    pub compensated_distrib: *const PiecewiseConstant2D,
+    pub compensated_distrib: Ptr<PiecewiseConstant2D>,
     pub scale: Float,
     pub scene_radius: Float,
     pub scene_center: Point3f,
@@ -139,26 +135,16 @@ unsafe impl Send for InfiniteImageLight {}
 unsafe impl Sync for InfiniteImageLight {}
 
 impl InfiniteImageLight {
-    #[inline(always)]
-    fn color_space(&self) -> &RGBColorSpace {
-        unsafe { &*self.image_color_space }
-    }
-
-    #[inline(always)]
-    fn image(&self) -> &Image {
-        unsafe { &*self.image }
-    }
-
     fn image_le(&self, uv: Point2f, lambda: &SampledWavelengths) -> SampledSpectrum {
         let mut rgb = RGB::default();
         for c in 0..3 {
             rgb[c] = self.image.lookup_nearest_channel_with_wrap(
                 uv,
-                c,
+                c as i32,
                 WrapMode::OctahedralSphere.into(),
             );
         }
-        let spec = RGBIlluminantSpectrum::new(self.color_space(), RGB::clamp_zero(rgb));
+        let spec = RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero());
         self.scale * spec.sample(lambda)
     }
 }
@@ -190,13 +176,14 @@ impl LightTrait for InfiniteImageLight {
         let pdf = map_pdf / (4. * PI);
 
         // Return radiance value for infinite light direction
-        let mut simple_intr = SimpleInteraction::new_interface(
+        let base = InteractionBase::new_boundary(
             ctx.p() + wi * (2. * self.scene_radius),
-            Some(MediumInterface::default()),
+            0.,
+            self.base.medium_interface,
         );
-
+        let simple_intr = SimpleInteraction::new(base);
-        simple_intr.common.medium_interface = Some(self.base.medium_interface.clone());
         let intr = Interaction::Simple(simple_intr);
+
         Some(LightLiSample::new(self.image_le(uv, lambda), wi, pdf, intr))
     }
 
@@ -243,14 +230,11 @@ impl LightTrait for InfiniteImageLight {
                 for c in 0..3 {
                     rgb[c] = self.image.get_channel_with_wrap(
                         Point2i::new(u, v),
-                        c,
+                        c as i32,
                         WrapMode::OctahedralSphere.into(),
                     );
                 }
-                sum_l += RGBIlluminantSpectrum::new(
+                sum_l += RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero())
-                    self.image_color_space.as_ref(),
-                    RGB::clamp_zero(rgb),
-                )
                     .sample(&lambda);
             }
         }
@@ -274,8 +258,8 @@ impl LightTrait for InfiniteImageLight {
 #[derive(Debug, Copy, Clone)]
 pub struct InfinitePortalLight {
     pub base: LightBase,
-    pub image: Image,
+    pub image: Ptr<Image>,
-    pub image_color_space: RGBColorSpace,
+    pub image_color_space: Ptr<RGBColorSpace>,
     pub scale: Float,
     pub portal: [Point3f; 4],
     pub portal_frame: Frame,
@@ -288,10 +272,9 @@ impl InfinitePortalLight {
     pub fn image_lookup(&self, uv: Point2f, lambda: &SampledWavelengths) -> SampledSpectrum {
         let mut rgb = RGB::default();
         for c in 0..3 {
-            rgb[c] = self.image.lookup_nearest_channel(uv, c)
+            rgb[c] = self.image.lookup_nearest_channel(uv, c as i32)
         }
-        let spec =
+        let spec = RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero());
-            RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb));
         self.scale * spec.sample(lambda)
     }
 
@@ -363,7 +346,8 @@ impl LightTrait for InfinitePortalLight {
         let pdf = map_pdf / duv_dw;
         let l = self.image_lookup(uv, lambda);
         let pl = ctx.p() + 2. * self.scene_radius * wi;
-        let sintr = SimpleInteraction::new_interface(pl, Some(self.base.medium_interface.clone()));
+        let base = InteractionBase::new_boundary(pl, 0., self.base.medium_interface);
+        let sintr = SimpleInteraction::new(base);
         let intr = Interaction::Simple(sintr);
         Some(LightLiSample::new(l, wi, pdf, intr))
     }

shared/src/lights/point.rs

@@ -1,9 +1,13 @@
-use crate::core::geometry::{Bounds3f, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector3f};
+use crate::core::geometry::{
-use crate::core::interaction::{Interaction, SimpleInteraction};
+    Bounds3f, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector3f, VectorLike,
+};
+use crate::core::interaction::{Interaction, InteractionBase, SimpleInteraction};
 use crate::core::light::{
     Light, LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType,
 };
+use crate::core::spectrum::SpectrumTrait;
 use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, SampledWavelengths};
+use crate::utils::ptr::Ptr;
 use crate::{Float, PI};
 
 #[repr(C)]
@@ -11,25 +15,7 @@ use crate::{Float, PI};
 pub struct PointLight {
     pub base: LightBase,
     pub scale: Float,
-    pub i: *const DenselySampledSpectrum,
+    pub i: Ptr<DenselySampledSpectrum>,
-}
-
-impl PointLight {
-    fn sample_li_base(
-        &self,
-        ctx_p: Point3f,
-        lambda: &SampledWavelengths,
-    ) -> (SampledSpectrum, Vector3f, Float, Point3fi) {
-        let pi = self
-            .base
-            .render_from_light
-            .apply_to_interval(&Point3fi::default());
-        let p: Point3f = pi.into();
-        let wi = (p - ctx_p).normalize();
-        let spectrum = DenselySampledSpectrum::from_array(&self.i_coeffs);
-        let li = self.scale * spectrum.sample(lambda) / p.distance_squared(ctx_p);
-        (li, wi, 1.0, pi)
-    }
 }
 
 impl LightTrait for PointLight {
@@ -51,7 +37,8 @@ impl LightTrait for PointLight {
         let p: Point3f = pi.into();
         let wi = (p - ctx.p()).normalize();
         let li = self.scale * self.i.sample(lambda) / p.distance_squared(ctx.p());
-        let intr = SimpleInteraction::new(pi, 0.0, Some(self.base.medium_interface.clone()));
+        let base = InteractionBase::new_boundary(p, 0., self.base.medium_interface);
+        let intr = SimpleInteraction::new(base);
         Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
     }
 

shared/src/lights/projection.rs

@@ -8,8 +8,10 @@ use crate::core::light::{
     LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType,
 };
 use crate::core::medium::MediumInterface;
+use crate::core::spectrum::SpectrumTrait;
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::utils::math::{radians, square};
+use crate::utils::ptr::Ptr;
 use crate::{
     spectra::{RGBColorSpace, RGBIlluminantSpectrum},
     utils::{Transform, sampling::PiecewiseConstant2D},
@@ -25,14 +27,12 @@ pub struct ProjectionLight {
     pub screen_from_light: Transform,
     pub light_from_screen: Transform,
     pub a: Float,
-    pub image: *const Image,
+    pub image: Ptr<Image>,
-    pub distrib: *const PiecewiseConstant2D,
+    pub distrib: Ptr<PiecewiseConstant2D>,
-    pub image_color_space: *const RGBColorSpace,
+    pub image_color_space: Ptr<RGBColorSpace>,
 }
 
 impl ProjectionLight {
-    #[cfg(not(target_os = "cuda"))]
-
     pub fn i(&self, w: Vector3f, lambda: SampledWavelengths) -> SampledSpectrum {
         if w.z() < self.hither {
             return SampledSpectrum::new(0.);
@@ -44,9 +44,9 @@ impl ProjectionLight {
         let uv = Point2f::from(self.screen_bounds.offset(&Point2f::new(ps.x(), ps.y())));
         let mut rgb = RGB::default();
         for c in 0..3 {
-            rgb[c] = self.image.lookup_nearest_channel(uv, c);
+            rgb[c] = self.image.lookup_nearest_channel(uv, c as i32);
         }
-        let s = RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb));
+        let s = RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero());
         self.scale * s.sample(&lambda)
     }
 }
@@ -107,15 +107,10 @@ impl LightTrait for ProjectionLight {
                 let dwda = cos_theta(w).powi(3);
                 let mut rgb = RGB::default();
                 for c in 0..3 {
-                    rgb[c] = self.image.get_channel(Point2i::new(x, y), c);
+                    rgb[c] = self.image.get_channel(Point2i::new(x, y), c as i32);
                 }
 
-                let s = unsafe {
+                let s = RGBIlluminantSpectrum::new(&*self.image_color_space, rgb.clamp_zero());
-                    RGBIlluminantSpectrum::new(
-                        self.image_color_space.as_ref(),
-                        RGB::clamp_zero(rgb),
-                    );
-                };
                 sum += s.sample(&lambda) * dwda;
             }
         }

shared/src/lights/sampler.rs

@@ -160,7 +160,10 @@ pub struct SampledLight {
 
 impl SampledLight {
     pub fn new(light: Light, p: Float) -> Self {
-        Self { light, p }
+        Self {
+            light: Ptr::from(&light),
+            p,
+        }
     }
 }
 
@@ -211,7 +214,7 @@ impl LightSamplerTrait for UniformLightSampler {
 
         let light_index = (u as u32 * self.lights_len).min(self.lights_len - 1) as usize;
         Some(SampledLight {
-            light: self.light(light_index),
+            light: Ptr::from(&self.light(light_index)),
             p: 1. / self.lights_len as Float,
         })
     }

shared/src/lights/spot.rs

@@ -1,7 +1,7 @@
 use crate::core::geometry::{
     Bounds3f, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector3f, VectorLike,
 };
-use crate::core::interaction::{Interaction, InteractionTrait, SimpleInteraction};
+use crate::core::interaction::{Interaction, InteractionBase, InteractionTrait, SimpleInteraction};
 use crate::core::light::{LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait};
 use crate::core::spectrum::SpectrumTrait;
 use crate::spectra::{DenselySampledSpectrum, SampledSpectrum, SampledWavelengths};
@@ -51,7 +51,8 @@ impl LightTrait for SpotLight {
         let w_light = self.base.render_from_light.apply_inverse_vector(-wi);
         let li = self.i(w_light, lambda) / p.distance_squared(ctx.p());
 
-        let intr = SimpleInteraction::new(pi, 0.0, Ptr::from(&self.base.medium_interface));
+        let base = InteractionBase::new_boundary(p, 0., self.base.medium_interface);
+        let intr = SimpleInteraction::new(base);
         Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr)))
     }
 

shared/src/materials/coated.rs

@@ -1,28 +1,30 @@
-use crate::core::bssrdf::BSSRDF;
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::BSSRDF;
+use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::math::clamp;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct CoatedDiffuseMaterial {
-    pub normal_map: *const Image,
+    pub normal_map: Ptr<Image>,
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub reflectance: RelPtr<GPUSpectrumTexture>,
+    pub reflectance: Ptr<GPUSpectrumTexture>,
-    pub albedo: RelPtr<GPUSpectrumTexture>,
+    pub albedo: Ptr<GPUSpectrumTexture>,
-    pub u_roughness: RelPtr<GPUFloatTexture>,
+    pub u_roughness: Ptr<GPUFloatTexture>,
-    pub v_roughness: RelPtr<GPUFloatTexture>,
+    pub v_roughness: Ptr<GPUFloatTexture>,
-    pub thickness: RelPtr<GPUFloatTexture>,
+    pub thickness: Ptr<GPUFloatTexture>,
-    pub g: RelPtr<GPUFloatTexture>,
+    pub g: Ptr<GPUFloatTexture>,
-    pub eta: RelPtr<Spectrum>,
+    pub eta: Ptr<Spectrum>,
     pub remap_roughness: bool,
     pub max_depth: usize,
     pub n_samples: usize,
@@ -32,29 +34,29 @@ impl CoatedDiffuseMaterial {
     #[allow(clippy::too_many_arguments)]
     #[cfg(not(target_os = "cuda"))]
     pub fn new(
-        reflectance: GPUSpectrumTexture,
+        reflectance: &GPUSpectrumTexture,
-        u_roughness: GPUFloatTexture,
+        u_roughness: &GPUFloatTexture,
-        v_roughness: GPUFloatTexture,
+        v_roughness: &GPUFloatTexture,
-        thickness: GPUFloatTexture,
+        thickness: &GPUFloatTexture,
-        albedo: GPUSpectrumTexture,
+        albedo: &GPUSpectrumTexture,
-        g: GPUFloatTexture,
+        g: &GPUFloatTexture,
-        eta: Spectrum,
+        eta: &Spectrum,
-        displacement: GPUFloatTexture,
+        displacement: &GPUFloatTexture,
-        normal_map: *const Image,
+        normal_map: &Image,
         remap_roughness: bool,
         max_depth: usize,
         n_samples: usize,
     ) -> Self {
         Self {
-            displacement,
+            displacement: Ptr::from(displacement),
-            normal_map,
+            normal_map: Ptr::from(normal_map),
-            reflectance,
+            reflectance: Ptr::from(reflectance),
-            albedo,
+            albedo: Ptr::from(albedo),
-            u_roughness,
+            u_roughness: Ptr::from(u_roughness),
-            v_roughness,
+            v_roughness: Ptr::from(v_roughness),
-            thickness,
+            thickness: Ptr::from(thickness),
-            g,
+            g: Ptr::from(g),
-            eta,
+            eta: Ptr::from(eta),
             remap_roughness,
             max_depth,
             n_samples,
@@ -113,7 +115,7 @@ impl MaterialTrait for CoatedDiffuseMaterial {
             self.n_samples,
         ));
 
-        BSDF::new(ctx.ns, ctx.dpdus, Some(bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
     }
 
     fn get_bssrdf<T>(
@@ -132,11 +134,11 @@ impl MaterialTrait for CoatedDiffuseMaterial {
         )
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.normal_map)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 
@@ -148,19 +150,19 @@ impl MaterialTrait for CoatedDiffuseMaterial {
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct CoatedConductorMaterial {
-    normal_map: *const Image,
+    normal_map: Ptr<Image>,
-    displacement: RelPtr<GPUFloatTexture>,
+    displacement: Ptr<GPUFloatTexture>,
-    interface_uroughness: RelPtr<GPUFloatTexture>,
+    interface_uroughness: Ptr<GPUFloatTexture>,
-    interface_vroughness: RelPtr<GPUFloatTexture>,
+    interface_vroughness: Ptr<GPUFloatTexture>,
-    thickness: RelPtr<GPUFloatTexture>,
+    thickness: Ptr<GPUFloatTexture>,
-    interface_eta: RelPtr<Spectrum>,
+    interface_eta: Ptr<Spectrum>,
-    g: RelPtr<GPUFloatTexture>,
+    g: Ptr<GPUFloatTexture>,
-    albedo: RelPtr<GPUSpectrumTexture>,
+    albedo: Ptr<GPUSpectrumTexture>,
-    conductor_uroughness: RelPtr<GPUFloatTexture>,
+    conductor_uroughness: Ptr<GPUFloatTexture>,
-    conductor_vroughness: RelPtr<GPUFloatTexture>,
+    conductor_vroughness: Ptr<GPUFloatTexture>,
-    conductor_eta: RelPtr<GPUSpectrumTexture>,
+    conductor_eta: Ptr<GPUSpectrumTexture>,
-    k: RelPtr<GPUSpectrumTexture>,
+    k: Ptr<GPUSpectrumTexture>,
-    reflectance: RelPtr<GPUSpectrumTexture>,
+    reflectance: Ptr<GPUSpectrumTexture>,
     remap_roughness: bool,
     max_depth: u32,
     n_samples: u32,
@@ -170,37 +172,37 @@ impl CoatedConductorMaterial {
     #[allow(clippy::too_many_arguments)]
     #[cfg(not(target_os = "cuda"))]
     pub fn new(
-        displacement: GPUFloatTexture,
+        normal_map: &Image,
-        normal_map: *const Image,
+        displacement: &GPUFloatTexture,
-        interface_uroughness: GPUFloatTexture,
+        interface_uroughness: &GPUFloatTexture,
-        interface_vroughness: GPUFloatTexture,
+        interface_vroughness: &GPUFloatTexture,
-        thickness: GPUFloatTexture,
+        thickness: &GPUFloatTexture,
-        interface_eta: Spectrum,
+        interface_eta: &Spectrum,
-        g: GPUFloatTexture,
+        g: &GPUFloatTexture,
-        albedo: GPUSpectrumTexture,
+        albedo: &GPUSpectrumTexture,
-        conductor_uroughness: GPUFloatTexture,
+        conductor_uroughness: &GPUFloatTexture,
-        conductor_vroughness: GPUFloatTexture,
+        conductor_vroughness: &GPUFloatTexture,
-        conductor_eta: Option<GPUSpectrumTexture>,
+        conductor_eta: &GPUSpectrumTexture,
-        k: Option<GPUSpectrumTexture>,
+        k: &GPUSpectrumTexture,
-        reflectance: GPUSpectrumTexture,
+        reflectance: &GPUSpectrumTexture,
         remap_roughness: bool,
-        max_depth: usize,
+        max_depth: u32,
-        n_samples: usize,
+        n_samples: u32,
     ) -> Self {
         Self {
-            displacement,
+            displacement: Ptr::from(displacement),
-            normal_map,
+            normal_map: Ptr::from(normal_map),
-            interface_uroughness,
+            interface_uroughness: Ptr::from(interface_uroughness),
-            interface_vroughness,
+            interface_vroughness: Ptr::from(interface_vroughness),
-            thickness,
+            thickness: Ptr::from(thickness),
-            interface_eta,
+            interface_eta: Ptr::from(interface_eta),
-            g,
+            g: Ptr::from(g),
-            albedo,
+            albedo: Ptr::from(albedo),
-            conductor_uroughness,
+            conductor_uroughness: Ptr::from(conductor_uroughness),
-            conductor_vroughness,
+            conductor_vroughness: Ptr::from(conductor_vroughness),
-            conductor_eta,
+            conductor_eta: Ptr::from(conductor_eta),
-            k,
+            k: Ptr::from(k),
-            reflectance,
+            reflectance: Ptr::from(reflectance),
             remap_roughness,
             max_depth,
             n_samples,
@@ -235,12 +237,12 @@ impl MaterialTrait for CoatedConductorMaterial {
             ieta = 1.;
         }
 
-        let (mut ce, mut ck) = if let Some(eta_tex) = &self.conductor_eta {
+        let (mut ce, mut ck) = if !self.conductor_eta.is_null() {
             let k_tex = self
                 .k
                 .as_ref()
                 .expect("CoatedConductor: 'k' must be provided if 'conductor_eta' is present");
-            let ce = tex_eval.evaluate_spectrum(eta_tex, ctx, lambda);
+            let ce = tex_eval.evaluate_spectrum(&self.conductor_eta, ctx, lambda);
             let ck = tex_eval.evaluate_spectrum(k_tex, ctx, lambda);
             (ce, ck)
         } else {
@@ -280,10 +282,10 @@ impl MaterialTrait for CoatedConductorMaterial {
             thick,
             a,
             gg,
-            self.max_depth,
+            self.max_depth as usize,
-            self.n_samples,
+            self.n_samples as usize,
         ));
-        BSDF::new(ctx.ns, ctx.dpdus, Some(bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
     }
 
     fn get_bssrdf<T>(
@@ -307,27 +309,28 @@ impl MaterialTrait for CoatedConductorMaterial {
 
         let mut spectrum_textures = Vec::with_capacity(4);
 
-        spectrum_textures.push(&self.albedo);
+        spectrum_textures.push(self.albedo);
 
-        if let Some(eta) = &self.conductor_eta {
+        if !self.conductor_eta.is_null() {
-            spectrum_textures.push(eta);
+            spectrum_textures.push(self.conductor_eta);
-        }
-        if let Some(k) = &self.k {
-            spectrum_textures.push(k);
         }
 
-        if self.conductor_eta.is_none() {
+        if !self.k.is_null() {
+            spectrum_textures.push(self.k);
+        }
+
+        if !self.conductor_eta.is_null() {
             spectrum_textures.push(self.reflectance);
         }
 
         tex_eval.can_evaluate(&float_textures, &spectrum_textures)
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.normal_map)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 

shared/src/materials/complex.rs

@@ -1,9 +1,11 @@
 use crate::Float;
-use crate::core::bssrdf::{BSSRDF, BSSRDFTable};
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF, HairBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
+    MeasuredBxDF, MeasuredBxDFData,
-    HairBxDF, MeasuredBxDF, MeasuredBxDFData,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::{BSSRDF, BSSRDFTable};
+use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
@@ -11,33 +13,33 @@ use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
 use crate::textures::GPUSpectrumMixTexture;
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::math::clamp;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct HairMaterial {
-    pub sigma_a: RelPtr<GPUSpectrumTexture>,
+    pub sigma_a: Ptr<GPUSpectrumTexture>,
-    pub color: RelPtr<GPUSpectrumTexture>,
+    pub color: Ptr<GPUSpectrumTexture>,
-    pub eumelanin: RelPtr<GPUFloatTexture>,
+    pub eumelanin: Ptr<GPUFloatTexture>,
-    pub pheomelanin: RelPtr<GPUFloatTexture>,
+    pub pheomelanin: Ptr<GPUFloatTexture>,
-    pub eta: RelPtr<GPUFloatTexture>,
+    pub eta: Ptr<GPUFloatTexture>,
-    pub beta_m: RelPtr<GPUFloatTexture>,
+    pub beta_m: Ptr<GPUFloatTexture>,
-    pub beta_n: RelPtr<GPUFloatTexture>,
+    pub beta_n: Ptr<GPUFloatTexture>,
-    pub alpha: RelPtr<GPUFloatTexture>,
+    pub alpha: Ptr<GPUFloatTexture>,
 }
 
 impl HairMaterial {
     #[cfg(not(target_os = "cuda"))]
     pub fn new(
-        sigma_a: RelPtr<GPUSpectrumTexture>,
+        sigma_a: Ptr<GPUSpectrumTexture>,
-        color: RelPtr<GPUSpectrumTexture>,
+        color: Ptr<GPUSpectrumTexture>,
-        eumelanin: RelPtr<GPUFloatTexture>,
+        eumelanin: Ptr<GPUFloatTexture>,
-        pheomelanin: RelPtr<GPUFloatTexture>,
+        pheomelanin: Ptr<GPUFloatTexture>,
-        eta: RelPtr<GPUFloatTexture>,
+        eta: Ptr<GPUFloatTexture>,
-        beta_m: RelPtr<GPUFloatTexture>,
+        beta_m: Ptr<GPUFloatTexture>,
-        beta_n: RelPtr<GPUFloatTexture>,
+        beta_n: Ptr<GPUFloatTexture>,
-        alpha: RelPtr<GPUFloatTexture>,
+        alpha: Ptr<GPUFloatTexture>,
     ) -> Self {
         Self {
             sigma_a,
@@ -74,12 +76,12 @@ impl MaterialTrait for HairMaterial {
         todo!()
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
         todo!()
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
-        RelPtr::null()
+        Ptr::null()
     }
 
     fn has_subsurface_scattering(&self) -> bool {
@@ -90,9 +92,9 @@ impl MaterialTrait for HairMaterial {
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct MeasuredMaterial {
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub normal_map: *const Image,
+    pub normal_map: Ptr<Image>,
-    pub brdf: *const MeasuredBxDFData,
+    pub brdf: Ptr<MeasuredBxDFData>,
 }
 
 impl MaterialTrait for MeasuredMaterial {
@@ -100,9 +102,10 @@ impl MaterialTrait for MeasuredMaterial {
         &self,
         _tex_eval: &T,
         _ctx: &MaterialEvalContext,
-        lambda: &SampledWavelengths,
+        _lambda: &SampledWavelengths,
     ) -> BSDF {
-        MeasuredBxDF::new(self.brdf, lambda)
+        // MeasuredBxDF::new(&self.brdf, lambda)
+        todo!()
     }
 
     fn get_bssrdf<T>(
@@ -118,11 +121,11 @@ impl MaterialTrait for MeasuredMaterial {
         true
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.normal_map)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 
@@ -134,16 +137,16 @@ impl MaterialTrait for MeasuredMaterial {
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct SubsurfaceMaterial {
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub normal_map: Ptr<Image>,
-    pub normal_map: *const Image,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub sigma_a: RelPtr<GPUSpectrumTexture>,
+    pub sigma_a: Ptr<GPUSpectrumTexture>,
-    pub sigma_s: RelPtr<GPUSpectrumMixTexture>,
+    pub sigma_s: Ptr<GPUSpectrumMixTexture>,
-    pub reflectance: RelPtr<GPUSpectrumMixTexture>,
+    pub reflectance: Ptr<GPUSpectrumMixTexture>,
-    pub mfp: RelPtr<GPUSpectrumMixTexture>,
+    pub mfp: Ptr<GPUSpectrumMixTexture>,
     pub eta: Float,
     pub scale: Float,
-    pub u_roughness: RelPtr<GPUFloatTexture>,
+    pub u_roughness: Ptr<GPUFloatTexture>,
-    pub v_roughness: RelPtr<GPUFloatTexture>,
+    pub v_roughness: Ptr<GPUFloatTexture>,
     pub remap_roughness: bool,
     pub table: BSSRDFTable,
 }
@@ -169,12 +172,15 @@ impl MaterialTrait for SubsurfaceMaterial {
     fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool {
         todo!()
     }
-    fn get_normal_map(&self) -> *const Image {
+
+    fn get_normal_map(&self) -> Option<&Image> {
         todo!()
     }
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         todo!()
     }
+
     fn has_subsurface_scattering(&self) -> bool {
         true
     }

shared/src/materials/conductor.rs

@@ -1,28 +1,29 @@
-use crate::core::bssrdf::BSSRDF;
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF, HairBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
-    HairBxDF,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::BSSRDF;
+use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::math::clamp;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct ConductorMaterial {
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub eta: RelPtr<GPUSpectrumTexture>,
+    pub eta: Ptr<GPUSpectrumTexture>,
-    pub k: RelPtr<GPUSpectrumTexture>,
+    pub k: Ptr<GPUSpectrumTexture>,
-    pub reflectance: RelPtr<GPUSpectrumTexture>,
+    pub reflectance: Ptr<GPUSpectrumTexture>,
-    pub u_roughness: RelPtr<GPUFloatTexture>,
+    pub u_roughness: Ptr<GPUFloatTexture>,
-    pub v_roughness: RelPtr<GPUFloatTexture>,
+    pub v_roughness: Ptr<GPUFloatTexture>,
     pub remap_roughness: bool,
-    pub normal_map: *const Image,
+    pub normal_map: Ptr<Image>,
 }
 
 impl MaterialTrait for ConductorMaterial {
@@ -49,11 +50,11 @@ impl MaterialTrait for ConductorMaterial {
         )
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
         todo!()
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         todo!()
     }
 

shared/src/materials/dielectric.rs

@@ -1,26 +1,27 @@
-use crate::core::bssrdf::BSSRDF;
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF, HairBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
-    HairBxDF,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::BSSRDF;
+use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::math::clamp;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct DielectricMaterial {
-    normal_map: *const Image,
+    normal_map: Ptr<Image>,
-    displacement: GPUFloatTexture,
+    displacement: Ptr<GPUFloatTexture>,
-    u_roughness: GPUFloatTexture,
+    u_roughness: Ptr<GPUFloatTexture>,
-    v_roughness: GPUFloatTexture,
+    v_roughness: Ptr<GPUFloatTexture>,
+    eta: Ptr<Spectrum>,
     remap_roughness: bool,
-    eta: Spectrum,
 }
 
 impl MaterialTrait for DielectricMaterial {
@@ -50,7 +51,7 @@ impl MaterialTrait for DielectricMaterial {
         let distrib = TrowbridgeReitzDistribution::new(u_rough, v_rough);
         let bxdf = BxDF::Dielectric(DielectricBxDF::new(sampled_eta, distrib));
 
-        BSDF::new(ctx.ns, ctx.dpdus, Some(bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
     }
 
     fn get_bssrdf<T>(
@@ -66,11 +67,11 @@ impl MaterialTrait for DielectricMaterial {
         tex_eval.can_evaluate(&[self.u_roughness, self.v_roughness], &[])
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.normal_map)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 
@@ -82,9 +83,9 @@ impl MaterialTrait for DielectricMaterial {
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct ThinDielectricMaterial {
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub normal_map: *const Image,
+    pub normal_map: Ptr<Image>,
-    pub eta: RelPtr<Spectrum>,
+    pub eta: Ptr<Spectrum>,
 }
 impl MaterialTrait for ThinDielectricMaterial {
     fn get_bsdf<T: TextureEvaluator>(
@@ -108,11 +109,11 @@ impl MaterialTrait for ThinDielectricMaterial {
         true
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.normal_map)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 

shared/src/materials/diffuse.rs

@@ -1,24 +1,25 @@
 use crate::Float;
-use crate::core::bssrdf::BSSRDF;
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF, HairBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
-    HairBxDF,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::BSSRDF;
+use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::math::clamp;
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct DiffuseMaterial {
-    pub normal_map: *const Image,
+    pub normal_map: Ptr<Image>,
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub reflectance: RelPtr<GPUSpectrumTexture>,
+    pub reflectance: Ptr<GPUSpectrumTexture>,
 }
 
 impl MaterialTrait for DiffuseMaterial {
@@ -30,7 +31,7 @@ impl MaterialTrait for DiffuseMaterial {
     ) -> BSDF {
         let r = tex_eval.evaluate_spectrum(&self.reflectance, ctx, lambda);
         let bxdf = BxDF::Diffuse(DiffuseBxDF::new(r));
-        BSDF::new(ctx.ns, ctx.dpdus, Some(bxdf))
+        BSDF::new(ctx.ns, ctx.dpdus, Ptr::from(&bxdf))
     }
 
     fn get_bssrdf<T>(
@@ -46,11 +47,11 @@ impl MaterialTrait for DiffuseMaterial {
         tex_eval.can_evaluate(&[], &[self.reflectance])
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.normal_map)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 
@@ -62,10 +63,10 @@ impl MaterialTrait for DiffuseMaterial {
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct DiffuseTransmissionMaterial {
-    pub displacement: RelPtr<GPUFloatTexture>,
+    pub image: Ptr<Image>,
-    pub image: *const Image,
+    pub displacement: Ptr<GPUFloatTexture>,
-    pub reflectance: RelPtr<GPUFloatTexture>,
+    pub reflectance: Ptr<GPUFloatTexture>,
-    pub transmittance: RelPtr<GPUFloatTexture>,
+    pub transmittance: Ptr<GPUFloatTexture>,
     pub scale: Float,
 }
 
@@ -88,14 +89,14 @@ impl MaterialTrait for DiffuseTransmissionMaterial {
     }
 
     fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool {
-        tex_eval.can_evaluate(&[], &[self.reflectance, self.transmittance])
+        tex_eval.can_evaluate(&[self.reflectance, self.transmittance], &[])
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        self.normal_map
+        Some(&*self.image)
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         self.displacement
     }
 

shared/src/materials/mix.rs

@@ -1,23 +1,24 @@
-use crate::core::bssrdf::BSSRDF;
+use crate::bxdfs::{
-use crate::core::bxdf::{
+    CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF, HairBxDF,
-    BSDF, BxDF, CoatedConductorBxDF, CoatedDiffuseBxDF, ConductorBxDF, DielectricBxDF, DiffuseBxDF,
-    HairBxDF,
 };
+use crate::core::bsdf::BSDF;
+use crate::core::bssrdf::BSSRDF;
+use crate::core::bxdf::BxDF;
 use crate::core::image::Image;
 use crate::core::material::{Material, MaterialEvalContext, MaterialTrait};
 use crate::core::scattering::TrowbridgeReitzDistribution;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvaluator};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
 use crate::utils::hash::hash_float;
 use crate::utils::math::clamp;
+use crate::utils::{ArenaPtr, Ptr};
 
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct MixMaterial {
-    pub amount: RelPtr<GPUFloatTexture>,
+    pub amount: Ptr<GPUFloatTexture>,
-    pub materials: [RelPtr<Material>; 2],
+    pub materials: [ArenaPtr<Material>; 2],
 }
 
 impl MixMaterial {
@@ -51,7 +52,7 @@ impl MaterialTrait for MixMaterial {
         if let Some(mat) = self.choose_material(tex_eval, ctx) {
             mat.get_bsdf(tex_eval, ctx, lambda)
         } else {
-            BSDF::empty()
+            BSDF::default()
         }
     }
 
@@ -68,11 +69,11 @@ impl MaterialTrait for MixMaterial {
         tex_eval.can_evaluate(&[self.amount], &[])
     }
 
-    fn get_normal_map(&self) -> *const Image {
+    fn get_normal_map(&self) -> Option<&Image> {
-        core::ptr::null()
+        None
     }
 
-    fn get_displacement(&self) -> RelPtr<GPUFloatTexture> {
+    fn get_displacement(&self) -> Ptr<GPUFloatTexture> {
         panic!(
             "MixMaterial::get_displacement() shouldn't be called. \
                Displacement is not supported on Mix materials directly."

shared/src/shapes/bilinear.rs

@@ -11,6 +11,7 @@ use crate::utils::mesh::BilinearPatchMesh;
 use crate::utils::sampling::{
     bilinear_pdf, invert_spherical_rectangle_sample, sample_bilinear, sample_spherical_rectangle,
 };
+use core::ops::Add;
 
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
@@ -60,7 +61,11 @@ impl BilinearPatchShape {
     #[inline(always)]
     fn get_vertex_indices(&self) -> [usize; 4] {
         unsafe {
-            let base_ptr = self.mesh.vertex_indices.add((self.blp_index as usize) * 4);
+            let base_ptr = self
+                .mesh
+                .vertex_indices
+                .0
+                .add((self.blp_index as usize) * 4);
             [
                 *base_ptr.add(0) as usize,
                 *base_ptr.add(1) as usize,
@@ -317,7 +322,7 @@ impl BilinearPatchShape {
         dpdv: &mut Vector3f,
     ) -> (Point2f, Option<TextureDerivative>) {
         let Some(uvs) = patch_uvs else {
-            return (uv, TextureDerivative::default());
+            return (uv, Some(TextureDerivative::default()));
         };
         let uv00 = uvs[0];
         let uv01 = uvs[1];
@@ -485,7 +490,7 @@ impl BilinearPatchShape {
         u: Point2f,
         corner_dirs: &[Vector3f; 4],
     ) -> Option<ShapeSample> {
-        let (p00, p10, p01, p11) = (corners[0], corners[1], corners[2], corners[3]);
+        let (p00, p10, p01, _p11) = (corners[0], corners[1], corners[2], corners[3]);
         let mut pdf = 1.;
         if ctx.ns != Normal3f::zero() {
             let w = [

shared/src/shapes/cylinder.rs

@@ -266,7 +266,7 @@ impl ShapeTrait for CylinderShape {
             (p_obj.z() - self.z_min) / (self.z_max - self.z_min),
         );
         Some(ShapeSample {
-            intr: SurfaceInteraction::new_simple(pi, n, uv),
+            intr: Interaction::Surface(SurfaceInteraction::new_simple(pi, n, uv)),
             pdf: 1. / self.area(),
         })
     }

shared/src/spectra/colorspace.rs

@@ -8,7 +8,7 @@ use crate::utils::ptr::Ptr;
 use std::cmp::{Eq, PartialEq};
 
 #[repr(C)]
-#[derive(Copy, Clone)]
+#[derive(Copy, Debug, Clone)]
 pub struct StandardColorSpaces {
     pub srgb: Ptr<RGBColorSpace>,
     pub dci_p3: Ptr<RGBColorSpace>,

shared/src/spectra/sampled.rs

@@ -1,5 +1,5 @@
 use crate::core::pbrt::Float;
-use crate::core::spectrum::StandardSpectra;
+use crate::core::spectrum::{SpectrumTrait, StandardSpectra};
 use crate::utils::math::{clamp, lerp};
 use std::ops::{
     Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
@@ -118,7 +118,7 @@ impl SampledSpectrum {
     }
 
     pub fn y(&self, lambda: &SampledWavelengths, std: &StandardSpectra) -> Float {
-        let ys = std.cie_y().sample(lambda);
+        let ys = std.y.sample(lambda);
         let pdf = lambda.pdf();
         SampledSpectrum::safe_div(&(ys * *self), &pdf).average() / CIE_Y_INTEGRAL
     }
@@ -305,7 +305,7 @@ pub struct SampledWavelengths {
 
 impl SampledWavelengths {
     pub fn pdf(&self) -> SampledSpectrum {
-        SampledSpectrum::from_vector(self.pdf.to_vec())
+        SampledSpectrum::from_array(&self.pdf)
     }
 
     pub fn secondary_terminated(&self) -> bool {

shared/src/spectra/simple.rs

@@ -1,8 +1,9 @@
 use super::cie::*;
 use super::sampled::{LAMBDA_MAX, LAMBDA_MIN};
-use crate::Float;
 use crate::core::spectrum::{Spectrum, SpectrumTrait};
 use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths};
+use crate::utils::ptr::Ptr;
+use crate::{Float, find_interval};
 use core::slice;
 use std::hash::{Hash, Hasher};
 use std::sync::LazyLock;
@@ -34,7 +35,7 @@ impl SpectrumTrait for ConstantSpectrum {
 pub struct DenselySampledSpectrum {
     pub lambda_min: i32,
     pub lambda_max: i32,
-    pub values: *const Float,
+    pub values: Ptr<Float>,
 }
 
 unsafe impl Send for DenselySampledSpectrum {}
@@ -42,148 +43,155 @@ unsafe impl Sync for DenselySampledSpectrum {}
 
 impl DenselySampledSpectrum {
     #[inline(always)]
-    fn as_slice(&self) -> &[Float] {
+    pub fn count(&self) -> usize {
         if self.values.is_null() {
-            return &[];
+            0
-        }
-        let len = (self.lambda_max - self.lambda_min + 1).max(0) as usize;
-        unsafe { slice::from_raw_parts(self.values, len) }
-    }
-
-    pub fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
-        let mut s = SampledSpectrum::default();
-
-        for i in 0..N_SPECTRUM_SAMPLES {
-            let offset = lambda[i].round() as i32 - self.lambda_min;
-            let len = (self.lambda_max - self.lambda_min + 1) as i32;
-
-            if offset < 0 || offset >= len {
-                s[i] = 0.0;
         } else {
-                unsafe { s[i] = *self.values.add(offset as usize) };
+            (self.lambda_max - self.lambda_min + 1) as usize
         }
     }
-        s
-    }
 
-    pub fn min_component_value(&self) -> Float {
+    #[inline(always)]
-        self.as_slice()
+    fn get(&self, idx: u32) -> Float {
-            .iter()
+        unsafe { *self.values.0.add(idx as usize) }
-            .fold(Float::INFINITY, |a, &b| a.min(b))
-    }
-
-    pub fn max_component_value(&self) -> Float {
-        self.as_slice()
-            .iter()
-            .fold(Float::NEG_INFINITY, |a, &b| a.max(b))
-    }
-
-    pub fn average(&self) -> Float {
-        let slice = self.as_slice();
-        if slice.is_empty() {
-            return 0.0;
-        }
-        slice.iter().sum::<Float>() / (slice.len() as Float)
-    }
-
-    pub fn safe_div(&self, rhs: SampledSpectrum) -> Self {
-        let mut r = Self::new(1, 1);
-        for i in 0..N_SPECTRUM_SAMPLES {
-            r.values[i] = if rhs[i] != 0.0 {
-                self.values[i] / rhs.values[i]
-            } else {
-                0.0
-            }
-        }
-        r
     }
 }
 
 impl PartialEq for DenselySampledSpectrum {
     fn eq(&self, other: &Self) -> bool {
-        if self.lambda_min != other.lambda_min
+        self.lambda_min == other.lambda_min
-            || self.lambda_max != other.lambda_max
+            && self.lambda_max == other.lambda_max
-            || self.values.len() != other.values.len()
+            && self.values.0 == other.values.0
-        {
-            return false;
-        }
-
-        self.values
-            .iter()
-            .zip(&other.values)
-            .all(|(a, b)| a.to_bits() == b.to_bits())
     }
 }
 
 impl Eq for DenselySampledSpectrum {}
 
-impl Hash for DenselySampledSpectrum {
+// impl Hash for DenselySampledSpectrum {
-    fn hash<H: Hasher>(&self, state: &mut H) {
+//     fn hash<H: Hasher>(&self, state: &mut H) {
-        self.lambda_min.hash(state);
+//         self.lambda_min.hash(state);
-        self.lambda_max.hash(state);
+//         self.lambda_max.hash(state);
-
+//
-        for v in &self.values {
+//         for v in &self.values {
-            v.to_bits().hash(state);
+//             v.to_bits().hash(state);
-        }
+//         }
-    }
+//     }
-}
+// }
 
 impl SpectrumTrait for DenselySampledSpectrum {
+    fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum {
+        let mut s = SampledSpectrum::default();
+        let n = self.count() as i32;
+
+        for i in 0..N_SPECTRUM_SAMPLES {
+            let offset = lambda[i].round() as i32 - self.lambda_min;
+
+            if offset < 0 || offset >= n {
+                s[i] = 0.0;
+            } else {
+                unsafe {
+                    s[i] = *self.values.0.add(offset as usize);
+                }
+            }
+        }
+        s
+    }
+
     fn evaluate(&self, lambda: Float) -> Float {
         let offset = (lambda.round() as i32) - self.lambda_min;
-        if offset < 0 || offset as usize >= self.values.len() {
+        let n = self.count() as i32;
+        if offset < 0 || offset >= n {
             0.0
         } else {
-            self.values[offset as usize]
+            unsafe { *self.values.0.add(offset as usize) }
         }
     }
 
     fn max_value(&self) -> Float {
-        self.values.iter().fold(Float::MIN, |a, b| a.max(*b))
+        if self.values.is_null() {
+            return 0.;
+        }
+
+        let n = self.count();
+        let mut max_val = Float::NEG_INFINITY;
+
+        for i in 0..n {
+            unsafe {
+                let val = *self.values.0.add(i);
+                if val > max_val {
+                    max_val = val;
+                }
+            }
+        }
+        max_val
     }
 }
 
 #[repr(C)]
 #[derive(Debug, Clone, Copy)]
 pub struct PiecewiseLinearSpectrum {
-    pub lambdas: *const Float,
+    pub lambdas: Ptr<Float>,
-    pub values: *const Float,
+    pub values: Ptr<Float>,
     pub count: u32,
 }
 
+impl PiecewiseLinearSpectrum {
+    #[inline(always)]
+    fn lambda(&self, i: u32) -> Float {
+        unsafe { *self.lambdas.0.add(i as usize) }
+    }
+
+    #[inline(always)]
+    fn value(&self, i: u32) -> Float {
+        unsafe { *self.values.0.add(i as usize) }
+    }
+}
+
 unsafe impl Send for PiecewiseLinearSpectrum {}
 unsafe impl Sync for PiecewiseLinearSpectrum {}
 
 impl SpectrumTrait for PiecewiseLinearSpectrum {
     fn evaluate(&self, lambda: Float) -> Float {
-        if self.lambdas.is_empty() {
+        if self.lambdas.is_null() {
             return 0.0;
         }
 
-        if lambda <= self.lambdas[0] {
+        if lambda <= self.lambda(0) {
-            return self.values[0];
+            return self.value(0);
         }
 
-        if lambda >= *self.lambdas.last().unwrap() {
+        if lambda >= self.lambda(self.count - 1) {
-            return *self.values.last().unwrap();
+            return self.value(self.count - 1);
         }
 
-        let i = self.lambdas.partition_point(|&l| l < lambda);
+        let i = find_interval(self.count, |idx| self.lambda(idx) <= lambda);
-        let l0 = self.lambdas[i - 1];
+
-        let l1 = self.lambdas[i];
+        let l0 = self.lambda(i);
-        let v0 = self.values[i - 1];
+        let l1 = self.lambda(i + 1);
-        let v1 = self.values[i];
+        let v0 = self.value(i);
+        let v1 = self.value(i + 1);
 
         let t = (lambda - l0) / (l1 - l0);
-
         v0 + t * (v1 - v0)
     }
 
     fn max_value(&self) -> Float {
-        if self.values.is_empty() {
+        if self.values.is_null() {
-            return 0.0;
+            return 0.;
         }
-        self.values.iter().fold(0.0, |acc, &v| acc.max(v))
+
+        let n = self.count;
+        let mut max_val = Float::NEG_INFINITY;
+
+        for i in 0..n {
+            unsafe {
+                let val = *self.values.0.add(i as usize);
+                if val > max_val {
+                    max_val = val;
+                }
+            }
+        }
+        max_val
     }
 }
 

shared/src/textures/checkerboard.rs

@@ -4,7 +4,7 @@ use crate::core::texture::{
     TextureMapping3DTrait,
 };
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::{Ptr, RelPtr, math::square};
+use crate::utils::{ArenaPtr, Ptr, math::square};
 
 fn checkerboard(
     ctx: &TextureEvalContext,
@@ -45,7 +45,7 @@ fn checkerboard(
 pub struct FloatCheckerboardTexture {
     pub map2d: Ptr<TextureMapping2D>,
     pub map3d: Ptr<TextureMapping3D>,
-    pub tex: [RelPtr<GPUFloatTexture>; 2],
+    pub tex: [ArenaPtr<GPUFloatTexture>; 2],
 }
 
 impl FloatCheckerboardTexture {
@@ -76,7 +76,7 @@ impl FloatCheckerboardTexture {
 pub struct SpectrumCheckerboardTexture {
     pub map2d: Ptr<TextureMapping2D>,
     pub map3d: Ptr<TextureMapping3D>,
-    pub tex: [RelPtr<GPUSpectrumTexture>; 2],
+    pub tex: [ArenaPtr<GPUSpectrumTexture>; 2],
 }
 
 impl SpectrumCheckerboardTexture {

shared/src/textures/dots.rs

@@ -4,7 +4,7 @@ use crate::core::texture::{
     GPUFloatTexture, GPUSpectrumTexture, TextureEvalContext, TextureMapping2D,
 };
 use crate::spectra::sampled::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::Ptr;
 use crate::utils::math::square;
 use crate::utils::noise::noise_2d;
 
@@ -28,21 +28,23 @@ fn inside_polka_dot(st: Point2f) -> bool {
 #[derive(Debug, Clone, Copy)]
 pub struct FloatDotsTexture {
     pub mapping: TextureMapping2D,
-    pub outside_dot: RelPtr<GPUFloatTexture>,
+    pub outside_dot: Ptr<GPUFloatTexture>,
-    pub inside_dot: RelPtr<GPUFloatTexture>,
+    pub inside_dot: Ptr<GPUFloatTexture>,
 }
 
 impl FloatDotsTexture {
     pub fn evaluate(&self, ctx: &TextureEvalContext) -> Float {
         let c = self.mapping.map(ctx);
-        if inside_polka_dot(c.st) {
+        let target_texture = if inside_polka_dot(c.st) {
-            if let Some(tex) = self.inside_dot.get() {
+            self.inside_dot
-                tex.evaluate(ctx)
-            }
         } else {
-            if let Some(tex) = self.outside_dot.get() {
+            self.outside_dot
-                tex.evaluate(ctx)
+        };
-            }
+
+        if !target_texture.is_null() {
+            target_texture.evaluate(ctx)
+        } else {
+            0.0
         }
     }
 }
@@ -51,8 +53,8 @@ impl FloatDotsTexture {
 #[derive(Clone, Copy, Debug)]
 pub struct SpectrumDotsTexture {
     pub mapping: TextureMapping2D,
-    pub outside_dot: RelPtr<GPUSpectrumTexture>,
+    pub outside_dot: Ptr<GPUSpectrumTexture>,
-    pub inside_dot: RelPtr<GPUSpectrumTexture>,
+    pub inside_dot: Ptr<GPUSpectrumTexture>,
 }
 
 impl SpectrumDotsTexture {
@@ -62,14 +64,17 @@ impl SpectrumDotsTexture {
         lambda: &SampledWavelengths,
     ) -> SampledSpectrum {
         let c = self.mapping.map(ctx);
-        if inside_polka_dot(c.st) {
+
-            if let Some(tex) = self.inside_dot.get() {
+        let target_texture = if inside_polka_dot(c.st) {
-                tex.evaluate(ctx, &lambda)
+            self.inside_dot
-            }
         } else {
-            if let Some(tex) = self.outside_dot.get() {
+            self.outside_dot
-                tex.evaluate(ctx, &lambda)
+        };
-            }
+
+        if !target_texture.is_null() {
+            target_texture.evaluate(ctx, lambda)
+        } else {
+            SampledSpectrum::new(0.0)
         }
     }
 }

shared/src/textures/marble.rs

@@ -6,6 +6,7 @@ use crate::core::texture::{TextureEvalContext, TextureMapping3D};
 use crate::spectra::{RGBAlbedoSpectrum, RGBColorSpace, SampledSpectrum, SampledWavelengths};
 use crate::utils::math::clamp;
 use crate::utils::noise::fbm;
+use crate::utils::ptr::Ptr;
 use crate::utils::splines::evaluate_cubic_bezier;
 
 #[repr(C)]
@@ -17,7 +18,7 @@ pub struct MarbleTexture {
     pub scale: Float,
     pub variation: Float,
     // TODO: DO not forget to pass StandardColorSpace here!!
-    pub colorspace: *const RGBColorSpace,
+    pub colorspace: Ptr<RGBColorSpace>,
 }
 
 unsafe impl Send for MarbleTexture {}
@@ -31,19 +32,26 @@ impl MarbleTexture {
     ) -> SampledSpectrum {
         let mut c = self.mapping.map(ctx);
         c.p = Point3f::from(Vector3f::from(c.p) * self.scale);
-        let marble =
+        let marble = c.p.y()
-            c.p.y() + self.variation * fbm(c.p, self.scale, c.dpdy, self.omega, self.octaves);
+            + self.variation
+                * fbm(
+                    c.p,
+                    self.scale * c.dpdx,
+                    self.scale * c.dpdy,
+                    self.omega,
+                    self.octaves,
+                );
         let t = 0.5 + 0.5 * marble.sin();
-        let colors: [RGB; 9] = [
+        let colors: [Point3f; 9] = [
-            RGB::new(0.58, 0.58, 0.6),
+            Point3f::new(0.58, 0.58, 0.6),
-            RGB::new(0.58, 0.58, 0.6),
+            Point3f::new(0.58, 0.58, 0.6),
-            RGB::new(0.58, 0.58, 0.6),
+            Point3f::new(0.58, 0.58, 0.6),
-            RGB::new(0.5, 0.5, 0.5),
+            Point3f::new(0.5, 0.5, 0.5),
-            RGB::new(0.6, 0.59, 0.58),
+            Point3f::new(0.6, 0.59, 0.58),
-            RGB::new(0.58, 0.58, 0.6),
+            Point3f::new(0.58, 0.58, 0.6),
-            RGB::new(0.58, 0.58, 0.6),
+            Point3f::new(0.58, 0.58, 0.6),
-            RGB::new(0.2, 0.2, 0.33),
+            Point3f::new(0.2, 0.2, 0.33),
-            RGB::new(0.58, 0.58, 0.6),
+            Point3f::new(0.58, 0.58, 0.6),
         ];
 
         const N_SEG: i32 = 6; // (9 - 3)
@@ -51,9 +59,9 @@ impl MarbleTexture {
         let first = ((t_clamped * N_SEG as Float).floor() as i32).clamp(0, N_SEG - 1);
         let t_segment = t_clamped * N_SEG as Float - first as Float;
         let first_idx = first as usize;
-        let rgb = evaluate_cubic_bezier(&colors[first_idx..first_idx + 4], t_segment) * 1.5;
+        let (rgb_vec, _) = evaluate_cubic_bezier(&colors[first_idx..first_idx + 4], t_segment);
+        let rgb = RGB::new(rgb_vec.x() * 1.5, rgb_vec.y() * 1.5, rgb_vec.z() * 1.5);
 
-        let color_space = unsafe { &*self.colorspace };
+        RGBAlbedoSpectrum::new(&*self.colorspace, rgb).sample(lambda)
-        RGBAlbedoSpectrum::new(color_space, rgb).sample(lambda)
     }
 }

shared/src/textures/mix.rs

@@ -2,14 +2,14 @@ use crate::Float;
 use crate::core::geometry::{Vector3f, VectorLike};
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvalContext};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::ArenaPtr;
 
 #[repr(C)]
 #[derive(Copy, Clone, Debug)]
 pub struct GPUFloatMixTexture {
-    pub tex1: RelPtr<GPUFloatTexture>,
+    pub tex1: ArenaPtr<GPUFloatTexture>,
-    pub tex2: RelPtr<GPUFloatTexture>,
+    pub tex2: ArenaPtr<GPUFloatTexture>,
-    pub amount: RelPtr<GPUFloatTexture>,
+    pub amount: ArenaPtr<GPUFloatTexture>,
 }
 
 impl GPUFloatMixTexture {
@@ -34,8 +34,8 @@ impl GPUFloatMixTexture {
 #[repr(C)]
 #[derive(Copy, Clone, Debug)]
 pub struct GPUFloatDirectionMixTexture {
-    pub tex1: RelPtr<GPUFloatTexture>,
+    pub tex1: ArenaPtr<GPUFloatTexture>,
-    pub tex2: RelPtr<GPUFloatTexture>,
+    pub tex2: ArenaPtr<GPUFloatTexture>,
     pub dir: Vector3f,
 }
 
@@ -61,9 +61,9 @@ impl GPUFloatDirectionMixTexture {
 #[repr(C)]
 #[derive(Copy, Clone, Debug)]
 pub struct GPUSpectrumMixTexture {
-    pub tex1: RelPtr<GPUSpectrumTexture>,
+    pub tex1: ArenaPtr<GPUSpectrumTexture>,
-    pub tex2: RelPtr<GPUSpectrumTexture>,
+    pub tex2: ArenaPtr<GPUSpectrumTexture>,
-    pub amount: RelPtr<GPUFloatTexture>,
+    pub amount: ArenaPtr<GPUFloatTexture>,
 }
 
 impl GPUSpectrumMixTexture {
@@ -98,8 +98,8 @@ impl GPUSpectrumMixTexture {
 #[repr(C)]
 #[derive(Copy, Clone, Debug)]
 pub struct GPUSpectrumDirectionMixTexture {
-    pub tex1: RelPtr<GPUSpectrumTexture>,
+    pub tex1: ArenaPtr<GPUSpectrumTexture>,
-    pub tex2: RelPtr<GPUSpectrumTexture>,
+    pub tex2: ArenaPtr<GPUSpectrumTexture>,
     pub dir: Vector3f,
 }
 

shared/src/textures/ptex.rs

@@ -21,7 +21,7 @@ impl GPUFloatPtexTexture {
 }
 
 #[repr(C)]
-#[derive(Clone, Copy)]
+#[derive(Clone, Debug, Copy)]
 pub struct GPUSpectrumPtexTexture {
     pub face_values: Slice<RGB>,
     pub n_faces: u32,
@@ -35,10 +35,8 @@ impl GPUSpectrumPtexTexture {
         ctx: &TextureEvalContext,
         lambda: &SampledWavelengths,
     ) -> SampledSpectrum {
-        let index = ctx
+        let index = ctx.face_index.clamp(0, self.n_faces.saturating_sub(1));
-            .face_index
+        let rgb = self.face_values[index as usize];
-            .clamp(0, self.n_faces.saturating_sub(1) as usize);
-        let rgb = self.face_values[index];
         let s_rgb = self.colorspaces.srgb;
 
         match self.spectrum_type {

shared/src/textures/scaled.rs

@@ -1,13 +1,13 @@
 use crate::Float;
 use crate::core::texture::{GPUFloatTexture, GPUSpectrumTexture, TextureEvalContext};
 use crate::spectra::{SampledSpectrum, SampledWavelengths};
-use crate::utils::RelPtr;
+use crate::utils::ArenaPtr;
 
 #[repr(C)]
 #[derive(Debug, Clone, Copy)]
 pub struct GPUFloatScaledTexture {
-    pub tex: RelPtr<GPUFloatTexture>,
+    pub tex: ArenaPtr<GPUFloatTexture>,
-    pub scale: RelPtr<GPUFloatTexture>,
+    pub scale: ArenaPtr<GPUFloatTexture>,
 }
 
 impl GPUFloatScaledTexture {
@@ -23,8 +23,8 @@ impl GPUFloatScaledTexture {
 #[repr(C)]
 #[derive(Debug, Clone, Copy)]
 pub struct GPUSpectrumScaledTexture {
-    pub tex: RelPtr<GPUSpectrumTexture>,
+    pub tex: ArenaPtr<GPUSpectrumTexture>,
-    pub scale: RelPtr<GPUFloatTexture>,
+    pub scale: ArenaPtr<GPUFloatTexture>,
 }
 
 impl GPUSpectrumScaledTexture {

shared/src/utils/math.rs

@@ -5,6 +5,7 @@ use crate::core::pbrt::{Float, FloatBitOps, FloatBits, ONE_MINUS_EPSILON, PI, PI
 use crate::utils::hash::{hash_buffer, mix_bits};
 use crate::utils::sobol::{SOBOL_MATRICES_32, VDC_SOBOL_MATRICES, VDC_SOBOL_MATRICES_INV};
 
+use crate::utils::Ptr;
 use half::f16;
 use num_traits::{Float as NumFloat, Num, One, Signed, Zero};
 use std::error::Error;
@@ -341,7 +342,7 @@ pub fn wrap_equal_area_square(uv: &mut Point2f) -> Point2f {
     *uv
 }
 
-pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(usize, [Float; 4])> {
+pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(u32, [Float; 4])> {
     if nodes.len() < 4 {
         return None;
     }
@@ -398,7 +399,7 @@ pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(usize, [Float;
         weights[3] = 0.0;
     }
 
-    Some((offset, weights))
+    Some((offset as u32, weights))
 }
 
 pub fn equal_area_sphere_to_square(d: Vector3f) -> Point2f {
@@ -416,8 +417,7 @@ pub fn equal_area_sphere_to_square(d: Vector3f) -> Point2f {
     let t5 = 0.881770664775316294736387951347e-1;
     let t6 = 0.419038818029165735901852432784e-1;
     let t7 = -0.251390972343483509333252996350e-1;
-    let mut phi = evaluate_polynomial(b, &[t1, t2, t3, t4, t5, t6, t7])
+    let mut phi = evaluate_polynomial(b, &[t1, t2, t3, t4, t5, t6, t7]);
-        .expect("Could not evaluate polynomial");
 
     if x < y {
         phi = 1. - phi;
@@ -749,45 +749,14 @@ pub fn inverse_radical_inverse(mut inverse: u64, base: u64, n_digits: u64) -> u6
 
 // Digit scrambling
 #[repr(C)]
-#[derive(Default, Debug, Clone)]
+#[derive(Default, Debug, Copy, Clone)]
 pub struct DigitPermutation {
-    base: u32,
+    pub base: u32,
-    n_digits: u32,
+    pub n_digits: u32,
-    permutations: Vec<u16>,
+    pub permutations: Ptr<u16>,
 }
 
 impl DigitPermutation {
-    pub fn new(base: u32, seed: u64) -> Self {
-        let mut n_digits: u32 = 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 as usize * base as usize];
-
-        for digit_index in 0..n_digits {
-            let hash_input = [base as u64, digit_index as u64, seed];
-            let dseed = hash_buffer(&hash_input, 0);
-
-            for digit_value in 0..base {
-                let index = (digit_index * base + digit_value) as usize;
-
-                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;
@@ -795,15 +764,8 @@ impl DigitPermutation {
     }
 }
 
-pub fn compute_radical_inverse_permutations(seed: u64) -> Vec<DigitPermutation> {
+pub fn scrambled_radical_inverse(base_index: u32, mut a: u64, perm: &DigitPermutation) -> Float {
-    PRIMES
+    let base = PRIMES[base_index as usize] as u64;
-        .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 {
-    let base = PRIMES[base_index] as u64;
 
     let limit = (u64::MAX / base).saturating_sub(base);
 
@@ -828,8 +790,8 @@ pub fn scrambled_radical_inverse(base_index: usize, mut a: u64, perm: &DigitPerm
     (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 {
+pub fn owen_scrambled_radical_inverse(base_index: u32, mut a: u64, hash: u32) -> Float {
-    let base = PRIMES[base_index] as u64;
+    let base = PRIMES[base_index as usize] as u64;
 
     let limit = (u64::MAX / base).saturating_sub(base);
     let inv_base = 1.0 / (base as Float);
@@ -1030,8 +992,8 @@ impl<F: Fn(u32) -> u32> Scrambler for F {
     }
 }
 
-const N_SOBOL_DIMENSIONS: usize = 1024;
+const N_SOBOL_DIMENSIONS: u32 = 1024;
-const SOBOL_MATRIX_SIZE: usize = 52;
+const SOBOL_MATRIX_SIZE: u32 = 52;
 #[inline]
 pub fn sobol_sample<S: Scrambler>(mut a: u64, dimension: u32, randomizer: S) -> Float {
     debug_assert!(
@@ -1046,7 +1008,7 @@ pub fn sobol_sample<S: Scrambler>(mut a: u64, dimension: u32, randomizer: S) ->
 
     while a != 0 {
         if (a & 1) != 0 {
-            v ^= SOBOL_MATRICES_32[i];
+            v ^= SOBOL_MATRICES_32[i as usize];
         }
         a >>= 1;
         i += 1;

shared/src/utils/mod.rs

@@ -15,7 +15,7 @@ pub mod sobol;
 pub mod splines;
 pub mod transform;
 
-pub use ptr::{Ptr, RelPtr};
+pub use ptr::{ArenaPtr, Ptr};
 pub use transform::{AnimatedTransform, Transform, TransformGeneric};
 
 #[inline]

shared/src/utils/ptr.rs

@@ -3,41 +3,43 @@ use core::ops::Index;
 
 #[repr(C)]
 #[derive(Debug)]
-pub struct RelPtr<T: ?Sized> {
+pub struct ArenaPtr<T: ?Sized> {
     offset: i32,
-    _phantom: PhantomData<T>,
+    _marker: PhantomData<T>,
 }
 
-impl<T: ?Sized> Clone for RelPtr<T> {
+impl<T: ?Sized> Clone for ArenaPtr<T> {
     fn clone(&self) -> Self {
         *self
     }
 }
 
-impl<T: ?Sized> Copy for RelPtr<T> {}
+impl<T: ?Sized> Copy for ArenaPtr<T> {}
 
-impl<T> RelPtr<T> {
+impl<T> ArenaPtr<T> {
     pub fn null() -> Self {
         Self {
-            offset: 0,
+            offset: 0xFFFFFFFF,
-            _phantom: PhantomData,
+            _marker: PhantomData,
         }
     }
 
     pub fn is_null(&self) -> bool {
-        self.offset == 0
+        self.offset == 0xFFFFFFFF
     }
 
-    pub fn get(&self) -> Option<&T> {
+    #[inline(always)]
-        if self.offset == 0 {
+    pub unsafe fn as_ptr(&self, base: *const u8) -> *const T {
-            None
+        if self.is_null() {
+            core::ptr::null()
         } else {
-            unsafe {
+            unsafe { base.add(self.offset as usize) as *const T }
-                let base = self as *const _ as *const u8;
-                let target = base.offset(self.offset as isize) as *const T;
-                target.as_ref()
         }
     }
+
+    #[inline(always)]
+    pub unsafe fn as_ref<'a>(&self, base: *const u8) -> &'a T {
+        unsafe { &*self.as_ptr(base) }
     }
 
     #[cfg(not(target_os = "cuda"))]
@@ -54,7 +56,7 @@ impl<T> RelPtr<T> {
 
         Self {
             offset: diff as i32,
-            _phantom: PhantomData,
+            _marker: PhantomData,
         }
     }
 }
@@ -69,6 +71,15 @@ impl<T> Default for Ptr<T> {
     }
 }
 
+impl<T> PartialEq for Ptr<T> {
+    #[inline(always)]
+    fn eq(&self, other: &Self) -> bool {
+        self.0 == other.0
+    }
+}
+
+impl<T> Eq for Ptr<T> {}
+
 impl<T> Ptr<T> {
     pub fn null() -> Self {
         Self::default()
@@ -83,6 +94,17 @@ impl<T> Ptr<T> {
     pub fn as_ref(&self) -> Option<&T> {
         unsafe { self.0.as_ref() }
     }
+
+    /// UNSTABLE: Casts the const pointer to mutable and returns a mutable reference.
+    /// THIS IS VERY DANGEROUS
+    /// The underlying data is not currently borrowed or accessed by any other thread/kernel.
+    /// The memory is actually writable.
+    /// No other mutable references exist to this data.
+    #[inline(always)]
+    pub unsafe fn as_mut(&mut self) -> &mut T {
+        debug_assert!(!self.is_null());
+        unsafe { &mut *(self.0 as *mut T) }
+    }
 }
 
 unsafe impl<T: Sync> Send for Ptr<T> {}
@@ -109,6 +131,26 @@ impl<T> From<&mut T> for Ptr<T> {
     }
 }
 
+impl<T> Index<usize> for Ptr<T> {
+    type Output = T;
+
+    #[inline(always)]
+    fn index(&self, index: usize) -> &Self::Output {
+        // There is no bounds checking because we dont know the length.
+        // It is host responsbility to check bounds
+        unsafe { &*self.0.add(index) }
+    }
+}
+
+impl<T> Index<u32> for Ptr<T> {
+    type Output = T;
+
+    #[inline(always)]
+    fn index(&self, index: u32) -> &Self::Output {
+        unsafe { &*self.0.add(index as usize) }
+    }
+}
+
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct Slice<T> {

shared/src/utils/sampling.rs

@@ -476,7 +476,7 @@ pub fn sample_catmull_rom(
     assert_eq!(f.len(), big_f.len());
 
     u *= big_f.last().copied().unwrap_or(0.);
-    let i = find_interval(big_f.len(), |i| big_f[i] <= u);
+    let i = find_interval(big_f.len() as u32, |i| big_f[i as usize] <= u) as usize;
     let x0 = nodes[i];
     let x1 = nodes[i + 1];
     let f0 = f[i];
@@ -518,8 +518,8 @@ pub fn sample_catmull_rom(
     let mut big_fhat = 0.;
 
     let eval = |t: Float| -> (Float, Float) {
-        big_fhat = evaluate_polynomial(t, big_fhat_coeffs_ref).unwrap_or(0.);
+        big_fhat = evaluate_polynomial(t, big_fhat_coeffs_ref);
-        fhat = evaluate_polynomial(t, fhat_coeffs_ref).unwrap_or(0.);
+        fhat = evaluate_polynomial(t, fhat_coeffs_ref);
         (big_fhat - u, fhat)
     };
     let t = newton_bisection(0., 1., eval);
@@ -543,12 +543,13 @@ pub fn sample_catmull_rom_2d(
         None => return (0., 0., 0.),
     };
 
-    let n2 = nodes2.len();
+    let n2 = nodes2.len() as u32;
-    let interpolate = |array: &[Float], idx: usize| -> Float {
+    let interpolate = |array: &[Float], idx: u32| -> Float {
         let mut v = 0.;
         for i in 0..4 {
             if weights[i] != 0. {
-                v += array[(offset + i) * n2 + idx] * weights[i];
+                let ind = (offset + i as u32) * n2 + idx;
+                v += array[ind as usize] * weights[i];
             }
         }
         v
@@ -558,23 +559,25 @@ pub fn sample_catmull_rom_2d(
         return (0., 0., 0.);
     }
     u *= maximum;
-    let idx = find_interval(n2, |i| interpolate(cdf, i) <= u);
+    // TODO: Make find_interval(binary_search) integer agnostic, this is a PITA
-    let f0 = interpolate(values, idx);
+    let id = find_interval(n2 as u32, |i| interpolate(cdf, i) <= u);
-    let f1 = interpolate(values, idx + 1);
+    let f0 = interpolate(values, id);
+    let f1 = interpolate(values, id + 1);
+    let idx = id as usize;
     let x0 = nodes2[idx];
     let x1 = nodes2[idx + 1];
     let width = x1 - x0;
     let d0 = if idx > 0 {
-        width * (f1 - interpolate(values, idx - 1)) / (x1 - nodes2[idx - 1])
+        width * (f1 - interpolate(values, id - 1)) / (x1 - nodes2[idx - 1])
     } else {
         f1 - f0
     };
-    let d1 = if idx + 2 < n2 {
+    let d1 = if id + 2 < n2 {
-        width * (interpolate(values, idx + 2) - f0) / (nodes2[idx + 2] - x0)
+        width * (interpolate(values, id + 2) - f0) / (nodes2[idx + 2] - x0)
     } else {
         f1 - f0
     };
-    u = (u - interpolate(cdf, idx)) / width;
+    u = (u - interpolate(cdf, id)) / width;
 
     let fhat_coeffs = [
         f0,
@@ -597,8 +600,8 @@ pub fn sample_catmull_rom_2d(
     let mut fhat = 0.0;
 
     let eval = |t: Float| -> (Float, Float) {
-        big_fhat = evaluate_polynomial(t, big_fhat_coeffs_ref).unwrap_or(0.);
+        big_fhat = evaluate_polynomial(t, big_fhat_coeffs_ref);
-        fhat = evaluate_polynomial(t, fhat_coeffs_ref).unwrap_or(0.);
+        fhat = evaluate_polynomial(t, fhat_coeffs_ref);
         (big_fhat - u, fhat)
     };
 
@@ -700,67 +703,28 @@ pub struct PLSample {
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
 pub struct PiecewiseConstant1D {
-    pub func: *mut Float,
+    pub func: Ptr<Float>,
-    pub cdf: *mut Float,
+    pub cdf: Ptr<Float>,
     pub min: Float,
     pub max: Float,
-    pub n: usize,
+    pub n: u32,
     pub func_integral: Float,
 }
 
+unsafe impl Send for PiecewiseConstant1D {}
+unsafe impl Sync for PiecewiseConstant1D {}
+
 impl PiecewiseConstant1D {
-    #[cfg(not(target_os = "cuda"))]
-    pub fn new_with_bounds(f: &[Float], min: Float, max: Float) -> Self {
-        let n = f.len();
-        let mut func_vec = f.to_vec();
-        let mut cdf_vec = vec![0.0; n + 1];
-
-        cdf_vec[0] = 0.0;
-        for i in 1..=n {
-            cdf_vec[i] = cdf_vec[i - 1] + func_vec[i - 1] / n as Float;
-        }
-
-        let func_integral = cdf_vec[n];
-        if func_integral > 0.0 {
-            for i in 1..=n {
-                cdf_vec[i] /= func_integral;
-            }
-        } else {
-            for i in 1..=n {
-                cdf_vec[i] = i as Float / n as Float;
-            }
-        }
-
-        let func = func_vec.as_mut_ptr();
-        let cdf = cdf_vec.as_mut_ptr();
-        std::mem::forget(func_vec);
-        std::mem::forget(cdf_vec);
-
-        Self {
-            func,
-            cdf,
-            min,
-            max,
-            n,
-            func_integral,
-        }
-    }
-
-    #[cfg(not(target_os = "cuda"))]
-    pub fn new(f: &[Float]) -> Self {
-        Self::new_with_bounds(f, 0., 1.)
-    }
-
     pub fn integral(&self) -> Float {
         self.func_integral
     }
 
-    pub fn size(&self) -> usize {
+    pub fn size(&self) -> u32 {
         self.n
     }
 
-    pub fn sample(&self, u: Float) -> (Float, Float, usize) {
+    pub fn sample(&self, u: Float) -> (Float, Float, u32) {
-        let o = find_interval(self.cdf.len(), |idx| self.cdf[idx] <= u);
+        let o = find_interval(self.size(), |idx| self.cdf[idx] <= u) as usize;
         let mut du = u - self.cdf[o];
         if self.cdf[o + 1] - self.cdf[o] > 0. {
             du /= self.cdf[o + 1] - self.cdf[o];
@@ -772,7 +736,7 @@ impl PiecewiseConstant1D {
         } else {
             0.
         };
-        (value, pdf_val, o)
+        (value, pdf_val, o as u32)
     }
 }
 
@@ -780,62 +744,16 @@ impl PiecewiseConstant1D {
 #[derive(Debug, Copy, Clone)]
 pub struct PiecewiseConstant2D {
     pub domain: Bounds2f,
-    pub p_conditional_v: Ptr<PiecewiseConstant1D>,
     pub p_marginal: PiecewiseConstant1D,
     pub n_conditionals: usize,
+    pub p_conditional_v: Ptr<PiecewiseConstant1D>,
 }
 
 impl PiecewiseConstant2D {
-    #[cfg(not(target_os = "cuda"))]
-    pub fn new(data: &Array2D<Float>, x_size: u32, y_size: u32, domain: Bounds2f) -> Self {
-        let nu = x_size as usize;
-        let nv = y_size as usize;
-        let mut conditionals = Vec::with_capacity(nv);
-        for v in 0..nv {
-            let row = unsafe { core::slice::from_raw_parts(data.values.add(v * nu), nu) };
-            conditionals.push(PiecewiseConstant1D::new_with_bounds(
-                row,
-                domain.p_min.x(),
-                domain.p_max.x(),
-            ));
-        }
-
-        let marginal_funcs: Vec<Float> = conditionals.iter().map(|c| c.func_integral).collect();
-        let p_marginal = PiecewiseConstant1D::new_with_bounds(
-            &marginal_funcs,
-            domain.p_min.y(),
-            domain.p_max.y(),
-        );
-
-        let p_conditional_v = conditionals.as_mut_ptr();
-        std::mem::forget(conditionals);
-
-        Self {
-            p_conditional_v,
-            p_marginal,
-            domain,
-            n_conditionals: nv,
-        }
-    }
-
-    #[cfg(not(target_os = "cuda"))]
-    pub fn new_with_bounds(data: &Array2D<Float>, domain: Bounds2f) -> Self {
-        Self::new(data, data.x_size(), data.y_size(), domain)
-    }
-
-    #[cfg(not(target_os = "cuda"))]
-    pub fn new_with_data(data: &Array2D<Float>) -> Self {
-        let nx = data.x_size();
-        let ny = data.y_size();
-        let domain = Bounds2f::new(Point2f::new(0.0, 0.0), Point2f::new(1.0, 1.0));
-
-        Self::new(data, nx, ny, domain)
-    }
-
     pub fn resolution(&self) -> Point2i {
         Point2i::new(
-            self.p_conditional_v[0].size() as i32,
+            self.p_conditional_v[0u32].size() as i32,
-            self.p_conditional_v[1].size() as i32,
+            self.p_conditional_v[1u32].size() as i32,
         )
     }
 
@@ -853,7 +771,7 @@ impl PiecewiseConstant2D {
 
     pub fn pdf(&self, p: Point2f) -> f32 {
         let p_offset = self.domain.offset(&p);
-        let nu = self.p_conditional_v[0].size();
+        let nu = self.p_conditional_v[0u32].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;
@@ -875,35 +793,35 @@ pub struct SummedAreaTable {
 }
 
 impl SummedAreaTable {
-    pub fn new(values: &Array2D<Float>) -> Self {
+    // pub fn new(values: &Array2D<Float>) -> Self {
-        let width = values.x_size();
+    //     let width = values.x_size();
-        let height = values.y_size();
+    //     let height = values.y_size();
-
+    //
-        let mut sum = Array2D::<f64>::new_with_dims(width, height);
+    //     let mut sum = Array2D::<f64>::new_with_dims(width, height);
-        sum[(0, 0)] = values[(0, 0)] as f64;
+    //     sum[(0, 0)] = values[(0, 0)] as f64;
-
+    //
-        for x in 1..width {
+    //     for x in 1..width {
-            sum[(x, 0)] = values[(x, 0)] as f64 + sum[(x - 1, 0)];
+    //         sum[(x, 0)] = values[(x as i32, 0)] as f64 + sum[(x - 1, 0)];
-        }
+    //     }
-
+    //
-        for y in 1..height {
+    //     for y in 1..height {
-            sum[(0, y)] = values[(0, y)] as f64 + sum[(0, y - 1)];
+    //         sum[(0, y)] = values[(0, y as i32)] as f64 + sum[(0, y - 1)];
-        }
+    //     }
-
+    //
-        for y in 1..height {
+    //     for y in 1..height {
-            for x in 1..width {
+    //         for x in 1..width {
-                let term = values[(x, y)] as f64;
+    //             let term = values[(x as i32, y as i32)] as f64;
-                let left = sum[(x - 1, y)];
+    //             let left = sum[(x - 1, y)];
-                let up = sum[(x, y - 1)];
+    //             let up = sum[(x, y - 1)];
-                let diag = sum[(x - 1, y - 1)];
+    //             let diag = sum[(x - 1, y - 1)];
-
+    //
-                sum[(x, y)] = term + left + up - diag;
+    //             sum[(x, y)] = term + left + up - diag;
-            }
+    //         }
-        }
+    //     }
-
+    //
-        Self { sum }
+    //     Self { sum }
-    }
+    // }
-
+    //
     pub fn integral(&self, extent: Bounds2f) -> Float {
         let s = self.lookup(extent.p_max.x(), extent.p_max.y())
             - self.lookup(extent.p_min.x(), extent.p_max.y())
@@ -941,8 +859,8 @@ impl SummedAreaTable {
             return 0.0;
         }
 
-        let ix = (x - 1).min(self.sum.x_size() as i32 - 1) as usize;
+        let ix = (x - 1).min(self.sum.x_size() as i32 - 1);
-        let iy = (y - 1).min(self.sum.y_size() as i32 - 1) as usize;
+        let iy = (y - 1).min(self.sum.y_size() as i32 - 1);
 
         self.sum[(ix, iy)]
     }
@@ -956,10 +874,10 @@ pub struct WindowedPiecewiseConstant2D {
 }
 
 impl WindowedPiecewiseConstant2D {
-    pub fn new(func: Array2D<Float>) -> Self {
+    // pub fn new(func: Array2D<Float>) -> Self {
-        let sat = SummedAreaTable::new(&func);
+    //     let sat = SummedAreaTable::new(&func);
-        Self { sat, func }
+    //     Self { sat, func }
-    }
+    // }
 
     pub fn sample(&self, u: Point2f, b: Bounds2f) -> Option<(Point2f, Float)> {
         let b_int = self.sat.integral(b);
@@ -1023,13 +941,13 @@ impl WindowedPiecewiseConstant2D {
         let nx = self.func.x_size();
         let ny = self.func.y_size();
 
-        let ix = ((p.x() * nx as Float) as i32).min(nx as i32 - 1).max(0) as usize;
+        let ix = ((p.x() * nx as Float) as i32).min(nx as i32 - 1).max(0);
-        let iy = ((p.y() * ny as Float) as i32).min(ny as i32 - 1).max(0) as usize;
+        let iy = ((p.y() * ny as Float) as i32).min(ny as i32 - 1).max(0);
 
         self.func[(ix, iy)]
     }
 
-    fn sample_bisection<F>(p_func: F, u: Float, mut min: Float, mut max: Float, n: usize) -> Float
+    fn sample_bisection<F>(p_func: F, u: Float, mut min: Float, mut max: Float, n: u32) -> Float
     where
         F: Fn(Float) -> Float,
     {
@@ -1368,8 +1286,8 @@ impl<const N: usize> PiecewiseLinear2D<N> {
                 continue;
             }
 
-            let param_index = find_interval(size as usize, |idx| {
+            let param_index = find_interval(size, |idx| {
-                self.get_param_value(dim, idx) <= params[dim]
+                self.get_param_value(dim, idx as usize) <= params[dim]
             }) as u32;
 
             let p0 = self.get_param_value(dim, param_index as usize);

shared/src/utils/splines.rs

@@ -1,6 +1,7 @@
 use crate::core::geometry::{Bounds3f, Lerp, Point3f, Vector3f, VectorLike};
 use crate::core::pbrt::Float;
 use crate::utils::math::lerp;
+use core::ops::Sub;
 use num_traits::Num;
 
 fn bounds_cubic_bezier(cp: &[Point3f]) -> Bounds3f {
@@ -39,7 +40,7 @@ where
 }
 
 pub fn subdivide_cubic_bezier(cp: &[Point3f]) -> [Point3f; 7] {
-    let v: Vec<Vector3f> = cp.iter().map(|&p| p.into()).collect();
+    let v: [Vector3f; 4] = core::array::from_fn(|i| Vector3f::from(cp[i]));
     let v01 = (v[0] + v[1]) / 2.0;
     let v12 = (v[1] + v[2]) / 2.0;
     let v23 = (v[2] + v[3]) / 2.0;
@@ -99,7 +100,9 @@ pub fn quadratic_bspline_to_bezier(cp: &[Point3f]) -> [Point3f; 3] {
     [p11, cp[1], p22]
 }
 
-pub fn evaluate_cubic_bezier(cp: &[Point3f], u: Float) -> (Point3f, Vector3f) {
+pub fn evaluate_cubic_bezier(cp: &[Point3f], u: Float) -> (Point3f, Vector3f)
+where
+{
     let cp1 = [
         lerp(u, cp[0], cp[1]),
         lerp(u, cp[1], cp[2]),

shared/src/utils/transform.rs

@@ -1,5 +1,4 @@
 use num_traits::Float as NumFloat;
-use std::error::Error;
 use std::iter::{Product, Sum};
 use std::ops::{Add, Div, Index, IndexMut, Mul};
 use std::sync::Arc;
@@ -167,7 +166,7 @@ impl TransformGeneric<Float> {
                 *t -= dt;
             }
         }
-        Ray::new(o.into(), r.d, Some(r.time), r.medium.clone())
+        Ray::new(o.into(), r.d, Some(r.time), &*r.medium)
     }
 
     pub fn apply_to_interval(&self, pi: &Point3fi) -> Point3fi {
@@ -262,7 +261,7 @@ impl TransformGeneric<Float> {
                 ret.shading.dndu = self.apply_to_normal(si.shading.dndu);
                 ret.shading.dndv = self.apply_to_normal(si.shading.dndv);
 
-                ret.common.n = n.normalize().face_forward(ret.shading.n.into());
+                ret.common.n = n.normalize().face_forward(ret.shading.n);
 
                 Interaction::Surface(ret)
             }
@@ -368,10 +367,7 @@ impl TransformGeneric<Float> {
                 t = t_max - dt;
             }
         }
-        (
+        (Ray::new(Point3f::from(o), d, Some(r.time), &*r.medium), t)
-            Ray::new(Point3f::from(o), d, Some(r.time), r.medium.clone()),
-            t,
-        )
     }
 
     pub fn to_quaternion(self) -> Quaternion {
@@ -819,7 +815,7 @@ impl AnimatedTransform {
                 actually_animated: false,
                 t: [Vector3f::default(); 2],
                 r: [Quaternion::default(); 2],
-                s: std::array::from_fn(|_| SquareMatrix::default()),
+                s: core::array::from_fn(|_| SquareMatrix::default()),
                 has_rotation: false,
                 c1: [DerivativeTerm::default(); 3],
                 c2: [DerivativeTerm::default(); 3],

src/core/light.rs

@@ -1,6 +1,11 @@
-use shared::core::light::LIghtBase;
+use shared::core::geometry::{Bounds3f, Point2i};
+use shared::core::image::Image;
+use shared::core::light::LightBase;
+use shared::core::medium::MediumInterface;
 use shared::core::spectrum::Spectrum;
-use shared::spectra::DenselySampledSpectrum;
+use shared::spectra::{DenselySampledSpectrum, SampledSpectrum, SampledWavelengths};
+use shared::utils::Transform;
+use shared::{Float, PI};
 
 use crate::core::spectrum::SPECTRUM_CACHE;
 use crate::utils::containers::InternCache;
@@ -14,3 +19,13 @@ pub trait LightBaseTrait {
 }
 
 impl LightBaseTrait for LightBase {}
+
+pub trait LightFactory {
+    fn new(
+        render_from_light: Transform,
+        medium_interface: MediumInterface,
+        scale: Float,
+        iemit: &Spectrum,
+        image: &Image,
+    ) -> Self;
+}

src/lights/goniometric.rs

@@ -0,0 +1,28 @@
+use crate::core::light::{LightBaseTrait, LightFactory};
+
+impl LightFactory for GoniometricLight {
+    fn new(
+        render_from_light: Transform,
+        medium_interface: MediumInterface,
+        scale: Float,
+        iemit: &Spectrum,
+        image: &Image,
+    ) -> Self {
+        let base = LightBase::new(
+            LightType::DeltaPosition,
+            render_from_light,
+            medium_interface,
+        );
+
+        let i_interned = LightBase::lookup_spectrum(&iemit);
+        let d = image.get_sampling_distribution_uniform();
+        let distrib = PiecewiseConstant2D::new_with_data(&d);
+        Self {
+            base,
+            iemit: i_interned,
+            scale,
+            image: Ptr::from(image),
+            distrib,
+        }
+    }
+}

src/lights/infinite.rs

@@ -8,6 +8,8 @@ use shared::utils::Transform;
 use shared::utils::sampling::PiecewiseConstant2D;
 use std::sync::Arc;
 
+use crate::core::light::{LightBaseTrait, LightFactory};
+
 #[derive(Debug)]
 struct InfiniteImageLightStorage {
     image: Image,
@@ -267,3 +269,21 @@ impl InfinitePortalLightHost {
         }
     }
 }
+
+impl LightFactory for InfiniteUniformLight {
+    fn new(render_from_light: Transform, le: Spectrum, scale: Float) -> Self {
+        let base = LightBase::new(
+            LightType::Infinite,
+            render_from_light,
+            MediumInterface::default(),
+        );
+        let lemit = LightBase::lookup_spectrum(&le);
+        Self {
+            base,
+            lemit,
+            scale,
+            scene_center: Point3f::default(),
+            scene_radius: 0.,
+        }
+    }
+}

src/lights/mod.rs

@@ -1,4 +1,5 @@
 pub mod diffuse;
+pub mod goniometric;
 pub mod infinite;
 pub mod projection;
 pub mod sampler;

src/utils/math.rs

@@ -0,0 +1,67 @@
+use half::f16;
+use shared::Float;
+use shared::utils::Ptr;
+use shared::utils::math::{DigitPermutation, PRIMES};
+
+pub fn new_digit_permutation(base: u32, seed: u64) -> Vec<u16> {
+    let mut n_digits: u32 = 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 as usize * base as usize];
+
+    for digit_index in 0..n_digits {
+        let hash_input = [base as u64, digit_index as u64, seed];
+        let dseed = hash_buffer(&hash_input, 0);
+
+        for digit_value in 0..base {
+            let index = (digit_index * base + digit_value) as usize;
+
+            permutations[index] =
+                permutation_element(digit_value as u32, base as u32, dseed as u32) as u16;
+        }
+    }
+
+    permutations
+}
+
+pub fn compute_radical_inverse_permutations(seed: u64) -> (Vec<u16>, Vec<DigitPermutation>) {
+    let temp_data: Vec<Vec<u16>> = PRIMES
+        .iter()
+        .map(|&base| new_digit_permutation(base as u32, seed))
+        .collect();
+    let mut storage: Vec<u16> = Vec::with_capacity(temp_data.iter().map(|v| v.len()).sum());
+
+    for vec in &temp_data {
+        storage.extend_from_slice(vec);
+    }
+
+    let mut views = Vec::with_capacity(PRIMES.len());
+    let mut current_offset = 0;
+
+    let storage_base_ptr = storage.as_ptr();
+
+    for (i, &base) in PRIMES.iter().enumerate() {
+        let len = temp_data[i].len();
+        let n_digits = len as u32 / base as u32;
+
+        unsafe {
+            let ptr_to_data = storage_base_ptr.add(current_offset);
+
+            views.push(DigitPermutation::new(
+                base as u32,
+                n_digits,
+                Ptr(ptr_to_data),
+            ));
+        }
+
+        current_offset += len;
+    }
+
+    (storage, views)
+}

src/utils/mod.rs

@@ -2,6 +2,7 @@ pub mod containers;
 pub mod error;
 pub mod file;
 pub mod io;
+pub mod math;
 pub mod mipmap;
 pub mod parallel;
 pub mod parameters;

src/utils/sampling.rs

@@ -1,7 +1,131 @@
 use shared::Float;
-use shared::utils::sampling::{AliasTable, PiecewiseLinear2D};
+use shared::core::geometry::{Bounds2f, Point2f};
+use shared::utils::Ptr;
+use shared::utils::sampling::{
+    AliasTable, PiecewiseConstant1D, PiecewiseConstant2D, PiecewiseLinear2D,
+};
 use std::sync::Arc;
 
+#[derive(Debug, Clone)]
+pub struct PiecewiseConstant1DHost {
+    pub view: PiecewiseConstant1D,
+    _func: Vec<Float>,
+    _cdf: Vec<Float>,
+}
+
+impl std::ops::Deref for PiecewiseConstant1DHost {
+    type Target = PiecewiseConstant1D;
+    fn deref(&self) -> &Self::Target {
+        &self.view
+    }
+}
+
+impl PiecewiseConstant1DHost {
+    pub fn new(f: &[Float]) -> Self {
+        Self::new_with_bounds(f, 0.0, 1.0)
+    }
+
+    pub fn new_with_bounds(f: &[Float], min: Float, max: Float) -> Self {
+        let n = f.len();
+
+        let mut func_vec = f.to_vec();
+        let mut cdf_vec = vec![0.0; n + 1];
+
+        cdf_vec[0] = 0.0;
+        for i in 1..=n {
+            cdf_vec[i] = cdf_vec[i - 1] + func_vec[i - 1] / n as Float;
+        }
+
+        let func_integral = cdf_vec[n];
+        if func_integral > 0.0 {
+            for i in 1..=n {
+                cdf_vec[i] /= func_integral;
+            }
+        } else {
+            for i in 1..=n {
+                cdf_vec[i] = i as Float / n as Float;
+            }
+        }
+
+        let view = PiecewiseConstant1D {
+            func: Ptr(func_vec.as_ptr()),
+            cdf: Ptr(cdf_vec.as_ptr()),
+            min,
+            max,
+            n: n as u32,
+            func_integral,
+        };
+
+        Self {
+            view,
+            _func: func_vec,
+            _cdf: cdf_vec,
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct PiecewiseConstant2DHost {
+    pub view: PiecewiseConstant2D,
+    _p_conditional_v: Vec<PiecewiseConstant1D>,
+}
+
+impl std::ops::Deref for PiecewiseConstant2DHost {
+    type Target = PiecewiseConstant2D;
+    fn deref(&self) -> &Self::Target {
+        &self.view
+    }
+}
+
+impl PiecewiseConstant2DHost {
+    pub fn new(data: &Array2D<Float>, x_size: u32, y_size: u32, domain: Bounds2f) -> Self {
+        let nu = x_size as usize;
+        let nv = y_size as usize;
+        let mut conditionals = Vec::with_capacity(nv);
+        for v in 0..nv {
+            let row = unsafe { core::slice::from_raw_parts(data.values.add(v * nu), nu) };
+            conditionals.push(PiecewiseConstant1D::new_with_bounds(
+                row,
+                domain.p_min.x(),
+                domain.p_max.x(),
+            ));
+        }
+
+        let marginal_funcs: Vec<Float> = conditionals.iter().map(|c| c.func_integral).collect();
+        let p_marginal = PiecewiseConstant1D::new_with_bounds(
+            &marginal_funcs,
+            domain.p_min.y(),
+            domain.p_max.y(),
+        );
+
+        let p_conditional_v = conditionals.as_mut_ptr();
+        std::mem::forget(conditionals);
+        let view = PiecewiseConstant2D {
+            domain,
+            p_marginal,
+            n_conditionals: nv,
+            p_conditional_v: Ptr(p_conditional_v),
+        };
+
+        Self {
+            view,
+            _p_conditional_v: p_conditional_v,
+        }
+    }
+
+    pub fn new_with_bounds(data: &Array2D<Float>, domain: Bounds2f) -> Self {
+        Self::new(data, data.x_size(), data.y_size(), domain)
+    }
+
+    pub fn new_with_data(data: &Array2D<Float>) -> Self {
+        let nx = data.x_size();
+        let ny = data.y_size();
+        let domain = Bounds2f::from_points(Point2f::new(0.0, 0.0), Point2f::new(1.0, 1.0));
+
+        Self::new(data, nx, ny, domain)
+    }
+}
+
 struct PiecewiseLinear2DStorage<const D: usize> {
     data: Vec<Float>,
     marginal_cdf: Vec<Float>,