• Use ordered-float crate that provides strong ordering on floats by ignoring IEEE.754 spec in terms of handling NaN (which is what we did anyway). This allows to use *_by_key functions for comparisons on iterables.
• Change public API to return Option, so that when triangulation doesn't exist, None would be returned instead of panicking (this is more idiomatic in Rust).

I was using this library for a project of mine, and I came across a bug caused by the algorithm generating a triangulation that badly fails to meet the Delaunay condition.

Here is code that demonstrates the bug:

use delaunator::{triangulate, Point};

fn main {
    let points_raw = [
        [-1.0849334460551594, -1.0849334460551594], // 0
        [-1.5265847189170323, -1.5265847189170323], // 1
        [-2.095815826893781, -1.274850699584578], // 2
        [-2.770865690566727, -0.9763199041819535], // 3
        [-3.6843898126113546, -0.5723274031100705], // 4
        [-4.403658177176129, -0.25424163117441645], // 5
        [-5.02567003534954, 0.020833892521026076], // 6
        [-5.701084620214019, 0.3195259804640507], // 7
        [-6.561463270870451, 0.7000156846325452], // 8
        [-7.31105511135829, 1.0315115642859167], // 9
        [-8.0, 1.336187228503853], // 10
        [-7.339371017948894, 1.1048861305058357], // 11
        [-6.616986689211032, 0.8519630935920505], // 12
        [-5.816767071935726, 0.5717881676590966], // 13
        [-5.121128245254447, 0.3282293338915143], // 14
        [-4.512948676142796, 0.11529195763725286], // 15
        [-3.850960067633096, -0.11648517623155441], // 16
        [-3.1594534122386113, -0.3585972436874558], // 17
        [-2.288934450277422, -0.663385554827794], // 18
        [-1.6751076145244035, -0.8783001664294665], // 19
    ];
    let mut points = Vec::with_capacity(20);
    for [x, y] in &points_raw {
        points.push(Point { x: *x, y: *y });
    }
    let tris = triangulate(&points).unwrap();
    println!("{:?}", tris.triangles);
}

The output is:

[5, 6, 15, 6, 14, 15, 14, 16, 15, 15, 16, 5, 6, 7, 14, 16, 4, 5, 5, 4, 6, 7, 13, 14, 16, 17, 4, 14, 17, 16, 9, 8, 7, 7, 8, 13, 17, 3, 4, 4, 3, 6, 8, 12, 13, 19, 18, 17, 17, 18, 3, 6, 9, 7, 8, 9, 12, 18, 2, 3, 9, 11, 12, 12, 11, 13, 13, 11, 14, 14, 19, 17, 18, 19, 2, 19, 1, 2, 2, 10, 3, 10, 0, 3]

Which contains the triangle [10, 0, 3]. By plotting these points we can see that this triangle definitely shouldn't be part of the Delaunay triangulation, as the circumcircle of these points contains most of the other points: It looks like there are other triangles in this triangulation that fail the condition. I have idea why this specific set of points causes a problem.

I'm trying to triangulate a detailed map of the world. On several "chunks" (think islands on the map) of my data, delaunator crashes. Here is one example:

extern crate delaunator;

use delaunator::{Point, triangulate};

fn main() {
    let points = vec!(
        Point { x: 11.484139442443848, y: 65.20500183105469 },
        Point { x: 11.484139442443848, y: 65.2066421508789 },
        Point { x: 11.484999656677246, y: 65.20708465576172 },
        Point { x: 11.491666793823242, y: 65.20708465576172 },
        Point { x: 11.492444038391113, y: 65.2066650390625 },
        Point { x: 11.492500305175781, y: 65.20580291748047 },
        Point { x: 11.491639137268066, y: 65.20541381835938 },
        Point { x: 11.489860534667969, y: 65.20541381835938 },
        Point { x: 11.488332748413086, y: 65.20622253417969 },
        Point { x: 11.487471580505371, y: 65.2058334350586 },
        Point { x: 11.487388610839844, y: 65.20500183105469 },
        Point { x: 11.486528396606445, y: 65.20455932617188 },
        Point { x: 11.484916687011719, y: 65.20458221435547 },
        Point { x: 11.484139442443848, y: 65.20500183105469 },
    );
    let result = triangulate(&points).expect("No triangulation exists.");
    println!("{:?}", result.triangles);
}

When I run this code, delaunator panics:

> rustc -V
rustc 1.30.1 (1433507eb 2018-11-07)
> RUST_BACKTRACE=1 cargo run
    Finished dev [unoptimized + debuginfo] target(s) in 0.00s                                                                                 
     Running `target/debug/tmprusttest`
thread 'main' panicked at 'index out of bounds: the len is 14 but the index is 18446744073709551615', libcore/slice/mod.rs:2046:10
stack backtrace:
   0: std::sys::unix::backtrace::tracing::imp::unwind_backtrace
             at libstd/sys/unix/backtrace/tracing/gcc_s.rs:49
   1: std::sys_common::backtrace::print
             at libstd/sys_common/backtrace.rs:71
             at libstd/sys_common/backtrace.rs:59
   2: std::panicking::default_hook::{{closure}}
             at libstd/panicking.rs:211
   3: std::panicking::default_hook
             at libstd/panicking.rs:227
   4: std::panicking::rust_panic_with_hook
             at libstd/panicking.rs:477
   5: std::panicking::continue_panic_fmt
             at libstd/panicking.rs:391
   6: rust_begin_unwind
             at libstd/panicking.rs:326
   7: core::panicking::panic_fmt
             at libcore/panicking.rs:77
   8: core::panicking::panic_bounds_check
             at libcore/panicking.rs:59
   9: <usize as core::slice::SliceIndex<[T]>>::index
             at libcore/slice/mod.rs:2046
  10: core::slice::<impl core::ops::index::Index<I> for [T]>::index
             at libcore/slice/mod.rs:1914
  11: <alloc::vec::Vec<T> as core::ops::index::Index<I>>::index
             at liballoc/vec.rs:1725
  12: delaunator::Triangulation::legalize
             at /home/trehn/.cargo/registry/src/github.com-1ecc6299db9ec823/delaunator-0.2.0/src/lib.rs:232
  13: delaunator::triangulate
             at /home/trehn/.cargo/registry/src/github.com-1ecc6299db9ec823/delaunator-0.2.0/src/lib.rs:484
  14: tmprusttest::main
             at src/main.rs:27
  15: std::rt::lang_start::{{closure}}
             at libstd/rt.rs:74
  16: std::panicking::try::do_call
             at libstd/rt.rs:59
             at libstd/panicking.rs:310
  17: __rust_maybe_catch_panic
             at libpanic_unwind/lib.rs:103
  18: std::rt::lang_start_internal
             at libstd/panicking.rs:289
             at libstd/panic.rs:392
             at libstd/rt.rs:58
  19: std::rt::lang_start
             at libstd/rt.rs:74
  20: main
  21: __libc_start_main
  22: _start

Any idea what's going on here?

This pull requests replaces all dependencies on libstd to dependencies on libcore and liballoc, so that this algorithm can be used on targets that do not have access to the standard library.

This PR addresses the robustness issue tracked in #15, based on mapbox/delaunator/pull/68

Changes

• Takes dependency on robust, a pure rust implementation of Shewchuk's robust predicates
• Uses robust::orient2d implementation for the orient() checks
• Split the test fixtures from a single test case into separate tests to make it easier to see errors when running cargo test
• Adds a example that generates a svg of the triangulation
• Brought in the additional fixture tests from the javascript implementation

Benchmark

The robust checks are slower by about ~8-10%

triangulate/100         time:   [18.330 us 18.350 us 18.392 us]
                        change: [+6.7385% +7.6582% +8.8075%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 1 outliers among 20 measurements (5.00%)
  1 (5.00%) high severe
triangulate/1000        time:   [275.33 us 276.51 us 277.46 us]
                        change: [+12.091% +12.617% +13.055%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 2 outliers among 20 measurements (10.00%)
  2 (10.00%) high mild
triangulate/10000       time:   [3.4929 ms 3.4974 ms 3.5028 ms]
                        change: [+8.8572% +9.3444% +9.7613%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 2 outliers among 20 measurements (10.00%)
  2 (10.00%) high severe
triangulate/100000      time:   [50.839 ms 50.982 ms 51.150 ms]
                        change: [+6.2803% +6.6616% +7.0349%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 1 outliers among 20 measurements (5.00%)
  1 (5.00%) high mild

Commentary

During testing, I was not able to reproduce issue mapbox/delaunator/issues/61 with or without this change.

For example, this is what the triangulation output looks on current master:

This is what it looks like with the new orient2d checks:

And this is what it looks as reported on the original issue:

While testing, I came across #18, so I added the checks for connectivity on the test case but left commented out for now, until we can address that issue (which is unrelated to this).

Additionally, this change may lead to supporting the ask in #4, at least for f32 if we move from our Point struct to the Coord struct defined in the robust crate.

Overall, the performance impact is significant but worth it in my mind to guarantee correctness. I am just unsure about these changes because I couldn't repro the original issue yet. Please let me know your thoughts.

This is the port for the fix to remove recursion in the legalize method. First issue tracked in #15

From benchmarks, it seems the loop approach is slightly worse than the recursion approach (numbers below).

My first commit was actually much worse, because I had the edge_stack in the stack, which was getting initialized on each method call. I had to move the edge_stack in the triangulation type to avoid a global variable, which would make cause problems for folks wanting to use this create with multiple threads. Still, it is not ideal to keep it there as it takes unnecessary memory. Another alternative is to use a unsafe array initializer but then it taints the create with unsafe.

I will try to dig a bit more later to see if I can tweak this; but my guess is that I won't be able to beat the compile and the recursion implementation was probably getting optimized out by the compiler quite well.

I am unsure if we should take this change, but I am sending the PR anyways with the details so we can review and decide.

Bench for master right now

triangulate/100         time:   [17.412 us 17.510 us 17.658 us]                            
Found 1 outliers among 20 measurements (5.00%)
  1 (5.00%) high mild
triangulate/1000        time:   [249.41 us 250.44 us 252.33 us]                            
Found 3 outliers among 20 measurements (15.00%)
  1 (5.00%) high mild
  2 (10.00%) high severe
triangulate/10000       time:   [3.2132 ms 3.2291 ms 3.2552 ms]                              
Found 4 outliers among 20 measurements (20.00%)
  4 (20.00%) high severe
triangulate/100000      time:   [47.804 ms 48.003 ms 48.293 ms]                              
Found 1 outliers among 20 measurements (5.00%)
  1 (5.00%) high severe

Bench after this change:

triangulate/100         time:   [19.405 us 19.498 us 19.591 us]                            
                        change: [+10.503% +11.068% +11.579%] (p = 0.00 < 0.05)
                        Performance has regressed.
Found 3 outliers among 20 measurements (15.00%)
  3 (15.00%) high severe
triangulate/1000        time:   [251.40 us 253.38 us 255.88 us]                            
                        change: [+0.3120% +0.9845% +1.6889%] (p = 0.01 < 0.05)
                        Change within noise threshold.
Found 1 outliers among 20 measurements (5.00%)
  1 (5.00%) high mild
triangulate/10000       time:   [3.2142 ms 3.2216 ms 3.2329 ms]                              
                        change: [-0.6539% -0.0132% +0.5873%] (p = 0.97 > 0.05)
                        No change in performance detected.
triangulate/100000      time:   [48.406 ms 48.505 ms 48.613 ms]                              
                        change: [+0.4143% +1.0451% +1.5393%] (p = 0.00 < 0.05)
                        Change within noise threshold.
Found 4 outliers among 20 measurements (20.00%)
  2 (10.00%) low mild
  1 (5.00%) high mild
  1 (5.00%) high severe

Thanks for the port!

The docs indicate that the indices returned by Triangulation.hull are counter-clockwise.

hull: Vec<usize>
[−]
A vector of indices that reference points on the convex hull of the triangulation, counter-clockwise.

However I see them being returned clockwise. Which one is the right expectation?

Sample test:

use delaunator::*;

fn main() {
    let p = vec![
        Point { x: 0.0, y: 0.0 },
        Point { x: 1.0, y: 1.0 },
        Point { x: -1.0, y: 1.0 },
    ];
    triangulate(&p)
        .unwrap().hull.iter()
        .for_each(|s| println!("Hull {}: {:?}", *s, p[*s]));
}

for which I get:

Hull 0: [0, 0]
Hull 2: [-1, 1]
Hull 1: [1, 1]

Points distributed in a square pattern now have an incorrect hull (left side is v1 delaunator-rs, and right side is the latest delaunator js):

Here's the result for v0.2.1 that shows the hull matched the js result

Likely a regression from the robust checks.

I'd like to use this algorithm in my breadx project, as a way of breaking shapes into triangles before sending them to the X11 server. As a test, I ran the algorithm on this shape:

    // create an arbitrary polygon and convert it to triangles
    let t: Vec<_> = tesselate_shape(&vec![
        Pointfix {
            x: double_to_fixed(100.0),
            y: double_to_fixed(100.0),
        },
        Pointfix {
            x: double_to_fixed(100.0),
            y: double_to_fixed(150.0),
        },
        Pointfix {
            x: double_to_fixed(150.0),
            y: double_to_fixed(250.0),
        },
        Pointfix {
            x: double_to_fixed(200.0),
            y: double_to_fixed(250.0),
        },
        Pointfix {
            x: double_to_fixed(200.0),
            y: double_to_fixed(150.0),
        },
        Pointfix {
            x: double_to_fixed(250.0),
            y: double_to_fixed(100.0),
        },
        Pointfix {
            x: double_to_fixed(300.0),
            y: double_to_fixed(150.0),
        },
        Pointfix {
            x: double_to_fixed(300.0),
            y: double_to_fixed(300.0),
        },
        Pointfix {
            x: double_to_fixed(150.0),
            y: double_to_fixed(450.0),
        },
        Pointfix {
            x: double_to_fixed(50.0),
            y: double_to_fixed(250.0),
        },
        Pointfix {
            x: double_to_fixed(50.0),
            y: double_to_fixed(150.0),
        },
    ])
    .collect();

I implemented tesselate_shape using the following:

/// From the given set of points, return an iterator over the triangles.
#[inline]
pub fn tesselate_shape<'a>(points: &'a [Pointfix]) -> impl Iterator<Item = Triangle> + 'a {
    let floating_points: Vec<Point> = points
        .iter()
        .copied()
        .map(|Pointfix { x, y }| Point {
            x: fixed_to_double(x),
            y: fixed_to_double(y),
        })
        .collect();

    let vector = match triangulate(&floating_points) {
        Some(t) => t.triangles ,
        None => vec![],
    };

    vector
        .into_iter()
        .map(move |index| &points[index])
        .copied()
        .scan(ArrayVec::<[Pointfix; 3]>::new(), |av, point| {
            av.push(point);
            if av.len() == 3 {
                let [p1, p2, p3] = mem::take(av).into_inner();
                Some(Some(Triangle { p1, p2, p3 }))
            } else {
                Some(None)
            }
        })
        .flatten()
}

Note: XRender does things in fixed-point 32-bit numbers by default, so I have to convert them to and from floating point notation.

When I run this, it produces radically different results than expected. Here's what I should see; this is what I get when I manually triangulated the shape:

However, this happens when I use delaunator:

I feel like this is probably a result of me misunderstanding how the crate is supposed to work; is there something I'm missing?

Depending on the input, a triangulation may be returned, however not all points are connected.

For example, for the input below,

[
    [0, 0],
    [0.0000000000000002220446049250313, 0.0000000000000002220446049250313],
    [0.0000000000000004440892098500626, 0.0000000000000002220446049250313],
    [-0.0000000000000004440892098500626, 0.0000000000000002220446049250313],
    [-0.0000000000000005551115123125783, 0.0000000000000004440892098500626]
]

The following triangulation is returned:

Triangulation {
    triangles: [
        0,
        1,
        2,
        0,
        3,
        1,
        1,
        3,
        2,
    ],
    halfedges: [
        5,
        8,
        18446744073709551615,
        18446744073709551615,
        6,
        0,
        4,
        18446744073709551615,
        1,
    ],
    hull: [
        0,
        3,
        2,
    ],
}

which does not include the last point.

During hull walk, the orientation check was missing taking into account collinear points ( roust == 0). This was my miss when I up-took the robustness check changes.

Fixes #24

• ~~Eliminate recursion in edge flipping to improve performance~~ https://github.com/mapbox/delaunator/pull/35
• [x] Return an ordered set of points if they are collinear https://github.com/mapbox/delaunator/pull/47 — done in #17
• [ ] Add an update method for iterative algorithms https://github.com/mapbox/delaunator/pull/48
• [x] Fix a race condition that caused an infinite loop https://github.com/mapbox/delaunator/pull/49 — done in #13
• [x] Validate that hull is convex in robustness tests https://github.com/mapbox/delaunator/pull/51
• [x] ❗ Fully robust orientation check based on Shewchuk predicates https://github.com/mapbox/delaunator/pull/68
• [x] Make sure it passes all the latest Delaunator fixtures/tests

cc @andreesteve

There is a panic on the following line when one of the values is NaN: https://github.com/mourner/delaunator-rs/blob/82ddbdb33651139414f2fe581f67dc12522dec32/src/lib.rs#L434

You can get around this by implementing a default ordering with something like

dists.sort_unstable_by(|&(_, da), &(_, db)| da.partial_cmp(&db).unwrap_or(std::cmp::Ordering::Less));

Or by checking if a given argument contains a NaN in the first place if NaN breaks further work. I think having one of these and returning None is preferable to a panic inside the method. At the very least this behavior should be documented.

Currently delaunator only works with f64. I need to learn how Rust generics and numeric traits work to enable using the library with other numeric types.