# How to check if a 2D point is in a polygon?

Background: I use code from An Efficient Test For A Point To Be In A Convex Polygon Wolfram Demonstration to check if a point ( mouse pointer ) is in a ( convex ) polygon. Clearly this code fails for non-convex polygons.

Question: I am looking for an efficient routine to check if a 2D-point is in a polygon.

• I don't believe Mathematica has a build-in function for this. You could build your own in which case this is a good starting point: point in polygon Aug 14, 2012 at 8:39
• Perhaps this answer by Heike?
– kglr
Aug 14, 2012 at 8:49
• It's pretty depressing that you already accepted an answer before I posted mine, but I posted anyway. :-/ Aug 14, 2012 at 10:02
• The best answer to this question depends on the intended use and other constraints. The most important determinants are (1) whether this will be a one-off test or if many points will be tested for a given polygon; (2) where the points are likely to fall; and (3) whether the test needs to be absolutely accurate. Except for the one-off accurate test, by far the fastest method--and one not yet offered in any answer--is to rasterize the polygon's interior, probe the raster at the point's location (a $O(1)$ operation), and revert to a more expensive test only if the probe is inconclusive. Aug 14, 2012 at 17:52
• The various responses are quite good. That said, if you have to test many points and the polygon has many vertices, the method indicated at this link should be fairly efficient. forums.wolfram.com/mathgroup/archive/2009/Feb/msg00519.html Aug 14, 2012 at 21:10

Using the function winding from Heike's answer to a related question

 winding[poly_, pt_] :=
Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly),
2 Pi, -Pi]/2/Pi)]


to modify the test function in this Wolfram Demonstration by R. Nowak to

testpoint[poly_, pt_] :=
Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly),
2 Pi, -Pi]/2/Pi)] != 0


gives Update: Full code:

Manipulate[With[{p = Rest@pts, pt = First@pts},
Graphics[{If[testpoint[p, pt], Pink, Orange], Polygon@p},
PlotRange -> 3 {{-1, 1}, {-1, 1}},
ImageSize -> {400, 475},
PlotLabel -> Text[Style[If[testpoint[p, pt], "True ", "False"], Bold, Italic]]]],
{{pts, {{0, 0}, {-2, -2}, {2, -2}, {0, 2}}},
Sequence @@ (3 {{-1, -1}, {1, 1}}), Locator, LocatorAutoCreate -> {4, Infinity}},
SaveDefinitions -> True,
Initialization :> {
(* test if point pt inside polygon poly *)
testpoint[poly_, pt_] :=
Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly),
2 Pi, -Pi]/2/Pi)] != 0 } ]


Update 2: An alternative point-in-polygon test using yet another undocumented function:

 testpoint2[poly_, pt_] := GraphicsMeshInPolygonQ[poly, pt]

testpoint2[{{-1, 0}, {0, 1}, {1, 0}}, {1/3, 1/3}]
(*True*)
testpoint2[{{-1, 0}, {0, 1}, {1, 0}}, {1, 1}]
(*False*)

• Looks great! I'll try to implement this in my app and let you know asap. Thank you very much. Aug 14, 2012 at 9:14
• The testpoint function is a beauty. I would give a bonus if possible, I'll do a fast accept, instead. Thanks kguler Aug 14, 2012 at 9:57
• This doesn't work for self-intersecting polygons, but you could use PolygonTessellation from the package linked in my answer to pre-process. +1 Aug 14, 2012 at 10:39
• @matheorem, i used ??**PointWindingNumber*
– kglr
Dec 3, 2014 at 2:02
• @matheorem There is no GraphicsMesh* context in v10. You must be seeing it because a symbol was created in that context when you tried kguler's answer or mine (which didn't work). Try again in a fresh session and you won't see it.
– rm -rf
Dec 3, 2014 at 14:07

The undocumented GraphicsPolygonUtilsPointWindingNumber (if you're on versions < 10, use GraphicsMeshPointWindingNumber) does exactly this — it gives you the winding number of a point. A point lies inside the polygon if and only if its winding number is non-zero.

Using this, you can create a Boolean function to test if a point is inside the polygon

inPolyQ[poly_, pt_] := GraphicsPolygonUtilsPointWindingNumber[poly, pt] =!= 0

(* Examples *)
inPolyQ[{{-1, 0}, {0, 1}, {1, 0}}, {1/3, 1/3}]
(* True *)
inPolyQ[{{-1, 0}, {0, 1}, {1, 0}}, {1, 1}]
(* False *)


Or, you can use the aptly named GraphicsPolygonUtilsInPolygonQ which has the same 2-argument syntax and is a predicate.

• One could also use GraphicsMeshInPolygonQ, which has the same 2 argument syntax as GraphicsMeshPointWindingNumber. The former also takes a Method option, which I haven't explored yet...
– rm -rf
Aug 14, 2012 at 16:37
• I implement windnig number myself but this is better, I bet faster. A newbie Q: how can I see how this functions is coded, its syntax? I can see there's also PolygonWindingNumber in ver. 9 and I'd like to see what it does if I can.
– BoLe
Apr 19, 2013 at 8:16
• @BoLe Unfortunately, it doesn't look like the code for this one is accessible... In general, the tighter a function's integration with the kernel, the harder it is (almost impossible) to read its definitions (e.g. Pick, Cases, etc.). You can, however, read them for several other functions that are implemented in top level Mathematica. Usually all that one needs to do is to clear the attributes, read the definitions and follow the rabbit hole of internal definitions. However, instead of that, I'll recommend this answer. Prepare to be amazed :)
– rm -rf
Apr 19, 2013 at 13:41
• @rm-rf As I tested it, this undocumented GraphicsMeshPointWindingNumber is not working in V10, so maybe undocumented function is dangerous? Dec 3, 2014 at 0:43
• @rm-rf Thank you! kguler told me using ??**PointWindingNumber* could see which package that PointWindingNumber belongs to. But I found V10 gives both Mesh and PolygonUtils, you can see here imgur.com/bCIANiR, then why Mesh didn't work any more? Dec 3, 2014 at 6:35

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 is 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[]) != (poly[[j, 2]] > pt[])) && (pt[[
1]] < (poly[[j, 1]] -
poly[[i, 1]])*(pt[] - 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= {0.062, Null}
Out= {0.016, Null}


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}];][]]
Print["testpoint: ",
Timing[Table[testpoint[poly, pts[[i]]], {i, npts}];][]]
Print["testpoint2: ",
Timing[Table[testpoint2[poly, pts[[i]]], {i, npts}];][]]
Print["inPolyQ: ",
Timing[Table[inPolyQ[poly, pts[[i]]], {i, npts}];][]]
Print["InsidePolygonQ: ",
Timing[Table[InsidePolygonQ[poly, pts[[i]]], {i, npts}];] []]
Print["inPolyQ2: ",
Timing[Table[
inPolyQ2[poly, pts[[i, 1]], pts[[i, 2]]], {i, npts}];][]]


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}];][]]
Print["testpoint: ",
Timing[Table[testpoint[poly, pts[[i]]], {i, npts}];][]]
Print["testpoint2: ",
Timing[Table[testpoint2[poly, pts[[i]]], {i, npts}];][]]
Print["inPolyQ: ",
Timing[Table[inPolyQ[poly, pts[[i]]], {i, npts}];][]]
Print["inPolyQ2: ",
Timing[Table[
inPolyQ2[poly, pts[[i, 1]], pts[[i, 2]]], {i, npts}];][]]


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.

• Mac, thanks for contributing. Looks interesting! I recommend that you localize i and j. Aug 14, 2012 at 11:26
• Thanks for the suggestion - have updated the code accordingly with i and j now localised.
– Mac
Aug 14, 2012 at 11:48
• I'd replace the For[] with a Do[] myself... Aug 14, 2012 at 13:21
• @Mac, you could use a slightly larger polygon (e.g. the polygon from a country from CountryData) and some Random points (with a SeedRandom) and see how that performs ;-)
– user21
Aug 14, 2012 at 14:21
• @mac p is not right, should be p = CountryData["Canada", "Polygon"][[1, 1,25]];. Because there are 25 parts in Canada. The main part is 25th. Tested in V10 Oct 11, 2015 at 6:59

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.

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


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}] • Damn, I was about to post this myself... :D Aug 14, 2012 at 10:22
• For speed I would try: Xip1 = RotateLeft@Xi; Yip1 = RotateLeft@Yi; and Total[u (1 - v) (1 - w) - (1 - u) v w] Aug 15, 2012 at 0:22
• I'll just note that inPolyQ2 is exactly what you'll end up with if you try to vectorize Sunday's implementation of the winding number algorithm. Jan 15, 2017 at 14:13

Another approach to this problem is computing the winding number by integrating $1/z$ centered on the point of interest along the contour of the polygon in the complex plane. Sure this isn't exactly efficient, but still I think it's nice to see this working in action. And since complex integration is feasible in Mathematica, I just tried :)

PointToComplex[{x_, y_}] := x + I y
Windingnumber[polygon_, point_] := Module[{wn,z},
Off[NIntegrate::ncvb, NIntegrate::slwcon];
wn = Round@
Re@Chop[1/(2 π I)
NIntegrate[1/(z - PointToComplex[point]),
Evaluate@{z, Sequence @@ (PointToComplex /@ Append[#, #[]]&[polygon])}]];
On[NIntegrate::ncvb, NIntegrate::slwcon];
wn
]
InsidePolygonQ[polygon_, point_] := Windingnumber[polygon, point] != 0

• You can use Complex @@@ r in place of PointToComplex /@ polygon, etc., I believe. Aug 14, 2012 at 13:40
• @Mr. Wizard, if they're explicitly numbers (i.e. passes NumberQ[] but not NumericQ[]). Complex @@ {1, Pi} won't work, for instance. Aug 14, 2012 at 15:34
• I think this is the first time that code has made me laugh (at least for positive reasons), so +1!
– acl
Aug 14, 2012 at 15:57
• @Mr.Wizard Thanks for the hint! Applying Complex didn't cross my mind when i wrote the code. @J.M. interesting detail, never looked at how Mathematica handles complex number representation! Aug 15, 2012 at 8:32
• @acl Ha, glad to hear that :) was definitely worth posting then ;) Aug 16, 2012 at 22:48

You could use this package to triangulate your polygon, and then use this barycentric formula on each of the triangles.

inside[{{x1_, y1_}, {x2_, y2_}, r3 : {x3_, y3_}}, r : {_, _}] :=
# >= 0 && #2 >= 0 && # + #2 < 1 & @@
LinearSolve[{{x1 - x3, x2 - x3}, {y1 - y3, y2 - y3}}, r - r3]


Example for a single triangle:

tri = {{13.2, 11.9}, {10.3, 12.3}, {9.5, 14.9}};

{
LocatorPane[Dynamic @ pt, Graphics @ {Orange, Polygon@tri}],
Dynamic @ inside[tri, pt]
} Example for a polygon:

<< PolygonTriangulationSimplePolygonTriangulation

poly = {{4.4, 14}, {6.7, 15.25}, {6.9, 12.8}, {9.5, 14.9}, {13.2,
11.9}, {10.3, 12.3}, {6.8, 9.5}, {13.3, 7.7}, {0.6, 1.1}, {1.3,
2.4}, {2.45, 4.7}};

tris = poly[[#]] & /@ SimplePolygonTriangulation[poly];

colors = MapIndexed[{ColorData @ #2[], Polygon@#} &, tris];

DynamicModule[{pt},
{LocatorPane[Dynamic[pt], colors // Graphics],
Or @@ (inside[#, pt] & /@ tris) // Dynamic}
] • Didn't notice until now, but Inverse[{{x1 - x3, x2 - x3}, {y1 - y3, y2 - y3}}].(r - r3) is better done as LinearSolve[{{x1 - x3, x2 - x3}, {y1 - y3, y2 - y3}}, r - r3]. Aug 31, 2012 at 15:10
• @J.M. good to know; thanks! Sep 1, 2012 at 19:57

As per Szabolcs's suggestion:

Version 10 alternatives are RegionMember and Element, but the latter is unreasonably slow.

### A drop in alternative

RegionMember[reg] returns a RegionMemberFunction[...] that can be applied repeatedly to different points.

(* Memoizing the RegionMemberFunction[...] for a given polygon *)
inPolyQHelper[poly_] := inPolyQHelper[poly] = RegionMember[Polygon@poly];
inPolyQ[poly_, pt_] := inPolyQHelper[poly]@pt


### A faster alternative

RegionMember also accepts a list of points to be tested!

RegionMember[Polygon@list, data]


### Benchmarks

data = Table[{RandomReal[{-10, 10}], RandomReal[{-10, 10}]}, {i, 1, 1000000}];
list = {{0.5735,5.274},{-4.961,2.333},<<10>>,{-1.662,-0.1829}};

(* Compiled version from @Simon Wood's answer *)
inPolyQSimonWoods[list, Sequence @@ #] & /@ data // AbsoluteTiming // First
(* 11.465298 *)

(* The drop-in RegionMember replacement *)
inPolyQ[list, #] & /@ data // AbsoluteTiming // First
(* 2.994139 *)

(*The fast replacement*))
RegionMember[Polygon@list, data] // AbsoluteTiming // First
(* 0.399948 *)


Just for the record, Element[#, Polygon @ list] /@ data takes 13 seconds with only 100 points.

• I wish I could vote more. I think RegionMember is by far the fastest based on my test. Who stands with me ! I wish you could add a "Canada test" as Daniel Lichtblau has done. That will better show the power of RegionMember when dealing with complex polygons and lots of points. With a pregenerated RegionFunction at first, RegionMember beats other method badly. Oct 11, 2015 at 10:20

Here is the code from a MathGroup post I had referenced. I have modified to compile to C and that speeds it further. The one-off preprocessing does take time but it seems not unreasonable. It takes a list of lists of polygons (so the "region" need not be connected). To account for this I slightly alter the setup from Mac's response.

Preprocessing the polygon:

getSegsC =
Compile[{{j, _Integer}, {minx, _Real}, {len, _Real}, {eps, _Real}, \
{segs, _Real, 3}}, Module[{lo, hi}, lo = minx + (j - 1)*len - eps;
hi = minx + j*len + eps;
Select[segs,
Module[{xlo, xhi}, {xlo, xhi} = Sort[{#[[1, 1]], #[[2, 1]]}];
lo <= xlo <= hi ||
lo <= xhi <= hi || (xlo <= lo && xhi >= hi)] &]]];

polyToSegmentList[poly_, nbins_] :=
Module[{xvals, yvals, minx, maxx, miny, maxy, segments, flatsegments,
segmentbins, xrange, len, eps}, {xvals, yvals} =
Transpose[Flatten[poly, 1]];
{minx, maxx} = {Min[xvals], Max[xvals]};
{miny, maxy} = {Min[yvals], Max[yvals]};
segments = Map[Partition[#, 2, 1, {1, 1}] &, poly];
flatsegments = Flatten[segments, 1];
xrange = maxx - minx;
eps = 1/nbins*len;
len = xrange/nbins;
segmentbins =
Table[getSegsC[j, minx, len, eps, flatsegments], {j, nbins}];
{{minx, maxx}, {miny, maxy}, segmentbins}]


The actual in-or-out code.

pointInPolygon[{x_, y_}, bins_, xmin_, xmax_, ymin_, ymax_] :=
Catch[Module[{nbins = Length[bins], bin},
If[x < xmin || x > xmax || y < ymin || y > ymax, Throw[False]];
bin = Ceiling[nbins*(x - xmin)/(xmax - xmin)];
If[EvenQ[countIntersectionsC[bins[[bin]], x, y, ymin - 1.]], False,
True]]]

countIntersectionsC =
Compile[{{segs, _Real, 3}, {x, _Real}, {yhi, _Real}, {ylo, _Real}},
Module[{tally = 0, yval, xlo, xhi, y1, y2},
Do[{{xlo, y1}, {xhi, y2}} = segs[[j]];
If[(x < xlo && x < xhi) || (x > xlo && x > xhi), Continue[]];
yval = y1 + (x - xlo)/(xhi - xlo)*(y2 - y1);
If[ylo < yval < yhi, tally++];, {j, Length[segs]}];
tally]];


The mainland of Canada will be the test again. As in Mac's example I rescale so coordinates are all between -1 and 1. This means I really don't need the x/ymin/max stuff but I opted to keep that in.

p = CountryData["Canada", "Polygon"][[1, 1]];
poly = {Transpose[{Rescale[
p[[All, 1]], {Min@#, Max@#} &@p[[All, 1]], {-1, 1}],
Rescale[p[[All, 2]], {Min@#, Max@#} &@p[[All, 2]], {-1, 1}]}]};


I'll use 1000 bins and do the preprocessing.

nbins = 1000;
Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} =
polyToSegmentList[poly, nbins];]

(* Out= {5.15, Null} *)


For testing I'll take 10000 points.

npts = 10000;
pts = Partition[RandomReal[{-1, 1}, 2*npts], 2];

Timing[
inout = Map[pointInPolygon[#, segmentbins, xmin, xmax, ymin, ymax] &,
pts];]

(* Out= {0.37, Null} *)


Visual check:

ListPlot[Pick[pts, inout], Joined -> False] The northeast reminds me a bit of the duck's head seen here. But then...I've always found the Baffin...to be bafflin'.

• For once Canada is mentioned in another context than exporting cold air down south...(Canadians will know what I mean).
– Mac
Aug 16, 2012 at 9:17
• I don't think we will find a better solution in terms of performance unless perhaps a GPU-solution can be implemented (way above my programming skills in any case). One minor point - you're actually testing 10000 pairs of coordinates (points) not 20000.
– Mac
Aug 16, 2012 at 9:27
• Thanks. I changed from my point generator to yours and failed to catch the effect of the partitioning. Will correct response accordingly.. Aug 16, 2012 at 14:33
• Worth pointing out too that the "Canada" polygon example assumes implicitly that each vertex defining the polygon can be connected by a straight line. This has to be a good assumption for such a densely sampled contour. However, for large geographical polygons with few points the connecting edges should be modelled as "great circles".
– Mac
Aug 20, 2012 at 7:25
• @maththeorem It is very difficult to future-proof code of this sort. It appears that CountryData has changed since this thread was in its infancy. Oct 11, 2015 at 16:22

In version 10 (now available through the Programming Cloud) it is now possible to simply use Element:

For example,

Element[{0,0}, Polygon[{{-1,-1},{-1,1},{1,1},{1,-1}}]]

(* True *)


This works for arbitrary regions in general. Most graphics primitives can be used as regions.

• Use RegionMember! It is much faster, specially if you pass the list of points to be tested instead of Map'ing it over them! Nov 9, 2014 at 17:28
• @Aisamu Sounds like a much better solution, why don't you post it as an answer? Nov 9, 2014 at 17:32
• I thought there were too many answers already... but I took your advice and just posted it! Nov 9, 2014 at 18:05

Another approach you could use is to draw a line (or define a vector) between a line guaranteed to be outside the polygon and the point you wish to test, then counting the number of line segments of the polygon that intersect with this line. If this number is odd, the point is inside the polygon.

To determine if two line segments intersect, you can use the vector algebra from this SO answer: How do you detect where two line segments intersect?. The short of it is that for any two vectors that intersect, there are two scalars that can be applied, one to each vector, to produce a parallel vector of the exact magnitude needed to reach the intersection. These scalars are a function of the cross product of the vectors. If both scalars are $0 < x < 1$ then this intersection happens within the magnitudes of the original vectors. If $x > 1$ or $x < 0$ for either scalar, they intersect beyond the bounds of the defined vectors, while if $x=0$ the vectors are parallel.

This test should be linear to the number of points defining the polygon (requiring a scan of all points to determine the max X-coord and y-coord to produce a point outside the polygon, and then a scan of all adjacent pairs of points to produce line segments followed by constant-time operations to determine intersection). And, it should work with any 2D polygon you can imagine, no matter how twisted.

• This is the gist of what I had here. In that case preprocessing, which does take real time, makes the individual queries closer to constant time for "typical" polygons. Aug 15, 2012 at 18:45
• the trick in real-world use of this algorithm is the boundary cases - what do you do when the imaginary ray to infinity exactly crosses the intersection of two lines, or is exactly coincident with one of the line segments. Aug 15, 2012 at 23:22
• @ddyer In actual practice I would probably choose a random direction to approach from. I might also change the binning to use a random orientation rather than horizontal. I might also add a check for hitting a vertex exactly. Depends on how concerned I was about getting a wrong result in a "small" set of cases. Aug 16, 2012 at 14:39
• sure, all those things are potential strategies. My point is that a conceptually elegant algorithm has lots of gritty details to be handled if you want to make it work reliably. Aug 16, 2012 at 18:19
• @ddyer I actually do this stuff for a living... Aug 22, 2012 at 16:36

Sorry to be late to the party. I'll throw in the following Mathematica implementation of an algorithm by W. Randolph Franklin which I wrote up here a while ago.

The implementation has a number of nice features:

• Polygon can be closed or not.
• A point will be inside exactly one member of a polygonal partitioning.
• No trigonometry, so it's blazing fast.

pnPoly[{testx_, testy_}, pts_List] := Xor @@ ((
Xor[#[[1, 2]] > testy, #[[2, 2]] > testy] &&
((testx - #[[2, 1]]) < (#[[1, 1]] - #[[2, 1]]) (testy - #[[2, 2]])/(#[[1, 2]] - #[[2, 2]]))
) & /@ Partition[pts, 2, 1, {2, 2}])

• The speed will depend linearly on the number of vertices. That is often fine but could be a problem if there are both many vertices and many query points. Jan 20, 2014 at 22:53
• @DanielLichtblau: Yes, you are right of course that for a large polygon you want to do something hierarchical along the lines of your answer to get decent scaling. One reason I keep coming back to this implementation is the partitioning guarantee which is critical in much of what I do. Jan 21, 2014 at 8:41
• I had a look at "Insignificance Galore" where it mentions the partitioning guarantee. But I still do not understand what it means. Is it for the case of multiple disconnected polygons? Self-intersecting? Or does it also have meaning in the case of one non-self-intersecting polygon. Jan 21, 2014 at 15:29
• A partitioning of a set S is a collection of disjoint subsets of S whose union is S mathworld.wolfram.com/SetPartition.html. The practical problem with partitioning (part of) the plane into polygons is to specify what happens to points on the edges and vertices: it's a lot of tedious details which are usually unimportant from a mathematical point of view (since the combined edges have 0 area), but still needs to be done right for some numerical algorithms to work. Jan 22, 2014 at 8:45
• Okay, thanks for the explanation. I will add that it is also critical, in polynomial irreducibility testing, to know if an exponent vector is or is not a vertex in the convex hull corresponding to a certain a Newton polytope. I can say that numerical convex hull methods have made such determination much more difficult than I would like. So there is at least one math algorithm where this does matter. Jan 22, 2014 at 15:51