MathGroup Archive 2007

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

Search the Archive

Re: Point inside polygon

  • To: mathgroup at smc.vnet.net
  • Subject: [mg73211] Re: Point inside polygon
  • From: Jens-Peer Kuska <kuska at informatik.uni-leipzig.de>
  • Date: Wed, 7 Feb 2007 06:14:46 -0500 (EST)
  • Organization: Uni Leipzig
  • References: <eq9btd$muu$1@smc.vnet.net>
  • Reply-to: kuska at informatik.uni-leipzig.de

Hi,

and the simple method to sum the angles of the triangles
build up from the point and all cyclic consecutive vertices
does not help. For all inner points it is 2Pi and for all outer
points zero or keeping the numerical errors in mind much less
than 2Pi

Regards
   Jens

Goyder, Hugh wrote:
> I have been reading mathgroup items for determining if a point is inside or
>  outside a general polygon. David Park gave a method which he acknowledged
> still has some flaws. With the help of an internet search (Kreylos and "point inside a polygon") I found an improvement on his method which works better than I expected.
> 
> The basic method is to draw a line from the test point to infinity and count the number of edge crossings. If the number of edge crossings is odd then the test point is inside the polygon. This basic method needs extending to cover the cases where the test line goes through a vertex, where a vertex is tangent to the test line and where an edge lies on the test line. An additional problem, which I take separately, is when the test point lines on the boundary.
> 
> The method can be simplified, from that given by David Park,  by always drawing a horizontal test line from the test point to infinity in the positive x direction. Edges are only counted when they cross the test line. A crossing occurs if one edge end-point is above or on the test line while the other edge end-point is below the test line. If an edge crosses, the crossing point must be to the right of the test point. The careful definition eliminates edges that lie on the test line because they have both end-points on the test line. Tangent vertices cause no difficulty because those that touch the test line from below are counted twice while those that touch from above are not counted at all. Both cases add no odd counts to the total count.  The passage of the test line through a vertex correctly gives just one count. The method can be coded in only a few lines.
> 
> The code, PointInPolygon, is given below together with some test cases. I have introduced a scale factor so that the test may be repeated with big or small numbers. They all seem to work. Am I lucky or does Mathematica's inequality testing save me? Try changing the scale factors to 10^100 or 10^-100 to try and break the code.
> 
> For the case of test points on the boundary I do an initial check and exit with False if the point lies on the boundary. I determine if a point is on the boundary by working out the area of the triangle formed by the test point and the two points defining a boundary edge. For these calculations I have converted to Integers using Rationalize. The method now works unless very big or small numbers (ie outside 10^14 to 10^-14) are used which I guess is down to Rationalize. The codes which include a boundary check are PointInsidePolygon and PointOnLine. I give examples below.
> 
> Are there examples which break this method?
> Can this be improved? How can it be made faster?
> 
> Hugh Goyder
> Cranfield University
> Tel: +44 (0) 1793 785122
> Mob: +44 (0) 7804 252770
> 
> ClearAll[PointInPolygon];
> PointInPolygon::usage == "PointInPolygon[p0,poly] gives True if the point p0 \
> is inside the polygon poly. The polygon is given as a list of lists where \
> each sublist corresponds to a list of external or internal vertices. If p0 is \
> on the boundary the behaviour is unpredictable.";
> PointInPolygon[{x0_, y0_}, poly_] :== Module[{edges, incount},
>    edges == Flatten[(Partition[Join[#1, {First[#1]}], 2, 1] & ) /@ poly, 1];
>     incount == 0;
>  Do[{{x1, y1}, {x2, y2}} == edges[[n]];
>       If[(y1 >== y0 && y2 < y0) || (y2 >== y0 && y1 < y0),
>        If[x1 + ((y1 - y0)*(x2 - x1))/(y1 - y2) > x0, ++incount]],
>  {n, Length[edges]}];
>  OddQ[incount]];
> 
> 
> size == 10;
> poly == N[{{{0, 0}, {1, -1}, {2, 0}, {2.5, -1.1}, {3, 0}, {4, 0}, {4, 1.5},
>        {3, 1.5}, {3, 1}, {2.5, 0}, {2, 1}, {2, 2}, {0.5, 2}, {0.5, 1}},
>       {{0.75, 1}, {1.75, 1}, {1.75, 1.75}, {0.75, 1.75}}}*size];
> inpoints == N[{{0.62, 1.5}, {1.88, 1.5}, {0.62, 1}, {1.88, 1}, {1, 0},
>       {2.25, 0}, {2.75, 0}, {3, 0.5}}*size];
> outpoints == N[{{1, 1.5}, {-0.5, 0}, {4.5, 0}, {0, 1}, {2.5, 1}}*size];
> midpoints == N[((1/2)*(RotateLeft[#1] + #1) & )[poly[[1]], 1]];
> 
> Show[Graphics[{(Line[#1 /. {a_, b__} :> {a, b, a}] & ) /@ poly,
>       {Green, PointSize[0.02], (Point[#1] & ) /@ inpoints},
>       {Red, PointSize[0.02], (Point[#1] & ) /@ outpoints},
>       {Blue, PointSize[0.02], (Point[#1] & ) /@ midpoints}},
>      AspectRatio -> Automatic, Frame -> True]];
> 
> (PointInPolygon[#1, poly] & ) /@ inpoints
> 
> (PointInPolygon[#1, poly] & ) /@ outpoints
> 
> Comment["The method fails if the vertices are tested"];
> (PointInPolygon[#1, poly] & ) /@ Partition[Flatten[poly], 2]
> 
> Comment["The method fails if points on the boundary are tested"];
> (PointInPolygon[#1, poly] & ) /@ midpoints
> 
> ClearAll[PointOnLine];
> PointOnLine::usage == "PointOnLine[p0,{p1,p2}] determines if the point p0 is \
> on the line between points p1 and p2. The method works by determining if the \
> area of the triangle formed by the three points is zero and if the point p0 \
> lies between p1 and p2 when projected onto the x-axis (or y-axis if the line \
> is vertical). The values are converted to integers.";
> PointOnLine[p0_, {p1_, p2_}] :== Module[{x0, y0, x1, y1, x2, y2},
>  {x0, y0} == Rationalize[p0, 0];
>  {x1, y1} == Rationalize[p1, 0];
>  {x2, y2} == Rationalize[p2, 0];
>   x2*(y0 - y1) + x0*(y1 - y2) + x1*(-y0 + y2) ==== 0 &&
>   If[x1 !== x2,
>  0 <== (x0 - x1)/(x2 - x1) <== 1,
>  0 <== (y0 - y1)/(y2 - y1) <== 1]];
> 
> ClearAll[PointInsidePolygon];
> PointInsidePolygon::usage == "PointInPolygon[p0,poly] gives True if the point \
> p0 is inside the polygon poly. The polygon is given as a list of lists where \
> each sublist corresponds to a list of external or internal vertices.";
> PointInsidePolygon[{x0_, y0_}, poly_] :==
>   Module[{edges, incount, x1, y1, x2, y2},
>    edges == Flatten[(Partition[Join[#1, {First[#1]}], 2, 1] & ) /@ poly, 1] ;
>     incount == 0;
> Catch[Do[{{x1, y1}, {x2, y2}} == edges[[n]];
>   If[PointOnLine[{x0, y0},{{x1, y1}, {x2, y2}}], Throw[False]];
>   If[(y1 >== y0 && y2 < y0) || (y2 >== y0 && y1 < y0),
>   If[x1 + ((y1 - y0)/(y1 - y2))*(x2 - x1) > x0,
>     ++incount]],
>   {n, Length[edges]}];
>    OddQ[incount]]];
> 
> (PointInsidePolygon[#1, poly] & ) /@ inpoints
> 
> (PointInsidePolygon[#1, poly] & ) /@ outpoints
> 
> Comment["The method works if the vertices are tested"];
> (PointInsidePolygon[#1, poly] & ) /@ Partition[Flatten[poly], 2]
> 
> Comment["The method works if points on the boundary are tested"];
> (PointInsidePolygon[#1, poly] & ) /@ midpoints
> 
> 
> 
> 


  • Prev by Date: Re: String to number
  • Next by Date: Re: String to number
  • Previous by thread: Re: Point inside polygon
  • Next by thread: BarChart and Ticks problem.