Document accuracy.

examples/co-spherical.rs

@@ -2,33 +2,102 @@ use miniball::{
  nalgebra::{Point3, Vector3},
  {Ball, Enclosing},
 };
-use rand::thread_rng;
+use rand::{seq::SliceRandom, thread_rng};
 use rand_distr::{Distribution, UnitSphere};
 use std::collections::VecDeque;
 
+type T = f64;
+
 fn main() {
+ // The epsilon chosen in [`Enclosing::contains()`] for [`Ball`].
+ //
+ // It should be large enough to prevent duplicate bounds.
+ let epsilon = T::EPSILON.sqrt();
+
+ let mut rng = thread_rng();
+
  let m = 100_000;
- println!("n = 3, m = {m}");
- println!();
+ let s = 8;
  let center = Vector3::new(-3.0, 7.0, 4.8);
- let radius = 3.0;
- let radius_squared = radius * radius;
- let mut points = UnitSphere
- .sample_iter(&mut thread_rng())
+ let inner_radius = 3.0 - epsilon;
+ let inner = UnitSphere
+ .sample_iter(&mut rng)
  .take(m)
- .map(Point3::<f64>::from)
- .map(|point| point * radius + center)
- .collect::<VecDeque<_>>();
- let ball = (0..8)
+ .map(Point3::from)
+ .map(|point| point * inner_radius + center)
+ .collect::<Vec<_>>();
+ let outer_radius = 3.0;
+ let outer = UnitSphere
+ .sample_iter(&mut rng)
+ .take(m)
+ .map(Point3::from)
+ .map(|point| point * outer_radius + center)
+ .collect::<Vec<_>>();
+ let mut inner_then_outer = VecDeque::new();
+ inner_then_outer.extend(inner.iter().cloned());
+ inner_then_outer.extend(outer.iter().cloned());
+ let mut outer_then_inner = VecDeque::new();
+ outer_then_inner.extend(outer.iter().cloned());
+ outer_then_inner.extend(inner.iter().cloned());
+ let mut random = Vec::new();
+ random.extend(inner.iter().cloned());
+ random.extend(outer.iter().cloned());
+ random.shuffle(&mut rng);
+ let mut random = VecDeque::from(random);
+
+ println!("Demonstrates the accuracy depending on the order of points.");
+ println!("Takes minimum of {s} samples permuted by move-to-front heuristic.");
+ println!("Rerun multiple times to see the accuracy changing for random points.");
+ println!();
+ println!("On-surface tolerance: 1{epsilon:+.1e}");
+ println!();
+ println!("n = 3, m = {m}, inner-then-outer (worst accuracy)");
+ println!();
+ let radius_squared = outer_radius * outer_radius;
+ let ball = (0..s)
+ .map(|_| {
+ let ball = Ball::enclosing_points(&mut inner_then_outer);
+ let epsilon = ball.radius_squared / radius_squared - 1.0;
+ println!("Sample with accuracy: 1{epsilon:+.1e}");
+ ball
+ })
+ .min()
+ .unwrap();
+ println!();
+ let epsilon = ball.radius_squared / radius_squared - 1.0;
+ println!("Result with accuracy: 1{epsilon:+.1e}");
+
+ println!();
+ println!("n = 3, m = {m}, outer-then-inner (best accuracy)");
+ println!();
+ let radius_squared = outer_radius * outer_radius;
+ let ball = (0..s)
+ .map(|_| {
+ let ball = Ball::enclosing_points(&mut outer_then_inner);
+ let epsilon = ball.radius_squared / radius_squared - 1.0;
+ println!("Sample with accuracy: 1{epsilon:+.1e}");
+ ball
+ })
+ .min()
+ .unwrap();
+ println!();
+ let epsilon = ball.radius_squared / radius_squared - 1.0;
+ println!("Result with accuracy: 1{epsilon:+.1e}");
+
+ println!();
+ println!("n = 3, m = {m}, random-inner/outer (from worst to best accuracy)");
+ println!();
+ let radius_squared = outer_radius * outer_radius;
+ let ball = (0..s)
  .map(|_| {
- let ball = Ball::enclosing_points(&mut points);
+ let ball = Ball::enclosing_points(&mut random);
  let epsilon = ball.radius_squared / radius_squared - 1.0;
- println!("sample with accuracy: 1{epsilon:+.1e}");
+ println!("Sample with accuracy: 1{epsilon:+.1e}");
  ball
  })
  .min()
  .unwrap();
  println!();
  let epsilon = ball.radius_squared / radius_squared - 1.0;
- println!("result with accuracy: 1{epsilon:+.1e}");
+ println!("Result with accuracy: 1{epsilon:+.1e}");
 }

src/enclosing.rs

@@ -94,7 +94,9 @@ where
  /// # Stability
  ///
  /// Due to floating-point inaccuracies, the returned ball might not exactly be the minimum for
- /// degenerate (e.g., co-spherical) `points`.
+ /// degenerate (e.g., co-spherical) `points`. The accuracy is depending on the shape and order
+ /// of `points` with an expected worst-case factor of `T::one() ± T::default_epsilon().sqrt()`
+ /// where `T::one()` is exact.
  ///
  /// # Example
  ///