diff --git a/src/core/bssrdf.rs b/src/core/bssrdf.rs index 74d1344..b247f9a 100644 --- a/src/core/bssrdf.rs +++ b/src/core/bssrdf.rs @@ -1,18 +1,20 @@ use crate::core::bxdf::BSDF; -use crate::core::interaction::SurfaceInteraction; +use crate::core::interaction::{InteractionData, ShadingGeometry, SurfaceInteraction}; use crate::core::pbrt::{Float, PI}; use crate::geometry::{Frame, Normal3f, Point2f, Point3f, Point3fi, Vector3f}; +use crate::shapes::Shape; use crate::spectra::{N_SPECTRUM_SAMPLES, SampledSpectrum}; use crate::utils::math::{catmull_rom_weights, square}; use crate::utils::sampling::sample_catmull_rom_2d; +use enum_dispatch::enum_dispatch; use std::sync::Arc; #[derive(Debug)] pub struct BSSRDFSample<'a> { - sp: SampledSpectrum, - pdf: SampledSpectrum, - sw: BSDF<'a>, - wo: Vector3f, + pub sp: SampledSpectrum, + pub pdf: SampledSpectrum, + pub sw: BSDF<'a>, + pub wo: Vector3f, } #[derive(Clone, Debug)] @@ -58,6 +60,43 @@ impl From for SubsurfaceInteraction { } } +impl From<&SubsurfaceInteraction> for SurfaceInteraction { + fn from(ssi: &SubsurfaceInteraction) -> SurfaceInteraction { + SurfaceInteraction { + common: InteractionData { + pi: ssi.pi, + n: ssi.n, + wo: Vector3f::zero(), + time: 0., + medium_interface: None, + medium: None, + }, + uv: Point2f::zero(), + dpdu: ssi.dpdu, + dpdv: ssi.dpdv, + dndu: Normal3f::zero(), + dndv: Normal3f::zero(), + shading: ShadingGeometry { + n: ssi.ns, + dpdu: ssi.dpdus, + dpdv: ssi.dpdvs, + dndu: Normal3f::zero(), + dndv: Normal3f::zero(), + }, + face_index: 0, + area_light: None, + material: None, + dpdx: Vector3f::zero(), + dpdy: Vector3f::zero(), + dudx: 0., + dvdx: 0., + dudy: 0., + dvdy: 0., + shape: Shape::default().into(), + } + } +} + #[derive(Clone, Debug)] pub struct BSSRDFTable { rho_samples: Vec, @@ -92,15 +131,18 @@ impl BSSRDFTable { #[derive(Clone, Default, Debug)] pub struct BSSRDFProbeSegment { - p0: Point3f, - p1: Point3f, + pub p0: Point3f, + pub p1: Point3f, } +#[enum_dispatch] pub trait BSSRDFTrait: Send + Sync + std::fmt::Debug { fn sample_sp(&self, u1: Float, u2: Point2f) -> Option; - fn probe_intersection_to_sample(si: &SubsurfaceInteraction) -> BSSRDFSample<'_>; + fn probe_intersection_to_sample(&self, si: &SubsurfaceInteraction) -> BSSRDFSample<'_>; } +#[enum_dispatch(BSSRDFTrait)] +#[derive(Debug, Clone)] pub enum BSSRDF<'a> { Tabulated(TabulatedBSSRDF<'a>), } @@ -279,7 +321,7 @@ impl<'a> BSSRDFTrait for TabulatedBSSRDF<'a> { }) } - fn probe_intersection_to_sample(_si: &SubsurfaceInteraction) -> BSSRDFSample<'_> { + fn probe_intersection_to_sample(&self, _si: &SubsurfaceInteraction) -> BSSRDFSample<'_> { todo!() } } diff --git a/src/core/bxdf.rs b/src/core/bxdf.rs index 1e88beb..86948ce 100644 --- a/src/core/bxdf.rs +++ b/src/core/bxdf.rs @@ -492,21 +492,12 @@ pub enum BxDF { // DiffuseTransmission(DiffuseTransmissionBxDF), } -#[derive(Debug)] +#[derive(Debug, Default)] pub struct BSDF<'a> { bxdf: Option<&'a mut dyn BxDFTrait>, shading_frame: Frame, } -impl<'a> Default for BSDF<'a> { - fn default() -> Self { - Self { - bxdf: None, - shading_frame: Frame::default(), - } - } -} - impl<'a> BSDF<'a> { pub fn new(ns: Normal3f, dpdus: Vector3f, bxdf: Option<&'a mut dyn BxDFTrait>) -> Self { Self { @@ -1375,8 +1366,7 @@ impl BxDFTrait for HairBxDF { 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)).expect("Could not sample from AP"); + 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], diff --git a/src/core/interaction.rs b/src/core/interaction.rs index 45352af..ce3688e 100644 --- a/src/core/interaction.rs +++ b/src/core/interaction.rs @@ -1,4 +1,5 @@ use crate::camera::{Camera, CameraTrait}; +use crate::core::bssrdf::BSSRDF; use crate::core::bxdf::{BSDF, BxDFFlags, DiffuseBxDF}; use crate::core::material::{ Material, MaterialEvalContext, MaterialTrait, NormalBumpEvalContext, bump_map, normal_map, @@ -343,6 +344,30 @@ impl SurfaceInteraction { Some(bsdf) } + pub fn get_bssrdf( + &self, + _ray: &Ray, + lambda: &SampledWavelengths, + _camera: &Camera, + _scratch: &Bump, + ) -> Option> { + let material = { + let root_mat = self.material.as_deref()?; + let mut active_mat: &Material = root_mat; + let tex_eval = UniversalTextureEvaluator; + while let Material::Mix(mix) = active_mat { + // We need a context to evaluate the 'amount' texture + let ctx = MaterialEvalContext::from(self); + active_mat = mix.choose_material(&tex_eval, &ctx); + } + active_mat.clone() + }; + + let ctx = MaterialEvalContext::from(self); + let tex_eval = UniversalTextureEvaluator; + material.get_bssrdf(&tex_eval, &ctx, lambda) + } + fn compute_bump_geometry( &mut self, tex_eval: &UniversalTextureEvaluator, @@ -628,6 +653,29 @@ pub struct MediumInteraction { pub phase: PhaseFunction, } +impl MediumInteraction { + pub fn new( + p: Point3f, + wo: Vector3f, + time: Float, + medium: Arc, + phase: PhaseFunction, + ) -> Self { + Self { + common: InteractionData { + pi: Point3fi::new_from_point(p), + n: Normal3f::default(), + time, + wo: wo.normalize(), + medium_interface: None, + medium: Some(medium.clone()), + }, + medium, + phase, + } + } +} + impl InteractionTrait for MediumInteraction { fn is_medium_interaction(&self) -> bool { true diff --git a/src/core/material.rs b/src/core/material.rs index 65c68ca..9c61e2d 100644 --- a/src/core/material.rs +++ b/src/core/material.rs @@ -160,7 +160,7 @@ pub trait MaterialTrait: Send + Sync + std::fmt::Debug { &self, tex_eval: &T, ctx: &MaterialEvalContext, - lambda: &mut SampledWavelengths, + lambda: &SampledWavelengths, ) -> Option>; fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool; @@ -201,7 +201,7 @@ impl MaterialTrait for CoatedDiffuseMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -238,7 +238,7 @@ impl MaterialTrait for CoatedConductorMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -272,7 +272,7 @@ impl MaterialTrait for ConductorMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -334,7 +334,7 @@ impl MaterialTrait for DielectricMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { None } @@ -379,7 +379,7 @@ impl MaterialTrait for DiffuseMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -414,7 +414,7 @@ impl MaterialTrait for DiffuseTransmissionMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -448,7 +448,7 @@ impl MaterialTrait for HairMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -483,7 +483,7 @@ impl MaterialTrait for MeasuredMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -517,7 +517,7 @@ impl MaterialTrait for SubsurfaceMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -551,7 +551,7 @@ impl MaterialTrait for ThinDielectricMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { todo!() } @@ -615,7 +615,7 @@ impl MaterialTrait for MixMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &mut SampledWavelengths, + _lambda: &SampledWavelengths, ) -> Option> { None } diff --git a/src/core/medium.rs b/src/core/medium.rs index bbb695e..9e78d75 100644 --- a/src/core/medium.rs +++ b/src/core/medium.rs @@ -314,10 +314,10 @@ impl Iterator for DDAMajorantIterator { } pub struct MediumProperties { - sigma_a: SampledSpectrum, - sigma_s: SampledSpectrum, - phase: PhaseFunction, - le: SampledSpectrum, + pub sigma_a: SampledSpectrum, + pub sigma_s: SampledSpectrum, + pub phase: PhaseFunction, + pub le: SampledSpectrum, } impl MediumProperties { @@ -343,12 +343,12 @@ pub trait MediumTrait: Send + Sync + std::fmt::Debug { mut ray: Ray, mut t_max: Float, u: Float, - _rng: &mut Rng, + rng: &mut Rng, lambda: &SampledWavelengths, mut callback: F, ) -> SampledSpectrum where - F: FnMut(Point3f, MediumProperties, SampledSpectrum, SampledSpectrum) -> bool, + F: FnMut(Point3f, MediumProperties, SampledSpectrum, SampledSpectrum, &mut Rng) -> bool, { let len = ray.d.norm(); t_max *= len; @@ -377,7 +377,7 @@ pub trait MediumTrait: Send + Sync + std::fmt::Debug { let p = ray.at(t); let mp = self.sample_point(p, lambda); - if !callback(p, mp, seg.sigma_maj, t_maj) { + if !callback(p, mp, seg.sigma_maj, t_maj, rng) { return SampledSpectrum::new(1.); } @@ -398,6 +398,7 @@ pub trait MediumTrait: Send + Sync + std::fmt::Debug { } #[derive(Debug, Clone)] +#[enum_dispatch(MediumTrait)] pub enum Medium { Homogeneous(HomogeneousMedium), Grid(GridMedium), @@ -754,8 +755,42 @@ impl MediumTrait for RGBGridMedium { #[derive(Debug, Clone)] pub struct CloudMedium; +impl MediumTrait for CloudMedium { + fn is_emissive(&self) -> bool { + todo!() + } + fn sample_point(&self, _p: Point3f, _lambda: &SampledWavelengths) -> MediumProperties { + todo!() + } + fn sample_ray( + &self, + _ray: &Ray, + _t_max: Float, + _lambda: &SampledWavelengths, + _buf: &Bump, + ) -> RayMajorantIterator { + todo!() + } +} #[derive(Debug, Clone)] pub struct NanoVDBMedium; +impl MediumTrait for NanoVDBMedium { + fn is_emissive(&self) -> bool { + todo!() + } + fn sample_point(&self, _p: Point3f, _lambda: &SampledWavelengths) -> MediumProperties { + todo!() + } + fn sample_ray( + &self, + _ray: &Ray, + _t_max: Float, + _lambda: &SampledWavelengths, + _buf: &Bump, + ) -> RayMajorantIterator { + todo!() + } +} #[derive(Debug, Default, Clone)] pub struct MediumInterface { diff --git a/src/core/texture.rs b/src/core/texture.rs index e19a763..59e24a1 100644 --- a/src/core/texture.rs +++ b/src/core/texture.rs @@ -725,7 +725,7 @@ impl ImageTextureBase { } let path = Path::new(&filename); - let mipmap_raw = MIPMap::create_from_file(&path, filter_options, wrap_mode, encoding) + let mipmap_raw = MIPMap::create_from_file(path, filter_options, wrap_mode, encoding) .expect("Failed to create MIPMap from file"); let mipmap_arc = Arc::new(mipmap_raw); diff --git a/src/integrators/mod.rs b/src/integrators/mod.rs index 6125927..a5b64d2 100644 --- a/src/integrators/mod.rs +++ b/src/integrators/mod.rs @@ -3,9 +3,13 @@ mod pipeline; use pipeline::*; use crate::camera::{Camera, CameraTrait}; -use crate::core::bxdf::{BSDF, BxDFFlags, BxDFTrait, FArgs, TransportMode}; +use crate::core::bssrdf::{BSSRDFTrait, SubsurfaceInteraction}; +use crate::core::bxdf::{BSDF, BxDFFlags, BxDFReflTransFlags, BxDFTrait, FArgs, TransportMode}; use crate::core::film::{FilmTrait, VisibleSurface}; -use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction}; +use crate::core::interaction::{ + self, Interaction, InteractionTrait, MediumInteraction, SimpleInteraction, SurfaceInteraction, +}; +use crate::core::medium::{MediumTrait, PhaseFunctionTrait}; use crate::core::options::get_options; use crate::core::pbrt::{Float, SHADOW_EPSILON}; use crate::core::primitive::{Primitive, PrimitiveTrait}; @@ -15,10 +19,12 @@ use crate::lights::sampler::{LightSampler, UniformLightSampler}; use crate::lights::{Light, LightSampleContext, LightTrait, sampler::LightSamplerTrait}; use crate::shapes::ShapeIntersection; use crate::spectra::{SampledSpectrum, SampledWavelengths}; -use crate::utils::math::square; +use crate::utils::hash::{hash_buffer, mix_bits}; +use crate::utils::math::{float_to_bits, sample_discrete, square}; +use crate::utils::rng::Rng; use crate::utils::sampling::{ - power_heuristic, sample_uniform_hemisphere, sample_uniform_sphere, uniform_hemisphere_pdf, - uniform_sphere_pdf, + WeightedReservoirSampler, power_heuristic, sample_uniform_hemisphere, sample_uniform_sphere, + uniform_hemisphere_pdf, uniform_sphere_pdf, }; use bumpalo::Bump; use rayon::prelude::*; @@ -244,10 +250,10 @@ impl RayIntegratorTrait for SimplePathIntegrator { if !self.sample_lights || specular_bounce { l += beta * isect.le(-ray.d, lambda); } - depth += 1; if depth == self.max_depth { break; } + depth += 1; let Some(bsdf) = isect.get_bsdf(&ray, lambda, self.camera.as_ref(), scratch, sampler) else { @@ -459,16 +465,16 @@ impl RayIntegratorTrait for PathIntegrator { bsdf.regularize(); } - depth += 1; - if depth == self.max_depth { break; } - if BxDFFlags::is_non_specular(&bsdf.flags()) { - if let Some(ld) = self.sample_ld(isect, &bsdf, lambda, sampler) { - l += beta * ld; - } + depth += 1; + + if bsdf.flags().is_non_specular() + && let Some(ld) = self.sample_ld(isect, &bsdf, lambda, sampler) + { + l += beta * ld; } let wo = -ray.d; @@ -507,6 +513,669 @@ impl RayIntegratorTrait for PathIntegrator { } } +#[derive(Clone, Debug)] +pub struct SimpleVolPathIntegrator { + base: IntegratorBase, + max_depth: usize, + camera: Arc, + sampler_prototype: Sampler, +} + +impl RayIntegratorTrait for SimpleVolPathIntegrator { + fn evaluate_pixel_sample( + &self, + p_pixel: Point2i, + sample_ind: usize, + sampler: &mut Sampler, + scratch: &Bump, + ) { + pipeline::evaluate_pixel_sample( + self, + self.camera.as_ref(), + sampler, + p_pixel, + sample_ind, + scratch, + ); + } + + fn li( + &self, + mut ray: Ray, + lambda: &SampledWavelengths, + sampler: &mut Sampler, + scratch: &Bump, + _visible_surface: bool, + ) -> (SampledSpectrum, Option) { + let mut l = SampledSpectrum::new(0.); + let mut beta = 1.; + let mut depth = 0; + + let lambda = lambda.terminate_secondary(); + + loop { + let si = self.base.intersect(&ray, None); + let mut scattered = false; + let mut terminated = false; + let current_medium = ray.medium.take(); + if let Some(medium) = ¤t_medium { + let hash0 = hash_buffer(&sampler.get1d(), 0); + let hash1 = hash_buffer(&sampler.get1d(), 0); + let mut rng = Rng::new_with_offset(hash0, hash1); + let t_max = if let Some(hit) = &si { + hit.t_hit + } else { + Float::INFINITY + }; + let u = sampler.get1d(); + let mut u_mode = sampler.get1d(); + medium.as_ref().sample_t_maj( + ray.clone(), + t_max, + u, + &mut rng, + &lambda, + |p, mp, sigma_maj, _t_maj, rng| -> bool { + let sig_maj = sigma_maj[0]; + // Compute medium event probabilities for interaction + let p_absorb: Float = mp.sigma_a[0] / sig_maj; + let p_scatter: Float = mp.sigma_s[0] / sig_maj; + let p_null = (1. - p_absorb - p_scatter).max(0_f32); + let mode = + sample_discrete(&[p_absorb, p_scatter, p_null], u_mode, None, None); + if mode == 0 { + // Handle absorption event for medium sample + l += beta * mp.le; + terminated = true; + false + } else if mode == 1 { + // Handle regular scattering event for medium sample + if depth >= self.max_depth { + return false; + } + depth += 1; + + let u_phase = Point2f::new(rng.uniform(), rng.uniform()); + + let Some(ps) = mp.phase.sample_p(-ray.d, u_phase) else { + terminated = true; + return false; + }; + + beta *= ps.p / ps.pdf; + ray.o = p; + ray.d = ps.wi; + scattered = true; + false + } else { + // Handle null-scattering event for medium sample + u_mode = rng.uniform(); + true + } + }, + ); + } + // Handle terminated and unscattered rays after medium sampling + if terminated { + return (l, None); + } + + if scattered { + return (l, None); + } + + let Some(mut hit) = si else { + for light in &self.base.infinite_lights { + l += beta * light.le(&ray, &lambda); + } + return (l, None); + }; + + // Handle surface intersection along ray path + if let Some(bsdf) = &hit + .intr + .get_bsdf(&ray, &lambda, &self.camera, scratch, sampler) + { + let uc = sampler.get1d(); + let u = sampler.get2d(); + if bsdf.sample_f(-ray.d, uc, u, FArgs::default()).is_some() { + panic!("SimpleVolPathIntegrator does not support surface scattering"); + } + } else { + hit.intr.skip_intersection(&mut ray, hit.t_hit); + } + } + } +} + +#[derive(Clone, Debug)] +pub struct VolPathIntegrator { + base: IntegratorBase, + sampler_prototype: Sampler, + camera: Arc, + light_sampler: LightSampler, + regularize: bool, + max_depth: usize, +} + +impl VolPathIntegrator { + fn sample_ld( + &self, + intr: &Interaction, + bsdf: Option<&BSDF>, + lambda: &SampledWavelengths, + sampler: &mut Sampler, + beta: SampledSpectrum, + r_p: SampledSpectrum, + ) -> SampledSpectrum { + let mut ctx = LightSampleContext::from(intr); + if let Some(bsdf) = bsdf { + let flags = bsdf.flags(); + let is_refl = flags.is_reflective(); + let is_trans = flags.is_transmissive(); + if is_refl ^ is_trans { + let vector = if is_refl { intr.wo() } else { -intr.wo() }; + ctx.pi = Point3fi::new_from_point(intr.offset_ray_vector(vector)); + } + } + let u = sampler.get1d(); + let Some(sampled_light) = self.light_sampler.sample_with_context(&ctx, u) else { + return SampledSpectrum::new(0.); + }; + let u_light = sampler.get2d(); + let light = sampled_light.light; + debug_assert!(sampled_light.p != 0.); + let Some(ls) = light + .sample_li(&ctx, u_light, lambda, true) + .filter(|ls| !ls.l.is_black() && ls.pdf > 0.0) + else { + return SampledSpectrum::new(0.0); + }; + let p_l = sampled_light.p * ls.pdf; + let wo = intr.wo(); + let wi = ls.wi; + let (f_hat, scatter_pdf) = match (bsdf, intr) { + (Some(bsdf), Interaction::Surface(si)) => { + let f = bsdf + .f(wo, wi, TransportMode::Radiance) + .map(|f_val| f_val * wi.abs_dot(si.shading.n.into())) + .unwrap_or(SampledSpectrum::new(0.0)); + + let pdf = bsdf.pdf(wo, wi, FArgs::default()); + (f, pdf) + } + + (None, Interaction::Medium(mi)) => { + let phase = &mi.phase; + let p = phase.p(wo, wi); + let pdf = phase.pdf(wo, wi); + + (SampledSpectrum::new(p), pdf) + } + + (Some(_), _) => panic!("BSDF provided for a non-Surface interaction"), + (None, _) => panic!("No BSDF provided for a non-Medium interaction"), + }; + + if f_hat.is_black() { + return SampledSpectrum::new(0.); + } + + let mut light_ray = intr.spawn_ray_to_interaction(ls.p_light.as_ref()); + let mut t_ray = SampledSpectrum::new(1.); + let mut r_l = SampledSpectrum::new(1.); + let mut r_u = SampledSpectrum::new(1.); + let hash0 = hash_buffer(&[light_ray.o.x(), light_ray.o.y(), light_ray.o.z()], 0); + let hash1 = hash_buffer(&[light_ray.d.x(), light_ray.d.y(), light_ray.d.z()], 0); + let mut rng = Rng::new_with_offset(hash0, hash1); + + while light_ray.d != Vector3f::zero() { + // Trace ray through media to estimate transmittance + let si = &self.base.intersect(&light_ray, Some(1. - SHADOW_EPSILON)); + if si.as_ref().filter(|s| s.intr.material.is_some()).is_some() { + return SampledSpectrum::new(0.); + } + + let current_medium = light_ray.medium.take(); + if let Some(medium) = ¤t_medium { + let t_max = if let Some(s) = si { + s.t_hit() + } else { + 1. - SHADOW_EPSILON + }; + + let u: Float = rng.uniform(); + let t_maj = medium.as_ref().sample_t_maj( + light_ray.clone(), + t_max, + u, + &mut rng, + lambda, + |_p, mp, sigma_maj, t_maj, rng| { + let sigma_n = + SampledSpectrum::clamp_zero(&(sigma_maj - mp.sigma_a - mp.sigma_s)); + let pdf = t_maj[0] * sigma_maj[0]; + t_ray *= t_maj * sigma_n / pdf; + r_l *= t_maj * sigma_maj / pdf; + r_u *= t_maj * sigma_n / pdf; + + let tr = t_ray / (r_l + r_u).average(); + if tr.max_component_value() < 0.05 { + let q = 0.75; + if rng.uniform::() < 1. { + t_ray = SampledSpectrum::new(0.); + } else { + t_ray /= 1. - q; + } + } + + if t_ray.is_black() { + return false; + } + true + }, + ); + t_ray *= t_maj / t_maj[0]; + r_l *= t_maj / t_maj[0]; + r_u *= t_maj / t_maj[0]; + } + + if t_ray.is_black() { + return SampledSpectrum::new(0.); + } + + if let Some(s) = si { + light_ray = s.intr.spawn_ray_to_interaction(ls.p_light.as_ref()); + } else { + break; + } + } + + r_l *= r_p * p_l; + r_u *= r_p * scatter_pdf; + + if light.light_type().is_delta_light() { + beta * f_hat * t_ray * ls.l / r_l.average() + } else { + beta * f_hat * t_ray * ls.l / (r_l + r_u).average() + } + } +} + +impl RayIntegratorTrait for VolPathIntegrator { + fn evaluate_pixel_sample( + &self, + p_pixel: Point2i, + sample_ind: usize, + sampler: &mut Sampler, + scratch: &Bump, + ) { + pipeline::evaluate_pixel_sample( + self, + self.camera.as_ref(), + sampler, + p_pixel, + sample_ind, + scratch, + ); + } + + fn li( + &self, + mut ray: Ray, + lambda: &SampledWavelengths, + sampler: &mut Sampler, + scratch: &Bump, + visible_surface: bool, + ) -> (SampledSpectrum, Option) { + let mut l = SampledSpectrum::new(0.); + let mut beta = SampledSpectrum::new(1.); + let mut r_u = SampledSpectrum::new(1.); + let mut r_l = SampledSpectrum::new(1.); + let mut specular_bounce = false; + let mut any_non_specular_bounces = false; + let mut depth = 0; + let mut eta_scale = 1.; + let mut prev_int_ctx = LightSampleContext::default(); + let mut visible: Option = None; + + loop { + let si = self.base.intersect(&ray, None); + let current_medium = ray.medium.take(); + if let Some(medium) = ¤t_medium { + let mut scattered = false; + let mut terminated = false; + let hash0 = hash_buffer(&sampler.get1d(), 0); + let hash1 = hash_buffer(&sampler.get1d(), 0); + let mut rng = Rng::new_with_offset(hash0, hash1); + + let t_max = if let Some(hit) = &si { + hit.t_hit + } else { + Float::INFINITY + }; + let t_maj = medium.as_ref().sample_t_maj( + ray.clone(), + t_max, + sampler.get1d(), + &mut rng, + lambda, + |p, mp, sigma_maj, t_maj, rng| -> bool { + // Handle medium scattering event for ray + if beta.is_black() { + terminated = true; + return false; + } + + if depth < self.max_depth && !mp.le.is_black() { + let pdf = sigma_maj[0] * t_maj[0]; + let betap = beta * t_maj / pdf; + let r_e = r_u * sigma_maj * t_maj / pdf; + if !r_e.is_black() { + l += betap * mp.sigma_a * mp.le / r_e.average(); + } + } + let p_absorb = mp.sigma_a[0] / sigma_maj[0]; + let p_scatter = mp.sigma_s[0] / sigma_maj[0]; + let p_null = (1. - p_absorb - p_scatter).max(0.); + assert!(1. - p_absorb - p_scatter > -1e-6); + // Sample medium scattering event type and update path + let um = rng.uniform(); + let mode = sample_discrete(&[p_absorb, p_scatter, p_null], um, None, None); + if mode == 0 { + terminated = true; + false + } else if mode == 1 { + // Handle scattering along ray path + // Stop path sampling if maximum depth has been reached + if depth >= self.max_depth { + terminated = true; + return false; + } + depth += 1; + + let pdf = t_maj[0] * mp.sigma_s[0]; + beta *= t_maj * mp.sigma_s / pdf; + r_u *= t_maj * mp.sigma_s / pdf; + + if !beta.is_black() && !r_u.is_black() { + let ray_medium = + ray.medium.as_ref().expect("Void scattering").clone(); + let intr = MediumInteraction::new( + p, -ray.d, ray.time, ray_medium, mp.phase, + ); + l += self.sample_ld( + &Interaction::Medium(intr.clone()), + None, + lambda, + sampler, + beta, + r_u, + ); + let u = sampler.get2d(); + let Some(ps) = + intr.phase.sample_p(-ray.d, u).filter(|ps| ps.pdf > 0.0) + else { + terminated = true; + // Stop traversal (Path died) + return false; + }; + + beta *= ps.p / ps.pdf; + r_l = r_u / ps.pdf; + prev_int_ctx = LightSampleContext::from(&intr); + scattered = true; + ray.o = p; + ray.d = ps.wi; + specular_bounce = false; + any_non_specular_bounces = false; + } + false + } else { + let sigma_n = + SampledSpectrum::clamp_zero(&(sigma_maj - mp.sigma_a - mp.sigma_s)); + let pdf = t_maj[0] * sigma_n[0]; + beta *= t_maj * sigma_n / pdf; + if pdf == 0. { + beta = SampledSpectrum::new(0.); + } + r_u *= t_maj * sigma_n / pdf; + r_l *= t_maj * sigma_maj / pdf; + !beta.is_black() && !r_u.is_black() + } + }, + ); + + if terminated || beta.is_black() || r_u.is_black() { + return (l, None); + } + if scattered { + continue; + } + beta *= t_maj / t_maj[0]; + r_u *= t_maj / t_maj[0]; + r_l *= t_maj / t_maj[0]; + } + + // Handle surviving unscattered rays + // Add emitted light at volume path vertex or from the environment + let Some(mut hit) = si else { + for light in &self.base.infinite_lights { + let le = light.le(&ray, lambda); + if !le.is_black() { + if depth == 0 || specular_bounce { + l += beta * light.le(&ray, lambda); + } + } else { + // Add infinite light contribution using both PDFs with MIS + let p_l = self.light_sampler.pmf_with_context(&prev_int_ctx, light) + * light.pdf_li(&prev_int_ctx, ray.d, true); + r_l *= p_l; + l += beta * le / (r_u + r_l).average(); + } + } + break; + }; + + let le = &hit.intr.le(-ray.d, lambda); + if !le.is_black() { + if depth == 0 || specular_bounce { + l += beta * *le / r_u.average(); + } + } else { + // Add surface light contribution using both PDFs with MIS + let area_light = hit.intr.area_light.as_ref().expect("No surface light"); + let p_l = self + .light_sampler + .pmf_with_context(&prev_int_ctx, area_light); + r_l *= p_l; + l += beta * *le / (r_u + r_l).average(); + } + + let isect = &mut hit.intr; + let Some(mut bsdf) = isect.get_bsdf(&ray, lambda, &self.camera, scratch, sampler) + else { + hit.intr.skip_intersection(&mut ray, hit.t_hit()); + continue; + }; + + if depth == 0 && visible_surface { + const N_RHO_SAMPLES: usize = 16; + const UC_RHO: [Float; N_RHO_SAMPLES] = [ + 0.75741637, + 0.37870818, + 0.7083487, + 0.18935409, + 0.9149363, + 0.35417435, + 0.5990858, + 0.09467703, + 0.8578725, + 0.45746812, + 0.686759, + 0.17708716, + 0.9674518, + 0.2995429, + 0.5083201, + 0.047338516, + ]; + let u_rho: [Point2f; N_RHO_SAMPLES] = [ + Point2f::new(0.855985, 0.570367), + Point2f::new(0.381823, 0.851844), + Point2f::new(0.285328, 0.764262), + Point2f::new(0.733380, 0.114073), + Point2f::new(0.542663, 0.344465), + Point2f::new(0.127274, 0.414848), + Point2f::new(0.964700, 0.947162), + Point2f::new(0.594089, 0.643463), + Point2f::new(0.095109, 0.170369), + Point2f::new(0.825444, 0.263359), + Point2f::new(0.429467, 0.454469), + Point2f::new(0.244460, 0.816459), + Point2f::new(0.756135, 0.731258), + Point2f::new(0.516165, 0.152852), + Point2f::new(0.180888, 0.214174), + Point2f::new(0.898579, 0.503897), + ]; + let albedo = bsdf.rho_wo(isect.wo(), &UC_RHO, &u_rho); + visible = Some(VisibleSurface::new(isect, &albedo, lambda)); + } + + if depth >= self.max_depth { + return (l, visible); + } + + depth += 1; + + if self.regularize && any_non_specular_bounces { + bsdf.regularize(); + } + + if BxDFFlags::is_non_specular(&bsdf.flags()) { + l += self.sample_ld( + &Interaction::Surface(isect.clone()), + Some(&bsdf), + lambda, + sampler, + beta, + r_u, + ); + debug_assert!(l.y(lambda).is_finite()); + } + let wo = isect.wo(); + let n = isect.shading.n; + prev_int_ctx = LightSampleContext::from(&*isect); + let u = sampler.get1d(); + let Some(bs) = bsdf.sample_f(wo, u, sampler.get2d(), FArgs::default()) else { + break; + }; + beta *= bs.f * bs.wi.abs_dot(n.into()) / bs.pdf; + if bs.pdf_is_proportional { + r_l = r_u / bsdf.pdf(wo, bs.wi, FArgs::default()); + } else { + r_l = r_u / bs.pdf; + } + debug_assert!(beta.y(lambda).is_finite()); + // Update volumetric integrator path state after surface scattering + specular_bounce = bs.is_specular(); + if bs.is_transmissive() { + eta_scale *= square(bs.eta); + } + ray = isect.spawn_ray_with_differentials(&ray, bs.wi, bs.flags, bs.eta); + + if let Some(bssrdf) = (*isect).get_bssrdf(&ray, lambda, self.camera.as_ref(), scratch) + && bs.is_transmissive() + { + let uc = sampler.get1d(); + let up = sampler.get2d(); + let Some(probe_seg) = bssrdf.sample_sp(uc, up) else { + break; + }; + let seed = mix_bits(float_to_bits(sampler.get1d()).into()); + let mut interaction_sampler = + WeightedReservoirSampler::::new(seed); + let base = + SimpleInteraction::new(Point3fi::new_from_point(probe_seg.p0), ray.time, None); + loop { + let r = base.spawn_ray_to_point(probe_seg.p1); + if r.d == Vector3f::zero() { + break; + } + let Some(si) = self.base.intersect(&r, Some(1.)) else { + break; + }; + + let materials_match = match (&si.intr.material, &isect.material) { + (Some(m1), Some(m2)) => Arc::ptr_eq(m1, m2), + _ => false, + }; + + if materials_match { + interaction_sampler.add(SubsurfaceInteraction::new(&si.intr), 1.); + } + } + + if !interaction_sampler.has_sample() { + break; + } + + if let Some(ssi) = interaction_sampler.get_sample() { + let bssrdf_sample = bssrdf.probe_intersection_to_sample(ssi); + if bssrdf_sample.sp.is_black() || bssrdf_sample.pdf.is_black() { + break; + } + let pdf = interaction_sampler.sample_probability() * bssrdf_sample.pdf[0]; + beta *= bssrdf_sample.sp / pdf; + r_u *= bssrdf_sample.pdf / bssrdf_sample.pdf[0]; + let mut pi = SurfaceInteraction::from(ssi); + pi.common.wo = bssrdf_sample.wo; + prev_int_ctx = LightSampleContext::from(&pi); + let mut sw = bssrdf_sample.sw; + if self.regularize { + sw.regularize(); + } + + l += self.sample_ld( + &Interaction::Surface(pi.clone()), + Some(&sw), + lambda, + sampler, + beta, + r_u, + ); + + let u = sampler.get1d(); + let Some(bs) = sw.sample_f(pi.wo(), u, sampler.get2d(), FArgs::default()) + else { + break; + }; + beta *= bs.f * bs.wi.abs_dot(pi.shading.n.into()) / bs.pdf; + r_l = r_u / bs.pdf; + debug_assert!(beta.y(lambda).is_finite()); + specular_bounce = bs.is_specular(); + ray = pi.spawn_ray(bs.wi); + } + } + + if beta.is_black() { + break; + } + + let rr_beta = beta * eta_scale / r_u.average(); + let u_rr = sampler.get1d(); + if rr_beta.max_component_value() < 1. && depth > 1 { + let q = (1. - rr_beta.max_component_value()).max(0.); + if u_rr < 1. { + break; + } + beta /= 1. - q; + } + } + + (l, visible) + } +} + pub enum Integrator { BPDT(BPDTIntegrator), MLT(MLTIntegrator), diff --git a/src/lights/diffuse.rs b/src/lights/diffuse.rs index a5b7183..7734a16 100644 --- a/src/lights/diffuse.rs +++ b/src/lights/diffuse.rs @@ -169,7 +169,7 @@ impl LightTrait for DiffuseAreaLight { } let wi = (p - ctx.p()).normalize(); - let le = self.l(p, n, uv, -wi, &lambda); + let le = self.l(p, n, uv, -wi, lambda); if le.is_black() { return None; diff --git a/src/lights/infinite.rs b/src/lights/infinite.rs index 0de25a7..c3299ec 100644 --- a/src/lights/infinite.rs +++ b/src/lights/infinite.rs @@ -78,7 +78,7 @@ impl LightTrait for InfiniteUniformLight { let intr = Interaction::Simple(intr_simple); Some(LightLiSample::new( - self.scale * self.lemit.sample(&lambda), + self.scale * self.lemit.sample(lambda), wi, pdf, intr, @@ -244,12 +244,7 @@ impl LightTrait for InfiniteImageLight { 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, - )) + Some(LightLiSample::new(self.image_le(uv, lambda), wi, pdf, intr)) } fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { @@ -541,7 +536,7 @@ impl LightTrait for InfinitePortalLight { return None; } let pdf = map_pdf / duv_dw; - let l = self.image_lookup(uv, &lambda); + 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 intr = Interaction::Simple(sintr); diff --git a/src/lights/mod.rs b/src/lights/mod.rs index e47fc2f..3e7a5c4 100644 --- a/src/lights/mod.rs +++ b/src/lights/mod.rs @@ -395,7 +395,7 @@ impl LightTrait for DistantLight { .normalize(); let p_outside = ctx.p() + wi * 2. * self.scene_radius; - let li = self.scale * self.lemit.sample(&lambda); + let li = self.scale * self.lemit.sample(lambda); let intr = SimpleInteraction::new( Point3fi::new_from_point(p_outside), 0.0, @@ -476,7 +476,7 @@ impl GoniometricLight { 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) + self.scale * self.iemit.sample(lambda) * self.image.lookup_nearest_channel(uv, 0) } } @@ -584,7 +584,7 @@ impl LightTrait for PointLight { .apply_to_interval(&Point3fi::default()); 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 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())); Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr))) } @@ -741,7 +741,7 @@ 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).into(); + rgb[c] = self.image.get_channel(Point2i::new(x, y), c); } let s = RGBIlluminantSpectrum::new( diff --git a/src/lights/sampler.rs b/src/lights/sampler.rs index c8a9ef0..1eefc1e 100644 --- a/src/lights/sampler.rs +++ b/src/lights/sampler.rs @@ -439,7 +439,7 @@ impl LightSamplerTrait for BVHLightSampler { } let mut node_pmf: Float = 0.; - let child = sample_discrete(&ci, u, Some(&mut node_pmf), Some(&mut u))?; + let child = sample_discrete(&ci, u, Some(&mut node_pmf), Some(&mut u)); pmf *= node_pmf; node_ind = if child == 0 { node_ind + 1 diff --git a/src/spectra/sampled.rs b/src/spectra/sampled.rs index 4355712..9ead79c 100644 --- a/src/spectra/sampled.rs +++ b/src/spectra/sampled.rs @@ -314,6 +314,8 @@ impl SampledWavelengths { } } + pub fn terminate_secondary_inplace(&mut self) {} + pub fn sample_uniform(u: Float, lambda_min: Float, lambda_max: Float) -> Self { let mut lambda = [0.0; N_SPECTRUM_SAMPLES]; lambda[0] = crate::core::pbrt::lerp(u, lambda_min, lambda_min); diff --git a/src/utils/containers.rs b/src/utils/containers.rs index 829577f..0572719 100644 --- a/src/utils/containers.rs +++ b/src/utils/containers.rs @@ -253,7 +253,7 @@ impl SampledGrid { where T: Interpolatable + Clone, { - self.lookup_convert(p, |v| v.clone()) + self.lookup_convert(p, |v| *v) } pub fn max_value_convert(&self, bounds: Bounds3f, convert: F) -> U diff --git a/src/utils/math.rs b/src/utils/math.rs index 1c4829e..22c2aaf 100644 --- a/src/utils/math.rs +++ b/src/utils/math.rs @@ -520,13 +520,13 @@ pub fn sample_discrete( u: f32, pmf: Option<&mut f32>, u_remapped: Option<&mut f32>, -) -> Option { +) -> usize { // Handle empty weights for discrete sampling. if weights.is_empty() { if let Some(p) = pmf { *p = 0.0; } - return None; + return 0; } // If the total weight is zero, sampling is not possible. @@ -535,7 +535,7 @@ pub fn sample_discrete( if let Some(p) = pmf { *p = 0.0; } - return None; + return 0; } let mut up = u * sum_weights; @@ -560,13 +560,12 @@ pub fn sample_discrete( *ur = ((up - sum) / weight).min(ONE_MINUS_EPSILON); } - Some(offset) + offset } pub fn sample_tent(u: Float, r: Float) -> Float { let mut u_remapped = 0.0; - let offset = sample_discrete(&[0.5, 0.5], u, None, Some(&mut u_remapped)) - .expect("Discrete sampling shouldn't fail"); + let offset = sample_discrete(&[0.5, 0.5], u, None, Some(&mut u_remapped)); if offset == 0 { -r + r * sample_linear(u, 0., 1.) } else { diff --git a/src/utils/rng.rs b/src/utils/rng.rs index 005e55a..4e5e558 100644 --- a/src/utils/rng.rs +++ b/src/utils/rng.rs @@ -23,6 +23,26 @@ impl Rng { self.uniform_u32(); } + pub fn set_sequence_with_offset(&mut self, sequence_index: u64, offset: u64) { + self.state = 0; + self.inc = (sequence_index << 1) | 1; + self.uniform_u32(); + self.state = self.state.wrapping_add(offset); + self.uniform_u32(); + } + + pub fn new(sequence_index: u64) -> Self { + let mut rng = Self { state: 0, inc: 0 }; + rng.set_sequence(sequence_index); + rng + } + + pub fn new_with_offset(sequence_index: u64, offset: u64) -> Self { + let mut rng = Self { state: 0, inc: 0 }; + rng.set_sequence_with_offset(sequence_index, offset); + rng + } + pub fn advance(&mut self, delta: u64) { // TODO: Implementation of PCG32 advance/seek for _ in 0..delta { diff --git a/src/utils/sampling.rs b/src/utils/sampling.rs index 6a890b4..369a746 100644 --- a/src/utils/sampling.rs +++ b/src/utils/sampling.rs @@ -13,6 +13,7 @@ use crate::utils::math::{ catmull_rom_weights, difference_of_products, logistic, newton_bisection, square, sum_of_products, }; +use crate::utils::rng::Rng; use std::sync::atomic::{AtomicU64, Ordering as SyncOrdering}; pub fn sample_linear(u: Float, a: Float, b: Float) -> Float { @@ -1533,3 +1534,118 @@ impl PiecewiseLinear2D { result } } + +#[derive(Clone, Debug)] +pub struct WeightedReservoirSampler { + rng: Rng, + weight_sum: Float, + reservoir_weight: Float, + reservoir: Option, +} + +impl Default for WeightedReservoirSampler { + fn default() -> Self { + Self { + rng: Rng::default(), + weight_sum: 0.0, + reservoir_weight: 0.0, + reservoir: None, + } + } +} + +impl WeightedReservoirSampler { + pub fn new(seed: u64) -> Self { + let mut rng = Rng::default(); + rng.set_sequence(seed); + Self { + rng, + weight_sum: 0.0, + reservoir_weight: 0.0, + reservoir: None, + } + } + + pub fn seed(&mut self, seed: u64) { + self.rng.set_sequence(seed); + } + + pub fn add(&mut self, sample: T, weight: Float) -> bool { + self.weight_sum += weight; + + let p = weight / self.weight_sum; + if self.rng.uniform::() < p { + self.reservoir = Some(sample); + self.reservoir_weight = weight; + return true; + } + + // debug_assert!(self.weight_sum < 1e80); + false + } + + pub fn add_closure(&mut self, func: F, weight: Float) -> bool + where + F: FnOnce() -> T, + { + self.weight_sum += weight; + + let p = weight / self.weight_sum; + if self.rng.uniform::() < p { + self.reservoir = Some(func()); + self.reservoir_weight = weight; + return true; + } + + // debug_assert!(self.weight_sum < 1e80); + false + } + + pub fn copy_from(&mut self, other: &Self) + where + T: Clone, + { + self.weight_sum = other.weight_sum; + self.reservoir = other.reservoir.clone(); + self.reservoir_weight = other.reservoir_weight; + } + + pub fn has_sample(&self) -> bool { + self.weight_sum > 0.0 + } + + pub fn get_sample(&self) -> Option<&T> { + self.reservoir.as_ref() + } + + pub fn sample_probability(&self) -> Float { + if self.weight_sum == 0.0 { + 0.0 + } else { + self.reservoir_weight / self.weight_sum + } + } + + pub fn weight_sum(&self) -> Float { + self.weight_sum + } + + pub fn reset(&mut self) { + self.reservoir_weight = 0.0; + self.weight_sum = 0.0; + self.reservoir = None; + } + + pub fn merge(&mut self, other: &WeightedReservoirSampler) + where + T: Clone, + { + // debug_assert!(self.weight_sum + other.weight_sum < 1e80); + + if let Some(other_sample) = &other.reservoir + && self.add(other_sample.clone(), other.weight_sum) + { + self.reservoir_weight = other.reservoir_weight; + } + } +}