Greiner–Hormann O(nm)

[mrehayden1]

I'm also thinking about tackling this too.

It's less general than the Vatti algorithm in that it can't handle holes, therefore SimplePolygon p r -> SimplePolygon p r -> [SimplePolygon p r] seems like the best fit I could come up with for a type signature.

GH can also handle self intersecting polygons but as I understand, that's not excluded by SimplePolygon, just unadvised.

What do you think?

[mrehayden1]

Oh, and maybe there should be parameters for the set theoretic operations to be performed between the two polygons, or better still separate functions.

Any suggestions on what to call them?

[noinia]

Some collection of remarks related to set operations on polygons:

• No matter what algorithm we implement, one question is how to represent the result. Intersecting two polygons does not necessarily result in one polygon (as you already observed). The 0.X series of hgeometry had a "PlanarSubdivision" type that could represent embedded planar graph graphs consisting of multiple connected components. I think we will want to include that into hgeometry again (probably together with some PlanarSubdivison_ typeclass, consistent with the other typeclasses).

  Given a "red" polygon, and a "blue" polygon, we could then have some functions of type roughly:

  overlay :: redPolygon -> bluePolygon -> PlanarSubdivision s v e FaceType

  in which FaceType is some enum with options: Outside, OnlyRed, OnlyBlue, RedAndBlue. (i.e. every face of the subdivision is tagged as outside, red, blue, or both). We could similarly tag the edges actually.

  To compute the actual intersection, we would keep only the faces with type 'both'. To compute the union we could keep all but Outside etc.

  Bringing back the Planarsubdivision type should not be too difficult I think. Most of the annoying bookkeeping code I have already written before :).

• Implementing the overlay (and/or intersections) function itself is what could be done using the Greiner–Hormann algorithm as you suggested, but I think there are also two other (maybe easier/better?) options:

  1. computing the overlay means computing the intersections; that is already done in the LineSegmentIntersection algorithm. What remains is "only" the face information. That should be doable using another sweep line very similar to the one in VerticalRayShooting implementation.

     In total this would give an O((n+m+k)\log (n+m)) time algorithm, where k is the output complexity. In the worst case that is indeed slightly slower than the O(nm) time algorithm from Greiner–Hormann, but for many practical applications k is probably much lower than nm. So I feel that the two sentences in their paper don't really do that algorithm justice.

  2. Furthermore, if you really do want a fast algorithm, even in the worst case, then Chan's bichromatic intersections algorithm actually runs in optimal O((n+m)\log(n+m) + k) time. (That is for computing the overlay, one would still have to be a bit careful in computing the face information I think). That again actually seems more promising to me than the O(nm) time algorithm.

  Implementing algorithm 1 should not be too much work given what we already have I think.

• The 'SimplePolygon' type should keep the guarantee that the boundary is non-self intersecting. We can add a separate type that represents possibly selfintersecting polygons if need be.

• In the 0.X range, there used to be another 'MultiPolygon' type representing polygons with holes. It used to be defined as s.t. like

  haskell data MultiPolygon = MultiPolygon { holes :: [SimplePolygon] , outerBoundary :: SimplePolygon }

  but I found it a not so super helpful type. I think simply representing MultiPolygons using the (to be ported) PlanarSubdivision type might be more convenient/consistent.

Having said all of that. Having any polygon intersection algorithm is better than none, so if you want to implement and add Greiner-Hormann I would certainly welcome it :).

[mrehayden1]

There's also additionally the degenerate cases to address in the intersection detection with GH. I think the gist of it is that you would obtain one or more polygons with two duplicated points where the clip/subject polygon enter and exit the other at a single point.

The authors suggest to perturb the points where such cases exist. While it might work fine for rasterization, it doesn't sit right with me for my purposes.

I'm thinking de-duplication would be a more favourable alternative, but I haven't investigated it any further.

Maybe another route would be better?

[noinia]

The line segment intersection algorithm that is in hgeometry should be able handle such degeneracies. So I think producing the overlay itself should be fine.

[mrehayden1]

The degeneracies, as I understand, stem from computing the faces.

Are your proposing the Bentley–Ottmann algorithm for i)? There's some discussion on wikipedia about how to deal with the degenerate cases there.

Do you have any material I can read to understand the PlanarSubdivision choice? I'm rather new to computational geometry I'm afraid.

[mrehayden1]

I'm also lost when it comes to processing data of type of CoRec. While I understand the dependent typing, but I'm also struggling with tracking down the instances for intersections of LineSegments.

Edit: Nevermind, I found the instances. I think I'll have to read up on vinyl to fully understand how to use intersect, but I'm getting there.

[mrehayden1]

Okay, I think I understand the motivation for PlanarSubdivision now.

I was hoping to implement and start using an implementation pretty quickly with the existing API. I can see the hgeometry_0 branch is tagged with 0.15, and it looks as if master is undergoing some heavy refactoring.

Can I branch this off 0.15 for it to be ported to the work going on in master later?

[noinia]

The degeneracies, as I understand, stem from computing the faces.

Are your proposing the Bentley–Ottmann algorithm for i)? There's some discussion on wikipedia about how to deal with the degenerate cases there.

Indeed. This is the algorithm implemented for LineSegmentIntersection, and like I said, it already takes care of degeneracies as is. (And if it isn't its a bug ;)).

Do you have any material I can read to understand the PlanarSubdivision choice? I'm rather new to computational geometry I'm afraid.

No worries, I think your best bet may be to look at chapter 2 of the book "Computational Geometry - Algorithms and Applications by de Berg, Cheong, van Kreveld, and Overmars, third edition, 2008." You can also have a look at the slides that go with that chapter. For example, For the LineSegment intersection algorithm and for the overlay algorithm itself.

Okay, I think I understand the motivation for PlanarSubdivision now.

I was hoping to implement and start using an implementation pretty quickly with the existing API. I can see the hgeometry_0 branch is tagged with 0.15, and it looks as if master is undergoing some heavy refactoring.

Can I branch this off 0.15 for it to be ported to the work going on in master later?

I think your better of working from master directly. As you observed there are quite a few changes, so basing on master might make your life easier. (For example: many of the functions now can take more general arguments (controlled by typeclasses), and the intersections are now represented using regular data types rather than vinyl's CoRecs).