From 4bb94877022d96b4a4585d1da79cd8b6ce47cddb Mon Sep 17 00:00:00 2001 From: Wentao Liu Date: Tue, 26 Dec 2023 10:55:10 +0800 Subject: [PATCH] feature: implement SimplePath integrator (sample lights) --- src/base/integrator.rs | 22 +++ src/base/interaction.rs | 12 +- src/base/shape.rs | 28 +++ src/base/texture.rs | 2 +- src/bin/main.rs | 2 +- src/euclidean_space/normal.rs | 53 ++++- src/euclidean_space/transform.rs | 2 +- src/euclidean_space/vector3.rs | 10 + src/integrators/ambient_occlusion.rs | 7 +- src/integrators/random_walk.rs | 7 +- src/integrators/simple_path.rs | 17 +- src/integrators/surface_normal.rs | 8 +- src/lib.rs | 1 + src/light_samplers/uniform_light_sampler.rs | 2 +- src/lights/diffuse_area.rs | 44 ++++- src/pbrt.rs | 10 +- src/scene/scene_builder.rs | 34 +++- src/shapes/sphere.rs | 20 +- src/shapes/triangle.rs | 204 +++++++++++++++++++- src/shapes/triangle_mesh.rs | 8 +- src/spectra/rgb_to_spectrum_table.rs | 2 +- src/util/math.rs | 37 ++-- src/util/sampling.rs | 161 +++++++++++++++ 23 files changed, 613 insertions(+), 80 deletions(-) diff --git a/src/base/integrator.rs b/src/base/integrator.rs index d7812cea..76a8efec 100644 --- a/src/base/integrator.rs +++ b/src/base/integrator.rs @@ -7,6 +7,28 @@ pub struct IntegratorBase { pub infinite_lights: Vec>, } +impl IntegratorBase { + pub fn new( + aggregate: Arc, + camera: Arc, + lights: Vec>, + ) -> Self { + let mut infinite_lights = vec![]; + for _light in &lights { + if _light.light_type() == LightType::Infinite { + infinite_lights.push(_light.clone()); + } + } + + return Self { + camera, + aggregate, + lights, + infinite_lights, + }; + } +} + pub trait Integrator: Send + Sync { fn fast_intersect(&self, ray: &Ray, t_max: f64) -> bool; diff --git a/src/base/interaction.rs b/src/base/interaction.rs index c0f4bc28..46dd25c0 100644 --- a/src/base/interaction.rs +++ b/src/base/interaction.rs @@ -127,11 +127,11 @@ impl SurfaceInteraction { // Compute auxiliary intersection points with plane, _px_ and _py_ let p = Point3f::from(self.interaction.pi); - let d = -n.dot(Vector3f::from(p)); - let tx = (-n.dot(Vector3::from(ray.rx_origin)) - d) / n.dot(ray.rx_direction); + let d = -n.dot(p.into()); + let tx = (-n.dot(ray.rx_origin.into()) - d) / n.dot(ray.rx_direction); let px = ray.rx_origin + tx * ray.rx_direction; - let ty = (-n.dot(Vector3f::from(ray.ry_origin)) - d) / n.dot(ray.ry_direction); + let ty = (-n.dot(ray.ry_origin.into()) - d) / n.dot(ray.ry_direction); let py = ray.ry_origin + ty * ray.ry_direction; @@ -140,7 +140,7 @@ impl SurfaceInteraction { } else { // Approximate screen-space change in $\pt{}$ based on camera projection (self.dpdx, self.dpdy) = camera.approximate_dp_dxy( - Point3f::from(self.interaction.pi), + self.interaction.pi.into(), self.interaction.n, samples_per_pixel, ); @@ -176,7 +176,7 @@ impl SurfaceInteraction { // Clamp derivatives of $u$ and $v$ to reasonable values let local_clamp = |x: f64| { if x.is_finite() { - clamp_float(x, -1e8, 1e8) + x.clamp(-1e8, 1e8) } else { 0.0 } @@ -223,7 +223,7 @@ impl SurfaceInteraction { return match &self.area_light { None => SampledSpectrum::same_value(0.0), Some(are_light) => are_light.l( - Point3f::from(self.interaction.pi), + self.interaction.pi.into(), self.interaction.n, self.interaction.uv, w, diff --git a/src/base/shape.rs b/src/base/shape.rs index 36c18ce8..7b1df8e2 100644 --- a/src/base/shape.rs +++ b/src/base/shape.rs @@ -1,5 +1,8 @@ use crate::pbrt::*; +pub const MIN_SPHERICAL_SAMPLE_AREA: f64 = 3e-4; +pub const MAX_SPHERICAL_SAMPLE_AREA: f64 = 6.22; + pub struct Shading { pub n: Normal3f, pub dpdu: Vector3f, @@ -20,6 +23,27 @@ impl Shading { } } +pub struct ShapeSampleContext { + pub pi: Point3fi, + pub n: Normal3f, + pub ns: Normal3f, +} + +pub struct ShapeSample { + pub interaction: Interaction, + pub pdf: f64, +} + +impl ShapeSampleContext { + pub fn offset_ray_origin(&self, w: Vector3f) -> Point3f { + panic!("ShapeSampleContext::offset_ray_origin() not implemented"); + } + + pub fn spawn_ray(&self, w: Vector3f) -> Ray { + panic!("ShapeSampleContext::spawn_ray() not implemented"); + } +} + pub struct ShapeIntersection { pub t_hit: f64, pub surface_interaction: SurfaceInteraction, @@ -39,4 +63,8 @@ pub trait Shape: Send + Sync { fn bounds(&self) -> Bounds3f; fn area(&self) -> f64; + + fn sample(&self, u: Point2f) -> Option; + + fn sample_with_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option; } diff --git a/src/base/texture.rs b/src/base/texture.rs index dca98a90..269e1eac 100644 --- a/src/base/texture.rs +++ b/src/base/texture.rs @@ -23,7 +23,7 @@ impl Display for TextureEvalContext { impl TextureEvalContext { pub fn new(si: &SurfaceInteraction) -> Self { return Self { - p: Point3f::from(si.interaction.pi), + p: si.interaction.pi.into(), dpdx: si.dpdx, dpdy: si.dpdy, n: si.interaction.n, diff --git a/src/bin/main.rs b/src/bin/main.rs index c7a3610e..4251562f 100644 --- a/src/bin/main.rs +++ b/src/bin/main.rs @@ -53,7 +53,7 @@ fn main() { let absolute_path = fs::canonicalize(args.scene_file).unwrap(); - let spp = args.spp.unwrap_or_else(|| 16); + let spp = args.spp.unwrap_or_else(|| 32); render(&absolute_path.display().to_string(), spp); } diff --git a/src/euclidean_space/normal.rs b/src/euclidean_space/normal.rs index 683289f6..8e5c350e 100644 --- a/src/euclidean_space/normal.rs +++ b/src/euclidean_space/normal.rs @@ -1,6 +1,6 @@ use crate::pbrt::*; -#[derive(Copy, Clone)] +#[derive(Copy, Clone, PartialEq)] pub struct Normal3f { pub x: f64, pub y: f64, @@ -20,6 +20,13 @@ impl Normal3f { }; } + pub fn is_valid(&self) -> bool { + return self.x.is_finite() + && self.y.is_finite() + && self.z.is_finite() + && (self.x * self.y * self.z != 0.0); + } + pub fn normalize(&self) -> Normal3f { let length = (sqr(self.x) + sqr(self.y) + sqr(self.z)).sqrt(); return Normal3f { @@ -75,3 +82,47 @@ impl Neg for Normal3f { }; } } + +impl Add for Normal3f { + type Output = Normal3f; + + fn add(self, rhs: Normal3f) -> Self::Output { + return Normal3f { + x: self.x + rhs.x, + y: self.y + rhs.y, + z: self.z + rhs.z, + }; + } +} + +impl Mul for Normal3f { + type Output = Normal3f; + + fn mul(self, rhs: f64) -> Self::Output { + return Self { + x: self.x * rhs, + y: self.y * rhs, + z: self.z * rhs, + }; + } +} + +impl Mul for f64 { + type Output = Normal3f; + + fn mul(self, rhs: Normal3f) -> Self::Output { + return Normal3f { + x: self * rhs.x, + y: self * rhs.y, + z: self * rhs.z, + }; + } +} + +impl MulAssign for Normal3f { + fn mul_assign(&mut self, rhs: f64) { + self.x *= rhs; + self.y *= rhs; + self.z *= rhs; + } +} diff --git a/src/euclidean_space/transform.rs b/src/euclidean_space/transform.rs index 922325b2..ed912dbf 100644 --- a/src/euclidean_space/transform.rs +++ b/src/euclidean_space/transform.rs @@ -382,7 +382,7 @@ impl Transform { let dt = d.abs().dot(o.error()) / length_squared; let offset_o = o + Vector3fi::from(d * dt); - return (Ray::new(Point3f::from(offset_o), d), dt); + return (Ray::new(offset_o.into(), d), dt); } pub fn on_differential_ray(&self, r: &DifferentialRay) -> (DifferentialRay, f64) { diff --git a/src/euclidean_space/vector3.rs b/src/euclidean_space/vector3.rs index a483c3d3..31d97cea 100644 --- a/src/euclidean_space/vector3.rs +++ b/src/euclidean_space/vector3.rs @@ -179,6 +179,16 @@ impl Vector3f { z, }; } + + // Equivalent to std::acos(Dot(a, b)), but more numerically stable. + // via http://www.plunk.org/~hatch/rightway.html + pub fn angle_between(&self, v: Vector3f) -> f64 { + if self.dot(v) < 0.0 { + return PI - 2.0 * safe_asin((*self + v).length() / 2.0); + } + + return 2.0 * safe_asin((v - *self).length() / 2.0); + } } impl Vector3 { diff --git a/src/integrators/ambient_occlusion.rs b/src/integrators/ambient_occlusion.rs index bea871bc..dc56e1db 100644 --- a/src/integrators/ambient_occlusion.rs +++ b/src/integrators/ambient_occlusion.rs @@ -15,12 +15,7 @@ impl AmbientOcclusion { let illuminant_scale = 1.0 / illuminant_spectrum.to_photometric(); return Self { - base: IntegratorBase { - aggregate, - camera, - lights: vec![], - infinite_lights: vec![], - }, + base: IntegratorBase::new(aggregate, camera, vec![]), illuminant_spectrum, illuminant_scale, }; diff --git a/src/integrators/random_walk.rs b/src/integrators/random_walk.rs index 2b8ffc53..84815d5f 100644 --- a/src/integrators/random_walk.rs +++ b/src/integrators/random_walk.rs @@ -23,12 +23,7 @@ impl RandomWalkIntegrator { } return Self { - base: IntegratorBase { - aggregate, - camera, - lights, - infinite_lights, - }, + base: IntegratorBase::new(aggregate, camera, lights), illuminant_spectrum, illuminant_scale, }; diff --git a/src/integrators/simple_path.rs b/src/integrators/simple_path.rs index aa4c8d30..0f86ed6d 100644 --- a/src/integrators/simple_path.rs +++ b/src/integrators/simple_path.rs @@ -117,16 +117,19 @@ impl Integrator for SimplePath { } impl SimplePath { - pub fn new(base: IntegratorBase) -> Self { - unreachable!(); + pub fn new( + aggregate: Arc, + camera: Arc, + lights: Vec>, + ) -> Self { + let light_sampler = UniformLightSampler { + lights: lights.clone(), + }; - /* return Self { - base, - sample_bsdf: true, - sample_light: true, + base: IntegratorBase::new(aggregate, camera, lights), max_depth: 5, + light_sampler, }; - */ } } diff --git a/src/integrators/surface_normal.rs b/src/integrators/surface_normal.rs index 5066f7ff..9f8d5dde 100644 --- a/src/integrators/surface_normal.rs +++ b/src/integrators/surface_normal.rs @@ -11,14 +11,8 @@ impl SurfaceNormal { camera: Arc, color_space: &RGBColorSpace, ) -> Self { - let val = 0.01; return SurfaceNormal { - base: IntegratorBase { - aggregate, - camera, - lights: vec![], - infinite_lights: vec![], - }, + base: IntegratorBase::new(aggregate, camera, vec![]), rgb: color_space.generate_albedo_rgb(), }; } diff --git a/src/lib.rs b/src/lib.rs index 05b489d2..9144c6ed 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,5 +1,6 @@ #![feature(const_fmt_arguments_new)] #![feature(const_fn_floating_point_arithmetic)] +#![feature(const_float_bits_conv)] #![feature(const_float_classify)] #![feature(const_trait_impl)] #![feature(generic_const_exprs)] diff --git a/src/light_samplers/uniform_light_sampler.rs b/src/light_samplers/uniform_light_sampler.rs index 74d22bd3..3703b3fc 100644 --- a/src/light_samplers/uniform_light_sampler.rs +++ b/src/light_samplers/uniform_light_sampler.rs @@ -1,7 +1,7 @@ use crate::pbrt::*; pub struct UniformLightSampler { - lights: Vec>, + pub lights: Vec>, } impl LightSampler for UniformLightSampler { diff --git a/src/lights/diffuse_area.rs b/src/lights/diffuse_area.rs index 3710de04..d5d92503 100644 --- a/src/lights/diffuse_area.rs +++ b/src/lights/diffuse_area.rs @@ -1,4 +1,5 @@ use crate::pbrt::*; +use crate::FilterFunction::Point; pub struct DiffuseAreaLight { base: LightBase, @@ -42,8 +43,47 @@ impl Light for DiffuseAreaLight { lambda: &SampledWavelengths, allow_incomplete_pdf: bool, ) -> Option { - //TODO: progress 2023/12/20 implement DiffuseAreaLight - panic!("DiffuseAreaLight::sample_li() not implemented"); + // Sample point on shape for _DiffuseAreaLight_ + let shape_ctx = ShapeSampleContext { + pi: ctx.pi, + n: ctx.n, + ns: ctx.ns, + }; + + let ss = match self.shape.sample_with_context(&shape_ctx, u) { + None => { + return None; + } + Some(_ss) => _ss, + }; + + if ss.pdf == 0.0 + || (Point3f::from(ss.interaction.pi) - Point3f::from(ctx.pi)).length_squared() == 0.0 + { + return None; + } + + // Return _LightLiSample_ for sampled point on shape + let wi = (Point3f::from(ss.interaction.pi) - Point3f::from(ctx.pi)).normalize(); + + let Le = self.l( + ss.interaction.pi.into(), + ss.interaction.n, + ss.interaction.uv, + -wi, + lambda, + ); + + if !Le.is_positive() { + return None; + } + + return Some(LightLiSample { + l: Le, + wi, + pdf: ss.pdf, + p_light: ss.interaction, + }); } } diff --git a/src/pbrt.rs b/src/pbrt.rs index fa42511c..d79530dd 100644 --- a/src/pbrt.rs +++ b/src/pbrt.rs @@ -2,6 +2,11 @@ pub struct GlobalVariable { pub rgb_color_space: Arc, } +pub const ONE_MINUS_EPSILON: f64 = { + let bits_of_one = 1.0f64.to_bits(); + f64::from_bits(bits_of_one - 1) +}; + pub use fma::fma; pub use rand::{rngs::StdRng, thread_rng, Rng, SeedableRng}; pub use rayon::prelude::*; @@ -65,7 +70,6 @@ pub use crate::{ pub type Point2f = Point2; pub type Point2i = Point2; -// TODO: rewrite Point2i to Point2 usize pub type Point3f = Point3; pub type Point3fi = Point3; @@ -75,7 +79,3 @@ pub type Vector3fi = Vector3; pub type Bounds2f = Bounds2; pub type Bounds3f = Bounds3; - -pub fn same_type() -> bool { - return TypeId::of::() == TypeId::of::(); -} diff --git a/src/scene/scene_builder.rs b/src/scene/scene_builder.rs index e72fb4eb..35cba16d 100644 --- a/src/scene/scene_builder.rs +++ b/src/scene/scene_builder.rs @@ -127,6 +127,8 @@ fn build_integrator( lights: Vec>, global_variable: &GlobalVariable, ) -> Arc { + println!("Integrator: `{}`", name); + return match name { "ambientocclusion" => Arc::new(AmbientOcclusion::new( global_variable.rgb_color_space.illuminant, @@ -140,11 +142,15 @@ fn build_integrator( camera, lights, )), + + "simplepath" => Arc::new(SimplePath::new(aggregate, camera, lights)), + "surfacenormal" => Arc::new(SurfaceNormal::new( aggregate, camera, global_variable.rgb_color_space.as_ref(), )), + _ => { panic!("unrecognized integrator: `{}`", name); } @@ -491,10 +497,21 @@ impl SceneBuilder { debug_assert!(tokens.len() == 2); debug_assert!(tokens[0].clone() == Token::Keyword("Transform".to_string())); - let floats = tokens[1..] - .into_iter() - .map(|t| t.convert_to_float()) - .collect::>(); + let floats = match tokens[1].clone() { + Token::List(numbers) => { + assert_eq!(tokens.len(), 2); + + numbers + .into_iter() + .map(|x| x.parse::().unwrap()) + .collect::>() + } + _ => tokens[1..] + .into_iter() + .map(|t| t.convert_to_float()) + .collect::>(), + }; + debug_assert!(floats.len() == 16); self.graphics_state.current_transform = @@ -568,7 +585,7 @@ impl SceneBuilder { } "Integrator" => { - self.integrator_name = tokens[1].convert_to_string(); + self.option_integrator(tokens); } "LightSource" => { @@ -687,6 +704,13 @@ impl SceneBuilder { }); } + fn option_integrator(&mut self, tokens: &[Token]) { + debug_assert!(tokens[0].clone() == Token::Keyword("Integrator".to_string())); + + self.integrator_name = tokens[1].convert_to_string(); + //TODO: parse Integrator: "integer pixelsamples" [ 32 ] + } + fn option_look_at(&mut self, tokens: &[Token]) { debug_assert!(tokens[0].clone() == Token::Keyword("LookAt".to_string())); diff --git a/src/shapes/sphere.rs b/src/shapes/sphere.rs index 8bf15e96..910b6a3d 100644 --- a/src/shapes/sphere.rs +++ b/src/shapes/sphere.rs @@ -22,11 +22,13 @@ impl Sphere { z_max: f64, phi_max: f64, ) -> Self { - let z_min = clamp_float(z_min, -radius, radius); - let z_max = clamp_float(z_max, -radius, radius); - let theta_z_min = clamp_float(z_min.min(z_max) / radius, -1.0, 1.0).acos(); - let theta_z_max = clamp_float(z_min.max(z_max) / radius, -1.0, 1.0).acos(); - let phi_max = degree_to_radian(clamp_float(phi_max, 0.0, 360.0)); + let z_min = z_min.clamp(-radius, radius); + let z_max = z_max.clamp(-radius, radius); + + let theta_z_min = (z_min.min(z_max) / radius).clamp(-1.0, 1.0).acos(); + let theta_z_max = (z_min.max(z_max) / radius).clamp(-1.0, 1.0).acos(); + + let phi_max = degree_to_radian(phi_max.clamp(0.0, 360.0)); return Sphere { render_from_object, @@ -250,4 +252,12 @@ impl Shape for Sphere { let r = self.radius; return 4.0 * PI * r * r; } + + fn sample(&self, u: Point2f) -> Option { + panic!("Sphere::sample() not implemented"); + } + + fn sample_with_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + panic!("Sphere::sample_with_context() not implemented"); + } } diff --git a/src/shapes/triangle.rs b/src/shapes/triangle.rs index 84e83f71..8199545f 100644 --- a/src/shapes/triangle.rs +++ b/src/shapes/triangle.rs @@ -53,9 +53,9 @@ fn intersect_triangle( p2t.x += sx * p2t.z; p2t.y += sy * p2t.z; - let mut e0 = difference_of_products(p1t.x, p2t.y, p1t.y, p2t.x); - let mut e1 = difference_of_products(p2t.x, p0t.y, p2t.y, p0t.x); - let mut e2 = difference_of_products(p0t.x, p1t.y, p0t.y, p1t.x); + let e0 = difference_of_products(p1t.x, p2t.y, p1t.y, p2t.x); + let e1 = difference_of_products(p2t.x, p0t.y, p2t.y, p0t.x); + let e2 = difference_of_products(p0t.x, p1t.y, p0t.y, p1t.x); if (e0 < 0.0 || e1 < 0.0 || e2 < 0.0) && (e0 > 0.0 || e1 > 0.0 || e2 > 0.0) { return None; @@ -133,10 +133,16 @@ impl Triangle { let v1 = self.mesh.indices[self.idx + 1]; let v2 = self.mesh.indices[self.idx + 2]; - return ( - self.mesh.points[v0], - self.mesh.points[v1], - self.mesh.points[v2], + return (self.mesh.p[v0], self.mesh.p[v1], self.mesh.p[v2]); + } + + fn solid_angle(&self, p: Point3f) -> f64 { + let (p0, p1, p2) = self.get_points(); + + return spherical_triangle_area( + (p0 - p).normalize(), + (p1 - p).normalize(), + (p2 - p).normalize(), ); } @@ -153,9 +159,9 @@ impl Triangle { let v1 = self.mesh.indices[self.idx + 1]; let v2 = self.mesh.indices[self.idx + 2]; - let p0 = self.mesh.points[v0]; - let p1 = self.mesh.points[v1]; - let p2 = self.mesh.points[v2]; + let p0 = self.mesh.p[v0]; + let p1 = self.mesh.p[v1]; + let p2 = self.mesh.p[v2]; let uv = if self.mesh.uv.len() > 0 { let _uv = &self.mesh.uv; @@ -262,4 +268,182 @@ impl Shape for Triangle { let (p0, p1, p2) = self.get_points(); return (p1 - p0).cross(p2 - p0).length() * 0.5; } + + fn sample(&self, u: Point2f) -> Option { + // Get triangle vertices in _p0_, _p1_, and _p2_ + + let v = [ + self.mesh.indices[self.idx + 0], + self.mesh.indices[self.idx + 1], + self.mesh.indices[self.idx + 2], + ]; + + let p0 = self.mesh.p[v[0]]; + let p1 = self.mesh.p[v[1]]; + let p2 = self.mesh.p[v[2]]; + + // Sample point on triangle uniformly by area + let b = sample_uniform_triangle(u); + let p = b[0] * p0 + b[1] * p1 + b[2] * p2; + + // Compute surface normal for sampled point on triangle + let n = { + let n = Normal3f::from((p1 - p0).cross(p2 - p0)).normalize(); + if self.mesh.n.len() > 0 { + let ns = b[0] * self.mesh.n[v[0]] + + b[1] * self.mesh.n[v[1]] + + (1.0 - b[0] - b[1]) * self.mesh.n[v[2]]; + + n.face_forward(ns.into()) + } else if self.mesh.reverse_orientation { + -n + } else { + n + } + }; + + // Compute $(u,v)$ for sampled point on triangle + // Get triangle texture coordinates in _uv_ array + let uv = if self.mesh.uv.len() > 0 { + [self.mesh.uv[v[0]], self.mesh.uv[v[1]], self.mesh.uv[v[2]]] + } else { + [ + Point2f::new(0.0, 0.0), + Point2f::new(1.0, 0.0), + Point2f::new(1.0, 1.0), + ] + }; + + let uv_sample = b[0] * uv[0] + b[1] * uv[1] + b[2] * uv[2]; + // Compute error bounds _pError_ for sampled point on triangle + + let p_abs_sum = (b[0] * p0).abs() + (b[1] * p1).abs() + ((1.0 - b[0] - b[1]) * p2).abs(); + let p_error = gamma(6) * p_abs_sum; + + let shape_sample = ShapeSample { + interaction: Interaction { + pi: Point3fi::from_value_and_error(p, p_error.into()), + n, + wo: Vector3f::nan(), + uv: uv_sample, + }, + pdf: 1.0 / self.area(), + }; + + return Some(shape_sample); + } + + fn sample_with_context(&self, ctx: &ShapeSampleContext, u: Point2f) -> Option { + let v = [ + self.mesh.indices[self.idx + 0], + self.mesh.indices[self.idx + 1], + self.mesh.indices[self.idx + 2], + ]; + let (p0, p1, p2) = (self.mesh.p[v[0]], self.mesh.p[v[1]], self.mesh.p[v[2]]); + + // Use uniform area sampling for numerically unstable cases + let solid_angle = self.solid_angle(ctx.pi.into()); + if solid_angle < MIN_SPHERICAL_SAMPLE_AREA || solid_angle > MAX_SPHERICAL_SAMPLE_AREA { + // Sample shape by area and compute incident direction _wi_ + let mut ss = self.sample(u).unwrap(); + let _wi = Point3f::from(ss.interaction.pi) - Point3f::from(ctx.pi); + + if _wi.length_squared() == 0.0 { + return None; + } + let wi = _wi.normalize(); + + ss.pdf /= ss.interaction.n.dot(-wi).abs() + / (Point3f::from(ctx.pi) - Point3f::from(ss.interaction.pi)).length_squared(); + if ss.pdf.is_infinite() { + return None; + } + + return Some(ss); + } + + // Sample spherical triangle from reference point + // Apply warp product sampling for cosine factor at reference point + + let mut pdf = 1.0; + let mut u = u; + if ctx.ns.is_valid() { + // This part is slightly different with PBRT-v4 + // Compute $\cos\theta$-based weights _w_ at sample domain corners + + let rp = Point3f::from(ctx.pi); + let wi = [ + (p0 - rp).normalize(), + (p1 - rp).normalize(), + (p2 - rp).normalize(), + ]; + + let w = [ + ctx.ns.dot(wi[1]).abs().max(0.01), + ctx.ns.dot(wi[1]).abs().max(0.01), + ctx.ns.dot(wi[0]).abs().max(0.01), + ctx.ns.dot(wi[2]).abs().max(0.01), + ]; + + u = sample_bilinear(u, &w); + pdf = bilinear_pdf(u, &w); + } + + let (b, tri_pdf) = sample_spherical_triangle(&[p0, p1, p2], ctx.pi.into(), u); + if tri_pdf == 0.0 { + return None; + } + pdf *= tri_pdf; + + // Compute error bounds _pError_ for sampled point on triangle + let p_abs_sum = (b[0] * p0).abs() + (b[1] * p1).abs() + ((1.0 - b[0] - b[1]) * p2).abs(); + let p_error = Vector3f::from(gamma(6) * p_abs_sum); + + // Return _ShapeSample_ for solid angle sampled point on triangle + let p = b[0] * p0 + b[1] * p1 + b[2] * p2; + + // Compute surface normal for sampled point on triangle + + let n = { + let n = Normal3f::from((p1 - p0).cross(p2 - p0)).normalize(); + + if self.mesh.n.len() > 0 { + let ns = b[0] * self.mesh.n[v[0]] + + b[1] * self.mesh.n[v[1]] + + (1.0 - b[0] - b[1]) * self.mesh.n[v[2]]; + + n.face_forward(ns.into()) + } else if self.mesh.reverse_orientation { + n * -1.0 + } else { + n + } + }; + + // Compute $(u,v)$ for sampled point on triangle + // Get triangle texture coordinates in _uv_ array + + let uv = if self.mesh.uv.len() > 0 { + [self.mesh.uv[v[0]], self.mesh.uv[v[1]], self.mesh.uv[v[2]]] + } else { + [ + Point2f::new(0.0, 0.0), + Point2f::new(1.0, 0.0), + Point2f::new(1.0, 1.0), + ] + }; + + let uv_sample = b[0] * uv[0] + b[1] * uv[1] + b[2] * uv[2]; + let shape_sample = ShapeSample { + interaction: Interaction { + pi: Point3fi::from_value_and_error(p, p_error), + n, + wo: Vector3::nan(), + uv: uv_sample, + }, + pdf, + }; + + return Some(shape_sample); + } } diff --git a/src/shapes/triangle_mesh.rs b/src/shapes/triangle_mesh.rs index 157fd6de..b7edabfe 100644 --- a/src/shapes/triangle_mesh.rs +++ b/src/shapes/triangle_mesh.rs @@ -3,8 +3,8 @@ use crate::pbrt::*; pub struct TriangleMesh { pub reverse_orientation: bool, pub indices: Vec, - pub points: Vec, - pub normals: Vec, + pub p: Vec, + pub n: Vec, pub uv: Vec, } @@ -38,9 +38,9 @@ impl TriangleMesh { return TriangleMesh { reverse_orientation, - points: transformed_points, + p: transformed_points, indices, - normals: transformed_normals, + n: transformed_normals, uv, }; } diff --git a/src/spectra/rgb_to_spectrum_table.rs b/src/spectra/rgb_to_spectrum_table.rs index e790290b..324ef50b 100644 --- a/src/spectra/rgb_to_spectrum_table.rs +++ b/src/spectra/rgb_to_spectrum_table.rs @@ -74,7 +74,7 @@ impl RGBtoSpectrumTable { }; } - return clamp_usize(first - 1, 0, sz - 2); + return (first - 1).clamp(0, sz - 2); } pub fn eval(&self, rgb: RGB) -> RGBSigmoidPolynomial { diff --git a/src/util/math.rs b/src/util/math.rs index bd50a65a..310c456a 100644 --- a/src/util/math.rs +++ b/src/util/math.rs @@ -21,16 +21,6 @@ pub fn mod_i32(a: i32, b: i32) -> i32 { return if result < 0 { result + b } else { result }; } -pub const fn clamp_float(val: f64, low: f64, high: f64) -> f64 { - if val < low { - return low; - } - if val > high { - return high; - } - return val; -} - pub const fn clamp_usize(val: usize, low: usize, high: usize) -> usize { if val < low { return low; @@ -84,6 +74,13 @@ pub fn difference_of_products_vec3(a: f64, b: Vector3f, c: f64, d: Vector3f) -> return difference_of_prod + error; } +pub fn sum_of_products(a: f64, b: f64, c: f64, d: f64) -> f64 { + let cd = c * d; + let sum = fma(a, b, cd); + let error = fma(c, d, -cd); + return sum + error; +} + pub const fn is_power_of_2(v: i32) -> bool { return v > 0 && !(v & (v - 1) > 0); } @@ -99,10 +96,20 @@ pub const fn round_up_pow_2(v: i32) -> i32 { return x + 1; } +pub fn safe_asin(x: f64) -> f64 { + debug_assert!(x >= -1.0001 && x <= 1.0001); + + return x.clamp(-1.0, 1.0).asin(); +} + pub fn safe_acos(x: f64) -> f64 { debug_assert!(x >= -1.0001 && x <= 1.0001); - return clamp_float(x, -1.0, 1.0).acos(); + return x.clamp(-1.0, 1.0).acos(); +} + +pub fn gram_schmidt(v: Vector3f, w: Vector3f) -> Vector3f { + return v - v.dot(w) * w; } pub fn sqr + Copy>(x: T) -> T { @@ -133,3 +140,11 @@ pub fn windowed_sinc(x: f64, radius: f64, tau: f64) -> f64 { } return sinc(x) * sinc(x / tau); } + +pub fn spherical_triangle_area(a: Vector3f, b: Vector3f, c: Vector3f) -> f64 { + return a + .dot(b.cross(c)) + .atan2(1.0 + a.dot(b) + a.dot(c) + b.dot(c)) + .abs() + * 2.0; +} diff --git a/src/util/sampling.rs b/src/util/sampling.rs index 150c8d70..562c79be 100644 --- a/src/util/sampling.rs +++ b/src/util/sampling.rs @@ -1,5 +1,35 @@ use crate::pbrt::*; +/* +PBRT_CPU_GPU inline pstd::array SampleUniformTriangle(Point2f u) { + Float b0, b1; + if (u[0] < u[1]) { + b0 = u[0] / 2; + b1 = u[1] - b0; + } else { + b1 = u[1] / 2; + b0 = u[0] - b1; + } + return {b0, b1, 1 - b0 - b1}; +} +*/ + +pub fn sample_uniform_triangle(u: Point2f) -> [f64; 3] { + let (b0, b1) = { + if u[0] < u[1] { + let b0 = u[0] / 2.0; + let b1 = u[1] - b0; + (b0, b1) + } else { + let b1 = u[1] / 2.0; + let b0 = u[0] - b1; + (b0, b1) + } + }; + + return [b0, b1, 1.0 - b0 - b1]; +} + pub fn sample_uniform_disk_concentric(u: Point2f) -> Point2f { // Map _u_ to $[-1,1]^2$ and handle degeneracy at the origin @@ -43,3 +73,134 @@ pub fn visible_wavelengths_pdf(lambda: f64) -> f64 { return 0.0039398042 / sqr((0.0072 * (lambda - 538.0)).cosh()); } + +pub fn sample_linear(u: f64, a: f64, b: f64) -> f64 { + if u == 0.0 && a == 0.0 { + return 0.0; + } + + let x = u * (a + b) / (a + lerp(u, a * a, b * b).sqrt()); + return x.min(ONE_MINUS_EPSILON); +} + +pub fn sample_bilinear(u: Point2f, w: &[f64; 4]) -> Point2f { + let y = sample_linear(u[1], w[0] + w[1], w[2] + w[3]); + let x = sample_linear(u[0], lerp(y, w[0], w[2]), lerp(y, w[1], w[3])); + + return Point2f::new(x, y); +} + +pub fn bilinear_pdf(p: Point2f, w: &[f64; 4]) -> f64 { + let sum = w[0] + w[1] + w[2] + w[3]; + + if sum == 0.0 { + return 1.0; + } + + return 4.0 + * ((1.0 - p[0]) * (1.0 - p[1]) * w[0] + + p[0] * (1.0 - p[1]) * w[1] + + (1.0 - p[0]) * p[1] * w[2] + + p[0] * p[1] * w[3]) + / sum; +} + +pub fn sample_spherical_triangle(v: &[Point3f; 3], p: Point3f, u: Point2f) -> ([f64; 3], f64) { + // Compute vectors _a_, _b_, and _c_ to spherical triangle vertices + + let a = (v[0] - p).normalize(); + let b = (v[1] - p).normalize(); + let c = (v[2] - p).normalize(); + + // Compute normalized cross products of all direction pairs + let n_ab = a.cross(b); + let n_bc = b.cross(c); + let n_ca = c.cross(a); + + if n_ab.length_squared() == 0.0 || n_bc.length_squared() == 0.0 || n_ca.length_squared() == 0.0 + { + return ([f64::NAN, f64::NAN, f64::NAN], 0.0); + } + + let n_ab = n_ab.normalize(); + let n_bc = n_bc.normalize(); + let n_ca = n_ca.normalize(); + + // Find angles $\alpha$, $\beta$, and $\gamma$ at spherical triangle vertices + let alpha = n_ab.angle_between(-n_ca); + let beta = n_bc.angle_between(-n_ab); + let gamma = n_ca.angle_between(-n_bc); + + // Uniformly sample triangle area $A$ to compute $A'$ + let A_pi = alpha + beta + gamma; + let Ap_pi = lerp(u[0], PI, A_pi); + + let pdf = { + let A = A_pi - PI; + if A <= 0.0 { + 0.0 + } else { + 1.0 / A + } + }; + + // Find $\cos\beta'$ for point along _b_ for sampled area + let cos_alpha = alpha.cos(); + let sin_alpha = alpha.sin(); + let sin_phi = Ap_pi.sin() * cos_alpha - Ap_pi.cos() * sin_alpha; + let cos_phi = Ap_pi.cos() * cos_alpha + Ap_pi.sin() * sin_alpha; + + let k1 = cos_phi + cos_alpha; + let k2 = sin_phi - sin_alpha * a.dot(b); + + let cos_bp = (k2 + difference_of_products(k2, cos_phi, k1, sin_phi) * cos_alpha) + / ((sum_of_products(k2, sin_phi, k1, cos_phi)) * sin_alpha); + + // Happens if the triangle basically covers the entire hemisphere. + // We currently depend on calling code to detect this case, which + // is sort of ugly/unfortunate. + let cos_bp = cos_bp.clamp(-1.0, 1.0); + + // Sample $c'$ along the arc between $b'$ and $a$ + let sin_bp = safe_sqrt(1.0 - sqr(cos_bp)); + let cp = cos_bp * a + sin_bp * gram_schmidt(c, a).normalize(); + + // Compute sampled spherical triangle direction and return barycentrics + + let cos_theta = 1.0 - u[1] * (1.0 - cp.dot(b)); + let sin_theta = safe_sqrt(1.0 - sqr(cos_theta)); + + let w = cos_theta * b + sin_theta * gram_schmidt(cp, b).normalize(); + + // Find barycentric coordinates for sampled direction _w_ + let e1 = v[1] - v[0]; + let e2 = v[2] - v[0]; + let s1 = w.cross(e2); + let divisor = s1.dot(e1); + + if divisor == 0.0 { + // This happens with triangles that cover (nearly) the whole + // hemisphere. + return ([1.0 / 3.0; 3], pdf); + } + + let inv_divisor = 1.0 / divisor; + let s = p - v[0]; + let b1 = s.dot(s1) * inv_divisor; + let b2 = w.dot(s.cross(e1)) * inv_divisor; + + // Return clamped barycentrics for sampled direction + let (b1, b2) = { + let _b1 = b1.clamp(0.0, 1.0); + let _b2 = b2.clamp(0.0, 1.0); + + let sum = _b1 + _b2; + if sum > 1.0 { + (_b1 / sum, _b2 / sum) + } else { + (_b1, _b2) + } + }; + + return ([1.0 - b1 - b2, b1, b2], pdf); +}