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

MathGroup Archive 2009

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

Search the Archive

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





  • 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