31

This seems non-trivial (it gets asked quite a lot on various forums), but I absolutely need this as a building block for a more complex algorithm.

Input: 2 polygons (A and B) in 2D, given as a list of edges [(x0, y0, x1, y2), ...] each. The points are represented by pairs of doubles. I do not know if they are given clockwise, counter-clockwise or in any direction at all. I do know that they are not necessarily convex.

Output: 3 polygons representing A, B and the intersecting polygon AB. Either of which may be an empty (?) polygon, e.g. null.

Hint for optimization: These polygons represent room and floor boundaries. So the room boundary will normally fully intersect with the floor boundary, unless it belongs to another floor on the same plane (argh!).

I'm kind of hoping someone has already done this in c# and will let me use their strategy/code, as what I have found so far on this problem is rather daunting.

EDIT: So it seems I'm not entirely chicken for feiling faint at the prospect of doing this. I would like to restate the desired output here, as this is a special case and might make computation simpler:

Output: First polygon minus all the intersecting bits, intersection polygons (plural is ok). I'm not really interested in the second polygon, just its intersection with the first.

EDIT2: I am currently using the GPC (General Polygon Clipper) library that makes this really easy!

mvw
  • 5,075
  • 1
  • 28
  • 34
Daren Thomas
  • 67,947
  • 40
  • 154
  • 200
  • 2
    This is more complicated than you might think. How do you plan to represent the resulting shape? Bear in mind that the two inputs could (a) not intersect at all, (b) intersect partially, resulting in a convex or concave closed shape, (c) intersect wholly, resulting in a shape with TWO boundaries (i.e., donut) where you must find a way to represent the inside and outside of the shape. – Jon Seigel Oct 06 '09 at 15:47
  • 3
    Indeed, Jon is right. Your problem is not well-stated; the resulting set *might not be a polygon*, and therefore the output of the function ought not to be a polygon. What do you want to do in the case where the resulting shape is not connected? Imagine for example a C-shaped poly intersecting an I shaped poly, resulting in a colon. – Eric Lippert Oct 06 '09 at 17:51
  • thanks, will have to think hard about this and come up with a solution. Didn't think of it that way before... – Daren Thomas Oct 06 '09 at 18:04
  • And as DreamWalker pointed out, a non-convex polygon cannot be defined in the way you describe. Put four points into a square, and then a single point in the center: you can create four different Pacman-shaped polygons from this configuration. Your problem is not solvable as it stands. – StriplingWarrior Oct 06 '09 at 19:36
  • StriplingWarrior, are you aware that I specify the edges and not only the points? – Daren Thomas Oct 07 '09 at 06:53
  • 1
    Most algorithms to do what you describe rely on the winding rule to make it possible, so your first step should probably be to connect all the edges into an ordered set of points with known winding (clockwise is most common, but I've seen counterclockwise as well). Once you have an ordered set of points, you can use dot products and the right-hand rule to quickly (well in O(m*n)) find if any point from polygon A is inside polygon B. This is a necessary precondition to determine what kind of output geometry you can get. – Daniel Pryden Oct 08 '09 at 19:36
  • 1
    It may also be useful for you to read up on Point Set Theory. This page describes the JTS topology suite's implementation: http://docs.codehaus.org/display/GEOTDOC/Point+Set+Theory+and+the+DE-9IM+Matrix . JTS does what you want, but is written in Java. You might want to check out GEOS (a C++ port of JTS), or the NetTopologySuite that is mentioned below: http://stackoverflow.com/questions/1526352/how-to-intersect-two-polygons/1527105#1527105 – Daniel Pryden Oct 08 '09 at 19:40
  • This academic paper explains how to do this for both convex and nonconvext polygons: http://www.gvu.gatech.edu/~jarek/graphics/papers/04PolygonBooleansMargalit.pdf – Charles Bretana Oct 05 '11 at 00:23

9 Answers9

17

Arash Partow's FastGEO library contains implementations of many interesting algorithms in computational geometry. Polygon intersection is one of them. It's written in Pascal, but it's only implementing math so it's pretty readable. Note that you will certainly need to preprocess your edges a little, to get them into clockwise or counterclockwise order.

ETA: But really, the best way to do this is to not do this. Find another way to approach your problem that doesn't involve arbitrary polygon intersections.

David Seiler
  • 9,557
  • 2
  • 31
  • 39
  • 7
    If you do this, be extremely careful and don't change anything of the algorithm from the original code. If possible, get a Pascal-to-C compiler or compile it into a library you can use, rather than trying to translate the code yourself. – jprete Oct 06 '09 at 18:13
14

If you are programming in .NET Framework, you may want to take a look at SqlGeometry class available in .NET assemblies shipped as Microsoft SQL Server System CLR Types

The SqlGeometry class provides STIntersection method

SqlGeometry g1 = SqlGeometry.Parse("POLYGON ((...))");
SqlGeometry g2 = SqlGeometry.Parse("POLYGON ((...))");
SqlGeometry intersection = g1.STIntersection(g2);
mloskot
  • 37,086
  • 11
  • 109
  • 136
  • 3
    Is this answer incorrect, so it got downvote? If it is, please point where is a bug. – mloskot Feb 04 '10 at 15:11
  • 1
    i agree! what's wrong with this solution? I personally do all my spatial stuff in the DB ... but if u need to do it in the code, leverage the same code :) – Pure.Krome Oct 03 '10 at 23:54
12

What I think you should do

Do not attempt to do this yourself if you can possibly help it. Instead, use one of the many available polygon intersection algorithms that already exist.

I was strongly considering the following codebase on the strength of their demonstration code and the fact that they mentioned their handling of most/all of the weird cases. You would need to donate an amount (of you/your company's choice) if you use it commercially, but it's worth it to get a robust version of this kind of code.

http://www.cs.man.ac.uk/~toby/gpc/

What I actually did was to use a polygon-intersection algorithm that is part of the Java2D libraries. You can possibly find something similar in MS's own C# libraries to use.

There are other options out there as well; look for "polygon clipper" or "polygon clipping", since the same basic algorithms that handle polygon intersection also tend to be usable for the general clipping cases.

Once you actually have a polygon clipping library, you just need to subtract polygon B from polygon A to get your first piece of output, and intersect polygons A and B to get your second piece of output.

How to roll your own, for the hopelessly masochistic

When I was considering rolling my own, I found the Weiler-Atherton algorithm to have the most potential for general polygon-cutting. I used the following as a reference:

http://cs1.bradley.edu/public/jcm/weileratherton.html

http://en.wikipedia.org/wiki/Weiler-Atherton

The details, as they say, are too dense to include here, but I have no doubt that you'll be able to find references on Weiler-Atherton for years to come. Essentially, you split all the points into those that are entering the final polygon or exiting the final polygon, then you form a graph out of all the points, and then walk the graph in the appropriate directions in order to extract all the polygon pieces you want. By changing the way you define and treat the "entering" and "exiting" polygons, you can achieve several possible polygon intersections (AND, OR, XOR, etc.).

It's actually fairly implementable, but like with any computational geometry code, the devil is in the degeneracies.

yurez
  • 2,826
  • 1
  • 28
  • 22
jprete
  • 3,769
  • 21
  • 24
3

You may also want to have a look at the NetTopologySuite or even try importing it into Sql server 2008 & it's spatial tools.

oharab
  • 4,405
  • 3
  • 19
  • 15
3

Clipper is an open source freeware polygon clipping library (written in Delphi and C++) that does exactly what you're asking - http://sourceforge.net/projects/polyclipping/

In my testing, Clipper is both significantly faster and far less prone to error than GPC (see more detailed comparisons here - http://www.angusj.com/delphi/clipper.php#features). Also, while there's source code for both Delphi and C++, the Clipper library also includes a compiled DLL to make it very easy to use the clipping functions in other (Windows) languages too.

Angus Johnson
  • 4,565
  • 2
  • 26
  • 28
2

If you dare to take a look to other people's GPL C++ code, you can see how do they do it in Inkscape:

http://bazaar.launchpad.net/~inkscape.dev/inkscape/trunk/view/head:/src/2geom/shape.cpp (line 126)

mloskot
  • 37,086
  • 11
  • 109
  • 136
fortran
  • 74,053
  • 25
  • 135
  • 175
2

A polygon is fully described by an ordered list of points (P1, P2, ..., Pn). The edges are (P1 - P2), (P2 - P3), ..., (Pn - P1). If you have two polygons A and B which overlaps, there will be a point An (from the list on points describing polygon A) which lies within the area surrounded by polygon B or vice versa (a point of B lies in A). If no such point is found, then the polygons does not overlap. If such a point is found (i.e. Ai) check the adjacent points of the polygon A(i-1) and A(i+1). Repeat until you find a point outside the area or all points are checked (then the first polygon lies completly within the second polygon). If you found a point outside then you can calculate the crossing point. Find the corresponding edge of polygon B and follow it with resersed roles = now check if the points of polygon B lie within polygon A. This way you can build a list of points which describe the overlapping polygon. If needed you should check if the polygons are identical, (P1, P2, P3) === (P2, P3, P1).

This is just an idea and there maybe better ways. If you find a working and tested solution I would recommend that you use it...

narozed

narozed
  • 21
  • 1
2

Try to use GIS tools for that, such as ArcObjects, TopologySuite, GEOS, OGR, etc. I'm not sure if all of these distributions are availuable to .net, but they all do the same.

George Silva
  • 3,454
  • 10
  • 39
  • 64
  • 1
    FYI, OGR (from GDAL/OGR) uses GEOS (http://trac.osgeo.org/geos/) library, so there is no difference regarding implementation of computational geometry between these two. – mloskot Feb 03 '10 at 19:10
2

This academic paper explains how to do this.

Charles Bretana
  • 143,358
  • 22
  • 150
  • 216