Add self intersection for polygon and polyline

[mpizenberg]

Computing self intersection for a polygon or a polyline is not obvious, even if an intersection function for line segments is provided. I recently had to do this so I wonder if that's something that would be interesting for OpenSolid. I have implemented a function with this interface:

selfIntersection : Polygon2d -> Maybe Point2d

It "stops" as soon as one intersection point is found. This was because I only wanted to know if the polygon was simple or not. I used (almost) the Shamos-Hoey algorithm which is a simplification of Bentley-Ottmann intersection algorithm stopping at first intersection.

Let me know if that's something that you'd like to add to OpenSolid.

[ianmackenzie]

It would certainly be a useful addition, although I'd like to take some time to understand the implementation a bit (especially stuff like numerical stability). Can you share your current code and/or send me a link to a description of the algorithm?

Signature-wise, I think

selfIntersection : Polygon2d -> Maybe Point2d

could be slightly misleading since it kind of implies the result is "the" self-intersection if it exists; perhaps

findSelfIntersection : Polygon2d -> Maybe Point2d

since I think there's a fair bit of precedent for 'find' functions returning only the first result found. Or, for maximum flexibility (find whatever you want, without paying for anything you don't care about):

hasSelfIntersection : Polygon2d -> Bool
findSelfIntersection : Polygon2d -> Maybe Point2d
allSelfIntersections : Polygon2d -> List Point2d

[mpizenberg]

Yes findSelfIntersection is way better ;)

The Bentley - Ottmann algorithm is also called the sweep line algorithm. You can find two interesting explanations at these references: http://www.webcitation.org/6ahkPQIsN http://jeffe.cs.illinois.edu/teaching/373/notes/x06-sweepline.pdf

My current implementation is accessible in file Polygon2d.elm. Please be aware of the following remarks:

• This package was initially designed to be available on elm package an was my second ever created elm package. Well, it should not have ended on elm package but it was to try the adventure ...
• I used the Pivot data structure for the sweep line. It does not have a sorted insertion and search in O log(n). n is here the number of segments currently crossing the sweep line. But since in my case n would never be large (usually < 6) it didn't really matter. For a polyline, where lines are usually well less organized than with regular polygones, n could be higher and this could be problematic.

To see a working use of the intersection algorithm, you can try running the examples/protoFeedback.elm file.

[ianmackenzie]

What are your thoughts on the various numerical issues described at https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm#Special_position_and_numerical_precision_issues? Polylines/polygons will certainly often have vertical segments and by definition will have segments with coincident endpoints. I think that findSelfIntersection would be a useful function, but I'd like its behaviour to be consistent with a brute-force 'iterate over every pair of segments and check that pair using LineSegment2d.intersection' (that is, you should never have a case where one of those methods finds an intersection point and the other one doesn't). Fortunately, fuzz testing with elm-test makes that sort of thing pretty easy to check =)

An alternative would be to build an R-tree of line segments and then query that R-tree with each individual line segment to look for possibly-intersecting other segments. This wouldn't be quite as fast as a specialized algorithm like Shamos-Hoey but should be reasonably efficient. I'm planning on building an OpenSolid R-tree package anyways at some point since it's a very useful data structure for a lot of cases.

[mpizenberg]

What are your thoughts on the various numerical issues described at https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm#Special_position_and_numerical_precision_issues? Polylines/polygons will certainly often have vertical segments and by definition will have segments with coincident endpoints.

The listed restrictions related to numerical issues that my happen are:

1. No two line segment endpoints or crossings have the same x-coordinate
2. No line segment endpoint lies upon another line segment
3. No three line segments intersect at a single point

The findSelfIntersection function only cares about finding at least one intersection point. Therefore, it uses the Shamos–Hoey algorithm and so issue 3 is not an issue anymore.

Issue 1 is basically a problem when we have coincident endpoints. However, the implementation I chose ignores successive segments in the polygon by keeping their index in the polygon. So it will only detect "true" intersection, not successive segments intersections.

I don't exactly know what the issue with issue 2 is (except for pure numerical instability). However, I feel that it will not be an issue with this implementation for the following reasons:

• if it is a >| configuration. Then I guess intersection or no intersection is just a matter of perception.
• if it is a -|- configuration (where 2 successive segments cross another). I'm confident that either the first or the second segment (at least, maybe both) will raise an intersection. If both are detected as an intersection, same as 3, we don't care. If however neither is detected, that's an issue. In practice this case should be fuzzy tested.

I'd like its behaviour to be consistent with a brute-force 'iterate over every pair of segments and check that pair using LineSegment2d.intersection'

Segments are not modified by the algorithm, and internally, each intersection computation is done with that very function. So the only reason for not having the same result would be if some pair comparison is not done. In this algorithm, pair comparisons are done only between "adjacent" segments in the sweep line data structure. Therefore a difference between those two algorithms would happen only if all "brute" intersection points are never appear as adjacent segments in the data structure. However the construction of the algorithm "should" avoid such an issue. This may also deserve some fuzzy testing.

If you can prepare some tests (I might be biased by my implementation), I can implement this and make sure it has a valid behavior.

I'm planning on building an OpenSolid R-tree package anyways at some point since it's a very useful data structure for a lot of cases.

Sounds cool :)