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