Re: Point inside a plygon? and PointInsidePolyhedronQ
- To: mathgroup at smc.vnet.net
- Subject: [mg25354] Re: [mg25239] Point inside a plygon? and [mg24992] PointInsidePolyhedronQ
- From: Andrzej Kozlowski <andrzej at tuins.ac.jp>
- Date: Sat, 23 Sep 2000 03:36:38 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
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/