RE: Point inside a plygon? and PointInsidePolyhedronQ
- To: mathgroup at smc.vnet.net
- Subject: [mg25350] RE: [mg25239] Point inside a plygon? and [mg24992] PointInsidePolyhedronQ
- From: "David Park" <djmp at earthlink.net>
- Date: Sat, 23 Sep 2000 03:36:25 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Dear MathGroup and All,
Here is another algorithm for determining if a point is within a polygon. It
works for any polygon defined by a sequence of points, and is approximately
100 times faster than the complex contour integral method. The method was
suggested to me by Richard I. Pelletier.
The idea is to draw a line from the point in question to a point that is
known to be outside the polygon. We then find out how many times that line
intersects the edge segments of the polygon. If it is an odd number of
times, the point is within the polygon, otherwise it is outside. We can
easily generate a point outside the polygon by finding the maximum and
minimum coordinates of the vertices.
The following routine tests if two line segments intersect. The first
segment would be from the test point to the outside point. The second
segment would be an edge of the polygon. The expressions were obtained by
parametrizing the two line segments so that their parameters run from 0 to 1
over the segment, and solving for the intersection point. d is a denominator
and is 0 for parallel lines. The test for the edges is designed to avoid
double counting intersections at vertices.
SegmentIntersection::usage = "SegmentIntersection[{p1,p2},{p3,p4}] \
determine if the line segment {p1,p2} intersects the line segment {p3,p4}."\
;
SegmentIntersection[{{x1_, y1_}, {x2_, y2_}}, {{x3_, y3_}, {x4_, y4_}}] :=
With[{d = x4*(y1 - y2) + x3*(-y1 + y2) + (x1 - x2)*(y3 - y4)},
d != 0 && 0 <= (x4*(y1 - y3) + x1*(y3 - y4) + x3*(-y1 + y4))/d <= 1 &&
Inequality[0, LessEqual, (x3*(-y1 + y2) + x2*(y1 - y3) + x1*(-y2 + y3))/
d, Less, 1]]
This is the routine that tests if a point is within a polygon.
IntersectionPointInPolygon::usage =
"IntersectionPointInPolygon[point, polygonvertices] will determine if a
\
point is within a polygon defined by its vertices. Its method is to define a
\
line segment between the point and an outside point, and count how many
times \
it intersects the edge segments of the polygon. An odd count puts the point
\
within the polygon.";
IntersectionPointInPolygon[pt : {_, _}, polypoints : {{_, _} ..}] :=
Module[{xmin, xmax, ymin, ymax, xpts, ypts, outpoint, edges, count},
{xpts, ypts} = Transpose[polypoints];
{xmax, ymax} = Max /@ {xpts, ypts};
{xmin, ymin} = Min /@ {xpts, ypts};
outpoint = {(3xmax - xmin)/2, (ymin + ymax)/2};
edges = Partition[Join[polypoints, {First[polypoints]}], 2, 1];
count = Count[SegmentIntersection[{pt, outpoint}, #] & /@ edges, True];
OddQ[count]
]
Needs["Graphics`Colors`"]
This tests the routine with a moderately complex polygon and random points.
Module[
{pt = {Random[Real, {-5, 5}], Random[Real, {-5, 5}]},
poly = {{-3.72826, -3.17076}, {-0.801402, -1.84671}, {1.70733, \
-3.79795}, {5.26138, -1.21952}, {2.4042, 0.801401}, {4.00701,
4.14638}, {0.243904, 3.24045}, {-3.65857, 5.05231}, {-5.88856,
1.56796}, {-4.0767, -0.174218}},
linepts},
linepts = Join[poly, {First[poly]}];
Show[Graphics[
{Line[linepts],
AbsolutePointSize[4], Point /@ poly,
Red, Point[pt],
Black,
Text[IntersectionPointInPolygon[pt, poly], {-9, 9}, {-1, 0}]}],
AspectRatio -> Automatic, PlotRange -> {{-10, 10}, {-10, 10}},
Background -> Linen]];
David Park
djmp at earthlink.net
http://home.earthlink.net/~djmp/