       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,
{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,
{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,
{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,
{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-
>
> 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];
inside += isicr;
];
(* CHECK SEGMENT FROM LAST VERTEX TO FIRST VERTEX *)
{px,py} = poly[[n]];
{pxx,pyy} = poly[];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
If[isicr == 4, Return];
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 (* 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;, (* no crossing *)
Return[-2]; (* downward crossing *)
];
,
(* CASE   Y1 < 0 < Y2 *)
If[x1 * y2 >= y1 * x2, Return;, (* no crossing *)
Return; (* 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; (* no crossing *)
,
(* UPWARD OR DOWNWARD? *)
If[y1 > 0.,
Return[-1]; (* Downward half crossing *)
,
Return;  (* Upward half crossing *)
];
];
,
(* HERE Y1==0   CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return;
,
(* UPWARD OR DOWNWARD? *)
If[y2 > 0,
Return;  (* Upward half crossing *)
,
Return[-1]; (* Downward half crossing *)
];
];
];

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

```

• Prev by Date: Re: Mathematicas simplifications
• Next by Date: Re: Log[x]//TraditionalForm
• Previous by thread: Re: testing if a point is inside a polygon
• Next by thread: Re: Re: Re: testing if a point is