Closed rmanoka closed 1 year ago
Nice work! I'd never heard of this concept before. The algorithm part is kind a bit dense, so I don't feel like I adequately reviewed it.
My use-case was a generalization of polygon rasterization, but I felt it's a useful algorithm to abstract out. I was getting performance comparable to gdal_rasterize for polygons. The algorithm essentially traces along edges, but yeah, the details are a bit dense; I'll try to flesh out some more docs in-code.
For validation, it looks like you're checking that the area of several input polygons is equal to it's monotonic decomposition, is that right?
Yeah, since it's integer arithmetic, we just check equality. We could add Relate
based assertions too: all components are within original polygon, all components are interior disjoint.
via a wrapper around a Vec
. Do you think we should introduce something like this so we can hang impls on it?
struct MonotonicPolygons(Vec<MonoPoly>);
As well as saving people from writing boilerplate, it seems like it'd make it more obvious to explain to people how to leverage this functionality, without necessarily having to explain the entire concept behind what a monotonic polygon is.
+ let mono_polys = MonotonicPolygons::from(polygon); for candidate in polygons { - polygon.intersects(&candidate) + mono_polys.intersects(&candidate) }
Yes absolutely, this makes sense; the for loop would move into the impl.
That part could always follow in another PR though.
Shall we merge this @michaelkirk @frewsxcv ? I've some bool ops ideas on top of this PR, so would be nice to have this in main.
Thanks for the reviews! Merging this.
bors r+
Build succeeded!
The publicly hosted instance of bors-ng is deprecated and will go away soon.
If you want to self-host your own instance, instructions are here. For more help, visit the forum.
If you want to switch to GitHub's built-in merge queue, visit their help page.
[x] I added an entry to
CHANGES.md
if knowledge of this change could be valuable to users.This PR provides a pre-processed representation for a polygon (port from my
geo-crossings
unpublished crate), that allows fast point-in-polygon queries. This implementation is based on these awesome lecture notes by David Mount. The results are promising: the pre-processed queries are faster across the board, about 1.5-2x in the worst-case, to upto 100x in the best cases.In this PR, we provide:
A
MonoPoly<T>
structure that represents a monotone polygon, supporting fast point-in-polygon queries.A builder algorithm:
monotone_subdivision(Polygon<T>)
that sub-divides a polygon into a collection ofMonoPoly<T>
s.Benchmarks that show that the pre-processed point-in-polygon queries are fast.
Benchmarks
We use three types of polygons:
A steppy polygon that zig-zags along the X-axis as it moves up along Y-axis, then again zig-zags down the Y-axis. This is the worst case for our algorithm as it gives many monotone polygons along the X-axis.
Same as above but along Y-axis. This is a good case as it is representible by a single monotone polygon.
Circular polygons. These also give many monotone-polygons on sub-division, and are a good test case.
With each polygon class, we generate polygons with different number of vertices, and measure query times for checking if a random point in the bounding box lies in the polygon.
Try out with
cargo criterion --bench monotone_subdiv
.Remarks / Caveats
This uses a simpler version of sweep used for the boolean ops. I've re-implemented the sweep to try to move away from the complicated float issues we've been facing. For instance, this whole algorithm works without requiring
T: Float
.The algorithm currently only supports a single polygon, and no two segments in it can intersect anywhere except at their endpoints. This is almost the SFS spec, but I think SFS also allows intersections at interior that are "tangential". I believe this can be handled without involving float ops as the intersection is "exact", but needs more work.
This data-structure is currently not a stand-in replacement for a
Polygon
forCoordinatePosition
algorithm as we don't currently remember the new edges added during subdivision. These edges actually lie inside the polygon, and points on these edges would be on the border of two of the sub-dividedMonoPoly
. This could be handled by some clevermod 2
counting but I wasn't sure of the logic. As it stands, we can get fast equivalentIntersects
implementation with this impl via a wrapper around aVec<MonoPoly<T>>
.