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 > >

**Re: Point inside a plygon? and PointInsidePolyhedronQ**

**Re: Point inside a plygon? and PointInsidePolyhedronQ**

**Re: Point inside a plygon? and PointInsidePolyhedronQ**

**Re: Point inside a plygon? and PointInsidePolyhedronQ**