How to compile Heike's winding number function?

10

2

Heike gave the following function for winding number:

winding[poly_, pt_] :=  
Round[(Total@Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 
   2 Pi, -Pi]/2./Pi)]

I attempted to compile it as follows:

winding2 := Compile[{{poly, _Real, 2}, {pt, _Real, 1}},   
Round[(Total@Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 
    2 Pi, -Pi]/2./Pi)]]

Applied to the following simple problem, the compiled version gives error messages:

poly = {{0., 0.}, {10., 0.}, {10., 6.}, {0, 6}, {0., 0.}};
pt = {5., 3.}; 
winding2[poly, pt]

The error messages include:

Compile::cpapot: Compilation of ArcTan@@(ptCompile`GetElement[poly,System`Private`CompileSymbol[0]]) 
is not supported for the function argument ArcTan. The only function arguments supported are 
Times, Plus, or List. Evaluation will use the uncompiled function. >>
CompiledFunction::cfse: Compiled expression 6.283185307179586` should be a machine-size integer. >>
CompiledFunction::cfex: Could not complete external evaluation at instruction 1;
proceeding with uncompiled evaluation. >>

Where am I going wrong?

Duns

Posted 2013-07-14T09:27:29.667

Reputation: 275

1In version 11.3, Graphics`Mesh`PointWindingNumber seems to have been moved to Graphics`PolygonUtils`PointWindingNumber. – Henrik Schumacher – 2018-05-07T08:58:04.887

1You can also use the undocumented Graphics`Mesh`PointWindingNumber which seems to be just as fast as the compiled version. – rm -rf – 2013-07-15T01:35:17.313

Answers

17

First, if you use := in your assignment, then the compilation is not done instantly but every time you call winding2. That's btw the reason why you get the error message when you try to call the function because it is not compiled until then and the error is a compilation error.

Secondly, as the error messages sais, @@ can only be used with Times, Plus or List, so you have to expand this part:

winding2 = 
 Compile[{{poly, _Real, 2}, {pt, _Real, 1}}, 
  Round[Total@Mod[# - RotateRight[#] &@(ArcTan[#[[1]], #[[2]]] &@
    (Transpose@poly - pt)), 2 Pi, -Pi]/(2. Pi)]
 ]

Seems to work pretty smoothly

enter image description here

And here the code:

winding2 = 
  Compile[{{poly, _Real, 2}, {pt, _Real, 1}}, 
   0 != Round[
     Total@Mod[# - 
           RotateRight[#] &@(ArcTan[#[[1]], #[[2]]] &@(Transpose@
              poly - pt)), 2 Pi, -Pi]/(2. Pi)], 
   CompilationTarget -> "C", RuntimeOptions -> "Speed"];

With[{gr = 
   RegionPlot[x^2 + y^3 < 2, {x, -2, 2}, {y, -2, 2}, Mesh -> All, 
    FrameTicks -> None, PlotPoints -> 3, MaxRecursion -> 4, 
    PlotStyle -> RGBColor[14/15, 232/255, 71/85], 
    MeshStyle -> RGBColor[88/255, 22/51, 39/85]]}, 
 DynamicModule[{pt = {0, 1}, polyPts}, 
  polyPts = 
   gr[[1, 1, #]] & /@ 
    Flatten[Cases[gr, Polygon[pts__] :> Sequence[pts], Infinity], 1];
  LocatorPane[Dynamic[pt], 
   Dynamic@Show[gr, 
     Graphics[{Opacity[.5], RGBColor[38/255, 139/255, 14/17], 
       Polygon[Pick[polyPts, winding2[#, pt] & /@ polyPts]]}]]]]]

halirutan

Posted 2013-07-14T09:27:29.667

Reputation: 109 574

1Nice animation (+1). A bit faster: 0 != Round[Total@Mod[# - RotateRight[#] &@(ArcTan[#[[1]], #[[2]]] &@(Transpose@poly - pt)), 2 Pi, -Pi]/(2. Pi)]. You might want 0 != instead of 1 ==, in case the winding number is -1. The OP can worry about an even/odd rule or whether to use the winding number itself, instead of inside/outside the polygon. – Michael E2 – 2013-07-15T01:14:21.780

@MichaelE2 Thanks, I edited my post. Honestly, I haven't barely looked at the implementation of Heike. I just saw the bug and fixed this. Thanks, for looking in more detail at this. – halirutan – 2013-07-15T01:28:30.707