From 63f4a36e69cfccf5629e12c044c93bd0aca8304f Mon Sep 17 00:00:00 2001 From: pingu Date: Tue, 9 Dec 2025 14:07:25 +0000 Subject: [PATCH] Updated light management and creation, refactored module. Started work on rendering and integrators, fixed up maths functions --- src/core/options.rs | 89 ++++ src/integrators/mod.rs | 520 ++++++++++++++++++++ src/integrators/pipeline.rs | 264 ++++++++++ src/lights/diffuse.rs | 260 +++++++++- src/lights/infinite.rs | 593 +++++++++++++++++++++++ src/lights/mod.rs | 937 +++++++++++++++++++++++++++++++++++- src/lights/sampler.rs | 521 ++++++++++++++++++++ 7 files changed, 3149 insertions(+), 35 deletions(-) create mode 100644 src/core/options.rs create mode 100644 src/integrators/mod.rs create mode 100644 src/integrators/pipeline.rs create mode 100644 src/lights/infinite.rs create mode 100644 src/lights/sampler.rs diff --git a/src/core/options.rs b/src/core/options.rs new file mode 100644 index 0000000..5a320b2 --- /dev/null +++ b/src/core/options.rs @@ -0,0 +1,89 @@ +use std::ops::Deref; +use std::sync::OnceLock; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum RenderingCoordinateSystem { + Camera, + CameraWorld, + World, +} + +#[derive(Debug, Clone, Copy)] +pub struct BasicPBRTOptions { + pub seed: i32, + pub quiet: bool, + pub disable_pixel_jitter: bool, + pub disable_wavelength_jitter: bool, + pub disable_texture_filtering: bool, + pub force_diffuse: bool, + pub use_gpu: bool, + pub wavefront: bool, + pub rendering_space: RenderingCoordinateSystem, +} + +impl Default for BasicPBRTOptions { + fn default() -> Self { + Self { + seed: 0, + quiet: false, + disable_pixel_jitter: false, + disable_wavelength_jitter: false, + disable_texture_filtering: false, + force_diffuse: false, + use_gpu: false, + wavefront: false, + rendering_space: RenderingCoordinateSystem::CameraWorld, + } + } +} + +#[derive(Debug, Clone)] +pub struct PBRTOptions { + pub basic: BasicPBRTOptions, + + pub n_threads: usize, + pub log_level: String, + pub image_file: String, +} + +impl Default for PBRTOptions { + fn default() -> Self { + Self { + basic: BasicPBRTOptions::default(), + n_threads: 0, // 0 usually implies "autodetect" + log_level: "info".to_string(), + image_file: "output.exr".to_string(), + } + } +} + +impl Deref for PBRTOptions { + type Target = BasicPBRTOptions; + + fn deref(&self) -> &Self::Target { + &self.basic + } +} + +// ----------------------------------------------------------------------- +// Global State Management +// ----------------------------------------------------------------------- + +static OPTIONS: OnceLock = OnceLock::new(); + +pub fn init_pbrt(options: PBRTOptions) { + OPTIONS + .set(options) + .expect("PBRT has already been initialized!"); +} + +pub fn cleanup_pbrt() { + todo!() +} + +pub fn get_options() -> &'static PBRTOptions { + OPTIONS.get().unwrap_or_else(|| { + static DEFAULT: OnceLock = OnceLock::new(); + DEFAULT.get_or_init(PBRTOptions::default) + }) +} diff --git a/src/integrators/mod.rs b/src/integrators/mod.rs new file mode 100644 index 0000000..6125927 --- /dev/null +++ b/src/integrators/mod.rs @@ -0,0 +1,520 @@ +mod pipeline; + +use pipeline::*; + +use crate::camera::{Camera, CameraTrait}; +use crate::core::bxdf::{BSDF, BxDFFlags, BxDFTrait, FArgs, TransportMode}; +use crate::core::film::{FilmTrait, VisibleSurface}; +use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction}; +use crate::core::options::get_options; +use crate::core::pbrt::{Float, SHADOW_EPSILON}; +use crate::core::primitive::{Primitive, PrimitiveTrait}; +use crate::core::sampler::{CameraSample, Sampler, SamplerTrait}; +use crate::geometry::{Bounds2i, Point2f, Point2i, Point3fi, Ray, Vector3f, VectorLike}; +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::sampling::{ + power_heuristic, sample_uniform_hemisphere, sample_uniform_sphere, uniform_hemisphere_pdf, + uniform_sphere_pdf, +}; +use bumpalo::Bump; +use rayon::prelude::*; + +use std::sync::Arc; + +#[derive(Clone, Debug)] +pub struct IntegratorBase { + pub aggregate: Arc, + pub lights: Vec>, + pub infinite_lights: Vec>, +} + +impl IntegratorBase { + pub fn new(aggregate: Arc, lights: Vec>) -> Self { + let infinite_lights = lights + .iter() + .filter(|light| light.light_type().is_infinite()) + .cloned() + .collect(); + + Self { + aggregate, + lights, + infinite_lights, + } + } + + pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + self.aggregate.intersect(ray, t_max) + } + + pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + self.aggregate.intersect_p(ray, t_max) + } + + pub fn unoccluded(&self, p0: &Interaction, p1: &Interaction) -> bool { + !self.intersect_p(&p0.spawn_ray_to_interaction(p1), Some(1. - SHADOW_EPSILON)) + } +} + +pub trait IntegratorTrait { + fn render(&self); +} + +pub trait RayIntegratorTrait: Send + Sync + std::fmt::Debug { + fn evaluate_pixel_sample( + &self, + p_pixel: Point2i, + sample_ind: usize, + sampler: &mut Sampler, + scratch: &Bump, + ); + + fn li( + &self, + ray: Ray, + lambda: &SampledWavelengths, + sampler: &mut Sampler, + scratch: &Bump, + visible_surface: bool, + ) -> (SampledSpectrum, Option); +} + +struct NextRaySample { + wi: Vector3f, + f: SampledSpectrum, + pdf: Float, + is_specular: bool, +} + +#[derive(Clone, Debug)] +pub struct SimplePathIntegrator { + base: IntegratorBase, + camera: Arc, + light_sampler: UniformLightSampler, + sample_lights: bool, + sample_bsdf: bool, + max_depth: usize, +} + +impl SimplePathIntegrator { + fn sample_direct_lighting( + &self, + isect: &SurfaceInteraction, + bsdf: &BSDF, + sampler: &mut Sampler, + lambda: &SampledWavelengths, + ) -> SampledSpectrum { + let Some(sampled_light) = self.light_sampler.sample(sampler.get1d()) else { + return SampledSpectrum::new(0.0); + }; + let u_light = sampler.get2d(); + let p0 = LightSampleContext::from(isect); + + let Some(ls) = sampled_light.light.sample_li(&p0, u_light, lambda, false) else { + return SampledSpectrum::new(0.0); + }; + + if ls.pdf == 0.0 || ls.l.is_black() { + return SampledSpectrum::new(0.0); + } + + let wi = ls.wi; + let wo = -isect.wo(); + + let Some(mut f) = bsdf.f(wo, wi, TransportMode::Radiance) else { + return SampledSpectrum::new(0.); + }; + + f *= wi.abs_dot(isect.shading.n.into()); + + if f.is_black() { + return SampledSpectrum::new(0.0); + } + + // Visibility check + let p0_surf = Interaction::Surface(isect.clone()); + if !self.base.unoccluded(&p0_surf, ls.p_light.as_ref()) { + return SampledSpectrum::new(0.0); + } + + f * ls.l / (sampled_light.p * ls.pdf) + } + + fn sample_bsdf_direction( + &self, + wo: Vector3f, + bsdf: &BSDF, + sampler: &mut Sampler, + _isect: &SurfaceInteraction, + ) -> Option { + let u = sampler.get1d(); + let bs = bsdf.sample_f(wo, u, sampler.get2d(), FArgs::default())?; + + Some(NextRaySample { + wi: bs.wi, + f: bs.f, + pdf: bs.pdf, + is_specular: bs.is_specular(), + }) + } + + fn sample_uniform_direction( + &self, + wo: Vector3f, + bsdf: &BSDF, + sampler: &mut Sampler, + isect: &SurfaceInteraction, + ) -> Option { + let flags = bsdf.flags(); + let mut wi; + let pdf; + + if BxDFFlags::is_reflective(&flags) && BxDFFlags::is_transmissive(&flags) { + wi = sample_uniform_sphere(sampler.get2d()); + pdf = uniform_sphere_pdf(); + } else { + wi = sample_uniform_hemisphere(sampler.get2d()); + pdf = uniform_hemisphere_pdf(); + + let same_hemi = wo.dot(isect.n().into()) * wi.dot(isect.n().into()) > 0.0; + if BxDFFlags::is_reflective(&flags) && !same_hemi + || BxDFFlags::is_transmissive(&flags) && same_hemi + { + wi = -wi; + } + } + + let f = bsdf.f(wo, wi, TransportMode::Radiance)?; + + Some(NextRaySample { + wi, + f, + pdf, + is_specular: false, + }) + } +} + +impl RayIntegratorTrait for SimplePathIntegrator { + 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.0); + let mut beta = SampledSpectrum::new(1.0); + let mut specular_bounce = true; + let mut depth = 0; + while !beta.is_black() { + let Some(mut si) = self.base.intersect(&ray, None) else { + if !self.sample_lights || specular_bounce { + for light in &self.base.infinite_lights { + l += light.le(&ray, lambda); + } + } + break; + }; + + let t_hit = si.t_hit(); + let isect = &mut si.intr; + if !self.sample_lights || specular_bounce { + l += beta * isect.le(-ray.d, lambda); + } + depth += 1; + if depth == self.max_depth { + break; + } + + let Some(bsdf) = isect.get_bsdf(&ray, lambda, self.camera.as_ref(), scratch, sampler) + else { + specular_bounce = false; + isect.skip_intersection(&mut ray, t_hit); + continue; + }; + + let wo = -ray.d; + if self.sample_lights { + l += beta * self.sample_direct_lighting(isect, &bsdf, sampler, lambda); + } + + let sample = if self.sample_bsdf { + self.sample_bsdf_direction(wo, &bsdf, sampler, isect) + } else { + self.sample_uniform_direction(wo, &bsdf, sampler, isect) + }; + + let Some(bs) = sample else { + break; + }; + + beta *= bs.f * bs.wi.abs_dot(isect.shading.n.into()) / bs.pdf; + specular_bounce = bs.is_specular; + ray = isect.spawn_ray(bs.wi); + } + assert!(beta.y(lambda) > 0.); + debug_assert!(beta.y(lambda).is_finite()); + (l, None) + } +} + +#[derive(Clone, Debug)] +pub struct PathIntegrator { + base: IntegratorBase, + camera: Arc, + sampler_prototype: Sampler, + light_sampler: LightSampler, + regularize: bool, + max_depth: usize, +} + +impl PathIntegrator { + fn sample_ld( + &self, + intr: &SurfaceInteraction, + bsdf: &BSDF, + lambda: &SampledWavelengths, + sampler: &mut Sampler, + ) -> Option { + let mut ctx: LightSampleContext = intr.into(); + let flags = bsdf.flags(); + if BxDFFlags::is_reflective(&flags) ^ !BxDFFlags::is_transmissive(&flags) { + ctx.pi = Point3fi::new_from_point(intr.offset_ray_vector(-intr.wo())); + } + + let u = sampler.get1d(); + let sampled_light = self.light_sampler.sample_with_context(&ctx, u)?; + + let u_light = sampler.get2d(); + let light = sampled_light.light; + let ls = light.sample_li(&ctx, u_light, lambda, true)?; + + if ls.l.is_black() || ls.pdf == 0.0 { + return None; + } + + let wo = intr.wo(); + let wi = ls.wi; + let f = bsdf.f(wo, wi, TransportMode::Radiance)? * wi.abs_dot(intr.shading.n.into()); + + if !self + .base + .unoccluded(&Interaction::Surface(intr.clone()), &ls.p_light) + { + return None; + } + + let pl = sampled_light.p * ls.pdf; + if light.light_type().is_delta_light() { + Some(ls.l * f / pl) + } else { + let pb = bsdf.pdf(wo, wi, FArgs::default()); + let wl = power_heuristic(1, pl, 1, pb); + Some(wl * ls.l * f / pl) + } + } +} + +impl RayIntegratorTrait for PathIntegrator { + 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.0); + let mut beta = SampledSpectrum::new(1.0); + let mut visible: Option = None; + let mut depth = 0; + let mut pb = 1.; + let mut eta_scale = 1.; + let mut specular_bounce = false; + let mut any_non_specular_bounces = false; + + let mut prev_int_ctx = LightSampleContext::default(); + loop { + let Some(mut si) = self.base.intersect(&ray, None) else { + for light in &self.base.infinite_lights { + let le = light.le(&ray, lambda); + if depth == 0 || specular_bounce { + l += beta * le; + } else { + let pl = self.light_sampler.pmf_with_context(&prev_int_ctx, light) + * light.pdf_li(&prev_int_ctx, ray.d, true); + let w_b = power_heuristic(1, pb, 1, pl); + l += beta * w_b * le; + } + } + break; + }; + + let t_hit = si.t_hit(); + let isect = &mut si.intr; + let le = isect.le(-ray.d, lambda); + if depth == 0 || specular_bounce { + l += beta * le; + } else if let Some(area_light) = &isect.area_light { + let pl = self + .light_sampler + .pmf_with_context(&prev_int_ctx, area_light) + + area_light.pdf_li(&prev_int_ctx, ray.d, true); + let wl = power_heuristic(1, pb, 1, pl); + l += beta * wl * le; + } + + let Some(mut bsdf) = + isect.get_bsdf(&ray, lambda, self.camera.as_ref(), scratch, sampler) + else { + specular_bounce = true; + isect.skip_intersection(&mut ray, 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 self.regularize && any_non_specular_bounces { + 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; + } + } + + let wo = -ray.d; + 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(isect.shading.n.into()) / bs.pdf; + pb = if bs.pdf_is_proportional { + bsdf.pdf(wo, bs.wi, FArgs::default()) + } else { + bs.pdf + }; + specular_bounce = bs.is_specular(); + any_non_specular_bounces |= !bs.is_specular(); + if bs.is_transmissive() { + eta_scale *= square(bs.eta); + } + prev_int_ctx = LightSampleContext::from(&si.intr); + ray = si + .intr + .spawn_ray_with_differentials(&ray, bs.wi, bs.flags, bs.eta); + + let rr_beta = beta * eta_scale; + if rr_beta.max_component_value() < 1. && depth > 1 { + let q = (1. - rr_beta.max_component_value()).max(0.); + if sampler.get1d() < q { + break; + } + beta /= 1. - q; + } + } + + (l, visible) + } +} + +pub enum Integrator { + BPDT(BPDTIntegrator), + MLT(MLTIntegrator), + SPPM(SPPMIntegrator), + Sampler(SamplerIntegrator), +} + +pub struct BPDTIntegrator; +pub struct MLTIntegrator; +pub struct SPPMIntegrator; +pub struct SamplerIntegrator; diff --git a/src/integrators/pipeline.rs b/src/integrators/pipeline.rs new file mode 100644 index 0000000..61fd54b --- /dev/null +++ b/src/integrators/pipeline.rs @@ -0,0 +1,264 @@ +use crate::core::{options::PBRTOptions, sampler::get_camera_sample}; +use crate::image::{Image, ImageMetadata}; +use indicatif::{ProgressBar, ProgressStyle}; +use std::io::Write; +use std::path::Path; + +use super::*; + +struct PbrtProgress { + bar: ProgressBar, +} + +impl PbrtProgress { + fn new(total_work: u64, description: &str, quiet: bool) -> Self { + if quiet { + return Self { + bar: ProgressBar::hidden(), + }; + } + + let bar = ProgressBar::new(total_work); + + bar.set_style( + ProgressStyle::default_bar() + .template("[{elapsed_precise}] {bar:40.cyan/blue} {pos:>7}/{len:7} {msg}") + .unwrap() + .progress_chars("=>-"), + ); + + bar.set_message(description.to_string()); + + Self { bar } + } + + fn update(&self, amount: u64) { + self.bar.inc(amount); + } + + fn done(&self) { + self.bar.finish_with_message("Done"); + } + + fn elapsed_seconds(&self) -> f32 { + self.bar.elapsed().as_secs_f32() + } +} + +fn generate_tiles(bounds: Bounds2i) -> Vec { + let mut tiles = Vec::new(); + const TILE_SIZE: i32 = 16; + for y in (bounds.p_min.y()..bounds.p_max.y()).step_by(TILE_SIZE as usize) { + for x in (bounds.p_min.x()..bounds.p_max.x()).step_by(TILE_SIZE as usize) { + let p_min = Point2i::new(x, y); + let p_max = Point2i::new( + (x + TILE_SIZE).min(bounds.p_max.x()), + (y + TILE_SIZE).min(bounds.p_max.y()), + ); + tiles.push(Bounds2i::from_points(p_min, p_max)); + } + } + tiles +} + +pub fn render( + integrator: &T, + _base: &IntegratorBase, + camera: &Camera, + sampler_prototype: &Sampler, +) where + T: RayIntegratorTrait, +{ + let options = get_options(); + if let Some((p_pixel, sample_index)) = options.debug_start { + let s_index = sample_index as usize; + let scratch = Bump::new(); + let mut tile_sampler = sampler_prototype.clone(); + + tile_sampler.start_pixel_sample(p_pixel, s_index, None); + + evaluate_pixel_sample( + integrator, + camera, + &mut tile_sampler, + p_pixel, + s_index, + &scratch, + ); + return; + } + + let pixel_bounds = camera.get_film().pixel_bounds(); + let spp = sampler_prototype.samples_per_pixel(); + let total_work = (pixel_bounds.area() as u64) * (spp as u64); + let progress = PbrtProgress::new(total_work, "Rendering", options.quiet); + let mut wave_start = 0; + let mut wave_end = 1; + let mut next_wave_size = 1; + + let mut reference_image: Option = None; + let mut mse_out_file: Option = None; + + if let Some(ref_path) = &options.mse_reference_image { + let image_and_metadata = + Image::read(Path::new(&ref_path), None).expect("Could not load image"); + let image = image_and_metadata.image; + let metadata = image_and_metadata.metadata; + let resolution = image.resolution(); + // reference_image = Some(image); + let mse_pixel_bounds = metadata + .pixel_bounds + .unwrap_or_else(|| Bounds2i::from_points(Point2i::new(0, 0), resolution)); + + if !mse_pixel_bounds.overlaps(&pixel_bounds) { + panic!("Output pixel bounds dont fit inside the reference image"); + } + + let crop_p_min = Point2i::from(pixel_bounds.p_min - mse_pixel_bounds.p_min); + let crop_p_max = Point2i::from(pixel_bounds.p_max - mse_pixel_bounds.p_min); + + let crop_bounds = Bounds2i::from_points(crop_p_min, crop_p_max); + + let cropped_image = image.crop(crop_bounds).clone(); + let cropped_resolution = cropped_image.resolution(); + + let expected_res = Point2i::new( + pixel_bounds.p_max.x() - pixel_bounds.p_min.x(), + pixel_bounds.p_max.y() - pixel_bounds.p_min.y(), + ); + + reference_image = Some(cropped_image); + + assert_eq!( + cropped_resolution, expected_res, + "Cropped reference image resolution mismatch" + ); + + if let Some(out_path) = &options.mse_reference_output { + mse_out_file = Some( + std::fs::File::create(out_path) + .expect(&format!("Failed to create MSE output file: {}", out_path)), + ); + } + } + + let tiles = generate_tiles(pixel_bounds); + while wave_start < spp { + tiles.par_iter().for_each(|tile_bounds| { + let mut arena = Bump::with_capacity(65 * 1024); + let mut sampler = sampler_prototype.clone(); + + for p_pixel in tile_bounds { + for sample_index in wave_start..wave_end { + sampler.start_pixel_sample(*p_pixel, sample_index, None); + + evaluate_pixel_sample( + integrator, + camera, + &mut sampler, + *p_pixel, + sample_index, + &arena, + ); + + arena.reset(); + } + } + + let work_done = (tile_bounds.area() as u64) * ((wave_end - wave_start) as u64); + progress.update(work_done); + }); + + wave_start = wave_end; + wave_end = (wave_end + next_wave_size).min(spp); + + if reference_image.is_none() { + next_wave_size = (2 * next_wave_size).min(64); + } + + if wave_start == spp { + progress.done(); + } + + if wave_start == spp || options.write_partial_images || reference_image.is_some() { + let mut metadata = ImageMetadata { + render_time_seconds: Some(progress.elapsed_seconds()), + samples_per_pixel: Some(wave_start as i32), + ..Default::default() + }; + + if wave_start == spp || options.write_partial_images { + camera.init_metadata(&mut metadata); + camera + .get_film() + .write_image(&metadata, 1.0 / wave_start as Float); + } + + if let Some(ref_img) = &reference_image { + let splat_scale = 1.0 / (wave_start as Float); + + let film_metadata = ImageMetadata::default(); + let film_image = camera.get_film().get_image(&film_metadata, splat_scale); + + let (mse_values, _mse_debug_img) = + film_image.mse(film_image.all_channels_desc(), ref_img, false); + let mse_avg = mse_values.average(); + + if let Some(file) = &mut mse_out_file { + writeln!(file, "{}, {:.9}", wave_start, mse_avg).ok(); + file.flush().ok(); + } + + metadata.mse = Some(mse_avg); + } + } + } +} + +pub fn evaluate_pixel_sample( + integrator: &T, + camera: &Camera, + sampler: &mut Sampler, + pixel: Point2i, + _sample_index: usize, + scratch: &Bump, +) { + let mut lu = sampler.get1d(); + if get_options().disable_wavelength_jitter { + lu = 0.5; + } + let lambda = camera.get_film().sample_wavelengths(lu); + let mut film = camera.get_film(); + let filter = film.get_filter(); + let camera_sample = get_camera_sample(sampler, pixel, filter); + if let Some(mut camera_ray) = camera.generate_ray_differential(camera_sample, &lambda) { + debug_assert!(camera_ray.ray.d.norm() > 0.999); + debug_assert!(camera_ray.ray.d.norm() < 1.001); + let ray_diff_scale = (sampler.samples_per_pixel() as Float).sqrt().max(0.125); + if get_options().disable_pixel_jitter { + camera_ray.ray.scale_differentials(ray_diff_scale); + } + + let initialize_visible_surface = film.uses_visible_surface(); + let (mut l, visible_surface) = integrator.li( + camera_ray.ray, + &lambda, + sampler, + scratch, + initialize_visible_surface, + ); + l *= camera_ray.weight; + + if l.has_nans() || l.y(&lambda).is_infinite() { + l = SampledSpectrum::new(0.); + } + + film.add_sample( + pixel, + l, + &lambda, + visible_surface.as_ref(), + camera_sample.filter_weight, + ); + } +} diff --git a/src/lights/diffuse.rs b/src/lights/diffuse.rs index a9c15f3..a5b7183 100644 --- a/src/lights/diffuse.rs +++ b/src/lights/diffuse.rs @@ -1,19 +1,255 @@ +use super::{ + DenselySampledSpectrum, LightBase, LightBounds, LightLiSample, LightSampleContext, LightTrait, + LightType, RGBIlluminantSpectrum, SampledSpectrum, SampledWavelengths, Spectrum, SpectrumTrait, +}; +use crate::core::interaction::{ + Interaction, InteractionTrait, SimpleInteraction, SurfaceInteraction, +}; use crate::core::medium::MediumInterface; -use crate::core::pbrt::Float; -use crate::shapes::Shape; -use crate::spectra::Spectrum; +use crate::core::pbrt::{Float, PI}; +use crate::core::texture::{ + FloatTexture, FloatTextureTrait, TextureEvalContext, TextureEvaluator, + UniversalTextureEvaluator, +}; +use crate::geometry::{ + Bounds3f, Normal3f, Point2f, Point2fi, Point2i, Point3f, Point3fi, Ray, Vector3f, VectorLike, +}; +use crate::image::Image; +use crate::shapes::{Shape, ShapeSample, ShapeSampleContext, ShapeTrait}; +use crate::utils::color::RGB; +use crate::utils::colorspace::RGBColorSpace; +use crate::utils::hash::hash_float; use crate::utils::transform::Transform; use std::sync::Arc; +#[derive(Clone, Debug)] pub struct DiffuseAreaLight { - pub l_emit: Spectrum, - pub shape: Arc, - pub two_sided: bool, - pub area: Float, - pub flags: u8, - pub n_samples: i32, - pub medium_interface: MediumInterface, - light_to_world: Transform, - world_to_light: Transform, + base: LightBase, + shape: Shape, + alpha: Option, + area: Float, + two_sided: bool, + lemit: Arc, + scale: Float, + image: Option, + image_color_space: Option>, +} + +impl DiffuseAreaLight { + #[allow(clippy::too_many_arguments)] + pub fn new( + render_from_light: Transform, + medium_interface: MediumInterface, + le: Spectrum, + scale: Float, + shape: Shape, + alpha: FloatTexture, + image: Option, + image_color_space: Option>, + two_sided: bool, + ) -> Self { + let is_constant_zero = match &alpha { + FloatTexture::FloatConstant(tex) => tex.evaluate(&TextureEvalContext::default()) == 0.0, + _ => false, + }; + + let (light_type, stored_alpha) = if is_constant_zero { + (LightType::DeltaPosition, None) + } else { + (LightType::Area, Some(alpha)) + }; + + let base = LightBase::new(light_type, &render_from_light, &medium_interface); + + let lemit = LightBase::lookup_spectrum(&le); + if let Some(im) = &image { + let desc = im + .get_channel_desc(&["R", "G", "B"]) + .expect("Image used for DiffuseAreaLight doesn't have R, G, B channels"); + + assert_eq!(3, desc.size(), "Image channel description size mismatch"); + assert!( + desc.is_identity(), + "Image channel description is not identity" + ); + + assert!( + image_color_space.is_some(), + "Image provided but ColorSpace is missing" + ); + } + let is_triangle_or_bilinear = matches!(shape, Shape::Triangle(_) | Shape::BilinearPatch(_)); + if render_from_light.has_scale(None) && !is_triangle_or_bilinear { + println!( + "Scaling detected in rendering to light space transformation! \ + The system has numerous assumptions, implicit and explicit, \ + that this transform will have no scale factors in it. \ + Proceed at your own risk; your image may have errors." + ); + } + + Self { + base, + area: shape.area(), + shape, + alpha: stored_alpha, + two_sided, + lemit, + scale, + image, + image_color_space, + } + } + + fn alpha_masked(&self, intr: &Interaction) -> bool { + let Some(alpha_tex) = &self.alpha else { + return false; + }; + let ctx = TextureEvalContext::from(intr); + let a = UniversalTextureEvaluator.evaluate_float(alpha_tex, &ctx); + if a >= 1.0 { + return false; + } + if a <= 0.0 { + return true; + } + hash_float(&intr.p()) > a + } +} + +impl LightTrait for DiffuseAreaLight { + fn base(&self) -> &LightBase { + &self.base + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + let mut l = SampledSpectrum::new(0.); + if let Some(image) = &self.image { + for y in 0..image.resolution().y() { + for x in 0..image.resolution().x() { + let mut rgb = RGB::default(); + for c in 0..3 { + rgb[c] = image.get_channel(Point2i::new(x, y), c); + } + l += RGBIlluminantSpectrum::new( + self.image_color_space.as_ref().unwrap(), + RGB::clamp_zero(rgb), + ) + .sample(&lambda); + } + } + l *= self.scale / (image.resolution().x() * image.resolution().y()) as Float; + } else { + l = self.lemit.sample(&lambda) * self.scale; + } + let two_side = if self.two_sided { 2. } else { 1. }; + PI * two_side * self.area * l + } + + fn sample_li( + &self, + ctx: &LightSampleContext, + u: Point2f, + lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + let shape_ctx = ShapeSampleContext::new(ctx.pi, ctx.n, ctx.ns, 0.0); + let ss = self.shape.sample_from_context(&shape_ctx, u)?; + let mut intr: SurfaceInteraction = ss.intr.as_ref().clone(); + + intr.common.medium_interface = Some(self.base.medium_interface.clone()); + let p = intr.p(); + let n = intr.n(); + let uv = intr.uv; + + let generic_intr = Interaction::Surface(intr); + if self.alpha_masked(&generic_intr) { + return None; + } + + let wi = (p - ctx.p()).normalize(); + let le = self.l(p, n, uv, -wi, &lambda); + + if le.is_black() { + return None; + } + + Some(LightLiSample::new(le, wi, ss.pdf, generic_intr)) + } + + fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, _allow_incomplete_pdf: bool) -> Float { + let shape_ctx = ShapeSampleContext::new(ctx.pi, ctx.n, ctx.ns, 0.); + self.shape.pdf_from_context(&shape_ctx, wi) + } + + fn l( + &self, + p: Point3f, + n: Normal3f, + mut uv: Point2f, + w: Vector3f, + lambda: &SampledWavelengths, + ) -> SampledSpectrum { + if self.two_sided && n.dot(w.into()) < 0. { + return SampledSpectrum::new(0.); + } + let intr = Interaction::Surface(SurfaceInteraction::new_minimal( + Point3fi::new_from_point(p), + uv, + )); + + if self.alpha_masked(&intr) { + return SampledSpectrum::new(0.); + } + if let Some(image) = &self.image { + let mut rgb = RGB::default(); + uv[1] = 1. - uv[1]; + for c in 0..3 { + rgb[c] = image.bilerp_channel(uv, c); + } + + let spec = RGBIlluminantSpectrum::new( + self.image_color_space.as_ref().unwrap(), + RGB::clamp_zero(rgb), + ); + + self.scale * spec.sample(lambda) + } else { + self.scale * self.lemit.sample(lambda) + } + } + + fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + todo!() + } + + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + unimplemented!() + } + + fn bounds(&self) -> Option { + let mut phi = 0.; + if let Some(image) = &self.image { + for y in 0..image.resolution.y() { + for x in 0..image.resolution.x() { + for c in 0..3 { + phi += image.get_channel(Point2i::new(x, y), c); + } + } + } + } else { + phi = self.lemit.max_value(); + } + + let nb = self.shape.normal_bounds(); + Some(LightBounds::new( + &self.shape.bounds(), + nb.w, + phi, + nb.cos_theta, + (PI / 2.).cos(), + self.two_sided, + )) + } } diff --git a/src/lights/infinite.rs b/src/lights/infinite.rs new file mode 100644 index 0000000..0de25a7 --- /dev/null +++ b/src/lights/infinite.rs @@ -0,0 +1,593 @@ +use crate::{ + core::medium::Medium, + geometry::Frame, + spectra::RGBIlluminantSpectrum, + utils::{ + color::RGB, + colorspace::RGBColorSpace, + math::equal_area_square_to_sphere, + sampling::{ + AliasTable, PiecewiseConstant2D, WindowedPiecewiseConstant2D, sample_uniform_sphere, + uniform_sphere_pdf, + }, + }, +}; + +use rayon::prelude::*; + +use crate::core::pbrt::clamp_t; +use crate::image::{PixelFormat, WrapMode}; + +use std::path::Path; + +use super::{ + Arc, Bounds2f, Bounds3f, DenselySampledSpectrum, Float, Image, Interaction, LightBase, + LightBounds, LightLiSample, LightSampleContext, LightTrait, LightType, MediumInterface, + Normal3f, PI, Point2f, Point2i, Point3f, Ray, SampledSpectrum, SampledWavelengths, + SimpleInteraction, Spectrum, Transform, Vector3f, VectorLike, equal_area_sphere_to_square, + square, +}; + +#[derive(Debug, Clone)] +pub struct InfiniteUniformLight { + base: LightBase, + lemit: Arc, + scale: Float, + scene_center: Point3f, + scene_radius: Float, +} + +impl InfiniteUniformLight { + pub fn new(render_from_light: Transform, le: Spectrum, scale: Float) -> Self { + let base = LightBase::new( + LightType::Infinite, + &render_from_light, + &MediumInterface::default(), + ); + let lemit = LightBase::lookup_spectrum(&le); + Self { + base, + lemit, + scale, + scene_center: Point3f::default(), + scene_radius: 0., + } + } +} + +impl LightTrait for InfiniteUniformLight { + fn base(&self) -> &LightBase { + &self.base + } + fn sample_li( + &self, + ctx: &LightSampleContext, + u: Point2f, + lambda: &SampledWavelengths, + allow_incomplete_pdf: bool, + ) -> Option { + if allow_incomplete_pdf { + return None; + } + let wi = sample_uniform_sphere(u); + let pdf = uniform_sphere_pdf(); + let intr_simple = SimpleInteraction::new_interface( + ctx.p() + wi * (2. * self.scene_radius), + Some(MediumInterface::default()), + ); + + let intr = Interaction::Simple(intr_simple); + Some(LightLiSample::new( + self.scale * self.lemit.sample(&lambda), + wi, + pdf, + intr, + )) + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + 4. * PI * PI * square(self.scene_radius) * self.scale * self.lemit.sample(&lambda) + } + + fn pdf_li( + &self, + _ctx: &LightSampleContext, + _wi: Vector3f, + allow_incomplete_pdf: bool, + ) -> Float { + if allow_incomplete_pdf { + return 0.; + } + uniform_sphere_pdf() + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + + fn le(&self, _ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum { + self.scale * self.lemit.sample(lambda) + } + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + todo!() + } + fn bounds(&self) -> Option { + todo!() + } +} + +#[derive(Clone, Debug)] +pub struct InfiniteImageLight { + base: LightBase, + image: Image, + image_color_space: Arc, + scale: Float, + scene_radius: Float, + scene_center: Point3f, + distrib: PiecewiseConstant2D, + compensated_distrib: PiecewiseConstant2D, +} + +impl InfiniteImageLight { + pub fn new( + render_from_light: Transform, + image: Image, + image_color_space: Arc, + scale: Float, + filename: String, + ) -> Self { + let base = LightBase::new( + LightType::Infinite, + &render_from_light, + &MediumInterface::default(), + ); + + let desc = image + .get_channel_desc(&["R", "G", "B"]) + .expect("Image used for DiffuseAreaLight doesn't have R, G, B channels"); + + assert_eq!(3, desc.size()); + assert!(desc.is_identity()); + if image.resolution().x() != image.resolution().y() { + panic!( + "{}: image resolution ({}, {}) is non-square. It's unlikely this is an equal area environment map.", + filename, + image.resolution.x(), + image.resolution.y() + ); + } + let mut d = image.get_sampling_distribution_uniform(); + let domain = Bounds2f::from_points(Point2f::new(0., 0.), Point2f::new(1., 1.)); + let distrib = PiecewiseConstant2D::new_with_bounds(&d, domain); + let slice = &mut d.values; // or d.as_slice_mut() + let count = slice.len() as Float; + let sum: Float = slice.iter().sum(); + let average = sum / count; + + for v in slice.iter_mut() { + *v = (*v - average).max(0.0); + } + + let all_zero = slice.iter().all(|&v| v == 0.0); + if all_zero { + for v in slice.iter_mut() { + *v = 1.0; + } + } + + let compensated_distrib = PiecewiseConstant2D::new_with_bounds(&d, domain); + + Self { + base, + image, + image_color_space, + scene_center: Point3f::default(), + scene_radius: 0., + scale, + distrib, + compensated_distrib, + } + } + + fn image_le(&self, uv: Point2f, lambda: &SampledWavelengths) -> SampledSpectrum { + let mut rgb = RGB::default(); + for c in 0..3 { + rgb[c] = self.image.lookup_nearest_channel_with_wrap( + uv, + c, + WrapMode::OctahedralSphere.into(), + ); + } + let spec = + RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb)); + self.scale * spec.sample(lambda) + } +} + +impl LightTrait for InfiniteImageLight { + fn base(&self) -> &LightBase { + &self.base + } + fn sample_li( + &self, + ctx: &LightSampleContext, + u: Point2f, + lambda: &SampledWavelengths, + allow_incomplete_pdf: bool, + ) -> Option { + let (uv, map_pdf, _) = if allow_incomplete_pdf { + self.compensated_distrib.sample(u) + } else { + self.distrib.sample(u) + }; + + if map_pdf == 0. { + return None; + } + // Convert infinite light sample point to direction + let w_light = equal_area_square_to_sphere(uv); + let wi = self.base.render_from_light.apply_to_vector(w_light); + let pdf = map_pdf / (4. * PI); + + // Return radiance value for infinite light direction + let mut simple_intr = SimpleInteraction::new_interface( + ctx.p() + wi * (2. * self.scene_radius), + Some(MediumInterface::default()), + ); + + 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, + )) + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + let mut sum_l = SampledSpectrum::new(0.); + let width = self.image.resolution.x(); + let height = self.image.resolution.y(); + for v in 0..height { + for u in 0..width { + let mut rgb = RGB::default(); + for c in 0..3 { + rgb[c] = self.image.get_channel_with_wrap( + Point2i::new(u, v), + c, + WrapMode::OctahedralSphere.into(), + ); + } + sum_l += RGBIlluminantSpectrum::new( + self.image_color_space.as_ref(), + RGB::clamp_zero(rgb), + ) + .sample(&lambda); + } + } + 4. * PI * PI * square(self.scene_radius) * self.scale * sum_l / (width * height) as Float + } + + fn pdf_li(&self, _ctx: &LightSampleContext, wi: Vector3f, allow_incomplete_pdf: bool) -> Float { + let w_light = self.base.render_from_light.apply_inverse_vector(wi); + let uv = equal_area_sphere_to_square(w_light); + let pdf = if allow_incomplete_pdf { + self.compensated_distrib.pdf(uv) + } else { + self.distrib.pdf(uv) + }; + pdf / (4. * PI) + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + + fn le(&self, ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum { + let w_light = self + .base + .render_from_light + .apply_inverse_vector(ray.d) + .normalize(); + let uv = equal_area_sphere_to_square(w_light); + self.image_le(uv, lambda) + } + + fn preprocess(&mut self, scene_bounds: &Bounds3f) { + let (scene_center, scene_radius) = scene_bounds.bounding_sphere(); + self.scene_center = scene_center; + self.scene_radius = scene_radius; + } + + fn bounds(&self) -> Option { + None + } +} + +#[derive(Debug, Clone)] +pub struct InfinitePortalLight { + pub base: LightBase, + pub image: Image, + pub image_color_space: Arc, + pub scale: Float, + pub filename: String, + pub portal: [Point3f; 4], + pub portal_frame: Frame, + pub distribution: WindowedPiecewiseConstant2D, + pub scene_center: Point3f, + pub scene_radius: Float, +} + +impl InfinitePortalLight { + fn base(&self) -> &LightBase { + &self.base + } + pub fn new( + render_from_light: Transform, + equal_area_image: &Image, + image_color_space: Arc, + scale: Float, + filename: String, + points: Vec, + ) -> Self { + let base = LightBase::new( + LightType::Infinite, + &render_from_light, + &MediumInterface::default(), + ); + + let desc = equal_area_image + .get_channel_desc(&["R", "G", "B"]) + .unwrap_or_else(|_| { + panic!( + "{}: image used for PortalImageInfiniteLight doesn't have R, G, B channels.", + filename + ) + }); + + assert_eq!(3, desc.offset.len()); + let src_res = equal_area_image.resolution; + if src_res.x() != src_res.y() { + panic!( + "{}: image resolution ({}, {}) is non-square. It's unlikely this is an equal area environment map.", + filename, + src_res.x(), + src_res.y() + ); + } + + if points.len() != 4 { + panic!( + "Expected 4 vertices for infinite light portal but given {}", + points.len() + ); + } + + let portal: [Point3f; 4] = [points[0], points[1], points[2], points[3]]; + + let p01 = (portal[1] - portal[0]).normalize(); + let p12 = (portal[2] - portal[1]).normalize(); + let p32 = (portal[2] - portal[3]).normalize(); + let p03 = (portal[3] - portal[0]).normalize(); + + if (p01.dot(p32) - 1.0).abs() > 0.001 || (p12.dot(p03) - 1.0).abs() > 0.001 { + panic!("Infinite light portal isn't a planar quadrilateral (opposite edges)"); + } + + if p01.dot(p12).abs() > 0.001 + || p12.dot(p32).abs() > 0.001 + || p32.dot(p03).abs() > 0.001 + || p03.dot(p01).abs() > 0.001 + { + panic!("Infinite light portal isn't a planar quadrilateral (perpendicular edges)"); + } + + let portal_frame = Frame::from_xy(p03, p01); + + let width = src_res.x(); + let height = src_res.y(); + + let mut new_pixels = vec![0.0 as Float; (width * height * 3) as usize]; + + new_pixels + .par_chunks_mut((width * 3) as usize) + .enumerate() + .for_each(|(y, row_pixels)| { + let y = y as i32; + + for x in 0..width { + let uv = Point2f::new( + (x as Float + 0.5) / width as Float, + (y as Float + 0.5) / height as Float, + ); + + let (w_world, _) = Self::render_from_image(portal_frame, uv); + let w_local = render_from_light.apply_inverse_vector(w_world).normalize(); + let uv_equi = equal_area_sphere_to_square(w_local); + + let pixel_idx = (x * 3) as usize; + + for c in 0..3 { + let val = equal_area_image.bilerp_channel_with_wrap( + uv_equi, + c, + WrapMode::OctahedralSphere.into(), + ); + row_pixels[pixel_idx + c] = val; + } + } + }); + + let image = Image::new( + PixelFormat::F32, + src_res, + &["R", "G", "B"], + equal_area_image.encoding, + ); + + let duv_dw_closure = |p: Point2f| -> Float { + let (_, jacobian) = Self::render_from_image(portal_frame, p); + jacobian + }; + + let d = image.get_sampling_distribution( + duv_dw_closure, + Bounds2f::from_points(Point2f::new(0., 0.), Point2f::new(1., 1.)), + ); + + let distribution = WindowedPiecewiseConstant2D::new(d); + + Self { + base, + image, + image_color_space, + scale, + scene_center: Point3f::default(), + scene_radius: 0., + filename, + portal, + portal_frame, + distribution, + } + } + + pub fn image_lookup(&self, uv: Point2f, lambda: &SampledWavelengths) -> SampledSpectrum { + let mut rgb = RGB::default(); + for c in 0..3 { + rgb[c] = self.image.lookup_nearest_channel(uv, c) + } + let spec = + RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb)); + self.scale * spec.sample(lambda) + } + + pub fn image_from_render(&self, w_render: Vector3f) -> Option<(Point2f, Float)> { + let w = self.portal_frame.to_local(w_render); + + if w.z() <= 0.0 { + return None; + } + + let alpha = w.x().atan2(w.z()); + let beta = w.y().atan2(w.z()); + + let duv_dw = square(PI) * (1. - square(w.x())) * (1. - square(w.y())) / w.z(); + + Some(( + Point2f::new( + clamp_t((alpha + PI / 2.0) / PI, 0.0, 1.0), + clamp_t((beta + PI / 2.0) / PI, 0.0, 1.0), + ), + duv_dw, + )) + } + + pub fn image_bounds(&self, p: Point3f) -> Option { + let (p0, _) = self.image_from_render((self.portal[0] - p).normalize())?; + let (p1, _) = self.image_from_render((self.portal[2] - p).normalize())?; + + Some(Bounds2f::from_points(p0, p1)) + } + + pub fn area(&self) -> Float { + (self.portal[1] - self.portal[0]).norm() * (self.portal[3] - self.portal[0]).norm() + } + + pub fn render_from_image(portal_frame: Frame, uv: Point2f) -> (Vector3f, Float) { + let alpha = -PI / 2.0 + uv.x() * PI; + let beta = -PI / 2.0 + uv.y() * PI; + + let x = alpha.tan(); + let y = beta.tan(); + + let w = Vector3f::new(x, y, 1.0).normalize(); + + let duv_dw = square(PI) * (1.0 - square(w.x())) * (1.0 - square(w.y())) / w.z(); + + (portal_frame.from_local(w), duv_dw) + } +} + +impl LightTrait for InfinitePortalLight { + fn base(&self) -> &LightBase { + &self.base + } + fn sample_li( + &self, + ctx: &LightSampleContext, + u: Point2f, + lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + let b = self.image_bounds(ctx.p())?; + let (uv, map_pdf) = self.distribution.sample(u, b)?; + let (wi, duv_dw) = Self::render_from_image(self.portal_frame, uv); + if duv_dw == 0. { + return None; + } + let pdf = map_pdf / duv_dw; + 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); + Some(LightLiSample::new(l, wi, pdf, intr)) + } + + fn phi(&self, _lambda: SampledWavelengths) -> SampledSpectrum { + todo!() + } + + fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, _allow_incomplete_pdf: bool) -> Float { + let Some((uv, duv_dw)) = self.image_from_render(wi) else { + return 0.; + }; + let Some(b) = self.image_bounds(ctx.p()) else { + return 0.; + }; + let pdf = self.distribution.pdf(uv, b); + pdf / duv_dw + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + + fn le(&self, ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum { + let uv = self.image_from_render(ray.d.normalize()); + let b = self.image_bounds(ray.o); + match (uv, b) { + (Some((p, duv_dw)), Some(bounds)) if bounds.contains(p) => self.image_lookup(p, lambda), + _ => SampledSpectrum::new(0.0), + } + } + + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + todo!() + } + + fn bounds(&self) -> Option { + None + } +} diff --git a/src/lights/mod.rs b/src/lights/mod.rs index ce202ae..e47fc2f 100644 --- a/src/lights/mod.rs +++ b/src/lights/mod.rs @@ -1,27 +1,918 @@ pub mod diffuse; +pub mod infinite; +pub mod sampler; -#[derive(Debug, Clone)] -pub struct DiffuseAreaLight; -#[derive(Debug, Clone)] -pub struct DistantLight; -#[derive(Debug, Clone)] -pub struct GonioPhotometricLight; -#[derive(Debug, Clone)] -pub struct InfiniteAreaLight; -#[derive(Debug, Clone)] -pub struct PointLight; -#[derive(Debug, Clone)] -pub struct ProjectionLight; -#[derive(Debug, Clone)] -pub struct SpotLight; +use crate::core::interaction::{ + Interaction, InteractionTrait, MediumInteraction, SimpleInteraction, SurfaceInteraction, +}; +use crate::core::medium::MediumInterface; +use crate::core::pbrt::{Float, InternCache, PI}; +use crate::geometry::{ + Bounds2f, Bounds3f, DirectionCone, Normal3f, Point2f, Point2i, Point3f, Point3fi, Ray, + Vector3f, VectorLike, cos_theta, +}; +use crate::image::Image; +use crate::spectra::{ + DenselySampledSpectrum, LAMBDA_MAX, LAMBDA_MIN, RGBIlluminantSpectrum, SampledSpectrum, + SampledWavelengths, Spectrum, SpectrumTrait, +}; +use crate::utils::color::RGB; +use crate::utils::colorspace::RGBColorSpace; +use crate::utils::math::{equal_area_sphere_to_square, radians, safe_sqrt, smooth_step, square}; +use crate::utils::sampling::PiecewiseConstant2D; +use crate::utils::transform::Transform; +use bitflags::bitflags; -#[derive(Debug, Clone)] -pub enum Light { - DiffuseArea(Box), - Distant(Box), - GonioPhotometric(Box), - InfiniteArea(Box), - Point(Box), - Projection(Box), - Spot(Box), +use enum_dispatch::enum_dispatch; +use std::sync::{Arc, OnceLock}; + +use diffuse::DiffuseAreaLight; +use infinite::{InfiniteImageLight, InfinitePortalLight, InfiniteUniformLight}; + +static SPECTRUM_CACHE: OnceLock> = OnceLock::new(); + +fn get_spectrum_cache() -> &'static InternCache { + SPECTRUM_CACHE.get_or_init(InternCache::new) +} + +bitflags! { + #[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] + pub struct LightType: u32 { + const DeltaPosition = 1; + const DeltaDirection = 2; + const Area = 4; + const Infinite = 8; + } +} + +impl LightType { + pub fn is_infinite(&self) -> bool { + self.contains(LightType::Infinite) + } + + pub fn is_delta_light(&self) -> bool { + self.contains(LightType::DeltaPosition) || self.contains(LightType::DeltaDirection) + } +} + +#[derive(Debug, Clone)] +pub struct LightLeSample { + l: SampledSpectrum, + ray: Ray, + intr: Option, + pdf_pos: Float, + pdf_dir: Float, +} + +#[derive(Debug, Default, Clone)] +pub struct LightSampleContext { + pub pi: Point3fi, + pub n: Normal3f, + pub ns: Normal3f, +} + +impl LightSampleContext { + pub fn new(pi: Point3fi, n: Normal3f, ns: Normal3f) -> Self { + Self { pi, n, ns } + } + + pub fn p(&self) -> Point3f { + self.pi.into() + } +} + +impl From<&SurfaceInteraction> for LightSampleContext { + fn from(si: &SurfaceInteraction) -> Self { + Self { + pi: si.common.pi, + n: si.common.n, + ns: si.shading.n, + } + } +} + +impl From<&MediumInteraction> for LightSampleContext { + fn from(mi: &MediumInteraction) -> Self { + Self { + pi: mi.common.pi, + n: Normal3f::default(), + ns: Normal3f::default(), + } + } +} + +impl From<&Interaction> for LightSampleContext { + fn from(intr: &Interaction) -> Self { + match intr { + Interaction::Surface(si) => Self { + pi: si.common.pi, + n: si.common.n, + ns: si.shading.n, + }, + + Interaction::Medium(mi) => Self { + pi: mi.common.pi, + n: mi.common.n, + ns: mi.common.n, + }, + + Interaction::Simple(sim) => Self { + pi: sim.common.pi, + n: sim.common.n, + ns: sim.common.n, + }, + } + } +} + +#[derive(Debug, Clone)] +pub struct LightLiSample { + pub l: SampledSpectrum, + pub wi: Vector3f, + pub pdf: Float, + pub p_light: Arc, +} + +impl LightLiSample { + pub fn new(l: SampledSpectrum, wi: Vector3f, pdf: Float, p_light: Interaction) -> Self { + Self { + l, + wi, + pdf, + p_light: Arc::new(p_light), + } + } +} + +#[derive(Debug, Clone)] +pub struct LightBase { + pub light_type: LightType, + pub render_from_light: Transform, + pub medium_interface: MediumInterface, +} + +impl LightBase { + pub fn new( + light_type: LightType, + render_from_light: &Transform, + medium_interface: &MediumInterface, + ) -> Self { + Self { + light_type, + render_from_light: *render_from_light, + medium_interface: medium_interface.clone(), + } + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + SampledSpectrum::default() + } + + pub fn light_type(&self) -> LightType { + self.light_type + } + + pub fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + SampledSpectrum::default() + } + + pub fn lookup_spectrum(s: &Spectrum) -> Arc { + let cache = SPECTRUM_CACHE.get_or_init(InternCache::new); + let dense_spectrum = DenselySampledSpectrum::from_spectrum(s, LAMBDA_MIN, LAMBDA_MAX); + cache.lookup(dense_spectrum) + } +} + +#[derive(Debug, Clone)] +pub struct LightBounds { + bounds: Bounds3f, + phi: Float, + w: Vector3f, + cos_theta_o: Float, + cos_theta_e: Float, + two_sided: bool, +} + +impl LightBounds { + pub fn new( + bounds: &Bounds3f, + w: Vector3f, + phi: Float, + cos_theta_o: Float, + cos_theta_e: Float, + two_sided: bool, + ) -> Self { + Self { + bounds: *bounds, + phi, + w, + cos_theta_o, + cos_theta_e, + two_sided, + } + } + + pub fn centroid(&self) -> Point3f { + self.bounds.p_min + Vector3f::from(self.bounds.p_max) / 2. + } + + pub fn importance(&self, p: Point3f, n: Normal3f) -> Float { + // Compute clamped squared distance to reference point + let pc = self.centroid(); + let d2_raw = p.distance_squared(pc); + let d2 = d2_raw.max(self.bounds.diagonal().norm()) / 2.; + let cos_sub_clamped = |sin_theta_a: Float, + cos_theta_a: Float, + sin_theta_b: Float, + cos_theta_b: Float| + -> Float { + if cos_theta_a > cos_theta_b { + return 1.; + } + cos_theta_a * cos_theta_b + sin_theta_a * sin_theta_b + }; + + let sin_sub_clamped = |sin_theta_a: Float, + cos_theta_a: Float, + sin_theta_b: Float, + cos_theta_b: Float| + -> Float { + if cos_theta_a > cos_theta_b { + return 1.; + } + sin_theta_a * cos_theta_b - cos_theta_a * sin_theta_b + }; + + let wi = (p - pc).normalize(); + let mut cos_theta_w = self.w.dot(wi); + if self.two_sided { + cos_theta_w = cos_theta_w.abs(); + } + let sin_theta_w = safe_sqrt(1. - square(cos_theta_w)); + let cos_theta_b = DirectionCone::bound_subtended_directions(&self.bounds, p).cos_theta; + let sin_theta_b = safe_sqrt(1. - square(cos_theta_b)); + let sin_theta_o = safe_sqrt(1. - square(self.cos_theta_o)); + let cos_theta_x = cos_sub_clamped(sin_theta_w, cos_theta_w, sin_theta_o, self.cos_theta_o); + let sin_theta_x = sin_sub_clamped(sin_theta_w, cos_theta_w, sin_theta_o, self.cos_theta_o); + let cos_theta_p = cos_sub_clamped(sin_theta_x, cos_theta_x, sin_theta_b, cos_theta_b); + + if cos_theta_p <= self.cos_theta_e { + return 0.; + } + + let mut importance = self.phi * cos_theta_p / d2; + if n != Normal3f::new(0., 0., 0.) { + let cos_theta_i = wi.abs_dot(n.into()); + let sin_theta_i = safe_sqrt(1. - square(cos_theta_i)); + let cos_thetap_i = cos_sub_clamped(sin_theta_i, cos_theta_i, sin_theta_b, cos_theta_b); + importance *= cos_thetap_i; + } + importance + } + + pub fn union(a: &Self, b: &Self) -> Self { + if a.phi == 0. { + return a.clone(); + } + if b.phi == 0. { + return b.clone(); + } + + let a_cone = DirectionCone::new(a.w, a.cos_theta_o); + let b_cone = DirectionCone::new(b.w, b.cos_theta_o); + let cone = DirectionCone::union(&a_cone, &b_cone); + let cos_theta_o = cone.cos_theta; + let cos_theta_e = a.cos_theta_e.min(b.cos_theta_e); + LightBounds::new( + &a.bounds.union(b.bounds), + cone.w, + a.phi + b.phi, + cos_theta_o, + cos_theta_e, + a.two_sided || b.two_sided, + ) + } +} + +#[enum_dispatch] +pub trait LightTrait: Send + Sync + std::fmt::Debug { + fn base(&self) -> &LightBase; + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum; + fn sample_li( + &self, + ctx: &LightSampleContext, + u: Point2f, + lambda: &SampledWavelengths, + allow_incomplete_pdf: bool, + ) -> Option; + fn pdf_li(&self, ctx: &LightSampleContext, wi: Vector3f, allow_incomplete_pdf: bool) -> Float; + fn l( + &self, + p: Point3f, + n: Normal3f, + uv: Point2f, + w: Vector3f, + lambda: &SampledWavelengths, + ) -> SampledSpectrum; + fn le(&self, ray: &Ray, lambda: &SampledWavelengths) -> SampledSpectrum; + fn preprocess(&mut self, scene_bounds: &Bounds3f); + fn bounds(&self) -> Option; + fn light_type(&self) -> LightType { + self.base().light_type() + } +} + +#[derive(Debug, Clone)] +#[enum_dispatch(LightTrait)] +#[allow(clippy::large_enum_variant)] +pub enum Light { + DiffuseArea(DiffuseAreaLight), + Distant(DistantLight), + Goniometric(GoniometricLight), + InfiniteUniform(InfiniteUniformLight), + InfiniteImage(InfiniteImageLight), + InfinitePortal(InfinitePortalLight), + Point(PointLight), + Projection(ProjectionLight), + Spot(SpotLight), +} + +#[derive(Debug, Clone)] +pub struct DistantLight { + pub base: LightBase, + lemit: Arc, + scale: Float, + scene_center: Point3f, + scene_radius: Float, +} + +impl DistantLight { + pub fn new(render_from_light: &Transform, lemit: Spectrum, scale: Float) -> Self { + let base = LightBase::new( + LightType::DeltaDirection, + render_from_light, + &MediumInterface::default(), + ); + + let l_interned = LightBase::lookup_spectrum(&lemit); + Self { + base, + lemit: l_interned, + scale, + scene_center: Point3f::default(), + scene_radius: 0., + } + } +} + +impl LightTrait for DistantLight { + fn base(&self) -> &LightBase { + &self.base + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + self.scale * self.lemit.sample(&lambda) * PI * self.scene_radius.sqrt() + } + + fn sample_li( + &self, + ctx: &LightSampleContext, + _u: Point2f, + lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + let wi = self + .base + .render_from_light + .apply_to_vector(Vector3f::new(0., 0., 1.)) + .normalize(); + let p_outside = ctx.p() + wi * 2. * self.scene_radius; + + let li = self.scale * self.lemit.sample(&lambda); + let intr = SimpleInteraction::new( + Point3fi::new_from_point(p_outside), + 0.0, + Some(self.base.medium_interface.clone()), + ); + + Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr))) + } + + fn pdf_li( + &self, + _ctx: &LightSampleContext, + _wi: Vector3f, + _allow_incomplete_pdf: bool, + ) -> Float { + 0. + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + todo!() + } + + fn preprocess(&mut self, scene_bounds: &Bounds3f) { + let (center, radius) = scene_bounds.bounding_sphere(); + self.scene_center = center; + self.scene_radius = radius; + } + + fn bounds(&self) -> Option { + None + } +} + +#[derive(Debug, Clone)] +pub struct GoniometricLight { + pub base: LightBase, + iemit: Arc, + scale: Float, + image: Image, + distrib: PiecewiseConstant2D, +} + +impl GoniometricLight { + pub fn new( + render_from_light: &Transform, + medium_interface: &MediumInterface, + iemit: Spectrum, + scale: Float, + image: Image, + ) -> Self { + let base = LightBase::new( + LightType::DeltaPosition, + render_from_light, + medium_interface, + ); + + let i_interned = LightBase::lookup_spectrum(&iemit); + let d = image.get_sampling_distribution_uniform(); + let distrib = PiecewiseConstant2D::new_with_data(&d); + Self { + base, + iemit: i_interned, + scale, + image, + distrib, + } + } + + 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) + } +} + +impl LightTrait for GoniometricLight { + fn base(&self) -> &LightBase { + &self.base + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + let mut sum_y = 0.; + for y in 0..self.image.resolution.y() { + for x in 0..self.image.resolution.x() { + sum_y += self.image.get_channel(Point2i::new(x, y), 0); + } + } + self.scale * self.iemit.sample(&lambda) * 4. * PI * sum_y + / (self.image.resolution.x() * self.image.resolution.y()) as Float + } + + fn sample_li( + &self, + _ctx: &LightSampleContext, + _u: Point2f, + _lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + todo!() + } + + fn pdf_li( + &self, + _ctx: &LightSampleContext, + _wi: Vector3f, + _allow_incomplete_pdf: bool, + ) -> Float { + 0. + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + todo!() + } + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + todo!() + } + fn bounds(&self) -> Option { + todo!() + } +} + +#[derive(Debug, Clone)] +pub struct PointLight { + base: LightBase, + i: Arc, + scale: Float, +} + +impl PointLight { + pub fn new( + render_from_light: Transform, + medium_interface: MediumInterface, + i: &Spectrum, + scale: Float, + ) -> Self { + let base = LightBase::new( + LightType::DeltaPosition, + &render_from_light, + &medium_interface, + ); + + let i_interned = LightBase::lookup_spectrum(i); + + Self { + base, + i: i_interned, + scale, + } + } +} + +impl LightTrait for PointLight { + fn base(&self) -> &LightBase { + &self.base + } + + fn sample_li( + &self, + ctx: &LightSampleContext, + _u: Point2f, + lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + let pi = self + .base + .render_from_light + .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 intr = SimpleInteraction::new(pi, 0.0, Some(self.base.medium_interface.clone())); + Some(LightLiSample::new(li, wi, 1., Interaction::Simple(intr))) + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + 4. * PI * self.scale * self.i.sample(&lambda) + } + + fn pdf_li( + &self, + _ctx: &LightSampleContext, + _wi: Vector3f, + _allow_incomplete_pdf: bool, + ) -> Float { + 0. + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + todo!() + } + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + todo!() + } + fn bounds(&self) -> Option { + let p = self + .base + .render_from_light + .apply_to_point(Point3f::new(0., 0., 0.)); + let phi = 4. * PI * self.scale * self.i.max_value(); + Some(LightBounds::new( + &Bounds3f::from_points(p, p), + Vector3f::new(0., 0., 1.), + phi, + PI.cos(), + (PI / 2.).cos(), + false, + )) + } +} + +#[derive(Debug, Clone)] +pub struct ProjectionLight { + base: LightBase, + image: Image, + image_color_space: Arc, + scale: Float, + screen_bounds: Bounds2f, + hither: Float, + screen_from_light: Transform, + light_from_screen: Transform, + a: Float, + distrib: PiecewiseConstant2D, +} + +impl ProjectionLight { + pub fn new( + render_from_light: Transform, + medium_interface: MediumInterface, + image: Image, + image_color_space: Arc, + scale: Float, + fov: Float, + ) -> Self { + let base = LightBase::new( + LightType::DeltaPosition, + &render_from_light, + &medium_interface, + ); + let aspect = image.resolution().x() as Float / image.resolution().y() as Float; + let screen_bounds = if aspect > 1. { + Bounds2f::from_points(Point2f::new(-aspect, -1.), Point2f::new(aspect, 1.)) + } else { + Bounds2f::from_points( + Point2f::new(-1., 1. / aspect), + Point2f::new(1., 1. / aspect), + ) + }; + + let hither = 1e-3; + let screen_from_light = Transform::perspective(fov, hither, 1e30).unwrap(); + let light_from_screen = screen_from_light.inverse(); + let opposite = (radians(fov) / 2.).tan(); + let aspect_ratio = if aspect > 1. { aspect } else { 1. / aspect }; + let a = 4. * square(opposite) * aspect_ratio; + let dwda = |p: Point2f| { + let w = + Vector3f::from(light_from_screen.apply_to_point(Point3f::new(p.x(), p.y(), 0.))); + cos_theta(w.normalize()).powi(3) + }; + + let d = image.get_sampling_distribution(dwda, screen_bounds); + let distrib = PiecewiseConstant2D::new_with_bounds(&d, screen_bounds); + + Self { + base, + image, + image_color_space, + screen_bounds, + screen_from_light, + light_from_screen, + scale, + hither, + a, + distrib, + } + } + + pub fn i(&self, w: Vector3f, lambda: SampledWavelengths) -> SampledSpectrum { + if w.z() < self.hither { + return SampledSpectrum::new(0.); + } + let ps = self.screen_from_light.apply_to_point(w.into()); + if !self.screen_bounds.contains(Point2f::new(ps.x(), ps.y())) { + return SampledSpectrum::new(0.); + } + let uv = Point2f::from(self.screen_bounds.offset(&Point2f::new(ps.x(), ps.y()))); + let mut rgb = RGB::default(); + for c in 0..3 { + rgb[c] = self.image.lookup_nearest_channel(uv, c); + } + let s = RGBIlluminantSpectrum::new(self.image_color_space.as_ref(), RGB::clamp_zero(rgb)); + self.scale * s.sample(&lambda) + } +} + +impl LightTrait for ProjectionLight { + fn base(&self) -> &LightBase { + &self.base + } + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + let mut sum = SampledSpectrum::new(0.); + for y in 0..self.image.resolution.y() { + for x in 0..self.image.resolution.x() { + let ps = self.screen_bounds.lerp(Point2f::new( + (x as Float + 0.5) / self.image.resolution.x() as Float, + (y as Float + 0.5) / self.image.resolution.y() as Float, + )); + let w_raw = Vector3f::from(self.light_from_screen.apply_to_point(Point3f::new( + ps.x(), + ps.y(), + 0., + ))); + let w = w_raw.normalize(); + 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(); + } + + let s = RGBIlluminantSpectrum::new( + self.image_color_space.as_ref(), + RGB::clamp_zero(rgb), + ); + sum += s.sample(&lambda) * dwda; + } + } + self.scale * self.a * sum / (self.image.resolution.x() * self.image.resolution.y()) as Float + } + + fn sample_li( + &self, + _ctx: &LightSampleContext, + _u: Point2f, + _lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + todo!() + } + fn pdf_li( + &self, + _ctx: &LightSampleContext, + _wi: Vector3f, + _allow_incomplete_pdf: bool, + ) -> Float { + todo!() + } + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + todo!() + } + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + todo!() + } + fn bounds(&self) -> Option { + todo!() + } +} + +#[derive(Debug, Clone)] +pub struct SpotLight { + base: LightBase, + iemit: Arc, + scale: Float, + cos_fallof_start: Float, + cos_fallof_end: Float, +} + +impl SpotLight { + pub fn new( + render_from_light: &Transform, + medium_interface: &MediumInterface, + iemit: Spectrum, + scale: Float, + total_width: Float, + fallof_start: Float, + ) -> Self { + let base = LightBase::new( + LightType::DeltaPosition, + render_from_light, + medium_interface, + ); + + let i_interned = LightBase::lookup_spectrum(&iemit); + let cos_fallof_end = radians(total_width).cos(); + let cos_fallof_start = radians(fallof_start).cos(); + assert!(fallof_start < total_width); + Self { + base, + iemit: i_interned, + scale, + cos_fallof_start, + cos_fallof_end, + } + } + + pub fn i(&self, w: Vector3f, lambda: SampledWavelengths) -> SampledSpectrum { + smooth_step(cos_theta(w), self.cos_fallof_end, self.cos_fallof_start) + * self.scale + * self.iemit.sample(&lambda) + } +} + +impl LightTrait for SpotLight { + fn base(&self) -> &LightBase { + &self.base + } + fn sample_li( + &self, + ctx: &LightSampleContext, + _u: Point2f, + lambda: &SampledWavelengths, + _allow_incomplete_pdf: bool, + ) -> Option { + let pi = self + .base + .render_from_light + .apply_to_interval(&Point3fi::default()); + let p: Point3f = pi.into(); + let wi = (p - ctx.p()).normalize(); + let w_light = self.base.render_from_light.apply_inverse_vector(-wi); + let li = self.i(w_light, *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))) + } + + fn phi(&self, lambda: SampledWavelengths) -> SampledSpectrum { + self.scale + * self.iemit.sample(&lambda) + * 2. + * PI + * ((1. - self.cos_fallof_start) + (self.cos_fallof_start - self.cos_fallof_end) / 2.) + } + + fn pdf_li( + &self, + _ctx: &LightSampleContext, + _wi: Vector3f, + _allow_incomplete_pdf: bool, + ) -> Float { + 0. + } + + fn l( + &self, + _p: Point3f, + _n: Normal3f, + _uv: Point2f, + _w: Vector3f, + _lambda: &SampledWavelengths, + ) -> SampledSpectrum { + todo!() + } + fn le(&self, _ray: &Ray, _lambda: &SampledWavelengths) -> SampledSpectrum { + todo!() + } + fn preprocess(&mut self, _scene_bounds: &Bounds3f) { + todo!() + } + + fn bounds(&self) -> Option { + let p = self + .base + .render_from_light + .apply_to_point(Point3f::default()); + let w = self + .base + .render_from_light + .apply_to_vector(Vector3f::new(0., 0., 1.)) + .normalize(); + let phi = self.scale * self.iemit.max_value() * 4. * PI; + let cos_theta_e = (self.cos_fallof_end.acos() - self.cos_fallof_start.acos()).cos(); + Some(LightBounds::new( + &Bounds3f::from_points(p, p), + w, + phi, + self.cos_fallof_start, + cos_theta_e, + false, + )) + } } diff --git a/src/lights/sampler.rs b/src/lights/sampler.rs new file mode 100644 index 0000000..c8a9ef0 --- /dev/null +++ b/src/lights/sampler.rs @@ -0,0 +1,521 @@ +use super::{ + Bounds3f, Float, Light, LightBounds, LightSampleContext, LightTrait, Normal3f, PI, Point3f, + SampledSpectrum, SampledWavelengths, Vector3f, VectorLike, safe_sqrt, square, +}; +use crate::geometry::primitives::OctahedralVector; +use crate::geometry::{DirectionCone, Normal}; +use crate::utils::math::sample_discrete; +use std::collections::HashMap; +use std::sync::Arc; + +use crate::core::pbrt::{ONE_MINUS_EPSILON, clamp_t, lerp}; +use crate::utils::sampling::AliasTable; +use enum_dispatch::enum_dispatch; + +#[derive(Clone, Copy, Debug, Default)] +#[repr(C)] +pub struct CompactLightBounds { + pub w: OctahedralVector, + pub phi: Float, + + // [0..15] = qCosTheta_o + // [15..30] = qCosTheta_e + // [30..31] = twoSided + // [31..32] = Unused/Padding + packed_info: u32, + + pub qb: [[u16; 3]; 2], +} + +const _: () = assert!(std::mem::size_of::() == 24); + +impl CompactLightBounds { + pub fn new(lb: &LightBounds, all_b: &Bounds3f) -> Self { + let q_cos_o = Self::quantize_cos(lb.cos_theta_o); + let q_cos_e = Self::quantize_cos(lb.cos_theta_e); + let two_sided = if lb.two_sided { 1 } else { 0 }; + + // | - twoSided (1) - | - qCosTheta_e (15) - | - qCosTheta_o (15) - | + let packed_info = (q_cos_o & 0x7FFF) | ((q_cos_e & 0x7FFF) << 15) | (two_sided << 30); + + let mut qb = [[0u16; 3]; 2]; + for i in 0..3 { + qb[0][i] = Self::quantize_bounds(lb.bounds.p_min[i], all_b.p_min[i], all_b.p_max[i]) + .floor() as u16; + qb[1][i] = Self::quantize_bounds(lb.bounds.p_max[i], all_b.p_min[i], all_b.p_max[i]) + .ceil() as u16; + } + + Self { + w: OctahedralVector::new(lb.w.normalize()), + phi: lb.phi, + packed_info, + qb, + } + } + + #[inline(always)] + pub fn two_sided(&self) -> bool { + (self.packed_info >> 30) & 1 == 1 + } + + #[inline(always)] + fn q_cos_theta_o(&self) -> u32 { + self.packed_info & 0x7FFF + } + + #[inline(always)] + fn q_cos_theta_e(&self) -> u32 { + (self.packed_info >> 15) & 0x7FFF + } + + #[inline] + pub fn cos_theta_o(&self) -> Float { + 2.0 * (self.q_cos_theta_o() as Float / 32767.0) - 1.0 + } + + #[inline] + pub fn cos_theta_e(&self) -> Float { + 2.0 * (self.q_cos_theta_e() as Float / 32767.0) - 1.0 + } + + pub fn importance(&self, p: Point3f, n: Normal3f, all_b: &Bounds3f) -> Float { + let bounds = self.bounds(all_b); + let cos_o = self.cos_theta_o(); + let cos_e = self.cos_theta_e(); + + let pc = bounds.centroid(); + let d2 = p.distance_squared(pc).max(bounds.diagonal().norm() * 0.5); + + let cos_sub_clamped = |sin_a: Float, cos_a: Float, sin_b: Float, cos_b: Float| { + if cos_a > cos_b { + 1.0 + } else { + cos_a * cos_b + sin_a * sin_b + } + }; + let sin_sub_clamped = |sin_a: Float, cos_a: Float, sin_b: Float, cos_b: Float| { + if cos_a > cos_b { + 0.0 + } else { + sin_a * cos_b - cos_a * sin_b + } + }; + + let wi = (p - pc).normalize(); + let w_vec = self.w.to_vector(); + let mut cos_w = w_vec.dot(wi); + if self.two_sided() { + cos_w = cos_w.abs(); + } + let sin_w = safe_sqrt(1.0 - square(cos_w)); + let cos_b = DirectionCone::bound_subtended_directions(&bounds, p).cos_theta; + let sin_b = safe_sqrt(1. - square(cos_b)); + let sin_o = safe_sqrt(1. - square(cos_o)); + let cos_x = cos_sub_clamped(sin_w, cos_w, sin_o, cos_o); + let sin_x = sin_sub_clamped(sin_w, cos_w, sin_o, cos_o); + let cos_p = cos_sub_clamped(sin_x, cos_x, sin_b, cos_b); + if cos_p <= cos_e { + return 0.; + } + let mut importance = self.phi * cos_p / d2; + if n != Normal3f::zero() { + let cos_i = wi.abs_dot(n.into()); + let sin_i = safe_sqrt(1. - square(cos_i)); + let cos_pi = cos_sub_clamped(sin_i, cos_i, sin_b, cos_b); + importance *= cos_pi; + } + importance + } + + pub fn bounds(&self, all_b: &Bounds3f) -> Bounds3f { + let mut p_min = Point3f::default(); + let mut p_max = Point3f::default(); + + for i in 0..3 { + let t_min = self.qb[0][i] as Float / 65535.0; + let t_max = self.qb[1][i] as Float / 65535.0; + p_min[i] = lerp(t_min, all_b.p_min[i], all_b.p_max[i]); + p_max[i] = lerp(t_max, all_b.p_min[i], all_b.p_max[i]); + } + Bounds3f::from_points(p_min, p_max) + } + + fn quantize_cos(c: Float) -> u32 { + (32767.0 * ((c + 1.0) * 0.5)).floor() as u32 + } + + fn quantize_bounds(c: Float, min: Float, max: Float) -> Float { + if min == max { + return 0.0; + } + 65535.0 * clamp_t((c - min) / (max - min), 0.0, 1.0) + } +} + +#[derive(Debug, Clone)] +pub struct SampledLight { + pub light: Arc, + pub p: Float, +} + +impl SampledLight { + pub fn new(light: Arc, p: Float) -> Self { + Self { light, p } + } +} + +#[enum_dispatch] +pub trait LightSamplerTrait: Send + Sync + std::fmt::Debug { + fn sample_with_context(&self, ctx: &LightSampleContext, u: Float) -> Option; + fn pmf_with_context(&self, ctx: &LightSampleContext, light: &Arc) -> Float; + fn sample(&self, u: Float) -> Option; + fn pmf(&self, light: &Arc) -> Float; +} + +#[derive(Clone, Debug)] +#[enum_dispatch(LightSamplerTrait)] +pub enum LightSampler { + Uniform(UniformLightSampler), + Power(PowerLightSampler), + BVH(BVHLightSampler), +} + +#[derive(Clone, Debug)] +pub struct UniformLightSampler { + lights: Vec>, +} + +impl UniformLightSampler { + pub fn new(lights: &[Arc]) -> Self { + Self { + lights: lights.to_vec(), + } + } +} + +impl LightSamplerTrait for UniformLightSampler { + fn sample_with_context(&self, _ctx: &LightSampleContext, u: Float) -> Option { + self.sample(u) + } + fn pmf_with_context(&self, _ctx: &LightSampleContext, light: &Arc) -> Float { + self.pmf(light) + } + fn sample(&self, u: Float) -> Option { + if self.lights.is_empty() { + return None; + } + + let light_index = (u as usize * self.lights.len()).min(self.lights.len() - 1); + Some(SampledLight { + light: self.lights[light_index].clone(), + p: 1. / self.lights.len() as Float, + }) + } + fn pmf(&self, _light: &Arc) -> Float { + if self.lights.is_empty() { + return 0.; + } + 1. / self.lights.len() as Float + } +} + +#[derive(Clone, Debug)] +pub struct PowerLightSampler { + lights: Vec>, + light_to_index: HashMap, + alias_table: AliasTable, +} + +impl PowerLightSampler { + pub fn new(lights: &[Arc]) -> Self { + if lights.is_empty() { + return Self { + lights: Vec::new(), + light_to_index: HashMap::new(), + alias_table: AliasTable::new(&[]), + }; + } + + let mut lights_vec = Vec::with_capacity(lights.len()); + let mut light_to_index = HashMap::with_capacity(lights.len()); + let mut light_power = Vec::with_capacity(lights.len()); + + let lambda = SampledWavelengths::sample_visible(0.5); + + for (i, light) in lights.iter().enumerate() { + lights_vec.push(light.clone()); + + let ptr = Arc::as_ptr(light) as usize; + light_to_index.insert(ptr, i); + + let phi = SampledSpectrum::safe_div(&light.phi(lambda), &lambda.pdf()); + light_power.push(phi.average()); + } + + let alias_table = AliasTable::new(&light_power); + + Self { + lights: lights_vec, + light_to_index, + alias_table, + } + } +} + +impl LightSamplerTrait for PowerLightSampler { + fn sample_with_context(&self, _ctx: &LightSampleContext, u: Float) -> Option { + self.sample(u) + } + + fn pmf_with_context(&self, _ctx: &LightSampleContext, light: &Arc) -> Float { + self.pmf(light) + } + + fn sample(&self, u: Float) -> Option { + if self.alias_table.size() == 0 { + return None; + } + + let (light_index, pmf, _) = self.alias_table.sample(u); + + Some(SampledLight { + light: self.lights[light_index].clone(), + p: pmf, + }) + } + + fn pmf(&self, light: &Arc) -> Float { + if self.alias_table.size() == 0 { + return 0.; + } + + let ptr = Arc::as_ptr(light) as usize; + + if let Some(&index) = self.light_to_index.get(&ptr) { + self.alias_table.pmf(index) + } else { + 0.0 + } + } +} + +#[derive(Clone, Copy, Debug, Default)] +#[repr(C, align(32))] +pub struct LightBVHNode { + pub light_bounds: CompactLightBounds, + + // Bit 31 (MSB) : isLeaf (1 bit) + // Bits 0..31 : childOrLightIndex (31 bits) + packed_data: u32, +} + +const _: () = assert!(std::mem::size_of::() == 32); + +impl LightBVHNode { + /// Mask to isolate the Leaf Flag (Bit 31) + const LEAF_MASK: u32 = 0x8000_0000; + /// Mask to isolate the Index (Bits 0-30) + const INDEX_MASK: u32 = 0x7FFF_FFFF; + + pub fn make_leaf(light_index: u32, cb: CompactLightBounds) -> Self { + debug_assert!( + (light_index & Self::LEAF_MASK) == 0, + "Light index too large" + ); + + Self { + light_bounds: cb, + // Set index and flip the MSB to 1 + packed_data: light_index | Self::LEAF_MASK, + } + } + + pub fn make_interior(child_index: u32, cb: CompactLightBounds) -> Self { + debug_assert!( + (child_index & Self::LEAF_MASK) == 0, + "Child index too large" + ); + + Self { + light_bounds: cb, + // Set index, MSB remains 0 + packed_data: child_index, + } + } + + #[inline(always)] + pub fn is_leaf(&self) -> bool { + (self.packed_data & Self::LEAF_MASK) != 0 + } + + #[inline(always)] + pub fn light_index(&self) -> u32 { + debug_assert!(self.is_leaf()); + self.packed_data & Self::INDEX_MASK + } + + #[inline(always)] + pub fn child_index(&self) -> u32 { + debug_assert!(!self.is_leaf()); + self.packed_data & Self::INDEX_MASK + } + + #[inline(always)] + pub fn child_or_light_index(&self) -> u32 { + self.packed_data & Self::INDEX_MASK + } + + pub fn sample(&self, _ctx: &LightSampleContext, _u: Float) -> Option { + todo!("Implement LightBVHNode::Sample logic") + } +} + +#[derive(Clone, Debug)] +pub struct BVHLightSampler { + lights: Vec>, + infinite_lights: Vec>, + all_light_bounds: Bounds3f, + nodes: Vec, + light_to_bit_trail: HashMap, +} + +impl BVHLightSampler { + fn evaluate_cost(&self, b: &LightBounds, bounds: &Bounds3f, dim: usize) -> Float { + let theta_o = b.cos_theta_o.acos(); + let theta_e = b.cos_theta_e.acos(); + let theta_w = (theta_o + theta_e).min(PI); + let sin_o = safe_sqrt(1. - square(b.cos_theta_o)); + let m_omega = 2. * PI * (1. - b.cos_theta_o) + + PI / 2. + * (2. * theta_w - (theta_o - 2. * theta_w).cos() - 2. * theta_o * sin_o + + b.cos_theta_o); + let kr = bounds.diagonal().max_component_value() / bounds.diagonal()[dim]; + b.phi * m_omega * kr * b.bounds.surface_area() + } +} + +impl LightSamplerTrait for BVHLightSampler { + fn sample_with_context(&self, ctx: &LightSampleContext, mut u: Float) -> Option { + let empty_nodes = if self.nodes.is_empty() { 0. } else { 1. }; + let inf_size = self.infinite_lights.len() as Float; + let light_size = self.lights.len() as Float; + + let p_inf = inf_size / (inf_size + empty_nodes); + + if u < p_inf { + u /= p_inf; + let ind = (u * light_size).min(light_size - 1.) as usize; + let pmf = p_inf / inf_size; + Some(SampledLight::new(self.infinite_lights[ind].clone(), pmf)) + } else { + if self.nodes.is_empty() { + return None; + } + let p = ctx.p(); + let n = ctx.ns; + u = ((u - p_inf) / (1. - p_inf)).min(ONE_MINUS_EPSILON); + let mut node_ind = 0; + let mut pmf = 1. - p_inf; + + loop { + let node = self.nodes[node_ind]; + if !node.is_leaf() { + let children: [LightBVHNode; 2] = [ + self.nodes[node_ind + 1], + self.nodes[node.child_or_light_index() as usize], + ]; + let ci: [Float; 2] = [ + children[0] + .light_bounds + .importance(p, n, &self.all_light_bounds), + children[1] + .light_bounds + .importance(p, n, &self.all_light_bounds), + ]; + + if ci[0] == 0. && ci[1] == 0. { + return None; + } + + let mut node_pmf: Float = 0.; + let child = sample_discrete(&ci, u, Some(&mut node_pmf), Some(&mut u))?; + pmf *= node_pmf; + node_ind = if child == 0 { + node_ind + 1 + } else { + node.child_or_light_index() as usize + }; + } else { + if node_ind > 0 + || node.light_bounds.importance(p, n, &self.all_light_bounds) > 0. + { + return Some(SampledLight::new( + self.lights[node.child_or_light_index() as usize].clone(), + pmf, + )); + } + return None; + } + } + } + } + + fn pmf_with_context(&self, ctx: &LightSampleContext, light: &Arc) -> Float { + let ptr = Arc::as_ptr(light) as usize; + let empty_nodes = if self.nodes.is_empty() { 0. } else { 1. }; + if self.light_to_bit_trail.contains_key(&ptr) { + return 1. / (self.infinite_lights.len() as Float + empty_nodes); + } + + let mut bit_trail = self.light_to_bit_trail[&ptr]; + let p = ctx.p(); + let n = ctx.ns; + let p_inf = self.infinite_lights.len() as Float + / (self.infinite_lights.len() as Float + empty_nodes); + let mut pmf = 1. - p_inf; + let mut node_ind = 0; + + loop { + let node = self.nodes[node_ind]; + if node.is_leaf() { + return pmf; + } + let child0 = self.nodes[node_ind + 1]; + let child1 = self.nodes[node.child_or_light_index() as usize]; + let ci = [ + child0.light_bounds.importance(p, n, &self.all_light_bounds), + child1.light_bounds.importance(p, n, &self.all_light_bounds), + ]; + pmf *= ci[bit_trail & 1] / (ci[0] + ci[1]); + node_ind = if (bit_trail & 1) != 0 { + node.child_or_light_index() as usize + } else { + node_ind + 1 + }; + bit_trail >>= 1; + } + } + + fn sample(&self, u: Float) -> Option { + if self.lights.is_empty() { + return None; + } + + let light_ind = + (u * self.lights.len() as Float).min(self.lights.len() as Float - 1.) as usize; + + Some(SampledLight::new( + self.lights[light_ind].clone(), + 1. / self.lights.len() as Float, + )) + } + + fn pmf(&self, _light: &Arc) -> Float { + if self.lights.is_empty() { + return 0.; + } + + 1. / self.lights.len() as Float + } +}