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

87

56

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.

Notebook from Mathematica-users.org PointInsidePolygon.nb contains additional solutions.

– Alexey Popkov – 2013-10-11T00:14:01.300

1

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

– jVincent – 2012-08-14T08:39:00.740

– kglr – 2012-08-14T08:49:56.780

@jVincent, I was hoping someone was willing to share their Mathematica code. – nilo de roock – 2012-08-14T08:52:17.610

@kguler Interesting stuff in the answers. Thanks, I'll have a closer look. – nilo de roock – 2012-08-14T08:53:31.913

2It's pretty depressing that you already accepted an answer before I posted mine, but I posted anyway. :-/ – Mr.Wizard – 2012-08-14T10:02:55.690

Anyway: there are a number of methods to choose from here (also here). Now, adapting them to Mathematica is a different kettle of fish...

– J. M.'s ennui – 2012-08-14T13:27:58.140

4The 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. – whuber – 2012-08-14T17:52:25.653

4

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. http://forums.wolfram.com/mathgroup/archive/2009/Feb/msg00519.html

– Daniel Lichtblau – 2012-08-14T21:10:18.463

52

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*)


@kguler Hello! I tested GraphicsMeshInPolygonQ and unfortunately found it didn't work in V10. What is more, GraphicsMeshPointWindingNumber proposed by rm-rf also didn't work in V10 – matheorem – 2014-12-03T01:07:24.113

@kguler And for the first version of testpoint, sometimes it encounters error message "ArcTan::indet: "Indeterminate expression ArcTan[0.,0.] encountered"", if the edge of polygon passes points. – matheorem – 2014-12-03T01:10:48.680

@matheorem, I think with the new version 10 built-ins we no-longer need these methods; the new built-in RegionMember is much faster and Element is much more elegant. – kglr – 2014-12-03T01:23:42.357

@kguler Wow, you're so familiar with new features in v10. I was wondering if you already known RegionMember before I posted my question? – matheorem – 2014-12-03T01:28:15.487

@matheorem, i learned about Element and RegionMember from the answers posted here sometime back: Aisamu posted this answer below using RegionMember and Szabolcs had posted this one using Element.

– kglr – 2014-12-03T01:34:05.887

@kguler Ok. By the way, I just tested RegionMember, but I found it much slower than your original testpoint version, and also much slower than GraphicsMeshPointWindingNumber when I used it in v9 – matheorem – 2014-12-03T01:37:46.413

@kguler Another question, where is GraphicsMeshPointWindingNumber now in V10? Deleted? – matheorem – 2014-12-03T01:39:06.630

@matheorem, I don't have v10. Using the free version of v10 on wolfram cloud, it seems that it moved to a new package GraphicsPolygonUtilsPointWindingNumber. – kglr – 2014-12-03T01:45:05.210

@kguler Wow, thank you! Much faster now, though GraphicsPolygonUtilsPointWindingNumber doesn't include points threaded by polygon edge. Could please show me how can you find the package moved to PolygonUtils? Because if I search PointWindingNumber in help documentation, it doesn't give any result. – matheorem – 2014-12-03T01:59:58.543

1@matheorem, i used ??**PointWindingNumber* – kglr – 2014-12-03T02:02:55.187

@kguler Please forgive my verbose comment, I know ?? is for Information, but what does those "" mean here? And when I type ??*PointWindingNumber*, v10 gives both "GraphicsMesh" and "GraphicsPolygonUtils", so why "GraphicsMesh" didn't work any more? – matheorem – 2014-12-03T02:11:14.547

@matheorem, * is a the usual wildcard. I only get GraphicsPolygonUtilsPointWindingNumber on the cloud version. No idea why GraphicsMeshPointWindingNumber does not work. – kglr – 2014-12-03T02:15:51.037

@kguler Oh, wildcard! I understand! And this is my v10 version gives http://imgur.com/bCIANiR, it is quite odd that mesh didn't work.

– matheorem – 2014-12-03T02:27:08.967

1@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 – 2014-12-03T14:07:11.273

@rm-rf Thank you! You're right. – matheorem – 2014-12-03T23:19:47.680

Looks great! I'll try to implement this in my app and let you know asap. Thank you very much. – nilo de roock – 2012-08-14T09:14:12.927

The testpoint function is a beauty. I would give a bonus if possible, I'll do a fast accept, instead. Thanks kguler – nilo de roock – 2012-08-14T09:57:23.247

1This doesn't work for self-intersecting polygons, but you could use PolygonTessellation from the package linked in my answer to pre-process. +1 – Mr.Wizard – 2012-08-14T10:39:54.600

@Mr.Wizard, good point. Will update with the needed pre-processing. ndroock1, pls re-consider your accept. – kglr – 2012-08-14T11:09:26.887

83

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.

@rm-rf As I tested it, this undocumented GraphicsMeshPointWindingNumber is not working in V10, so maybe undocumented function is dangerous? – matheorem – 2014-12-03T00:43:08.707

@matheorem The inner context was changed from Mesh to PolygonUtils. I knew about this before, but Mathematica 10 was in private beta and I couldn't post it then. I've updated the post now. – rm -rf – 2014-12-03T03:12:01.710

1@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? – matheorem – 2014-12-03T06:35:38.290

7One 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 – 2012-08-14T16:37:15.990

Depending on how you define the interior of the polygon, you may want to check if the winding number is odd instead. – Solomon Ucko – 2020-03-11T21:13:45.040

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 – 2013-04-19T08:16:50.263

2

@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 – 2013-04-19T13:41:39.610

27

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 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 – matheorem – 2015-10-11T06:59:46.563

Mac, thanks for contributing. Looks interesting! I recommend that you localize i and j. – Mr.Wizard – 2012-08-14T11:26:48.813

Thanks for the suggestion - have updated the code accordingly with i and j now localised. – Mac – 2012-08-14T11:48:02.267

I'd replace the For[] with a Do[] myself... – J. M.'s ennui – 2012-08-14T13:21:06.977

2@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 ;-) – None – 2012-08-14T14:21:46.277

26

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 – J. M.'s ennui – 2012-08-14T10:22:14.673

For speed I would try: Xip1 = RotateLeft@Xi; Yip1 = RotateLeft@Yi; and Total[u (1 - v) (1 - w) - (1 - u) v w] – Mr.Wizard – 2012-08-15T00:22:19.250

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.

– J. M.'s ennui – 2017-01-15T14:13:22.247

24

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


1You can use Complex @@@ r in place of PointToComplex /@ polygon, etc., I believe. – Mr.Wizard – 2012-08-14T13:40:59.240

2@Mr. Wizard, if they're explicitly numbers (i.e. passes NumberQ[] but not NumericQ[]). Complex @@ {1, Pi} won't work, for instance. – J. M.'s ennui – 2012-08-14T15:34:13.607

1I think this is the first time that code has made me laugh (at least for positive reasons), so +1! – acl – 2012-08-14T15:57:05.107

@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! – Thies Heidecke – 2012-08-15T08:32:03.307

@acl Ha, glad to hear that :) was definitely worth posting then ;) – Thies Heidecke – 2012-08-16T22:48:38.413

22

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}
] 3Didn'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]. – J. M.'s ennui – 2012-08-31T15:10:14.900

@J.M. good to know; thanks! – Mr.Wizard – 2012-09-01T19:57:00.060

18

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. – matheorem – 2015-10-11T10:20:48.913

17

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'.

@DanielLichtblau p should be p = CountryData["Canada", "Polygon"][[1, 1,25]];. The main part of Canada is 25th list. – matheorem – 2015-10-11T07:03:12.440

1@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. – Daniel Lichtblau – 2015-10-11T16:22:06.910

For once Canada is mentioned in another context than exporting cold air down south...(Canadians will know what I mean). – Mac – 2012-08-16T09:17:13.783

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 – 2012-08-16T09:27:18.513

Thanks. I changed from my point generator to yours and failed to catch the effect of the partitioning. Will correct response accordingly.. – Daniel Lichtblau – 2012-08-16T14:33:00.280

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 – 2012-08-20T07:25:56.713

15

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.

3Use RegionMember! It is much faster, specially if you pass the list of points to be tested instead of Map'ing it over them! – Aisamu – 2014-11-09T17:28:30.830

@Aisamu Sounds like a much better solution, why don't you post it as an answer? – Szabolcs – 2014-11-09T17:32:06.663

11

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.

– Daniel Lichtblau – 2012-08-15T18:45:53.533

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. – ddyer – 2012-08-15T23:22:12.087

@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. – Daniel Lichtblau – 2012-08-16T14:39:42.987

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. – ddyer – 2012-08-16T18:19:06.727

3@ddyer I actually do this stuff for a living... – Daniel Lichtblau – 2012-08-22T16:36:37.620

8

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. – Daniel Lichtblau – 2014-01-20T22:53:02.007

@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. – Janus – 2014-01-21T08:41:19.350

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. – Daniel Lichtblau – 2014-01-21T15:29:40.867

A partitioning of a set S is a collection of disjoint subsets of S whose union is S http://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.

– Janus – 2014-01-22T08:45:20.287

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. – Daniel Lichtblau – 2014-01-22T15:51:07.843