       Re: testing if a point is inside a polygon

• To: mathgroup at smc.vnet.net
• Subject: [mg96233] Re: [mg96189] testing if a point is inside a polygon
• From: Frank Scherbaum <Frank.Scherbaum at geo.uni-potsdam.de>
• Date: Tue, 10 Feb 2009 05:48:09 -0500 (EST)
• References: <200902091032.FAA12225@smc.vnet.net>

```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: "mapping" functions over lists, again!
• Next by Date: Re: Re: FourierTransform
• Previous by thread: testing if a point is inside a polygon
• Next by thread: Re: testing if a point is inside a polygon