From 11a731247dd7cbf4cfaf7be7186d070e3953fd91 Mon Sep 17 00:00:00 2001 From: pingu Date: Tue, 9 Dec 2025 14:07:46 +0000 Subject: [PATCH] Updated light management and creation, refactored module. Started work on rendering and integrators, fixed up maths functions --- Cargo.toml | 1 + src/camera/mod.rs | 86 +++-- src/camera/orthographic.rs | 8 +- src/camera/perspective.rs | 8 +- src/camera/realistic.rs | 7 +- src/camera/spherical.rs | 9 +- src/core/aggregates.rs | 12 +- src/core/bssrdf.rs | 11 +- src/core/bxdf.rs | 8 +- src/core/film.rs | 127 ++++--- src/core/filter.rs | 18 +- src/core/integrator.rs | 13 - src/core/interaction.rs | 334 +++++++++------- src/core/material.rs | 393 ++++++++++++++----- src/core/medium.rs | 753 ++++++++++++++++++++++++++++++++++++- src/core/mod.rs | 2 +- src/core/options.rs | 17 +- src/core/pbrt.rs | 30 ++ src/core/primitive.rs | 75 ++-- src/core/sampler.rs | 3 - src/core/texture.rs | 44 ++- src/geometry/bounds.rs | 19 +- src/geometry/cone.rs | 110 +++--- src/geometry/primitives.rs | 109 +++++- src/geometry/ray.rs | 2 +- src/image/io.rs | 27 +- src/image/metadata.rs | 62 ++- src/image/mod.rs | 277 +++++++++++++- src/image/ops.rs | 14 +- src/lib.rs | 1 + src/shapes/bilinear.rs | 34 +- src/shapes/curves.rs | 34 +- src/shapes/cylinder.rs | 23 +- src/shapes/disk.rs | 23 +- src/shapes/mod.rs | 180 +++------ src/shapes/sphere.rs | 25 +- src/shapes/triangle.rs | 30 +- src/spectra/mod.rs | 12 +- src/spectra/sampled.rs | 28 +- src/spectra/simple.rs | 42 +++ src/utils/color.rs | 2 +- src/utils/containers.rs | 215 ++++++++++- src/utils/hash.rs | 4 +- src/utils/math.rs | 58 ++- src/utils/mipmap.rs | 25 +- src/utils/mod.rs | 48 +++ src/utils/rng.rs | 11 +- src/utils/sampling.rs | 327 +++++++++++++++- src/utils/transform.rs | 130 +++++-- 49 files changed, 3008 insertions(+), 823 deletions(-) delete mode 100644 src/core/integrator.rs diff --git a/Cargo.toml b/Cargo.toml index c6de8ac..9d4483f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,6 +11,7 @@ enum_dispatch = "0.3.13" exr = "1.73.0" half = "2.7.1" image_rs = { package = "image", version = "0.25.8" } +indicatif = "0.18.3" num = "0.4.3" num-integer = "0.1.46" num-traits = "0.2.19" diff --git a/src/camera/mod.rs b/src/camera/mod.rs index 533e653..7e2293d 100644 --- a/src/camera/mod.rs +++ b/src/camera/mod.rs @@ -8,14 +8,20 @@ pub use perspective::PerspectiveCamera; pub use realistic::RealisticCamera; pub use spherical::SphericalCamera; -use crate::core::film::Film; +use crate::core::film::{Film, FilmTrait}; +use crate::core::interaction::Interaction; use crate::core::medium::Medium; use crate::core::pbrt::{Float, RenderingCoordinateSystem, lerp}; use crate::core::sampler::CameraSample; -use crate::geometry::{Normal3f, Point3f, Ray, RayDifferential, Vector3f, VectorLike}; +use crate::geometry::{ + Normal3f, Point2f, Point2i, Point3f, Ray, RayDifferential, Vector3f, VectorLike, +}; +use crate::image::ImageMetadata; use crate::spectra::{SampledSpectrum, SampledWavelengths}; use crate::utils::transform::{AnimatedTransform, Transform}; +use enum_dispatch::enum_dispatch; + use std::sync::Arc; #[derive(Debug, Clone)] @@ -24,6 +30,17 @@ pub struct CameraRay { pub weight: SampledSpectrum, } +#[derive(Debug, Clone)] +pub struct CameraWiSample { + wi_spec: SampledSpectrum, + wi: Vector3f, + pdf: Float, + p_raster: Point2f, + p_ref: Interaction, + p_lens: Interaction, +} + +#[derive(Debug)] pub struct CameraTransform { render_from_camera: AnimatedTransform, world_from_render: Transform, @@ -63,6 +80,10 @@ impl CameraTransform { } } + pub fn camera_from_world(&self, time: Float) -> Transform { + (self.world_from_render * self.render_from_camera.interpolate(time)).inverse() + } + pub fn render_from_camera(&self, p: Point3f, time: Float) -> Point3f { self.render_from_camera.apply_point(p, time) } @@ -88,6 +109,7 @@ impl CameraTransform { } } +#[derive(Debug)] pub struct CameraBase { pub camera_transform: CameraTransform, pub shutter_open: Float, @@ -100,8 +122,28 @@ pub struct CameraBase { pub min_dir_differential_y: Vector3f, } +impl CameraBase { + pub fn init_metadata(&self, metadata: &mut ImageMetadata) { + let camera_from_world: Transform = + self.camera_transform.camera_from_world(self.shutter_open); + + metadata.camera_from_world = Some(camera_from_world.get_matrix()); + } +} + +#[enum_dispatch] pub trait CameraTrait { fn base(&self) -> &CameraBase; + fn get_film(&self) -> Film { + self.base().film.clone() + } + + fn resolution(&self) -> Point2i { + self.base().film.full_resolution() + } + + fn init_metadata(&self, metadata: &mut ImageMetadata); + fn generate_ray(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option; fn generate_ray_differential( &self, @@ -187,8 +229,8 @@ pub trait CameraTrait { let n_down = Vector3f::from(n_down_z); let tx = -(n_down.dot(y_ray.o.into())) / n_down.dot(x_ray.d); let ty = -(n_down.dot(x_ray.o.into()) - d) / n_down.dot(y_ray.d); - let px = x_ray.evaluate(tx); - let py = y_ray.evaluate(ty); + let px = x_ray.at(tx); + let py = y_ray.at(ty); let spp_scale = 0.125_f32.max((samples_per_pixel as Float).sqrt()); *dpdx = spp_scale * self.base().camera_transform.render_from_camera_vector( @@ -209,6 +251,8 @@ pub trait CameraTrait { } } +#[enum_dispatch(CameraTrait)] +#[derive(Debug)] pub enum Camera { Perspective(PerspectiveCamera), Orthographic(OrthographicCamera), @@ -216,39 +260,7 @@ pub enum Camera { Realistic(RealisticCamera), } -impl CameraTrait for Camera { - fn base(&self) -> &CameraBase { - match self { - Camera::Perspective(c) => c.base(), - Camera::Orthographic(c) => c.base(), - Camera::Spherical(c) => c.base(), - Camera::Realistic(c) => c.base(), - } - } - - fn generate_ray(&self, sample: CameraSample, lambda: &SampledWavelengths) -> Option { - match self { - Camera::Perspective(c) => c.generate_ray(sample, lambda), - Camera::Orthographic(c) => c.generate_ray(sample, lambda), - Camera::Spherical(c) => c.generate_ray(sample, lambda), - Camera::Realistic(c) => c.generate_ray(sample, lambda), - } - } - - fn generate_ray_differential( - &self, - sample: CameraSample, - lambda: &SampledWavelengths, - ) -> Option { - match self { - Camera::Perspective(c) => c.generate_ray_differential(sample, lambda), - Camera::Orthographic(c) => c.generate_ray_differential(sample, lambda), - Camera::Spherical(c) => c.generate_ray_differential(sample, lambda), - Camera::Realistic(c) => c.generate_ray_differential(sample, lambda), - } - } -} - +#[derive(Debug)] pub struct LensElementInterface { pub curvature_radius: Float, pub thickness: Float, diff --git a/src/camera/orthographic.rs b/src/camera/orthographic.rs index 15937e7..3ef514d 100644 --- a/src/camera/orthographic.rs +++ b/src/camera/orthographic.rs @@ -9,6 +9,7 @@ use crate::spectra::{SampledSpectrum, SampledWavelengths}; use crate::utils::sampling::sample_uniform_disk_concentric; use crate::utils::transform::Transform; +#[derive(Debug)] pub struct OrthographicCamera { pub base: CameraBase, pub screen_from_camera: Transform, @@ -71,6 +72,11 @@ impl CameraTrait for OrthographicCamera { fn base(&self) -> &CameraBase { &self.base } + + fn init_metadata(&self, metadata: &mut crate::image::ImageMetadata) { + self.base.init_metadata(metadata) + } + fn generate_ray( &self, sample: CameraSample, @@ -95,7 +101,7 @@ impl CameraTrait for OrthographicCamera { // Compute point on plane of focus let ft = self.focal_distance / ray.d.z(); - let p_focus = ray.evaluate(ft); + let p_focus = ray.at(ft); // Update ray for effect of lens ray.o = Point3f::new(p_lens.x(), p_lens.y(), 0.); diff --git a/src/camera/perspective.rs b/src/camera/perspective.rs index 456815c..120e21c 100644 --- a/src/camera/perspective.rs +++ b/src/camera/perspective.rs @@ -11,6 +11,7 @@ use crate::spectra::{SampledSpectrum, SampledWavelengths}; use crate::utils::sampling::sample_uniform_disk_concentric; use crate::utils::transform::Transform; +#[derive(Debug)] pub struct PerspectiveCamera { pub base: CameraBase, pub screen_from_camera: Transform, @@ -77,6 +78,11 @@ impl CameraTrait for PerspectiveCamera { fn base(&self) -> &CameraBase { &self.base } + + fn init_metadata(&self, metadata: &mut crate::image::ImageMetadata) { + self.base.init_metadata(metadata) + } + fn generate_ray( &self, sample: CameraSample, @@ -101,7 +107,7 @@ impl CameraTrait for PerspectiveCamera { // Compute point on plane of focus let ft = self.focal_distance / r.d.z(); - let p_focus = r.evaluate(ft); + let p_focus = r.at(ft); // Update ray for effect of lens r.o = Point3f::new(p_lens.x(), p_lens.y(), 0.); diff --git a/src/camera/realistic.rs b/src/camera/realistic.rs index 1188c8d..53676ea 100644 --- a/src/camera/realistic.rs +++ b/src/camera/realistic.rs @@ -10,6 +10,7 @@ use crate::spectra::{SampledSpectrum, SampledWavelengths}; use crate::utils::math::{quadratic, square}; use crate::utils::scattering::refract; +#[derive(Debug)] pub struct RealisticCamera { base: CameraBase, focus_distance: Float, @@ -216,7 +217,7 @@ impl RealisticCamera { } // Test intersection point against element aperture - let p_hit = r_lens.evaluate(t); + let p_hit = r_lens.at(t); if square(p_hit.x()) + square(p_hit.y()) > square(element.aperture_radius) { return None; } @@ -297,6 +298,10 @@ impl CameraTrait for RealisticCamera { &self.base } + fn init_metadata(&self, metadata: &mut crate::image::ImageMetadata) { + self.base.init_metadata(metadata) + } + fn generate_ray( &self, sample: CameraSample, diff --git a/src/camera/spherical.rs b/src/camera/spherical.rs index ad29039..5e6ae75 100644 --- a/src/camera/spherical.rs +++ b/src/camera/spherical.rs @@ -6,14 +6,15 @@ use crate::geometry::{Bounds2f, Point2f, Point3f, Ray, Vector3f, spherical_direc use crate::spectra::{SampledSpectrum, SampledWavelengths}; use crate::utils::math::{equal_area_square_to_sphere, wrap_equal_area_square}; -#[derive(PartialEq)] +#[derive(Debug, PartialEq)] pub struct EquiRectangularMapping; -#[derive(PartialEq)] +#[derive(Debug, PartialEq)] pub enum Mapping { EquiRectangular(EquiRectangularMapping), } +#[derive(Debug)] pub struct SphericalCamera { pub base: CameraBase, pub screen: Bounds2f, @@ -27,6 +28,10 @@ impl CameraTrait for SphericalCamera { &self.base } + fn init_metadata(&self, metadata: &mut crate::image::ImageMetadata) { + self.base.init_metadata(metadata) + } + fn generate_ray( &self, sample: CameraSample, diff --git a/src/core/aggregates.rs b/src/core/aggregates.rs index f8c99be..211f4d1 100644 --- a/src/core/aggregates.rs +++ b/src/core/aggregates.rs @@ -783,7 +783,11 @@ impl BVHAggregate { let node = &self.nodes[current_node_index]; // Check ray against BVH node bounds using the current closest hit_t - if node.bounds.intersect_p(r.o, hit_t, inv_dir, &dir_is_neg) { + if node + .bounds + .intersect_p(r.o, hit_t, inv_dir, &dir_is_neg) + .is_some() + { if node.n_primitives > 0 { // Intersect ray with all primitives in this leaf for i in 0..node.n_primitives { @@ -858,7 +862,11 @@ impl BVHAggregate { let node = &self.nodes[current_node_index]; // Check AABB - if node.bounds.intersect_p(r.o, t_max, inv_dir, &dir_is_neg) { + if node + .bounds + .intersect_p(r.o, t_max, inv_dir, &dir_is_neg) + .is_some() + { if node.n_primitives > 0 { for i in 0..node.n_primitives { let prim_idx = node.primitives_offset + i as usize; diff --git a/src/core/bssrdf.rs b/src/core/bssrdf.rs index df40db9..74d1344 100644 --- a/src/core/bssrdf.rs +++ b/src/core/bssrdf.rs @@ -7,11 +7,11 @@ use crate::utils::math::{catmull_rom_weights, square}; use crate::utils::sampling::sample_catmull_rom_2d; use std::sync::Arc; -#[derive(Clone, Debug)] -pub struct BSSRDFSample { +#[derive(Debug)] +pub struct BSSRDFSample<'a> { sp: SampledSpectrum, pdf: SampledSpectrum, - sw: BSDF, + sw: BSDF<'a>, wo: Vector3f, } @@ -98,7 +98,7 @@ pub struct BSSRDFProbeSegment { 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(si: &SubsurfaceInteraction) -> BSSRDFSample<'_>; } pub enum BSSRDF<'a> { @@ -214,7 +214,6 @@ impl<'a> TabulatedBSSRDF<'a> { for (j, rho_weight) in rho_weights.iter().enumerate() { if *rho_weight != 0. { // Update _rhoEff_ and _sr_ for wavelength - // We use 'j' for the offset calculation and 'rho_weight' for calculation rho_eff += self.table.rho_eff[rho_offset + j] * rho_weight; // Fix: Use .iter().enumerate() for 'k' @@ -280,7 +279,7 @@ impl<'a> BSSRDFTrait for TabulatedBSSRDF<'a> { }) } - fn probe_intersection_to_sample(_si: &SubsurfaceInteraction) -> BSSRDFSample { + fn probe_intersection_to_sample(_si: &SubsurfaceInteraction) -> BSSRDFSample<'_> { todo!() } } diff --git a/src/core/bxdf.rs b/src/core/bxdf.rs index c615298..1e88beb 100644 --- a/src/core/bxdf.rs +++ b/src/core/bxdf.rs @@ -539,20 +539,20 @@ impl<'a> BSDF<'a> { wo_render: Vector3f, wi_render: Vector3f, mode: TransportMode, - ) -> SampledSpectrum { + ) -> Option { let bxdf = match &self.bxdf { Some(b) => b, - None => return SampledSpectrum::default(), + None => return None, }; let wi = self.render_to_local(wi_render); let wo = self.render_to_local(wo_render); if wo.z() == 0.0 || wi.z() == 0.0 { - return SampledSpectrum::default(); + return None; } - bxdf.f(wo, wi, mode) + Some(bxdf.f(wo, wi, mode)) } pub fn sample_f( diff --git a/src/core/film.rs b/src/core/film.rs index 84dd6cb..998bd93 100644 --- a/src/core/film.rs +++ b/src/core/film.rs @@ -11,6 +11,7 @@ use crate::spectra::{ ConstantSpectrum, DenselySampledSpectrum, N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths, Spectrum, SpectrumTrait, cie_x, cie_y, cie_z, }; +use crate::utils::AtomicFloat; use crate::utils::color::{MatrixMulColor, RGB, SRGB, Triplet, XYZ, white_balance}; use crate::utils::colorspace::RGBColorSpace; use crate::utils::containers::Array2D; @@ -25,27 +26,37 @@ use once_cell::sync::Lazy; use std::error::Error; use std::sync::Mutex; +#[derive(Clone, Debug)] pub struct RGBFilm { pub base: FilmBase, pub max_component_value: Float, pub write_fp16: bool, pub filter_integral: Float, pub output_rgbf_from_sensor_rgb: SquareMatrix, - pub pixels: Array2D, + pub pixels: Arc>, } +#[derive(Debug)] pub struct RGBPixel { - rgb_sum: [f64; 3], - weight_sum: f64, - rgb_splat: Mutex<[f64; 3]>, + rgb_sum: [AtomicFloat; 3], + weight_sum: AtomicFloat, + rgb_splat: [AtomicFloat; 3], } impl Default for RGBPixel { fn default() -> Self { Self { - rgb_sum: [0., 0., 0.], - weight_sum: 0., - rgb_splat: Mutex::new([0., 0., 0.]), + rgb_sum: [ + AtomicFloat::new(0.0), + AtomicFloat::new(0.0), + AtomicFloat::new(0.0), + ], + weight_sum: AtomicFloat::new(0.0), + rgb_splat: [ + AtomicFloat::new(0.0), + AtomicFloat::new(0.0), + AtomicFloat::new(0.0), + ], } } } @@ -53,11 +64,10 @@ impl Default for RGBPixel { impl RGBFilm { pub fn new( base: FilmBase, - colorspace: RGBColorSpace, + colorspace: &RGBColorSpace, max_component_value: Float, write_fp16: bool, ) -> Self { - let pixels = Array2D::new(base.pixel_bounds); let filter_integral = base.filter.integral(); let sensor_matrix = base .sensor @@ -65,22 +75,35 @@ impl RGBFilm { .expect("Sensor must exist") .xyz_from_sensor_rgb; let output_rgbf_from_sensor_rgb = colorspace.rgb_from_xyz * sensor_matrix; + + let width = base.pixel_bounds.p_max.x() - base.pixel_bounds.p_min.x(); + let height = base.pixel_bounds.p_max.y() - base.pixel_bounds.p_min.y(); + let count = (width * height) as usize; + + let mut pixel_vec = Vec::with_capacity(count); + for _ in 0..count { + pixel_vec.push(RGBPixel::default()); + } + + let pixels_array = Array2D::new(base.pixel_bounds); + Self { base, max_component_value, write_fp16, filter_integral, output_rgbf_from_sensor_rgb, - pixels, + pixels: Arc::new(pixels_array), } } } +#[derive(Clone, Debug)] pub struct GBufferBFilm { pub base: FilmBase, output_from_render: AnimatedTransform, apply_inverse: bool, - pixels: Array2D, + pixels: Arc>, colorspace: RGBColorSpace, max_component_value: Float, write_fp16: bool, @@ -88,27 +111,30 @@ pub struct GBufferBFilm { output_rgbf_from_sensor_rgb: SquareMatrix, } +#[derive(Debug, Default)] struct GBufferPixel { - rgb_sum: [f64; 3], - weight_sum: f64, - g_bugger_weight_sum: f64, - rgb_splat: Mutex<[f64; 3]>, - p_sum: Point3f, - dz_dx_sum: Float, - dz_dy_sum: Float, - n_sum: Normal3f, - ns_sum: Normal3f, - uv_sum: Point2f, - rgb_albedo_sum: [f64; 3], - rgb_variance: VarianceEstimator, + pub rgb_sum: [AtomicFloat; 3], + pub weight_sum: AtomicFloat, + pub g_bugger_weight_sum: AtomicFloat, + pub rgb_splat: [AtomicFloat; 3], + pub p_sum: Point3f, + pub dz_dx_sum: AtomicFloat, + pub dz_dy_sum: Float, + pub n_sum: Normal3f, + pub ns_sum: Normal3f, + pub uv_sum: Point2f, + pub rgb_albedo_sum: [AtomicFloat; 3], + pub rgb_variance: VarianceEstimator, } +#[derive(Clone, Debug)] pub struct SpectralFilm { pub base: FilmBase, output_from_render: AnimatedTransform, pixels: Array2D, } +#[derive(Clone, Debug)] struct SpectralPixel; const N_SWATCH_REFLECTANCES: usize = 24; @@ -263,12 +289,12 @@ pub struct VisibleSurface { impl VisibleSurface { pub fn new( - _si: SurfaceInteraction, - albedo: SampledSpectrum, - _lambda: SampledWavelengths, + _si: &SurfaceInteraction, + albedo: &SampledSpectrum, + _lambda: &SampledWavelengths, ) -> Self { VisibleSurface { - albedo, + albedo: *albedo, ..Default::default() } } @@ -290,7 +316,7 @@ impl Default for VisibleSurface { } } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct FilmBase { pub full_resolution: Point2i, pub pixel_bounds: Bounds2i, @@ -352,7 +378,7 @@ pub trait FilmTrait: Sync { PixelFormat::F32 }; - let channel_names = vec!["R".to_string(), "G".to_string(), "B".to_string()]; + let channel_names = &["R", "G", "B"]; let pixel_bounds = self.base().pixel_bounds; let resolution = Point2i::from(pixel_bounds.diagonal()); @@ -461,6 +487,7 @@ pub trait FilmTrait: Sync { } } +#[derive(Clone, Debug)] pub enum Film { RGB(RGBFilm), GBuffer(GBufferBFilm), @@ -493,11 +520,11 @@ impl FilmTrait for RGBFilm { rgb *= self.max_component_value / m; } - let pixel = &mut self.pixels[p_film]; + let pixel = &self.pixels[p_film]; for c in 0..3 { - pixel.rgb_sum[c] += (weight * rgb[c]) as f64; + pixel.rgb_sum[c].add((weight * rgb[c]) as f64); } - pixel.weight_sum += weight as f64; + pixel.weight_sum.add(weight as f64); } fn add_splat(&mut self, p: Point2f, l: SampledSpectrum, lambda: &SampledWavelengths) { @@ -527,8 +554,7 @@ impl FilmTrait for RGBFilm { if wt != 0. { let pixel = &self.pixels[*pi]; for i in 0..3 { - let mut rgb_splat = pixel.rgb_splat.lock().unwrap(); - rgb_splat[i] += wt as f64 * rgb_splat[i]; + pixel.rgb_splat[i].add((wt * rgb[i]) as f64); } } } @@ -537,23 +563,24 @@ impl FilmTrait for RGBFilm { fn get_pixel_rgb(&self, p: Point2i, splat_scale: Option) -> RGB { let pixel = &self.pixels[p]; let mut rgb = RGB::new( - pixel.rgb_sum[0] as Float, - pixel.rgb_sum[1] as Float, - pixel.rgb_sum[2] as Float, + pixel.rgb_sum[0].load() as Float, + pixel.rgb_sum[1].load() as Float, + pixel.rgb_sum[2].load() as Float, ); - let weight_sum = pixel.weight_sum; + let weight_sum = pixel.weight_sum.load(); if weight_sum != 0. { rgb /= weight_sum as Float } - let rgb_splat = pixel.rgb_splat.lock().unwrap(); if let Some(splat) = splat_scale { for c in 0..3 { - rgb[c] += splat * rgb_splat[c] as Float / self.filter_integral; + let splat_val = pixel.rgb_splat[c].load(); + rgb[c] += splat * splat_val as Float / self.filter_integral; } } else { for c in 0..3 { - rgb[c] += rgb_splat[c] as Float / self.filter_integral; + let splat_val = pixel.rgb_splat[c].load(); + rgb[c] += splat_val as Float / self.filter_integral; } } self.output_rgbf_from_sensor_rgb.mul_rgb(rgb) @@ -615,8 +642,7 @@ impl FilmTrait for GBufferBFilm { if wt != 0. { let pixel = &self.pixels[*pi]; for i in 0..3 { - let mut rgb_splat = pixel.rgb_splat.lock().unwrap(); - rgb_splat[i] += wt as f64 * rgb_splat[i]; + pixel.rgb_splat[i].add((wt * rgb[i]) as f64); } } } @@ -633,23 +659,24 @@ impl FilmTrait for GBufferBFilm { fn get_pixel_rgb(&self, p: Point2i, splat_scale: Option) -> RGB { let pixel = &self.pixels[p]; let mut rgb = RGB::new( - pixel.rgb_sum[0] as Float, - pixel.rgb_sum[1] as Float, - pixel.rgb_sum[2] as Float, + pixel.rgb_sum[0].load() as Float, + pixel.rgb_sum[1].load() as Float, + pixel.rgb_sum[2].load() as Float, ); - let weight_sum = pixel.weight_sum; + let weight_sum = pixel.weight_sum.load(); if weight_sum != 0. { rgb /= weight_sum as Float } - let rgb_splat = pixel.rgb_splat.lock().unwrap(); if let Some(splat) = splat_scale { for c in 0..3 { - rgb[c] += splat * rgb_splat[c] as Float / self.filter_integral; + let splat_val = pixel.rgb_splat[c].load(); + rgb[c] += splat * splat_val as Float / self.filter_integral; } } else { for c in 0..3 { - rgb[c] += rgb_splat[c] as Float / self.filter_integral; + let splat_val = pixel.rgb_splat[c].load(); + rgb[c] += splat_val as Float / self.filter_integral; } } self.output_rgbf_from_sensor_rgb.mul_rgb(rgb) diff --git a/src/core/filter.rs b/src/core/filter.rs index e53f304..3bb3ba9 100644 --- a/src/core/filter.rs +++ b/src/core/filter.rs @@ -13,7 +13,7 @@ pub struct FilterSample { pub weight: Float, } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct FilterSampler { domain: Bounds2f, distrib: PiecewiseConstant2D, @@ -30,8 +30,8 @@ impl FilterSampler { Point2f::new(radius.x(), radius.y()), ); - let nx = (32.0 * radius.x()) as i32; - let ny = (32.0 * radius.y()) as i32; + let nx = (32.0 * radius.x()) as usize; + let ny = (32.0 * radius.y()) as usize; let mut f = Array2D::new_with_dims(nx, ny); for y in 0..f.y_size() { @@ -67,7 +67,7 @@ pub trait FilterTrait { } #[enum_dispatch(FilterTrait)] -#[derive(Debug)] +#[derive(Clone, Debug)] pub enum Filter { Box(BoxFilter), Gaussian(GaussianFilter), @@ -76,7 +76,7 @@ pub enum Filter { Triangle(TriangleFilter), } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct BoxFilter { pub radius: Vector2f, } @@ -113,7 +113,7 @@ impl FilterTrait for BoxFilter { } } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct GaussianFilter { pub radius: Vector2f, pub sigma: Float, @@ -165,7 +165,7 @@ impl FilterTrait for GaussianFilter { } } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct MitchellFilter { pub radius: Vector2f, pub b: Float, @@ -231,7 +231,7 @@ impl FilterTrait for MitchellFilter { } } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct LanczosSincFilter { pub radius: Vector2f, pub tau: Float, @@ -290,7 +290,7 @@ impl FilterTrait for LanczosSincFilter { } } -#[derive(Debug)] +#[derive(Clone, Debug)] pub struct TriangleFilter { pub radius: Vector2f, } diff --git a/src/core/integrator.rs b/src/core/integrator.rs deleted file mode 100644 index f16b913..0000000 --- a/src/core/integrator.rs +++ /dev/null @@ -1,13 +0,0 @@ -pub struct BPDTIntegrator; -pub struct MLTIntegrator; -pub struct SPPMIntegrator; -pub struct SamplerIntegrator; - -pub enum Integrator { - BPDT(BPDTIntegrator), - MLT(MLTIntegrator), - SPPM(SPPMIntegrator), - Sampler(SamplerIntegrator), -} - -impl Integrator {} diff --git a/src/core/interaction.rs b/src/core/interaction.rs index 3ed8cd1..45352af 100644 --- a/src/core/interaction.rs +++ b/src/core/interaction.rs @@ -1,15 +1,24 @@ use crate::camera::{Camera, CameraTrait}; -use crate::core::bxdf::{BSDF, BxDFFlags}; -use crate::core::material::Material; -use crate::core::medium::{Medium, MediumInterface}; +use crate::core::bxdf::{BSDF, BxDFFlags, DiffuseBxDF}; +use crate::core::material::{ + Material, MaterialEvalContext, MaterialTrait, NormalBumpEvalContext, bump_map, normal_map, +}; +use crate::core::medium::{Medium, MediumInterface, PhaseFunction}; +use crate::core::options::get_options; use crate::core::pbrt::{Float, clamp_t}; +use crate::core::sampler::{Sampler, SamplerTrait}; +use crate::core::texture::{FloatTexture, UniversalTextureEvaluator}; use crate::geometry::{ Normal3f, Point2f, Point3f, Point3fi, Ray, RayDifferential, Vector3f, VectorLike, }; -use crate::lights::Light; +use crate::image::Image; +use crate::lights::{Light, LightTrait}; use crate::shapes::Shape; +use crate::spectra::{SampledSpectrum, SampledWavelengths}; use crate::utils::math::{difference_of_products, square}; +use bumpalo::Bump; +use enum_dispatch::enum_dispatch; use std::any::Any; use std::sync::Arc; @@ -19,16 +28,14 @@ pub struct InteractionData { pub n: Normal3f, pub time: Float, pub wo: Vector3f, - pub medium_interface: Option>, + pub medium_interface: Option, pub medium: Option>, } -pub trait Interaction: Send + Sync { +#[enum_dispatch] +pub trait InteractionTrait: Send + Sync + std::fmt::Debug { fn get_common(&self) -> &InteractionData; fn get_common_mut(&mut self) -> &mut InteractionData; - fn as_any(&self) -> &dyn Any; - fn as_any_mut(&mut self) -> &mut dyn Any; - fn into_any(self: Box) -> Box; fn p(&self) -> Point3f { self.get_common().pi.into() @@ -48,31 +55,49 @@ pub trait Interaction: Send + Sync { } fn is_surface_interaction(&self) -> bool { - self.as_any().is::() + false } fn is_medium_interaction(&self) -> bool { - self.as_any().is::() + false } - fn get_medium(&self, w: Vector3f) -> Option>; - - fn spawn_ray(&self, d: Vector3f) -> Ray; - - fn spawn_ray_to_point(&self, p2: Point3f) -> Ray { - let origin = self.p(); - let direction = p2 - origin; - Ray { - o: origin, - d: direction, - time: self.time(), - medium: self.get_medium(direction), - differential: None, + fn get_medium(&self, w: Vector3f) -> Option> { + let data = self.get_common(); + if let Some(mi) = &data.medium_interface { + if w.dot(data.n.into()) > 0.0 { + mi.outside.clone() + } else { + mi.inside.clone() + } + } else { + data.medium.clone() } } - fn spawn_ray_to_interaction(&self, other: &dyn Interaction) -> Ray { - self.spawn_ray_to_point(other.p()) + fn spawn_ray(&self, d: Vector3f) -> Ray { + let data = self.get_common(); + let mut ray = Ray::spawn(&data.pi, &data.n, data.time, d); + + ray.medium = self.get_medium(d); + ray + } + + fn spawn_ray_to_point(&self, p2: Point3f) -> Ray { + let data = self.get_common(); + let mut ray = Ray::spawn_to_point(&data.pi, &data.n, data.time, p2); + ray.medium = self.get_medium(ray.d); + ray + } + + fn spawn_ray_to_interaction(&self, other: &dyn InteractionTrait) -> Ray { + let data = self.get_common(); + let other_data = other.get_common(); + + let mut ray = + Ray::spawn_to_interaction(&data.pi, &data.n, data.time, &other_data.pi, &other_data.n); + ray.medium = self.get_medium(ray.d); + ray } fn offset_ray_vector(&self, w: Vector3f) -> Point3f { @@ -84,6 +109,67 @@ pub trait Interaction: Send + Sync { } } +#[enum_dispatch(InteractionTrait)] +#[derive(Debug, Clone)] +pub enum Interaction { + Surface(SurfaceInteraction), + Medium(MediumInteraction), + Simple(SimpleInteraction), +} + +impl Interaction { + pub fn set_medium_interface(&mut self, mi: Option) { + match self { + Interaction::Surface(si) => si.common.medium_interface = mi, + Interaction::Simple(si) => si.common.medium_interface = mi, + Interaction::Medium(_) => {} // Medium interactions don't usually sit on boundaries + } + } +} + +#[derive(Debug, Clone)] +pub struct SimpleInteraction { + pub common: InteractionData, +} + +impl SimpleInteraction { + pub fn new(pi: Point3fi, time: Float, medium_interface: Option) -> Self { + Self { + common: InteractionData { + pi, + time, + medium_interface, + n: Normal3f::default(), + wo: Vector3f::default(), + medium: None, + }, + } + } + + pub fn new_interface(p: Point3f, medium_interface: Option) -> Self { + Self { + common: InteractionData { + pi: Point3fi::new_from_point(p), + n: Normal3f::zero(), + wo: Vector3f::zero(), + time: 0.0, + medium: None, + medium_interface, + }, + } + } +} + +impl InteractionTrait for SimpleInteraction { + fn get_common(&self) -> &InteractionData { + &self.common + } + + fn get_common_mut(&mut self) -> &mut InteractionData { + &mut self.common + } +} + #[derive(Default, Clone, Debug)] pub struct ShadingGeometry { pub n: Normal3f, @@ -93,8 +179,8 @@ pub struct ShadingGeometry { pub dndv: Normal3f, } -#[derive(Default, Debug)] -pub struct SurfaceInteraction<'a> { +#[derive(Debug, Default, Clone)] +pub struct SurfaceInteraction { pub common: InteractionData, pub uv: Point2f, pub dpdu: Vector3f, @@ -102,7 +188,6 @@ pub struct SurfaceInteraction<'a> { pub dndu: Normal3f, pub dndv: Normal3f, pub shading: ShadingGeometry, - pub medium_interface: Option>, pub face_index: usize, pub area_light: Option>, pub material: Option>, @@ -113,11 +198,18 @@ pub struct SurfaceInteraction<'a> { pub dudy: Float, pub dvdy: Float, pub shape: Arc, - // pub bsdf: Option>, } -impl<'a> SurfaceInteraction<'a> { - pub fn compute_differentials(&mut self, r: &Ray, camera: Camera, samples_per_pixel: i32) { +impl SurfaceInteraction { + pub fn le(&self, w: Vector3f, lambda: &SampledWavelengths) -> SampledSpectrum { + if let Some(area_light) = &self.area_light { + area_light.l(self.p(), self.n(), self.uv, w, lambda) + } else { + SampledSpectrum::new(0.) + } + } + + pub fn compute_differentials(&mut self, r: &Ray, camera: &Camera, samples_per_pixel: i32) { let computed = if let Some(diff) = &r.differential { let dot_rx = self.common.n.dot(diff.rx_direction.into()); let dot_ry = self.common.n.dot(diff.ry_direction.into()); @@ -211,25 +303,68 @@ impl<'a> SurfaceInteraction<'a> { } } - pub fn spawn_ray(&self, d: Vector3f) -> Ray { - let mut ray = Ray::spawn(&self.pi(), &self.n(), self.time(), d); - ray.medium = self.get_medium(d); - ray + pub fn get_bsdf<'a>( + &mut self, + r: &Ray, + lambda: &SampledWavelengths, + camera: &Camera, + scratch: &'a Bump, + sampler: &mut Sampler, + ) -> Option> { + self.compute_differentials(r, camera, sampler.samples_per_pixel() as i32); + + 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; + let displacement = material.get_displacement(); + let normal_map = material.get_normal_map(); + if displacement.is_some() || normal_map.is_some() { + // This calls the function defined above + self.compute_bump_geometry(&tex_eval, displacement, normal_map); + } + + let mut bsdf = material.get_bxdf(&tex_eval, &ctx, lambda, scratch); + if get_options().force_diffuse { + let r = bsdf.rho_wo(self.common.wo, &[sampler.get1d()], &[sampler.get2d()]); + let diff_bxdf = scratch.alloc(DiffuseBxDF::new(r)); + bsdf = BSDF::new(self.shading.n, self.shading.dpdu, Some(diff_bxdf)); + } + Some(bsdf) } - pub fn spawn_ray_to(&self, p2: Point3f) -> Ray { - let mut ray = Ray::spawn_to_point(&self.pi(), &self.n(), self.time(), p2); - ray.medium = self.get_medium(ray.d); - ray - } + fn compute_bump_geometry( + &mut self, + tex_eval: &UniversalTextureEvaluator, + displacement: Option, + normal_image: Option<&Image>, + ) { + let ctx = NormalBumpEvalContext::from(&*self); + let (dpdu, dpdv) = if let Some(disp) = displacement { + bump_map(tex_eval, &disp, &ctx) + } else if let Some(map) = normal_image { + normal_map(map, &ctx) + } else { + (self.shading.dpdu, self.shading.dpdv) + }; - // fn get_medium(&self, w: &Vector3f) -> Option> { - // if w.dot(self.n().into()) > 0.0 { - // self.medium_interface.outside.clone() - // } else { - // self.medium_interface.inside.clone() - // } - // } + let mut ns = Normal3f::from(dpdu.cross(dpdv).normalize()); + if ns.dot(self.n()) < 0.0 { + ns = -ns; + } + + self.set_shading_geometry(ns, dpdu, dpdv, self.shading.dndu, self.shading.dndv, false); + } pub fn spawn_ray_with_differentials( &self, @@ -326,28 +461,7 @@ impl<'a> SurfaceInteraction<'a> { } } -#[derive(Default, Debug, Clone)] -pub struct PhaseFunction; - -pub struct MediumInteraction { - pub common: InteractionData, - pub medium: Arc, - pub phase: PhaseFunction, -} - -impl Interaction for SurfaceInteraction { - fn as_any(&self) -> &dyn Any { - self - } - - fn as_any_mut(&mut self) -> &mut dyn Any { - self - } - - fn into_any(self: Box) -> Box { - self - } - +impl InteractionTrait for SurfaceInteraction { fn get_common(&self) -> &InteractionData { &self.common } @@ -357,7 +471,7 @@ impl Interaction for SurfaceInteraction { } fn get_medium(&self, w: Vector3f) -> Option> { - self.medium_interface.as_ref().and_then(|interface| { + self.common.medium_interface.as_ref().and_then(|interface| { if self.n().dot(w.into()) > 0.0 { interface.outside.clone() } else { @@ -366,34 +480,8 @@ impl Interaction for SurfaceInteraction { }) } - fn spawn_ray(&self, d: Vector3f) -> Ray { - let mut ray = Ray::spawn(&self.pi(), &self.n(), self.time(), d); - ray.medium = self.get_medium(d); - ray - } - - fn spawn_ray_to_point(&self, p2: Point3f) -> Ray { - let mut ray = Ray::spawn_to_point(&self.pi(), &self.n(), self.time(), p2); - ray.medium = self.get_medium(ray.d); - ray - } - - fn spawn_ray_to_interaction(&self, other: &dyn Interaction) -> Ray { - // Check if the other interaction is a surface to use the robust spawn method - if let Some(si_to) = other.as_any().downcast_ref::() { - let mut ray = Ray::spawn_to_interaction( - &self.pi(), - &self.n(), - self.time(), - &si_to.pi(), - &si_to.n(), - ); - ray.medium = self.get_medium(ray.d); - ray - } else { - // Fallback for non-surface interactions - self.spawn_ray_to_point(other.p()) - } + fn is_surface_interaction(&self) -> bool { + true } } @@ -438,7 +526,6 @@ impl SurfaceInteraction { dndu, dndv, }, - medium_interface: None, material: None, face_index: 0, area_light: None, @@ -449,7 +536,6 @@ impl SurfaceInteraction { dvdx: 0.0, dvdy: 0.0, shape: Arc::new(Shape::default()), - bsdf: None, } } @@ -504,11 +590,22 @@ impl SurfaceInteraction { } } + pub fn new_minimal(pi: Point3fi, uv: Point2f) -> Self { + Self { + common: InteractionData { + pi, + ..Default::default() + }, + uv, + ..Default::default() + } + } + pub fn set_intersection_properties( &mut self, mtl: Arc, area: Arc, - prim_medium_interface: Option>, + prim_medium_interface: Option, ray_medium: Arc, ) { self.material = Some(mtl); @@ -524,17 +621,16 @@ impl SurfaceInteraction { } } -impl Interaction for MediumInteraction { - fn as_any(&self) -> &dyn Any { - self - } +#[derive(Clone, Debug)] +pub struct MediumInteraction { + pub common: InteractionData, + pub medium: Arc, + pub phase: PhaseFunction, +} - fn as_any_mut(&mut self) -> &mut dyn Any { - self - } - - fn into_any(self: Box) -> Box { - self +impl InteractionTrait for MediumInteraction { + fn is_medium_interaction(&self) -> bool { + true } fn get_common(&self) -> &InteractionData { @@ -544,18 +640,4 @@ impl Interaction for MediumInteraction { fn get_common_mut(&mut self) -> &mut InteractionData { &mut self.common } - - fn get_medium(&self, _w: Vector3f) -> Option> { - Some(self.medium.clone()) - } - - fn spawn_ray(&self, d: Vector3f) -> Ray { - Ray { - o: self.p(), - d, - time: self.time(), - medium: Some(self.medium.clone()), - differential: None, - } - } } diff --git a/src/core/material.rs b/src/core/material.rs index b3e3ae4..65c68ca 100644 --- a/src/core/material.rs +++ b/src/core/material.rs @@ -4,18 +4,23 @@ use std::ops::Deref; use std::sync::Arc; use crate::core::bssrdf::BSSRDF; -use crate::core::bxdf::{BSDF, DiffuseBxDF}; +use crate::core::bxdf::{BSDF, DielectricBxDF, DiffuseBxDF}; +use crate::core::interaction::InteractionTrait; +use crate::core::interaction::{Interaction, ShadingGeometry, SurfaceInteraction}; +use crate::core::pbrt::Float; use crate::core::texture::{FloatTexture, SpectrumTexture, TextureEvalContext, TextureEvaluator}; -use crate::geometry::{Normal3f, Vector3f}; -use crate::image::Image; -use crate::spectra::SampledWavelengths; +use crate::geometry::{Frame, Normal3f, Point2f, Point3f, Vector2f, Vector3f, VectorLike}; +use crate::image::{Image, WrapMode, WrapMode2D}; +use crate::spectra::{SampledWavelengths, Spectrum, SpectrumTrait}; +use crate::utils::hash::hash_float; +use crate::utils::scattering::TrowbridgeReitzDistribution; #[derive(Clone, Debug)] pub struct MaterialEvalContext { - texture: TextureEvalContext, - wo: Vector3f, - ns: Normal3f, - dpdus: Vector3f, + pub texture: TextureEvalContext, + pub wo: Vector3f, + pub ns: Normal3f, + pub dpdus: Vector3f, } impl Deref for MaterialEvalContext { @@ -26,9 +31,124 @@ impl Deref for MaterialEvalContext { } } +impl From<&SurfaceInteraction> for MaterialEvalContext { + fn from(si: &SurfaceInteraction) -> Self { + Self { + texture: TextureEvalContext::from(si), + wo: si.common.wo, + ns: si.shading.n, + dpdus: si.shading.dpdu, + } + } +} + +#[derive(Clone, Debug, Default)] +pub struct NormalBumpEvalContext { + p: Point3f, + uv: Point2f, + n: Normal3f, + shading: ShadingGeometry, + dpdx: Vector3f, + dpdy: Vector3f, + // All 0 + dudx: Float, + dudy: Float, + dvdx: Float, + dvdy: Float, + face_index: usize, +} + +impl From<&SurfaceInteraction> for NormalBumpEvalContext { + fn from(si: &SurfaceInteraction) -> Self { + Self { + p: si.p(), + uv: si.uv, + n: si.n(), + shading: si.shading.clone(), + dudx: si.dudx, + dudy: si.dudy, + dvdx: si.dvdx, + dvdy: si.dvdy, + dpdx: si.dpdx, + dpdy: si.dpdy, + face_index: si.face_index, + } + } +} + +impl From<&NormalBumpEvalContext> for TextureEvalContext { + fn from(ctx: &NormalBumpEvalContext) -> Self { + Self { + p: ctx.p, + uv: ctx.uv, + n: ctx.n, + dpdx: ctx.dpdx, + dpdy: ctx.dpdy, + dudx: ctx.dudx, + dudy: ctx.dudy, + dvdx: ctx.dvdx, + dvdy: ctx.dvdy, + face_index: ctx.face_index, + } + } +} + +pub fn normal_map(normal_map: &Image, ctx: &NormalBumpEvalContext) -> (Vector3f, Vector3f) { + let wrap = WrapMode2D::from(WrapMode::Repeat); + let uv = Point2f::new(ctx.uv[0], 1. - ctx.uv[1]); + let r = normal_map.bilerp_channel_with_wrap(uv, 0, wrap); + let g = normal_map.bilerp_channel_with_wrap(uv, 1, wrap); + let b = normal_map.bilerp_channel_with_wrap(uv, 2, wrap); + let mut ns = Vector3f::new(2.0 * r - 1.0, 2.0 * g - 1.0, 2.0 * b - 1.0); + ns = ns.normalize(); + let frame = Frame::from_xz(ctx.shading.dpdu.normalize(), Vector3f::from(ctx.shading.n)); + ns = frame.from_local(ns); + let ulen = ctx.shading.dpdu.norm(); + let vlen = ctx.shading.dpdv.norm(); + let dpdu = ctx.shading.dpdu.gram_schmidt(ns).normalize() * ulen; + let dpdv = ctx.shading.dpdu.cross(dpdu).normalize() * vlen; + (dpdu, dpdv) +} + +pub fn bump_map( + tex_eval: &T, + displacement: &FloatTexture, + ctx: &NormalBumpEvalContext, +) -> (Vector3f, Vector3f) { + debug_assert!(tex_eval.can_evaluate(&[displacement], &[])); + let mut du = 0.5 * (ctx.dudx.abs() + ctx.dudy.abs()); + if du == 0.0 { + du = 0.0005; + } + + let mut dv = 0.5 * (ctx.dvdx.abs() + ctx.dvdy.abs()); + if dv == 0.0 { + dv = 0.0005; + } + let mut shifted_ctx = TextureEvalContext::from(ctx); + shifted_ctx.p = ctx.p + ctx.shading.dpdu * du; + shifted_ctx.uv = ctx.uv + Vector2f::new(du, 0.0); + let u_displace = tex_eval.evaluate_float(displacement, &shifted_ctx); + shifted_ctx.p = ctx.p + ctx.shading.dpdv * dv; + shifted_ctx.uv = ctx.uv + Vector2f::new(0.0, dv); + let v_displace = tex_eval.evaluate_float(displacement, &shifted_ctx); + let center_ctx = TextureEvalContext::from(ctx); + let displace = tex_eval.evaluate_float(displacement, ¢er_ctx); + + let d_displace_du = (u_displace - displace) / du; + let d_displace_dv = (v_displace - displace) / dv; + let n_vec = Vector3f::from(ctx.shading.n); + let dndu_vec = Vector3f::from(ctx.shading.dndu); + let dndv_vec = Vector3f::from(ctx.shading.dndv); + + let dpdu = ctx.shading.dpdu + n_vec * d_displace_du + dndu_vec * displace; + let dpdv = ctx.shading.dpdv + n_vec * d_displace_dv + dndv_vec * displace; + (dpdu, dpdv) +} + #[enum_dispatch] pub trait MaterialTrait: Send + Sync + std::fmt::Debug { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, tex_eval: &T, ctx: &MaterialEvalContext, @@ -40,12 +160,12 @@ pub trait MaterialTrait: Send + Sync + std::fmt::Debug { &self, tex_eval: &T, ctx: &MaterialEvalContext, - lambda: &SampledWavelengths, - ) -> BSSRDF<'a>; + lambda: &mut SampledWavelengths, + ) -> Option>; fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool; - fn get_normal_map(&self) -> Image; - fn get_displacement(&self) -> FloatTexture; + fn get_normal_map(&self) -> Option<&Image>; + fn get_displacement(&self) -> Option; fn has_surface_scattering(&self) -> bool; } @@ -68,7 +188,7 @@ pub enum Material { #[derive(Clone, Debug)] pub struct CoatedDiffuseMaterial; impl MaterialTrait for CoatedDiffuseMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -81,27 +201,31 @@ impl MaterialTrait for CoatedDiffuseMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } + fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { todo!() } } + #[derive(Clone, Debug)] pub struct CoatedConductorMaterial; impl MaterialTrait for CoatedConductorMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -114,17 +238,18 @@ impl MaterialTrait for CoatedConductorMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -134,7 +259,7 @@ impl MaterialTrait for CoatedConductorMaterial { #[derive(Clone, Debug)] pub struct ConductorMaterial; impl MaterialTrait for ConductorMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -147,17 +272,17 @@ impl MaterialTrait for ConductorMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -165,36 +290,69 @@ impl MaterialTrait for ConductorMaterial { } } #[derive(Clone, Debug)] -pub struct DielectricMaterial; +pub struct DielectricMaterial { + normal_map: Option>, + displacement: FloatTexture, + u_roughness: FloatTexture, + v_roughness: FloatTexture, + remap_roughness: bool, + eta: Spectrum, +} + impl MaterialTrait for DielectricMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, - _tex_eval: &T, - _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - _scratch: &'a Bump, + tex_eval: &T, + ctx: &MaterialEvalContext, + lambda: &SampledWavelengths, + scratch: &'a Bump, ) -> BSDF<'a> { - todo!() + let mut sampled_eta = self.eta.evaluate(lambda[0]); + if !self.eta.is_constant() { + lambda.terminate_secondary(); + } + + if sampled_eta == 0.0 { + sampled_eta = 1.0; + } + + let mut u_rough = tex_eval.evaluate_float(&self.u_roughness, ctx); + let mut v_rough = tex_eval.evaluate_float(&self.v_roughness, ctx); + + if self.remap_roughness { + u_rough = TrowbridgeReitzDistribution::roughness_to_alpha(u_rough); + v_rough = TrowbridgeReitzDistribution::roughness_to_alpha(v_rough); + } + + let distrib = TrowbridgeReitzDistribution::new(u_rough, v_rough); + let bxdf = scratch.alloc(DielectricBxDF::new(sampled_eta, distrib)); + + BSDF::new(ctx.ns, ctx.dpdus, Some(bxdf)) } + fn get_bssrdf<'a, T>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { - todo!() + _lambda: &mut SampledWavelengths, + ) -> Option> { + None } - fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { - todo!() + + fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool { + tex_eval.can_evaluate(&[&self.u_roughness, &self.v_roughness], &[]) } - fn get_normal_map(&self) -> Image { - todo!() + + fn get_normal_map(&self) -> Option<&Image> { + self.normal_map.as_deref() } - fn get_displacement(&self) -> FloatTexture { - todo!() + + fn get_displacement(&self) -> Option { + Some(self.displacement.clone()) } + fn has_surface_scattering(&self) -> bool { - todo!() + false } } #[derive(Clone, Debug)] @@ -205,7 +363,7 @@ pub struct DiffuseMaterial { } impl MaterialTrait for DiffuseMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, tex_eval: &T, ctx: &MaterialEvalContext, @@ -221,28 +379,29 @@ impl MaterialTrait for DiffuseMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } - fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { - todo!() + fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool { + tex_eval.can_evaluate(&[], &[&self.reflectance]) } - fn get_normal_map(&self) -> Image { - todo!() + + fn get_normal_map(&self) -> Option<&Image> { + self.normal_map.as_deref() } - fn get_displacement(&self) -> FloatTexture { - todo!() + fn get_displacement(&self) -> Option { + Some(self.displacement.clone()) } fn has_surface_scattering(&self) -> bool { - todo!() + false } } #[derive(Clone, Debug)] pub struct DiffuseTransmissionMaterial; impl MaterialTrait for DiffuseTransmissionMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -255,18 +414,18 @@ impl MaterialTrait for DiffuseTransmissionMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -276,7 +435,7 @@ impl MaterialTrait for DiffuseTransmissionMaterial { #[derive(Clone, Debug)] pub struct HairMaterial; impl MaterialTrait for HairMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -289,18 +448,18 @@ impl MaterialTrait for HairMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -311,7 +470,7 @@ impl MaterialTrait for HairMaterial { #[derive(Clone, Debug)] pub struct MeasuredMaterial; impl MaterialTrait for MeasuredMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -324,18 +483,18 @@ impl MaterialTrait for MeasuredMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -345,7 +504,7 @@ impl MaterialTrait for MeasuredMaterial { #[derive(Clone, Debug)] pub struct SubsurfaceMaterial; impl MaterialTrait for SubsurfaceMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -358,18 +517,18 @@ impl MaterialTrait for SubsurfaceMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -379,7 +538,7 @@ impl MaterialTrait for SubsurfaceMaterial { #[derive(Clone, Debug)] pub struct ThinDielectricMaterial; impl MaterialTrait for ThinDielectricMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( + fn get_bxdf<'a, T: TextureEvaluator>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, @@ -392,17 +551,17 @@ impl MaterialTrait for ThinDielectricMaterial { &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { + _lambda: &mut SampledWavelengths, + ) -> Option> { todo!() } fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { todo!() } - fn get_normal_map(&self) -> Image { + fn get_normal_map(&self) -> Option<&Image> { todo!() } - fn get_displacement(&self) -> FloatTexture { + fn get_displacement(&self) -> Option { todo!() } fn has_surface_scattering(&self) -> bool { @@ -410,36 +569,74 @@ impl MaterialTrait for ThinDielectricMaterial { } } #[derive(Clone, Debug)] -pub struct MixMaterial; -impl MaterialTrait for MixMaterial { - fn get_bsdf<'a, T: TextureEvaluator>( +pub struct MixMaterial { + pub amount: FloatTexture, + pub materials: [Box; 2], +} + +impl MixMaterial { + pub fn choose_material( &self, - _tex_eval: &T, - _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - _scratch: &'a Bump, + tex_eval: &T, + ctx: &MaterialEvalContext, + ) -> &Material { + let amt = tex_eval.evaluate_float(&self.amount, ctx); + + if amt <= 0.0 { + return &self.materials[0]; + } + if amt >= 1.0 { + return &self.materials[1]; + } + + let u = hash_float(&(ctx.p, ctx.wo)); + + if amt < u { + &self.materials[0] + } else { + &self.materials[1] + } + } +} + +impl MaterialTrait for MixMaterial { + fn get_bxdf<'a, T: TextureEvaluator>( + &self, + tex_eval: &T, + ctx: &MaterialEvalContext, + lambda: &SampledWavelengths, + scratch: &'a Bump, ) -> BSDF<'a> { - todo!() + let chosen_mat = self.choose_material(tex_eval, ctx); + chosen_mat.get_bxdf(tex_eval, ctx, lambda, scratch) } fn get_bssrdf<'a, T>( &self, _tex_eval: &T, _ctx: &MaterialEvalContext, - _lambda: &SampledWavelengths, - ) -> BSSRDF<'a> { - todo!() + _lambda: &mut SampledWavelengths, + ) -> Option> { + None } - fn can_evaluate_textures(&self, _tex_eval: &dyn TextureEvaluator) -> bool { - todo!() + + fn can_evaluate_textures(&self, tex_eval: &dyn TextureEvaluator) -> bool { + tex_eval.can_evaluate(&[&self.amount], &[]) } - fn get_normal_map(&self) -> Image { - todo!() + + fn get_normal_map(&self) -> Option<&Image> { + None } - fn get_displacement(&self) -> FloatTexture { - todo!() + + fn get_displacement(&self) -> Option { + None + // panic!( + // "MixMaterial::get_displacement() shouldn't be called. \ + // Displacement is not supported on Mix materials directly." + // ); } + fn has_surface_scattering(&self) -> bool { - todo!() + false } } diff --git a/src/core/medium.rs b/src/core/medium.rs index bfc0715..bbb695e 100644 --- a/src/core/medium.rs +++ b/src/core/medium.rs @@ -1,15 +1,401 @@ +use bumpalo::Bump; +use enum_dispatch::enum_dispatch; use std::sync::Arc; +use crate::core::pbrt::{Float, INV_4_PI, PI, clamp_t}; +use crate::geometry::{ + Bounds3f, Frame, Point2f, Point3f, Point3i, Ray, Vector3f, VectorLike, spherical_direction, +}; +use crate::spectra::{ + BlackbodySpectrum, DenselySampledSpectrum, LAMBDA_MAX, LAMBDA_MIN, RGBIlluminantSpectrum, + RGBUnboundedSpectrum, SampledSpectrum, SampledWavelengths, Spectrum, SpectrumTrait, +}; +use crate::utils::containers::SampledGrid; +use crate::utils::math::square; +use crate::utils::rng::Rng; +use crate::utils::transform::Transform; + +#[derive(Debug, Clone, Copy)] +pub struct PhaseFunctionSample { + pub p: Float, + pub wi: Vector3f, + pub pdf: Float, +} + +impl PhaseFunctionSample { + pub fn new(p: Float, wi: Vector3f, pdf: Float) -> Self { + Self { p, wi, pdf } + } +} + +#[enum_dispatch] +pub trait PhaseFunctionTrait { + fn p(&self, wo: Vector3f, wi: Vector3f) -> Float; + fn sample_p(&self, wo: Vector3f, u: Point2f) -> Option; + fn pdf(&self, wo: Vector3f, wi: Vector3f) -> Float; +} + +#[enum_dispatch(PhaseFunctionTrait)] #[derive(Debug, Clone)] -pub struct HomogeneousMedium; +pub enum PhaseFunction { + HenyeyGreenstein(HGPhaseFunction), +} + +#[derive(Debug, Clone, Copy)] +pub struct HGPhaseFunction { + g: Float, +} + +impl HGPhaseFunction { + pub fn new(g: Float) -> Self { + Self { g } + } + + fn phase_hg(&self, cos_theta: Float) -> Float { + let denom = 1.0 + square(self.g) + 2.0 * self.g * cos_theta; + INV_4_PI * (1.0 - square(self.g)) / (denom * denom.sqrt()) + } +} + +impl PhaseFunctionTrait for HGPhaseFunction { + fn p(&self, wo: Vector3f, wi: Vector3f) -> Float { + self.phase_hg(-(wo.dot(wi))) + } + + fn sample_p(&self, wo: Vector3f, u: Point2f) -> Option { + let cos_theta = if self.g.abs() < 1e-3 { + 1.0 - 2.0 * u[0] + } else { + let square_term = (1.0 - square(self.g)) / (1.0 + self.g * (1.0 - 2.0 * u[0])); + -(1.0 + square(self.g) - square(square_term)) / (2.0 * self.g) + }; + + let sin_theta = (1.0 - square(cos_theta)).max(0.0).sqrt(); + let phi = 2.0 * PI * u[1]; + let w_frame = Frame::from_z(wo); + let wi = w_frame.from_local(spherical_direction(sin_theta, cos_theta, phi)); + + let p_val = self.phase_hg(cos_theta); + + Some(PhaseFunctionSample::new(p_val, wi, p_val)) + } + + fn pdf(&self, wo: Vector3f, wi: Vector3f) -> Float { + self.p(wo, wi) + } +} + #[derive(Debug, Clone)] -pub struct GridMedium; -#[derive(Debug, Clone)] -pub struct RGBGridMedium; -#[derive(Debug, Clone)] -pub struct CloudMedium; -#[derive(Debug, Clone)] -pub struct NanoVDBMedium; +pub struct MajorantGrid { + pub bounds: Bounds3f, + pub res: Point3i, + pub voxels: Vec, +} + +impl MajorantGrid { + pub fn new(bounds: Bounds3f, res: Point3i) -> Self { + Self { + bounds, + res, + voxels: Vec::with_capacity((res.x() * res.y() * res.z()) as usize), + } + } + pub fn lookup(&self, x: i32, y: i32, z: i32) -> Float { + let idx = z * self.res.x() * self.res.y() + y * self.res.x() + x; + if idx >= 0 && (idx as usize) < self.voxels.len() { + self.voxels[idx as usize] + } else { + 0.0 + } + } + pub fn set(&mut self, x: i32, y: i32, z: i32, v: Float) { + let idx = z * self.res.x() * self.res.y() + y * self.res.x() + x; + if idx >= 0 && (idx as usize) < self.voxels.len() { + self.voxels[idx as usize] = v; + } + } + + pub fn voxel_bounds(&self, x: i32, y: i32, z: i32) -> Bounds3f { + let p0 = Point3f::new( + x as Float / self.res.x() as Float, + y as Float / self.res.y() as Float, + z as Float / self.res.z() as Float, + ); + let p1 = Point3f::new( + (x + 1) as Float / self.res.x() as Float, + (y + 1) as Float / self.res.y() as Float, + (z + 1) as Float / self.res.z() as Float, + ); + Bounds3f::from_points(p0, p1) + } +} + +#[derive(Clone, Copy, Debug)] +pub struct RayMajorantSegment { + t_min: Float, + t_max: Float, + sigma_maj: SampledSpectrum, +} + +pub enum RayMajorantIterator { + Homogeneous(HomogeneousMajorantIterator), + DDA(DDAMajorantIterator), + // Grid(GridMajorantIterator<'a>), + Void, +} + +impl Iterator for RayMajorantIterator { + type Item = RayMajorantSegment; + + fn next(&mut self) -> Option { + match self { + RayMajorantIterator::Homogeneous(iter) => iter.next(), + RayMajorantIterator::DDA(iter) => iter.next(), + RayMajorantIterator::Void => None, + } + } +} +pub struct HomogeneousMajorantIterator { + called: bool, + seg: RayMajorantSegment, +} + +impl HomogeneousMajorantIterator { + pub fn new(t_min: Float, t_max: Float, sigma_maj: SampledSpectrum) -> Self { + Self { + seg: RayMajorantSegment { + t_min, + t_max, + sigma_maj, + }, + called: false, + } + } +} + +impl Iterator for HomogeneousMajorantIterator { + type Item = RayMajorantSegment; + + fn next(&mut self) -> Option { + if self.called { + return None; + } + + self.called = true; + Some(self.seg) + } +} + +pub struct DDAMajorantIterator { + sigma_t: SampledSpectrum, + t_min: Float, + t_max: Float, + grid: MajorantGrid, + next_crossing_t: [Float; 3], + delta_t: [Float; 3], + step: [i32; 3], + voxel_limit: [i32; 3], + voxel: [i32; 3], +} +impl DDAMajorantIterator { + pub fn new( + ray: &Ray, + t_min: Float, + t_max: Float, + grid: &MajorantGrid, + sigma_t: &SampledSpectrum, + ) -> Self { + let mut iter = Self { + t_min, + t_max, + sigma_t: *sigma_t, + grid: grid.clone(), + next_crossing_t: [0.0; 3], + delta_t: [0.0; 3], + step: [0; 3], + voxel_limit: [0; 3], + voxel: [0; 3], + }; + + let diag = grid.bounds.diagonal(); + + let mut ray_grid_d = Vector3f::new( + ray.d.x() / diag.x(), + ray.d.y() / diag.y(), + ray.d.z() / diag.z(), + ); + + let p_grid_start = grid.bounds.offset(&ray.at(t_min)); + let grid_intersect = Vector3f::from(p_grid_start); + + for axis in 0..3 { + iter.voxel[axis] = clamp_t( + (grid_intersect[axis] * grid.res[axis] as Float) as i32, + 0, + grid.res[axis] - 1, + ); + + iter.delta_t[axis] = 1.0 / (ray_grid_d[axis].abs() * grid.res[axis] as Float); + + if ray_grid_d[axis] == -0.0 { + ray_grid_d[axis] = 0.0; + } + + if ray_grid_d[axis] >= 0.0 { + let next_voxel_pos = (iter.voxel[axis] + 1) as Float / grid.res[axis] as Float; + iter.next_crossing_t[axis] = + t_min + (next_voxel_pos - grid_intersect[axis]) / ray_grid_d[axis]; + iter.step[axis] = 1; + iter.voxel_limit[axis] = grid.res[axis]; + } else { + let next_voxel_pos = (iter.voxel[axis]) as Float / grid.res[axis] as Float; + iter.next_crossing_t[axis] = + t_min + (next_voxel_pos - grid_intersect[axis]) / ray_grid_d[axis]; + iter.step[axis] = -1; + iter.voxel_limit[axis] = -1; + } + } + + iter + } +} + +impl Iterator for DDAMajorantIterator { + type Item = RayMajorantSegment; + + fn next(&mut self) -> Option { + if self.t_min >= self.t_max { + return None; + } + + // Find stepAxis for stepping to next voxel and exit point tVoxelExit + let d0 = self.next_crossing_t[0]; + let d1 = self.next_crossing_t[1]; + let d2 = self.next_crossing_t[2]; + + let bits = ((d0 < d1) as usize) << 2 | ((d0 < d2) as usize) << 1 | ((d1 < d2) as usize); + + const CMP_TO_AXIS: [usize; 8] = [2, 1, 2, 1, 2, 2, 0, 0]; + let step_axis = CMP_TO_AXIS[bits]; + + let t_voxel_exit = self.t_max.min(self.next_crossing_t[step_axis]); + // Get maxDensity for current voxel and initialize RayMajorantSegment, seg + + let density_scale = self + .grid + .lookup(self.voxel[0], self.voxel[1], self.voxel[2]); + + let sigma_maj = self.sigma_t * density_scale; + + let seg = RayMajorantSegment { + t_min: self.t_min, + t_max: t_voxel_exit, + sigma_maj, + }; + + // Advance to next voxel in maximum density grid + self.t_min = t_voxel_exit; + + if self.next_crossing_t[step_axis] > self.t_max { + self.t_min = self.t_max; + } + + self.voxel[step_axis] += self.step[step_axis]; + + if self.voxel[step_axis] == self.voxel_limit[step_axis] { + self.t_min = self.t_max; + } + + // Increment the crossing time for this axis + self.next_crossing_t[step_axis] += self.delta_t[step_axis]; + + Some(seg) + } +} + +pub struct MediumProperties { + sigma_a: SampledSpectrum, + sigma_s: SampledSpectrum, + phase: PhaseFunction, + le: SampledSpectrum, +} + +impl MediumProperties { + pub fn sigma_t(&self) -> SampledSpectrum { + self.sigma_a * self.sigma_s + } +} + +#[enum_dispatch] +pub trait MediumTrait: Send + Sync + std::fmt::Debug { + fn is_emissive(&self) -> bool; + fn sample_point(&self, p: Point3f, lambda: &SampledWavelengths) -> MediumProperties; + fn sample_ray( + &self, + ray: &Ray, + t_max: Float, + lambda: &SampledWavelengths, + buf: &Bump, + ) -> RayMajorantIterator; + + fn sample_t_maj( + &self, + mut ray: Ray, + mut t_max: Float, + u: Float, + _rng: &mut Rng, + lambda: &SampledWavelengths, + mut callback: F, + ) -> SampledSpectrum + where + F: FnMut(Point3f, MediumProperties, SampledSpectrum, SampledSpectrum) -> bool, + { + let len = ray.d.norm(); + t_max *= len; + ray.d /= len; + + let buf = Bump::new(); + let mut iter = self.sample_ray(&ray, t_max, lambda, &buf); + let mut t_maj = SampledSpectrum::new(1.0); + + while let Some(seg) = iter.next() { + if seg.sigma_maj[0] == 0. { + let dt = seg.t_max - seg.t_min; + let dt = if dt.is_infinite() { Float::MAX } else { dt }; + + t_maj *= (-dt * seg.sigma_maj).exp(); + continue; + } + + let mut t_min = seg.t_min; + loop { + let dist = -(1. - u.ln()) / seg.sigma_maj[0]; + let t = t_min + dist; + + if t < seg.t_max { + t_maj *= (-(t - t_min) * seg.sigma_maj).exp(); + let p = ray.at(t); + let mp = self.sample_point(p, lambda); + + if !callback(p, mp, seg.sigma_maj, t_maj) { + return SampledSpectrum::new(1.); + } + + t_maj = SampledSpectrum::new(1.); + t_min = t; + } else { + let dt = seg.t_max - t_min; + let dt = if dt.is_infinite() { Float::MAX } else { dt }; + + t_maj *= (-dt * seg.sigma_maj).exp(); + break; + } + } + } + + SampledSpectrum::new(1.) + } +} #[derive(Debug, Clone)] pub enum Medium { @@ -20,6 +406,357 @@ pub enum Medium { NanoVDB(NanoVDBMedium), } +#[derive(Debug, Clone)] +pub struct HomogeneousMedium { + sigma_a_spec: DenselySampledSpectrum, + sigma_s_spec: DenselySampledSpectrum, + le_spec: DenselySampledSpectrum, + phase: HGPhaseFunction, +} + +impl HomogeneousMedium { + pub fn new( + sigma_a: Spectrum, + sigma_s: Spectrum, + sigma_scale: Float, + le: Spectrum, + le_scale: Float, + g: Float, + ) -> Self { + let mut sigma_a_spec = + DenselySampledSpectrum::from_spectrum(&sigma_a, LAMBDA_MIN, LAMBDA_MAX); + let mut sigma_s_spec = + DenselySampledSpectrum::from_spectrum(&sigma_s, LAMBDA_MIN, LAMBDA_MAX); + let mut le_spec = DenselySampledSpectrum::from_spectrum(&le, LAMBDA_MIN, LAMBDA_MAX); + + sigma_a_spec.scale(sigma_scale); + sigma_s_spec.scale(sigma_scale); + le_spec.scale(le_scale); + + Self { + sigma_a_spec, + sigma_s_spec, + le_spec, + phase: HGPhaseFunction::new(g), + } + } +} + +impl MediumTrait for HomogeneousMedium { + fn is_emissive(&self) -> bool { + self.le_spec.max_value() > 0. + } + + fn sample_point(&self, _p: Point3f, lambda: &SampledWavelengths) -> MediumProperties { + let sigma_a = self.sigma_a_spec.sample(lambda); + let sigma_s = self.sigma_s_spec.sample(lambda); + let le = self.le_spec.sample(lambda); + MediumProperties { + sigma_a, + sigma_s, + phase: self.phase.into(), + le, + } + } + + fn sample_ray( + &self, + _ray: &Ray, + t_max: Float, + lambda: &SampledWavelengths, + _scratch: &Bump, + ) -> RayMajorantIterator { + let sigma_a = self.sigma_a_spec.sample(lambda); + let sigma_s = self.sigma_s_spec.sample(lambda); + let sigma_maj = sigma_a + sigma_s; + + let iter = HomogeneousMajorantIterator::new(0.0, t_max, sigma_maj); + + RayMajorantIterator::Homogeneous(iter) + } +} + +#[derive(Debug, Clone)] +pub struct GridMedium { + bounds: Bounds3f, + render_from_medium: Transform, + sigma_a_spec: DenselySampledSpectrum, + sigma_s_spec: DenselySampledSpectrum, + density_grid: SampledGrid, + phase: HGPhaseFunction, + temperature_grid: Option>, + le_spec: DenselySampledSpectrum, + le_scale: SampledGrid, + is_emissive: bool, + majorant_grid: MajorantGrid, +} + +impl GridMedium { + #[allow(clippy::too_many_arguments)] + pub fn new( + bounds: &Bounds3f, + render_from_medium: &Transform, + sigma_a: &Spectrum, + sigma_s: &Spectrum, + sigma_scale: Float, + g: Float, + density_grid: SampledGrid, + temperature_grid: Option>, + le: &Spectrum, + le_scale: SampledGrid, + ) -> Self { + let mut sigma_a_spec = + DenselySampledSpectrum::from_spectrum(sigma_a, LAMBDA_MIN, LAMBDA_MAX); + let mut sigma_s_spec = + DenselySampledSpectrum::from_spectrum(sigma_s, LAMBDA_MIN, LAMBDA_MAX); + let le_spec = DenselySampledSpectrum::from_spectrum(le, LAMBDA_MIN, LAMBDA_MAX); + sigma_a_spec.scale(sigma_scale); + sigma_s_spec.scale(sigma_scale); + + let mut majorant_grid = MajorantGrid::new(*bounds, Point3i::new(16, 16, 16)); + let is_emissive = if temperature_grid.is_some() { + true + } else { + le_spec.max_value() > 0. + }; + + for z in 0..majorant_grid.res.z() { + for y in 0..majorant_grid.res.y() { + for x in 0..majorant_grid.res.x() { + let bounds = majorant_grid.voxel_bounds(x, y, z); + majorant_grid.set(x, y, z, density_grid.max_value(bounds)); + } + } + } + + Self { + bounds: *bounds, + render_from_medium: *render_from_medium, + sigma_a_spec, + sigma_s_spec, + density_grid, + phase: HGPhaseFunction::new(g), + temperature_grid, + le_spec, + le_scale, + is_emissive, + majorant_grid, + } + } +} + +impl MediumTrait for GridMedium { + fn is_emissive(&self) -> bool { + self.is_emissive + } + + fn sample_point(&self, p: Point3f, lambda: &SampledWavelengths) -> MediumProperties { + let mut sigma_a = self.sigma_a_spec.sample(lambda); + let mut sigma_s = self.sigma_s_spec.sample(lambda); + let p_transform = self.render_from_medium.apply_inverse(p); + let p = Point3f::from(self.bounds.offset(&p_transform)); + let d = self.density_grid.lookup(p); + sigma_a *= d; + sigma_s *= d; + let scale = if self.is_emissive { + self.le_scale.lookup(p) + } else { + 0.0 + }; + + let le = if scale > 0.0 { + let raw_emission = match &self.temperature_grid { + Some(grid) => { + let temp = grid.lookup(p); + BlackbodySpectrum::new(temp).sample(lambda) + } + None => self.le_spec.sample(lambda), + }; + + raw_emission * scale + } else { + SampledSpectrum::new(0.0) + }; + + MediumProperties { + sigma_a, + sigma_s, + phase: self.phase.into(), + le, + } + } + + fn sample_ray( + &self, + ray: &Ray, + t_max: Float, + lambda: &SampledWavelengths, + _buf: &Bump, + ) -> RayMajorantIterator { + let (local_ray, local_t_max) = self.render_from_medium.apply_inverse_ray(ray, Some(t_max)); + + let inv_dir = Vector3f::new( + 1.0 / local_ray.d.x(), + 1.0 / local_ray.d.y(), + 1.0 / local_ray.d.z(), + ); + let dir_is_neg = [ + if local_ray.d.x() < 0.0 { 1 } else { 0 }, + if local_ray.d.y() < 0.0 { 1 } else { 0 }, + if local_ray.d.z() < 0.0 { 1 } else { 0 }, + ]; + + let (t_min, t_max) = + match self + .bounds + .intersect_p(local_ray.o, local_t_max, inv_dir, &dir_is_neg) + { + Some((t0, t1)) => (t0, t1), + None => return RayMajorantIterator::Void, // Missed the medium bounds + }; + + let sigma_t = self.sigma_a_spec.sample(lambda) + self.sigma_s_spec.sample(lambda); + + let iter = + DDAMajorantIterator::new(&local_ray, t_min, t_max, &self.majorant_grid, &sigma_t); + + RayMajorantIterator::DDA(iter) + } +} + +#[derive(Debug, Clone)] +pub struct RGBGridMedium { + bounds: Bounds3f, + render_from_medium: Transform, + le_grid: Option>, + le_scale: Float, + phase: HGPhaseFunction, + sigma_a_grid: Option>, + sigma_s_grid: Option>, + sigma_scale: Float, + majorant_grid: MajorantGrid, +} + +impl RGBGridMedium { + #[allow(clippy::too_many_arguments)] + pub fn new( + bounds: &Bounds3f, + render_from_medium: &Transform, + g: Float, + sigma_a_grid: Option>, + sigma_s_grid: Option>, + sigma_scale: Float, + le_grid: Option>, + le_scale: Float, + ) -> Self { + let mut majorant_grid = MajorantGrid::new(*bounds, Point3i::new(16, 16, 16)); + for z in 0..majorant_grid.res.x() { + for y in 0..majorant_grid.res.y() { + for x in 0..majorant_grid.res.x() { + let bounds = majorant_grid.voxel_bounds(x, y, z); + let convert = |s: &RGBUnboundedSpectrum| s.max_value(); + let max_sigma_t = sigma_a_grid + .as_ref() + .map_or(1.0, |g| g.max_value_convert(bounds, convert)) + + sigma_s_grid + .as_ref() + .map_or(1.0, |g| g.max_value_convert(bounds, convert)); + majorant_grid.set(x, y, z, sigma_scale * max_sigma_t); + } + } + } + + Self { + bounds: *bounds, + render_from_medium: *render_from_medium, + le_grid, + le_scale, + phase: HGPhaseFunction::new(g), + sigma_a_grid, + sigma_s_grid, + sigma_scale, + majorant_grid, + } + } +} + +impl MediumTrait for RGBGridMedium { + fn is_emissive(&self) -> bool { + self.le_grid.is_some() && self.le_scale > 0. + } + + fn sample_point(&self, p: Point3f, lambda: &SampledWavelengths) -> MediumProperties { + let p_transform = self.render_from_medium.apply_inverse(p); + let p = Point3f::from(self.bounds.offset(&p_transform)); + let convert = |s: &RGBUnboundedSpectrum| s.sample(lambda); + let sigma_a = self.sigma_scale + * self + .sigma_a_grid + .as_ref() + .map_or(SampledSpectrum::new(1.0), |g| g.lookup_convert(p, convert)); + + let sigma_s = self.sigma_scale + * self + .sigma_s_grid + .as_ref() + .map_or(SampledSpectrum::new(1.0), |g| g.lookup_convert(p, convert)); + + let le = self + .le_grid + .as_ref() + .filter(|_| self.le_scale > 0.0) + .map(|g| g.lookup_convert(p, |s| s.sample(lambda)) * self.le_scale) + .unwrap_or_else(|| SampledSpectrum::new(0.0)); + + MediumProperties { + sigma_a, + sigma_s, + phase: self.phase.into(), + le, + } + } + + fn sample_ray( + &self, + ray: &Ray, + t_max: Float, + _lambda: &SampledWavelengths, + _buf: &Bump, + ) -> RayMajorantIterator { + let (local_ray, local_t_max) = self.render_from_medium.apply_inverse_ray(ray, Some(t_max)); + + let inv_dir = Vector3f::new( + 1.0 / local_ray.d.x(), + 1.0 / local_ray.d.y(), + 1.0 / local_ray.d.z(), + ); + let dir_is_neg = [ + if local_ray.d.x() < 0.0 { 1 } else { 0 }, + if local_ray.d.y() < 0.0 { 1 } else { 0 }, + if local_ray.d.z() < 0.0 { 1 } else { 0 }, + ]; + let (t_min, t_max) = + match self + .bounds + .intersect_p(local_ray.o, local_t_max, inv_dir, &dir_is_neg) + { + Some((t0, t1)) => (t0, t1), + None => return RayMajorantIterator::Void, // Missed the medium bounds + }; + + let sigma_t = SampledSpectrum::new(1.); + + let iter = + DDAMajorantIterator::new(&local_ray, t_min, t_max, &self.majorant_grid, &sigma_t); + RayMajorantIterator::DDA(iter) + } +} + +#[derive(Debug, Clone)] +pub struct CloudMedium; +#[derive(Debug, Clone)] +pub struct NanoVDBMedium; + #[derive(Debug, Default, Clone)] pub struct MediumInterface { pub inside: Option>, diff --git a/src/core/mod.rs b/src/core/mod.rs index 67e631f..b8e0d9e 100644 --- a/src/core/mod.rs +++ b/src/core/mod.rs @@ -4,10 +4,10 @@ pub mod bxdf; pub mod cie; pub mod film; pub mod filter; -pub mod integrator; pub mod interaction; pub mod material; pub mod medium; +pub mod options; pub mod pbrt; pub mod primitive; pub mod sampler; diff --git a/src/core/options.rs b/src/core/options.rs index 5a320b2..b23ea8e 100644 --- a/src/core/options.rs +++ b/src/core/options.rs @@ -1,3 +1,4 @@ +use crate::geometry::Point2i; use std::ops::Deref; use std::sync::OnceLock; @@ -8,7 +9,7 @@ pub enum RenderingCoordinateSystem { World, } -#[derive(Debug, Clone, Copy)] +#[derive(Debug, Clone)] pub struct BasicPBRTOptions { pub seed: i32, pub quiet: bool, @@ -19,6 +20,10 @@ pub struct BasicPBRTOptions { pub use_gpu: bool, pub wavefront: bool, pub rendering_space: RenderingCoordinateSystem, + pub debug_start: Option<(Point2i, i32)>, + pub write_partial_images: bool, + pub mse_reference_image: Option, + pub mse_reference_output: Option, } impl Default for BasicPBRTOptions { @@ -33,6 +38,10 @@ impl Default for BasicPBRTOptions { use_gpu: false, wavefront: false, rendering_space: RenderingCoordinateSystem::CameraWorld, + debug_start: Some((Point2i::default(), 0)), + write_partial_images: false, + mse_reference_image: None, + mse_reference_output: None, } } } @@ -50,7 +59,7 @@ impl Default for PBRTOptions { fn default() -> Self { Self { basic: BasicPBRTOptions::default(), - n_threads: 0, // 0 usually implies "autodetect" + n_threads: 0, log_level: "info".to_string(), image_file: "output.exr".to_string(), } @@ -65,10 +74,6 @@ impl Deref for PBRTOptions { } } -// ----------------------------------------------------------------------- -// Global State Management -// ----------------------------------------------------------------------- - static OPTIONS: OnceLock = OnceLock::new(); pub fn init_pbrt(options: PBRTOptions) { diff --git a/src/core/pbrt.rs b/src/core/pbrt.rs index 9be4d18..2711de5 100644 --- a/src/core/pbrt.rs +++ b/src/core/pbrt.rs @@ -1,7 +1,10 @@ use crate::geometry::{Lerp, Point2f, Vector2f, Vector3f}; use num_traits::{Num, PrimInt}; +use std::collections::HashSet; +use std::hash::Hash; use std::ops::{Add, Mul}; use std::sync::atomic::{AtomicU64, Ordering as SyncOrdering}; +use std::sync::{Arc, Mutex}; pub type Float = f32; @@ -237,3 +240,30 @@ macro_rules! check_rare { } }; } + +pub struct InternCache { + cache: Mutex>>, +} + +impl InternCache +where + T: Eq + Hash + Clone, +{ + pub fn new() -> Self { + Self { + cache: Mutex::new(HashSet::new()), + } + } + + pub fn lookup(&self, value: T) -> Arc { + let mut lock = self.cache.lock().unwrap(); + + if let Some(existing) = lock.get(&value) { + return existing.clone(); // Returns a cheap Arc copy + } + + let new_item = Arc::new(value); + lock.insert(new_item.clone()); + new_item + } +} diff --git a/src/core/primitive.rs b/src/core/primitive.rs index 9007723..9b58594 100644 --- a/src/core/primitive.rs +++ b/src/core/primitive.rs @@ -1,17 +1,19 @@ use crate::core::aggregates::LinearBVHNode; -use crate::core::interaction::{Interaction, SurfaceInteraction}; +use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction}; use crate::core::material::Material; use crate::core::medium::{Medium, MediumInterface}; use crate::core::pbrt::Float; use crate::core::texture::{FloatTextureTrait, TextureEvalContext}; use crate::geometry::{Bounds3f, Ray}; use crate::lights::Light; -use crate::shapes::{Shape, ShapeIntersection}; +use crate::shapes::{Shape, ShapeIntersection, ShapeTrait}; use crate::utils::hash::hash_float; use crate::utils::transform::{AnimatedTransform, Transform}; +use enum_dispatch::enum_dispatch; use std::sync::Arc; +#[enum_dispatch] pub trait PrimitiveTrait: Send + Sync + std::fmt::Debug { fn bounds(&self) -> Bounds3f; fn intersect(&self, r: &Ray, t_max: Option) -> Option; @@ -23,7 +25,7 @@ pub struct GeometricPrimitive { shape: Arc, material: Arc, area_light: Arc, - medium_interface: Arc, + medium_interface: MediumInterface, alpha: Option>, } @@ -35,17 +37,17 @@ impl PrimitiveTrait for GeometricPrimitive { fn intersect(&self, r: &Ray, t_max: Option) -> Option { let mut si = self.shape.intersect(r, t_max)?; if let Some(ref alpha) = self.alpha { - let ctx: TextureEvalContext = (&*si.intr_mut()).into(); + let ctx = TextureEvalContext::from(&si.intr); let a = alpha.evaluate(&ctx); if a < 1.0 { let u = if a <= 0.0 { 1.0 } else { - hash_float((r.o, r.d)) + hash_float(&(r.o, r.d)) }; if u > a { - let r_next = si.intr().spawn_ray(r.d); + let r_next = si.intr.spawn_ray(r.d); let new_t_max = t_max.map(|t| t - si.t_hit()); @@ -59,11 +61,11 @@ impl PrimitiveTrait for GeometricPrimitive { } } - si.intr_mut().set_intersection_properties( + si.set_intersection_properties( self.material.clone(), self.area_light.clone(), Some(self.medium_interface.clone()), - r.medium.clone().expect("Medium not set"), + Some(r.medium.clone().expect("Medium not set")), ); Some(si) @@ -99,17 +101,17 @@ impl PrimitiveTrait for TransformedPrimitive { fn intersect(&self, r: &Ray, t_max: Option) -> Option { let (ray, t_max) = self.render_from_primitive.apply_inverse_ray(r, t_max); let mut si = self.primitive.intersect(&ray, Some(t_max))?; - let new_render: Box = self + let intr_wrapper = Interaction::Surface(si.intr); + let new_render_enum = self .render_from_primitive - .apply_to_interaction(si.intr.as_ref()); - let new_render_surf = new_render - .into_any() - .downcast::() - .expect( - "Failed to downcast Interaction to SurfaceInteraction. This should not happen.", - ); + .apply_to_interaction(&intr_wrapper); + + if let Interaction::Surface(new_surf) = new_render_enum { + si.intr = new_surf; + } else { + panic!("TransformedPrimitive: Transform changed interaction type (impossible)"); + } - si.intr = new_render_surf; Some(si) } @@ -134,16 +136,15 @@ impl PrimitiveTrait for AnimatedPrimitive { let interp_render_from_primitive = self.render_from_primitive.interpolate(r.time); let (ray, t_max) = interp_render_from_primitive.apply_inverse_ray(r, t_max); let mut si = self.primitive.intersect(&ray, Some(t_max))?; - let new_render: Box = - interp_render_from_primitive.apply_to_interaction(si.intr.as_ref()); - let new_render_surf = new_render - .into_any() - .downcast::() - .expect( - "Failed to downcast Interaction to SurfaceInteraction. This should not happen.", - ); + let wrapper = si.intr.into(); + let new_wrapper = interp_render_from_primitive.apply_to_interaction(&wrapper); + + if let Interaction::Surface(new_surf) = new_wrapper { + si.intr = new_surf; + } else { + panic!("TransformedPrimitive: Interaction type changed unexpectedly!"); + } - si.intr = new_render_surf; Some(si) } @@ -159,7 +160,7 @@ pub struct BVHAggregatePrimitive { nodes: Vec, } -impl BVHAggregatePrimitive { +impl PrimitiveTrait for BVHAggregatePrimitive { fn bounds(&self) -> Bounds3f { if !self.nodes.is_empty() { self.nodes[0].bounds @@ -182,12 +183,30 @@ impl BVHAggregatePrimitive { self.intersect_p(r, t_max) } } + +#[derive(Debug, Clone)] pub struct KdTreeAggregate; +impl PrimitiveTrait for KdTreeAggregate { + fn bounds(&self) -> Bounds3f { + todo!() + } + + fn intersect(&self, _r: &Ray, _t_max: Option) -> Option { + todo!() + } + + fn intersect_p(&self, _r: &Ray, _t_max: Option) -> bool { + todo!() + } +} + +#[derive(Clone, Debug)] +#[enum_dispatch(PrimitiveTrait)] pub enum Primitive { Geometric(GeometricPrimitive), Transformed(TransformedPrimitive), - Animated(Box), + Animated(AnimatedPrimitive), BVH(BVHAggregatePrimitive), KdTree(KdTreeAggregate), } diff --git a/src/core/sampler.rs b/src/core/sampler.rs index 8893254..1391578 100644 --- a/src/core/sampler.rs +++ b/src/core/sampler.rs @@ -745,9 +745,6 @@ pub trait SamplerTrait { fn get1d(&mut self) -> Float; fn get2d(&mut self) -> Point2f; fn get_pixel2d(&mut self) -> Point2f; - // fn clone_sampler(&self) -> Box { - // Box::new(self.clone()) - // } } #[enum_dispatch(SamplerTrait)] diff --git a/src/core/texture.rs b/src/core/texture.rs index 641d158..e19a763 100644 --- a/src/core/texture.rs +++ b/src/core/texture.rs @@ -1,4 +1,4 @@ -use crate::core::interaction::SurfaceInteraction; +use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction}; use crate::core::pbrt::Float; use crate::core::pbrt::{INV_2_PI, INV_PI, PI}; use crate::geometry::{Normal3f, Point2f, Point3f, Vector2f, Vector3f, VectorLike}; @@ -243,19 +243,18 @@ impl TextureMapping3DTrait for PointTransformMapping { } } -#[derive(Clone, Debug)] +#[derive(Clone, Default, Debug)] pub struct TextureEvalContext { - p: Point3f, - dpdx: Vector3f, - dpdy: Vector3f, - n: Normal3f, - uv: Point2f, - // All 0 - dudx: Float, - dudy: Float, - dvdx: Float, - dvdy: Float, - face_index: usize, + pub p: Point3f, + pub dpdx: Vector3f, + pub dpdy: Vector3f, + pub n: Normal3f, + pub uv: Point2f, + pub dudx: Float, + pub dudy: Float, + pub dvdx: Float, + pub dvdy: Float, + pub face_index: usize, } impl TextureEvalContext { @@ -290,7 +289,7 @@ impl TextureEvalContext { impl From<&SurfaceInteraction> for TextureEvalContext { fn from(si: &SurfaceInteraction) -> Self { Self { - p: si.common.pi.into(), + p: si.p(), dpdx: si.dpdx, dpdy: si.dpdy, n: si.common.n, @@ -304,6 +303,23 @@ impl From<&SurfaceInteraction> for TextureEvalContext { } } +impl From<&Interaction> for TextureEvalContext { + fn from(intr: &Interaction) -> Self { + match intr { + Interaction::Surface(si) => TextureEvalContext::from(si), + Interaction::Medium(mi) => TextureEvalContext { + p: mi.p(), + ..Default::default() + }, + + Interaction::Simple(si) => TextureEvalContext { + p: si.p(), + ..Default::default() + }, + } + } +} + #[enum_dispatch] pub trait FloatTextureTrait: Send + Sync + std::fmt::Debug { fn evaluate(&self, _ctx: &TextureEvalContext) -> Float { diff --git a/src/geometry/bounds.rs b/src/geometry/bounds.rs index 8e946e9..cecd97a 100644 --- a/src/geometry/bounds.rs +++ b/src/geometry/bounds.rs @@ -258,9 +258,7 @@ impl Bounds3f { ray_t_max: Float, inv_dir: Vector3f, dir_is_neg: &[usize; 3], - ) -> bool { - // We create a temporary array of references to access p_min/p_max by index - // The compiler optimizes this away entirely. + ) -> Option<(Float, Float)> { let bounds = [&self.p_min, &self.p_max]; // Check X @@ -271,11 +269,8 @@ impl Bounds3f { let ty_min = (bounds[dir_is_neg[1]].y() - o.y()) * inv_dir.y(); let ty_max = (bounds[1 - dir_is_neg[1]].y() - o.y()) * inv_dir.y(); - // Optional: PBRT Gamma correction for float precision (t_max *= 1.00000024) - // t_max *= 1.00000024; - if t_min > ty_max || ty_min > t_max { - return false; + return None; } if ty_min > t_min { t_min = ty_min; @@ -288,10 +283,8 @@ impl Bounds3f { let tz_min = (bounds[dir_is_neg[2]].z() - o.z()) * inv_dir.z(); let tz_max = (bounds[1 - dir_is_neg[2]].z() - o.z()) * inv_dir.z(); - // t_max *= 1.00000024; - if t_min > tz_max || tz_min > t_max { - return false; + return None; } if tz_min > t_min { t_min = tz_min; @@ -300,7 +293,11 @@ impl Bounds3f { t_max = tz_max; } - (t_min < ray_t_max) && (t_max > 0.0) + if (t_min < ray_t_max) && (t_max > 0.0) { + Some((t_min, t_max)) + } else { + None + } } pub fn intersect_with_inverse(&self, o: Point3f, d: Vector3f, ray_t_max: Float) -> bool { diff --git a/src/geometry/cone.rs b/src/geometry/cone.rs index 09e5c33..3fb84f3 100644 --- a/src/geometry/cone.rs +++ b/src/geometry/cone.rs @@ -4,8 +4,8 @@ use crate::utils::transform::Transform; #[derive(Debug, Clone)] pub struct DirectionCone { - w: Vector3f, - cos_theta: Float, + pub w: Vector3f, + pub cos_theta: Float, } impl Default for DirectionCone { @@ -60,57 +60,57 @@ impl DirectionCone { - wp.z() * (square(w.x() + square(w.y())))), ) } -} - -pub fn inside(d: &DirectionCone, w: Vector3f) -> bool { - !d.is_empty() && d.w.dot(w.normalize()) > d.cos_theta -} - -pub fn bound_subtended_directions(b: &Bounds3f, p: Point3f) -> DirectionCone { - let (p_center, radius) = b.bounding_sphere(); - if p.distance_squared(p_center) < square(radius) { - return DirectionCone::entire_sphere(); - } - - let w = (p_center - p).normalize(); - let sin2_theta_max = square(radius) / p_center.distance_squared(p); - let cos_theta_max = safe_sqrt(1. - sin2_theta_max); - DirectionCone::new(w, cos_theta_max) -} - -pub fn union(a: &DirectionCone, b: &DirectionCone) -> DirectionCone { - if a.is_empty() { - return b.clone(); - } - if b.is_empty() { - return a.clone(); - } - - // Handle the cases where one cone is inside the other - let theta_a = safe_acos(a.cos_theta); - let theta_b = safe_acos(b.cos_theta); - let theta_d = a.w.angle_between(b.w); - - if (theta_d + theta_b).min(PI) <= theta_b { - return a.clone(); - } - if (theta_d + theta_a).min(PI) <= theta_a { - return b.clone(); - } - - // Compute the spread angle of the merged cone, $\theta_o$ - let theta_o = (theta_a + theta_d + theta_b) / 2.; - if theta_o >= PI { - return DirectionCone::entire_sphere(); - } - - // Find the merged cone's axis and return cone union - let theta_r = theta_o - theta_a; - let wr = a.w.cross(b.w); - if wr.norm_squared() >= 0. { - return DirectionCone::entire_sphere(); - } - - let w = Transform::rotate_around_axis(degrees(theta_r), wr).apply_to_vector(a.w); - DirectionCone::new(w, theta_o.cos()) + + pub fn inside(d: &DirectionCone, w: Vector3f) -> bool { + !d.is_empty() && d.w.dot(w.normalize()) > d.cos_theta + } + + pub fn bound_subtended_directions(b: &Bounds3f, p: Point3f) -> DirectionCone { + let (p_center, radius) = b.bounding_sphere(); + if p.distance_squared(p_center) < square(radius) { + return DirectionCone::entire_sphere(); + } + + let w = (p_center - p).normalize(); + let sin2_theta_max = square(radius) / p_center.distance_squared(p); + let cos_theta_max = safe_sqrt(1. - sin2_theta_max); + DirectionCone::new(w, cos_theta_max) + } + + pub fn union(a: &DirectionCone, b: &DirectionCone) -> DirectionCone { + if a.is_empty() { + return b.clone(); + } + if b.is_empty() { + return a.clone(); + } + + // Handle the cases where one cone is inside the other + let theta_a = safe_acos(a.cos_theta); + let theta_b = safe_acos(b.cos_theta); + let theta_d = a.w.angle_between(b.w); + + if (theta_d + theta_b).min(PI) <= theta_b { + return a.clone(); + } + if (theta_d + theta_a).min(PI) <= theta_a { + return b.clone(); + } + + // Compute the spread angle of the merged cone, $\theta_o$ + let theta_o = (theta_a + theta_d + theta_b) / 2.; + if theta_o >= PI { + return DirectionCone::entire_sphere(); + } + + // Find the merged cone's axis and return cone union + let theta_r = theta_o - theta_a; + let wr = a.w.cross(b.w); + if wr.norm_squared() >= 0. { + return DirectionCone::entire_sphere(); + } + + let w = Transform::rotate_around_axis(degrees(theta_r), wr).apply_to_vector(a.w); + DirectionCone::new(w, theta_o.cos()) + } } diff --git a/src/geometry/primitives.rs b/src/geometry/primitives.rs index 0f130f1..e178b0b 100644 --- a/src/geometry/primitives.rs +++ b/src/geometry/primitives.rs @@ -1,5 +1,5 @@ use super::traits::{Sqrt, Tuple, VectorLike}; -use super::{Float, NumFloat, PI}; +use super::{Float, NumFloat, PI, clamp_t}; use crate::utils::interval::Interval; use crate::utils::math::{difference_of_products, quadratic, safe_asin}; use num_traits::{AsPrimitive, FloatConst, Num, Signed, Zero}; @@ -60,6 +60,50 @@ macro_rules! impl_tuple_core { pub fn floor(&self) -> $Struct { $Struct(self.0.map(|v| v.floor() as i32)) } + + #[inline] + pub fn average(&self) -> f32 { + let sum: f32 = self.0.iter().sum(); + sum / (N as f32) + } + } + + impl $Struct + where + T: Copy + PartialOrd, + { + #[inline] + pub fn min(&self, other: Self) -> Self { + let mut out = self.0; + for i in 0..N { + if other.0[i] < out[i] { + out[i] = other.0[i]; + } + } + Self(out) + } + + #[inline] + pub fn max(&self, other: Self) -> Self { + let mut out = self.0; + for i in 0..N { + if other.0[i] > out[i] { + out[i] = other.0[i] + } + } + Self(out) + } + + #[inline] + pub fn max_component_value(&self) -> T { + let mut m = self.0[0]; + for i in 1..N { + if self.0[i] > m { + m = self.0[i]; + } + } + m + } } impl $Struct @@ -695,6 +739,69 @@ where } } +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +#[repr(C)] +pub struct OctahedralVector { + x: u16, + y: u16, +} + +impl OctahedralVector { + pub fn new(mut v: Vector3f) -> Self { + v /= v.x().abs() + v.y().abs() + v.z().abs(); + + let (x_enc, y_enc) = if v.z() >= 0.0 { + (Self::encode(v.x()), Self::encode(v.y())) + } else { + ( + Self::encode((1.0 - v.y().abs()) * Self::sign(v.x())), + Self::encode((1.0 - v.x().abs()) * Self::sign(v.y())), + ) + }; + + Self { x: x_enc, y: y_enc } + } + + pub fn to_vector(self) -> Vector3f { + let mut v = Vector3f::default(); + + // Map [0, 65535] back to [-1, 1] + v[0] = -1.0 + 2.0 * (self.x as Float / 65535.0); + v[1] = -1.0 + 2.0 * (self.y as Float / 65535.0); + v[2] = 1.0 - (v.x().abs() + v.y().abs()); + + if v.z() < 0.0 { + let xo = v.x(); + v[0] = (1.0 - v.y().abs()) * Self::sign(xo); + v[1] = (1.0 - xo.abs()) * Self::sign(v.y()); + } + + v.normalize() + } + + #[inline] + pub fn sign(v: Float) -> Float { + 1.0.copysign(v) + } + + #[inline] + pub fn encode(f: Float) -> u16 { + (clamp_t((f + 1.0) / 2.0, 0.0, 1.0) * 65535.0).round() as u16 + } +} + +impl From for OctahedralVector { + fn from(v: Vector3f) -> Self { + Self::new(v) + } +} + +impl From for Vector3f { + fn from(ov: OctahedralVector) -> Self { + ov.to_vector() + } +} + #[derive(Copy, Clone, Debug, Default, PartialEq)] pub struct Frame { pub x: Vector3f, diff --git a/src/geometry/ray.rs b/src/geometry/ray.rs index 805b817..d6780c8 100644 --- a/src/geometry/ray.rs +++ b/src/geometry/ray.rs @@ -37,7 +37,7 @@ impl Ray { } } - pub fn evaluate(&self, t: Float) -> Point3f { + pub fn at(&self, t: Float) -> Point3f { self.o + self.d * t } diff --git a/src/image/io.rs b/src/image/io.rs index 2d50b11..cc9d9ff 100644 --- a/src/image/io.rs +++ b/src/image/io.rs @@ -110,15 +110,20 @@ impl Image { } fn write_exr(&self, path: &Path, _metadata: &ImageMetadata) -> Result<()> { - // EXR requires F32. If we have others, we rely on the float accessors or convert. + // EXR requires F32 let w = self.resolution.x() as usize; let h = self.resolution.y() as usize; let c = self.n_channels(); write_rgba_file(path, w, h, |x, y| { // Helper to get float value regardless of internal storage - let get = - |ch| self.get_channel(Point2i::new(x as i32, y as i32), ch, WrapMode::Clamp.into()); + let get = |ch| { + self.get_channel_with_wrap( + Point2i::new(x as i32, y as i32), + ch, + WrapMode::Clamp.into(), + ) + }; if c == 1 { let v = get(0); @@ -145,7 +150,6 @@ impl Image { // Header writeln!(writer, "PF")?; writeln!(writer, "{} {}", self.resolution.x(), self.resolution.y())?; - // Scale: Negative means Little Endian let scale = if cfg!(target_endian = "little") { -1.0 } else { @@ -153,11 +157,12 @@ impl Image { }; writeln!(writer, "{}", scale)?; - // PFM stores bottom-to-top. PBRT stores top-to-bottom. + // PBRT stores top-to-bottom. for y in (0..self.resolution.y()).rev() { for x in 0..self.resolution.x() { for c in 0..3 { - let val = self.get_channel(Point2i::new(x, y), c, WrapMode::Clamp.into()); + let val = + self.get_channel_with_wrap(Point2i::new(x, y), c, WrapMode::Clamp.into()); writer.write_all(&val.to_le_bytes())?; } } @@ -182,7 +187,6 @@ impl Image { } fn read_generic(path: &Path, encoding: Option) -> Result { - // 1. Load via 'image' crate let dyn_img = ImageReader::open(path) .with_context(|| format!("Failed to open image: {:?}", path))? .decode()?; @@ -191,7 +195,7 @@ fn read_generic(path: &Path, encoding: Option) -> Result Image { format: PixelFormat::F32, @@ -208,7 +212,7 @@ fn read_generic(path: &Path, encoding: Option) -> Result { - // Default to RGB8 for everything else (PNG, JPG) + // Default to RGB8 for everything else if dyn_img.color().has_alpha() { let buf = dyn_img.to_rgba8(); Image { @@ -314,9 +318,10 @@ fn read_pfm(path: &Path) -> Result { let mut pixels = vec![0.0 as Float; (w * h * channels) as usize]; - // PFM is Bottom-to-Top. We must flip Y while reading. + // PFM is Bottom-to-Top for y in 0..h { - let src_y = h - 1 - y; // Flip logic + // Flippety-do + let src_y = h - 1 - y; for x in 0..w { for c in 0..channels { let src_idx = ((src_y * w + x) * channels + c) as usize * 4; diff --git a/src/image/metadata.rs b/src/image/metadata.rs index 12c864a..07a005b 100644 --- a/src/image/metadata.rs +++ b/src/image/metadata.rs @@ -5,7 +5,55 @@ use crate::utils::math::SquareMatrix; use smallvec::SmallVec; use std::collections::HashMap; -pub type ImageChannelValues = SmallVec<[Float; 4]>; +use std::ops::{Deref, DerefMut}; + +#[derive(Clone, Debug, Default)] +pub struct ImageChannelValues(pub SmallVec<[Float; 4]>); + +impl ImageChannelValues { + pub fn average(&self) -> Float { + if self.0.is_empty() { + return 0.0; + } + let sum: Float = self.0.iter().sum(); + sum / (self.0.len() as Float) + } + + pub fn max_value(&self) -> Float { + self.0.iter().fold(Float::MIN, |a, &b| a.max(b)) + } +} + +impl From<&[Float]> for ImageChannelValues { + fn from(slice: &[Float]) -> Self { + Self(SmallVec::from_slice(slice)) + } +} + +impl From> for ImageChannelValues { + fn from(vec: Vec) -> Self { + Self(SmallVec::from_vec(vec)) + } +} + +impl From<[Float; N]> for ImageChannelValues { + fn from(arr: [Float; N]) -> Self { + Self(SmallVec::from_slice(&arr)) + } +} + +impl Deref for ImageChannelValues { + type Target = SmallVec<[Float; 4]>; + fn deref(&self) -> &Self::Target { + &self.0 + } +} + +impl DerefMut for ImageChannelValues { + fn deref_mut(&mut self) -> &mut Self::Target { + &mut self.0 + } +} #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] pub enum WrapMode { @@ -37,12 +85,22 @@ impl ImageChannelDesc { offset: offset.into(), } } - pub fn len(&self) -> usize { + + pub fn size(&self) -> usize { self.offset.len() } + pub fn is_empty(&self) -> bool { self.offset.is_empty() } + pub fn is_identity(&self) -> bool { + for i in 0..self.size() { + if self.offset[i] != i { + return false; + } + } + true + } } #[derive(Debug, Default)] diff --git a/src/image/mod.rs b/src/image/mod.rs index e2bb823..bafa6ad 100644 --- a/src/image/mod.rs +++ b/src/image/mod.rs @@ -4,10 +4,16 @@ pub mod ops; pub mod pixel; use crate::core::pbrt::{Float, lerp}; -use crate::geometry::{Point2f, Point2i}; -use crate::utils::color::{ColorEncoding, ColorEncodingTrait}; +use crate::geometry::{Bounds2f, Point2f, Point2fi, Point2i}; +use crate::utils::color::{ColorEncoding, ColorEncodingTrait, LINEAR}; +use crate::utils::containers::Array2D; +use crate::utils::math::square; +use core::hash; use half::f16; use pixel::PixelStorage; +use rayon::prelude::*; +use smallvec::{SmallVec, smallvec}; +use std::ops::{Deref, DerefMut}; pub use metadata::{ImageChannelDesc, ImageChannelValues, ImageMetadata, WrapMode, WrapMode2D}; @@ -73,27 +79,39 @@ pub struct ImageAndMetadata { } impl Image { - pub fn new( + fn from_vector( format: PixelFormat, resolution: Point2i, channel_names: Vec, encoding: ColorEncoding, ) -> Self { let size = (resolution.x() * resolution.y()) as usize * channel_names.len(); + let pixels = match format { PixelFormat::U8 => PixelData::U8(vec![0; size]), PixelFormat::F16 => PixelData::F16(vec![f16::ZERO; size]), PixelFormat::F32 => PixelData::F32(vec![0.0; size]), }; + Self { + format, resolution, channel_names, - format, - pixels, encoding, + pixels, } } + pub fn new( + format: PixelFormat, + resolution: Point2i, + channel_names: &[&str], + encoding: ColorEncoding, + ) -> Self { + let owned_names = channel_names.iter().map(|s| s.to_string()).collect(); + Self::from_vector(format, resolution, owned_names, encoding) + } + pub fn format(&self) -> PixelFormat { self.format } @@ -103,8 +121,14 @@ impl Image { pub fn n_channels(&self) -> usize { self.channel_names.len() } - pub fn channel_names(&self) -> &[String] { - &self.channel_names + pub fn channel_names(&self) -> Vec<&str> { + self.channel_names.iter().map(|s| s.as_str()).collect() + } + pub fn channel_names_from_desc(&self, desc: &ImageChannelDesc) -> Vec<&str> { + desc.offset + .iter() + .map(|&i| self.channel_names[i].as_str()) + .collect() } pub fn encoding(&self) -> ColorEncoding { self.encoding @@ -114,7 +138,11 @@ impl Image { (p.y() as usize * self.resolution.x() as usize + p.x() as usize) * self.n_channels() } - pub fn get_channel(&self, p: Point2i, c: usize, wrap: WrapMode2D) -> Float { + pub fn get_channel(&self, p: Point2i, c: usize) -> Float { + self.get_channel_with_wrap(p, c, WrapMode::Clamp.into()) + } + + pub fn get_channel_with_wrap(&self, p: Point2i, c: usize, wrap: WrapMode2D) -> Float { let mut pp = p; if !self.remap_pixel_coords(&mut pp, wrap) { return 0.0; @@ -128,12 +156,114 @@ impl Image { } } + pub fn get_channels(&self, p: Point2i, wrap: WrapMode2D) -> ImageChannelValues { + let mut pp = p; + + if !self.remap_pixel_coords(&mut pp, wrap) { + return ImageChannelValues(smallvec![0.0; self.n_channels()]); + } + + let start_idx = self.pixel_offset(pp); + let n_channels = self.n_channels(); + + let mut values: SmallVec<[Float; 4]> = SmallVec::with_capacity(n_channels); + + match &self.pixels { + PixelData::U8(data) => { + let slice = &data[start_idx..start_idx + n_channels]; + for &v in slice { + values.push(u8::to_linear(v, self.encoding)); + } + } + PixelData::F16(data) => { + let slice = &data[start_idx..start_idx + n_channels]; + for &v in slice { + values.push(f16::to_linear(v, self.encoding)); + } + } + PixelData::F32(data) => { + let slice = &data[start_idx..start_idx + n_channels]; + for &v in slice { + values.push(f32::to_linear(v, self.encoding)); + } + } + } + + ImageChannelValues(values) + } + + pub fn get_channels_desc( + &self, + p: Point2i, + desc: &ImageChannelDesc, + wrap: WrapMode2D, + ) -> ImageChannelValues { + let mut pp = p; + if !self.remap_pixel_coords(&mut pp, wrap) { + return ImageChannelValues(smallvec![0.0; desc.offset.len()]); + } + + let pixel_offset = self.pixel_offset(pp); + + let mut values: SmallVec<[Float; 4]> = SmallVec::with_capacity(desc.offset.len()); + + match &self.pixels { + PixelData::U8(data) => { + for &channel_idx in &desc.offset { + let val = data[pixel_offset + channel_idx]; + values.push(u8::to_linear(val, self.encoding)); + } + } + PixelData::F16(data) => { + for &channel_idx in &desc.offset { + let val = data[pixel_offset + channel_idx]; + values.push(f16::to_linear(val, self.encoding)); + } + } + PixelData::F32(data) => { + for &channel_idx in &desc.offset { + let val = data[pixel_offset + channel_idx]; + values.push(f32::to_linear(val, self.encoding)); + } + } + } + + ImageChannelValues(values) + } + + pub fn get_channels_default(&self, p: Point2i) -> ImageChannelValues { + self.get_channels(p, WrapMode::Clamp.into()) + } + pub fn all_channels_desc(&self) -> ImageChannelDesc { ImageChannelDesc { offset: (0..self.n_channels()).collect(), } } + pub fn get_channel_desc( + &self, + requested_channels: &[&str], + ) -> Result { + let mut offset = Vec::with_capacity(requested_channels.len()); + + for &req in requested_channels.iter() { + match self.channel_names.iter().position(|n| n == req) { + Some(idx) => { + offset.push(idx); + } + None => { + return Err(format!( + "Image is missing requested channel '{}'. Available channels: {:?}", + req, self.channel_names + )); + } + } + } + + Ok(ImageChannelDesc { offset }) + } + pub fn set_channel(&mut self, p: Point2i, c: usize, value: Float) { let val_no_nan = if value.is_nan() { 0.0 } else { value }; let offset = self.pixel_offset(p) + c; @@ -154,8 +284,8 @@ impl Image { desc: &ImageChannelDesc, values: &ImageChannelValues, ) { - assert_eq!(desc.len(), values.len()); - for i in 0..desc.len() { + assert_eq!(desc.size(), values.len()); + for i in 0..desc.size() { self.set_channel(p, desc.offset[i], values[i]); } } @@ -181,17 +311,134 @@ impl Image { true } - pub fn bilerp_channel(&self, p: Point2f, c: usize, wrap_mode: WrapMode2D) -> Float { + pub fn bilerp_channel(&self, p: Point2f, c: usize) -> Float { + self.bilerp_channel_with_wrap(p, c, WrapMode::Clamp.into()) + } + + pub fn bilerp_channel_with_wrap(&self, p: Point2f, c: usize, wrap_mode: WrapMode2D) -> Float { let x = p.x() * self.resolution.x() as Float - 0.5; let y = p.y() * self.resolution.y() as Float - 0.5; let xi = x.floor() as i32; let yi = y.floor() as i32; let dx = x - xi as Float; let dy = y - yi as Float; - let v00 = self.get_channel(Point2i::new(xi, yi), c, wrap_mode); - let v10 = self.get_channel(Point2i::new(xi + 1, yi), c, wrap_mode); - let v01 = self.get_channel(Point2i::new(xi, yi + 1), c, wrap_mode); - let v11 = self.get_channel(Point2i::new(xi + 1, yi + 1), c, wrap_mode); + let v00 = self.get_channel_with_wrap(Point2i::new(xi, yi), c, wrap_mode); + let v10 = self.get_channel_with_wrap(Point2i::new(xi + 1, yi), c, wrap_mode); + let v01 = self.get_channel_with_wrap(Point2i::new(xi, yi + 1), c, wrap_mode); + let v11 = self.get_channel_with_wrap(Point2i::new(xi + 1, yi + 1), c, wrap_mode); lerp(dy, lerp(dx, v00, v10), lerp(dx, v01, v11)) } + + pub fn lookup_nearest_channel_with_wrap( + &self, + p: Point2f, + c: usize, + wrap_mode: WrapMode2D, + ) -> Float { + let pi = Point2i::new( + p.x() as i32 * self.resolution.x(), + p.y() as i32 * self.resolution.y(), + ); + + self.get_channel_with_wrap(pi, c, wrap_mode) + } + + pub fn lookup_nearest_channel(&self, p: Point2f, c: usize) -> Float { + self.lookup_nearest_channel_with_wrap(p, c, WrapMode::Clamp.into()) + } + + pub fn get_sampling_distribution(&self, dxd_a: F, domain: Bounds2f) -> Array2D + where + F: Fn(Point2f) -> Float + Sync + Send, + { + let width = self.resolution.x(); + let height = self.resolution.y(); + + let mut dist = Array2D::new_with_dims(width as usize, height as usize); + + dist.values + .par_chunks_mut(width as usize) + .enumerate() + .for_each(|(y, row)| { + let y = y as i32; + + for (x, out_val) in row.iter_mut().enumerate() { + let x = x as i32; + + let value = self.get_channels_default(Point2i::new(x, y)).average(); + + let u = (x as Float + 0.5) / width as Float; + let v = (y as Float + 0.5) / height as Float; + let p = domain.lerp(Point2f::new(u, v)); + *out_val = value * dxd_a(p); + } + }); + + dist + } + + pub fn get_sampling_distribution_uniform(&self) -> Array2D { + let default_domain = Bounds2f::from_points(Point2f::new(0.0, 0.0), Point2f::new(1.0, 1.0)); + + self.get_sampling_distribution(|_| 1.0, default_domain) + } + + pub fn mse( + &self, + desc: ImageChannelDesc, + ref_img: &Image, + generate_mse_image: bool, + ) -> (ImageChannelValues, Option) { + let mut sum_se: Vec = vec![0.; desc.size()]; + let names_ref = self.channel_names_from_desc(&desc); + let ref_desc = ref_img + .get_channel_desc(&self.channel_names_from_desc(&desc)) + .expect("Channels not found in image"); + assert_eq!(self.resolution(), ref_img.resolution()); + + let width = self.resolution.x() as usize; + let height = self.resolution.y() as usize; + let n_channels = desc.offset.len(); + let mut mse_pixels = if generate_mse_image { + vec![0.0f32; width * height * n_channels] + } else { + Vec::new() + }; + + for y in 0..self.resolution().y() { + for x in 0..self.resolution().x() { + let v = self.get_channels_desc(Point2i::new(x, y), &desc, WrapMode::Clamp.into()); + let v_ref = + self.get_channels_desc(Point2i::new(x, y), &ref_desc, WrapMode::Clamp.into()); + for c in 0..desc.size() { + let se = square(v[c] as f64 - v_ref[c] as f64); + if se.is_infinite() { + continue; + } + sum_se[c] += se; + if generate_mse_image { + let idx = (y as usize * width + x as usize) * n_channels + c; + mse_pixels[idx] = se as f32; + } + } + } + } + + let pixel_count = (self.resolution.x() * self.resolution.y()) as f64; + let mse_values: SmallVec<[Float; 4]> = + sum_se.iter().map(|&s| (s / pixel_count) as Float).collect(); + + let mse_image = if generate_mse_image { + Some(Image::new( + PixelFormat::F32, + self.resolution, + &names_ref, + LINEAR, + )) + } else { + None + }; + + (ImageChannelValues(mse_values), mse_image) + } } diff --git a/src/image/ops.rs b/src/image/ops.rs index 6440e02..cff0793 100644 --- a/src/image/ops.rs +++ b/src/image/ops.rs @@ -31,7 +31,7 @@ impl Image { bounds.p_max.y() - bounds.p_min.y(), ); - let mut new_image = Image::new( + let mut new_image = Image::from_vector( self.format, new_res, self.channel_names.clone(), @@ -87,24 +87,20 @@ impl Image { "ResizeUp requires Float format" ); - // Create destination image wrapped in Mutex for tile-based write access - let resampled_image = Arc::new(Mutex::new(Image::new( + let resampled_image = Arc::new(Mutex::new(Image::from_vector( PixelFormat::F32, // Force float output new_res, self.channel_names.clone(), self.encoding, ))); - // Precompute weights let x_weights = resample_weights(self.resolution.x() as usize, new_res.x() as usize); let y_weights = resample_weights(self.resolution.y() as usize, new_res.y() as usize); let n_channels = self.n_channels(); - // Generate Tiles let tile_size = 16; let tiles = generate_tiles(new_res, tile_size); - // Parallel Execution tiles.par_iter().for_each(|out_extent| { let x_start = x_weights[out_extent.p_min.x() as usize].first_pixel; let x_end = x_weights[(out_extent.p_max.x() - 1) as usize].first_pixel + 4; @@ -150,7 +146,7 @@ impl Image { } let new_res = Point2i::new((old.x() / 2).max(1), (old.y() / 2).max(1)); - let mut next = Image::new( + let mut next = Image::from_vector( prev.format, new_res, prev.channel_names.clone(), @@ -235,7 +231,7 @@ fn copy_rect_out_kernel( // We fall back to get_channel which handles the wrapping math. let p = Point2i::new(x, y); for c in 0..channels { - row_buf[x_rel * channels + c] = image.get_channel(p, c, wrap); + row_buf[x_rel * channels + c] = image.get_channel_with_wrap(p, c, wrap); } } } @@ -301,7 +297,7 @@ fn downsample_kernel( let sx = src_x as i32 + dx; let sy = src_y as i32 + dy; if sx < old_res.x() && sy < old_res.y() { - sum += prev.get_channel(Point2i::new(sx, sy), c, wrap); + sum += prev.get_channel_with_wrap(Point2i::new(sx, sy), c, wrap); count += 1.0; } } diff --git a/src/lib.rs b/src/lib.rs index 682736b..127924c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -6,6 +6,7 @@ mod camera; mod core; mod geometry; mod image; +mod integrators; mod lights; mod shapes; mod spectra; diff --git a/src/shapes/bilinear.rs b/src/shapes/bilinear.rs index d903492..3d18b95 100644 --- a/src/shapes/bilinear.rs +++ b/src/shapes/bilinear.rs @@ -1,8 +1,9 @@ use super::{ BilinearIntersection, BilinearPatchShape, Bounds3f, DirectionCone, Interaction, Normal3f, Point2f, Point3f, Point3fi, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, - SurfaceInteraction, Vector3f, + ShapeTrait, SurfaceInteraction, Vector3f, }; +use crate::core::interaction::InteractionTrait; use crate::core::pbrt::{Float, clamp_t, gamma, lerp}; use crate::geometry::{Tuple, VectorLike, spherical_quad_area}; use crate::utils::math::{SquareMatrix, difference_of_products, quadratic}; @@ -501,12 +502,14 @@ impl BilinearPatchShape { } n } +} - pub fn area(&self) -> Float { +impl ShapeTrait for BilinearPatchShape { + fn area(&self) -> Float { self.area } - pub fn normal_bounds(&self) -> DirectionCone { + fn normal_bounds(&self) -> DirectionCone { let data = self.get_data(); if data.p00 == data.p10 || data.p10 == data.p11 @@ -553,19 +556,19 @@ impl BilinearPatchShape { DirectionCone::new(n_avg.into(), clamp_t(cos_theta, -1., 1.)) } - pub fn bounds(&self) -> Bounds3f { + fn bounds(&self) -> Bounds3f { let data = self.get_data(); Bounds3f::from_points(data.p00, data.p01).union(Bounds3f::from_points(data.p10, data.p11)) } - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + fn intersect(&self, ray: &Ray, t_max: Option) -> Option { let t_max_val = t_max?; let data = self.get_data(); if let Some(bilinear_hit) = self.intersect_bilinear_patch(ray, t_max_val, &data) { let intr = self.interaction_from_intersection(&data, bilinear_hit.uv, ray.time, -ray.d); Some(ShapeIntersection { - intr: Box::new(intr), + intr, t_hit: bilinear_hit.t, }) } else { @@ -573,14 +576,14 @@ impl BilinearPatchShape { } } - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { let t_max_val = t_max.unwrap_or(Float::INFINITY); let data = self.get_data(); self.intersect_bilinear_patch(ray, t_max_val, &data) .is_some() } - pub fn sample(&self, u: Point2f) -> Option { + fn sample(&self, u: Point2f) -> Option { let data = self.get_data(); // Sample bilinear patch parametric coordinate (u, v) let (uv, pdf) = self.sample_parametric_coords(&data, u); @@ -612,7 +615,7 @@ impl BilinearPatchShape { }) } - pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { let data = self.get_data(); let v00 = (data.p00 - ctx.p()).normalize(); let v10 = (data.p10 - ctx.p()).normalize(); @@ -629,10 +632,11 @@ impl BilinearPatchShape { } } - pub fn pdf(&self, intr: &dyn Interaction) -> Float { - let Some(si) = intr.as_any().downcast_ref::() else { - return 0.; + fn pdf(&self, intr: &Interaction) -> Float { + let Interaction::Surface(si) = intr else { + return 0.0; }; + let data = self.get_data(); let uv = if let Some(uvs) = &data.mesh.uv { Point2f::invert_bilinear(si.uv, uvs) @@ -659,7 +663,7 @@ impl BilinearPatchShape { if cross == 0. { 0. } else { param_pdf / cross } } - pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { + fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { let ray = ctx.spawn_ray(wi); let Some(isect) = self.intersect(&ray, None) else { return 0.; @@ -675,8 +679,10 @@ impl BilinearPatchShape { let use_area_sampling = !self.rectangle || data.mesh.image_distribution.is_some() || spherical_quad_area(v00, v10, v01, v11) <= Self::MIN_SPHERICAL_SAMPLE_AREA; + if use_area_sampling { - let isect_pdf = self.pdf(isect.intr.as_ref()); + let intr_wrapper = Interaction::Surface(isect.intr.clone()); + let isect_pdf = self.pdf(&intr_wrapper); let distsq = ctx.p().distance_squared(isect.intr.p()); let absdot = Vector3f::from(isect.intr.n()).abs_dot(-wi); if absdot == 0. { diff --git a/src/shapes/curves.rs b/src/shapes/curves.rs index 60eed9d..6e6a585 100644 --- a/src/shapes/curves.rs +++ b/src/shapes/curves.rs @@ -1,3 +1,4 @@ +use crate::core::interaction::InteractionTrait; use crate::core::pbrt::{clamp_t, lerp}; use crate::utils::math::square; use crate::utils::splines::{ @@ -8,7 +9,7 @@ use crate::utils::transform::look_at; use super::{ Bounds3f, CurveCommon, CurveShape, CurveType, DirectionCone, Float, Interaction, Normal3f, Point2f, Point3f, Point3fi, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, - SurfaceInteraction, Transform, Vector2f, Vector3f, VectorLike, + ShapeTrait, SurfaceInteraction, Transform, Vector2f, Vector3f, VectorLike, }; use std::sync::Arc; @@ -245,7 +246,7 @@ impl CurveShape { let p_error = Vector3f::fill(hit_width); let flip_normal = self.common.reverse_orientation ^ self.common.transform_swap_handedness; - let pi = Point3fi::new_with_error(context.ray.evaluate(t_hit), p_error); + let pi = Point3fi::new_with_error(context.ray.at(t_hit), p_error); let intr = SurfaceInteraction::new( pi, Point2f::new(u, v), @@ -258,13 +259,12 @@ impl CurveShape { flip_normal, ); - Some(ShapeIntersection { - intr: Box::new(intr), - t_hit, - }) + Some(ShapeIntersection { intr, t_hit }) } +} - pub fn bounds(&self) -> Bounds3f { +impl ShapeTrait for CurveShape { + fn bounds(&self) -> Bounds3f { let cs_span = self.common.cp_obj; let obj_bounds = bound_cubic_bezier(&cs_span, self.u_min, self.u_max); let width0 = lerp(self.u_min, self.common.width[0], self.common.width[1]); @@ -275,11 +275,11 @@ impl CurveShape { .apply_to_bounds(obj_bounds_expand) } - pub fn normal_bounds(&self) -> DirectionCone { + fn normal_bounds(&self) -> DirectionCone { DirectionCone::entire_sphere() } - pub fn area(&self) -> Float { + fn area(&self) -> Float { let cp_obj = cubic_bezier_control_points(&self.common.cp_obj, self.u_min, self.u_max); let width0 = lerp(self.u_min, self.common.width[0], self.common.width[1]); let width1 = lerp(self.u_max, self.common.width[0], self.common.width[1]); @@ -291,32 +291,28 @@ impl CurveShape { approx_length * avg_width } - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY)) .is_some() } - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + fn intersect(&self, ray: &Ray, t_max: Option) -> Option { self.intersect_ray(ray, t_max.unwrap_or(Float::INFINITY)) } - pub fn pdf(&self, _interaction: &dyn Interaction) -> Float { + fn pdf(&self, _interaction: &Interaction) -> Float { todo!() } - pub fn pdf_from_context(&self, _ctx: &ShapeSampleContext, _wi: Vector3f) -> Float { + fn pdf_from_context(&self, _ctx: &ShapeSampleContext, _wi: Vector3f) -> Float { todo!() } - pub fn sample(&self, _u: Point2f) -> Option { + fn sample(&self, _u: Point2f) -> Option { todo!() } - pub fn sample_from_context( - &self, - _ctx: &ShapeSampleContext, - _u: Point2f, - ) -> Option { + fn sample_from_context(&self, _ctx: &ShapeSampleContext, _u: Point2f) -> Option { todo!() } } diff --git a/src/shapes/cylinder.rs b/src/shapes/cylinder.rs index 4f05673..795fa64 100644 --- a/src/shapes/cylinder.rs +++ b/src/shapes/cylinder.rs @@ -1,8 +1,9 @@ use super::{ Bounds3f, CylinderShape, DirectionCone, Float, Interaction, Normal3f, PI, Point2f, Point3f, Point3fi, QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, - SurfaceInteraction, Transform, Vector3f, Vector3fi, + ShapeTrait, SurfaceInteraction, Transform, Vector3f, Vector3fi, }; +use crate::core::interaction::InteractionTrait; use crate::core::pbrt::{gamma, lerp}; use crate::geometry::{Sqrt, Tuple, VectorLike}; use crate::utils::interval::Interval; @@ -166,12 +167,14 @@ impl CylinderShape { flip_normal, ) } +} - pub fn area(&self) -> Float { +impl ShapeTrait for CylinderShape { + fn area(&self) -> Float { (self.z_max - self.z_min) * self.radius * self.phi_max } - pub fn bounds(&self) -> Bounds3f { + fn bounds(&self) -> Bounds3f { self.render_from_object .apply_to_bounds(Bounds3f::from_points( Point3f::new(-self.radius, -self.radius, self.z_min), @@ -179,11 +182,11 @@ impl CylinderShape { )) } - pub fn normal_bounds(&self) -> DirectionCone { + fn normal_bounds(&self) -> DirectionCone { DirectionCone::entire_sphere() } - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + fn intersect(&self, ray: &Ray, t_max: Option) -> Option { let t = t_max.unwrap_or(Float::INFINITY); if let Some(isect) = self.basic_intersect(ray, t) { let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time); @@ -193,7 +196,7 @@ impl CylinderShape { } } - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { if let Some(t) = t_max { self.basic_intersect(ray, t).is_some() } else { @@ -201,11 +204,11 @@ impl CylinderShape { } } - pub fn pdf(&self, _interaction: &dyn Interaction) -> Float { + fn pdf(&self, _interaction: &Interaction) -> Float { 1. / self.area() } - pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { + fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { let ray = ctx.spawn_ray(wi); if let Some(isect) = self.intersect(&ray, None) { let n = isect.intr.n(); @@ -220,7 +223,7 @@ impl CylinderShape { } } - pub fn sample(&self, u: Point2f) -> Option { + fn sample(&self, u: Point2f) -> Option { let z = lerp(u[0], self.z_min, self.z_max); let phi = u[1] * self.phi_max; let mut p_obj = Point3f::new(self.radius * phi.cos(), self.radius * phi.sin(), z); @@ -249,7 +252,7 @@ impl CylinderShape { }) } - pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { let mut ss = self.sample(u)?; let intr = Arc::make_mut(&mut ss.intr); intr.get_common_mut().time = ctx.time; diff --git a/src/shapes/disk.rs b/src/shapes/disk.rs index aab20a0..8cc51ea 100644 --- a/src/shapes/disk.rs +++ b/src/shapes/disk.rs @@ -1,8 +1,9 @@ use super::{ Bounds3f, DirectionCone, DiskShape, Float, Interaction, Normal3f, PI, Point2f, Point3f, Point3fi, QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, - SurfaceInteraction, Transform, Vector3f, + ShapeTrait, SurfaceInteraction, Transform, Vector3f, }; +use crate::core::interaction::InteractionTrait; use crate::geometry::VectorLike; use crate::utils::math::square; use crate::utils::sampling::sample_uniform_disk_concentric; @@ -100,19 +101,21 @@ impl DiskShape { flip_normal, ) } +} - pub fn area(&self) -> Float { +impl ShapeTrait for DiskShape { + fn area(&self) -> Float { self.phi_max * 0.5 * (square(self.radius) - square(self.inner_radius)) } - pub fn bounds(&self) -> Bounds3f { + fn bounds(&self) -> Bounds3f { self.render_from_object .apply_to_bounds(Bounds3f::from_points( Point3f::new(-self.radius, -self.radius, self.height), Point3f::new(self.radius, self.radius, self.height), )) } - pub fn normal_bounds(&self) -> DirectionCone { + fn normal_bounds(&self) -> DirectionCone { let mut n = self .render_from_object .apply_to_normal(Normal3f::new(0., 0., 1.)); @@ -122,7 +125,7 @@ impl DiskShape { DirectionCone::new_from_vector(Vector3f::from(n)) } - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + fn intersect(&self, ray: &Ray, t_max: Option) -> Option { let t = t_max.unwrap_or(Float::INFINITY); if let Some(isect) = self.basic_intersect(ray, t) { let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time); @@ -132,7 +135,7 @@ impl DiskShape { } } - pub fn sample(&self, u: Point2f) -> Option { + fn sample(&self, u: Point2f) -> Option { let pd = sample_uniform_disk_concentric(u); let p_obj = Point3f::new(pd.x() * self.radius, pd.y() * self.radius, self.height); let pi = self @@ -160,7 +163,7 @@ impl DiskShape { }) } - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { if let Some(t) = t_max { self.basic_intersect(ray, t).is_some() } else { @@ -168,7 +171,7 @@ impl DiskShape { } } - pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { let mut ss = self.sample(u)?; let intr = Arc::make_mut(&mut ss.intr); intr.get_common_mut().time = ctx.time; @@ -184,11 +187,11 @@ impl DiskShape { Some(ss) } - pub fn pdf(&self, _interaction: &dyn Interaction) -> Float { + fn pdf(&self, _interaction: &Interaction) -> Float { 1. / self.area() } - pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { + fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { let ray = ctx.spawn_ray(wi); if let Some(isect) = self.intersect(&ray, None) { let n = isect.intr.n(); diff --git a/src/shapes/mod.rs b/src/shapes/mod.rs index d4a6572..0edb614 100644 --- a/src/shapes/mod.rs +++ b/src/shapes/mod.rs @@ -5,7 +5,9 @@ pub mod disk; pub mod sphere; pub mod triangle; -use crate::core::interaction::{Interaction, MediumInteraction, SurfaceInteraction}; +use crate::core::interaction::{ + Interaction, InteractionTrait, MediumInteraction, SurfaceInteraction, +}; use crate::core::material::Material; use crate::core::medium::{Medium, MediumInterface}; use crate::core::pbrt::{Float, PI}; @@ -16,6 +18,7 @@ use crate::geometry::{ use crate::lights::Light; use crate::utils::math::{next_float_down, next_float_up}; use crate::utils::transform::Transform; +use enum_dispatch::enum_dispatch; use std::sync::{Arc, Mutex}; #[derive(Debug, Clone)] @@ -32,6 +35,20 @@ pub struct SphereShape { transform_swap_handedness: bool, } +impl Default for SphereShape { + fn default() -> Self { + Self::new( + Transform::default().into(), + Transform::default().into(), + false, + 1.0, + -1.0, + 1.0, + 360.0, + ) + } +} + #[derive(Debug, Clone)] pub struct CylinderShape { radius: Float, @@ -143,28 +160,13 @@ pub struct CurveShape { // Define Intersection objects. This only varies for #[derive(Debug, Clone)] pub struct ShapeIntersection { - pub intr: Box, + pub intr: SurfaceInteraction, pub t_hit: Float, } impl ShapeIntersection { pub fn new(intr: SurfaceInteraction, t_hit: Float) -> Self { - Self { - intr: Box::new(intr), - t_hit, - } - } - - pub fn as_ref(&self) -> &dyn Interaction { - &*self.intr - } - - pub fn intr(&self) -> &dyn Interaction { - &*self.intr - } - - pub fn intr_mut(&mut self) -> &mut SurfaceInteraction { - &mut self.intr + Self { intr, t_hit } } pub fn t_hit(&self) -> Float { @@ -179,7 +181,7 @@ impl ShapeIntersection { &mut self, mtl: Arc, area: Arc, - prim_medium_interface: Option>, + prim_medium_interface: Option, ray_medium: Option>, ) { let ray_medium = ray_medium.expect("Ray medium must be defined for intersection"); @@ -229,16 +231,16 @@ impl BilinearIntersection { #[derive(Clone)] pub struct ShapeSample { - intr: Arc, - pdf: Float, + pub intr: Arc, + pub pdf: Float, } #[derive(Clone, Debug)] pub struct ShapeSampleContext { - pi: Point3fi, - n: Normal3f, - ns: Normal3f, - time: Float, + pub pi: Point3fi, + pub n: Normal3f, + pub ns: Normal3f, + pub time: Float, } impl ShapeSampleContext { @@ -286,10 +288,22 @@ impl ShapeSampleContext { } } -#[derive(Default, Debug, Clone)] +#[enum_dispatch] +pub trait ShapeTrait { + fn bounds(&self) -> Bounds3f; + fn normal_bounds(&self) -> DirectionCone; + fn intersect(&self, ray: &Ray, t_max: Option) -> Option; + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool; + fn area(&self) -> Float; + fn pdf(&self, interaction: &Interaction) -> Float; + fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float; + fn sample(&self, u: Point2f) -> Option; + fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option; +} + +#[derive(Debug, Clone)] +#[enum_dispatch(ShapeTrait)] pub enum Shape { - #[default] - None, Sphere(SphereShape), Cylinder(CylinderShape), Disk(DiskShape), @@ -298,112 +312,8 @@ pub enum Shape { Curve(CurveShape), } -impl Shape { - pub fn bounds(&self) -> Bounds3f { - match self { - Shape::None => Bounds3f::default(), - Shape::Sphere(s) => s.bounds(), - Shape::Cylinder(c) => c.bounds(), - Shape::Disk(d) => d.bounds(), - Shape::Triangle(t) => t.bounds(), - Shape::BilinearPatch(b) => b.bounds(), - Shape::Curve(c) => c.bounds(), - } - } - - pub fn normal_bounds(&self) -> DirectionCone { - match self { - Shape::None => DirectionCone::entire_sphere(), - Shape::Sphere(s) => s.normal_bounds(), - Shape::Cylinder(c) => c.normal_bounds(), - Shape::Disk(d) => d.normal_bounds(), - Shape::Triangle(t) => t.normal_bounds(), - Shape::BilinearPatch(b) => b.normal_bounds(), - Shape::Curve(c) => c.normal_bounds(), - } - } - - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { - match self { - Shape::None => None, - Shape::Sphere(s) => s.intersect(ray, t_max), - Shape::Cylinder(c) => c.intersect(ray, t_max), - Shape::Disk(d) => d.intersect(ray, t_max), - Shape::Triangle(t) => t.intersect(ray, t_max), - Shape::BilinearPatch(b) => b.intersect(ray, t_max), - Shape::Curve(c) => c.intersect(ray, t_max), - } - } - - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { - match self { - Shape::None => false, - Shape::Sphere(s) => s.intersect_p(ray, t_max), - Shape::Cylinder(c) => c.intersect_p(ray, t_max), - Shape::Disk(d) => d.intersect_p(ray, t_max), - Shape::Triangle(t) => t.intersect_p(ray, t_max), - Shape::BilinearPatch(b) => b.intersect_p(ray, t_max), - Shape::Curve(c) => c.intersect_p(ray, t_max), - } - } - - pub fn area(&self) -> Float { - match self { - Shape::None => 0., - Shape::Sphere(s) => s.area(), - Shape::Cylinder(c) => c.area(), - Shape::Disk(d) => d.area(), - Shape::Triangle(t) => t.area(), - Shape::BilinearPatch(b) => b.area(), - Shape::Curve(c) => c.area(), - } - } - - pub fn pdf(&self, interaction: &dyn Interaction) -> Float { - match self { - Shape::None => 0., - Shape::Sphere(s) => s.pdf(interaction), - Shape::Cylinder(c) => c.pdf(interaction), - Shape::Disk(d) => d.pdf(interaction), - Shape::Triangle(t) => t.pdf(interaction), - Shape::BilinearPatch(b) => b.pdf(interaction), - Shape::Curve(c) => c.pdf(interaction), - } - } - - pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { - match self { - Shape::None => 0., - Shape::Sphere(s) => s.pdf_from_context(ctx, wi), - Shape::Cylinder(c) => c.pdf_from_context(ctx, wi), - Shape::Disk(d) => d.pdf_from_context(ctx, wi), - Shape::Triangle(t) => t.pdf_from_context(ctx, wi), - Shape::BilinearPatch(b) => b.pdf_from_context(ctx, wi), - Shape::Curve(c) => c.pdf_from_context(ctx, wi), - } - } - - pub fn sample(&self, u: Point2f) -> Option { - match self { - Shape::None => None, - Shape::Sphere(s) => s.sample(u), - Shape::Cylinder(c) => c.sample(u), - Shape::Disk(d) => d.sample(u), - Shape::Triangle(t) => t.sample(u), - Shape::BilinearPatch(b) => b.sample(u), - Shape::Curve(c) => c.sample(u), - } - } - - pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { - match self { - Shape::None => None, - Shape::Sphere(s) => s.sample_from_context(ctx, u), - Shape::Cylinder(c) => c.sample_from_context(ctx, u), - Shape::Disk(d) => d.sample_from_context(ctx, u), - Shape::Triangle(t) => t.sample_from_context(ctx, u), - Shape::BilinearPatch(b) => b.sample_from_context(ctx, u), - Shape::Curve(c) => c.sample_from_context(ctx, u), - } +impl Default for Shape { + fn default() -> Self { + Shape::Sphere(SphereShape::default()) } } diff --git a/src/shapes/sphere.rs b/src/shapes/sphere.rs index 828b379..d7ca001 100644 --- a/src/shapes/sphere.rs +++ b/src/shapes/sphere.rs @@ -1,8 +1,9 @@ use super::{ Bounds3f, DirectionCone, Float, Interaction, Normal3f, PI, Point2f, Point3f, Point3fi, - QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, SphereShape, - SurfaceInteraction, Transform, Vector3f, Vector3fi, + QuadricIntersection, Ray, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait, + SphereShape, SurfaceInteraction, Transform, Vector3f, Vector3fi, }; +use crate::core::interaction::InteractionTrait; use crate::core::pbrt::{clamp_t, gamma}; use crate::geometry::{Frame, Sqrt, VectorLike, spherical_direction}; use crate::utils::interval::Interval; @@ -190,8 +191,10 @@ impl SphereShape { flip_normal, ) } +} - pub fn bounds(&self) -> Bounds3f { +impl ShapeTrait for SphereShape { + fn bounds(&self) -> Bounds3f { self.render_from_object .apply_to_bounds(Bounds3f::from_points( Point3f::new(-self.radius, -self.radius, self.z_min), @@ -199,19 +202,19 @@ impl SphereShape { )) } - pub fn normal_bounds(&self) -> DirectionCone { + fn normal_bounds(&self) -> DirectionCone { DirectionCone::entire_sphere() } - pub fn area(&self) -> Float { + fn area(&self) -> Float { self.phi_max * self.radius * (self.z_max - self.z_min) } - pub fn pdf(&self, _interaction: &dyn Interaction) -> Float { + fn pdf(&self, _interaction: &Interaction) -> Float { 1. / self.area() } - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + fn intersect(&self, ray: &Ray, t_max: Option) -> Option { let t = t_max.unwrap_or(Float::INFINITY); if let Some(isect) = self.basic_intersect(ray, t) { let intr = self.interaction_from_intersection(isect.clone(), -ray.d, ray.time); @@ -221,7 +224,7 @@ impl SphereShape { } } - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { if let Some(t) = t_max { self.basic_intersect(ray, t).is_some() } else { @@ -229,7 +232,7 @@ impl SphereShape { } } - pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { + fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { let p_center = self .object_from_render .apply_to_point(Point3f::new(0., 0., 0.)); @@ -257,7 +260,7 @@ impl SphereShape { 1. / (2. * PI * one_minus_cos_theta_max) } - pub fn sample(&self, u: Point2f) -> Option { + fn sample(&self, u: Point2f) -> Option { let p_obj = Point3f::new(0., 0., 0.) + self.radius * sample_uniform_sphere(u); let mut p_obj_vec = Vector3f::from(p_obj); p_obj_vec *= self.radius / p_obj.distance(Point3f::zero()); @@ -289,7 +292,7 @@ impl SphereShape { }) } - pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { let p_center = self.render_from_object.apply_to_point(Point3f::zero()); let p_origin = ctx.offset_ray_origin_from_point(p_center); if p_origin.distance_squared(p_center) <= square(self.radius) { diff --git a/src/shapes/triangle.rs b/src/shapes/triangle.rs index 13c3550..b1751da 100644 --- a/src/shapes/triangle.rs +++ b/src/shapes/triangle.rs @@ -1,8 +1,9 @@ use super::{ Bounds3f, DirectionCone, Float, Interaction, Normal3f, Point2f, Point3f, Point3fi, Ray, - ShapeIntersection, ShapeSample, ShapeSampleContext, SurfaceInteraction, TriangleIntersection, - TriangleShape, Vector2f, Vector3f, + ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait, SurfaceInteraction, + TriangleIntersection, TriangleShape, Vector2f, Vector3f, }; +use crate::core::interaction::InteractionTrait; use crate::core::pbrt::gamma; use crate::geometry::{Sqrt, Tuple, VectorLike, spherical_triangle_area}; use crate::utils::math::{difference_of_products, square}; @@ -323,13 +324,15 @@ impl TriangleShape { isect.dndu = dndu; isect.dndv = dndv; } +} - pub fn bounds(&self) -> Bounds3f { +impl ShapeTrait for TriangleShape { + fn bounds(&self) -> Bounds3f { let [p0, p1, p2] = self.get_data().vertices; Bounds3f::from_points(p0, p1).union_point(p2) } - pub fn normal_bounds(&self) -> DirectionCone { + fn normal_bounds(&self) -> DirectionCone { let data = self.get_data(); let mut n = data.normal; if let Some([n0, n1, n2]) = data.normals { @@ -340,15 +343,15 @@ impl TriangleShape { DirectionCone::new_from_vector(Vector3f::from(n)) } - pub fn area(&self) -> Float { + fn area(&self) -> Float { self.get_data().area } - pub fn pdf(&self, _interaction: &dyn Interaction) -> Float { + fn pdf(&self, _interaction: &Interaction) -> Float { 1. / self.area() } - pub fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { + fn pdf_from_context(&self, ctx: &ShapeSampleContext, wi: Vector3f) -> Float { let solid_angle = self.solid_angle(ctx.p()); if (Self::MIN_SPHERICAL_SAMPLE_AREA..Self::MAX_SPHERICAL_SAMPLE_AREA).contains(&solid_angle) { @@ -384,7 +387,7 @@ impl TriangleShape { pdf } - pub fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + fn sample_from_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { let data = self.get_data(); let [p0, p1, p2] = data.vertices; let solid_angle = self.solid_angle(ctx.p()); @@ -463,7 +466,7 @@ impl TriangleShape { }) } - pub fn sample(&self, u: Point2f) -> Option { + fn sample(&self, u: Point2f) -> Option { let data = self.get_data(); let [p0, p1, p2] = data.vertices; let [uv0, uv1, uv2] = data.uvs; @@ -494,18 +497,15 @@ impl TriangleShape { }) } - pub fn intersect(&self, ray: &Ray, t_max: Option) -> Option { + fn intersect(&self, ray: &Ray, t_max: Option) -> Option { self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY)) .map(|ti| { let intr = self.interaction_from_intersection(ti, ray.time, -ray.d); - ShapeIntersection { - intr: Box::new(intr), - t_hit: ti.t, - } + ShapeIntersection { intr, t_hit: ti.t } }) } - pub fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { + fn intersect_p(&self, ray: &Ray, t_max: Option) -> bool { self.intersect_triangle(ray, t_max.unwrap_or(Float::INFINITY)) .is_some() } diff --git a/src/spectra/mod.rs b/src/spectra/mod.rs index 14bdf57..e6ee8b8 100644 --- a/src/spectra/mod.rs +++ b/src/spectra/mod.rs @@ -10,7 +10,7 @@ use enum_dispatch::enum_dispatch; pub use data::*; pub use rgb::*; -use sampled::{CIE_Y_INTEGRAL, LAMBDA_MAX, LAMBDA_MIN}; +pub use sampled::{CIE_Y_INTEGRAL, LAMBDA_MAX, LAMBDA_MIN}; pub use sampled::{N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths}; pub use simple::*; // CIE_X, etc // @@ -49,11 +49,7 @@ impl Spectrum { } pub fn sample(&self, wavelengths: &SampledWavelengths) -> SampledSpectrum { - let mut s = SampledSpectrum::default(); - for i in 0..N_SPECTRUM_SAMPLES { - s[i] = self.evaluate(wavelengths[i]); - } - s + SampledSpectrum::from_fn(|i| self.evaluate(wavelengths[i])) } pub fn inner_product(&self, other: &Spectrum) -> Float { @@ -65,4 +61,8 @@ impl Spectrum { } integral } + + pub fn is_constant(&self) -> bool { + matches!(self, Spectrum::Constant(_)) + } } diff --git a/src/spectra/sampled.rs b/src/spectra/sampled.rs index 1ba088c..4355712 100644 --- a/src/spectra/sampled.rs +++ b/src/spectra/sampled.rs @@ -3,6 +3,8 @@ use std::ops::{ Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign, }; +use super::{cie_x, cie_y, cie_z}; + pub const CIE_Y_INTEGRAL: Float = 106.856895; pub const N_SPECTRUM_SAMPLES: usize = 1200; @@ -29,7 +31,7 @@ impl SampledSpectrum { } } - #[inline] + #[inline(always)] pub fn from_fn(cb: F) -> Self where F: FnMut(usize) -> Float, @@ -97,6 +99,12 @@ impl SampledSpectrum { pub fn safe_div(a: &SampledSpectrum, b: &SampledSpectrum) -> SampledSpectrum { SampledSpectrum::from_fn(|i| if b[i] != 0. { a[i] / b[i] } else { 0. }) } + + pub fn y(&self, lambda: &SampledWavelengths) -> Float { + let ys = cie_y().sample(lambda); + let pdf = lambda.pdf(); + SampledSpectrum::safe_div(&(ys * *self), &pdf).average() / CIE_Y_INTEGRAL + } } impl<'a> IntoIterator for &'a SampledSpectrum { @@ -291,12 +299,18 @@ impl SampledWavelengths { true } - pub fn terminate_secondary(&mut self) { - if !self.secondary_terminated() { - for i in 1..N_SPECTRUM_SAMPLES { - self.pdf[i] = 0.0; - } - self.pdf[0] /= N_SPECTRUM_SAMPLES as Float; + pub fn terminate_secondary(self) -> Self { + if self.secondary_terminated() { + return self; + } + + let mut new_pdf = [0.0; N_SPECTRUM_SAMPLES]; + + new_pdf[0] = self.pdf[0] / (N_SPECTRUM_SAMPLES as Float); + + Self { + pdf: new_pdf, + ..self } } diff --git a/src/spectra/simple.rs b/src/spectra/simple.rs index a60ab2c..9ff3e8b 100644 --- a/src/spectra/simple.rs +++ b/src/spectra/simple.rs @@ -3,6 +3,7 @@ use crate::core::pbrt::Float; use crate::spectra::{ N_SPECTRUM_SAMPLES, SampledSpectrum, SampledWavelengths, Spectrum, SpectrumTrait, }; +use std::hash::{Hash, Hasher}; #[derive(Debug, Clone, Copy)] pub struct ConstantSpectrum { @@ -108,6 +109,41 @@ impl DenselySampledSpectrum { } r } + + pub fn scale(&mut self, factor: Float) { + for v in &mut self.values { + *v *= factor; + } + } +} + +impl PartialEq for DenselySampledSpectrum { + fn eq(&self, other: &Self) -> bool { + if self.lambda_min != other.lambda_min + || self.lambda_max != other.lambda_max + || self.values.len() != other.values.len() + { + return false; + } + + self.values + .iter() + .zip(&other.values) + .all(|(a, b)| a.to_bits() == b.to_bits()) + } +} + +impl Eq for DenselySampledSpectrum {} + +impl Hash for DenselySampledSpectrum { + fn hash(&self, state: &mut H) { + self.lambda_min.hash(state); + self.lambda_max.hash(state); + + for v in &self.values { + v.to_bits().hash(state); + } + } } impl SpectrumTrait for DenselySampledSpectrum { @@ -216,6 +252,12 @@ impl BlackbodySpectrum { numerator / denominator } } + + pub fn sample(&self, lambda: &SampledWavelengths) -> SampledSpectrum { + SampledSpectrum::from_fn(|i| { + Self::planck_law(lambda[i], self.temperature) * self.normalization_factor + }) + } } impl SpectrumTrait for BlackbodySpectrum { diff --git a/src/utils/color.rs b/src/utils/color.rs index 2953e59..06f03f9 100644 --- a/src/utils/color.rs +++ b/src/utils/color.rs @@ -252,7 +252,7 @@ impl fmt::Display for XYZ { } } -#[derive(Debug, Copy, Clone)] +#[derive(Debug, Default, Copy, Clone)] pub struct RGB { pub r: Float, pub g: Float, diff --git a/src/utils/containers.rs b/src/utils/containers.rs index f4a2751..829577f 100644 --- a/src/utils/containers.rs +++ b/src/utils/containers.rs @@ -1,16 +1,27 @@ +use crate::core::pbrt::{Float, lerp}; use std::collections::hash_map::RandomState; use std::hash::{BuildHasher, Hash, Hasher}; -use std::ops::{Index, IndexMut}; +use std::ops::{Add, Index, IndexMut, Mul, Sub}; use std::sync::RwLock; use crate::geometry::{ - Bounds2i, Bounds3i, Point2i, Point3f, Point3i, Vector2i, Vector3f, Vector3i, + Bounds2i, Bounds3f, Bounds3i, Point2i, Point3f, Point3i, Vector2i, Vector3f, Vector3i, }; -#[derive(Debug)] +pub trait Interpolatable: + Copy + Default + Add + Sub + Mul +{ +} + +impl Interpolatable for T where + T: Copy + Default + Add + Sub + Mul +{ +} + +#[derive(Debug, Clone)] pub struct Array2D { - values: Vec, - extent: Bounds2i, + pub values: Vec, + pub extent: Bounds2i, } impl Array2D { @@ -24,16 +35,27 @@ impl Array2D { Self { values, extent } } - pub fn new_with_dims(nx: i32, ny: i32) -> Self + pub fn new_with_dims(nx: usize, ny: usize) -> Self where T: Default, { Self::new(Bounds2i::from_points( Point2i::new(0, 0), - Point2i::new(nx, ny), + Point2i::new(nx as i32, ny as i32), )) } + pub fn new_filled(width: usize, height: usize, value: T) -> Self + where + T: Clone, + { + let extent = Bounds2i::from_points( + Point2i::new(0, 0), + Point2i::new(width as i32, height as i32), + ); + Self::new_from_bounds(extent, value) + } + pub fn new_from_bounds(extent: Bounds2i, default_val: T) -> Self where T: Clone, @@ -92,6 +114,20 @@ impl IndexMut for Array2D { } } +impl Index<(usize, usize)> for Array2D { + type Output = T; + + fn index(&self, (x, y): (usize, usize)) -> &Self::Output { + &self[(x as i32, y as i32)] + } +} + +impl IndexMut<(usize, usize)> for Array2D { + fn index_mut(&mut self, (x, y): (usize, usize)) -> &mut Self::Output { + &mut self[(x as i32, y as i32)] + } +} + impl Index<(i32, i32)> for Array2D { type Output = T; fn index(&self, index: (i32, i32)) -> &Self::Output { @@ -104,3 +140,168 @@ impl IndexMut<(i32, i32)> for Array2D { self.index_mut(Point2i::new(index.0, index.1)) } } + +#[derive(Debug, Clone)] +pub struct SampledGrid { + values: Vec, + nx: i32, + ny: i32, + nz: i32, +} + +impl SampledGrid { + pub fn new(values: Vec, nx: i32, ny: i32, nz: i32) -> Self { + assert_eq!( + values.len(), + (nx * ny * nz) as usize, + "Grid dimensions do not match data size" + ); + Self { values, nx, ny, nz } + } + + pub fn empty() -> Self { + Self { + values: Vec::new(), + nx: 0, + ny: 0, + nz: 0, + } + } + + pub fn bytes_allocated(&self) -> usize { + self.values.len() * std::mem::size_of::() + } + + pub fn x_size(&self) -> i32 { + self.nx + } + pub fn y_size(&self) -> i32 { + self.ny + } + pub fn z_size(&self) -> i32 { + self.nz + } + + pub fn lookup_int_convert(&self, p: Point3i, convert: F) -> U + where + F: Fn(&T) -> U, + U: Default, + { + let sample_bounds = Bounds3i::from_points( + Point3i::new(0, 0, 0), + Point3i::new(self.nx, self.ny, self.nz), + ); + + if !sample_bounds.contains_exclusive(p) { + return U::default(); + } + + let idx = (p.z() * self.ny + p.y()) * self.nx + p.x(); + convert(&self.values[idx as usize]) + } + + pub fn lookup_int(&self, p: Point3i) -> T + where + T: Clone + Default, + { + self.lookup_int_convert(p, |v| v.clone()) + } + + pub fn lookup_convert(&self, p: Point3f, convert: F) -> U + where + F: Fn(&T) -> U + Copy, + U: Interpolatable + Default, + { + let p_samples = Point3f::new( + p.x() * self.nx as Float - 0.5, + p.y() * self.ny as Float - 0.5, + p.z() * self.nz as Float - 0.5, + ); + + let pi = Point3i::from(p_samples.floor()); + let d = p_samples - Point3f::from(pi); + + // Helper to retrieve corners with conversion + let lk = |offset: Vector3i| -> U { self.lookup_int_convert(pi + offset, convert) }; + + // Trilinear interpolation + let d00 = lerp( + d.x(), + lk(Vector3i::new(0, 0, 0)), + lk(Vector3i::new(1, 0, 0)), + ); + let d10 = lerp( + d.x(), + lk(Vector3i::new(0, 1, 0)), + lk(Vector3i::new(1, 1, 0)), + ); + let d01 = lerp( + d.x(), + lk(Vector3i::new(0, 0, 1)), + lk(Vector3i::new(1, 0, 1)), + ); + let d11 = lerp( + d.x(), + lk(Vector3i::new(0, 1, 1)), + lk(Vector3i::new(1, 1, 1)), + ); + + lerp(d.z(), lerp(d.y(), d00, d10), lerp(d.y(), d01, d11)) + } + + pub fn lookup(&self, p: Point3f) -> T + where + T: Interpolatable + Clone, + { + self.lookup_convert(p, |v| v.clone()) + } + + pub fn max_value_convert(&self, bounds: Bounds3f, convert: F) -> U + where + F: Fn(&T) -> U, + U: PartialOrd + Default, + { + let ps = [ + Point3f::new( + bounds.p_min.x() * self.nx as Float - 0.5, + bounds.p_min.y() * self.ny as Float - 0.5, + bounds.p_min.z() * self.nz as Float - 0.5, + ), + Point3f::new( + bounds.p_max.x() * self.nx as Float - 0.5, + bounds.p_max.y() * self.ny as Float - 0.5, + bounds.p_max.z() * self.nz as Float - 0.5, + ), + ]; + + let pi_min = Point3i::from(ps[0].floor()).max(Point3i::new(0, 0, 0)); + let pi_max = (Point3i::from(ps[1].floor()) + Vector3i::new(1, 1, 1)).min(Point3i::new( + self.nx - 1, + self.ny - 1, + self.nz - 1, + )); + + // Initialize with the first voxel + let mut max_value = self.lookup_int_convert(pi_min, &convert); + + for z in pi_min.z()..=pi_max.z() { + for y in pi_min.y()..=pi_max.y() { + for x in pi_min.x()..=pi_max.x() { + let val = self.lookup_int_convert(Point3i::new(x, y, z), &convert); + if val > max_value { + max_value = val; + } + } + } + } + + max_value + } + + pub fn max_value(&self, bounds: Bounds3f) -> T + where + T: PartialOrd + Default + Clone, + { + self.max_value_convert(bounds, |v| v.clone()) + } +} diff --git a/src/utils/hash.rs b/src/utils/hash.rs index 3c2b7f2..e48893f 100644 --- a/src/utils/hash.rs +++ b/src/utils/hash.rs @@ -93,9 +93,9 @@ macro_rules! hash_values { } } -pub fn hash_float(data: T) -> Float +pub fn hash_float(data: &T) -> Float where - T: Copy, // Ensure it's plain old data + T: Copy + ?Sized, // Ensure it's plain old data { let h = hash_buffer(&data, 0); (h as u32) as Float * 2.3283064365386963e-10 // 0x1p-32f diff --git a/src/utils/math.rs b/src/utils/math.rs index 5fa2b0b..1c4829e 100644 --- a/src/utils/math.rs +++ b/src/utils/math.rs @@ -4,7 +4,7 @@ use crate::core::pbrt::{ Float, FloatBitOps, FloatBits, ONE_MINUS_EPSILON, PI, PI_OVER_4, clamp_t, evaluate_polynomial, lerp, }; -use crate::geometry::{Point, Point2f, Point2i, Vector, Vector3f}; +use crate::geometry::{Point, Point2f, Point2i, Vector, Vector3f, VectorLike}; use crate::utils::hash::{hash_buffer, mix_bits}; use crate::utils::sobol::{SOBOL_MATRICES_32, VDC_SOBOL_MATRICES, VDC_SOBOL_MATRICES_INV}; @@ -231,6 +231,15 @@ pub fn quadratic(a: Float, b: Float, c: Float) -> Option<(Float, Float)> { Some((t0, t1)) } +pub fn smooth_step(x: Float, a: Float, b: Float) -> Float { + if a == b { + if x < a { return 0. } else { return 1. } + } + + let t = clamp_t((x - a) / (b - a), 0., 1.); + t * t * (3. - 2. * t) +} + pub fn linear_least_squares( a: [[Float; R]; N], b: [[Float; R]; N], @@ -362,6 +371,43 @@ pub fn catmull_rom_weights(nodes: &[Float], x: Float) -> Option<(usize, [Float; Some((offset, weights)) } +pub fn equal_area_sphere_to_square(d: Vector3f) -> Point2f { + debug_assert!(d.norm_squared() > 0.999 && d.norm_squared() < 1.001); + let x = d.x().abs(); + let y = d.y().abs(); + let z = d.z().abs(); + let r = safe_sqrt(1. - z); + let a = x.max(y); + let b = if a == 0. { 0. } else { x.min(y) / a }; + let t1 = 0.406758566246788489601959989e-5; + let t2 = 0.636226545274016134946890922156; + let t3 = 0.61572017898280213493197203466e-2; + let t4 = -0.247333733281268944196501420480; + let t5 = 0.881770664775316294736387951347e-1; + let t6 = 0.419038818029165735901852432784e-1; + let t7 = -0.251390972343483509333252996350e-1; + let mut phi = evaluate_polynomial(b, &[t1, t2, t3, t4, t5, t6, t7]) + .expect("Could not evaluate polynomial"); + + if x < y { + phi = 1. - phi; + } + + let mut v = phi * r; + let mut u = r - v; + + if d.z() < 0. { + mem::swap(&mut u, &mut v); + u = 1. - u; + v = 1. - v; + } + + u = u.copysign(d.x()); + v = v.copysign(d.y()); + + Point2f::new(0.5 * (u + 1.), 0.5 * (v + 1.)) +} + pub fn equal_area_square_to_sphere(p: Point2f) -> Vector3f { assert!(p.x() >= 0. && p.x() <= 1. && p.y() >= 0. && p.y() <= 1.); @@ -727,7 +773,6 @@ pub fn compute_radical_inverse_permutations(seed: u64) -> Vec } pub fn scrambled_radical_inverse(base_index: usize, mut a: u64, perm: &DigitPermutation) -> Float { - // Access global PRIMES array (or pass base explicitly if preferred) let base = PRIMES[base_index] as u64; let limit = (u64::MAX / base).saturating_sub(base); @@ -737,7 +782,6 @@ pub fn scrambled_radical_inverse(base_index: usize, mut a: u64, perm: &DigitPerm let mut reversed_digits = 0u64; let mut digit_index = 0; - // Note: 'a' is not part of the condition; we pad with 0s if a runs out. while 1.0 - ((base as Float - 1.0) * inv_base_m) < 1.0 && reversed_digits < limit { let next = a / base; let digit_value = (a - next * base) as i32; @@ -784,7 +828,6 @@ pub fn owen_scrambled_radical_inverse(base_index: usize, mut a: u64, hash: u32) } pub fn permutation_element(mut i: u32, l: u32, p: u32) -> u32 { - // 1. Compute the mask 'w' (next power of 2 minus 1) let mut w = l - 1; w |= w >> 1; w |= w >> 2; @@ -792,9 +835,7 @@ pub fn permutation_element(mut i: u32, l: u32, p: u32) -> u32 { w |= w >> 8; w |= w >> 16; - // 2. Cycle Walking loop loop { - // We use wrapping_mul to replicate C++ unsigned overflow behavior i ^= p; i = i.wrapping_mul(0xe170893d); i ^= p >> 16; @@ -815,14 +856,13 @@ pub fn permutation_element(mut i: u32, l: u32, p: u32) -> u32 { i &= w; i ^= i >> 5; - // If the resulting index is within range [0, l), we are done. - // Otherwise, we loop again (cycle walking). + // If index is within range [0, l), we are done. + // Otherwise, we loop again. if i < l { break; } } - // 3. Final shift by p (i.wrapping_add(p)) % l } diff --git a/src/utils/mipmap.rs b/src/utils/mipmap.rs index b1013fd..b06bcbe 100644 --- a/src/utils/mipmap.rs +++ b/src/utils/mipmap.rs @@ -76,11 +76,11 @@ impl MIPMapSample for Float { } fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self { - image.bilerp_channel(st, 0, wrap) + image.bilerp_channel_with_wrap(st, 0, wrap) } fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self { - image.get_channel(st, 0, wrap) + image.get_channel_with_wrap(st, 0, wrap) } } @@ -92,12 +92,12 @@ impl MIPMapSample for RGB { fn sample_bilerp(image: &Image, st: Point2f, wrap: WrapMode2D) -> Self { let nc = image.n_channels(); if nc >= 3 { - let r = image.bilerp_channel(st, 0, wrap); - let g = image.bilerp_channel(st, 1, wrap); - let b = image.bilerp_channel(st, 2, wrap); + let r = image.bilerp_channel_with_wrap(st, 0, wrap); + let g = image.bilerp_channel_with_wrap(st, 1, wrap); + let b = image.bilerp_channel_with_wrap(st, 2, wrap); RGB::new(r, g, b) } else { - let v = image.bilerp_channel(st, 0, wrap); + let v = image.bilerp_channel_with_wrap(st, 0, wrap); RGB::new(v, v, v) } } @@ -105,12 +105,12 @@ impl MIPMapSample for RGB { fn sample_texel(image: &Image, st: Point2i, wrap: WrapMode2D) -> Self { let nc = image.n_channels(); if nc >= 3 { - let r = image.get_channel(st, 0, wrap); - let g = image.get_channel(st, 1, wrap); - let b = image.get_channel(st, 2, wrap); + let r = image.get_channel_with_wrap(st, 0, wrap); + let g = image.get_channel_with_wrap(st, 1, wrap); + let b = image.get_channel_with_wrap(st, 2, wrap); RGB::new(r, g, b) } else { - let v = image.get_channel(st, 0, wrap); + let v = image.get_channel_with_wrap(st, 0, wrap); RGB::new(v, v, v) } } @@ -214,8 +214,8 @@ impl MIPMap { let longer_len = dst0.norm(); let mut shorter_len = dst1.norm(); - // If the ellipse is too thin, we fatten the minor axis to limit the number - // of texels we have to filter, trading accuracy for performance/cache hits. + // If ellipse is too thin, fatten the minor axis to limit the number + // of texels if shorter_len * self.options.max_anisotropy < longer_len && shorter_len > 0.0 { let scale = longer_len / (shorter_len * self.options.max_anisotropy); dst1 *= scale; @@ -226,7 +226,6 @@ impl MIPMap { return self.bilerp(0, st); } - // EWA also interpolates between two MIP levels, but does complex filtering on both. let lod = (self.levels() as Float - 1.0 + shorter_len.log2()).max(0.0); let ilod = lod.floor() as usize; diff --git a/src/utils/mod.rs b/src/utils/mod.rs index 98715c5..bb5fc5e 100644 --- a/src/utils/mod.rs +++ b/src/utils/mod.rs @@ -1,3 +1,5 @@ +use std::sync::atomic::{AtomicU64, Ordering}; + pub mod color; pub mod colorspace; pub mod containers; @@ -29,3 +31,49 @@ where } i } + +#[derive(Debug)] +pub struct AtomicFloat { + bits: AtomicU64, +} + +impl AtomicFloat { + pub fn new(value: f64) -> Self { + Self { + bits: AtomicU64::new(value.to_bits()), + } + } + + pub fn load(&self) -> f64 { + f64::from_bits(self.bits.load(Ordering::Relaxed)) + } + + pub fn store(&self, value: f64) { + self.bits.store(value.to_bits(), Ordering::Relaxed); + } + + pub fn add(&self, value: f64) { + let mut current_bits = self.bits.load(Ordering::Relaxed); + loop { + let current_val = f64::from_bits(current_bits); + let new_val = current_val + value; + let new_bits = new_val.to_bits(); + + match self.bits.compare_exchange_weak( + current_bits, + new_bits, + Ordering::Relaxed, + Ordering::Relaxed, + ) { + Ok(_) => break, + Err(x) => current_bits = x, + } + } + } +} + +impl Default for AtomicFloat { + fn default() -> Self { + Self::new(0.0) + } +} diff --git a/src/utils/rng.rs b/src/utils/rng.rs index 4d2e083..005e55a 100644 --- a/src/utils/rng.rs +++ b/src/utils/rng.rs @@ -1,13 +1,11 @@ #[derive(Debug, Clone)] pub struct Rng { - // Internal state (PCG32 usually has state and sequence) state: u64, inc: u64, } impl Default for Rng { fn default() -> Self { - // Initialize with default PBRT constants Self { state: 0x853c49e6748fea9b, inc: 0xda3e39cb94b95bdb, @@ -20,14 +18,13 @@ impl Rng { self.state = 0; self.inc = (init_seq << 1) | 1; self.uniform_u32(); // Warm up - self.state = self.state.wrapping_add(0x853c49e6748fea9b); // PCG32_DEFAULT_STATE + // PCG32_DEFAULT_STATE + self.state = self.state.wrapping_add(0x853c49e6748fea9b); self.uniform_u32(); } pub fn advance(&mut self, delta: u64) { - // Implementation of PCG32 advance/seek - // ... (Math heavy, requires modular exponentiation for matrix mult) ... - // Simplest approximation for non-PCG: + // TODO: Implementation of PCG32 advance/seek for _ in 0..delta { self.uniform_u32(); } @@ -35,7 +32,7 @@ impl Rng { pub fn uniform(&mut self) -> T where - T: From, // Simplified trait bound + T: From, { // Convert u32 to float [0,1) let v = self.uniform_u32(); diff --git a/src/utils/sampling.rs b/src/utils/sampling.rs index 646ed61..6a890b4 100644 --- a/src/utils/sampling.rs +++ b/src/utils/sampling.rs @@ -1,7 +1,7 @@ use super::math::safe_sqrt; use crate::check_rare; use crate::core::pbrt::{ - Float, INV_2_PI, INV_PI, ONE_MINUS_EPSILON, PI, PI_OVER_2, PI_OVER_4, clamp_t, + Float, INV_2_PI, INV_4_PI, INV_PI, ONE_MINUS_EPSILON, PI, PI_OVER_2, PI_OVER_4, clamp_t, evaluate_polynomial, find_interval, lerp, }; use crate::core::pbrt::{RARE_EVENT_CONDITION_MET, RARE_EVENT_TOTAL_CALLS}; @@ -41,6 +41,16 @@ pub fn invert_bilinear_sample(p: Point2f, w: &[Float]) -> Point2f { ) } +pub fn power_heuristic(nf: i32, f_pdf: Float, ng: i32, g_pdf: Float) -> Float { + let f = nf as Float * f_pdf; + let g = ng as Float * g_pdf; + square(f) / (square(f) + square(g)) +} + +pub fn uniform_sphere_pdf() -> Float { + INV_4_PI +} + pub fn bilinear_pdf(p: Point2f, w: &[Float]) -> Float { if p.x() < 0. || p.x() > 1. || p.y() < 0. || p.y() > 1. { return 0.; @@ -779,7 +789,9 @@ impl PiecewiseConstant2D { Self::new(data, data.x_size(), data.y_size(), domain) } - pub fn new_with_data(data: &Array2D, nx: usize, ny: usize) -> Self { + pub fn new_with_data(data: &Array2D) -> Self { + let nx = data.x_size(); + let ny = data.y_size(); Self::new( data, nx, @@ -824,6 +836,317 @@ impl PiecewiseConstant2D { } } +#[derive(Debug, Clone)] +pub struct SummedAreaTable { + sum: Array2D, +} + +impl SummedAreaTable { + pub fn new(values: &Array2D) -> Self { + let width = values.x_size(); + let height = values.y_size(); + + let mut sum = Array2D::::new_with_dims(width, height); + sum[(0, 0)] = values[(0, 0)] as f64; + + for x in 1..width { + sum[(x, 0)] = values[(x, 0)] as f64 + sum[(x - 1, 0)]; + } + + for y in 1..height { + sum[(0, y)] = values[(0, y)] as f64 + sum[(0, y - 1)]; + } + + for y in 1..height { + for x in 1..width { + let term = values[(x, y)] as f64; + let left = sum[(x - 1, y)]; + let up = sum[(x, y - 1)]; + let diag = sum[(x - 1, y - 1)]; + + sum[(x, y)] = term + left + up - diag; + } + } + + Self { sum } + } + + pub fn integral(&self, extent: Bounds2f) -> Float { + let s = self.lookup(extent.p_max.x(), extent.p_max.y()) + - self.lookup(extent.p_min.x(), extent.p_max.y()) + + self.lookup(extent.p_min.x(), extent.p_min.y()) + - self.lookup(extent.p_max.x(), extent.p_min.y()); + + let total_area = (self.sum.x_size() * self.sum.y_size()) as f64; + + (s / total_area).max(0.0) as Float + } + + fn lookup(&self, mut x: Float, mut y: Float) -> f64 { + x *= self.sum.x_size() as Float; + y *= self.sum.y_size() as Float; + + let x0 = x as i32; + let y0 = y as i32; + + let v00 = self.lookup_int(x0, y0); + let v10 = self.lookup_int(x0 + 1, y0); + let v01 = self.lookup_int(x0, y0 + 1); + let v11 = self.lookup_int(x0 + 1, y0 + 1); + + let dx = (x - x0 as Float) as f64; + let dy = (y - y0 as Float) as f64; + + (1.0 - dx) * (1.0 - dy) * v00 + + dx * (1.0 - dy) * v10 + + (1.0 - dx) * dy * v01 + + dx * dy * v11 + } + + fn lookup_int(&self, x: i32, y: i32) -> f64 { + if x == 0 || y == 0 { + return 0.0; + } + + let ix = (x - 1).min(self.sum.x_size() as i32 - 1) as usize; + let iy = (y - 1).min(self.sum.y_size() as i32 - 1) as usize; + + self.sum[(ix, iy)] + } +} + +#[derive(Debug, Clone)] +pub struct WindowedPiecewiseConstant2D { + sat: SummedAreaTable, + func: Array2D, +} + +impl WindowedPiecewiseConstant2D { + pub fn new(func: Array2D) -> Self { + let sat = SummedAreaTable::new(&func); + Self { sat, func } + } + + pub fn sample(&self, u: Point2f, b: Bounds2f) -> Option<(Point2f, Float)> { + let b_int = self.sat.integral(b); + if b_int == 0.0 { + return None; + } + + let px = |x: Float| -> Float { + let mut bx = b; + bx.p_max[0] = x; + self.sat.integral(bx) / b_int + }; + + let nx = self.func.x_size(); + let px_val = Self::sample_bisection(px, u.x(), b.p_min.x(), b.p_max.x(), nx); + + let nx_f = nx as Float; + + let x_start = (px_val * nx_f).floor() / nx_f; + let x_end = (px_val * nx_f).ceil() / nx_f; + + let mut b_cond = Bounds2f::from_points( + Point2f::new(x_start, b.p_min.y()), + Point2f::new(x_end, b.p_max.y()), + ); + + if b_cond.p_min.x() == b_cond.p_max.x() { + b_cond.p_max[0] += 1.0 / nx_f; + } + + let cond_integral = self.sat.integral(b_cond); + if cond_integral == 0.0 { + return None; + } + + let py = |y: Float| -> Float { + let mut by = b_cond; + by.p_max[1] = y; + self.sat.integral(by) / cond_integral + }; + + let ny = self.func.y_size(); + let py_val = Self::sample_bisection(py, u.y(), b.p_min.y(), b.p_max.y(), ny); + + let p = Point2f::new(px_val, py_val); + + let pdf = self.eval(p) / b_int; + + Some((p, pdf)) + } + + pub fn pdf(&self, p: Point2f, b: Bounds2f) -> Float { + let func_int = self.sat.integral(b); + if func_int == 0.0 { + return 0.0; + } + self.eval(p) / func_int + } + + fn eval(&self, p: Point2f) -> Float { + let nx = self.func.x_size(); + let ny = self.func.y_size(); + + let ix = ((p.x() * nx as Float) as i32).min(nx as i32 - 1).max(0) as usize; + let iy = ((p.y() * ny as Float) as i32).min(ny as i32 - 1).max(0) as usize; + + self.func[(ix, iy)] + } + + fn sample_bisection(p_func: F, u: Float, mut min: Float, mut max: Float, n: usize) -> Float + where + F: Fn(Float) -> Float, + { + let n_f = n as Float; + + while (n_f * max).ceil() - (n_f * min).floor() > 1.0 { + debug_assert!(p_func(min) <= u + 1e-5); + debug_assert!(p_func(max) >= u - 1e-5); + + let mid = (min + max) / 2.0; + if p_func(mid) > u { + max = mid; + } else { + min = mid; + } + } + + let numerator = u - p_func(min); + let denominator = p_func(max) - p_func(min); + + let t = if denominator == 0.0 { + 0.5 + } else { + numerator / denominator + }; + + let res = crate::core::pbrt::lerp(t, min, max); + res.clamp(min, max) + } +} + +#[derive(Debug, Clone, Copy)] +pub struct Bin { + q: Float, + p: Float, + alias: usize, +} + +#[derive(Debug, Clone)] +pub struct AliasTable { + bins: Vec, +} + +impl AliasTable { + pub fn new(weights: &[Float]) -> Self { + let n = weights.len(); + if n == 0 { + return Self { bins: Vec::new() }; + } + let sum: f64 = weights.iter().map(|&w| w as f64).sum(); + assert!(sum > 0.0, "Sum of weights must be positive"); + let mut bins = Vec::with_capacity(n); + for &w in weights { + bins.push(Bin { + p: (w as f64 / sum) as Float, + q: 0.0, + alias: 0, + }); + } + + struct Outcome { + p_hat: f64, + index: usize, + } + + let mut under = Vec::with_capacity(n); + let mut over = Vec::with_capacity(n); + + for (i, bin) in bins.iter().enumerate() { + let p_hat = (bin.p as f64) * (n as f64); + if p_hat < 1.0 { + under.push(Outcome { p_hat, index: i }); + } else { + over.push(Outcome { p_hat, index: i }); + } + } + + while !under.is_empty() && !over.is_empty() { + let un = under.pop().unwrap(); + let ov = over.pop().unwrap(); + + bins[un.index].q = un.p_hat as Float; + bins[un.index].alias = ov.index; + + let p_excess = un.p_hat + ov.p_hat - 1.0; + + if p_excess < 1.0 { + under.push(Outcome { + p_hat: p_excess, + index: ov.index, + }); + } else { + over.push(Outcome { + p_hat: p_excess, + index: ov.index, + }); + } + } + + while let Some(ov) = over.pop() { + bins[ov.index].q = 1.0; + bins[ov.index].alias = ov.index; + } + + while let Some(un) = under.pop() { + bins[un.index].q = 1.0; + bins[un.index].alias = un.index; + } + + Self { bins } + } + + pub fn sample(&self, u: Float) -> (usize, Float, Float) { + let n = self.bins.len(); + + let val = u * (n as Float); + let offset = std::cmp::min(val as usize, n - 1); + + let up = (val - (offset as Float)).min(ONE_MINUS_EPSILON); + + let bin = &self.bins[offset]; + + if up < bin.q { + debug_assert!(bin.p > 0.0); + + let pmf = bin.p; + + let u_remapped = (up / bin.q).min(ONE_MINUS_EPSILON); + + (offset, pmf, u_remapped) + } else { + let alias_idx = bin.alias; + + let alias_p = self.bins[alias_idx].p; + debug_assert!(alias_p > 0.0); + + let u_remapped = ((up - bin.q) / (1.0 - bin.q)).min(ONE_MINUS_EPSILON); + + (alias_idx, alias_p, u_remapped) + } + } + + pub fn size(&self) -> usize { + self.bins.len() + } + + pub fn pmf(&self, index: usize) -> Float { + self.bins[index].p + } +} + #[derive(Debug, Clone)] pub struct PiecewiseLinear2D { size: Vector2i, diff --git a/src/utils/transform.rs b/src/utils/transform.rs index 79e56b4..4ddc812 100644 --- a/src/utils/transform.rs +++ b/src/utils/transform.rs @@ -6,10 +6,10 @@ use std::ops::{Add, Div, Index, IndexMut, Mul}; use std::sync::Arc; use super::color::{RGB, XYZ}; -use super::math::{SquareMatrix, safe_acos}; +use super::math::{SquareMatrix, radians, safe_acos}; use super::quaternion::Quaternion; use crate::core::interaction::{ - Interaction, InteractionData, MediumInteraction, SurfaceInteraction, + Interaction, InteractionData, InteractionTrait, MediumInteraction, SurfaceInteraction, }; use crate::core::pbrt::{Float, gamma}; use crate::geometry::{ @@ -29,7 +29,7 @@ impl Transform { Self { m, m_inv } } - pub fn from_matrix(m: SquareMatrix) -> Result> { + pub fn from_matrix(m: SquareMatrix) -> Result { let inv = m.inverse()?; Ok(Self { m, m_inv: inv }) } @@ -77,6 +77,16 @@ impl Transform { ]); s.determinant() < T::zero() } + + pub fn get_matrix(&self) -> SquareMatrix { + self.m + } +} + +impl Default for Transform { + fn default() -> Self { + Self::identity() + } } impl Transform { @@ -219,46 +229,60 @@ impl Transform { } } - pub fn apply_to_interaction(&self, inter: &dyn Interaction) -> Box { - if let Some(si) = inter.as_any().downcast_ref::() { - let mut ret = si.clone(); - ret.common.pi = self.apply_to_interval(&si.common.pi); - let n = self.apply_to_normal(si.common.n); - ret.common.wo = self.apply_to_vector(si.common.wo).normalize(); + pub fn apply_to_interaction(&self, inter: &Interaction) -> Interaction { + match inter { + Interaction::Surface(si) => { + let mut ret = si.clone(); - ret.dpdu = self.apply_to_vector(si.dpdu); - ret.dpdv = self.apply_to_vector(si.dpdv); - ret.dndu = self.apply_to_normal(si.dndu); - ret.dndv = self.apply_to_normal(si.dndv); - ret.dpdx = self.apply_to_vector(si.dpdx); - ret.dpdy = self.apply_to_vector(si.dpdy); + ret.common.pi = self.apply_to_interval(&si.common.pi); - let shading_n = self.apply_to_normal(si.shading.n); - ret.shading.n = shading_n.normalize(); - ret.shading.dpdu = self.apply_to_vector(si.shading.dpdu); - ret.shading.dpdv = self.apply_to_vector(si.shading.dpdv); - ret.shading.dndu = self.apply_to_normal(si.shading.dndu); - ret.shading.dndv = self.apply_to_normal(si.shading.dndv); + let n = self.apply_to_normal(si.common.n); + ret.common.wo = self.apply_to_vector(si.common.wo).normalize(); - ret.common.n = n.normalize().face_forward(shading_n.into()); + ret.dpdu = self.apply_to_vector(si.dpdu); + ret.dpdv = self.apply_to_vector(si.dpdv); + ret.dndu = self.apply_to_normal(si.dndu); + ret.dndv = self.apply_to_normal(si.dndv); - Box::new(ret) - } else if let Some(mi) = inter.as_any().downcast_ref::() { - let ret = MediumInteraction { - common: InteractionData { - pi: self.apply_to_interval(&mi.common.pi), - n: self.apply_to_normal(mi.common.n).normalize(), - wo: self.apply_to_vector(mi.common.wo).normalize(), - time: mi.common.time, - medium_interface: mi.common.medium_interface.clone(), - medium: mi.common.medium.clone(), - }, - medium: mi.medium.clone(), - phase: mi.phase.clone(), - }; - Box::new(ret) - } else { - panic!("Unhandled Interaction type in Transform::apply_to_interaction"); + ret.dpdx = self.apply_to_vector(si.dpdx); + ret.dpdy = self.apply_to_vector(si.dpdy); + + let shading_n = self.apply_to_normal(si.shading.n); + ret.shading.n = shading_n.normalize(); + ret.shading.dpdu = self.apply_to_vector(si.shading.dpdu); + ret.shading.dpdv = self.apply_to_vector(si.shading.dpdv); + ret.shading.dndu = self.apply_to_normal(si.shading.dndu); + ret.shading.dndv = self.apply_to_normal(si.shading.dndv); + + ret.common.n = n.normalize().face_forward(ret.shading.n.into()); + + Interaction::Surface(ret) + } + + Interaction::Medium(mi) => { + let mut ret = mi.clone(); + + ret.common.pi = self.apply_to_interval(&mi.common.pi); + ret.common.n = self.apply_to_normal(mi.common.n).normalize(); + ret.common.wo = self.apply_to_vector(mi.common.wo).normalize(); + + Interaction::Medium(ret) + } + + Interaction::Simple(sim) => { + let mut ret = sim.clone(); + + ret.common.pi = self.apply_to_interval(&sim.common.pi); + + if sim.common.n != Default::default() { + ret.common.n = self.apply_to_normal(sim.common.n).normalize(); + } + if sim.common.wo != Default::default() { + ret.common.wo = self.apply_to_vector(sim.common.wo).normalize(); + } + + Interaction::Simple(ret) + } } } @@ -416,6 +440,18 @@ impl Transform { Self { m, m_inv } } + pub fn perspective(fov: Float, n: Float, f: Float) -> Result, InversionError> { + let persp: SquareMatrix = SquareMatrix::new([ + [1., 0., 0., 0.], + [0., 1., 0., 0.], + [0., 0., f / (f - n), -f * n / (f - n)], + [0., 0., 1., 0.], + ]); + let inv_tan_ang = 1. / (radians(fov) / 2.).tan(); + let persp_transform = Transform::from_matrix(persp)?; + Ok(Transform::scale(inv_tan_ang, inv_tan_ang, 1.) * persp_transform) + } + pub fn orthographic(z_near: Float, z_far: Float) -> Self { Self::scale(1., 1., 1. / (z_far - z_near)) * Self::translate(Vector3f::new(0., 0., -z_near)) } @@ -557,6 +593,20 @@ impl Transform { * m; (t, r, s) } + + pub fn has_scale(&self, tolerance: Option) -> bool { + let la2 = self + .apply_to_vector(Vector3f::new(1., 0., 0.)) + .norm_squared(); + let lb2 = self + .apply_to_vector(Vector3f::new(0., 1., 0.)) + .norm_squared(); + let lc2 = self + .apply_to_vector(Vector3f::new(0., 0., 1.)) + .norm_squared(); + let tol = tolerance.unwrap_or(1e-3); + (la2 - 1.).abs() > tol || (lb2 - 1.).abs() > tol || (lc2 - 1.).abs() > tol + } } impl PartialEq for Transform { @@ -1935,7 +1985,7 @@ impl AnimatedTransform { t.apply_to_ray(r, t_max) } - pub fn apply_interaction(&self, si: &dyn Interaction) -> Box { + pub fn apply_interaction(&self, si: &Interaction) -> Interaction { if !self.actually_animated { return self.start_transform.apply_to_interaction(si); }