5
$\begingroup$

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

$\endgroup$
6
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics`PolygonUtils`SimplePolygonPartition No response. $\endgroup$ – GalAster May 7 '18 at 8:30
  • 1
    $\begingroup$ @GalAster This works: Block[{Region`Mesh`IntersectionFreeSegmentsQ},Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics`PolygonUtils`SimplePolygonPartition] $\endgroup$ – chyanog Jul 27 '18 at 10:47
  • $\begingroup$ wow, that's cool, but I can't understand the Block here. $\endgroup$ – GalAster Jul 27 '18 at 11:11
4
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ The syntax is reversed, but you could also use the internal function Graphics`PolygonUtils`PointWindingNumber to find the winding number. $\endgroup$ – Carl Woll May 6 '18 at 16:16
  • $\begingroup$ But mine is faster... ;) $\endgroup$ – Henrik Schumacher May 6 '18 at 16:22
  • $\begingroup$ On OSX with M11.3, I find that yours is slower? What is your configuration? $\endgroup$ – Carl Woll May 6 '18 at 16:40
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher May 6 '18 at 16:43
  • 1
    $\begingroup$ 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 $\endgroup$ – C. E. May 6 '18 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.