How to decompose a self intersecting polygon into simple polygons?

Question

I want to divide a self intersecting polygon into simple polygons with winding number.

For example poly = Polygon[{{0, 1}, {1, -1}, {-3/10, 0}, {3/10, 0}, {-1, -1}}]

It should be divided into two parts: the gray part(winding=1) and the darker part(winding=2).

So I want to get a list like this:

{
    {2,Polygon[{{-3/10,0},{0,-3/13},{3/10,0}}]},
    {1,Polygon[{{3/10,0},{0,-3/13},{1,-1},{0,1},{-1,-1},{0,-3/13},{-3/10,0}}]}
}

In other words, I would like a function polyDecompose which takes as input a Polygon and returns a list of {windingnumber,Polygon[...]}.

Another example, a pentagram can be divided into 5 parts with winding number 1 and 1 part with winding number 2.

Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics

Answer 1

If you don't mind using undocumented, internal functions, then the following returns the desired polygons without the winding numbers:

Graphics`PolygonUtils`SimplePolygonPartition[
    Polygon@{{0,1},{1,-1},{-3/10,0},{3/10,0},{-1,-1}}
]

{Polygon[{{-0.3, 0.}, {2.77556*10^-17, -0.230769}, {0.3, 0.}}], Polygon[{{0.3, 0.}, {2.77556*10^-17, -0.230769}, {1., -1.}, {0., 1.}, {-1., -1.}, {2.77556*10^-17, -0.230769}, {-0.3, 0.}}]}

Answer 2

Not a complete answer, but it might be helpful to have a function that computes the winding number of a point x in the plane with respect to the closed polygon poly (represented by a sequence of points with first and last point coinciding).

windingnumber = Compile[{{x, _Real, 1}, {poly, _Real, 2}},
   Block[{wn, p1, p2, q1, q2, x1, x2},
    wn = 0;
    x1 = Compile`GetElement[x, 1];
    x2 = Compile`GetElement[x, 2];

    Do[
     p1 = Compile`GetElement[poly, i, 1];
     p2 = Compile`GetElement[poly, i, 2];
     q1 = Compile`GetElement[poly, i + 1, 1];
     q2 = Compile`GetElement[poly, i + 1, 2];
     If[p2 <= x2,
      If[q2 > x2,(* possible upcrossing *)

       If[(q1 - p1) (x2 - p2) - (q2 - p2) (x1 - p1) > 0, wn++]
       ],
      If[q2 <= x2,(* possible downcrossing *)

       If[(q1 - p1) (x2 - p2) - (q2 - p2) (x1 - p1) < 0, wn--]]
      ],
     {i, 1, Length[poly] - 1}];
    wn
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Here a simple usage example:

γ = t \[Function] {Cos[t], Sin[t]} + 0.5 {Cos[-6 t], Sin[-6 t]};
poly = γ /@ Subdivide[-Pi, Pi, 200];
pts = Outer[List, Subdivide[-1.5, 1.5, 600],Subdivide[-1.5, 1.5, 600]];
ArrayPlot@Rescale@windingnumber[pts, poly]