Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
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 1:35
add comment

1 Answer

up vote 15 down vote accepted

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]]}]]]]]
share|improve this answer
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 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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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