Skip to main content
code edit for better performance
Source Link
Simon Woods
  • 85.4k
  • 8
  • 180
  • 326

Edit##

##Edit##

As Mr Wizard discovered, the function above does not work unless the last point in the polygon is the same as the first. Here is a version which doesn't have that limitation, and as a bonus is slightly faster.

Edit 2 : Code tweaked for more speed (thanks again to Mr. Wizard)

inPolyQ2 = Compile[{{poly, _Real, 2}, {x, _Real}, {y, _Real}},
   Block[{Xi, Yi, Xip1, Yip1, u, v, w},
    {Xi, Yi} = Transpose@poly;
    {Xip1, = RotateLeft@Xi;
    Yip1} = Transpose@RotateLeft@poly;RotateLeft@Yi;
    u = UnitStep[y - Yi];
    v = UnitStep[y - Yip1];RotateLeft@u;
    w = UnitStep[-((Xip1 - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi))];
    Plus @@ Total[(u (1 - v) (1 - w) - (1 - u) v w)] != 0]];

Edit##

As Mr Wizard discovered, the function above does not work unless the last point in the polygon is the same as the first. Here is a version which doesn't have that limitation, and as a bonus is slightly faster.

inPolyQ2 = Compile[{{poly, _Real, 2}, {x, _Real}, {y, _Real}},
   Block[{Xi, Yi, Xip1, Yip1, u, v, w},
    {Xi, Yi} = Transpose@poly;
    {Xip1, Yip1} = Transpose@RotateLeft@poly;
    u = UnitStep[y - Yi];
    v = UnitStep[y - Yip1];
    w = UnitStep[-((Xip1 - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi))];
    Plus @@ (u (1 - v) (1 - w) - (1 - u) v w) != 0]];

##Edit##

As Mr Wizard discovered, the function above does not work unless the last point in the polygon is the same as the first. Here is a version which doesn't have that limitation, and as a bonus is slightly faster.

Edit 2 : Code tweaked for more speed (thanks again to Mr. Wizard)

inPolyQ2 = Compile[{{poly, _Real, 2}, {x, _Real}, {y, _Real}},
   Block[{Xi, Yi, Xip1, Yip1, u, v, w},
    {Xi, Yi} = Transpose@poly;
    Xip1 = RotateLeft@Xi;
    Yip1 = RotateLeft@Yi;
    u = UnitStep[y - Yi];
    v = RotateLeft@u;
    w = UnitStep[-((Xip1 - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi))];
    Total[(u (1 - v) (1 - w) - (1 - u) v w)] != 0]];
added 1016 characters in body
Source Link
Simon Woods
  • 85.4k
  • 8
  • 180
  • 326

The second "Neat Example" in the documentation for SmoothKernelDistribution contains this compiled function:

(* A region function for a bounding polygon using winding numbers: *)

inPolyQ = 
  Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}}, 
   Block[{polySides = Length[polygon], X = polygon[[All, 1]], 
     Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = 0, i = 1}, 
    While[i < polySides, Yi = Y[[i]]; Yip1 = Y[[i + 1]]; 
     If[Yi <= y, 
      If[Yip1 > y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) > 0, 
         wn++;];];, 
      If[Yip1 <= y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) < 0, 
         wn--;];];]; i++]; ! wn == 0]];

Edit##

As Mr Wizard discovered, the function above does not work unless the last point in the polygon is the same as the first. Here is a version which doesn't have that limitation, and as a bonus is slightly faster.

inPolyQ2 = Compile[{{poly, _Real, 2}, {x, _Real}, {y, _Real}},
   Block[{Xi, Yi, Xip1, Yip1, u, v, w},
    {Xi, Yi} = Transpose@poly;
    {Xip1, Yip1} = Transpose@RotateLeft@poly;
    u = UnitStep[y - Yi];
    v = UnitStep[y - Yip1];
    w = UnitStep[-((Xip1 - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi))];
    Plus @@ (u (1 - v) (1 - w) - (1 - u) v w) != 0]];

Comparison showing that the defect in the original is not present in the new code:

poly = Table[RandomReal[{7, 10}] {Sin[th], Cos[th]}, {th, 2 Pi/100, 2 Pi, 2 Pi/100}];

Grid[Timing[RegionPlot[#[poly, x, y], {x, -15, 15}, {y, -15, 15}, 
     PlotPoints -> 100]] & /@ {inPolyQ, inPolyQ2}]

enter image description here

The second "Neat Example" in the documentation for SmoothKernelDistribution contains this compiled function:

(* A region function for a bounding polygon using winding numbers: *)

inPolyQ = 
  Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}}, 
   Block[{polySides = Length[polygon], X = polygon[[All, 1]], 
     Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = 0, i = 1}, 
    While[i < polySides, Yi = Y[[i]]; Yip1 = Y[[i + 1]]; 
     If[Yi <= y, 
      If[Yip1 > y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) > 0, 
         wn++;];];, 
      If[Yip1 <= y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) < 0, 
         wn--;];];]; i++]; ! wn == 0]];

The second "Neat Example" in the documentation for SmoothKernelDistribution contains this compiled function:

(* A region function for a bounding polygon using winding numbers: *)

inPolyQ = 
  Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}}, 
   Block[{polySides = Length[polygon], X = polygon[[All, 1]], 
     Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = 0, i = 1}, 
    While[i < polySides, Yi = Y[[i]]; Yip1 = Y[[i + 1]]; 
     If[Yi <= y, 
      If[Yip1 > y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) > 0, 
         wn++;];];, 
      If[Yip1 <= y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) < 0, 
         wn--;];];]; i++]; ! wn == 0]];

Edit##

As Mr Wizard discovered, the function above does not work unless the last point in the polygon is the same as the first. Here is a version which doesn't have that limitation, and as a bonus is slightly faster.

inPolyQ2 = Compile[{{poly, _Real, 2}, {x, _Real}, {y, _Real}},
   Block[{Xi, Yi, Xip1, Yip1, u, v, w},
    {Xi, Yi} = Transpose@poly;
    {Xip1, Yip1} = Transpose@RotateLeft@poly;
    u = UnitStep[y - Yi];
    v = UnitStep[y - Yip1];
    w = UnitStep[-((Xip1 - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi))];
    Plus @@ (u (1 - v) (1 - w) - (1 - u) v w) != 0]];

Comparison showing that the defect in the original is not present in the new code:

poly = Table[RandomReal[{7, 10}] {Sin[th], Cos[th]}, {th, 2 Pi/100, 2 Pi, 2 Pi/100}];

Grid[Timing[RegionPlot[#[poly, x, y], {x, -15, 15}, {y, -15, 15}, 
     PlotPoints -> 100]] & /@ {inPolyQ, inPolyQ2}]

enter image description here

added 9 characters in body
Source Link

The second "Neat Example" in the documentation for SmoothKernelDistributionSmoothKernelDistribution contains this compiled function:

(* A region function for a bounding polygon using winding numbers: *)

inPolyQ = 
  Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}}, 
   Block[{polySides = Length[polygon], X = polygon[[All, 1]], 
     Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = 0, i = 1}, 
    While[i < polySides, Yi = Y[[i]]; Yip1 = Y[[i + 1]]; 
     If[Yi <= y, 
      If[Yip1 > y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) > 0, 
         wn++;];];, 
      If[Yip1 <= y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) < 0, 
         wn--;];];]; i++]; ! wn == 0]];

The second "Neat Example" in the documentation for SmoothKernelDistribution contains this compiled function:

(* A region function for a bounding polygon using winding numbers: *)

inPolyQ = 
  Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}}, 
   Block[{polySides = Length[polygon], X = polygon[[All, 1]], 
     Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = 0, i = 1}, 
    While[i < polySides, Yi = Y[[i]]; Yip1 = Y[[i + 1]]; 
     If[Yi <= y, 
      If[Yip1 > y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) > 0, 
         wn++;];];, 
      If[Yip1 <= y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) < 0, 
         wn--;];];]; i++]; ! wn == 0]];

The second "Neat Example" in the documentation for SmoothKernelDistribution contains this compiled function:

(* A region function for a bounding polygon using winding numbers: *)

inPolyQ = 
  Compile[{{polygon, _Real, 2}, {x, _Real}, {y, _Real}}, 
   Block[{polySides = Length[polygon], X = polygon[[All, 1]], 
     Y = polygon[[All, 2]], Xi, Yi, Yip1, wn = 0, i = 1}, 
    While[i < polySides, Yi = Y[[i]]; Yip1 = Y[[i + 1]]; 
     If[Yi <= y, 
      If[Yip1 > y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) > 0, 
         wn++;];];, 
      If[Yip1 <= y, Xi = X[[i]]; 
        If[(X[[i + 1]] - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi) < 0, 
         wn--;];];]; i++]; ! wn == 0]];
Source Link
Simon Woods
  • 85.4k
  • 8
  • 180
  • 326
Loading