Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to quickly calculate the intersection of polygons with more than 6,000 points. A compiled solution would be preferable.

Here is one example of the problem:

o = First[
       First[ImportString[
         ExportString[
          Style["O", Italic, FontSize -> 24, FontFamily -> "Times"], 
          "PDF"], "PDF", "TextMode" -> "Outlines"]]];

p = First[
   First[ImportString[
     ExportString[
      Style["P", Italic, FontSize -> 24, FontFamily -> "Times"], 
      "PDF"], "PDF", "TextMode" -> "Outlines"]]];

Graphics[{EdgeForm[Black], ColorData["Crayola", "Sunglow"], {o, p}}]

Another:

    Module[{a = FilledCurve[{{Line[{{2, 3}, {0.8125, 0.625}}],
      BezierCurve[{{0.6875, 0.375}, {0.375, 0.25}, {1.125, 0.25}}, 
       SplineDegree -> 2], 
      BezierCurve[{{0.8125, 0.375}, {0.9375, 0.625}}],
      Line[{{1.3125, 1.375}, {2.4375, 1.375}, {2.8125, 0.625}}],
      BezierCurve[{{2.9375, 0.375}, {2.625, 0.25}, {3.625, 0.25}}, 
       SplineDegree -> 2], 
      BezierCurve[{{3.3125, 0.375}, {3.1875, 0.625}}]},
     {Line[{{1.875, 2.5}, {1.375, 1.5}, {2.375, 1.5}}]}}]}, 
 Graphics[Table[{EdgeForm[Black], Hue[RandomReal[]], 
    Translate[Rotate[Scale[a, RandomReal[5]], RandomReal[2 Pi]], 
     RandomReal[20, {2}]]}, {30}]]]

So in the problem is: how to intersect of multiple (or two at a time) polygons (with or without) holes (and with up to 6000 points in their triangulations)?

I have tried using the Weiler–Atherton clipping algorithm but my implementation was too slow (anything relying on bitmaps is too slow). Perhaps there is a solution that uses LibraryLink to harness a standard library? I found one here http://www.cs.man.ac.uk/~toby/alan/software/gpc.html

Updates

It was suggested in the comments that GraphicsMesh be used but this is way too slow in even 100 points and doesn't handle holes:

a = Polygon@RandomReal[1, {100, 2}];
b = Polygon@RandomReal[1, {100, 2}];
AbsoluteTiming[c = Graphics`Mesh`PolygonIntersection[a, b]]
Graphics[{Blue, a, Red, b, Yellow, Polygon /@ List @@ c}]

enter image description here

Here is a sample input polygon for the letter G:

 G = Polygon[{{-0.466796, -0.0328696}, {-0.466336, 0.0186753}, {-0.463089, 
   0.0682893}, {-0.459379, 0.100181}, {-0.451495, 
   0.146241}, {-0.444693, 0.175763}, {-0.432171, 0.21827}, {-0.422278,
    0.245423}, {-0.411147, 0.271628}, {-0.39878, 0.296886}, {-0.37791,
    0.332996}, {-0.362451, 0.355885}, {-0.345756, 
   0.377826}, {-0.327823, 0.398819}, {-0.300607, 
   0.426707}, {-0.281598, 0.443668}, {-0.251786, 
   0.466663}, {-0.231046, 0.480362}, {-0.198638, 
   0.498465}, {-0.176168, 0.508902}, {-0.141165, 
   0.522112}, {-0.116964, 0.529288}, {-0.0793649, 
   0.537605}, {-0.0534337, 0.541519}, {-0.013239, 
   0.544944}, {0.0144228, 0.545596}, {0.0521431, 
   0.544603}, {0.0973423, 0.540566}, {0.140374, 0.533423}, {0.181237, 
   0.523175}, {0.212366, 0.512741}, {0.249327, 0.496903}, {0.283288, 
   0.47845}, {0.310993, 0.459561}, {0.336333, 0.438254}, {0.359309, 
   0.414528}, {0.379921, 0.388384}, {0.398168, 0.359821}, {0.414052, 
   0.328839}, {0.427571, 0.295438}, {0.438725, 0.259618}, {0.444849, 
   0.234395}, {0.449921, 0.208096}, {0.309296, 0.208096}, {0.299494, 
   0.244597}, {0.29219, 0.264834}, {0.283819, 0.283824}, {0.267496, 
   0.312703}, {0.256279, 0.328367}, {0.243996, 0.342784}, {0.235215, 
   0.351703}, {0.221153, 0.364041}, {0.206024, 0.375133}, {0.189953, 
   0.385077}, {0.161346, 0.399325}, {0.14309, 0.406478}, {0.124015, 
   0.412584}, {0.0904026, 0.420435}, {0.0691431, 0.42375}, {0.0470644,
    0.426018}, {0.00844619, 0.427471}, {-0.0284694, 
   0.425533}, {-0.0551123, 0.421535}, {-0.0892443, 
   0.412813}, {-0.121786, 0.400214}, {-0.152737, 
   0.383739}, {-0.182097, 0.363387}, {-0.209866, 
   0.339158}, {-0.235987, 0.311056}, {-0.248019, 
   0.295569}, {-0.259191, 0.279133}, {-0.269503, 
   0.261747}, {-0.278956, 0.243411}, {-0.291525, 
   0.214127}, {-0.298829, 0.193418}, {-0.305275, 
   0.171759}, {-0.313331, 0.137489}, {-0.321066, 
   0.0884734}, {-0.325362, 
   0.0356589}, {-0.32624, -0.0188127}, {-0.324105, -0.0667413}, \
{-0.319123, -0.112199}, {-0.311294, -0.155187}, {-0.300618, \
-0.195705}, {-0.294213, -0.215037}, {-0.28327, -0.242878}, \
{-0.270726, -0.269329}, {-0.25151, -0.302436}, {-0.235136, \
-0.325433}, {-0.216766, -0.346138}, {-0.196376, -0.3645}, {-0.173967, \
-0.380518}, {-0.149538, -0.394191}, {-0.123089, -0.405521}, \
{-0.104334, -0.411771}, {-0.0846818, -0.41698}, {-0.0535202, \
-0.422841}, {-0.0203389, -0.426357}, {0.0148622, -0.427529}, \
{0.0509528, -0.426065}, {0.0853497, -0.421672}, {0.118053, \
-0.414352}, {0.149062, -0.404104}, {0.171208, -0.394496}, {0.1924, \
-0.383241}, {0.212641, -0.370339}, {0.231928, -0.35579}, {0.250217, \
-0.339572}, {0.266916, -0.321406}, {0.281843, -0.301207}, {0.290811, \
-0.286611}, {0.302785, -0.263022}, {0.312987, -0.237399}, {0.321418, \
-0.209743}, {0.326054, -0.190175}, {0.331531, -0.159129}, {0.334198, \
-0.137302}, {0.336722, -0.102866}, {0.337421, -0.0787786}, {0.174296, \
-0.0787786}, {0.0111708, -0.0787786}, {0.0111708, 
   0.0393464}, {0.238983, 0.0393464}, {0.466796, 
   0.0393464}, {0.466796, -0.517529}, {0.377235, -0.517529}, \
{0.343661, -0.386923}, {0.320697, -0.411362}, {0.294314, -0.43722}, \
{0.26921, -0.459296}, {0.245385, -0.477589}, {0.222841, -0.492099}, \
{0.200503, -0.503935}, {0.174026, -0.515381}, {0.145794, -0.524995}, \
{0.108033, -0.534437}, {0.0758485, -0.539931}, {0.0331483, \
-0.544223}, {-0.0029876, -0.545596}, {-0.0362581, -0.545124}, \
{-0.0713561, -0.542341}, {-0.10544, -0.537172}, {-0.138509, \
-0.529619}, {-0.170564, -0.51968}, {-0.191371, -0.511729}, \
{-0.211727, -0.502717}, {-0.231632, -0.492646}, {-0.260644, \
-0.475551}, {-0.279422, -0.462829}, {-0.29775, -0.449047}, \
{-0.315626, -0.434205}, {-0.324396, -0.426386}, {-0.34277, \
-0.405934}, {-0.367955, -0.373398}, {-0.390289, -0.338633}, \
{-0.403595, -0.314218}, {-0.415633, -0.288813}, {-0.426404, \
-0.262417}, {-0.440185, -0.220965}, {-0.447788, -0.192092}, \
{-0.456817, -0.146926}, {-0.461252, -0.115577}, {-0.465529, \
-0.0666954}, {-0.466796, -0.0328696}}];
share|improve this question
    
Great question! I've often wondered this myself... –  M.R. Sep 27 '12 at 1:08
    
mathematica.stackexchange.com/questions/7174/… only solved unions –  M.R. Sep 27 '12 at 1:09
    
The library you point to appears to be an outstanding option, particularly if you think you'll need to solve this problem many times and quickly each time. As a library, one could probably link to it relatively easily via a MathLink program. Also, the parent page points to a java port, which could be linked to via JLink. Both of those options would take a bit of time, though. (By parent page, I mean here: web.archive.org/web/20090213122910/http://www.seisw.com/GPCJ/… –  Mark McClure Sep 27 '12 at 1:56
    
1  
Theres ``GraphicsMesh`PolygonIntersection`` to check out perhaps? –  Rojo Sep 27 '12 at 4:02

2 Answers 2

Graphics`Mesh`PolygonIntersection[] is not documented; it builds full polygon triangulations. To handle holes, you can use:

PolygonIntersection[a, b, FillingMethod -> "OddEvenRule"]

or

PolygonIntersection[a, b, FillingMethod -> "WindingRule"]

To create the visualization:

Graphics`Mesh`MeshInit[];
a = Polygon@RandomReal[1, {100, 2}];
b = Polygon@RandomReal[1, {100, 2}];
c = Graphics`Mesh`PolygonIntersection[a, b, FillingMethod -> "OddEvenRule"];
Graphics[{Blue, Opacity[.5], a, Red, Opacity[.5], b, EdgeForm[Black], FaceForm[Yellow], 
          Polygon /@ List @@ c}]

intersections of two self-intersecting polygons

share|improve this answer
3  
@rm-rf Alternatively, executing Graphics`Mesh`MeshInit[] loads the full context, including any dependencies. –  Mark McClure Sep 27 '12 at 19:16
    
Hey Ulises! Welcome to the party! Got any more hidden gems for us? –  Yves Klett Sep 27 '12 at 20:55
1  
Right, PolygonIntersection[] is defined in the Graphics`Mesh` context. Just wanted to show which option needed to be set to make it more useful. As pointed out, it is not as fast as we want it to be, that is why is undocumented. Hi Yves!, if you are in town for the technology conference next month, we can discuss hidden and not so hidden features ... –  Ulises Cervantes Sep 27 '12 at 21:56
    
@UlisesCervantes If I try Graphics`Mesh`PolygonIntersection[a, b, FillingMethod -> "OddEvenRule"] I get no result where a and b is defined by Polygon@RandomReal[1, {100, 2}]!! –  PlatoManiac Sep 28 '12 at 8:48
    
In Mathematica V8.0.4. (or even 8.0): GraphicsMeshMeshInit[]; a = Polygon@RandomReal[1, {100, 2}]; b = Polygon@RandomReal[1, {100, 2}]; c = PolygonIntersection[a, b, FillingMethod -> "OddEvenRule"]; Graphics[{Blue, a, Red, b, Yellow, Polygon /@ List @@ c}] !Valid XHTML. –  Ulises Cervantes Sep 28 '12 at 18:09

I worked through this problem some years ago (Sorry I don't have mathematica code to share.) It turns out the most efficient algorithm is the fairly obvious direct approach. (assuming you want a 'precise' result..so rastering wont do)

1) compute all possible edge-edge intersections

2) decompose original polygons into edge sets on those boundaries

3) order edge sets by angle at each intersection point

4) construct sub-regions by following smallest angle at each intersection point

finally depending on what you want as output do an insideness test on each sub-region against the original input polygons.

The approaches discussed so far (meshing, etc) address the last parts, while for the problem posed with thousands of edges step 1 is by far the most expensive part of the calculation. Meshing actually makes that worse as now you have a least 2-3 x as many edges to deal with.

Now for the original example (Not the random table thing) we have thousands of mostly small segments. For this case you can dramatically reduce the O(N^2) edge-intersection problem by binning edges onto a regular grid and only check those edges that are near enough each other to possibly intersect.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.