Merge pull request #44 from Andlon/various-changes

Various changes

Author: Andlon

Cargo.toml

@@ -23,7 +23,7 @@ proptest-support = [ "proptest", "fenris-geometry/proptest-support", "nalgebra/p
 [dependencies]
 nalgebra = { version = "0.31", features = [ "serde-serialize" ] }
 nalgebra-sparse = { version = "0.7", features = ["compare"] }
-davenport = "0.1.0"
+davenport = "0.1.1"
 vtkio = "0.6"
 num = "0.4"
 numeric_literals = "0.2.0"

fenris-geometry/src/lib.rs

@@ -11,6 +11,7 @@ use std::fmt::Debug;
 mod polygon;
 mod polytope;
 mod primitives;
+use crate::util::index_set_nth_power_iter;
 pub use polygon::*;
 pub use polytope::*;
 pub use primitives::*;
@@ -86,25 +87,25 @@ where
 
 #[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
 #[serde(bound(
-    serialize = "OVector<T, D>: Serialize",
-    deserialize = "OVector<T, D>: Deserialize<'de>"
+    serialize = "OPoint<T, D>: Serialize",
+    deserialize = "OPoint<T, D>: Deserialize<'de>"
 ))]
 pub struct AxisAlignedBoundingBox<T, D>
 where
     T: Scalar,
     D: DimName,
     DefaultAllocator: Allocator<T, D>,
 {
-    min: OVector<T, D>,
-    max: OVector<T, D>,
+    min: OPoint<T, D>,
+    max: OPoint<T, D>,
 }
 
 impl<T, D> Copy for AxisAlignedBoundingBox<T, D>
 where
     T: Scalar,
     D: DimName,
     DefaultAllocator: Allocator<T, D>,
-    OVector<T, D>: Copy,
+    OPoint<T, D>: Copy,
 {
 }
 
@@ -117,18 +118,18 @@ where
     D: DimName,
     DefaultAllocator: Allocator<T, D>,
 {
-    pub fn new(min: OVector<T, D>, max: OVector<T, D>) -> Self {
+    pub fn new(min: OPoint<T, D>, max: OPoint<T, D>) -> Self {
         for i in 0..D::dim() {
             assert!(min[i] <= max[i]);
         }
         Self { min, max }
     }
 
-    pub fn min(&self) -> &OVector<T, D> {
+    pub fn min(&self) -> &OPoint<T, D> {
         &self.min
     }
 
-    pub fn max(&self) -> &OVector<T, D> {
+    pub fn max(&self) -> &OPoint<T, D> {
         &self.max
     }
 }
@@ -140,7 +141,7 @@ where
     DefaultAllocator: Allocator<T, D>,
 {
     fn from(point: OPoint<T, D>) -> Self {
-        AxisAlignedBoundingBox::new(point.coords.clone(), point.coords)
+        AxisAlignedBoundingBox::new(point.clone(), point)
     }
 }
 
@@ -152,13 +153,21 @@ where
 {
     /// Computes the minimal bounding box which encloses both `this` and `other`.
     pub fn enclose(&self, other: &AxisAlignedBoundingBox<T, D>) -> Self {
-        let min = self.min.iter().zip(&other.min).map(|(a, b)| T::min(*a, *b));
+        let min = self
+            .min
+            .iter()
+            .zip(&other.min.coords)
+            .map(|(a, b)| T::min(*a, *b));
         let min = OVector::<T, D>::from_iterator(min);
 
-        let max = self.max.iter().zip(&other.max).map(|(a, b)| T::max(*a, *b));
+        let max = self
+            .max
+            .iter()
+            .zip(&other.max.coords)
+            .map(|(a, b)| T::max(*a, *b));
         let max = OVector::<T, D>::from_iterator(max);
 
-        AxisAlignedBoundingBox::new(min, max)
+        AxisAlignedBoundingBox::new(min.into(), max.into())
     }
 
     pub fn from_points<'a>(points: impl IntoIterator<Item = &'a OPoint<T, D>>) -> Option<Self> {
@@ -179,37 +188,37 @@ where
     }
 
     pub fn center(&self) -> OPoint<T, D> {
-        OPoint::from((self.max() + self.min()) / T::from_f64(2.0).unwrap())
+        OPoint::from((&self.max().coords + &self.min().coords) / T::from_f64(2.0).unwrap())
     }
 
     /// Uniformly scales each axis by the given scale amount, with respect to the center of
     /// the box.
     ///
     /// ```rust
     /// # use fenris_geometry::AxisAlignedBoundingBox;
-    /// use nalgebra::vector;
+    /// use nalgebra::{point, vector};
     /// use matrixcompare::assert_matrix_eq;
     ///
-    /// let aabb = AxisAlignedBoundingBox::new(vector![0.0, 0.0], vector![1.0, 1.0]);
+    /// let aabb = AxisAlignedBoundingBox::new(point![0.0, 0.0], point![1.0, 1.0]);
     /// let scaled = aabb.uniformly_scale(0.5);
     ///
-    /// assert_matrix_eq!(scaled.min(), vector![0.25, 0.25], comp = float);
-    /// assert_matrix_eq!(scaled.max(), vector![0.75, 0.75], comp = float);
+    /// assert_matrix_eq!(scaled.min().coords, vector![0.25, 0.25], comp = float);
+    /// assert_matrix_eq!(scaled.max().coords, vector![0.75, 0.75], comp = float);
     /// ```
     #[replace_float_literals(T::from_f64(literal).unwrap())]
     pub fn uniformly_scale(&self, scale: T) -> Self {
         assert!(scale >= T::zero());
         let s = scale;
         let (a, b) = (&self.min, &self.max);
-        let ref c = self.center().coords;
+        let ref c = self.center();
         Self {
             min: c + (a - c) * s,
             max: c + (b - c) * s,
         }
     }
 
     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 {
@@ -227,18 +236,91 @@ where
     ///
     /// ```rust
     /// # use fenris_geometry::AxisAlignedBoundingBox;
-    /// # use nalgebra::vector;
-    /// let aabb = AxisAlignedBoundingBox::new(vector![0.0, 0.0], vector![1.0, 1.0]);
+    /// # use nalgebra::point;
+    /// let aabb = AxisAlignedBoundingBox::new(point![0.0, 0.0], point![1.0, 1.0]);
     /// let grown = aabb.grow_uniformly(1.0);
-    /// assert_eq!(grown.min(), &vector![-1.0, -1.0]);
-    /// assert_eq!(grown.max(), &vector![2.0, 2.0]);
+    /// assert_eq!(grown.min(), &point![-1.0, -1.0]);
+    /// assert_eq!(grown.max(), &point![2.0, 2.0]);
     /// ```
     ///
     pub fn grow_uniformly(&self, distance: T) -> Self {
         let min = self.min().map(|b_i| b_i - distance);
         let max = self.max().map(|b_i| b_i + distance);
         Self::new(min, max)
     }
+
+    /// Creates an iterator over the corners of the bounding box.
+    pub fn corners_iter<'a>(&'a self) -> impl 'a + Iterator<Item = OPoint<T, D>>
+    where
+        DefaultAllocator: Allocator<usize, D>,
+    {
+        // We can enumerate the corners by looking at {0, 1}^D, i.e. the D-th power of the
+        // set {0, 1}, and associating 0 and 1 with min and max coordinates for the i-th axis.
+        index_set_nth_power_iter::<D>(2).map(move |multi_idx| {
+            OVector::<T, D>::from_fn(|idx, _| match multi_idx[idx] {
+                0 => self.min[idx].clone(),
+                1 => self.max[idx].clone(),
+                _ => unreachable!(),
+            })
+            .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,38 @@
-// 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::allocator::Allocator;
 use nalgebra::proptest::vector;
-use nalgebra::{DimName, Point2, Point3, Unit, Vector2, Vector3, U2, U3};
+use nalgebra::{DefaultAllocator, DimName, Point2, Point3, RealField, Unit, Vector2, Vector3, U2, U3};
 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/src/sdf.rs

@@ -1,5 +1,5 @@
 use fenris_traits::Real;
-use nalgebra::{Point2, Scalar, Vector2, U2};
+use nalgebra::{OPoint, Point2, Scalar, Vector2, U2};
 
 use crate::{AxisAlignedBoundingBox2d, BoundedGeometry};
 use numeric_literals::replace_float_literals;
@@ -68,8 +68,8 @@ where
     fn bounding_box(&self) -> AxisAlignedBoundingBox2d<T> {
         let eps = self.radius * T::from_f64(0.01).unwrap();
         AxisAlignedBoundingBox2d::new(
-            self.center - Vector2::repeat(T::one()) * (self.radius + eps),
-            self.center + Vector2::repeat(T::one()) * (self.radius + eps),
+            OPoint::from(self.center - Vector2::repeat(T::one()) * (self.radius + eps)),
+            OPoint::from(self.center + Vector2::repeat(T::one()) * (self.radius + eps)),
         )
     }
 }

fenris-geometry/src/util.rs

@@ -1,5 +1,6 @@
 use fenris_traits::Real;
-use nalgebra::{UnitVector3, Vector3};
+use nalgebra::allocator::Allocator;
+use nalgebra::{DefaultAllocator, DimName, OVector, UnitVector3, Vector3};
 
 pub fn compute_orthonormal_vectors_3d<T: Real>(vector: &UnitVector3<T>) -> [UnitVector3<T>; 2] {
     // Ported from
@@ -82,3 +83,32 @@ macro_rules! assert_line_segments_approx_equal {
         $crate::assert_line_segments_approx_equal_base!(msg_handler, $segment1, $segment2, abstol = $tol);
     }};
 }
+
+/// Returns an iterator that iterates over the elements of $\{ 0, \dots, n - 1 \}^N$, where $n$ is
+/// the number of indices in the index set $\{ 0, \dots, n - 1 \}$.
+///
+/// TODO: Currently only implicitly tested through [`AxisAlignedBoundingBox::corners_iter`].
+pub(crate) fn index_set_nth_power_iter<D: DimName>(n: usize) -> impl Iterator<Item = OVector<usize, D>>
+where
+    DefaultAllocator: Allocator<usize, D>,
+{
+    let mut multi_idx = OVector::from_element(0);
+    let mut i = 0;
+    std::iter::from_fn(move || {
+        if i < n.pow(D::dim() as u32) {
+            let item = Some(multi_idx.clone());
+            for j in 0..D::dim() {
+                if multi_idx[j] + 1 < n {
+                    multi_idx[j] += 1;
+                    break;
+                } else {
+                    multi_idx[j] = 0;
+                }
+            }
+            i += 1;
+            item
+        } else {
+            None
+        }
+    })
+}