Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2009

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

Search the Archive

Re: Re: Re: testing if a point is

  • To: mathgroup at smc.vnet.net
  • Subject: [mg96387] Re: [mg96312] Re: [mg96233] Re: [mg96189] testing if a point is
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Thu, 12 Feb 2009 06:41:59 -0500 (EST)
  • References: <200902091032.FAA12225@smc.vnet.net>
  • Reply-to: drmajorbob at longhorns.com

> 1) None of these solutions will work for a country that straddles the  
> Greenwich meridian or a continent that straddles a pole. Antarctica does  
> both.

Oops, I meant the meridian 180 degrees from Greenwich.

Bobby

On Wed, 11 Feb 2009 17:38:29 -0600, DrMajorBob <btreat1 at austin.rr.com>  
wrote:

> Just a few caveats:
>
> 1) None of these solutions will work for a country that straddles the  
> Greenwich meridian or a continent that straddles a pole. Antarctica does  
> both.
>
> CountryData["Antarctica", "Shape"]
>
> To handle this in general, we'd need to rotate the globe so that the  
> country/continent in question doesn't have that problem... and test  
> points must be on the same hemisphere, too, to be inside a country.
>
> (This is why I was asking about RotationMatrix.)
>
> 2) The poster mentioned Canada, which is made up of disjoint polygons:
>
> Dimensions /@ First@canada;
> Length@%
> %%[[All, 1]] // Total
>
> 25
>
> 30458
>
> Twenty-five polygons, more than 30K points... and that's the BASIC  
> version. Here's the other:
>
> canada = CountryData["Canada", "FullPolygon"];
> Dimensions /@ First@canada;
> Length@%
> %%[[All, 1]] // Total
>
> 971
>
> 43262
>
> 971 polygons, 43 thousand points! For the OP's application, he may need  
> VERY fast code.
>
> 3) In the PointInsidePolygon code, I'd use MemberQ and Append, not FreeQ  
> and Join, and I NEVER, EVER use Return. (But I found no difference in  
> speed.)
>
> pointInsidePolygonQ[point_, polygon_] /; MemberQ[polygon, point] = True;
> pointInsidePolygonQ[point_, polygon_] :=
>    0 != Chop@
>     Total[angle /@
>       Partition[# - point & /@ Append[polygon, First@polygon], 2, 1]]
>
> 4) Points on edges of the polygon test as OUTSIDE:
>
> polypoints = {{-1.856, 3.293}, {1.257,
>      2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
> PointInsidePolygonQ[(
>   2 polypoints[[1]] + 8 polypoints[[2]])/10, polypoints]
>
> False
>
> One solution for this would be to stop calculating angles if one of them  
> happens to be Pi. For instance,
>
> Clear[angle, pointInsidePolygonQ]
> angle[{p1_, p2_}] :=
>   Module[{c1, c2, a}, {c1, c2} = #.{1, I} & /@ {p1, p2};
>    a = Arg[c2/c1];
>    Chop[Pi - a] == 0 && Throw[True];
>    a
>    ]
> pointInsidePolygonQ[point_, polygon_] /; MemberQ[polygon, point] =
>    True;
> pointInsidePolygonQ[point_, polygon_] :=
>
>   Catch[0 !=
>     Chop@Total[
>       angle /@
>        Partition[# - point & /@ Append[polygon, First@polygon], 2, 1]]]
>
> polypoints = {{-1.856, 3.293}, {1.257,
>      2.635}, {2.395, -0.6587}, {-1.018, -2.455}, {-3.293, -0.05988}};
> pointInsidePolygonQ[(
>   2 polypoints[[1]] + 8 polypoints[[2]])/10, polypoints]
>
> True
>
> This is slower for the test data, but it will be FASTER if many test  
> points are on the polygon.
>
> Bobby
>
> On Wed, 11 Feb 2009 04:22:02 -0600, David Park <djmpark at comcast.net>  
> wrote:
>
>> 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];
>>          ];
>> ];
>>
>>
>>
>>
>
>
>



-- 
DrMajorBob at longhorns.com


  • Prev by Date: New free introductory book on Mathematica programming, and a few
  • Next by Date: Re: LabeledListPlot
  • Previous by thread: Re: Re: testing if a point is inside a polygon
  • Next by thread: Re: Re: Re: testing if a point is