Fixing precision issues on vector operations, disk intersection

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<M = Self, A = Self> {
     type Output;
@@ -23,6 +23,14 @@ impl MulAdd<Float, Float> for Float {
     }
 }
 
+impl MulAdd<f64, f64> 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<T, const N: usize> num_traits::Zero for $Struct<T, N>
+        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<T: Copy> Vector4<T> {
 // Vector operations
 impl<T> Vector3<T>
 where
-    T: Num + Copy + Neg<Output = T>,
+    T: Num + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn cross(self, rhs: Self) -> Self {
         Self([
-            self[1] * rhs[2] - self[2] * rhs[1],
+            difference_of_products(self[1], rhs[2], self[2], rhs[1]),
-            self[2] * rhs[0] - self[0] * rhs[2],
+            difference_of_products(self[2], rhs[0], self[0], rhs[2]),
-            self[0] * rhs[1] - self[1] * rhs[0],
+            difference_of_products(self[0], rhs[1], self[1], rhs[0]),
         ])
     }
 }
 
 impl<T> Normal3<T>
 where
-    T: Num + Copy + Neg<Output = T>,
+    T: Num + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn cross(self, rhs: Self) -> Self {
         Self([
-            self[1] * rhs[2] - self[2] * rhs[1],
+            difference_of_products(self[1], rhs[2], self[2], rhs[1]),
-            self[2] * rhs[0] - self[0] * rhs[2],
+            difference_of_products(self[2], rhs[0], self[0], rhs[2]),
-            self[0] * rhs[1] - self[1] * rhs[0],
+            difference_of_products(self[0], rhs[1], self[1], rhs[0]),
         ])
     }
 }
 
 impl<T> Vector3<T>
 where
-    T: Num + NumFloat + Copy + Neg<Output = T>,
+    T: Num + NumFloat + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn coordinate_system(&self) -> (Self, Self)
     where
@@ -663,7 +692,7 @@ where
 
 impl<T> Normal3<T>
 where
-    T: Num + NumFloat + Copy + Neg<Output = T>,
+    T: Num + NumFloat + Copy + Neg<Output = T> + Zero + MulAdd<T, T, Output = T>,
 {
     pub fn coordinate_system(&self) -> (Self, Self)
     where

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<QuadricIntersection> {
-        let oi = self.object_from_render.apply_to_interval(Point3fi::new_from_point(r.o));
+        let oi = self
-        let di = self.object_from_render.apply_to_vector_interval(Vector3fi::new_from_vector(r.d));
+            .object_from_render
-        // Reject disk intersections for rays parallel to the disk’s plane
+            .apply_to_interval(&Point3fi::new_from_point(r.o));
-        if di.z() == 0. {
+        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();
+        let oi_f = Point3f::from(oi);
-        if t_shape_hit == 0. || t_shape_hit >= t_max {
+        let di_f = Vector3f::from(di);
-            return None;
+        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,
         })

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;
         }

shared/src/utils/math.rs

@@ -71,15 +71,16 @@ pub fn evaluate_polynomial(t: Float, coeffs: &[Float]) -> Float {
     result
 }
 
-pub fn difference_of_products<T>(a: Float, b: T, c: Float, d: T) -> T
+pub fn difference_of_products<T, U, V>(a: T, b: U, c: T, d: U) -> V
 where
-    T: Copy + Neg<Output = T> + Mul<Float, Output = T> + Add<Output = T>,
+    T: Copy + Neg<Output = T>,
-    T: MulAdd<Float, T, Output = T>,
+    U: Copy + MulAdd<T, V, Output = V>,
+    V: Copy + Zero + Neg<Output = V> + Add<Output = V>,
 {
-    let cd = d * c;
+    let cd: V = d.mul_add(c, V::zero());
-    let diff = b.mul_add(a, -cd);
+    let diff: V = b.mul_add(a, -cd);
-    let error = d.mul_add(-c, cd);
+    let err:  V = d.mul_add(-c, cd);
-    diff + error
+    diff + err
 }
 
 #[inline]