5

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

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
Aster
  • 3,836
  • 1
  • 18
  • 44

2 Answers2

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
  • 130,679
  • 6
  • 243
  • 355
  • Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics`PolygonUtils`SimplePolygonPartition No response. – Aster May 07 '18 at 08:30
  • 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 Jul 27 '18 at 10:47
  • wow, that's cool, but I can't understand the Block here. – Aster Jul 27 '18 at 11:11
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
  • 106,770
  • 7
  • 179
  • 309
  • 1
    The syntax is reversed, but you could also use the internal function Graphics`PolygonUtils`PointWindingNumber to find the winding number. – Carl Woll May 06 '18 at 16:16
  • But mine is faster... ;) – Henrik Schumacher May 06 '18 at 16:22
  • On OSX with M11.3, I find that yours is slower? What is your configuration? – Carl Woll May 06 '18 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. – Henrik Schumacher May 06 '18 at 16:43
  • 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. May 06 '18 at 16:50
  • @C.E. In version 11.3, Graphics`Mesh`PointWindingNumber seems to be nonexistent. It has probably been moved to Graphics`PolygonUtils`PointWindingNumber. – Henrik Schumacher May 07 '18 at 09:00
  • @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 May 07 '18 at 09:04