Create Intersection algorithm trait

[frewsxcv]

It might be the case (in the case of intersecting two polygons) that the operation will return a polygon or a linestring or a point or nothing, the return type will need to reflect this. Maybe something like this?

trait Intersection<...> {
    fn intersection(&self, other) -> Option<IntersectionResult> {
        ...
    }
}

// not a fan of "Result" here since it sounds similar to std::result::Result
enum IntersectionResult {   
    Polygon(Polygon),
    LineString(LineString),
    Point(Point),
}

Note that this is different from the Intersects trait which simply returns a bool indicating if two geometries intersect at any point

[ambaxter]

Yup. That would be perfect.

[mbattifarano]

I think Vatti's Algorithm would work well here. My understanding is that it can be parameterized to produce union, intersection, xor and difference.

https://en.wikipedia.org/wiki/Vatti_clipping_algorithm
extras.springer.com/2005/978-1-84628-108-2/VattiClip.pdf

[shterrett]

@mbattifarano and I were looking into this the other day. I am planning on implementing the algorithm that is outlined in de Berg, Cheong, van Kreveld, Overmars Computational Geometry. Briefly

• Represent the polygons in a doubly-linked edge list so that each edge has a reference to its successor and predecessor, and each face has a reference to an incident edge
• Combine the two polygons into a single list, ignoring edge intersections
• Use a scanning line algorithm to detect all of the intersections, and add a node at each point, aligning the references as necessary, including a reference to the originating polygon for each face (which may be both)
• Use boolean operations on the faces + originating polygons to keep only those desired (intersection/union/etc...)

This calls for a few new data structures:

• Doubly linked polygon edge list
• Line segment (only two points; not a LineString)
• Priority Queue that allows testing for containment and maintains order of the "event points" by upper-left-ness
• Binary search tree that returns the adjacent segments to given segment at any point

In addition, we've run in to one very concrete problem regarding Float. To use a BTreeSet as the queue, we need to order Float, which is impossible because of NaN. Our current thinking is to simply let the library panic if NaN is encountered (or infinity), but I'm not sure this is the right way to go.

Also, there isn't a binary search tree in the stdlib as far as I can tell, so I'm not sure what the best data thing to use there is. We can start with a simple vector, which isn't the best performance. If there's a binary search tree create, we can pull that in too.

Anyway - tl;dr - this is shaping up to be a huge change, and I wanted some guidance before writing a ton of code.

[urschrei]

I implemented a priority queue using BinaryHeap in distance.rs if it's any use – if I recall, I came to the same conclusion regarding NaN; intuitively, it feels like it shouldn't be possible to construct any of the geometry primitives with NaN values (I'm not sure whether GEOS allows them), though that's a separate discussion.

[shterrett]

@urschrei The book says that the priority queue cannot be implemented using a BinaryHeap because of a need to test for inclusion; however one option is to keep a secondary Set, and use it as a proxy for inclusion. Extra memory though.

I played with something similar for the Ord instance. I like yours, and I assume using the single-member struct is better form than impl-ing directly on Float?

[urschrei]

I've just come across @notriddle's QuickerSort crate, which includes a well-tested non-allocating float sort, which could be used in an Ord implementation on Point (lowest X followed by lowest Y?), which would give us a BTreeSet. Any thoughts, @shterrett & @mbattifarano?

[mbattifarano]

Hey, sorry for the radio silence, I've been really busy for the past few months but now have some free time on my hands and would love to start working on this again.

I went ahead and implemented a line segment struct (PR coming soon) as a first step. In addition to being useful for the algorithm @shterrett outlined above, I think it would be a helpful abstraction for several existing LineString algorithms -- basically anywhere there is a loop over linestring.0.windows(2) (here for example).

Over the next few days, I'll take a look at getting a priority queue working.

[notriddle]

BTW, instead of using quickersort, you'll want to use float-ord, which builds on top of the faster pdqsort crate.

[urschrei]

Amazing. Also, thanks for the heads-up, @notriddle. I'm currently stuck for time too, but I'll take a look at the segment PR – looks great at first glance. Also, it occurs to me that if you implement a line-scanning algorithm (Bentley-Ottmann?) we'll also get a check for self-intersection for free \o/

[Dushistov]

There is clipping crate with intersection, union, difference operations on polygons Greiner-Hormann algorithm. What do you think about porting it to usage of rust-geo types?

[urschrei]

Hmm interesting. It looks quite simple, once you have the points in a suitable structure (is it a DCEL?). I note that it can't deal with degeneracies (according to WP: common edges or intersections exactly at a vertex. It's not clear whether it can deal with holes), which limits practical use a bit.
On the subject of DCEL: it was a major missing component for us, because Vatti needs it, but Spade provides a really implementation, even though it doesn't have a public interface atm, and we've already impl'd the traits for points (and the line segment impl will merge at some point in the not-too-distant future).

[Indy2222]

It might be the case (in the case of intersecting two polygons) that the operation will return a polygon or a linestring or a point or nothing, the return type will need to reflect this.

Intersection of two polygons may result in any combination of the following geometries (there might be multiple instances of each in the result):

• Point, MultiPoint
• Line, LineString, MultiLineString
• Polygon, MultiPolygon

Popular Python geometry package Shapely represents such results as a GeometryCollection (or simpler objects like Point if they sufficiently represents the intersection). So IntersectionResult must have an option for basically any geometry type similary to types::Geometry (https://github.com/georust/rust-geo/blob/master/src/types.rs#L877)` why not to reuse the enum?

In the following example intersection of green and polygons consists of one point and two polygons (black areas).

[urschrei]

I came across the Martinez-Rueda Polygon Clipping Algorithm the other day. It looks like a faster replacement for Vatti, and there are some extant JS implementations, and even a Rust version (which isn't fully working, but got quite far).

It may well be worth exploring whether we could implement this without resorting to a load of unsafe with @fschutt's help.

[fschutt]

@urschrei There is a full Rust implementation of the same algorithm that someone used for clipping DOOM WAD files: https://eev.ee/blog/2018/03/30/a-geometric-rust-adventure/ - code here: https://github.com/eevee/idchoppers/blob/master/src/shapeops.rs - it might be better to start with that.

But I'd first just get it working before talking about the API design. Refactoring is the easy part here, the hard part is ownership, numerical stability and efficiency. For GIS operations, Vatti is AFAIK slower than Martinez-Rueda, see this paper: http://www.cs.ucr.edu/~vbz/cs230papers/martinez_boolean.pdf

[urschrei]

@fschutt Ohh yeah I remember looking at eevee's version a few months ago, and being a bit intimidated by the amount of it, and lots of comments like:

// hey what's up, this is a catastrophic mess ported from a paper

But maybe it's worth figuring out what we could use, and don't need to port, and what we can do without. There's no sweep-line functionality (let alone anything that deals with collinear segments) in Geo, and it seems to be really difficult to get right, so that would probably be a priority.
I read the Martinez paper yesterday, and it was the efficiency improvement over Vatti that got me interested. The problem is that there a quite a few high-quality Vatti implementations across several languages, and almost no Martinez-Rueda implementations…

[untoldwind]

FYI: I needed this feature and spend some time re-implementing the JavaScript version (https://github.com/w8r/martinez) to rust. There result could be found here: https://github.com/21re/rust-geo-booleanop

There certainly is still a bit of refactoring required (and potentially some optimization), but essentially it is producing the same results ... apart from some minor rounding issues for the extreme cases.

[fschutt]

@untoldwind Would you mind if I create a seperate version of your crate without the geo dependencies? It's just that I wouldn't want to pull in the entire geo library to intersect two polygons (esp. when these polygons are non-geographic, so I don't need the geo crate). Of course, with proper attribution and same license. Or you could do it and make a geo binding crate. But for example, in order to create intersections in a vector editor, I don't need any geo algorithms other than boolean operations, so for me it'd be kinda pointless having to compile the entire geo library just to get boolean operations.

[untoldwind]

Good point, I've reduced the dependency to "geo-types".
"geo" and "geojson" are still required for development/testing.

If you want/need a different implementation of Polygon/Coordindate/..., well this needs some refactoring.

As for the License: It's MIT like the original JavaScript implementation, so: Sure go ahead ;)

Be aware though that the current version is still somewhat basic. More specifically: I'm not very happy about the use of Rc and Weak and the RefCell, which might lead to runtime panics. I'd really like to make these obsolete somehow ...

[fschutt]

@untoldwind If you can make serde an optional feature and disable spade, then the only dependency you have is num-traits and I don't think a fork would be worth it in that case. geo-types depends on spade, which has a huge number of dependencies. However, spade can be turned off, since it's an optional feature. The library compiles fine without it, so I think spade isn't actually used or needed. The difference in dependencies is pretty huge:

geo-booleanop v0.1.1
├── geo-types v0.2.0
│   └── num-traits v0.2.6
└── num-traits v0.2.6 (*)

vs:

geo-booleanop v0.1.1
├── geo-types v0.2.0
│   ├── num-traits v0.2.6
│   └── spade v1.5.1
│       ├── cgmath v0.16.1
│       │   ├── approx v0.1.1
│       │   ├── num-traits v0.1.43
│       │   │   └── num-traits v0.2.6 (*)
│       │   └── rand v0.4.3
│       │       └── libc v0.2.43
│       ├── clamp v0.1.0
│       ├── nalgebra v0.15.3
│       │   ├── alga v0.6.0
│       │   │   ├── approx v0.2.1
│       │   │   │   └── num-traits v0.2.6 (*)
│       │   │   ├── num-complex v0.2.1
│       │   │   │   └── num-traits v0.2.6 (*)
│       │   │   └── num-traits v0.2.6 (*)
│       │   ├── approx v0.2.1 (*)
│       │   ├── generic-array v0.11.1
│       │   │   └── typenum v1.10.0
│       │   ├── matrixmultiply v0.1.14
│       │   │   └── rawpointer v0.1.0
│       │   ├── num-complex v0.2.1 (*)
│       │   ├── num-traits v0.2.6 (*)
│       │   ├── rand v0.5.5
│       │   │   ├── libc v0.2.43 (*)
│       │   │   └── rand_core v0.2.2
│       │   │       └── rand_core v0.3.0
│       │   └── typenum v1.10.0 (*)
│       ├── num v0.1.42
│       │   ├── num-bigint v0.1.44
│       │   │   ├── num-integer v0.1.39
│       │   │   │   └── num-traits v0.2.6 (*)
│       │   │   ├── num-traits v0.2.6 (*)
│       │   │   ├── rand v0.4.3 (*)
│       │   │   └── rustc-serialize v0.3.24
│       │   │   [dev-dependencies]
│       │   │   └── rand v0.4.3 (*)
│       │   ├── num-complex v0.1.43
│       │   │   ├── num-traits v0.2.6 (*)
│       │   │   └── rustc-serialize v0.3.24 (*)
│       │   ├── num-integer v0.1.39 (*)
│       │   ├── num-iter v0.1.37
│       │   │   ├── num-integer v0.1.39 (*)
│       │   │   └── num-traits v0.2.6 (*)
│       │   ├── num-rational v0.1.42
│       │   │   ├── num-bigint v0.1.44 (*)
│       │   │   ├── num-integer v0.1.39 (*)
│       │   │   ├── num-traits v0.2.6 (*)
│       │   │   └── rustc-serialize v0.3.24 (*)
│       │   └── num-traits v0.2.6 (*)
│       ├── pdqselect v0.1.0
│       └── smallvec v0.6.5
│           └── unreachable v1.0.0
│               └── void v1.0.2
│       [dev-dependencies]
│       └── cgmath v0.16.1 (*)
└── num-traits v0.2.6 (*)

However, I have to note that there doesn't seem to be a way to turn off spade via cargo yet, because there is no feature in the Cargo.toml for it (here).

[untoldwind]

I'd love to remove spade, but as noted there is no (obvious) feature to disable.

On the bright side: The core algorithm just requires geo_types::Coordinate and geo_types::Rect, both are pretty easy to replace. So separating the algorithm from the geo-binding should not be that big of a deal.

[fschutt]

@untoldwind yeah, I opened #307 because of this, once that's merged and re-published it should be possible to disable spade. In my demo above I did it via { git = "https://github.com/georust/geo", default-features = false } - for some reason, this works and disables spade, but { version = "0.2", default-features = false } doesn't.

[fschutt]

@untoldwind You can just turn off the unnecessary dependencies by doing this:

[dependencies]
geo-types = { version = "0.2.0", default-features = false }

That should turn off spade + serde.

[untoldwind]

Did that, but the dependency is still there. Somehow the dev-dependencies (which implimicitly require serde on geo-types) seem to interfere with the regular dependencies.

Maybe this could be fixed by moving all the tests to a sub-package.

[untoldwind]

Edit: I split the project into lib and tests, whereas tests contains all the tests requiring geojson. This should eventually cleanup the dependencies.

[frewsxcv]

FYI: I needed this feature and spend some time re-implementing the JavaScript version (https://github.com/w8r/martinez) to rust. There result could be found here: https://github.com/21re/rust-geo-booleanop

@untoldwind Are you planning to publish your crate to crates.io? I'm in favor of adding your crate as a dependency to geo and exposing that functionality. Looks like a great crate btw!

[untoldwind]

Was on vacation and somehow found no time to work on this :)

Anyhow: There is now https://crates.io/crates/geo-booleanop

[bluenote10]

@untoldwind Nice work! Unfortunately the reference implementation of the Martinez algorithm seems to have quite a lot of hard to solve issues (unmaintained?). So perhaps it's better to go for Vatti, because of easier handling of corner cases? Edit: I just came across another JS implementation of the Martinez algorithm -- maybe it's worth checking if this implementation has solved some of these issues...

In the C# world, a famous library that does this is the Clipper library -- old, but to my knowledge fairly stable (Vatti based). Maybe its implementation in C++ or C# could be used as inspiration. Some benchmarks (quite old, but probably still relevant):

• From the clipper docs
• Independent benchmark from 2011.

[mthh]

There is clipping crate with intersection, union, difference operations on polygons Greiner-Hormann algorithm. What do you think about porting it to usage of rust-geo types?

Hmm interesting. It looks quite simple, once you have the points in a suitable structure (is it a DCEL?). I note that it can't deal with degeneracies (according to WP: common edges or intersections exactly at a vertex. It's not clear whether it can deal with holes), which limits practical use a bit.

There has been a recent proposal in the scientific literature to extend the Greiner-Hormann algorithm in order to manage degenerate intersections (vertex on edge, overlapping edges, etc.) : https://doi.org/10.1016/j.cagx.2019.100007 - the paper is accompanied by a c++ implementation and data which allow the examples of the paper to be reproduced.
It seems a limitation subsists in their proposal ("if a self-intersection of one polygon lies on an edge or coincides with a vertex of the other polygon, then the algorithm may fail", discussed in part 6. Discussion and conclusions) but it could still be interesting !

[urschrei]

It seems a limitation subsists in their proposal ("if a self-intersection of one polygon lies on an edge or coincides with a vertex of the other polygon, then the algorithm may fail", discussed in part 6.

Doesn't this directly imply that the polygon is invalid (in the DE-9IM sense) anyway, so we wouldn't have to worry about it…

[mthh]

Yes, that's why I thought that the approach proposed in this paper might be a good candidate for the case discussed here (however I didn't investigate if it will be difficult to implement compared to the original algorithm).

[urschrei]

Ah, sorry!

[bluenote10]

Another aspect that comes to mind regarding the algorithmic interface:

When I was last touching this topic many years ago, my use case was to have a large number of polygons each with a low vertex count, which I wanted to union(*). I can't remember all the details of my benchmarks of various libraries at that time, but what seemed to make a rather a big difference is whether a library supports only pair-wise operations, or an n-way operation. With pair-wise boolean operations one has to fold all input polygons, and apparently repeating the algorithmic overhead can have a big effect. If I remember correctly, Angus Johnson's Clipper library had this feature, which might be why I kept it in good memory.

(*) in fact that's also my use case this time, so I would argue it's a fairly common use case ;).

[agersant]

@bluenote10 Your memories are correct, here is an example of what the API looks like for this in Clipper: https://stackoverflow.com/a/7548508

[velvia]

Hi folks, what is the latest on this issue? It was first opened in 2016, but other than the third party crate there is no purely Rust intersection / boolean op available, which is unfortunate.

[urschrei]

The current best solution is https://github.com/21re/rust-geo-booleanop, which may or may not work for your use case.

[velvia]

@urschrei I find that when I use the BooleanOp trait somehow rustc/cargo cannot find the intersection method. Very puzzling.

[lnicola]

That's usually a version mismatch (one crate implements a trait for a struct defined in an older version of a crate you are using).

[frewsxcv]

Very relevant to this is this GitHub issue which is a discussion about incorporating geo-booleanop algorithms directly into geo

[urschrei]

Completed by #835!