# Proper Filling of Self-Intersecting Polygon

I have a spiral pattern that I would like to fill, which self-intersects. My desired behavior is to have filling between neighboring lines only, and to represent this figure as a single polygon. The default behavior is somewhat different, producing an alternating color pattern:

The outer parts of the spiral are filling properly (only between the neighboring line), but the inner parts are not. There is probably a particular ordering of points that can be passed to Polygon that would produce the desired behavior, but I am not sure how to accomplish that.

The data is here, here's the code I use to currently generate the figure:

bhpl = Graphics[{Black, Circle[{0, 0}, 3.22823*10^11]}];
streampl = Graphics[{Red, Polygon[polydat]}];
Show[streampl, bhpl,
PlotRange -> {{-9.68468*10^11, 9.68468*10^11}, {-6.45645*10^11, 3.22823*10^12}}

-

The problem you face here is, that you draw one big polygon by going around your spiral with the coordinates. The solution is as simple as it sounds: Don't do that.

So how can you not go around? Let's start by looking from both ends on your list. Here I import and rescale your data and then I draw the first 3 points from the start in blue and the first 3 points from the end in red

Import["http://pastebin.com/raw.php?i=0ux8ep4W", "Package"];
data = Rescale[polydat];

Graphics[{PointSize[.02], Blue, Point[data[[1 ;; 3]]], Red,
Point[Take[data, -3]]}, PlotRange -> {{0.1, .14}, {0.05, .1}},
Frame -> True]


Now one possible solution should jump directly into your eye. Instead of drawing one big polygon around, we draw many which look like rectangles. The advantage of this is that when in a later stage a polygon is drawn over some other it just works because they are not the same. It's like drawing two things over another. Let's try this

getCoord[i_] :=
With[{p1 = Take[data, {i, i + 1}], p2 = Take[data, {-i - 1, -i}]},
Flatten[{p1, p2}, 1]]

Manipulate[
Graphics[{Red, Opacity[.2], EdgeForm[Black],
Polygon[getCoord /@ Range[i]]},
PlotRange -> {{0, .21}, {0, .21}}],
{i, 1, 200, 1}]


Looks like we nailed it. The work is done inside getCoord which takes an index i and extracts coordinates from your data and orders them in a way you can directly use it with Polygon.

Final solution is then

Graphics[{Red, EdgeForm[Red],
Polygon[getCoord /@ Range[Length[data]/2-1]]}]


et voila

## Building a GraphicsComplex

Building a GraphicsComplex where the coordinates are only used one time and in the Polygon primitive you work with indices is not hard. Just create a table of integers the same way you take coordinates. Two from the front and two from the back and you are done..

With[{l = Length[data]},
Graphics[{Red, EdgeForm[Red], GraphicsComplex[data,
Polygon[
Table[{i, i + 1, Mod[-i - 1, l] + 1, Mod[-i, l] + 1},{i,l/2 - 1}]]]}]
]

-
OK, that's what I suspected the solution would be like. However, I would like to "merge" the shapes together such that they are one polygon in the end. This is mostly because when the individual polygons are exported to PDF, the gaps between the shapes are noticeable, and it makes things look a bit ugly. How would I go about doing that? –  Guillochon Oct 4 '12 at 1:00
I imagine the solution to what I ask for above is going to be similar to the solution presented here: mathematica.stackexchange.com/questions/644/… –  Guillochon Oct 4 '12 at 1:02
Have you tried to just use my solution as it is and export it to pdf? ;-) –  halirutan Oct 4 '12 at 1:20
Thinking about it, this might be a good time to use GraphicsComplex[], so you don't have to keep repeating coordinates for each Polygon[]... –  Ｊ. Ｍ. Oct 4 '12 at 2:35
I just tried now, it looks OK in PDF. I guess my primary concern is the number of points, many of which are unnecessary as some of the polygons are completely hidden from view. –  Guillochon Oct 4 '12 at 3:53

Following on from Halirutan's answer, you can reduce the number of polygons using the undocumented functionality in GraphicsMesh

The code below uses Halirutan's getCoord function to obtain the polygons. The last polygon is dropped as it causes problems by being infinitely thin.

The general idea is to partition the polygon list into pairs and apply PolygonCombine to each pair, provided the pair don't intersect at more than two points. This means that adjacent polygons will be merged into one, but self-intersecting polygons won't be created. The process is repeated until the number of polygons no longer changes.

The final result contains just 9 polygons:

I suspect it's possible to go all the way down to a single polygon using the right functions, but it's a slow process to figure out all the undocumented functionality.

Here is the code:

data = Rescale[polydat];

getCoord[i_] :=
With[{p1 = Take[data, {i, i + 1}], p2 = Take[data, {-i - 1, -i}]},
Flatten[{p1, p2}, 1]];

p = Polygon /@ getCoord /@ Range[Length[data]/2];

p = Most@p;

GraphicsMeshMeshInit[];

pcomb[polys_] :=
If[Length[FindIntersections[polys]] > 2, polys, PolygonCombine[polys]]

p2 = FixedPoint[Flatten[pcomb /@ Partition[#, 2, 2, {1, 1}, {}]] &, p,
SameTest -> (Length[#1] == Length[#2] &)];

Graphics[{Yellow, EdgeForm[Black], p2}]

-
This is very nice. Too bad I can't accept both answers! –  Guillochon Oct 4 '12 at 18:26
This is indeed nice. The last polygon is infinitely thin because when I use i+1 in getCoord you only have to iterate until Length[data]/2-1`. I fixed this in my first edit. +1. –  halirutan Oct 5 '12 at 1:15