Spade 1.4.0 works correctly, but I'm seeing test failures on nearest-neighbour queries as of Spade 1.5.0: I'm getting a None value when querying. The worrying thing is that I can't reproduce them using Spade's native types.

This fails for our types: https://github.com/urschrei/rust-geo/blob/5d888080885f52dc656ebf7af3ae4d423bea5aed/geo/src/algorithm/euclidean_distance.rs#L945-L966

Source for edges which are bulk-loaded into the tree (edges built using windows(2)) Source for points which are used to retrieve the nearest-neighbour edge from the tree

The test fails when querying using the tree using the point (6.0001320971764285, -6.483226510722902)

The Spade trait impls are as follows:

For Point (same as the Point2 type):

https://github.com/georust/rust-geo/blob/b8354152381f735ab765847b656e12afc5117869/geo-types/src/point.rs#L297-L336

For Line (same as the SimpleEdge type):

https://github.com/georust/rust-geo/blob/b8354152381f735ab765847b656e12afc5117869/geo-types/src/line.rs#L145-L167]

Things I've checked:

• That the point values making up the Lines in the RTree are the same as those which make up the SimpleEdges when running the test using Spade types ✅
• That the mbr values are the same for our types and for Spade Types ✅
• That the distance2 types are the same for our types and for Spade types (in some cases, they vary in the least significant digit, so I don't think that's the problem)

Correct (using Spade Types):

95.29292557779627
100.53487738887094
105.67537257174098
110.66490534771486
115.4554237886039
120.00079258305857
124.25723734483688
128.1837661840736
131.74256448158783
134.8993590643375
137.6237482748112
64.54457621139609
69.33509465228511
74.324627428259
79.46512261112903
84.70707442220369
89.99999999999997
139.88949475560952
141.6747781295393
142.96240514177447
143.73997524029863
143.73997524029866
45.100640935662554
48.257435518412194
51.81623381592642
55.74276265516312
59.999207416941445
142.96240514177447
141.6747781295393
139.88949475560952
137.62374827481125
134.89935906433752
131.74256448158786
42.37625172518882
40.110505244390524
38.32522187046073
37.03759485822556
36.26002475970139
36.00000000000001
110.66490534771505
115.45542378860407
120.00079258305868
124.257237344837
128.18376618407368
36.00000000000001
36.26002475970138
37.03759485822556
38.32522187046073
40.11050524439051
79.46512261112939
84.70707442220402
90.00000000000028
95.29292557779655
100.53487738887117
105.67537257174118
42.376251725188816
45.10064093566255
48.257435518412194
51.81623381592641
55.742762655163105
59.99920741694144
59.999207416941815
64.54457621139646
69.33509465228548
74.32462742825936

Using our types:

95.29292557779625
100.53487738887094
105.67537257174098
110.66490534771485
115.45542378860391
120.00079258305858
124.2572373448369
128.1837661840736
131.74256448158786
134.8993590643375
137.6237482748112
64.5445762113961
69.3350946522851
74.324627428259
79.46512261112905
84.70707442220369
89.99999999999997
139.88949475560952
141.6747781295393
142.96240514177447
143.7399752402986
143.73997524029866
45.100640935662554
48.25743551841219
51.81623381592642
55.74276265516311
59.999207416941445
142.96240514177447
141.6747781295393
139.88949475560952
137.62374827481125
134.89935906433752
131.74256448158786
42.37625172518882
40.110505244390524
38.32522187046073
37.03759485822555
36.260024759701395
36.000000000000014
110.66490534771503
115.45542378860405
120.0007925830587
124.25723734483701
128.18376618407368
36.000000000000014
36.26002475970138
37.03759485822555
38.32522187046073
40.11050524439051
79.46512261112939
84.70707442220402
90.00000000000027
95.29292557779657
100.53487738887115
105.67537257174116
42.37625172518881
45.100640935662554
48.25743551841219
51.8162338159264
55.742762655163105
59.99920741694143
59.999207416941815
64.54457621139647
69.33509465228549
74.32462742825935

Has something in the nearest_neighbour implementation changed that could cause this?

I'm wondering whether you've benchmarked insertion performance when populating the R* Tree; We're currently using Spade in a couple of places in rust-geo (it's working perfectly, thank you!), and our primary applications involve one-time bulk-inserting line segments into the tree. It's hard to know without comparing to other libraries, but performance seems slower than expected; I see times around 90 ms for inserting segments built from ~9k points. If you have any ideas for improving this, or ideas for improving insertion performance in the library, I'd be interested in helping out with development.

Even though the project is set to passively maintained, I'm currently working on a large refactoring in the background. Here's an issue to discuss what Spade's current shortcomings are and how I plan to overcome them. Please, do feel free to start a discussion! That said, if you are interested in using spade, please don't be dismissed by this issue as it deliberately focuses only on the negative parts of Spade: It is still a feature rich and easy to use computational geometry library (let me know if you disagree!) that has proven itself in my personal projects. My experiences culminated in the ideas discussed here.

Current problems

Too many features

These features should probably be split into their own crates if required or be removed entirely

• Integer coordinates (along with the adaptive int kernel)
• f32 support (it is still possible to use f32 for storage, but it will be casted to f64 for calculation purposes)
• locate strategies. Those can probably be implemented on top of the triangulation.
• integration into cgmath and nalgebra, if the integration is kept, only feature gated and as a re-export
• farin interpolation
• Calculation Kernels should probably become an implementation detail, I don't see much sense for supporting the TrivialKernel at the moment.
• rtrees (These have already been extracted into the rstar crate)

Some missing features

• Bulk loading for triangulations
• Custom edge and face types
• Performance. Spade was not really focused on performance so far, I want to make this more of a priority.

Design changes

I believe that spade should be designed to be a "backend" for Delaunay related algorithms. Interpolation is a good example: Even today, interpolation algorithms can be implemented efficiently by using the current public API. Interpolation should probably be moved to a separate module or even its own crate. It may (I'm investigating that one) even be possible to implement the CDT based on the public Delaunay API. Finding the right design would allow to create custom tailored algorithms fitting many niches (e.g., alpha shapes).

Also, I plan to change an algorithmic technicality (although an influential one): The introduction of the "infinite point". This is an imaginary point "infinitely far away" from any other point, it is connected to every point of the convex hull. The main motivation is simplicity, I hope that the introduction will simplify the handling of the convex hull significantly and improving performance. However, I'm still unsure how and if this point should be part of the public API.

Open questions

• I'm not 100% sure if a 2.0 release or a split off as a new package is more apt, but I'm in favour of creating a new package (similar to the rstar outsourcing)
• Should the infinite point become an implementation detail or be part of the public API?
• If it becomes part of the public API, how should it be handled? Is it a regular vertex whose position is never evaluated? Can VertexHandle point to the infinite point?

Current Status

This reflects the status of my local prototype, I plan to upload it onto a branch once the dust settles a bit. The steps are roughly ordered according to when I want to work on them:

• [x] Removal of most of the above "unwanted features"
• [x] Locate in triangulation queries
• [x] Insertion
• [x] Infinite point
• [x] Vertex removal
• [x] CDT
• [x] Performance benchmarks
• [x] Documentation updates
• [x] Bulk loading

Final notes

Please, do feel free to open a discussion - there are many moving parts that need to be settled. If you even want to dive deeper, I'll happily answer questions or will mentor any kind of contribution.

Hello and thank you for all the hard work on Spade. It is a very useful library and I've had great success using it.

Context

I am using Spade to generate navigation meshes in a 2D video game. The inputs I give to Spade are the four corners of the playable space, and constraint edges representing all the polygonal obstacles within it. The resulting triangulation is a usable navigation mesh, PLUS unwanted faces within the obstacles.

My current approach to get rid of these unwanted faces is to inspect every face in the triangulation, and store which ones happen to be within obstacles. Whenever I make use the triangulation to compute paths, I simply ignore all these unwanted faces.

This image shows a playable space with red polygonal obstacles and the corresponding triangulation in cyan. Notice the many triangles overlapping the red obstacles:

This image shows the same triangulation, excluding the unwanted faces:

Sidenote: you may have noticed that the vertices in the triangulation are slightly offset from the obstacle vertices - that is intentional and not important for this discussion

Problem

The limitation of this approach is that calls to locate(), which are very useful to rein in stray game characters, have no way to ignore the extra faces.

I did manage to work around this by first projecting my arbitrary point onto the playable space without using Spade, and then calling locate() using that point. However, this is a fairly finicky and error-prone process.

Suggestion

In the past, I used Triangle to generate these navigation meshes. A very convenient feature of Triangle is that it supports the concept of 'holes', which can be used to eliminate unwanted faces from the output. The first example on this page illustrates the idea, with the inside of the CMU letters not being filled with faces:

An important subtlety is that holes can have holes themselves. In Triangle, this is implemented by specifying the triangulation input in three parts:

• Vertices (same as Spade)
• Segments (same as constraint edges in Spade)
• Holes (one point per hole)

Here is how they describe the implementation of the holes:

The third section lists holes (and concavities, if -c is selected) in the triangulation. Holes are specified by identifying a point inside each hole. After the triangulation is formed, Triangle creates holes by eating triangles, spreading out from each hole point until its progress is blocked by PSLG segments; you must be careful to enclose each hole in segments, or your whole triangulation might be eaten away. If the two triangles abutting a segment are eaten, the segment itself is also eaten. Do not place a hole directly on a segment; if you do, Triangle will choose one side of the segment arbitrarily.

Many thanks in advance!

I have to use nightly because of rocket. And I also need to use geo 0.10 which depends on spade 1.5.1. But spade fails to compile:

   Compiling spade v1.5.1
error[E0277]: the trait bound `num::BigInt: num::Signed` is not satisfied
  --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\traits.rs:39:6
   |
39 | impl SpadeNum for BigInt { }
   |      ^^^^^^^^ the trait `num::Signed` is not implemented for `num::BigInt`

error[E0277]: the trait bound `num::rational::Ratio<num::BigInt>: num::Signed` is not satisfied
  --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\traits.rs:40:6
   |
40 | impl SpadeNum for BigRational { }
   |      ^^^^^^^^ the trait `num::Signed` is not implemented for `num::rational::Ratio<num::BigInt>`

error[E0277]: the trait bound `num::rational::Ratio<bigvec::AdaptiveInt>: num::Signed` is not satisfied
  --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\traits.rs:42:6
   |
42 | impl SpadeNum for Ratio<AdaptiveInt> { }
   |      ^^^^^^^^ the trait `num::Signed` is not implemented for `num::rational::Ratio<bigvec::AdaptiveInt>`

error[E0277]: the trait bound `num::BigInt: num::Num` is not satisfied
  --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:99:10
   |
99 | impl <I> From<cg::Point2<I>> for BigVec2<BigInt> where I: cg::BaseNum + ToBigInt {
   |          ^^^^^^^^^^^^^^^^^^^ the trait `num::Num` is not implemented for `num::BigInt`
   |
note: required by `bigvec::BigVec2`
  --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:22:1
   |
22 | pub struct BigVec2<N: Num> {
   | ^^^^^^^^^^^^^^^^^^^^^^^^^^

error[E0277]: the trait bound `num::BigInt: num::Num` is not satisfied
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:105:10
    |
105 | impl <I> From<na::Point2<I>> for BigVec2<BigInt> where I: ToBigInt + na::Scalar {
    |          ^^^^^^^^^^^^^^^^^^^ the trait `num::Num` is not implemented for `num::BigInt`

    |
note: required by `bigvec::BigVec2`
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:22:1
    |
22  | pub struct BigVec2<N: Num> {
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^

error[E0277]: the trait bound `bigvec::AdaptiveInt: num_traits::Num` is not satisfied
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:149:6
    |
149 | impl Integer for AdaptiveInt {
    |      ^^^^^^^ the trait `num_traits::Num` is not implemented for `bigvec::AdaptiveInt`

error[E0277]: the trait bound `num::BigInt: num::Num` is not satisfied
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:100:5
    |
100 | /     fn from(v: cg::Point2<I>) -> Self {
101 | |         BigVec2::new(v[0].to_bigint().unwrap(), v[1].to_bigint().unwrap())
102 | |     }
    | |_____^ the trait `num::Num` is not implemented for `num::BigInt`
    |
note: required by `bigvec::BigVec2`
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:22:1
    |
22  | pub struct BigVec2<N: Num> {
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^

error[E0277]: the trait bound `num::BigInt: num::Num` is not satisfied
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:106:5
    |
106 | /     fn from(v: na::Point2<I>) -> Self {
107 | |         BigVec2::new(v[0].to_bigint().unwrap(), v[1].to_bigint().unwrap())
108 | |     }
    | |_____^ the trait `num::Num` is not implemented for `num::BigInt`
    |
note: required by `bigvec::BigVec2`
   --> C:\Users\me\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-1.5.1\src\bigvec.rs:22:1
    |
22  | pub struct BigVec2<N: Num> {
    | ^^^^^^^^^^^^^^^^^^^^^^^^^^

error: aborting due to 8 previous errors

For more information about this error, try `rustc --explain E0277`.
error: Could not compile `spade`.

rustc 1.30.0-nightly (90d36fb59 2018-09-13)

At first I got compilation errors in num-bigint and num-rational but i cargo updated those, but I can't update num to 0.2 because spade requires 0.1.

I want to perform interpolation on a sphere, with (theta, phi) coordinates.

I find that if I implement PointN for a SphPoint{theta:f64, phi:f64} then the SpatialObject is automatically implemented, then the distance2 is determined. However the distance of two points on a sphere is not calculated as theta^2+phi^2, but need to first convert into (x,y,z), then calculate the inner products.

How can I properly represent a point on a sphere?

I'm currently building on Spade to compute Voronoi diagrams for some research code, and I'm wondering if you'd be interested in getting a pull request to add this as a module. I understand that (c.f. #45) you want to avoid feature bloat, but this seems like it should be a fairly common use case of a Delaunay triangulation library (and small enough of an addition not to merit another crate).

The PR makes the following changes:

• change to PhantomData<fn() -> K> for storing Kernels
• added compile-time tests to verify common variants are Send+Sync
• fix minor lints

In addition, we may also tweak / duplicate the existing test_nearest_neighbor using crossbeam to test run-time multi-threaded behavior. Pls. let me know if that makes sense, and I can add that functionality to this PR

closes #42

While verifying the exact predicates implementation in Rust, I found that your incircle code produces results that differ from the original predicates code in C.

The problem lies in the following lines in exactpred.rs (lines 337-343):

        let mut aytcc = [0f64; 8];
        let aytcclen = scale_expansion_zeroelim(&cc, adytail, &mut aytcc);
        let temp16blen = scale_expansion_zeroelim(&aytcc[..aytcclen], cdx, &mut temp16b);

        let mut aytbb = [0f64; 8];
        let aytbblen = scale_expansion_zeroelim(&bb, adytail, &mut aytbb);
        let temp16clen = scale_expansion_zeroelim(&aytbb[..aytbblen], -bdx, &mut temp16c);

These lines should be changed as follows (variables named 'cc' have been swapped with 'bb' and vice versa):

        let mut aytbb = [0f64; 8];
        let aytbblen = scale_expansion_zeroelim(&bb, adytail, &mut aytbb);
        let temp16blen = scale_expansion_zeroelim(&aytbb[..aytbblen], cdx, &mut temp16b);

        let mut aytcc = [0f64; 8];
        let aytcclen = scale_expansion_zeroelim(&cc, adytail, &mut aytcc);
        let temp16clen = scale_expansion_zeroelim(&aytcc[..aytcclen], -bdx, &mut temp16c);

Note, the following are the original lines in predicates.c (lines 2839-2843):

    aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
    temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);

    aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
    temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);

After this change, the Rust incircle code produces identical results to the original C code.

I'm currently working on splitting off the r-tree into its own package. One of Rusts strengths has always been Cargo, which allows to have many, small packages with a single use. R-Trees would actually be perfect in such a scenario.

I would also be looking to improve performance and to implement some more advanced algorithms - I still have some ideas for this in my mind!

I'm still looking for a good name, just in case anyone has a good idea.

Can you tell us something about the quality of your CDT algorithm (speed and robustness) compared to CGAL lib, and maybe which papers you use to implement it?

Origin for my question is, that I would need a CDT for Nim language, and I would prefer a native implementation to restricted CGAL bindings (I made some CGAL bindings for Ruby some years ago.) Of course coding a robust and fast incremental constrained delaunay triangulation from scratch using research papers is a big task, so I am considering converting an existing implementation. Maybe from C++ or Rust. I have not yet found a nice standalone C++ implementation -- of course CGAL has one, but I guess extracting CDT code from lib is not easy. Your code is MIT licensed, same as Nim libs generally are, so it may a good candidate. I would have to learn some Rust before, but that is already on my todo list. (Funny fact is, that some weeks before I saw your spade package, I did a Nim RTree implementation myself from scratch following the original papers, see https://github.com/StefanSalewski/RTree.)

I'm getting an error for code that I believe should work. Here's my minimal reproduction:

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

fn main() {
    let vertices = vec![
        Point2::new(60.0, 20.0),
        Point2::new(84.0, 20.0),
        Point2::new(236.0, 20.0),
        Point2::new(60.0, 36.0),
        Point2::new(84.0, 36.0),
        Point2::new(92.0, 28.0),
        Point2::new(60.0, 44.0),
        Point2::new(92.0, 44.0),
        Point2::new(212.0, 44.0),
        Point2::new(20.0, 60.0),
        Point2::new(156.0, 68.0),
    ];

    let edges = vec![(1, 2), (0, 3)];

    let mut triangulation =
        ConstrainedDelaunayTriangulation::<Point2<f32>>::bulk_load(vertices.clone()).unwrap();
    let vertex_handles = triangulation.fixed_vertices().collect::<Vec<_>>();

    for (from, to) in edges {
        println!(
            "{:?} ({:?}), {:?} ({:?})",
            vertices[from], vertex_handles[from], vertices[to], vertex_handles[to]
        );
        triangulation.add_constraint(vertex_handles[from], vertex_handles[to]);
    }
}

Point2 { x: 84.0, y: 20.0 } (FixedHandle { index: 1 }), Point2 { x: 236.0, y: 20.0 } (FixedHandle { index: 2 })
Point2 { x: 60.0, y: 20.0 } (FixedHandle { index: 0 }), Point2 { x: 60.0, y: 36.0 } (FixedHandle { index: 3 })
thread 'main' panicked at 'Error - constraint edges must not intersect each other', C:\Users\Hayde\.cargo\registry\src\github.com-1ecc6299db9ec823\spade-2.0.0\src\cdt.rs:363:17 

I've removed as many points and edges as I can while preserving the error. Note that the edges don't actually intersect: Thanks!

Spade 2.0 looks great API-wise, and is almost perfect for my triangulation needs.

One thing I found missing is that, when iterating over inner faces, I need to get their 3 neighbour faces, which currently seems pretty difficult.

Would it be possible to add an adjacent_faces method similar to adjacent_edges to the face handle?

Hi!

I needed to do CDT where I define the outer contour as a constraint (as well as zero or more inner cutouts). Adding constraints along the outer vertices didn't help with this though, as the CDT implementation in spade will find convex-hull triangles outside my outer contour if they exist.

My solution thus far has been to check all triangle-centers using a vector_inside_polygon_test(..), and throwing away any triangle which is outside my outer contour.

Is there a better way of doing this using the library?

A lot of the time, when you call intersects_constraint, you're interested in which constraint you're intersecting, e.g. so you can get rid of it, or alter it, etc.

Do you think it would make sense for intersects_constraint to return Option<DirectedEdgeHandle> or something?

Can write a patch if you are interested.

You sometimes represent a mesh with an array of vertices, and an array of triangle indices. I would construct such an array with something like:

let verts = triangulation.vertices().map(|pt| pt.data).collect();
let tris = triangulation.inner_faces().map(|f| f.vertices()).flat_map(|v| v.index());

Except, in spade 2.0, the index method is private. Is there another way?

See #66.

This is the initial implementation which works sufficiently well already.

To use this, create a CDT, insert some constraint edges and then call

my_cdt.refine(RefinementParameters::new().exclude_outer_faces())

which will refine the cdt inplace. Note that, in order for exclude_outer_faces to work, the constraint edges should form a loop somewhere - otherwise, all faces will be considered to be outer faces and won't be refined.

Remaining TODOs (just mentioning for my own overview):

• [x] Documentation of AngleLimit
• [x] Documentation of RefinementParameters
• [x] Documentation of fn refine
• [x] Internal documentation by adding a few comments at appropriate places
• [ ] At least some improved test coverage, e.g. for testing special cases
• [ ] Expand fuzz tests (those are more or less just stubs at the moment)

TODOs for some follow-up PRs:

• Try to enable refinement algorithm for regular triangulations (if easy)
• Investigate if the refinement quality can be increased (though it's not too bad already)
• Add some images detailing how exclude_outer_faces works