From ac7fdd748611bcf34366c9e6f3ae194329f3c83b Mon Sep 17 00:00:00 2001 From: Wito Wiala Date: Fri, 5 Jun 2026 16:04:19 +0100 Subject: [PATCH] Fixing precision issues on vector operations, disk intersection --- shared/src/core/geometry/primitives.rs | 51 ++++++++++++++++++++------ shared/src/shapes/disk.rs | 32 +++++++++------- shared/src/shapes/sphere.rs | 8 ++-- shared/src/utils/math.rs | 15 ++++---- 4 files changed, 72 insertions(+), 34 deletions(-) diff --git a/shared/src/core/geometry/primitives.rs b/shared/src/core/geometry/primitives.rs index 88136e6..1daacc6 100644 --- a/shared/src/core/geometry/primitives.rs +++ b/shared/src/core/geometry/primitives.rs @@ -2,13 +2,13 @@ use super::traits::{SqrtExt, Tuple, VectorLike}; use super::{Float, NumFloat, PI}; use crate::utils::interval::Interval; use crate::utils::math::{clamp, difference_of_products, quadratic, safe_asin}; +use core::fmt; use core::hash::{Hash, Hasher}; use core::iter::Sum; use core::ops::{ Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign, }; use num_traits::{AsPrimitive, FloatConst, Num, Signed, Zero}; -use core::fmt; pub trait MulAdd { type Output; @@ -23,6 +23,14 @@ impl MulAdd for Float { } } +impl MulAdd for f64 { + type Output = f64; + #[inline(always)] + fn mul_add(self, multiplier: f64, addend: f64) -> Self::Output { + num_traits::Float::mul_add(self, multiplier, addend) + } +} + // N-dimensional displacement #[repr(C)] #[derive(Debug, Copy, Clone, PartialEq, Eq)] @@ -217,6 +225,27 @@ macro_rules! impl_tuple_core { }; } +#[macro_export] +macro_rules! impl_num_zero { + ($Struct:ident) => { + impl num_traits::Zero for $Struct + where + T: num_traits::Zero + Copy + PartialEq, + { + #[inline] + fn zero() -> Self { + Self([T::zero(); N]) + } + #[inline] + fn is_zero(&self) -> bool { + self.0.iter().all(|c| c.is_zero()) + } + } + }; +} +impl_num_zero!(Vector); +impl_num_zero!(Normal); + #[macro_export] macro_rules! impl_scalar_ops { ($Struct:ident) => { @@ -607,33 +636,33 @@ impl Vector4 { // Vector operations impl Vector3 where - T: Num + Copy + Neg, + T: Num + Copy + Neg + Zero + MulAdd, { pub fn cross(self, rhs: Self) -> Self { Self([ - self[1] * rhs[2] - self[2] * rhs[1], - self[2] * rhs[0] - self[0] * rhs[2], - self[0] * rhs[1] - self[1] * rhs[0], + difference_of_products(self[1], rhs[2], self[2], rhs[1]), + difference_of_products(self[2], rhs[0], self[0], rhs[2]), + difference_of_products(self[0], rhs[1], self[1], rhs[0]), ]) } } impl Normal3 where - T: Num + Copy + Neg, + T: Num + Copy + Neg + Zero + MulAdd, { pub fn cross(self, rhs: Self) -> Self { Self([ - self[1] * rhs[2] - self[2] * rhs[1], - self[2] * rhs[0] - self[0] * rhs[2], - self[0] * rhs[1] - self[1] * rhs[0], + difference_of_products(self[1], rhs[2], self[2], rhs[1]), + difference_of_products(self[2], rhs[0], self[0], rhs[2]), + difference_of_products(self[0], rhs[1], self[1], rhs[0]), ]) } } impl Vector3 where - T: Num + NumFloat + Copy + Neg, + T: Num + NumFloat + Copy + Neg + Zero + MulAdd, { pub fn coordinate_system(&self) -> (Self, Self) where @@ -663,7 +692,7 @@ where impl Normal3 where - T: Num + NumFloat + Copy + Neg, + T: Num + NumFloat + Copy + Neg + Zero + MulAdd, { pub fn coordinate_system(&self) -> (Self, Self) where diff --git a/shared/src/shapes/disk.rs b/shared/src/shapes/disk.rs index fe6c1f8..588a639 100644 --- a/shared/src/shapes/disk.rs +++ b/shared/src/shapes/disk.rs @@ -6,9 +6,10 @@ use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction use crate::core::shape::{ QuadricIntersection, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait, }; -use crate::utils::Transform; +use crate::utils::interval::Interval; use crate::utils::math::square; use crate::utils::sampling::sample_uniform_disk_concentric; +use crate::utils::Transform; use crate::{Float, PI}; use num_traits::Float as NumFloat; @@ -48,25 +49,30 @@ impl DiskShape { } fn basic_intersect(&self, r: &Ray, t_max: Float) -> Option { - let oi = self.object_from_render.apply_to_interval(Point3fi::new_from_point(r.o)); - let di = self.object_from_render.apply_to_vector_interval(Vector3fi::new_from_vector(r.d)); - // Reject disk intersections for rays parallel to the disk’s plane - if di.z() == 0. { + let oi = self + .object_from_render + .apply_to_interval(&Point3fi::new_from_point(r.o)); + let di = self + .object_from_render + .apply_to_vector_interval(&Vector3fi::new_from_vector(r.d)); + + if Float::from(di.z()) == 0. { + return None; + } + let t_shape_hit: Interval = (self.height - oi.z()) / di.z(); + if t_shape_hit.high <= 0. || t_shape_hit.low >= t_max { return None; } - let t_shape_hit = (self.height - oi.z()) / di.z(); - if t_shape_hit == 0. || t_shape_hit >= t_max { - return None; - } + let oi_f = Point3f::from(oi); + let di_f = Vector3f::from(di); + let t = Float::from(t_shape_hit); + let p_hit: Point3f = oi_f + di_f * t; - // See if hit point is inside disk radii and phi_max - let p_hit: Point3fi = oi + t_shape_hit * di; let dist2 = square(p_hit.x()) + square(p_hit.y()); if dist2 > square(self.radius) || dist2 < square(self.inner_radius) { return None; } - let mut phi = p_hit.y().atan2(p_hit.x()); if phi < 0. { phi += 2. * PI; @@ -76,7 +82,7 @@ impl DiskShape { } Some(QuadricIntersection { - t_hit: t_shape_hit, + t_hit: t, p_obj: p_hit, phi, }) diff --git a/shared/src/shapes/sphere.rs b/shared/src/shapes/sphere.rs index c7bba39..21d156e 100644 --- a/shared/src/shapes/sphere.rs +++ b/shared/src/shapes/sphere.rs @@ -1,17 +1,17 @@ +use crate::core::geometry::{spherical_direction, Frame, SqrtExt}; use crate::core::geometry::{ Bounds3f, DirectionCone, Normal3f, Point2f, Point3f, Point3fi, Ray, Vector2f, Vector3f, Vector3fi, VectorLike, }; -use crate::core::geometry::{Frame, SqrtExt, spherical_direction}; use crate::core::interaction::{Interaction, InteractionTrait, SurfaceInteraction}; use crate::core::pbrt::gamma; use crate::core::shape::{ QuadricIntersection, ShapeIntersection, ShapeSample, ShapeSampleContext, ShapeTrait, }; -use crate::utils::Transform; use crate::utils::interval::Interval; use crate::utils::math::{clamp, difference_of_products, radians, safe_acos, safe_sqrt, square}; use crate::utils::sampling::sample_uniform_sphere; +use crate::utils::Transform; use crate::{Float, PI}; use num_traits::Float as NumFloat; @@ -120,7 +120,9 @@ impl SphereShape { } let mut p_hit = Point3f::from(oi) + Float::from(t_shape_hit) * Vector3f::from(di); - p_hit *= self.radius / p_hit.distance(Point3f::new(0., 0., 0.)); + let scale = self.radius / p_hit.distance(Point3f::new(0., 0., 0.)); + p_hit = Point3f::from(Vector3f::from(p_hit) * scale); + if p_hit.x() == 0. && p_hit.y() == 0. { p_hit[0] = 1e-5 * self.radius; } diff --git a/shared/src/utils/math.rs b/shared/src/utils/math.rs index 53ab582..f906714 100644 --- a/shared/src/utils/math.rs +++ b/shared/src/utils/math.rs @@ -71,15 +71,16 @@ pub fn evaluate_polynomial(t: Float, coeffs: &[Float]) -> Float { result } -pub fn difference_of_products(a: Float, b: T, c: Float, d: T) -> T +pub fn difference_of_products(a: T, b: U, c: T, d: U) -> V where - T: Copy + Neg + Mul + Add, - T: MulAdd, + T: Copy + Neg, + U: Copy + MulAdd, + V: Copy + Zero + Neg + Add, { - let cd = d * c; - let diff = b.mul_add(a, -cd); - let error = d.mul_add(-c, cd); - diff + error + let cd: V = d.mul_add(c, V::zero()); + let diff: V = b.mul_add(a, -cd); + let err: V = d.mul_add(-c, cd); + diff + err } #[inline]