Just to be safe

src/core/bxdf.rs

@@ -7,7 +7,9 @@ use std::sync::{Arc, RwLock};
 
 use enum_dispatch::enum_dispatch;
 
-use crate::core::pbrt::{Float, INV_2_PI, INV_PI, PI, PI_OVER_2};
+use crate::core::medium::{HGPhaseFunction, PhaseFunctionTrait};
+use crate::core::options::get_options;
+use crate::core::pbrt::{Float, INV_2_PI, INV_PI, ONE_MINUS_EPSILON, PI, PI_OVER_2, clamp_t};
 use crate::geometry::{
     Frame, Normal3f, Point2f, Vector3f, VectorLike, abs_cos_theta, cos_theta, same_hemisphere,
     spherical_direction, spherical_theta,
@@ -17,13 +19,15 @@ use crate::spectra::{
 };
 use crate::utils::color::RGB;
 use crate::utils::colorspace::RGBColorSpace;
+use crate::utils::hash::hash_buffer;
 use crate::utils::math::{
     fast_exp, i0, 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, sample_cosine_hemisphere, sample_trimmed_logistic,
+    PiecewiseLinear2D, cosine_hemisphere_pdf, power_heuristic, sample_cosine_hemisphere,
-    sample_uniform_hemisphere, uniform_hemisphere_pdf,
+    sample_exponential, sample_trimmed_logistic, sample_uniform_hemisphere, uniform_hemisphere_pdf,
 };
 use crate::utils::scattering::{
     TrowbridgeReitzDistribution, fr_complex, fr_complex_from_spectrum, fr_dielectric, reflect,
@@ -395,6 +399,7 @@ impl ConductorBxDF {
 #[derive(Debug, Clone)]
 pub struct NormalizedFresnelBxDF;
 
+#[derive(Debug, Copy, Clone)]
 pub struct FArgs {
     pub mode: TransportMode,
     pub sample_flags: BxDFReflTransFlags,
@@ -1515,3 +1520,524 @@ impl BxDFTrait for HairBxDF {
         self
     }
 }
+
+#[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(),
+        }
+    }
+}
+
+#[derive(Clone, 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_t(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 mut r = || -> Float { rng.uniform::<Float>().min(ONE_MINUS_EPSILON) };
+        let f_args = FArgs {
+            mode,
+            sample_flags: BxDFReflTransFlags::TRANSMISSION,
+        };
+
+        let refl_args = FArgs {
+            mode,
+            sample_flags: BxDFReflTransFlags::REFLECTION,
+        };
+
+        for _ in 0..self.n_samples {
+            // Sample random walk through layers to estimate BSDF value
+            let mut uc = r();
+            let Some(wos) = enter_interface
+                .sample_f(wo, uc, Point2f::new(r(), r()), f_args)
+                .filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)
+            else {
+                continue;
+            };
+
+            uc = r();
+            let Some(wis) = exit_interface
+                .sample_f(wo, uc, Point2f::new(r(), r()), f_args)
+                .filter(|s| !s.f.is_black() && s.pdf > 0.0 && s.wi.z() != 0.0)
+            else {
+                continue;
+            };
+
+            // Declare state for random walk through BSDF layers
+            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 {
+                if depth > 3 && beta.max_component_value() < 0.25 {
+                    let q = (1. - beta.max_component_value()).max(0.);
+                    if r() < q {
+                        break;
+                    }
+                    beta /= 1. - q;
+                }
+                // Account for media between layers and possibly scatter
+                if !self.albedo.is_black() {
+                    z = if z == self.thickness {
+                        0.
+                    } else {
+                        self.thickness
+                    };
+                    beta *= self.tr(self.thickness, w);
+                } else {
+                    let sigma_t = 1.;
+                    let dz = sample_exponential(r(), sigma_t / w.z().abs());
+                    let zp = if w.z() > 0. { z + dz } else { z - dz };
+                    if 0. < zp && zp < self.thickness {
+                        // Handle scattering event in layered BSDF medium
+                        let mut wt = 1.;
+                        if exit_interface.flags().is_specular() {
+                            wt = power_heuristic(1, wis.pdf, 1, phase.pdf(-w, wis.wi));
+                        }
+                        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 u = Point2f::new(r(), r());
+                        let Some(ps) = phase
+                            .sample_p(-w, u)
+                            .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;
+
+                        if z < exit_z && w.z() > 0.
+                            || z > exit_z && w.z() < 0. && exit_interface.flags().is_specular()
+                        {
+                            let f_exit = exit_interface.f(-w, -wi, mode);
+                            if !f_exit.is_black() {
+                                let exit_pdf = exit_interface.pdf(-w, wi, f_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_t(zp, 0., self.thickness);
+                }
+
+                // Account for scattering at appropriate interface
+                if z == exit_z {
+                    let uc = r();
+                    let u = Point2f::new(r(), r());
+                    let Some(bs) = exit_interface
+                        .sample_f(-w, uc, u, 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;
+                } else {
+                    if !non_exit_interface.flags().is_specular() {
+                        let mut wt = 1.;
+                        if exit_interface.flags().is_specular() {
+                            wt = power_heuristic(
+                                1,
+                                wis.pdf,
+                                1,
+                                non_exit_interface.pdf(-w, -wis.wi, refl_args),
+                            );
+                        }
+                        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 using BSDF at _nonExitInterface_
+                    let uc = r();
+                    let u = Point2f::new(r(), r());
+                    let Some(bs) = non_exit_interface
+                        .sample_f(-w, uc, u, 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;
+
+                    if !exit_interface.flags().is_specular() {
+                        let f_exit = exit_interface.f(-w, wi, mode);
+                        if !f_exit.is_black() {
+                            let mut wt = 1.;
+                            if non_exit_interface.flags().is_specular() {
+                                let exit_pdf = exit_interface.pdf(-w, wi, f_args);
+                                wt = power_heuristic(1, bs.pdf, 1, exit_pdf);
+                            }
+                            f += beta * self.tr(self.thickness, bs.wi) * f_exit * wt;
+                        }
+                    }
+                }
+            }
+        }
+
+        f / self.n_samples as Float
+    }
+
+    fn sample_f(
+        &self,
+        _wo: Vector3f,
+        _uc: Float,
+        _u: Point2f,
+        _f_args: FArgs,
+    ) -> Option<BSDFSample> {
+        todo!()
+    }
+
+    fn pdf(&self, _wo: Vector3f, _wi: Vector3f, _f_args: FArgs) -> Float {
+        todo!()
+    }
+
+    fn regularize(&mut self) {
+        self.top.regularize();
+        self.bottom.regularize();
+    }
+
+    fn as_any(&self) -> &dyn Any {
+        todo!()
+    }
+}

src/utils/hash.rs

@@ -58,9 +58,6 @@ pub fn murmur_hash_64a(key: &[u8], seed: u64) -> u64 {
 }
 
 pub fn hash_buffer<T: ?Sized>(data: &T, seed: u64) -> u64 {
-    // In Rust, we need to treat the type as a byte slice.
-    // The safest way to do this generically without 'bytemuck' or 'unsafe' in signature
-    // is to require T to be castable, but usually in PBRT we just use unsafe for the hasher.
     let len = std::mem::size_of_val(data);
     let ptr = data as *const T as *const u8;
     let bytes = unsafe { std::slice::from_raw_parts(ptr, len) };

src/utils/sampling.rs

@@ -310,6 +310,10 @@ pub fn invert_spherical_rectangle_sample(
     }
 }
 
+pub fn sample_exponential(u: Float, a: Float) -> Float {
+    (1. - u).ln() / a
+}
+
 pub fn sample_spherical_triangle(
     v: &[Point3f; 3],
     p: Point3f,