10

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?

Community
  • 1
Duns
  • 275
  • 1
  • 5
  • 1
    You can also use the undocumented Graphics`Mesh`PointWindingNumber which seems to be just as fast as the compiled version. – rm -rf Jul 15 '13 at 01:35
  • 1
    In version 11.3, Graphics`Mesh`PointWindingNumber seems to have been moved to Graphics`PolygonUtils`PointWindingNumber. – Henrik Schumacher May 07 '18 at 08:58

1 Answers1

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
  • 112,764
  • 7
  • 263
  • 474
  • 1
    Nice 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 Jul 15 '13 at 01:14
  • @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 Jul 15 '13 at 01:28