I have 6894 polygons describing zones in a state plane (i.e., cartesian) coordinate system. Most are not very complex, and I believe none have holes. A random example is

Polygon[{{351633., 564236.}, {351681., 562612.}, {351663., 562555.}, {351699., 561341.}, {348792., 561268.}, {348733., 564168.}, {351633., 564236.}}]

  1. How can I test whether these polygons have holes?

  2. I want a new set of polygons that describe all areas in which one or more of these polygons intersect, and I need to know which polygons cover each of these regions. So, in the new set the region where polygon A is unintersected is one region, the region where A and B intersect is another region, the region where A, B, and C intersect is another, and so on. Any of these "regions" in the new set is not necessarily connected; for example, one could imagine a wide polygon being split in half by a thin one, producing two unintersected regions. In such cases I need a list of closed polygons (associated, as the other polygons, with the appropriate members of the domain set).

I have no idea how to do this mathematically, and it would be very cool to see a solution if someone has one. But on a more practical level I suspect that I can do this at high precision using image processing.

Here is an example of what is going on:

polygons = Get["https://s3.amazonaws.com/tblackburn/stack_exchange/polygons-sample"];
(* my polygons are benign, but beware what you "Get" from the internet! *)
n = 100;
Riffle[ColorData["Rainbow"] /@ Rescale[Range[n + 1]], RandomChoice[polygons, n]];
Graphics[Prepend[%, [email protected]], ImageSize -> Large]

enter image description here

Is there a way to scale a set of colors such that I can calculate what color each intersection would be, and such that the color is guaranteed to be unique? And then use the morphology suite of functions on a rasterized image to generate new polygons? I need the new polygons expressed in the same coordinate system, and I need the vertices to be as accurate as possible -- so the rasters would have to be quite large I think.

Or, is there some other way to obtain the result I am after?

I have 32 CPUs and 244 GB of memory to throw at this.

  • $\begingroup$ I love "such that I can calculate what color each intersection would be". Sounds like magic :) $\endgroup$ Nov 19, 2014 at 4:04
  • $\begingroup$ @belisarius I have seen image processing magic from you before :) $\endgroup$
    – mfvonh
    Nov 19, 2014 at 4:49
  • 3
    $\begingroup$ Oh! But I was sober that day! I almost remember ... $\endgroup$ Nov 19, 2014 at 4:58
  • $\begingroup$ You could use a binary encoding: assign each shape a colour that has exactly 1 nonzero bit, rasterize each on a black background, then use ImageAdd. But you can only distinguish up to 24 shapes. $\endgroup$
    – user484
    Nov 19, 2014 at 6:40
  • $\begingroup$ Another problem with an image processing approach is that if two edges intersect at a small angle you will get a chain of disconnected pixels: i.sstatic.net/TmEtp.png $\endgroup$
    – user484
    Nov 19, 2014 at 6:45

1 Answer 1


Well, since you bragged about computing power...

findIntersections[polyGroup__] := Module[{reg},
  reg = RegionIntersection[
    Sequence @@ 
       DiscretizeRegion[#, PerformanceGoal -> "Quality", 
        Method -> "DiscretizeGraphics", AccuracyGoal -> 10, 
        PrecisionGoal -> 10] & /@ polygons[[polyGroup]]];
  If[RegionDimension@reg > 0, 
    intPolygons, {Fold[regDiffPatch, reg, intPolygons[[All, 1]]], 

...you can just enumerate all the possible intersections and try them one by one, in reverse order:

intPolygons = {};
res = findIntersections /@ Most@Reverse@Subsets[Range[3]] // Quiet

Row[Column[{#, Graphics@polygons[[#]]}] & /@ {1, 2, 3}]
Column /@ intPolygons // Row
Row[{Show[intPolygons[[All, 1]]], Graphics@polygons[[{1, 2, 3}]]}]

Polygon list

The quality makes me cry and if I add one more polygon MMA just won't calculate the last intersection, and I don't even know why, because whatever the output is, it takes so long to appear that I never waited enough to see it.

But, hey, with that computing power you can just crank the Accuracy and Precision goals to over 9000 and you're probably good to go!

Here's the best a 2012 notebook could do:

Polygon list with 4

Here's the missing piece of code form the main function. It's just a bunch of silly heuristics to avoid some of the RegionDifference limitations (in many cases it won't compute and won't tell you why):

regDiffPatch["", reg2_] := "";
regDiffPatch[reg1_, ""] := reg1;
regDiffPatch[reg1_, reg2_] := 
 If[(Abs[Area@RegionUnion[reg1, reg2] - (Area@reg1 + Area@reg2)]/
     Area@RegionUnion[reg1, reg2]) < 0.01
  , reg1
  , Check[If[Area@RegionIntersection[reg1, reg2] == Area@reg1
    , ""
    , If[(Area@RegionIntersection[reg1, reg2]/Area@reg1) < 0.01, reg1,
      RegionDifference[reg1, reg2]]], RegionDifference[reg1, reg2]]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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