Merge pull request #149 from finnbear/bugfix_140_part_2

Fix #140 part 2 - detect NaN when intersection ray with AABB

CHANGELOG.md

@@ -24,6 +24,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/).
 - Fix intersection between rays and AABBs that lack depth. [#131](https://github.com/svenstaro/bvh/pull/131) (thanks @finnbear)
 - Document the limitations of distance traversal best-effort ordering. [#135](https://github.com/svenstaro/bvh/pull/135) (thanks @finnbear)
 - Add `Bvh::nearest_to`, which returns the nearest shape to a point. [#108](https://github.com/svenstaro/bvh/pull/108) (thanks @Azkellas)
+- Fix ray-AABB intersection such that a ray in the plane of an AABB face is never
+  considered intersecting, rather than returning an arbitrary answer in the case of
+  `Ray::intersects_aabb` or an erroneous answer in the case of
+  `Ray::intersection_slice_for_aabb`. [#149](https://github.com/svenstaro/bvh/pull/149) (thanks @finnbear)
 
 ## 0.10.0 - 2024-07-06
 - Don't panic when traversing empty BVH [#106](https://github.com/svenstaro/bvh/pull/106) (thanks @finnbear)

fuzz/fuzz_targets/fuzz.rs

@@ -383,13 +383,27 @@ impl<const D: usize> Workload<D> {
         let ray = self.ray.ray();
         let aabb = self.aabb.aabb();
 
-        // Make sure these functions agree and that the latter doesn't return NaN.
+        // Make sure these functions agree and that the latter doesn't return NaN/infinity.
         let intersects = ray.intersects_aabb(&aabb);
         if let Some((entry, exit)) = ray.intersection_slice_for_aabb(&aabb) {
-            assert!(intersects);
-            assert!(entry <= exit, "{entry} {exit}");
+            assert!(
+                intersects,
+                "intersection_slice_for_aabb returned Some for non-intersection"
+            );
+            assert!(entry >= 0.0, "entry ({entry}) must be non-negative");
+            assert!(
+                entry <= exit,
+                "entry ({entry}) must not exceed exit ({exit})"
+            );
+            assert!(
+                entry.is_finite() && exit.is_finite(),
+                "neither entry ({entry}) nor exit ({exit}) may be infinite or NaN for AABB's and rays with finite coordinates"
+            );
         } else {
-            assert!(!intersects);
+            assert!(
+                !intersects,
+                "intersection_slice_for_aabb returned None for intersection"
+            );
         }
 
         // If this environment variable is present, only test limited-size BVH's without mutations.

src/ray/intersect_default.rs

@@ -1,7 +1,11 @@
 //! This file contains the generic implementation of [`RayIntersection`]
 
 use super::Ray;
-use crate::{aabb::Aabb, bounding_hierarchy::BHValue, utils::fast_max};
+use crate::{
+    aabb::Aabb,
+    bounding_hierarchy::BHValue,
+    utils::{fast_max, has_nan},
+};
 
 /// The [`RayIntersection`] trait allows for generic implementation of ray intersection
 /// useful for our SIMD optimizations.
@@ -15,6 +19,14 @@ impl<T: BHValue, const D: usize> RayIntersection<T, D> for Ray<T, D> {
         let lbr = (aabb[0].coords - self.origin.coords).component_mul(&self.inv_direction);
         let rtr = (aabb[1].coords - self.origin.coords).component_mul(&self.inv_direction);
 
+        if has_nan(&lbr) | has_nan(&rtr) {
+            // Assumption: the ray is in the plane of an AABB face. Be consistent and
+            // consider this a non-intersection. This avoids making the result depend
+            // on which axis/axes have NaN (min/max in the code that follows are not
+            // commutative).
+            return false;
+        }
+
         let (inf, sup) = lbr.inf_sup(&rtr);
 
         let tmin = inf.max();
@@ -30,6 +42,10 @@ impl<T: BHValue, const D: usize> RayIntersection<T, D> for Ray<T, D> {
         let lbr = (aabb[0].coords - self.origin.coords).component_mul(&self.inv_direction);
         let rtr = (aabb[1].coords - self.origin.coords).component_mul(&self.inv_direction);
 
+        if has_nan(&lbr) | has_nan(&rtr) {
+            return false;
+        }
+
         let (inf, sup) = lbr.inf_sup(&rtr);
 
         let tmin = inf.max();

src/ray/intersect_x86_64.rs

@@ -5,7 +5,10 @@ use std::arch::x86_64::*;
 
 use nalgebra::SVector;
 
-use crate::{aabb::Aabb, utils::fast_max};
+use crate::{
+    aabb::Aabb,
+    utils::{fast_max, has_nan},
+};
 
 use super::{intersect_default::RayIntersection, Ray};
 
@@ -72,6 +75,33 @@ fn min_elem_m128(mm: __m128) -> f32 {
     }
 }
 
+#[inline(always)]
+fn has_nan_m128(mm: __m128) -> bool {
+    unsafe {
+        let mut data = std::mem::MaybeUninit::<[f32; 4]>::uninit();
+        _mm_store_ps(data.as_mut_ptr() as *mut f32, mm);
+        has_nan(data.assume_init_ref())
+    }
+}
+
+#[inline(always)]
+fn has_nan_m128d(mm: __m128d) -> bool {
+    unsafe {
+        let mut data = std::mem::MaybeUninit::<[f64; 2]>::uninit();
+        _mm_store_pd(data.as_mut_ptr() as *mut f64, mm);
+        has_nan(data.assume_init_ref())
+    }
+}
+
+#[inline(always)]
+fn has_nan_m256d(mm: __m256d) -> bool {
+    unsafe {
+        let mut data = std::mem::MaybeUninit::<[f64; 4]>::uninit();
+        _mm256_store_pd(data.as_mut_ptr() as *mut f64, mm);
+        has_nan(data.assume_init_ref())
+    }
+}
+
 #[inline(always)]
 fn ray_intersects_aabb_m128(
     ray_origin: __m128,
@@ -83,6 +113,10 @@ fn ray_intersects_aabb_m128(
         let v1 = _mm_mul_ps(_mm_sub_ps(aabb_0, ray_origin), ray_inv_dir);
         let v2 = _mm_mul_ps(_mm_sub_ps(aabb_1, ray_origin), ray_inv_dir);
 
+        if has_nan_m128(v1) | has_nan_m128(v2) {
+            return false;
+        }
+
         let inf = _mm_min_ps(v1, v2);
         let sup = _mm_max_ps(v1, v2);
 
@@ -176,6 +210,10 @@ fn ray_intersects_aabb_m128d(
         let v1 = _mm_mul_pd(_mm_sub_pd(aabb_0, ray_origin), ray_inv_dir);
         let v2 = _mm_mul_pd(_mm_sub_pd(aabb_1, ray_origin), ray_inv_dir);
 
+        if has_nan_m128d(v1) | has_nan_m128d(v2) {
+            return false;
+        }
+
         let inf = _mm_min_pd(v1, v2);
         let sup = _mm_max_pd(v1, v2);
 
@@ -257,6 +295,10 @@ fn ray_intersects_aabb_m256d(
         let v1 = _mm256_mul_pd(_mm256_sub_pd(aabb_0, ray_origin), ray_inv_dir);
         let v2 = _mm256_mul_pd(_mm256_sub_pd(aabb_1, ray_origin), ray_inv_dir);
 
+        if has_nan_m256d(v1) | has_nan_m256d(v2) {
+            return false;
+        }
+
         let inf = _mm256_min_pd(v1, v2);
         let sup = _mm256_max_pd(v1, v2);
 

src/ray/ray_impl.rs

@@ -4,7 +4,7 @@
 use std::cmp::Ordering;
 
 use crate::aabb::IntersectsAabb;
-use crate::utils::fast_max;
+use crate::utils::{fast_max, has_nan};
 use crate::{aabb::Aabb, bounding_hierarchy::BHValue};
 use nalgebra::{
     ClosedAddAssign, ClosedMulAssign, ClosedSubAssign, ComplexField, Point, SVector, SimdPartialOrd,
@@ -124,6 +124,14 @@ impl<T: BHValue, const D: usize> Ray<T, D> {
         let lbr = (aabb[0].coords - self.origin.coords).component_mul(&self.inv_direction);
         let rtr = (aabb[1].coords - self.origin.coords).component_mul(&self.inv_direction);
 
+        if has_nan(&lbr) | has_nan(&rtr) {
+            // Assumption: the ray is in the plane of an AABB face. Be consistent and
+            // consider this a non-intersection. This avoids making the result depend
+            // on which axis/axes have NaN (min/max in the code that follows are not
+            // commutative).
+            return None;
+        }
+
         let (inf, sup) = lbr.inf_sup(&rtr);
 
         let tmin = fast_max(inf.max(), T::zero());
@@ -277,6 +285,22 @@ mod tests {
         assert!(ray.intersection_slice_for_aabb(&aabb.aabb()).is_none());
     }
 
+    /// Ensure no slice is returned when the ray is in the plane of an AABB
+    /// face, which is a special case due to NaN's in the computation.
+    #[test]
+    fn test_in_plane_ray_slice() {
+        let aabb = UnitBox::new(0, TPoint3::new(0.0, 0.0, 0.0)).aabb();
+        let ray = TRay3::new(TPoint3::new(0.0, 0.0, -0.5), TVector3::new(1.0, 0.0, 0.0));
+        assert!(!ray.intersects_aabb(&aabb));
+        assert!(ray.intersection_slice_for_aabb(&aabb).is_none());
+
+        // Test a different ray direction, to ensure that order of `fast_max` (relevant
+        // to the result when NaN's are involved) doesn't matter.
+        let ray = TRay3::new(TPoint3::new(0.0, 0.5, 0.0), TVector3::new(0.0, 0.0, 1.0));
+        assert!(!ray.intersects_aabb(&aabb));
+        assert!(ray.intersection_slice_for_aabb(&aabb).is_none());
+    }
+
     proptest! {
         // Test whether a `Ray` which points at the center of an `Aabb` intersects it.
         #[test]

src/utils.rs

@@ -5,6 +5,7 @@ use crate::bvh::ShapeIndex;
 use crate::{aabb::Aabb, bvh::Shapes};
 
 use nalgebra::Scalar;
+use num::Float;
 
 /// Fast floating point minimum.  This function matches the semantics of
 ///
@@ -114,3 +115,9 @@ pub(crate) fn joint_aabb_of_shapes<T: BHValue, const D: usize, Shape: BHShape<T,
     }
     (aabb, centroid)
 }
+
+/// Returns `true` if and only if any of the floats returned by `iter` are NaN.
+#[inline(always)]
+pub(crate) fn has_nan<'a, T: Float + 'a>(iter: impl IntoIterator<Item = &'a T>) -> bool {
+    iter.into_iter().any(|f| f.is_nan())
+}