Return to Answer

Below my Mathematica implementation of the point in polygon problem which appears to be roughly 5x faster than the inPolyQ[] algorithm posted above.

with the following results

inPoly2: 0.202

testpoint: 0.25

testpoint2: 0.016

inPolyQ: 0.015

InsidePolygonQ: 12.277

inPolyQ2: 0.032

with the following results

inPoly2: 8.237

testpoint: 11.295

testpoint2: 0.156

inPolyQ: 0.436

inPolyQ2: 0.078

My verdict: There is an astonishing 3 orders of magnitude difference in the performance of the different routines. InsidePolygonQ[] while mathematically elegant is very slow. It pays to use either the undocumented routine for point in polygon in Mathematica, in this case testpoint2[] (with the usual caveats), or the compiled routine inPolyQ2[] which both had excellent performance for both simple and complex test polygons.

Below is my Mathematica implementation of the point in polygon problem which appears to be roughly 5x faster than the inPolyQ[] algorithm posted above.

with the following results

inPoly2: 0.202
testpoint: 0.25
testpoint2: 0.016
inPolyQ: 0.015
InsidePolygonQ: 12.277
inPolyQ2: 0.032

with the following results

inPoly2: 8.237
testpoint: 11.295
testpoint2: 0.156
inPolyQ: 0.436
inPolyQ2: 0.078

My verdict: There is an astonishing 3 orders of magnitude difference in the performance of the different routines. InsidePolygonQ[] while mathematically elegant, is very slow. It pays to use either the undocumented routine for point in polygon in Mathematica, in this case testpoint2[] (with the usual caveats), or the compiled routine inPolyQ2[] which both had excellent performance for both simple and complex test polygons.

Update Following a suggestion from ruebenko I've now investigated the actual performance of all the different point-in-polygon routines for two specific cases.

Test No1: Simple triangle polyon and testing using 5000 random test points

poly = {{-1, 0}, {0, 1}, {1, 0}};
pts = Partition[RandomReal[{-1, 1}, 10000], 2];
npts = Length@pts;
Print["inPoly2: ", 
 Timing[Table[inPoly2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint: ", 
 Timing[Table[testpoint[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint2: ", 
 Timing[Table[testpoint2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["inPolyQ: ", 
 Timing[Table[inPolyQ[poly, pts[[i]]], {i, npts}];][[1]]]
Print["InsidePolygonQ: ", 
 Timing[Table[InsidePolygonQ[poly, pts[[i]]], {i, npts}];] [[1]]]
Print["inPolyQ2: ", 
 Timing[Table[
     inPolyQ2[poly, pts[[i, 1]], pts[[i, 2]]], {i, npts}];][[1]]]

with the following results

inPoly2: 0.202

testpoint: 0.25

testpoint2: 0.016

inPolyQ: 0.015

InsidePolygonQ: 12.277

inPolyQ2: 0.032

Test No2: Very complicated polygon. The main CountryData[] polygon for Canada has over 10 000 vertices and a fairly complex shape. I've focused on the fastest routines and excluded the InsidePolygonQ[] routine in this case and used 200 test points.

p = CountryData["Canada", "Polygon"][[1, 1]];
poly = {Rescale[p[[All, 1]], {Min@#, Max@#} &@p[[All, 1]], {-1, 1}],
    Rescale[p[[All, 2]], {Min@#, Max@#} &@p[[All, 2]], {-1, 1}]} // 
   Transpose;
pts = Partition[RandomReal[{-1, 1}, 400], 2];
npts = Length@pts;
Print["inPoly2: ", 
 Timing[Table[inPoly2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint: ", 
 Timing[Table[testpoint[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint2: ", 
 Timing[Table[testpoint2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["inPolyQ: ", 
 Timing[Table[inPolyQ[poly, pts[[i]]], {i, npts}];][[1]]]
Print["inPolyQ2: ", 
 Timing[Table[
 inPolyQ2[poly, pts[[i, 1]], pts[[i, 2]]], {i, npts}];][[1]]]

with the following results

inPoly2: 8.237

testpoint: 11.295

testpoint2: 0.156

inPolyQ: 0.436

inPolyQ2: 0.078

My verdict: There is an astonishing 3 orders of magnitude difference in the performance of the different routines. InsidePolygonQ[] while mathematically elegant is very slow. It pays to use either the undocumented routine for point in polygon in Mathematica, in this case testpoint2[] (with the usual caveats), or the compiled routine inPolyQ2[] which both had excellent performance for both simple and complex test polygons.

Sometimes speed is an issue if there are many polygons and or many points to check. There is an excellent reference on this issue under http://erich.realtimerendering.com/ptinpoly/ with the main conclusion that the angle summation algorithm should be avoided if speed is the objective.

Below my Mathematica implementation of the point in polygon problem which appears to be roughly 5x faster than the inPolyQ[] algorithm posted above.

Test case - use triangle

In[8]:= poly = {{-1, 0}, {0, 1}, {1, 0}};

My code implementation

In[82]:= inPoly2[poly_, pt_] := Module[{c, nvert},
   nvert = Length[poly];
   c = False;
   For[i = 1, i <= nvert, i++,
    If[i != 1, j = i - 1, j = nvert];
    If[(
      ((poly[[i, 2]] > pt[[2]]) != (poly[[j, 2]] > pt[[2]])) && (pt[[
      1]] < (poly[[j, 1]] - 
         poly[[i, 1]])*(pt[[2]] - poly[[i, 2]])/(poly[[j, 2]] - 
          poly[[i, 2]]) + poly[[i, 1]])), c = ! c];
    ];
   c
   ];

An the timing output testing on point {0,0.99}

Timing[t1 = Table[inPolyQ[poly, 0, 0.99], {10000}];]
Timing[t2 = Table[inPoly2[poly, 0, 0.99], {10000}];]

Out[115]= {0.062, Null}
Out[116]= {0.016, Null}

Sometimes speed is an issue if there are many polygons and or many points to check. There is an excellent reference on this issue under http://erich.realtimerendering.com/ptinpoly/ with the main conclusion that the angle summation algorithm should be avoided if speed is the objective.

Below my Mathematica implementation of the point in polygon problem which appears to be roughly 5x faster than the inPolyQ[] algorithm posted above.

Test case - use triangle

poly = {{-1, 0}, {0, 1}, {1, 0}};

My code implementation

inPoly2[poly_, pt_] := Module[{c, nvert,i,j},
   nvert = Length[poly];
   c = False;
   For[i = 1, i <= nvert, i++,
    If[i != 1, j = i - 1, j = nvert];
    If[(
      ((poly[[i, 2]] > pt[[2]]) != (poly[[j, 2]] > pt[[2]])) && (pt[[
      1]] < (poly[[j, 1]] - 
         poly[[i, 1]])*(pt[[2]] - poly[[i, 2]])/(poly[[j, 2]] - 
          poly[[i, 2]]) + poly[[i, 1]])), c = ! c];
    ];
   c
   ];

An the timing output testing on point {0,0.99}

Timing[t1 = Table[inPolyQ[poly, 0, 0.99], {10000}];]
Timing[t2 = Table[inPoly2[poly, 0, 0.99], {10000}];]

Out[115]= {0.062, Null}
Out[116]= {0.016, Null}