Re: Re: testing if a point is inside a polygon
- To: mathgroup at smc.vnet.net
- Subject: [mg96312] Re: [mg96233] Re: [mg96189] testing if a point is inside a polygon
- From: "David Park" <djmpark at comcast.net>
- Date: Wed, 11 Feb 2009 05:22:02 -0500 (EST)
- References: <200902091032.FAA12225@smc.vnet.net> <25096359.1234267282082.JavaMail.root@m02>
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.
David Park
djmpark at comcast.net
http://home.comcast.net/~djmpark/
From: Frank Scherbaum [mailto:Frank.Scherbaum at 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-
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];
];
];
- References:
- testing if a point is inside a polygon
- From: Mitch Murphy <mitch@lemma.ca>
- testing if a point is inside a polygon