How to decompose a self intersecting polygon into simple polygons?

5

1

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

enter image description here

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

GalAster

Posted 2018-05-06T13:47:20.267

Reputation: 3 440

Answers

6

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

Carl Woll

Posted 2018-05-06T13:47:20.267

Reputation: 112 778

Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics`PolygonUtils`SimplePolygonPartition No response. – GalAster – 2018-05-07T08:30:21.940

1@GalAster This works: Block[{Region`Mesh`IntersectionFreeSegmentsQ},Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics`PolygonUtils`SimplePolygonPartition] – chyanog – 2018-07-27T10:47:29.373

wow, that's cool, but I can't understand the Block here. – GalAster – 2018-07-27T11:11:38.830

4

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]

enter image description here

Henrik Schumacher

Posted 2018-05-06T13:47:20.267

Reputation: 85 430

1The syntax is reversed, but you could also use the internal function Graphics`PolygonUtils`PointWindingNumber to find the winding number. – Carl Woll – 2018-05-06T16:16:26.553

But mine is faster... ;) – Henrik Schumacher – 2018-05-06T16:22:09.607

On OSX with M11.3, I find that yours is slower? What is your configuration? – Carl Woll – 2018-05-06T16:40:37.517

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. – Henrik Schumacher – 2018-05-06T16:43:57.770

1

Also related: Heike's algorithm here, which was written in compiled form by halirutan here. In the latter, rm -rf also points to Graphics`Mesh`PointWindingNumber. It may be the same as Graphics`PolygonUtils`PointWindingNumber. +1

– C. E. – 2018-05-06T16:50:37.107

@C.E. In version 11.3, Graphics`Mesh`PointWindingNumber seems to be nonexistent. It has probably been moved to Graphics`PolygonUtils`PointWindingNumber. – Henrik Schumacher – 2018-05-07T09:00:00.323

@CarlWoll Did you try that with a long list of points in a listable fashion like Partition[ GraphicsPolygonUtilsPointWindingNumber[poly, Flatten[pts, 1]], Dimensions[pts][[2]]]? – Henrik Schumacher – 2018-05-07T09:04:30.250