Added scene creator for SubSurface Material

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::{Float, GVec, PI, Ptr, gvec_with_capacity};
-use crate::{gvec_with_capacity, Float, GVec, PI, Ptr};
 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<Float> = gvec_with_capacity(n_rho);
+        let filled = |n: usize| {
-        let radius_samples: GVec<Float> = gvec_with_capacity(n_radius);
+            let mut v: GVec<Float> = gvec_with_capacity(n);
-        let profile: GVec<Float> = gvec_with_capacity(n_radius * n_rho);
+            v.resize(n, 0.);
-        let rho_eff: GVec<Float> = gvec_with_capacity(n_rho);
+            v
-        let profile_cdf: GVec<Float> = gvec_with_capacity(n_radius * n_rho);
+        };
         Self {
-            n_rho: n_rho.try_into().unwrap(),
+            n_rho: n_rho as u32,
-            n_radius: n_radius.try_into().unwrap(),
+            n_radius: n_radius as u32,
-            rho_samples,
+            rho_samples: filled(n_rho),
-            radius_samples,
+            radius_samples: filled(n_radius),
-            profile,
+            profile: filled(n_rho * n_radius),
-            rho_eff,
+            rho_eff: filled(n_rho),
-            profile_cdf,
+            profile_cdf: filled(n_rho * n_radius),
-
         }
     }
 

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,
+    Normal3f, Point2f, Vector2f, Vector3f, VectorLike, abs_cos_theta, cos_phi, cos2_theta, sin_phi,
-    VectorLike,
+    tan2_theta,
 };
-use crate::core::pbrt::{Float, PI};
+use crate::core::pbrt::{Float, INV_4_PI, PI};
-use crate::spectra::{SampledSpectrum, N_SPECTRUM_SAMPLES};
+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.45966 - 1.73965 * eta + 3.37668 * eta2 - 3.904945 * eta3 + 2.49277 * eta4 - 0.68441 * eta5
-            - 0.68441 * eta5;
     } else {
-        return -4.61686 + 11.1136 * eta - 10.4646 * eta2 + 5.11455 * eta3 - 1.27198 * eta4
+        -4.61686 + 11.1136 * eta - 10.4646 * eta2 + 5.11455 * eta3 - 1.27198 * eta4 + 0.12746 * eta5
-            + 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.27614 - 0.87350 * eta + 1.12077 * eta2 - 0.65095 * eta3 + 0.07883 * eta4 + 0.04860 * eta5
-            + 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))
+}

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<T: TextureEvaluator>(
         &self,
         _tex_eval: &T,
-        _ctx: &MaterialEvalContext,
+        ctx: &MaterialEvalContext,
-        _lambda: &mut SampledWavelengths,
+        lambda: &mut SampledWavelengths,
     ) -> BSDF {
-        // MeasuredBxDF::new(&self.brdf, lambda)
+        let bxdf = BxDF::Measured(MeasuredBxDF::new(&self.brdf, lambda));
-        todo!()
+        BSDF::new(ctx.ns, ctx.dpdus, bxdf)
     }
 
     fn get_bssrdf<T>(
@@ -180,15 +180,25 @@ impl MaterialTrait for MeasuredMaterial {
     }
 }
 
+#[repr(C)]
+#[derive(Copy, Clone, Debug)]
+pub enum SubsurfaceScattering {
+    Coefficients {
+        sigma_a: Ptr<SpectrumTexture>,
+        sigma_s: Ptr<SpectrumTexture>,
+    },
+    Reflectance {
+        reflectance: Ptr<SpectrumTexture>,
+        mfp: Ptr<SpectrumTexture>,
+    },
+}
+
 #[repr(C)]
 #[derive(Clone, Copy, Debug)]
 pub struct SubsurfaceMaterial {
     pub normal_map: Ptr<Image>,
     pub displacement: Ptr<FloatTexture>,
-    pub sigma_a: Ptr<SpectrumTexture>,
+    pub scattering: SubsurfaceScattering,
-    pub sigma_s: Ptr<SpectrumMixTexture>,
-    pub reflectance: Ptr<SpectrumMixTexture>,
-    pub mfp: Ptr<SpectrumMixTexture>,
     pub eta: Float,
     pub scale: Float,
     pub u_roughness: Ptr<FloatTexture>,

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;

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 =

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,
+        parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
-        _loc: &FileLoc,
+        loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Material> {
-        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,
+        parameters: &TextureParameterDictionary,
-        _normal_map: Option<Arc<HostImage>>,
+        normal_map: Option<Arc<HostImage>>,
         _named_materials: &HashMap<String, Material>,
-        _loc: &FileLoc,
+        loc: &FileLoc,
-        _arena: &Arena,
+        arena: &Arena,
     ) -> Result<Material> {
-        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)
+        }
     }
 }