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