AABB::furthest_point_to

Author: Andlon

fenris-geometry/src/lib.rs

@@ -210,7 +210,7 @@ where
     }
 
     pub fn contains_point(&self, point: &OPoint<T, D>) -> bool {
-        (0..D::dim()).all(|dim| point[dim] > self.min[dim] && point[dim] < self.max[dim])
+        (0..D::dim()).all(|dim| point[dim] >= self.min[dim] && point[dim] <= self.max[dim])
     }
 
     pub fn intersects(&self, other: &Self) -> bool {
@@ -257,6 +257,63 @@ where
                 }).into()
             })
     }
+
+    /// Compute the point in the bounding box furthest away from the given point.
+    ///
+    /// # Panics
+    ///
+    /// Panics if two distances cannot be ordered. This typically only happens if
+    /// one of the numbers is not a number (NaN) or the comparison is not sensible, such as
+    /// comparing two infinities. Since given finite coordinates no distance should be infinite,
+    /// this method will realistically only panic in cases where one of the points
+    /// --- either of the bounding box or the query point --- has components that are not
+    /// finite numbers.
+    #[replace_float_literals(T::from_f64(literal).unwrap())]
+    pub fn furthest_point_to(&self, point: &OPoint<T, D>) -> OPoint<T, D>
+    where
+        DefaultAllocator: Allocator<usize, D>
+    {
+        // It turns out that we can choose, along each dimension, the point in the interval
+        // [a_i, b_i] furthest away from p_i.
+        point.coords.zip_zip_map(
+            &self.min.coords,
+            &self.max.coords,
+            |p_i, a_i, b_i| {
+                let mid = (a_i + b_i) / 2.0;
+                if p_i < mid {
+                    b_i
+                } else {
+                    a_i
+                }
+            }
+        ).into()
+    }
+
+    /// The squared distance to the point in the bounding box furthest away from the given point.
+    ///
+    /// # Panics
+    ///
+    /// Panic behavior is identical to [`furthest_point_to`](Self::furthest_point_to).
+    pub fn max_dist2_to(&self, point: &OPoint<T, D>) -> T
+    where
+        // TODO: Use DimAllocator and SmallDim
+        DefaultAllocator: Allocator<usize, D>
+    {
+        (self.furthest_point_to(point) - point).norm_squared()
+    }
+
+    /// The distance to the point in the bounding box furthest away from the given point.
+    ///
+    /// # Panics
+    ///
+    /// Panic behavior is identical to [`max_dist2_to`](Self::max_dist2_to).
+    pub fn max_dist_to(&self, point: &OPoint<T, D>) -> T
+    where
+    // TODO: Use DimAllocator and SmallDim
+        DefaultAllocator: Allocator<usize, D>
+    {
+        self.max_dist2_to(point).sqrt()
+    }
 }
 
 fn intervals_intersect<T: Real>([l1, u1]: [T; 2], [l2, u2]: [T; 2]) -> bool {

fenris-geometry/src/proptest.rs

@@ -1,10 +1,32 @@
-// use crate::procedural::create_rectangular_uniform_quad_mesh_2d;
 use crate::Orientation::Counterclockwise;
-use crate::{HalfPlane, LineSegment2d, Orientation, Quad2d, Triangle, Triangle2d, Triangle3d};
+use crate::{AxisAlignedBoundingBox, AxisAlignedBoundingBox2d, AxisAlignedBoundingBox3d, HalfPlane, LineSegment2d, Orientation, Quad2d, Triangle, Triangle2d, Triangle3d};
 use nalgebra::proptest::vector;
-use nalgebra::{DimName, Point2, Point3, Unit, Vector2, Vector3, U2, U3};
+use nalgebra::{DimName, Point2, Point3, Unit, Vector2, Vector3, U2, U3, DefaultAllocator, RealField};
+use nalgebra::allocator::Allocator;
 use proptest::prelude::*;
 
+pub fn aabb<D>() -> impl Strategy<Value=AxisAlignedBoundingBox<f64, D>>
+where
+    D: DimName,
+    DefaultAllocator: Allocator<f64, D>
+{
+    let value_strategy = -10.0 .. 10.0;
+    (vector(value_strategy.clone(), D::name()), vector(value_strategy, D::name()))
+        .prop_map(|(a, b)| {
+            let min = a.zip_map(&b, |a_i, b_i| RealField::min(a_i, b_i));
+            let max = a.zip_map(&b, |a_i, b_i| RealField::max(a_i, b_i));
+            AxisAlignedBoundingBox::new(min.into(), max.into())
+        })
+}
+
+pub fn aabb2() -> impl Strategy<Value=AxisAlignedBoundingBox2d<f64>> {
+    aabb()
+}
+
+pub fn aabb3() -> impl Strategy<Value=AxisAlignedBoundingBox3d<f64>> {
+    aabb()
+}
+
 pub fn vector2() -> impl Strategy<Value = Vector2<f64>> {
     vector(-10.0..10.0, U2::name())
 }

fenris-geometry/tests/unit_tests/aabb.rs

@@ -1,7 +1,11 @@
+use matrixcompare::assert_scalar_eq;
 use fenris_geometry::AxisAlignedBoundingBox;
 use nalgebra::{DefaultAllocator, DimName, OPoint, point, U2};
 use nalgebra::allocator::Allocator;
+use proptest::prelude::*;
 use fenris::allocators::DimAllocator;
+use fenris_geometry::proptest::{aabb2, aabb3, point2, point3};
+use nalgebra::distance_squared;
 
 #[test]
 fn aabb_intersects_2d() {
@@ -110,6 +114,85 @@ fn test_aabb_corners_iter() {
         ].map(OPoint::from);
         assert_unordered_eq!(corners(&aabb), expected);
     }
+}
+
+#[test]
+fn test_furthest_point_2d() {
+    let aabb = AxisAlignedBoundingBox::new(point![1.0, 1.0], point![2.0, 3.0]);
+
+    {
+        let p = point![0.0, 0.0];
+        let q = aabb.furthest_point_to(&p);
+        assert_eq!(q, point![2.0, 3.0]);
+        // We check all the convenience method results here just to have a unit test that does
+        // this, but the consistency is separately tested by proptests, which is why
+        // we don't do this for every test
+        let dist2: f64 = distance_squared(&q, &p);
+        assert_scalar_eq!(dist2, 13.0);
+        assert_scalar_eq!(dist2, aabb.max_dist2_to(&p));
+        assert_scalar_eq!(aabb.max_dist_to(&p), f64::sqrt(13.0));
+    }
+
+    {
+        let p = point![1.5, 2.0];
+        let q = aabb.furthest_point_to(&p);
+        // The exact point is not unique: any corner will be applicable
+        assert_scalar_eq!(distance_squared(&q, &p), 1.25);
+    }
+}
+
+proptest! {
 
+    #[test]
+    fn aabb_max_dists_agree_with_furthest_point_2d(point in point2(), aabb in aabb2()) {
+        let q = aabb.furthest_point_to(&point);
+        let dist2 = distance_squared(&q, &point);
+        prop_assert_eq!(aabb.max_dist2_to(&point), dist2);
+        prop_assert_eq!(aabb.max_dist_to(&point), dist2.sqrt());
+    }
+
+    #[test]
+    fn aabb_max_dists_agree_with_furthest_point_3d(point in point3(), aabb in aabb3()) {
+        let q = aabb.furthest_point_to(&point);
+        let dist2 = distance_squared(&q, &point);
+        prop_assert_eq!(aabb.max_dist2_to(&point), dist2);
+        prop_assert_eq!(aabb.max_dist_to(&point), dist2.sqrt());
+    }
+
+    #[test]
+    fn aabb_furthest_point_2d(p in point2(), aabb in aabb2()) {
+        // The furthest point in the AABB is *always* a corner, so we must satisfy
+        //  dist(p, q) <= dist(p, c)
+        // for all corners c and furthest point q
+        let q = aabb.furthest_point_to(&p);
+        let further_away_than_all_corners = aabb.corners_iter()
+            .all(|corner| distance_squared(&p, &q) >= distance_squared(&p, &corner));
+        prop_assert!(further_away_than_all_corners);
+
+        // The result should be exactly one of the corners, and since there are no floating
+        // point operations applied to the result (all numbers are just copied),
+        // there should also be no round-off error in the result, so we should
+        // be safe to check if the point is contained in the AABB, despite the fact that it
+        // resides exactly on the boundary!
+        prop_assert!(aabb.contains_point(&q));
+    }
+
+    #[test]
+    fn aabb_furthest_point_3d(p in point3(), aabb in aabb3()) {
+        // The furthest point in the AABB is *always* a corner, so we must satisfy
+        //  dist(p, q) <= dist(p, c)
+        // for all corners c and furthest point q
+        let q = aabb.furthest_point_to(&p);
+        let further_away_than_all_corners = aabb.corners_iter()
+            .all(|corner| distance_squared(&p, &q) >= distance_squared(&p, &corner));
+        prop_assert!(further_away_than_all_corners);
+
+        // The result should be exactly one of the corners, and since there are no floating
+        // point operations applied to the result (all numbers are just copied),
+        // there should also be no round-off error in the result, so we should
+        // be safe to check if the point is contained in the AABB, despite the fact that it
+        // resides exactly on the boundary!
+        prop_assert!(aabb.contains_point(&q));
+    }
 
 }