# How to decompose a self intersecting polygon into simple polygons?

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

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

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

• Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//GraphicsPolygonUtilsSimplePolygonPartition No response. May 7, 2018 at 8:30
• @GalAster This works: Block[{RegionMeshIntersectionFreeSegmentsQ},Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//GraphicsPolygonUtilsSimplePolygonPartition] Jul 27, 2018 at 10:47
• wow, that's cool, but I can't understand the Block here. Jul 27, 2018 at 11:11

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 = CompileGetElement[x, 1];
x2 = CompileGetElement[x, 2];

Do[
p1 = CompileGetElement[poly, i, 1];
p2 = CompileGetElement[poly, i, 2];
q1 = CompileGetElement[poly, i + 1, 1];
q2 = CompileGetElement[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] • The syntax is reversed, but you could also use the internal function GraphicsPolygonUtilsPointWindingNumber to find the winding number. May 6, 2018 at 16:16
• But mine is faster... ;) May 6, 2018 at 16:22
• On OSX with M11.3, I find that yours is slower? What is your configuration? May 6, 2018 at 16:40
• M11.3, macOS 10.13.4, i7 4980HQ @ 2,8 GHz, 16GB DDR3 RAM @1600 Mhz. Code runs 5 times faster here. This is not the first time that we get quite different timings. May 6, 2018 at 16:43
• Also related: Heike's algorithm here, which was written in compiled form by halirutan here. In the latter, rm -rf also points to GraphicsMeshPointWindingNumber. It may be the same as GraphicsPolygonUtilsPointWindingNumber. +1 May 6, 2018 at 16:50