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
- References:
- testing if a point is inside a polygon
- From: Mitch Murphy <mitch@lemma.ca>
- testing if a point is inside a polygon