From 7fcc42886a8d8b9bc30139a06adb872cd729e7f8 Mon Sep 17 00:00:00 2001 From: pingu Date: Wed, 10 Dec 2025 01:16:25 +0000 Subject: [PATCH] Just to be safe --- src/core/bxdf.rs | 532 +++++++++++++++++++++++++++++++++++++++++- src/utils/hash.rs | 3 - src/utils/sampling.rs | 4 + 3 files changed, 533 insertions(+), 6 deletions(-) diff --git a/src/core/bxdf.rs b/src/core/bxdf.rs index 86948ce..94ff82d 100644 --- a/src/core/bxdf.rs +++ b/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, - sample_uniform_hemisphere, uniform_hemisphere_pdf, + PiecewiseLinear2D, cosine_hemisphere_pdf, power_heuristic, sample_cosine_hemisphere, + 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 { + 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 +where + T: BxDFTrait, + B: BxDFTrait, +{ + top: T, + bottom: B, + thickness: Float, + g: Float, + albedo: SampledSpectrum, + max_depth: usize, + n_samples: usize, +} + +impl LayeredBxDF +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, TopOrBottom, TopOrBottom), + 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::().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 BxDFTrait for LayeredBxDF +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::().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 { + 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!() + } +} diff --git a/src/utils/hash.rs b/src/utils/hash.rs index e48893f..56f339f 100644 --- a/src/utils/hash.rs +++ b/src/utils/hash.rs @@ -58,9 +58,6 @@ pub fn murmur_hash_64a(key: &[u8], seed: u64) -> u64 { } pub fn hash_buffer(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) }; diff --git a/src/utils/sampling.rs b/src/utils/sampling.rs index 369a746..8ae1435 100644 --- a/src/utils/sampling.rs +++ b/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,