Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

testing if a point is inside a polygon

339 views
Skip to first unread message

Mitch Murphy

unread,
Feb 9, 2009, 5:31:57 AM2/9/09
to

is there a way to test whether a point is inside a polygon? ie.

PointInsidePolygonQ[point_,polygon_] -> True or False

i'm trying to do something like ...

ListContourPlot[table,RegionFunction->CountryData["Canada","Polygon"]]

to create what would be called a "clipping mask" in photoshop.

cheers,
Mitch

Frank Scherbaum

unread,
Feb 10, 2009, 5:47:33 AM2/10/09
to
Mitch,

please find below what I use for his purpose. I hope it is useful.
Best regards,
Frank

PointInPolygonQ::usage="PointInPolygonQ[pt, poly] uses the winding-
number algorithm (Godkin and Pulli, 1984) to check, if point pt is
inside the closed polygon poly, which is given as list of its
vertices."
(*
checks, if a point is inside a polygon
pt: point as {lat [deg], lon [deg]} to test
poly: list of polygon vertices coordinates

GODKIN,C.B. AND J.J.PULLI: APPLICATION OF THE "WINDING-NUMBER
ALGORITHM" TO THE SPATIAL SORTING OF CATALOGED EARTHQUAKE DATA.
Bull. Seismol. Soc. Am. 74, 5, PP. 1845-1848, OCTOBER 1984

RETURN VALUE: 0 IF POINT OUTSIDE
+/-1 IF POINT INSIDE
2 IF POINT IS ON AN EDGE OR VERTEX

*)
PointInPolygonQ[pt_,poly_] := Module[
{ i,n,isicr,inside,px,py,pxx,pyy,x0,y0 },
n = Length[poly];
(* ACCUMULATE SIGNED CROSSING NUMBERS WITH INSIDE *)
inside = 0;
{x0,y0}=pt;

For[i=1,i < n,i++,
(* PROCEED AROUND POLYGON CHECKING EACH SEGMENT TO SEE IF
NEGATIVE X-AXIS WAS CROSSED
TRANSLATE COORDINATES OF POLYGON TO PUT TEST POINT AT
ORIGIN *)
{px,py} = poly[[i]];
{pxx,pyy} = poly[[i+1]];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
(* STOP COUNTING IF POINT IS ON EDGE *)
If[isicr == 4, Return[2]];
inside += isicr;
];
(* CHECK SEGMENT FROM LAST VERTEX TO FIRST VERTEX *)
{px,py} = poly[[n]];
{pxx,pyy} = poly[[1]];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
If[isicr == 4, Return[2]];
inside = (inside + isicr)/2;
Return[inside];
];
(*
COMPUTE SIGNED CROSSING NUMBER

A "SIGNED CROSSING NUMBER" TELLS WETHER A SEGMENT
(IN THIS CASE THE SEGMENT FROM (X1,Y1) TO (X2,Y2))
CROSSES THE NEGATIVE X-AXIS OR GOES THROUGH THE ORIGIN

THE RETURN VALUES ARE:
+2 IF SEGMENT CROSSES FROM BELOW
+1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM BELOW OR STARTS
UPWARDS FROM -X-AXIS ("HALF CROSSING")
0 IF THERE IS NO CROSSING
-1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM ABOVE OR STARTS
DOWNWARDS FROM -X-AXIS ("HALF CROSSING")
-2 IF SEGMENT CROSSES FROM ABOVE
+4 IF SEGMENT CROSSES THROUGH THE ORIGIN

*)

ksicr[x1_,y1_,x2_,y2_] := Module[
{
},
(* IF BOTH POINTS ARE ON THE SAME SIDE OF X-AXIS, RETURN 0 *)
If[N[y1*y2 > 0.], Return[0] (* no crossing *)];

(* CHECK IF SEGMENT CROSSES THROUGH THE ORIGIN *)
If[x1*y2 != x2*y1 || x1*x2 > 0.,
If[y1 * y2 < 0,
(*
COMPLETE CROSSING OF -X-AXIS?
BREAK INTO CASES ACCORDING TO CROSSING DIRECTION
*)
If[y1 > 0,
(* CASE Y1 > 0 > Y2 *)
If[y1 * x2 >= x1 * y2, Return[0];, (* no crossing *)
Return[-2]; (* downward crossing *)
];
,
(* CASE Y1 < 0 < Y2 *)
If[x1 * y2 >= y1 * x2, Return[0];, (* no crossing *)
Return[2]; (* upward crossing *)
];
];
,
(*
HALF CROSSING?
ONE END OF SEGMENT TOUCHES X-AXIS! WHICH END?
*)
If[y2 == 0,
(* HERE Y2==0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[y1 == 0 || x2 > 0,
Return[0]; (* no crossing *)
,
(* UPWARD OR DOWNWARD? *)
If[y1 > 0.,
Return[-1]; (* Downward half crossing *)
,
Return[1]; (* Upward half crossing *)
];
];
,
(* HERE Y1==0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return[0];
,
(* UPWARD OR DOWNWARD? *)
If[y2 > 0,
Return[1]; (* Upward half crossing *)
,
Return[-1]; (* Downward half crossing *)
];
];
];

(* HERE Y1=0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return[0]; (* no crossing *)
];
(* UPWARD OR DOWNWARD? *)
If[y2 > 0.,
Return[-1]; (* Downward half crossing *),
Return[1]; (* Upward half crossing *)
];
];
,
Return[4];
];
];


Jens-Peer Kuska

unread,
Feb 10, 2009, 5:48:26 AM2/10/09
to
Hi,

the web is full of this test as well as any graphics
book. Look at
http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/
http://tog.acm.org/editors/erich/ptinpoly/

Regards
Jens

David Park

unread,
Feb 10, 2009, 5:50:35 AM2/10/09
to
Mitch,

Here is my attempt. The idea is to move around the polygon and add up the
angles made with the point. If the point is outside this will be zero. The
steps are basically:

1) Check if the point is one of the vertices of the polygon.
2) Add the first polygon point to the end of the list of points.
3) Subtract the test point from each of the polygon vertices.
4) Partition the points into pairs.
5) Rotate each pair so the first one lies along the x axis. (To prevent
crossing the branch line)
6) Find the ArcTan of the second point. This will have a sign.
7) Sum the angles and take the absolute value.
8) Test if this is zero or some multiple of 2\[Pi].

PointInsidePolygonQ::usage =
"PointInsidePolygonQ[point,polygon] will return True if the point \
is on the boundary or inside the polygon and False otherwise.";

PointInsidePolygonQ[point_, polygon_] :=
Module[{work},
work = If[Head[polygon] === Polygon, First[polygon], polygon];
If[ \[Not] FreeQ[work, point], Return[True]];
work = # - point & /@ Join[work, {First[work]}];
work = Partition[work, 2, 1];
work = RotationTransform[{First[#], {1, 0}}]@# & /@ work;
work = Abs@Total[ArcTan @@ Last[#] & /@ work] // Chop // Rationalize;
TrueQ[work != 0]
]

Here is a graphical test for a simple polygon:

testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[4],
{If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
testpoints},
PlotRange -> 10,
Frame -> True]

Here is a test for a more complex polygon.

testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,
1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677,
0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],
AbsolutePointSize[4],
{If[PointInsidePolygonQ[#, Polygon[polypoints]], Black, Red],
Point[#]} & /@ testpoints},
PlotRange -> 10,
Frame -> True]

There might be simpler methods.

David Park
djm...@comcast.net
http://home.comcast.net/~djmpark/

Jean-Marc Gulliet

unread,
Feb 10, 2009, 5:49:09 AM2/10/09
to
In article <gmp0mt$btn$1...@smc.vnet.net>, Mitch Murphy <mi...@lemma.ca>
wrote:

> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False

[snip]

You could have a look at the different algorithms presented in Paul
Bourke's paper [1], pick up one of your liking, then code it from the
provided C source code to Mathematica syntax. How the algorithms work,
their benefits and draw backs, as well as tested C code and diagrams,
are neatly presented.

Depending of your level of familiarity with Mathematica programming, it
might be easier to translate the C syntax into equivalent Mathematica
procedural constructs, test the code, and then convert it into fast
functional Mathematica code.

If you have performance issues or troubles converting procedural code
(fast in C, but slow in Mathematica) into functional code (fast in
Mathematica), people in MathGroup will help you to improve performances
or overcome coding issues.

Regards,
--Jean-Marc

[1] Paul Bourke, "Determining If A Point Lies On The Interior Of A
Polygon",
http://local.wasp.uwa.edu.au/~pbourke/geometry/insidepoly/

dh

unread,
Feb 10, 2009, 5:52:12 AM2/10/09
to

Hi Mitch,

To write a routines that checks if a point is inside a polygon you must

realize that if you walk around the polygon you see a inner point

always on the same side. This is false for a outside point. Therefore,

the test is to check if the point in question: p0 is always on the same

side of the polygon sides.

Assume {p1,p1,..} are the polygon points, arrange in a definite order

(clockwise or anti-clockwise). How can we check on which side of {pi+1,

pi} p0 is? We may get a vector perpendicular to the side by rotating by

Pi/2: {{0,1},{-1,0}}.(pi+1 - pi). If we take the dot product with (p0

-pi) the sign will tell us the side. Here is an example:

pts = {{0, 0}, {1, 0}, {1, 1}, {0, 1}};

p0 = {0.5, 0.5};

sides= Subtract @@ # & /@ Partition[Append[pts, pts[[1]]], 2, 1];

perp= {{0,1},{-1,0}}.#& /@ sides;

Equal @@ Sign @ MapThread[Dot, {p0 - # & /@ pts, perp}]

The third lines calculates (pi+1 - pi), note that we added to first

point to the end of the list to close the polygon. The fourth turns each

vector by Pi/2. The fifth calculates p0-pi = {p0 - # & /@ pts and takes

the Dot product of corresponding vectors. Finally it takes the sign and

checks if all signs are the same.

hope this helps, Daniel

Sjoerd C. de Vries

unread,
Feb 10, 2009, 5:52:55 AM2/10/09
to
Here is a function based on the code in Robert Sedgewick's Algorithms
book:

PointInPoly[{x_, y_}, poly_List] :=
Module[{i, j, c = False, npol = Length[poly]},
For[i = 1; j = npol, i <= npol, j = i++,
If[((((poly[[i, 2]] <= y) && (y <
poly[[j, 2]])) || ((poly[[j, 2]] <= y) && (y <
poly[[i, 2]])))
&& (x < (poly[[j, 1]] -
poly[[i, 1]])*(y - poly[[i, 2]])/(poly[[j, 2]] -
poly[[i, 2]]) + poly[[i, 1]])), c = \[Not] c
];
];
Return [c];
]

Cheers -- Sjoerd

On Feb 9, 12:31 pm, Mitch Murphy <mi...@lemma.ca> wrote:
> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> i'm trying to do something like ...
>

> ListContourPlot[table,RegionFunction->CountryData["Canada=

Sjoerd C. de Vries

unread,
Feb 11, 2009, 5:20:15 AM2/11/09
to
Nice one David, but it's about 50 times slower than using Sedgewick's
algorithm...

Cheers -- Sjoerd

On Feb 10, 12:50 pm, "David Park" <djmp...@comcast.net> wrote:
> Mitch,
>
> Here is my attempt. The idea is to move around the polygon and add up the

> angles made with the point. If the point is outside this will be zero. Th=


e
> steps are basically:
>
> 1) Check if the point is one of the vertices of the polygon.
> 2) Add the first polygon point to the end of the list of points.
> 3) Subtract the test point from each of the polygon vertices.
> 4) Partition the points into pairs.
> 5) Rotate each pair so the first one lies along the x axis. (To prevent
> crossing the branch line)
> 6) Find the ArcTan of the second point. This will have a sign.
> 7) Sum the angles and take the absolute value.
> 8) Test if this is zero or some multiple of 2\[Pi].
>
> PointInsidePolygonQ::usage =
> "PointInsidePolygonQ[point,polygon] will return True if the point \
> is on the boundary or inside the polygon and False otherwise.";
>
> PointInsidePolygonQ[point_, polygon_] :=
> Module[{work},

> work = If[Head[polygon] === Polygon, First[polygon], polygon]=


;
> If[ \[Not] FreeQ[work, point], Return[True]];
> work = # - point & /@ Join[work, {First[work]}];
> work = Partition[work, 2, 1];
> work = RotationTransform[{First[#], {1, 0}}]@# & /@ work;

> work = Abs@Total[ArcTan @@ Last[#] & /@ work] // Chop // Rationaliz=


e;
> TrueQ[work != 0]
> ]
>
> Here is a graphical test for a simple polygon:
>
> testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
> polypoints = {{-1.856, 3.293}, {1.257,
> 2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
> Graphics[
> {Lighter[Green, .8], Polygon[polypoints],
> AbsolutePointSize[4],
> {If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
> testpoints},
> PlotRange -> 10,
> Frame -> True]
>
> Here is a test for a more complex polygon.
>
> testpoints = testpoints = RandomReal[{-9, 9}, {100, 2}];
> polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,

> 1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677=


,
> 0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
> Graphics[
> {Lighter[Green, .8], Polygon[polypoints],
> AbsolutePointSize[4],
> {If[PointInsidePolygonQ[#, Polygon[polypoints]], Black, Red],
> Point[#]} & /@ testpoints},
> PlotRange -> 10,
> Frame -> True]
>
> There might be simpler methods.
>
> David Park

> djmp...@comcast.nethttp://home.comcast.net/~djmpark/

Adriano Pascoletti

unread,
Feb 11, 2009, 5:16:20 AM2/11/09
to
Mitch,

the problem was discussed by the mathgroup in 2000 and Luc Barthelet has
collected some results at
http://www.mathematica-users.org/webMathematica/wiki/wiki.jsp?pageName=Notebook:PointInsidePolygon.nb


Adriano Pascoletti


2009/2/9 Mitch Murphy <mi...@lemma.ca>

David Park

unread,
Feb 11, 2009, 5:21:20 AM2/11/09
to
Frank,

The routine that you submitted was much faster than the one I submitted.
That was basically because I was using a RotationTransform to obtain the
proper signed rotations for each side of the polygon. However, this can be
done much faster using complex arithmetic. So here is a new routine:

angle[{p1_, p2_}] :=
Module[{c1, c2},
{c1, c2} = #.{1, I} & /@ {p1, p2};
Arg[c2/c1]]

PointInsidePolygonQ::usage =
"PointInsidePolygonQ[point,polygon] will return True if the point \
is on the boundary or inside the polygon and False otherwise.";

SyntaxInformation[
PointInsidePolygonQ] = {"ArgumentsPattern" -> {_, _}};
PointInsidePolygonQ[point_, polygon_] :=
Module[{work = Join[polygon, {First[polygon]}]},


If[ \[Not] FreeQ[work, point], Return[True]];

work = # - point & /@ work;
(Total[angle /@ Partition[work, 2, 1]] // Chop) != 0
]

Here are graphical test routines for a simple polygon and one that folds on
itself.

testpoints = RandomReal[{-9, 9}, {5000, 2}];


polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],

AbsolutePointSize[2],


{If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
testpoints},
PlotRange -> 10,

Frame -> True,
ImageSize -> 400] // Timing

testpoints = testpoints = RandomReal[{-9, 9}, {5000, 2}];


polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,

1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677,


0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],

AbsolutePointSize[2],


{If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
testpoints},
PlotRange -> 10,

Frame -> True,
ImageSize -> 400] // Timing


The following are the comparable graphical routines for the Godkin, Pulli
algorithm below.

testpoints = testpoints = RandomReal[{-9, 9}, {5000, 2}];


polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],

AbsolutePointSize[2],
{If[Abs@PointInPolygonQ[#, polypoints] == 1, Black, Red],

Point[#]} & /@ testpoints},
PlotRange -> 10,

Frame -> True,
ImageSize -> 400] // Timing

testpoints = testpoints = RandomReal[{-9, 9}, {5000, 2}];


polypoints = {{-3.653, 5.329}, {0.2395, 6.168}, {-0.8982,

1.138}, {-0.6587, 1.138}, {5.569, 3.234}, {6.527, -2.036}, {1.677,


0.479}, {-6.407, -1.976}, {-5.21, 2.635}, {1.856, -3.713}};
Graphics[
{Lighter[Green, .8], Polygon[polypoints],

AbsolutePointSize[2],
{If[Abs@PointInPolygonQ[#, polypoints] == 1, Black, Red],

Point[#]} & /@ testpoints},
PlotRange -> 10,

Frame -> True,
ImageSize -> 400] // Timing

The timings appear to be comparable, although the Godkin, Pulli routine
might be slightly faster.

From: Frank Scherbaum [mailto:Frank.S...@geo.uni-potsdam.de]


Mitch,

please find below what I use for his purpose. I hope it is useful.
Best regards,
Frank

Am Feb 9, 2009 um 11:32 AM schrieb Mitch Murphy:

>
> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> i'm trying to do something like ...
>
> ListContourPlot[table,RegionFunction-
> >CountryData["Canada","Polygon"]]
>
> to create what would be called a "clipping mask" in photoshop.
>
> cheers,
> Mitch
>

PointInPolygonQ::usage="PointInPolygonQ[pt, poly] uses the winding-

DrMajorBob

unread,
Feb 11, 2009, 5:22:58 AM2/11/09
to
If I calculate the rotation matrix from one vector to another in spherical
coordinates, the result is HUGELY complicated:

Needs["VectorAnalysis`"]

SetCoordinates[Spherical];
mean = {1, t0, p0};
pt = {1, t, p};
xy = CoordinatesToCartesian /@ {mean, pt};
rot = RotationMatrix@xy;
LeafCount@rot

32798

I don't know what to tell Simplify about this, but it seem there are MANY
unnecessary Conjugate mentions:

Cases[rot, Conjugate[z_], Infinity] // Length

1575

and MANY unnecessary cases like this, too:

Cases[rot, Power[Abs[z_], 2], Infinity] // Length

1980

Each of these could be simply z^2, I think.

Trying to eliminate these with rules doesn't seem to help much, if any.

The first thing I thought of was PowerExpand, but

Map[PowerExpand, rot, {2}] // LeafCount

32798

What am I missing? How can it be that complicated in the first place?

Bobby

DrMajorBob

unread,
Feb 12, 2009, 6:37:57 AM2/12/09
to
I decided a specific "mean" would do, and simplifications gave me the
following:

Needs["VectorAnalysis`"]

SetCoordinates[Spherical];
center = {1, 0, 0};


pt = {1, t, p};

xy = CoordinatesToCartesian /@ {center, pt};
rot = RotationMatrix@xy;
LeafCount@rot

1059

rot //. {Conjugate -> Identity, Abs[z_]^2 :> z^2} // Simplify
% // LeafCount

{{Cos[p]^2 (Cos[t] + Tan[p]^2), -Sin[2 p] Sin[t/2]^2,
Cos[p] Sin[t]}, {-Sin[2 p] Sin[t/2]^2, Cos[p]^2 + Cos[t] Sin[p]^2,
Sin[p] Sin[t]}, {-Cos[p] Sin[t], -Sin[p] Sin[t], Cos[t]}}

80

Bobby

On Wed, 11 Feb 2009 04:23:42 -0600, DrMajorBob <btr...@austin.rr.com>
wrote:

--
DrMaj...@longhorns.com

Sjoerd C. de Vries

unread,
Feb 12, 2009, 6:38:29 AM2/12/09
to
Since timing is often important for functions like this one, I've
compared the various algorithms posted here using David's two tests.
The results are below:

{
{"Daniel", 0.828, 1.219},
{"Daniel/Sjoerd", 0.688, 1.078},
{"Sedgewick/Sjoerd", 0.688, 0.954},
{"David1", 41.766, 81.89},
{"David2", 1.578, 2.531},
{"Godkin/Frank", 1.485, 2.359}
}

In the "Daniel/Sjoerd" method I have used Daniel's code and taken out
every intermediate step. This saves an additional 20%. The code is now
completely gibberish, but it still works ;-)

PointInsidePolygonQ[pt_, poly_] :=
Module[{},
Equal @@
Sign@MapThread[
Dot, {pt - # & /@
poly, {{0, 1}, {-1, 0}}.# & /@ (Subtract @@ # & /@
Partition[Append[poly, poly[[1]]], 2, 1])}]
]

Cheers --Sjoerd

On Feb 9, 12:31 pm, Mitch Murphy <mi...@lemma.ca> wrote:

> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> i'm trying to do something like ...
>

> ListContourPlot[table,RegionFunction->CountryData["Canada=

DrMajorBob

unread,
Feb 12, 2009, 6:40:16 AM2/12/09
to
Just a few caveats:

1) None of these solutions will work for a country that straddles the
Greenwich meridian or a continent that straddles a pole. Antarctica does
both.

CountryData["Antarctica", "Shape"]

To handle this in general, we'd need to rotate the globe so that the
country/continent in question doesn't have that problem... and test points
must be on the same hemisphere, too, to be inside a country.

(This is why I was asking about RotationMatrix.)

2) The poster mentioned Canada, which is made up of disjoint polygons:

Dimensions /@ First@canada;
Length@%
%%[[All, 1]] // Total

25

30458

Twenty-five polygons, more than 30K points... and that's the BASIC
version. Here's the other:

canada = CountryData["Canada", "FullPolygon"];
Dimensions /@ First@canada;
Length@%
%%[[All, 1]] // Total

971

43262

971 polygons, 43 thousand points! For the OP's application, he may need
VERY fast code.

3) In the PointInsidePolygon code, I'd use MemberQ and Append, not FreeQ
and Join, and I NEVER, EVER use Return. (But I found no difference in
speed.)

pointInsidePolygonQ[point_, polygon_] /; MemberQ[polygon, point] = True;
pointInsidePolygonQ[point_, polygon_] :=
0 != Chop@
Total[angle /@
Partition[# - point & /@ Append[polygon, First@polygon], 2, 1]]

4) Points on edges of the polygon test as OUTSIDE:

polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};

PointInsidePolygonQ[(
2 polypoints[[1]] + 8 polypoints[[2]])/10, polypoints]

False

One solution for this would be to stop calculating angles if one of them
happens to be Pi. For instance,

Clear[angle, pointInsidePolygonQ]
angle[{p1_, p2_}] :=
Module[{c1, c2, a}, {c1, c2} = #.{1, I} & /@ {p1, p2};
a = Arg[c2/c1];
Chop[Pi - a] == 0 && Throw[True];
a
]
pointInsidePolygonQ[point_, polygon_] /; MemberQ[polygon, point] =
True;
pointInsidePolygonQ[point_, polygon_] :=

Catch[0 !=
Chop@Total[
angle /@
Partition[# - point & /@ Append[polygon, First@polygon], 2, 1]]]

polypoints = {{-1.856, 3.293}, {1.257,
2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};

pointInsidePolygonQ[(
2 polypoints[[1]] + 8 polypoints[[2]])/10, polypoints]

True

This is slower for the test data, but it will be FASTER if many test
points are on the polygon.

Bobby

>> is there a way to test whether a point is inside a polygon? ie.
>>
>> PointInsidePolygonQ[point_,polygon_] -> True or False
>>
>> i'm trying to do something like ...
>>
>> ListContourPlot[table,RegionFunction-

>> >CountryData["Canada","Polygon"]]


>>
>> to create what would be called a "clipping mask" in photoshop.
>>
>> cheers,
>> Mitch
>>
>

--
DrMaj...@longhorns.com

DrMajorBob

unread,
Feb 12, 2009, 6:41:43 AM2/12/09
to
> 1) None of these solutions will work for a country that straddles the
> Greenwich meridian or a continent that straddles a pole. Antarctica does
> both.

Oops, I meant the meridian 180 degrees from Greenwich.

Bobby

On Wed, 11 Feb 2009 17:38:29 -0600, DrMajorBob <btr...@austin.rr.com>
wrote:

mark mcclure

unread,
Feb 13, 2009, 3:44:45 AM2/13/09
to
On Feb 9, 5:31 am, Mitch Murphy <mi...@lemma.ca> wrote:
> is there a way to test whether a point is inside apolygon? ie.

> PointInsidePolygonQ[point_,polygon_] -> True or False

Of course, we can write our own program to do just about
anything we want. If you want to know whether there is
a built in function "InPolygon", the answer is no. If
you want to if there are freely available packages that
include this functionality, the answer is yes:
CompG:
http://library.wolfram.com/infocenter/MathSource/5712/
and IMTEK:
http://www.imtek.de/simulation//mathematica/IMSweb/

Curiously, the functionality is clearly "in Mathematica",
it's just not exposed. You can see that a point in
polygon is implemented by playing with the following:
Graphics[Tooltip[Polygon[
{{0, 0}, {2, 0}, {2, 1}, {1, 1}, {1, 2}, {0, 2}}],
"in polygon"]]

In fact, there must be a large computational geometry
library contained in the Mathematica kernel to support
the extensive graphics routines. Hopefully, that
functionality will be exposed in the not too distant
future.

Mark McClure

Andrzej Kozlowski

unread,
Feb 13, 2009, 3:40:57 AM2/13/09
to
I know it's only a matter of "taste" but to me using Module without
any local variables is like eating rice with knife and fork. Can be
done but why?

Andrzej Kozlowski

On 12 Feb 2009, at 11:38, Sjoerd C. de Vries wrote:

> Since timing is often important for functions like this one, I've
> compared the various algorithms posted here using David's two tests.
> The results are below:
>
> {
> {"Daniel", 0.828, 1.219},
> {"Daniel/Sjoerd", 0.688, 1.078},
> {"Sedgewick/Sjoerd", 0.688, 0.954},
> {"David1", 41.766, 81.89},
> {"David2", 1.578, 2.531},
> {"Godkin/Frank", 1.485, 2.359}
> }
>
> In the "Daniel/Sjoerd" method I have used Daniel's code and taken out
> every intermediate step. This saves an additional 20%. The code is now
> completely gibberish, but it still works ;-)
>
> PointInsidePolygonQ[pt_, poly_] :=
> Module[{},
> Equal @@
> Sign@MapThread[
> Dot, {pt - # & /@
> poly, {{0, 1}, {-1, 0}}.# & /@ (Subtract @@ # & /@
> Partition[Append[poly, poly[[1]]], 2, 1])}]
> ]
>
> Cheers --Sjoerd
>
> On Feb 9, 12:31 pm, Mitch Murphy <mi...@lemma.ca> wrote:

>> is there a way to test whether a point is inside a polygon? ie.
>>
>> PointInsidePolygonQ[point_,polygon_] -> True or False
>>
>> i'm trying to do something like ...
>>

>> ListContourPlot[table,RegionFunction->CountryData["Canada=

Daniel Lichtblau

unread,
Feb 13, 2009, 3:45:29 AM2/13/09
to
On Feb 9, 4:31 am, Mitch Murphy <mi...@lemma.ca> wrote:
> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> i'm trying to do something like ...
>
> ListContourPlot[table,RegionFunction->CountryData["Canada=

","Polygon"]]
>
> to create what would be called a "clipping mask" in photoshop.
>
> cheers,
> Mitch

The application seems to involve many such tests rather than a one-off
sort of query, so it is probably acceptable to spend time in
preprocessing, if it allows us to avoid O(n) work per query (n being
the number of segments). I'll show an implementation that bins the
segments according to x coordinates. Specifically, a bin will contain
a segment if either starting or terminating x coordinate of the
segment corresponds to x values in the bin, or else the entire bin's
range is between the starting and terminating values. Note that a
given segment might be in multiple bins. So long as we bin "sensibly",
all bins should have but a smallish fraction of the total number of
segments.

The preprocessing is a bit on the slow side because I do nothing
fancy. I guess one could use Sort and some smarts to make it faster
(which might be a good idea, if the plan is to do this for many
different objects, e.g. all countries).

Anyway, we create the bins of segments and also keep track of min and
max x and y values (these we use as a cheap way of ruling out points
not in an obvious bounding rectangle).

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[
lo = minx + (j - 1)*len - eps;
hi = minx + j*len + eps;
Select[flatsegments,
Module[{xlo, xhi}, {xlo, xhi} = Sort[{#[[1, 1]], #[[2, 1]]}];
lo <= xlo <= hi || lo <= xhi <= hi || (xlo <= lo && xhi >=
= hi)
] &],
{j, nbins}];
{{minx, maxx}, {miny, maxy}, segmentbins}
]

We preprocess for Canada.

In[125]:= canpoly = First[CountryData["Canada", "Polygon"]];

In[256]:= nbins = 100;
Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} =
polyToSegmentList[canpoly, nbins];]

Out[257]= {47.0648, Null}

A bit slow, but it's only done once. And it covers all 25 pieces, not
just the main body.

To test a point that is not outside the bounding rectangle, we will
drop a vertical from it to some point below the minimal y value, then
count intersections between that vertical and all segments in the bin
corresponding to the point's x coordinate.

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

countIntersections[segs_, {x_, yhi_, ylo_}] := Module[
{tally = 0, seg, yval, xlo, xhi, y1, y2},
Do[
seg = segs[[j]];
{{xlo, y1}, {xhi, y2}} = seg;
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
]

Here are a few examples.

In[260]:= Timing[
pointInPolygon[{-86, 60}, segmentbins, xmin, xmax, ymin, ymax]]
Out[260]= {0.005, False}

In[261]:= Timing[
pointInPolygon[{-70, 55}, segmentbins, xmin, xmax, ymin, ymax]]
Out[261]= {0.003999, True}

In[262]:= Timing[
pointInPolygon[{-137, 64}, segmentbins, xmin, xmax, ymin, ymax]]
Out[262]= {0., True}

In[263]:= Timing[
pointInPolygon[{-122, 61}, segmentbins, xmin, xmax, ymin, ymax]]
Out[263]= {0.001, True}

In[266]:= Timing[
pointInPolygon[{-73, 70}, segmentbins, xmin, xmax, ymin, ymax]]
Out[266]= {0.002999, True}

Daniel Lichtblau
Wolfram Research

Daniel Lichtblau

unread,
Feb 14, 2009, 3:10:45 AM2/14/09
to
On Feb 13, 2:45 am, Daniel Lichtblau <d...@wolfram.com> wrote:
> On Feb 9, 4:31 am, Mitch Murphy <mi...@lemma.ca> wrote:
>
> > is there a way to test whether a point is inside a polygon? ie.
>
> > PointInsidePolygonQ[point_,polygon_] -> True or False
>
> > i'm trying to do something like ...
>
> > ListContourPlot[table,RegionFunction->CountryData["Cana=

da=
> ","Polygon"]]
>
> > to create what would be called a "clipping mask" in photoshop.
>
> > cheers,
> > Mitch
>
> The application seems to involve many such tests rather than a one-off
> sort of query, so it is probably acceptable to spend time in
> preprocessing, if it allows us to avoid O(n) work per query (n being
> the number of segments). I'll show an implementation that bins the
> segments according to x coordinates. Specifically, a bin will contain
> a segment if either starting or terminating x coordinate of the
> segment corresponds to x values in the bin, or else the entire bin's
> range is between the starting and terminating values. Note that a
> given segment might be in multiple bins. So long as we bin "sensibly",
> all bins should have but a smallish fraction of the total number of
> segments.
>
> The preprocessing is a bit on the slow side because I do nothing
> fancy. I guess one could use Sort and some smarts to make it faster
> (which might be a good idea, if the plan is to do this for many
> different objects, e.g. all countries).
>
> Anyway, we create the bins of segments and also keep track of min and
> max x and y values (these we use as a cheap way of ruling out points
> not in an obvious bounding rectangle).
> [...]

I missed a couple of useful optimizations. We can use Compile both in
preprocessing and in the predicate query itself.

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

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)

] &]]];

With this we can preprocess the polygons for Canada, using 1000 bins,
in around a minute on my machine.

In[346]:= canpoly = First[CountryData["Canada", "Polygon"]];

In[347]:= nbins = 1000;


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

Out[337]= {55.3256, Null}

To repeat from the last note, there are almost certainly smarter ways
to do the preprocessing, so as to make it faster. For most
applications I can think of that would be more trouble than it is
worth, so it's not something I have attempted.

For the predicate evaluation we can do:

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

]];

With this we can now process around 10,000 points in a second. I
selected the region so that most would be inside; this was so that we
would not gain speed due to a high percentage failing the basic
rectangle test.

In[352]:= n = 10000;
pts = Transpose[{RandomReal[{-115, -55}, {n}],
RandomReal[{45, 75}, {n}]}];

In[354]:= Timing[
inout = Map[pointInPolygon[#, segmentbins, xmin, xmax, ymin, ymax] &,
pts];]
Out[354]= {1.04284, Null}

In[355]:= Take[inout, 20]
Out[355]= {True, True, True, True, False, True, True, True, True, \
True, False, True, True, False, False, False, False, True, True, True}

I should make a few remarks about the speed/complexity. If we split
the n polygon segments into m bins, for m somewhat smaller than n,
then for "good" geographies we expect something like O(n/m). Figure
most segments appear in no more than two bins, most appear in only one
bin, and no bin has more than, say, three times as many segments as
any other.

To achieve better algorithmic complexity I believe one would need to
use something more complicated, along the lines of a triangulation.
With such a data structure, efficiently implemented, it might be only O
(log(n)) to see whether a new point falls into a triangle, and, if so,
whether that triangle is inside or outside the boundary. I doubt such
a tactic would be applicable for a polygonal set of the size in
question here.

If the goal is, say, to color all pixels inside Canada differently
than those outside (in effect, to query every pixel in the bounding
rectangle), then it would be faster simply to find all boundary points
for a given vertical. This transforms from a problem of query on
points to one of splitting segments (a query on lines, so to speak).
Much more efficient, I would guess.

Daniel Lichtblau
Wolfram Research

Sjoerd C. de Vries

unread,
Feb 14, 2009, 3:14:36 AM2/14/09
to
You're right. It's a left-over that I forgot and Daniel's code could
have been reduced even more.

> >> is there a way to test whether a point is inside a polygon? ie.
>
> >> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> >> i'm trying to do something like ...
>

> >> ListContourPlot[table,RegionFunction->CountryData["Cana=
da=

0 new messages