geometry: remove unwrap in BoundingBox::from_points

Author: hhirtz

src/algorithms/hilbert_curve.rs

@@ -22,6 +22,9 @@ fn hilbert_curve_partition(
     part_count: usize,
     order: usize,
 ) {
+    debug_assert!(!partition.is_empty()); // to compute the bounding box
+    debug_assert_eq!(partition.len(), points.len());
+    debug_assert_eq!(partition.len(), weights.len());
     debug_assert!(order < 64);
 
     let compute_hilbert_index = hilbert_index_computer(points, order);
@@ -97,8 +100,9 @@ fn segment_to_segment(min: f64, max: f64, order: usize) -> impl Fn(f64) -> u64 {
     }
 }
 
+/// Panics if `points` is empty.
 fn hilbert_index_computer(points: &[Point2D], order: usize) -> impl Fn(&Point2D) -> u64 {
-    let mbr = OrientedBoundingBox::from_points(points);
+    let mbr = OrientedBoundingBox::from_points(points).unwrap();
     let aabb = mbr.aabb();
     let x_mapping = segment_to_segment(aabb.p_min.x, aabb.p_max.x, order);
     let y_mapping = segment_to_segment(aabb.p_min.y, aabb.p_max.y, order);
@@ -326,6 +330,9 @@ where
                 actual: self.order,
             });
         }
+        if part_ids.is_empty() {
+            return Ok(());
+        }
         hilbert_curve_partition(
             part_ids,
             points.as_ref(),

src/algorithms/k_means.rs

@@ -286,12 +286,16 @@ fn balanced_k_means_iter<const D: usize>(
     }
 }
 
-// This is the main load balance routine. It handles:
-//   - reordering the clusters according to their distance to a bounding box of all the points
-//   - assigning each point to the closest cluster according to the effective distance
-//   - checking partitions imbalance
-//   - increasing of diminishing clusters influence based on their imbalance
-//   - relaxing upper and lower bounds
+/// This is the main load balance routine. It handles:
+///   - reordering the clusters according to their distance to a bounding box of all the points
+///   - assigning each point to the closest cluster according to the effective distance
+///   - checking partitions imbalance
+///   - increasing of diminishing clusters influence based on their imbalance
+///   - relaxing upper and lower bounds
+///
+/// # Panics
+///
+/// Panics if `points` is empty.
 fn assign_and_balance<const D: usize>(
     points: &[PointND<D>],
     weights: &[f64],
@@ -316,11 +320,11 @@ fn assign_and_balance<const D: usize>(
     } = clusters;
     // compute the distances from each cluster center to the minimal
     // bounding rectangle of the set of points
-    let mbr = OrientedBoundingBox::from_points(points);
+    let obb = OrientedBoundingBox::from_points(points).unwrap();
     let distances_to_mbr = centers
         .par_iter()
         .zip(influences.par_iter())
-        .map(|(center, influence)| mbr.distance_to_point(center) * influence)
+        .map(|(center, influence)| obb.distance_to_point(center) * influence)
         .collect::<Vec<_>>();
 
     let mut zipped = centers

src/algorithms/recursive_bisection.rs

@@ -892,8 +892,11 @@ where
     W::Item: RcbWeight,
     W::Iter: rayon::iter::IndexedParallelIterator,
 {
-    let mbr = OrientedBoundingBox::from_points(points);
-    let points = points.par_iter().map(|p| mbr.obb_to_aabb(p));
+    let obb = match OrientedBoundingBox::from_points(points) {
+        Some(v) => v,
+        None => return Ok(()),
+    };
+    let points = points.par_iter().map(|p| obb.obb_to_aabb(p));
     // When the rotation is done, we just apply RCB
     simple_rcb(partition, points, weights, n_iter, tolerance)
 }

src/algorithms/z_curve.rs

@@ -48,20 +48,25 @@ fn z_curve_partition<const D: usize>(
     DefaultAllocator: Allocator<f64, Const<D>, Const<D>, Buffer = ArrayStorage<f64, D, D>>
         + Allocator<f64, DimDiff<Const<D>, Const<1>>>,
 {
-    let max_order = (HASH_TYPE_MAX as f64).log(f64::from(2u32.pow(D as u32))) as u32;
+    debug_assert_eq!(partition.len(), points.len());
+
+    let max_order = (HASH_TYPE_MAX as f64).log(f64::from(1 << D)) as u32;
     assert!(
         order <= max_order,
         "Cannot use the z-curve partition algorithm with an order > {} because it would currently overflow hashes capacity",
         max_order,
     );
 
     // Bounding box used to construct Point hashes
-    let mbr = OrientedBoundingBox::from_points(points);
+    let obb = match OrientedBoundingBox::from_points(points) {
+        Some(v) => v,
+        None => return,
+    };
 
     let mut permutation: Vec<_> = (0..points.len()).into_par_iter().collect();
 
     // reorder points
-    z_curve_partition_recurse(points, order, &mbr, &mut permutation);
+    z_curve_partition_recurse(points, order, &obb, &mut permutation);
 
     let points_per_partition = points.len() / part_count;
     let remainder = points.len() % part_count;
@@ -133,64 +138,6 @@ fn z_curve_partition_recurse<const D: usize>(
     })
 }
 
-// reorders a slice of Point3D in increasing z-curve order
-#[allow(unused)]
-pub(crate) fn z_curve_reorder<const D: usize>(points: &[PointND<D>], order: u32) -> Vec<usize>
-where
-    Const<D>: DimSub<Const<1>>,
-    DefaultAllocator: Allocator<f64, Const<D>, Const<D>, Buffer = ArrayStorage<f64, D, D>>
-        + Allocator<f64, DimDiff<Const<D>, Const<1>>>,
-{
-    let max_order = (HASH_TYPE_MAX as f64).log(f64::from(2u32.pow(D as u32))) as u32;
-    assert!(
-        order <= max_order,
-        "Cannot use the z-curve partition algorithm with an order > {} because it would currently overflow hashes capacity",
-        max_order,
-    );
-
-    let mut permu: Vec<_> = (0..points.len()).into_par_iter().collect();
-    z_curve_reorder_permu(points, permu.as_mut_slice(), order);
-    permu
-}
-
-// reorders a slice of indices such that the associated array of Point3D is sorted
-// by increasing z-order
-#[allow(unused)]
-pub(crate) fn z_curve_reorder_permu<const D: usize>(
-    points: &[PointND<D>],
-    permu: &mut [usize],
-    order: u32,
-) where
-    Const<D>: DimSub<Const<1>>,
-    DefaultAllocator: Allocator<f64, Const<D>, Const<D>, Buffer = ArrayStorage<f64, D, D>>
-        + Allocator<f64, DimDiff<Const<D>, Const<1>>>,
-{
-    let mbr = OrientedBoundingBox::from_points(points);
-    let hashes = permu
-        .par_iter()
-        .map(|idx| compute_hash(&points[*idx], order, &mbr))
-        .collect::<Vec<_>>();
-
-    permu.par_sort_by_key(|idx| hashes[*idx]);
-}
-
-fn compute_hash<const D: usize>(
-    point: &PointND<D>,
-    order: u32,
-    mbr: &OrientedBoundingBox<D>,
-) -> HashType {
-    let current_hash = mbr
-        .region(point)
-        .expect("Cannot compute the z-hash of a point outside of the current Mbr.");
-
-    if order == 0 {
-        HashType::from(current_hash)
-    } else {
-        (2_u128.pow(D as u32)).pow(order) * HashType::from(current_hash)
-            + compute_hash(point, order - 1, &mbr.sub_mbr(current_hash))
-    }
-}
-
 /// # Z space-filling curve algorithm
 ///
 /// The Z-curve uses space hashing to partition points. The points in the same part of a partition

src/geometry.rs

@@ -28,22 +28,25 @@ impl<const D: usize> BoundingBox<D> {
     ///
     /// This is the smallest rectangle (rectangular cuboid in 3D) that both
     /// contains all given points and is aligned with the axises.
-    pub fn from_points<P>(points: P) -> Self
+    ///
+    /// Returns `None` iff the given iterator is empty.
+    pub fn from_points<P>(points: P) -> Option<Self>
     where
         P: IntoParallelIterator<Item = PointND<D>>,
         P::Iter: IndexedParallelIterator,
     {
         let points = points.into_par_iter();
-        assert!(1 < points.len());
-
-        let (min, max) = points
+        if points.len() == 0 {
+            return None;
+        }
+        let (p_min, p_max) = points
             .fold_with(
                 (
                     PointND::<D>::from_element(std::f64::MAX),
                     PointND::<D>::from_element(std::f64::MIN),
                 ),
                 |(mut mins, mut maxs), vals| {
-                    for ((min, max), val) in mins.iter_mut().zip(maxs.iter_mut()).zip(vals.iter()) {
+                    for ((min, max), val) in mins.iter_mut().zip(maxs.iter_mut()).zip(&vals) {
                         if *val < *min {
                             *min = *val;
                         }
@@ -59,23 +62,19 @@ impl<const D: usize> BoundingBox<D> {
                     PointND::<D>::from_iterator(
                         mins_left
                             .into_iter()
-                            .zip(mins_right.into_iter())
+                            .zip(&mins_right)
                             .map(|(left, right)| left.min(*right)),
                     ),
                     PointND::<D>::from_iterator(
                         maxs_left
                             .into_iter()
-                            .zip(maxs_right.into_iter())
+                            .zip(&maxs_right)
                             .map(|(left, right)| left.max(*right)),
                     ),
                 )
             })
-            .unwrap();
-
-        Self {
-            p_min: min,
-            p_max: max,
-        }
+            .unwrap(); // fold_with yields at least one element.
+        Some(Self { p_min, p_max })
     }
 
     fn center(&self) -> PointND<D> {
@@ -213,7 +212,7 @@ impl<const D: usize> OrientedBoundingBox<D> {
     /// The arbitrarily-oriented *minimum* bounding box.
     ///
     /// The smallest box that contains all given points.
-    pub fn from_points(points: &[PointND<D>]) -> Self
+    pub fn from_points(points: &[PointND<D>]) -> Option<Self>
     where
         Const<D>: DimSub<Const<1>>,
         DefaultAllocator: Allocator<f64, Const<D>, Const<D>, Buffer = ArrayStorage<f64, D, D>>
@@ -225,13 +224,13 @@ impl<const D: usize> OrientedBoundingBox<D> {
         let base_change = householder_reflection(&vec);
         let obb_to_aabb = base_change.try_inverse().unwrap();
         let mapped = points.par_iter().map(|p| obb_to_aabb * p);
-        let aabb = BoundingBox::from_points(mapped);
+        let aabb = BoundingBox::from_points(mapped)?;
 
-        Self {
+        Some(Self {
             aabb,
             aabb_to_obb: base_change,
             obb_to_aabb,
-        }
+        })
     }
 
     /// Constructs a new Mbr that is a sub-Mbr of the current one.
@@ -365,23 +364,16 @@ mod tests {
             Point2D::from([4., 5.]),
         ];
 
-        let aabb = BoundingBox::from_points(points);
+        let aabb = BoundingBox::from_points(points).unwrap();
 
         assert_ulps_eq!(aabb.p_min, Point2D::from([0., 0.]));
         assert_ulps_eq!(aabb.p_max, Point2D::from([5., 5.]));
     }
 
     #[test]
     #[should_panic]
-    fn test_aabb2d_invalid_input_1() {
-        let _aabb = BoundingBox::<2>::from_points([]);
-    }
-
-    #[test]
-    #[should_panic]
-    fn test_aabb2d_invalid_input_2() {
-        let points = [Point2D::from([5., -9.2])];
-        let _aabb = BoundingBox::from_points(points);
+    fn test_aabb2d_invalid_input() {
+        BoundingBox::<2>::from_points([]).unwrap();
     }
 
     #[test]
@@ -422,29 +414,29 @@ mod tests {
 
     #[test]
     fn test_obb_center() {
-        let points = vec![
+        let points = [
             Point2D::from([0., 1.]),
             Point2D::from([1., 0.]),
             Point2D::from([5., 6.]),
             Point2D::from([6., 5.]),
         ];
 
-        let obb = OrientedBoundingBox::from_points(&points);
+        let obb = OrientedBoundingBox::from_points(&points).unwrap();
 
         let center = obb.center();
         assert_ulps_eq!(center, Point2D::from([3., 3.]))
     }
 
     #[test]
     fn test_obb_contains() {
-        let points = vec![
+        let points = [
             Point2D::from([0., 1.]),
             Point2D::from([1., 0.]),
             Point2D::from([5., 6.]),
             Point2D::from([6., 5.]),
         ];
 
-        let obb = OrientedBoundingBox::from_points(&points);
+        let obb = OrientedBoundingBox::from_points(&points).unwrap();
 
         assert!(!obb.contains(&Point2D::from([0., 0.])));
         assert!(obb.contains(&obb.center()));
@@ -458,14 +450,14 @@ mod tests {
 
     #[test]
     fn test_obb_quadrant() {
-        let points = vec![
+        let points = [
             Point2D::from([0., 1.]),
             Point2D::from([1., 0.]),
             Point2D::from([5., 6.]),
             Point2D::from([6., 5.]),
         ];
 
-        let obb = OrientedBoundingBox::from_points(&points);
+        let obb = OrientedBoundingBox::from_points(&points).unwrap();
 
         let none = obb.region(&Point2D::from([0., 0.]));
         let q1 = obb.region(&Point2D::from([1.2, 1.2]));
@@ -524,7 +516,7 @@ mod tests {
             Point3D::from([4., 5., 3.]),
         ];
 
-        let aabb = BoundingBox::from_points(points);
+        let aabb = BoundingBox::from_points(points).unwrap();
 
         assert_ulps_eq!(aabb.p_min, Point3D::from([0., 0., -2.]));
         assert_ulps_eq!(aabb.p_max, Point3D::from([5., 5., 5.]));
@@ -551,7 +543,7 @@ mod tests {
 
     #[test]
     fn test_obb_octant() {
-        let points = vec![
+        let points = [
             Point3D::from([0., 1., 0.]),
             Point3D::from([1., 0., 0.]),
             Point3D::from([5., 6., 0.]),
@@ -562,8 +554,7 @@ mod tests {
             Point3D::from([6., 5., 1.]),
         ];
 
-        let obb = OrientedBoundingBox::from_points(&points);
-        eprintln!("obb = {:#?}", obb);
+        let obb = OrientedBoundingBox::from_points(&points).unwrap();
 
         let none = obb.region(&Point3D::from([0., 0., 0.]));
         let octants = vec![

tools/src/bin/medit2svg.rs

@@ -229,7 +229,10 @@ where
             mesh.node_count(),
         )
     };
-    let bb = coupe::BoundingBox::<2>::from_points(coordinates.par_iter().cloned());
+    let bb = match coupe::BoundingBox::<2>::from_points(coordinates.par_iter().cloned()) {
+        Some(v) => v,
+        None => return Ok(()),
+    };
     let xmin = bb.p_min[0];
     let xmax = bb.p_max[0];
     let ymin = bb.p_min[1];