Re: testing if a point is inside a polygon
- To: mathgroup at smc.vnet.net
- Subject: [mg96306] Re: testing if a point is inside a polygon
- From: "Sjoerd C. de Vries" <sjoerd.c.devries at gmail.com>
- Date: Wed, 11 Feb 2009 05:20:57 -0500 (EST)
- References: <23667197.1234176058052.JavaMail.root@m02> <gmrm5r$9nm$1@smc.vnet.net>
Nice one David, but it's about 50 times slower than using Sedgewick's
algorithm...
Cheers -- Sjoerd
On Feb 10, 12:50 pm, "David Park" <djmp... at comcast.net> wrote:
> Mitch,
>
> Here is my attempt. The idea is to move around the polygon and add up the
> angles made with the point. If the point is outside this will be zero. Th=
e
> steps are basically:
>
> 1) Check if the point is one of the vertices of the polygon.
> 2) Add the first polygon point to the end of the list of points.
> 3) Subtract the test point from each of the polygon vertices.
> 4) Partition the points into pairs.
> 5) Rotate each pair so the first one lies along the x axis. (To prevent
> crossing the branch line)
> 6) Find the ArcTan of the second point. This will have a sign.
> 7) Sum the angles and take the absolute value.
> 8) Test if this is zero or some multiple of 2\[Pi].
>
> PointInsidePolygonQ::usage =
> "PointInsidePolygonQ[point,polygon] will return True if the point \
> is on the boundary or inside the polygon and False otherwise.";
>
> PointInsidePolygonQ[point_, polygon_] :=
> Module[{work},
> work = If[Head[polygon] === Polygon, First[polygon], polygon]=
;
> If[ \[Not] FreeQ[work, point], Return[True]];
> work = # - point & /@ Join[work, {First[work]}];
> work = Partition[work, 2, 1];
> work = RotationTransform[{First[#], {1, 0}}]@# & /@ work;
> work = Abs@Total[ArcTan @@ Last[#] & /@ work] // Chop // Rationaliz=
e;
> TrueQ[work != 0]
> ]
>
> Here is a graphical test for a simple polygon:
>
> testpoints = testpoints = RandomReal[{-9, 9}, {100, 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[4],
> {If[PointInsidePolygonQ[#, polypoints], Black, Red], Point[#]} & /@
> testpoints},
> PlotRange -> 10,
> Frame -> True]
>
> Here is a test for a more complex polygon.
>
> testpoints = testpoints = RandomReal[{-9, 9}, {100, 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[4],
> {If[PointInsidePolygonQ[#, Polygon[polypoints]], Black, Red],
> Point[#]} & /@ testpoints},
> PlotRange -> 10,
> Frame -> True]
>
> There might be simpler methods.
>
> David Park
> djmp... at comcast.nethttp://home.comcast.net/~djmpark/
>
> From: Mitch Murphy [mailto:mi... at lemma.ca]
>
> 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