6

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[%, Opacity@.5], 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.

mfvonh
  • 8,460
  • 27
  • 42
  • I love "such that I can calculate what color each intersection would be". Sounds like magic :) – Dr. belisarius Nov 19 '14 at 04:04
  • @belisarius I have seen image processing magic from you before :) – mfvonh Nov 19 '14 at 04:49
  • 3
    Oh! But I was sober that day! I almost remember ... – Dr. belisarius Nov 19 '14 at 04:58
  • 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. –  Nov 19 '14 at 06:40
  • 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: http://i.stack.imgur.com/TmEtp.png –  Nov 19 '14 at 06:45
  • Just a note that the colour you will get for overlaps depends on the ordering of the polygons. Try evaluating Row[Graphics /@ Map[{Opacity[0.5], #, Disk[]} & {{Red, Yellow, Blue}, {Yellow, Blue, Red}}, {2}]] and you'll see what I mean. A way out would be to GatherBy color and then Sorting and look at the colour spectrum you get and whether it matches the topological properties you are after. The colour issue has been observed in this question – gpap Nov 19 '14 at 12:11
  • Related: http://stackoverflow.com/questions/6805333/determine-if-polygon-has-a-hole – alancalvitti Nov 22 '14 at 02:54
  • polygons[[9]] includes each perimeter point twice, so it is all hole. – Simon Woods Nov 22 '14 at 11:43
  • If a boundary has no self intersection then it cannot have any holes. That would require multiple boundaries. – Daniel Lichtblau Nov 22 '14 at 21:07

1 Answers1

3

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

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]]
  ]
Aisamu
  • 2,618
  • 14
  • 17