From 6e23698e2dd36dafcf0b18358ed613ca89d30fdb Mon Sep 17 00:00:00 2001 From: Wito Wiala Date: Wed, 2 Sep 2026 12:42:49 +0100 Subject: [PATCH] Added scene creator for SubSurface Material --- shared/src/core/bssrdf.rs | 31 ++-- shared/src/core/scattering.rs | 34 +++-- shared/src/materials/complex.rs | 28 ++-- shared/src/utils/math.rs | 30 ++++ src/core/medium.rs | 4 +- src/materials/complex.rs | 250 ++++++++++++++++++++++++++++++-- 6 files changed, 325 insertions(+), 52 deletions(-) diff --git a/shared/src/core/bssrdf.rs b/shared/src/core/bssrdf.rs index f1766c5..57a107d 100644 --- a/shared/src/core/bssrdf.rs +++ b/shared/src/core/bssrdf.rs @@ -3,11 +3,11 @@ use crate::core::bsdf::BSDF; use crate::core::geometry::{Frame, Normal3f, Point2f, Point3f, Point3fi, Vector3f}; use crate::core::interaction::{InteractionBase, ShadingGeom, SurfaceInteraction}; use crate::core::shape::Shape; -use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES}; +use crate::core::{LightIdx, MaterialIdx}; +use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum}; use crate::utils::math::{catmull_rom_weights, square}; use crate::utils::sampling::sample_catmull_rom_2d; -use crate::core::{LightIdx, MaterialIdx}; -use crate::{gvec_with_capacity, Float, GVec, PI, Ptr}; +use crate::{Float, GVec, PI, Ptr, gvec_with_capacity}; use enum_dispatch::enum_dispatch; use num_traits::Float as NumFloat; @@ -105,20 +105,19 @@ pub struct BSSRDFTable { impl BSSRDFTable { pub fn new(n_rho: usize, n_radius: usize) -> Self { - let rho_samples: GVec = gvec_with_capacity(n_rho); - let radius_samples: GVec = gvec_with_capacity(n_radius); - let profile: GVec = gvec_with_capacity(n_radius * n_rho); - let rho_eff: GVec = gvec_with_capacity(n_rho); - let profile_cdf: GVec = gvec_with_capacity(n_radius * n_rho); + let filled = |n: usize| { + let mut v: GVec = gvec_with_capacity(n); + v.resize(n, 0.); + v + }; Self { - n_rho: n_rho.try_into().unwrap(), - n_radius: n_radius.try_into().unwrap(), - rho_samples, - radius_samples, - profile, - rho_eff, - profile_cdf, - + n_rho: n_rho as u32, + n_radius: n_radius as u32, + rho_samples: filled(n_rho), + radius_samples: filled(n_radius), + profile: filled(n_rho * n_radius), + rho_eff: filled(n_rho), + profile_cdf: filled(n_rho * n_radius), } } diff --git a/shared/src/core/scattering.rs b/shared/src/core/scattering.rs index 5cd6c20..e1553f9 100644 --- a/shared/src/core/scattering.rs +++ b/shared/src/core/scattering.rs @@ -1,9 +1,9 @@ use crate::core::geometry::{ - abs_cos_theta, cos2_theta, cos_phi, sin_phi, tan2_theta, Normal3f, Point2f, Vector2f, Vector3f, - VectorLike, + Normal3f, Point2f, Vector2f, Vector3f, VectorLike, abs_cos_theta, cos_phi, cos2_theta, sin_phi, + tan2_theta, }; -use crate::core::pbrt::{Float, PI}; -use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES}; +use crate::core::pbrt::{Float, INV_4_PI, PI}; +use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum}; use crate::utils::math::{clamp, lerp, safe_sqrt, square}; use crate::utils::sampling::sample_uniform_disk_polar; use num_traits::Float as NumFloat; @@ -188,11 +188,9 @@ pub fn fresnel_moment1(eta: Float) -> Float { 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; + 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; + -4.61686 + 11.1136 * eta - 10.4646 * eta2 + 5.11455 * eta3 - 1.27198 * eta4 + 0.12746 * eta5 } } @@ -203,18 +201,28 @@ pub fn fresnel_moment2(eta: Float) -> Float { 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; + 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 + -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; + + 0.63942 * eta5 } } + +#[inline] +pub fn henyey_greenstein(cos_theta: Float, mut g: Float) -> Float { + // The Henyey-Greenstein phase function isn't suitable for |g| \approx + // 1 so we clamp it before it becomes numerically instable. (It's an + // analogous situation to BSDFs: if the BSDF is perfectly specular, one + // should use one based on a Dirac delta distribution rather than a + // very smooth microfacet distribution...) + g = g.clamp(-0.99, 0.99); + let denom = 1. + square(g) + 2. * g * cos_theta; + INV_4_PI * (1. - square(g)) / (denom * safe_sqrt(denom)) +} diff --git a/shared/src/materials/complex.rs b/shared/src/materials/complex.rs index 38a482c..70faffc 100644 --- a/shared/src/materials/complex.rs +++ b/shared/src/materials/complex.rs @@ -71,7 +71,7 @@ impl MaterialTrait for HairMaterial { let a = tex_eval.evaluate_float(&self.alpha, ctx); let e = tex_eval.evaluate_float(&self.eta, ctx); let sig_a = match self.hair_absorption { - // Absorption coefficient, not a reflectance, can be larger than 1 + // Absorption coefficient, not reflectance, can be larger than 1 HairAbsorption::SigmaA(sigma_a) => { SampledSpectrum::clamp_zero(&tex_eval.evaluate_spectrum(&sigma_a, ctx, lambda)) } @@ -147,11 +147,11 @@ impl MaterialTrait for MeasuredMaterial { fn get_bsdf( &self, _tex_eval: &T, - _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + ctx: &MaterialEvalContext, + lambda: &mut SampledWavelengths, ) -> BSDF { - // MeasuredBxDF::new(&self.brdf, lambda) - todo!() + let bxdf = BxDF::Measured(MeasuredBxDF::new(&self.brdf, lambda)); + BSDF::new(ctx.ns, ctx.dpdus, bxdf) } fn get_bssrdf( @@ -180,15 +180,25 @@ impl MaterialTrait for MeasuredMaterial { } } +#[repr(C)] +#[derive(Copy, Clone, Debug)] +pub enum SubsurfaceScattering { + Coefficients { + sigma_a: Ptr, + sigma_s: Ptr, + }, + Reflectance { + reflectance: Ptr, + mfp: Ptr, + }, +} + #[repr(C)] #[derive(Clone, Copy, Debug)] pub struct SubsurfaceMaterial { pub normal_map: Ptr, pub displacement: Ptr, - pub sigma_a: Ptr, - pub sigma_s: Ptr, - pub reflectance: Ptr, - pub mfp: Ptr, + pub scattering: SubsurfaceScattering, pub eta: Float, pub scale: Float, pub u_roughness: Ptr, diff --git a/shared/src/utils/math.rs b/shared/src/utils/math.rs index 928dbf4..6e9962a 100644 --- a/shared/src/utils/math.rs +++ b/shared/src/utils/math.rs @@ -349,6 +349,36 @@ pub fn wrap_equal_area_square(uv: &mut Point2f) -> Point2f { *uv } +pub fn integrate_catmull_rom(nodes: &[Float], f: &[Float], cdf: &mut [Float]) -> Float { + debug_assert_eq!(nodes.len(), f.len()); + let mut sum = 0.; + cdf[0] = 0.; + for i in 0..nodes.len() { + let x0 = nodes[i]; + let x1 = nodes[i + 1]; + let f0 = f[i]; + let f1 = f[i + 1]; + let width = x1 - x0; + + // Approximate derivatives using finite differences + let d0 = if i > 0 { + width * (f1 - f[i - 1]) / (x1 - nodes[i - 1]) + } else { + f1 - f0 + }; + let d1 = if i + 2 < nodes.len() { + width * (f[i + 2] - f0) / (nodes[i + 2] - x0) + } else { + f1 - f0 + }; + + // Keep a running sum and build a cumulative distribution function + sum += width * ((f0 + f1) / 2. + (d0 - d1) / 12.); + cdf[i + 1] = sum; + } + sum +} + pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(u32, [Float; 4])> { if nodes.len() < 4 { return None; diff --git a/src/core/medium.rs b/src/core/medium.rs index 036a96d..3508221 100644 --- a/src/core/medium.rs +++ b/src/core/medium.rs @@ -1,7 +1,7 @@ use crate::core::spectrum::spectrum_to_photometric; use crate::spectra::SRGB; use crate::{Arena, FileLoc, ParameterDictionary}; -use anyhow::{bail, Result}; +use anyhow::{Result, bail}; use shared::core::geometry::{Bounds3f, Point3f, Point3i}; use shared::core::medium::{ GridMedium, HGPhaseFunction, HomogeneousMedium, MajorantGrid, Medium, RGBGridMedium, @@ -262,7 +262,7 @@ static SUBSURFACE_TABLE: &[MeasuredSS] = &[ }, ]; -fn get_medium_scattering_properties(name: &str) -> Option<(Spectrum, Spectrum)> { +pub fn get_medium_scattering_properties(name: &str) -> Option<(Spectrum, Spectrum)> { SUBSURFACE_TABLE.iter().find(|m| m.name == name).map(|m| { let sigma_a = Spectrum::RGBUnbounded(RGBUnboundedSpectrum::new(&SRGB, m.sigma_a.into())); let sigma_s = diff --git a/src/materials/complex.rs b/src/materials/complex.rs index 609ed03..65c0519 100644 --- a/src/materials/complex.rs +++ b/src/materials/complex.rs @@ -1,16 +1,29 @@ use crate::core::image::HostImage; use crate::core::material::CreateMaterial; +use crate::core::medium::get_medium_scattering_properties; use crate::core::texture::SpectrumTexture; use crate::spectra::get_colorspace_device; -use crate::utils::TextureParameterDictionary; +use crate::utils::{TextureParameterDictionary, resolve_filename}; use crate::{Arena, ArenaUpload, FileLoc}; -use anyhow::Result; +use anyhow::{Result, anyhow, bail}; +use rayon::iter::{IndexedParallelIterator, IntoParallelRefMutIterator, ParallelIterator}; +use rayon::prelude::ParallelSliceMut; use shared::bxdfs::HairBxDF; +use shared::core::bssrdf::BSSRDFTable; +use shared::core::color::RGB; use shared::core::material::Material; +use shared::core::scattering::{ + fr_dielectric, fresnel_moment1, fresnel_moment2, henyey_greenstein, +}; use shared::core::spectrum::Spectrum; use shared::core::texture::SpectrumType; use shared::materials::complex::*; +use shared::spectra::{ConstantSpectrum, RGBUnboundedSpectrum}; use shared::textures::SpectrumConstantTexture; +use shared::utils::math::{fast_exp, integrate_catmull_rom, safe_sqrt, square}; +use shared::utils::sampling::sample_exponential; +use shared::{Float, INV_4_PI, PI}; + use std::collections::HashMap; use std::sync::Arc; @@ -65,26 +78,239 @@ impl CreateMaterial for HairMaterial { } } +fn cube(x: Float) -> Float { + x * x * x +} + +fn beam_diffusion_ms(sigma_s: Float, sigma_a: Float, g: Float, eta: Float, r: Float) -> Float { + const N_SAMPLES: usize = 100; + let mut e_d = 0.; + // Precompute information for dipole integrand + // Compute reduced scattering coefficients $\sigmaps, \sigmapt$ and albedo $\rhop$ + let sigmap_s = sigma_s * (1. - g); + let sigmap_t = sigma_a + sigmap_s; + let rhop = sigmap_s / sigmap_t; + + // Compute non-classical diffusion coefficient $D_\roman{G}$ using Equation + // $(\ref{eq:diffusion-coefficient-grosjean})$ + let d_g = (2. * sigma_a + sigmap_s) / (3. * sigmap_t * sigmap_t); + + // Compute effective transport coefficient $\sigmatr$ based on $D_\roman{G}$ + let sigma_tr = safe_sqrt(sigma_a / d_g); + + // Determine linear extrapolation distance $\depthextrapolation$ using Equation + // $(\ref{eq:dipole-boundary-condition})$ + let fm1 = fresnel_moment1(eta); + let fm2 = fresnel_moment2(eta); + let ze = -2. * d_g * (1. + 3. * fm2) / (1. - 2. * fm1); + + // Determine exitance scale factors using Equations $(\ref{eq:kp-exitance-phi})$ and + // $(\ref{eq:kp-exitance-e})$ + let c_phi = 0.25 * (1. - 2. * fm1); + let c_e = 0.5 * (1. - 3. * fm2); + + for i in 0..N_SAMPLES { + // Sample real point source depth $\depthreal$ + let zr = sample_exponential((i as Float + 0.5) / N_SAMPLES as Float, sigmap_t); + + // Evaluate dipole integrand $E_{\roman{d}}$ at $\depthreal$ and add to _Ed_ + let zv = -zr + 2. * ze; + let dr = (square(r) + square(zr)).sqrt(); + let dv = (square(r) + square(zv)).sqrt(); + // Compute dipole fluence rate $\dipole(r)$ using Equation + // $(\ref{eq:diffusion-dipole})$ + let phi_d = + INV_4_PI / d_g * (fast_exp(-sigma_tr * dr) / dr - fast_exp(-sigma_tr * dv) / dv); + + // Compute dipole vector irradiance $-\N{}\cdot\dipoleE(r)$ using Equation + // $(\ref{eq:diffusion-dipole-vector-irradiance-normal})$ + let e_dn = INV_4_PI + * (zr * (1. + sigma_tr * dr) * fast_exp(-sigma_tr * dr) / cube(dr) + - zv * (1. + sigma_tr * dv) * fast_exp(-sigma_tr * dv) / cube(dv)); + + // Add contribution from dipole for depth $\depthreal$ to _Ed_ + let e = phi_d * c_phi + e_dn * c_e; + let kappa = 1. - fast_exp(-2. * sigmap_t * (dr + zr)); + e_d += kappa * rhop * rhop * e; + } + return e_d / N_SAMPLES as Float; +} + +fn beam_diffusion_ss(sigma_s: Float, sigma_a: Float, g: Float, eta: Float, r: Float) -> Float { + // Compute material parameters and minimum $t$ below the critical angle + let sigma_t = sigma_a + sigma_s; + let rho = sigma_s / sigma_t; + let t_crit = r * safe_sqrt(square(eta) - 1.); + + let mut ess = 0.0; + const N_SAMPLES: usize = 100; + for i in 0..N_SAMPLES { + // Evaluate single-scattering integrand and add to _Ess_ + let ti = t_crit + sample_exponential((i as Float + 0.5) / N_SAMPLES as Float, sigma_t); + // Determine length $d$ of connecting segment and $\cos\theta_\roman{o}$ + let d = (square(r) + square(ti)).sqrt(); + let cos_theta_o = ti / d; + + // Add contribution of single scattering at depth $t$ + ess += rho * fast_exp(-sigma_t * (d + t_crit)) / square(d) + * henyey_greenstein(cos_theta_o, g) + * (1. - fr_dielectric(-cos_theta_o, eta)) + * cos_theta_o.abs(); + } + return ess / N_SAMPLES as Float; +} + +fn compute_beam_diffusion_bssrdf(g: Float, eta: Float, t: &mut BSSRDFTable) { + let n_rho = t.rho_samples.len(); + let n_radius = t.radius_samples.len(); + + t.radius_samples[0] = 0.; + t.radius_samples[1] = 2.5e-3; + for i in 2..n_radius { + t.radius_samples[i] = t.radius_samples[i - 1] * 1.2; + } + + for i in 0..n_rho { + t.rho_samples[i] = + (1. - fast_exp(-8. * i as Float / (n_rho - 1) as Float)) / (1. - fast_exp(-8.)); + } + + let rho_samples = &t.rho_samples; + let radius_samples = &t.radius_samples; + t.profile + .par_chunks_mut(n_radius) + .zip(t.profile_cdf.par_chunks_mut(n_radius)) + .zip(t.rho_eff.par_iter_mut()) + .enumerate() + .for_each(|(i, ((profile, cdf), rho_eff))| { + // Compute the diffusion profile for the _i_th albedo sample + // Compute scattering profile for chosen albedo $\rho$ + let rho = rho_samples[i]; + for j in 0..n_radius { + let r = radius_samples[j]; + profile[j] = 2. + * PI + * r + * (beam_diffusion_ss(rho, 1. - rho, g, eta, r) + + beam_diffusion_ms(rho, 1. - rho, g, eta, r)); + } + *rho_eff = integrate_catmull_rom(radius_samples, profile, cdf); + }); +} + impl CreateMaterial for SubsurfaceMaterial { fn create( - _parameters: &TextureParameterDictionary, - _normal_map: Option>, + parameters: &TextureParameterDictionary, + normal_map: Option>, _named_materials: &HashMap, - _loc: &FileLoc, - _arena: &Arena, + loc: &FileLoc, + arena: &Arena, ) -> Result { - todo!() + let mut g = parameters.get_one_float("g", 0.0)?; + let name = parameters.get_one_string("name", "")?; + let scattering = if !name.is_empty() { + let (sig_a, sig_s) = get_medium_scattering_properties(&name) + .ok_or_else(|| anyhow!("{loc}: named medium {name} not found"))?; + if g != 0. { + log::warn!("{loc}: non-zero \"g\" ignored with named scattering coefficients"); + } + g = 0.; + let sigma_a = SpectrumTexture::Constant(SpectrumConstantTexture::new(sig_a)); + let sigma_s = SpectrumTexture::Constant(SpectrumConstantTexture::new(sig_s)); + + SubsurfaceScattering::Coefficients { + sigma_a: arena.upload(&sigma_a), + sigma_s: arena.upload(&sigma_s), + } + } else { + let sigma_a = + parameters.get_spectrum_texture_or_null("sigma_a", SpectrumType::Unbounded); + let sigma_s = + parameters.get_spectrum_texture_or_null("sigma_s", SpectrumType::Unbounded); + match (sigma_a, sigma_s) { + (Some(a), Some(b)) => SubsurfaceScattering::Coefficients { + sigma_a: arena.upload(a), + sigma_s: arena.upload(b), + }, + (Some(_), None) => bail!("{loc}: provided \"sigma_a\" without \"sigma_s\""), + (None, Some(_)) => bail!("{loc}: provided \"sigma_s\" without \"sigma_a\""), + (None, None) => match parameters + .get_spectrum_texture_or_null("reflectance", SpectrumType::Albedo) + { + Some(r) => { + let one = Spectrum::Constant(ConstantSpectrum::new(1.)); + let mfp = parameters + .get_spectrum_texture("mfp", Some(one), SpectrumType::Unbounded) + .expect("default supplied"); + SubsurfaceScattering::Reflectance { + reflectance: arena.upload(r), + mfp: arena.upload(mfp), + } + } + None => { + let stdcs = get_colorspace_device(); + let default_sigma_a = + RGBUnboundedSpectrum::new(&stdcs.srgb, RGB::new(0.0011, 0.0024, 0.014)); + let default_sigma_s = + RGBUnboundedSpectrum::new(&stdcs.srgb, RGB::new(2.55, 3.21, 3.77)); + let sigma_a = SpectrumTexture::Constant(SpectrumConstantTexture::new( + Spectrum::RGBUnbounded(default_sigma_a), + )); + let sigma_s = SpectrumTexture::Constant(SpectrumConstantTexture::new( + Spectrum::RGBUnbounded(default_sigma_s), + )); + + SubsurfaceScattering::Coefficients { + sigma_a: arena.upload(&sigma_a), + sigma_s: arena.upload(&sigma_s), + } + } + }, + } + }; + + let scale = parameters.get_one_float("scale", 1.)?; + let eta = parameters.get_one_float("eta", 1.33)?; + + let u_roughness = + parameters.get_float_texture_with_fallback("uroughness", "roughness", 0.)?; + let v_roughness = + parameters.get_float_texture_with_fallback("vroughness", "roughness", 0.)?; + let displacement = parameters.get_float_texture_or_null("displacement")?; + let remap_roughness = parameters.get_one_bool("remaproughness", true)?; + let mut table = BSSRDFTable::new(100, 64); + compute_beam_diffusion_bssrdf(g, eta, &mut table); + + let mut ss_material = SubsurfaceMaterial { + scattering, + displacement: arena.upload(displacement), + normal_map: arena.upload(normal_map), + scale, + u_roughness: arena.upload(u_roughness), + v_roughness: arena.upload(v_roughness), + eta, + remap_roughness, + table: arena.alloc(table), + }; + + Ok(Material::Subsurface(ss_material)) } } impl CreateMaterial for MeasuredMaterial { fn create( - _parameters: &TextureParameterDictionary, - _normal_map: Option>, + parameters: &TextureParameterDictionary, + normal_map: Option>, _named_materials: &HashMap, - _loc: &FileLoc, - _arena: &Arena, + loc: &FileLoc, + arena: &Arena, ) -> Result { - todo!() + let filename = resolve_filename(parameters.get_one_string("filename", "")?); + let displacement = parameters.get_float_texture_or_null("displacement")?; + let brdf = MeasuredBxDF::brdf_data_from_file(filename); + let mat = MeasuredMaterial { + displacement: arena.upload(displacement), + normal_map: arena.upload(normal_map) + } } }