Failed to locate position. This is a bug in spade.

[mockersf]

I am randomly getting this panic when adding points to a Constrained Delaunay Triangulation.

I am running a physics simulation and extracting points of contact, which gives me a lot of points. I was able to get a small reproducer with only 4 edges:

use spade::{ConstrainedDelaunayTriangulation, Point2, Triangulation};

fn main() {
    let edges = [
        (
            Point2::new(17.064112, -17.96008),
            Point2::new(16.249594, -17.145563),
        ),
        (
            Point2::new(-25.290726, -24.435482),
            Point2::new(-5.6608872, -24.435482),
        ),
        (
            Point2::new(17.878626, -18.774595),
            Point2::new(15.435078, -16.331045),
        ),
        (
            Point2::new(16.331045, -21.951208),
            Point2::new(-3.29879, -24.435482),
        ),
    ];
    let mut cdt: ConstrainedDelaunayTriangulation<Point2<f32>> =
        ConstrainedDelaunayTriangulation::new();
    for edge in edges {
        println!("insert {:?}", edge.0);
        let point_a = cdt.insert(edge.0).unwrap();
        println!("insert {:?}", edge.1);
        let point_b = cdt.insert(edge.1).unwrap();
        println!("add constraint {:?}", edge);
        cdt.add_constraint_and_split(point_a, point_b, |v| {
            println!("inserted from split: {:?}", v);
            v
        });
    }
}

This logs:

insert Point2 { x: 17.064112, y: -17.96008 }
insert Point2 { x: 16.249594, y: -17.145563 }
add constraint (Point2 { x: 17.064112, y: -17.96008 }, Point2 { x: 16.249594, y: -17.145563 })
insert Point2 { x: -25.290726, y: -24.435482 }
insert Point2 { x: -5.6608872, y: -24.435482 }
add constraint (Point2 { x: -25.290726, y: -24.435482 }, Point2 { x: -5.6608872, y: -24.435482 })
insert Point2 { x: 17.878626, y: -18.774595 }
insert Point2 { x: 15.435078, y: -16.331045 }
add constraint (Point2 { x: 17.878626, y: -18.774595 }, Point2 { x: 15.435078, y: -16.331045 })
inserted from split: Point2 { x: 16.23077, y: -17.134615 }
insert Point2 { x: 16.331045, y: -21.951208 }
thread 'main' panicked at .cargo/registry/src/index.crates.io-6f17d22bba15001f/spade-2.10.0/src/delaunay_core/triangulation_ext.rs:505:17:
Failed to locate position. This is a bug in spade.
stack backtrace:
   0: rust_begin_unwind
             at /rustc/051478957371ee0084a7c0913941d2a8c4757bb9/library/std/src/panicking.rs:652:5
   1: core::panicking::panic_fmt
             at /rustc/051478957371ee0084a7c0913941d2a8c4757bb9/library/core/src/panicking.rs:72:14
   2: spade::delaunay_core::triangulation_ext::TriangulationExt::locate_with_hint_fixed_core
             at .cargo/registry/src/index.crates.io-6f17d22bba15001f/spade-2.10.0/src/delaunay_core/triangulation_ext.rs:505:17
   3: spade::delaunay_core::triangulation_ext::TriangulationExt::locate_with_hint_option_core
             at .cargo/registry/src/index.crates.io-6f17d22bba15001f/spade-2.10.0/src/delaunay_core/triangulation_ext.rs:221:9
   4: spade::delaunay_core::triangulation_ext::TriangulationExt::insert_with_hint_option_impl
             at .cargo/registry/src/index.crates.io-6f17d22bba15001f/spade-2.10.0/src/delaunay_core/triangulation_ext.rs:91:53
   5: spade::delaunay_core::triangulation_ext::TriangulationExt::insert_with_hint_option
             at .cargo/registry/src/index.crates.io-6f17d22bba15001f/spade-2.10.0/src/delaunay_core/triangulation_ext.rs:67:22
   6: spade::triangulation::Triangulation::insert
             at .cargo/registry/src/index.crates.io-6f17d22bba15001f/spade-2.10.0/src/triangulation.rs:536:9

If I manually insert the point from split (Point2::new(16.23077, -17.134615)) as start, it no longer panics.

[Stoeoef]

Interesting - thanks for reporting this!

It appears that this is an edge case related that can occur if some of the edges are nearly collinear:

Here is what happens:

• Spade inserts vertices 0, 1, 2, 3, 4, 5. Everything is fine.
• Then, it attempts to add a constraint edge going from 4 -> 5 . That edge just barely intersects edge 0 -> 1, generating vertex 6.
• For some reason, V6 creates an invalid topology - some of its surrounding faces are not ordered CCW .
• This invariant violation then leads to the crash when inserting an additional point.

Now, my memory tells me that there is some code that is supposed to handle this exact situation. I'll investigate why that doesn't fix the topology in this case...

As a workaround in the meantime you could consider to embed your input vertices in a big, outer triangle / rectangle that contains all other vertices (which can then later be removed) - I would assume the crash should only happen if the offending edge crosses the convex hull. edit: that likely won't help - my bad!

[Stoeoef]

I just realized that you are using an f32 based triangulation.

Spade calculates the intersection point using f32 vectors which makes a huge difference in this situation. I believe this issue should be fixed by using f64 for the edge intersection calculation.

[mockersf]

Yup, switching to f64 fixes it, and I didn't have other random panic in my app. My app uses f32 otherwise and now just move to f64 for the cdt, but I somewhat plan to support f64 all the way, would there be a risk of having the same kind of issues with numbers coming from the physics simulation in f64?

I was hoping to be able to keep using the same kind of float (f32 for now), I tried to look at spade code but didn't manage to find a fix.

[Stoeoef]

I guess the physics simulation "just" returns inaccurate results which might be okay, depending on the context. Delaunay triangulations unfortunately are reliant on high precision calculations, otherwise their internal invariants don't hold.

Anyway. Usually this is solved by using higher precision calculations internally. I'm still figuring out why this doesn't work in this case. My goal is that spade can be used with f32 in such situations, I'll keep this issue updated...

[mockersf]

thanks for looking into it!

at least I have a workaround for now 🙂

[Stoeoef]

Should be fixed with v2.11 (just released).

The fix should address a few additional edge cases as well that can arise due to imprecisions. Please keep in mind that the splitting result may not always match what you would expect - e.g. you wouldn't expect that an edge going from 4 to 5 (in the picture above) would create a split at all but rather be part of a straight constraint edge going from 4 -> 1 -> 0 -> 5. Yet, from a purely geometric point of view, that's unfortunately not a straight line...

With the fix, the example above contains the somewhat unintuitive (yet correct) constraint edges 0->6, 1->6, 2->3, 4->6, 5->6

Thanks again for the report! Please let me know if you find any other crashes, I'm interested to make this method as robust as it can reasonably get.

[hazelsparrow]

i'm getting the same error on an f64 triangulation consistently. does not happen on 2.6.0.

i'm not exactly sure what information you need but i will say that my triangulation is initialized from a voronoi diagram with 500k cells. i switched to using f64 a long time ago because with f32 i was running into precision issues--spade was complaining about invalid geometry or something like that, i don't remember exactly.

*edit: actually it might have been that it just hung indefinitely on trying to insert a new vertex. at any rate, it decidedly did not work with f32 but worked fine with f64s

i've downgraded to 2.6.0 for now but i'm curious to see if there is a better solution. thanks! :)

[mockersf]

trying to find a small reproducer could help! what I do is print points / edges in a way I can copy paste them in a simplified script like the one above, and then run it in a loop, but instead of using all edges I use choose_multiple from the rand crate to find a smaller case

[Stoeoef]

i'm getting the same error on an f64 triangulation consistently. does not happen on 2.6.0.

Thank you so much for the report! This bug is beginning to haunt me way more than anticipated :smile:

i'm not exactly sure what information you need

Most importantly, is this happening for a CDT or for a regular Delaunay triangulation? f32 / f64 shouldn't really affect this error (thus, it's equally likely with one or the other). It also shouldn't hang for f32 - that's strange.

As @mockersf mentioned, any data that lets us reproduce this is going to be helpful. Usually, only a small (<100) subset of vertices is responsible for issues like that. Finding those can be tricky though - randomly removing half of the input points and observing if the exception repeats might be the simplest way to to reduce this to a more manageable size.

That is, if you can share your input data, either as original or in a reduced form, I'd be very thankful - though I realize that this may not be an option for various reasons.

i've downgraded to 2.6.0 for now but i'm curious to see if there is a better solution. thanks! :)

The algorithm in 2.6 was different and works really well for regular Delaunay triangulations. I've changed that for CDT's due to other issues. You should be fine with v2.6 (though I'm interested in fixing this!)

[hazelsparrow]

let me see if I could get you a small reproducible example. but just for general information, i'm using your crate as follows. i start with an array of 500k points on a sphere (fibonacci points). i then use the stereographic projection to convert the points in the 3d space to a plane. then i use DelaunayTriangulation::bulk_load to load all the points into a delanay triangulation. it's at this point that i get the error. and it's at this point as well that it used to hang indefinitely with f32s, but i'm pretty sure it's because of precision issues -- the points in the stereographic projection were too cramped around the origin (0, 0), there wasn't enough precision to differentiate between them. i might be wrong about that one though -- it's been a while since i solved it.

let me try to get you a sample that reproduces the problem. it might not be possible -- maybe the issue only arises when the number of points exceeds a certain threshold though. but let me take a look

[hazelsparrow]

ok so i've done some digging and the problem i'm encountering is, in fact, with the triangulation.locate method, not the bulk_load at all. sorry for the confusion!

i wanted to highlight 2 things here which i'm not sure even qualify as bugs.

1) the behaviour of the locate method has changed when passing in points that lie far outside of the triangulation 2) the method blows up in certain scenarios in 2.11 but not in 2.6 (which is actually a good thing because it helped me uncover a bug in my code!)

here's the setup:

use std::f64::consts::PI;

use rand::{thread_rng, Rng};
use spade::{DelaunayTriangulation, Point2, Triangulation};

fn wrap(point: Point2<f64>) -> Point2<f64> {
  let quadrant = ((point.x.abs() / 90.) % 4.) as u8;
  let pole =  if point.x > 0. { 90. } else { -90. };
  let offset = point.x % 90.;
  let mut result = point.clone();

  match quadrant {
    0 => {
      result.x = offset;
    },
    1 => {
      result.x = pole - offset;
      result.y += 180.;
    },
    2 => {
      result.x = -offset;
      result.y += 180.;
    },
    3 => {
      result.x = offset - pole;
    },
    _ => {
      panic!("Invalid quadrant");
    }
  }

  if result.y > 180. || result.y < -180. {
    result.y -= ((result.y + 180.) / 360.).floor() * 360.;
  }

  result
}

fn generate_fibonacci_sphere(n: usize) -> Vec<Point2<f64>> {
  let mut rng = thread_rng();
  let jitter = 0.4;
  let nf = n as f64;
  let s = 3.6 / nf.sqrt();
  let dlong = PI * (3. - 5f64.sqrt());
  let dz = 2. / nf;

  let mut points = Vec::new();
  let mut z = 1. - dz / 2.;
  let mut long: f64 = 0.;

  for _ in 0..n {
    let r = (1. - z * z).sqrt();
    let mut lat_deg = z.asin().to_degrees();
    let mut long_deg = long.to_degrees();
    lat_deg += jitter * (rng.gen::<f64>() - rng.gen::<f64>()) * (lat_deg - (z - dz * 2. * PI * r / s).max(-1.).asin().to_degrees());
    long_deg += jitter * (rng.gen::<f64>() - rng.gen::<f64>()) * (s / r).to_degrees();

    points.push(wrap(Point2::new(lat_deg, long_deg)));
    long += dlong;
    z -= dz;
    if z < -1. { break; }
  }

  points
}

fn delanay(points: Vec<Point2<f64>>) -> DelaunayTriangulation<Point2<f64>> {
  DelaunayTriangulation::bulk_load(points).unwrap()
}

fn main() {
  let points = generate_fibonacci_sphere(1000);
  let triangulation = delanay(points);

  let distances = [5., 100_000_000., f64::NAN];

  for d in distances {
    let p = Point2::new(d, d);
    println!("locating {p:?}");
    let v = triangulation.locate(Point2::new(d, d));
    println!("found: {v:?}");
  }
}

passing in [5., 100_000_000., f64::NAN]

2.6

in 2.6, this prints out the following

locating Point2 { x: 5.0, y: 5.0 }
found: OnFace(FixedHandle { index: 14 })
locating Point2 { x: 100000000.0, y: 100000000.0 }
found: OutsideOfConvexHull(FixedHandle { index: 5825 })
locating Point2 { x: NaN, y: NaN }
found: OnFace(FixedHandle { index: 1884 }) <- this is the behaviour i was erroneously relying on previously

2.11

in 2.11, this prints out the following

locating Point2 { x: 5.0, y: 5.0 }
found: OnFace(FixedHandle { index: 12 })
locating Point2 { x: 100000000.0, y: 100000000.0 }
found: OutsideOfConvexHull(FixedHandle { index: 5909 })
locating Point2 { x: NaN, y: NaN }
found: OutsideOfConvexHull(FixedHandle { index: 5957 })

passing NaNs only

perhaps even more interesting is the fact that it seems to matter if NaNs are passed in with the very first locate call.

2.6

in 2.6, this:

fn main() {
  let points = generate_fibonacci_sphere(1000);
  let triangulation = delanay(points);

  let distances = [f64::NAN]; // note that the very first call will pass in NaNs

  for d in distances {
    let p = Point2::new(d, d);
    println!("locating {p:?}");
    let v = triangulation.locate(Point2::new(d, d));
    println!("found: {v:?}");
  }
}

produces the following output:

locating Point2 { x: NaN, y: NaN }
found: OnFace(FixedHandle { index: 9 })

2.11

in 2.11, this blows up:

locating Point2 { x: NaN, y: NaN }
thread 'main' panicked at C:\Users\grksparrow\.cargo\registry\src\index.crates.io-6f17d22bba15001f\spade-2.11.0\src\delaunay_core\triangulation_ext.rs:513:17:
Failed to locate position. This is a bug in spade.

[Stoeoef]

Oh, that's it - NaN. Indeed, locate does not do any check against NaN. Insertion, however, does. The behavior will be fairly wild and, indeed, depend on the previous locate call (that is expected).

The correct behavior is probably to panic - NaN are almost certainly a user error.

Thank you for providing all the information!

[Stoeoef]

"Fixed" with v2.12.0 (just released) - any locate call with NaN values will now always panic instead.