Obtain polygons describing all intersections of many polygons

Question

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]

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.

Answer 1

Well, since you bragged about computing power...

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

...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}]]}]

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:

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]]
  ]