# How to compile Heike's winding number function?

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@@(ptCompileGetElement[poly,SystemPrivateCompileSymbol[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?

-
You can also use the undocumented GraphicsMeshPointWindingNumber which seems to be just as fast as the compiled version. –  The Toad Jul 15 '13 at 1:35

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

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

-
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 1: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 1:28