"Empty point set" for non-empty point set (possible numerical stability issue)

[cormacrelf]

The following points cause a panic with "Empty point set", which is clearly a bit wrong because I'm counting six points there. I have a suspicion that if the circumscribing sphere part of the algorithm used https://crates.io/crates/robust, that this might not happen.

use nalgebra::{Point3, Vector3};
use std::collections::VecDeque;
type Ball = miniball::Ball<f32, nalgebra::Const<3>>;

#[test]
fn test_miniball_failure_case() {
    let mut buffer = [
        [0.8155192, -0.5787301, 0.0],
        [0.7634287, -0.6458921, 0.0],
        [0.7634287, -0.6458921, 8.751844e-16],
        [0.8993817, -0.43716416, -1.7844731e-15],
        [0.8155192, -0.5787301, -2.9261717e-15],
        [0.8993817, -0.43716416, 0.0],
    ]
    .into_iter()
    .map(|v| Vector3::<f32>::from_column_slice(&v[..]))
    .map(From::from)
    .collect::<VecDeque<Point3<f32>>>();
    let _ = Ball::enclosing_points(&mut buffer);
}

[cormacrelf]

The test passes if you delete the 3rd point, which is very close to being a duplicate of the second point.

[cormacrelf]

(A tricky thing about the robust crate is that there are only implementations for circles and 3d spheres, so using it would require specialisation of some kind for those cases somehow. So it is perhaps not a great solution. I'll have a proper crack at debugging why this is happening tomorrow.)

[n3vu0r]

Oh, good catch. Thanks! The panic message is wrong in case of numerical instability. I assume it happens when a point very close to the surface of the current ball is pushed onto the bounds vector resulting in an inaccurate decision of Self::with_bounds returning either None or Some(Ball) where Ball can have an unstable large radius.

The robust crate is interesting. It could help but I'm not sure if it's the best solution long term. There exists a variation of the Welzl algorithm which does use pivoting to find larger but stable balls with less recursion levels significantly improving on stability and performance. It requires more modifications for which I have currently no time but PRs are welcome. I can look up the paper if you are interested.

For now I've made a commit which simply retries the algorithm with the reordered set of points due to the move-to-front heuristic which should increase the chance of success and therefore hopefully succeed with a small number of retry attempts. Your test passes now with the second retry attempt and I guess it should work for a lot more problematic cases:

use miniball::Enclosing;
use nalgebra::{Point3, Vector3};
use std::collections::VecDeque;
type Ball = miniball::Ball<f32, nalgebra::Const<3>>;

fn main() {
    let mut buffer = [
        [0.8155192, -0.5787301, 0.0],
        [0.7634287, -0.6458921, 0.0],
        [0.7634287, -0.6458921, 8.751844e-16],
        [0.8993817, -0.43716416, -1.7844731e-15],
        [0.8155192, -0.5787301, -2.9261717e-15],
        [0.8993817, -0.43716416, 0.0],
    ]
    .into_iter()
    .map(|v| Vector3::<f32>::from_column_slice(&v[..]))
    .map(From::from)
    .collect::<VecDeque<Point3<f32>>>();
    let ball = Ball::enclosing_points(&mut buffer);
    println!("{ball:?}");
}

[cormacrelf]

Ah, you mean Bernd Gaertner's 1999 paper? The GPL code for it being here. If I were to port it, probably would be best to have it in a different crate so as not to force GPL onto your codebase. And then actually port the C++ that's been battle-tested and re-issue it under GPLv3, rather than attempt to divine the code from the paper, for example. I would think folks could purchase a commercial license for the ported Rust code from Gaertner/ETHZ under their existing scheme. Otherwise, if you only care about a specific dimension like 2 or 3, you can just compile some C++ and link it.

The use case for us is computing the bounding cone of some vectors. You compute the enclosing sphere of some unit vectors interpreted as points on the unit sphere, & then intersect it with the unit sphere to find a cone, somewhat similar to Shirmun & Abi-Ezzi's 1993 paper, although they went with a bounding box technique to approximate the enclosing sphere. We get a lot of degenerate cases due to where we're getting these vectors and how they're all normalised onto the unit sphere. So we are essentially a big fuzzer for miniball's stability WRT duplicate points, co-spherical points, zero- or near-zero-sized spheres, etc. I generated another test case pretty quickly, which is below. I don't expect you to play whack-a-mole with us! That could be a lot of work!

For the time being I think I will also compute the sphere from the diagonal of the bounding box, and pick the smaller of the two. Otherwise will look into the Gaertner things.

use nalgebra::{Point3, Vector3};
use std::collections::VecDeque;
type Ball = miniball::Ball<f32, nalgebra::Const<3>>;

#[test]
fn test_miniball_failure_case_duplicates_2() {
    let mut buffer = [
        [0.7822104, -0.6230143, 2.0600416e-16],
        [0.7473425, -0.66443896, 0.0],
        [0.7649825, -0.644051, 0.0],
        // This point duplicates the second point, above.
        // With this point included, you get a ball centred near the origin with radius ~ 1.0.
        // If you delete this point, the enclosing ball is correct
        // (a very small ball near the centroid of all these very similar points.)
        [0.7473425, -0.66443896, 0.0],
        [0.7822104, -0.6230143, 0.0],
        [0.7649825, -0.644051, 1.5816076e-15],
    ]
    .into_iter()
    .map(|v| Vector3::<f32>::from_column_slice(&v[..]))
    .map(From::from)
    .collect::<VecDeque<Point3<f32>>>();
    let ball = Ball::enclosing_points(&mut buffer);
    assert!(
        ball.radius_squared.sqrt() < 0.03,
        "expected radius < 0.03, got radius {:?}",
        ball.radius_squared.sqrt()
    );
}

[n3vu0r]

Yes, it's this paper. It describes an exact algorithm for non-exact numbers (floating-point) by motivating a threshold value (epsilon) to prevent inaccuracies from adding up but they mention bad things can still happen. You have an interesting use case, I have no idea if any exact algorithm will work well for your degenerate inputs. Alternatively, there are approximate algorithms which don't have this robustness problem at all but therefore their complexity also depends on epsilon, e.g., O(1/epsilon). This preprint from 2020 seems to be quite simple to implement.

I guess for this crate it's good to have both an exact and an approximate algorithm.

[n3vu0r]

As a last attempt on modifying the original algorithm, I've made a commit which encloses approximately co-spherical points. The epsilon is chosen roughly, such that it should also catch duplicates.