Re: Point inside polygon
- To: mathgroup at smc.vnet.net
- Subject: [mg73207] Re: Point inside polygon
- From: dh <dh at metrohm.ch>
- Date: Wed, 7 Feb 2007 05:58:41 -0500 (EST)
- References: <eq9btd$muu$1@smc.vnet.net>
Hi Goyder,
It seems to me that if you draw the test line from -infinity to
+infinity, the check on "x" becomes superfluous.
Daniel
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
>
>
>
>