elm-geometry icon indicating copy to clipboard operation
elm-geometry copied to clipboard

Simple polygon intersection

Open gampleman opened this issue 5 years ago • 2 comments

This is an odd one.

I've built this the simplest, dumbest way possible just to get something out of the door (I happen to need this for my elm-europe talk). It is most certainly not state-of-the-art algorithmically. Nor do I claim it to handle degenerate cases. It seems to do a good enough job for my kind of use-case (i.e. reasonably complex polygons without holes or kinks).

So I can foresee some things you might want to do with this PR:

  • Close it and pretend it never happened.
  • Ship it (after some of the usual tidy up) until we can build something better.
  • Keep it around for some of the following uses:
    • an oracle for fuzz testing
    • a baseline for benchmarking
    • scrapping it for parts (i.e. use the subdivide algorithm, but use something better for actually assembling the polygons, like Rivero & Feito, 2000).

Some details on how this works.

The algorithm works in 2 phases. First there is the subdivide phase. In this phase extra vertices are added into the intersections of the edges of the two polygons. This now means we have polygons with more edges, but we now have all the edges for the resulting polygon. I'm not entirely sure about the runtime here, but it must be at the minimum O(n*m) (where n is the number of vertices in the first polygon, and m is the number of vertices in the second polygon; so effectively quadratic performance).

Next we need to actually assemble the intersected polygon. We go through each edge in both the subdivided polygons, and check if the midpoint of that edge is in the other polygon. We then attempt to order the edges into proper loops. (Note: at the moment I haven't built support for holes, but conceptually this shouldn't be too hard - we could run buildChain until no more edges remain in edgeDict, then simply figuring out which contains which). This operation is at least O(n'*m + m'*n), where the primes are the vertex counts after subdivision.

gampleman avatar Jun 12 '19 09:06 gampleman

Thanks for the PR! I'd be very nervous about publishing this as-is, though, since I know that polygon intersection is a very tricky thing to get right and will probably require a pretty sophisticated algorithm to be appropriate for general-purpose use.

I've also tried very hard to avoid having any hidden 'epsilon' or 'tolerance' values in elm-geometry code, since it's usually impossible to find a single value that works for all cases: some people might be working in pixels, some in miles, some in nanometers. So things like Triangle2d.area (...) < 0.00001 and the funky rounding in toKey make me suspect that this is tuned for the kinds of coordinate value magnitudes that you see in your work, and could easily break down for people working at different scales. Passing in a tolerance value is an option, but I've found that even that is pretty error-prone - I try to keep the use of tolerances at all to a minimum, and only in places where the meaning of the tolerance value is very well-defined and geometrically meaningful (like Point2d.equalWithin).

I do agree that there's a lot of useful code in here, though - what would you think of moving the code into a non-exposed Polygon2d.Intersection module and renaming the intersection function to simple? Then we can call it from the tests, experiment with it from the REPL etc. without any additions to the public API. And potentially in the future we can run tests/benchmarks comparing Intersection.simple against Intersection.weilerAtherton or whatever.

(Side note: not sure what's going on with Travis, seems like it might be an issue with the Elm NPM package? If it doesn't go away in a few days I will investigate further...)

ianmackenzie avatar Jun 13 '19 19:06 ianmackenzie

Makes sense, will do in a couple of weeks.

Regarding the tolerances: It's possible neither may be necessary. The area check is just checking if all three points lie on a single line - there must be a better way to compute this than checking if the area of the triangle formed is zero. The distance filter is there just to remove extra zero-ish length edges, which is not necessary, but leads to simpler results.

It's possible the rounding in toKey may not be necessary anymore, I'll investigate that later (one of those things you add trying to get it to work and it stays on in the code).

gampleman avatar Jun 14 '19 10:06 gampleman