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$

2 Answers 2

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$
3
  • $\begingroup$ Polygon[Table[{Cos[t],Sin[t]},{t,0,4 Pi,(4 Pi)/5}]]//Graphics`PolygonUtils`SimplePolygonPartition No response. $\endgroup$
    – Aster
    May 7, 2018 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, 2018 at 10:47
  • $\begingroup$ wow, that's cool, but I can't understand the Block here. $\endgroup$
    – Aster
    Jul 27, 2018 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$
7
  • 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, 2018 at 16:16
  • $\begingroup$ But mine is faster... ;) $\endgroup$ May 6, 2018 at 16:22
  • $\begingroup$ On OSX with M11.3, I find that yours is slower? What is your configuration? $\endgroup$
    – Carl Woll
    May 6, 2018 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$ May 6, 2018 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, 2018 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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