Re: Point inside a plygon? and PointInsidePolyhedronQ
- To: mathgroup at smc.vnet.net
- Subject: [mg25355] Re: [mg25239] Point inside a plygon? and [mg24992] PointInsidePolyhedronQ
- From: Andrzej Kozlowski <andrzej at tuins.ac.jp>
- Date: Sat, 23 Sep 2000 03:36:40 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
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
on 00.9.22 1:26 PM, Andrzej Kozlowski at andrzej at tuins.ac.jp wrote:
> Sorry about this last message . It was the result of quickly reading two
> different messages, one from David Park, and one from Will Self and confusing
> them in my mind. David's method does not use triangulation's, so what I wrote
> below does not apply. However, I am also sceptical of any methods of deciding
> if a point is inside or outside a polygon as simple as the one David is
> describing. If it worked it would make it easy to use Mathematica to integrate
> (very fast) complex functions with isolated singularities over arbitrary
> closed contours (I mean using Integrate, not NIntegrate) in the complex plane,
> (since Mathematica can already find residues). I suspect this would be too
> good to be true.
>
>
> Andrzej
>
> on 00.9.22 8:06 AM, Andrzej Kozlowski at andrzej at tuins.ac.jp wrote:
>
>> I am really to busy at this time to do so myself but I hope other people will
>> test this. I must admit however that I am sceptical. The reason is this. Take
>> almost any book on combinatorial or piecewise linear topology e.g. Rourke and
>> Sanderson, "Piecewise-Linear Topology" and look at the proof of the theorem
>> (Theorem 2.11 of R&S) , which says that any compact polyhedron can be
>> triangulated. Admittedly this is a lot more general than what you are trying
>> to do. Still, you will see that the proof is non-constructive and much more
>> complicated than the one you are proposing. Of course I am not saying that it
>> is impossible to find simple and constructive proofs of classical results
>> that
>> earlier had complicated non-constructive ones, but it does not happen
>> terribly
>> often. On the other hand finding mistakes in other people proofs is generally
>> time consuming and not very rewarding, unless it concerns your own work. Of
>> course I wish you luck and I hope an interesting Mathematica package will
>> arise out of this but I suggest you write out the mathematical part of your
>> argument and get some specialists on combinatorial topology and geometry to
>> have a look at it.
>>
>>
>> Andrzej
>>
>>
>> on 00.9.22 5:32 AM, David Park at djmp at earthlink.net wrote:
>>
>>> 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/
>>>
--
Andrzej Kozlowski
Toyama International University
JAPAN
http://platon.c.u-tokyo.ac.jp/andrzej/
http://sigma.tuins.ac.jp/