RE: Point inside a plygon? and PointInsidePolyhedronQ
- To: mathgroup at smc.vnet.net
- Subject: [mg25321] RE: [mg25239] Point inside a plygon? and [mg24992]PointInsidePolyhedronQ
- From: "David Park" <djmp at earthlink.net>
- Date: Sat, 23 Sep 2000 03:35:46 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
Dear Andrzej, MathGroup and All,
First in regard to the so-called counter-example that Will Self provided: My
routine does calculate the correct answer. His example was one in which the
test line passed exactly through a vertex of the polygon. His assumption was
that the algorithm would count twice (or perhaps zero times). However the
algorithm uses a <= at the beginning of an edge, and a < at the end of an
edge and counts only once. I don't know if we could occasionally get into
trouble here with approximate numbers. But so far it works. The reliability
could be vastly increased by using a second outside test point to confirm.
It is not too difficult to modify the algorithm to calculate winding number
as suggested by Andrzej. The essential idea is to check each edge for an
intersection. If there is one, rotate the test line vector counterclockwise
by 90 degrees, and take the Sign of the dot product with the edge. Sum this
for all the edges to obtain the winding number.
Here is the code and a testing plot.
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]]
JOperator::usage = "JOperator[{x,y}] returns {-y,x}, rotating the \
vector counterclockwise by 90°.";
JOperator[{x_, y_}] := {-y, x}
WindingIncrement::usage = "WindingIncrement[testline, edge] calculates \
the increment in winding number for a testline, {testpoint, \
outsidepoint}, crossing the edge, {start, end}, of a polygon.";
WindingIncrement[testline_, edge_] :=
If[SegmentIntersection[testline, edge],
Sign[JOperator[Subtract @@ Reverse[testline]] .
Subtract @@ Reverse[edge]], 0]
PointWindingNumber::usage = "PointWindingNumber[point, \
polygonvertices] will calculate the winding number of a point in a \
polygon which may be self intersecting.";
PointWindingNumber[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 = {(3*xmax - xmin)/2, (ymin + ymax)/2};
edges = Partition[Join[polypoints, {First[polypoints]}], 2, 1];
count = (WindingIncrement[{pt, outpoint}, #1] & ) /@ edges;
Plus @@ count]
This tests with random points and makes a plot.
Module[{pt = {Random[Real, {-5, 5}], Random[Real, {-5, 5}]},
poly = {{-4.21607, -3.58889}, {1.08015, -3.5192},
{4.49482, -1.01046}, {4.49482, 3.17076}, {0.940775, 5.47044},
{-1.49827, 3.10108}, {-0.940776, -0.522654}, {1.98608, 1.28921},
{-1.98608, 6.09762}, {-5.81887, 4.28576}, {-5.81887, -1.01046}},
linepts}, linepts = Join[poly, {First[poly]}];
Show[Graphics[{Line[linepts], AbsolutePointSize[4], Point /@ poly,
Red, Point[pt], Black, Text[PointWindingNumber[pt, poly],
{-9, 9}, {-1, 0}]}], AspectRatio -> Automatic,
PlotRange -> {{-10, 10}, {-10, 10}}, Background -> Linen]];
This does the same thing with the polygon orientation reversed.
Module[{pt = {Random[Real, {-5, 5}], Random[Real, {-5, 5}]},
poly = Reverse[{{-4.21607, -3.58889}, {1.08015, -3.5192},
{4.49482, -1.01046}, {4.49482, 3.17076}, {0.940775, 5.47044},
{-1.49827, 3.10108}, {-0.940776, -0.522654},
{1.98608, 1.28921}, {-1.98608, 6.09762}, {-5.81887, 4.28576},
{-5.81887, -1.01046}}], linepts},
linepts = Join[poly, {First[poly]}];
Show[Graphics[{Line[linepts], AbsolutePointSize[4], Point /@ poly,
Red, Point[pt], Black, Text[PointWindingNumber[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/
> -----Original Message-----
> From: Andrzej Kozlowski [mailto:andrzej at tuins.ac.jp]
To: mathgroup at smc.vnet.net
>
> Well, I have at also found a bit of time to actually read David's message
> and must admit that the idea is quite correct (even if the implementation
> may not imperfect), and I was quite wrong to doubt it. It works fine for
> polygons without self-intersections. It can be improved to deal also with
> those with self-intersections (when a line connecting the point
> to one which
> lies outside passes through the self-intersection point). This seems
> actually to give a fast geometric method of computing winding numbers and
> hence can be used to write a package that would greatly improve
> Mathematica's abilities to compute certain integrals in the complex plane.
>
>
>
> Andrzej
>
>