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?
Answer
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]]}]]]]]
Comments
Post a Comment