104 views

Skip to first unread message

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

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

];

];

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

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/

Feb 10, 2009, 5:49:09 AM2/10/09

to

In article <gmp0mt$btn$1...@smc.vnet.net>, Mitch Murphy <mi...@lemma.ca>

wrote:

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/

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

Feb 10, 2009, 5:52:55 AM2/10/09

to

Here is a function based on the code in Robert Sedgewick's Algorithms

book:

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=

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

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/

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>

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-

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:

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

Feb 12, 2009, 6:37:57 AM2/12/09

to

I decided a specific "mean" would do, and simplifications gave me the

following:

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:

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:

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=

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

>>

>

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.

> 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:

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

> 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

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?

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=

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

>

","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

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

>

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

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.

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=

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu