Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • 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