benson@blake.acs.washington.edu (Dan Benson) (08/08/90)
I know this has been asked before but I missed the answers. I want to know the best way to determine whether a point is inside of a polygon where the polygon is represented as a set of ordered floating point numbers (its boundary). I'm actually using 3D points but the polygons are defined on 2D planes. Also, the polygons can be any shape, concave, convex. If anyone has some suggestions or real answers I would appreciate it. And if there is more interest I will post a summary. Thanks in advance, -- | Dan Benson benson@ee.washington.edu | | Dept. of Elec. Engr., FT-10 | | University of Washington (206) 685-7567 | | Seattle, WA 98195 |
metz@iam.unibe.ch (Igor Metz) (08/08/90)
In article <5990@milton.u.washington.edu> benson@blake.acs.washington.edu (Dan Benson) writes: >I know this has been asked before but I missed the answers. I want to >know the best way to determine whether a point is inside of a polygon >where the polygon is represented as a set of ordered floating point >numbers (its boundary). I'm actually using 3D points but the polygons >are defined on 2D planes. Also, the polygons can be any shape, concave, >convex. get THE book: F.P. Preparata/ M.I. Shamos: Computational Geometry: An Introduction. 2nd ed. New York: Springer, 1988. p. 41ff -- Igor Metz Institut fuer Informatik und angew. Mathematik, Universitaet Bern, Switzerland. domainNet: metz@iam.unibe.ch Phone: (0041) 31 65 49 90 ARPA: metz%iam.unibe.ch@relay.cs.net Fax: (0041) 31 65 39 65
phil@hpsmdca.HP.COM (Philip Walden) (08/08/90)
I use a simple routine that computes the sum of arc angles between the
test point and each vertex on the polygon. If the sum == 0, its outside
If the summ == 2pi its inside.
I have attached a brief explanation in pseudo-fortran.
In the actual algorithm I use, the polygon can actually be a set
of polygons (members) such that "holes" can be created. The key to it all
is the fortran intrinsic ATAN2(rise,run) which returns the 4 quandrant
angle (-pi to pi) of the specified vector. Don't use ATAN(rise/run) as
it returns only 2 quadrant angles (-pi/2 to pi/2).
Yes it works for anything.
Brief explanation:
*************************
Is a point inside a polygon?
Use sum of arc angles:
Given a point (x,y) and a polygon SUM(i=1,n) (xi,yi), total arc angle is
sum = 0
for i = 1,n ! assumes (x1,y1) = (xn,yn)
do
dx1 = (x - xi)
dy1 = (y - yi)
dx2 = (xi+1 - x)
dy2 = (yi+1 - y)
ds1 = SQRT(dx1**2 + dy1**2)
ds2 = SQRT(dx2**2 + dy2**2)
sini = (dx1*dy2 - dx2*dy1)/ds1*ds2 !from cross product
cosi = (dx1*dx2 + dy1*dy2)/ds1*ds2 !from dot product
sum = sum + atan2(sini,cosi)
end do
If |sum| > 0, point is inside
Note: you can skip the computation of ds1, ds2 and the divisions by ds1*ds2,
as they cancel out in the atan2 function. I just put them in to
clarify the algorithm.rusmin@twinsun.com (Rusmin Dirgantoro) (08/22/90)
In article <5990@milton.u.washington.edu> benson@blake.acs.washington.edu (Dan Benson) writes: >I know this has been asked before but I missed the answers. I want to >know the best way to determine whether a point is inside of a polygon >where the polygon is represented as a set of ordered floating point >numbers (its boundary). I'm actually using 3D points but the polygons >are defined on 2D planes. Also, the polygons can be any shape, concave, >convex. First of all, I missed the answers of previous postings as well. So.. this idea I'm posting might not be new to you guys. I was working on something like this some time ago for my Geometric Modeling class. The idea I was using was adapted from some good articles of some journals. I don't remember the exact name and date of the journals by heart. The idea is.. 0. Assume your polygon is PG and your point is PT. 1. Create a coordinate frame CF, where the XY plane is the plane of PG. We'll be working with this coordinate system from now on. 2. Transform PT's coordinate into CF coordinate, say we get (CF)PT. 3. If (CF)PTz != 0, PT IS OUTSIDE PG; DONE. If (CF)PTz == 0, PT may be inside PG; continue step 4. 4. Transform PG's boundary (the set of ordered points) into CF coordinates. The points should then have Z == 0, of course. 5. Build segment equations S1..Sn for each segment on the boundary. S1..Sn are line eqs with endpoints. They don't go from/to infinity. S1..Sn are based on CF coordinate system, so they are simply 2D-eqs since we don't care about Z anymore (Z is always 0). 6. Build a random (yes, random) ray equation R. R starts from (CF)PT to infinity. R is on the XY plane of CF. Again, we don't care about Z anymore. 7. Find intersections between R and S1..Sn. Say the number of intersections is NI. 8. If NI is even or zero, PT IS OUTSIDE PG; DONE. If NI is odd, PT IS INSIDE PG; DONE. Please note here that there's a very special case, that is.. Somehow R intersects any of the ordered points (the corners of PG), so step 8 won't necessary give the right answer. One of the solutions for this special case is.. Add step 7 with: When intersecting R with S1..Sn, if R intersects Si at an Si's endpoint, back to step 6 (try another random R). This solution may look dumb to you. But.. believe me.. it works just fine. Here, I assume we are all aware of computation errors, such as round-up errors, floating point errors, etc. So appropiate floating point comparisons should be applied. Of course there are many other ideas on this problem. I'd like to hear them as well. BTW, I'm sure I still have some piece of code doing this problem. And I still have hardcopies of the articles I mentioned above. If anyone is desperate, I'm willing to search those resources in my hayestacks. Just let me know. Ciao. -- -- Rusmin Dirgantoro rusmin@twinsun.com Twin Sun Inc. 1(213)640-6885 voice El Segundo, California ---
markv@gauss.Princeton.EDU (Mark VandeWettering) (08/22/90)
In article <1990Aug22.022743.4139@twinsun.com> rusmin@twinsun.com (Rusmin Dirgantoro) writes: >In article <5990@milton.u.washington.edu> benson@blake.acs.washington.edu (Dan Benson) writes: >> [ ask for points in polygon algorithms ] > [ responds with an "almost right" algorithm ] The solution proposed by rusmin@twinsun.com was based upon the Jordan Curve theorem: any ray from inside a polygon crosses an odd number of edges on its way to infinity. He chose a random ray that began at the point in question, and counted the number of intersections. Problem: if you intersected a vertex you were kind of clueless. Solution, fire another ray. You solve these problems by simplifying everything. The ray you shoot should go to the positive X axis. Assume without loss of generality that your point is the origin. Now: if you are going to intersect a vertex, its because the y component of an edge endpoint is == 0. Well, decide whether you want to count this as positive or negative. Assume positive (I always do). It turns out you get the right answer anyway. For example origin ^ | o o counts as one intersection | v origin ^ ^ |/ o o counts as zero or two intersections origin o o counts as zero or two intersections |\ v v Voila! C'est explique! mark
rusmin@twinsun.com (Rusmin Dirgantoro) (08/23/90)
In article <2018@idunno.Princeton.EDU> markv@gauss.Princeton.EDU (Mark VandeWettering) writes: >In article <1990Aug22.022743.4139@twinsun.com> rusmin@twinsun.com (Rusmin Dirgantoro) writes: >>In article <5990@milton.u.washington.edu> benson@blake.acs.washington.edu (Dan Benson) writes: >>> [ ask for points in polygon algorithms ] >> [ responds with an "almost right" algorithm ] > >You solve these problems by simplifying everything. The ray you shoot should >go to the positive X axis. Assume without loss of generality that your >point is the origin. Now: if you are going to intersect a vertex, its because >the y component of an edge endpoint is == 0. Well, decide whether you want >to count this as positive or negative. Assume positive (I always do). It >turns out you get the right answer anyway. For example >... Mark came up with a smart solution for the "ambiguous ray-edge intersection" (AREI) problem (anyone agrees with this name?). I like his article. But his solution has a little glitch. Consider this case.. origin ^ ^ | | o o---> or <---o or o---> or <---o | | v v (1) (2) (3) (4) Are we going to count (1), (2), (3), or (4) as one or zero? Consider this polygon... +--------------------------------+ | | | | | D E F | y=0 - - o - - o----o----o - - o - - o----o-----o - - o - - A B C | | G H I | | | | +-----------+ A..I are checkpoints for origin or possible positions of the point relative to the polygon. How are we going to do it? I believe there're solutions for this kind of glitch. Are they simple, though? Do we really want to do them? Ciao.. -- Rusmin Dirgantoro rusmin@twinsun.com Twin Sun Inc. 1(213)640-6885 office El Segundo, California 1(213)837-3941 home
markv@gauss.Princeton.EDU (Mark VandeWettering) (08/23/90)
>Mark came up with a smart solution for the "ambiguous ray-edge intersection" >(AREI) problem (anyone agrees with this name?). I like his article. >But his solution has a little glitch. Hardly my own idea. It was shown to me by Eric Haines originally, but I don't think he claims credit for it either. Any takers? Is it patented by Unisys as well :-)? But, it is stump the computer jock time, he gave a nice tough diagram, and its time to show that this really will work. >Consider this case.. > > origin ^ ^ > | | > o o---> or <---o or o---> or <---o > | | > v v > (1) (2) (3) (4) > >Are we going to count (1), (2), (3), or (4) as one or zero? Apply the rule I gave you blindly, lets assume we break ties by assuming they are in above the point in question. Therfore 1 & 2 count as zero intersections, and 3 & 4 count as 1 intersection. Now, on to the case you gave: >Consider this polygon... > > +--------------------------------+ > | | > | | > | D E F | >y=0 - - o - - o----o----o - - o - - o----o-----o - - o - - > A B C | | G H I > | | > | | > +-----------+ > A, B, C, E, G, H, are all counted as zero. D & F are counted as one. Therefore the total is two, and the point is outside. If we broke ties the other way, then B and H count one, and we get the same answer. Pretty damned amazing, isn't it? I really do thank Eric for showing me this. Its one of my favorite little snippets of code. Mark
jindak@surfside.sgi.com (Chris Schoeneman) (08/24/90)
In article <2034@idunno.Princeton.EDU> markv@gauss.Princeton.EDU (Mark VandeWettering) writes: >>Mark came up with a smart solution for the "ambiguous ray-edge intersection" >>(AREI) problem (anyone agrees with this name?). I like his article. >>But his solution has a little glitch. > >Apply the rule I gave you blindly, lets assume we break ties by assuming >they are in above the point in question. Therfore 1 & 2 count as zero >intersections, and 3 & 4 count as 1 intersection. Now, on to the case >you gave: > >>Consider this polygon... >> >> +--------------------------------+ >> | | >> | | >> | D E F | >>y=0 - - o - - o----o----o - - o - - o----o-----o - - o - - >> A B C | | G H I >> | | >> | | >> +-----------+ >> >Mark Mark is right, of course: the same code works. But I feel this solution can be stated more clearly, and make acceleration techniques more obvious. After all, speed is what it's all about, isn't it? The ray interesects an edge if and only if the edge crosses the x-axis. A vertex on the x-axis is considered above it, that is: zero is considered to be positive. * If we assume 0 is positive, the answer presents itself: * count an edge if and only if sign(y1)!=sign(y2). Now we can determine a crossing very fast. Two's-compliment (integer) numbers of n bits can be XOR'd and then AND'd with 2^(n-1). Non-zero means a crossing, zero means no crossing. But in 3D floating point is the norm. What do we do? The SAME thing! Many floating point storage schemes are arranged like so: +-+----------+----------------+ MSB |S| Exponent | Mantissa | LSB +-+----------+----------------+ The S is the sign bit, and it gives the sign of the mantissa, i.e. the sign of the overall number. Like two's-compliment, zero has a zero sign bit; it is in this sense positive. If your floating point is stored this way, then you can XOR and AND, just like before. This is NOT portable, some fp's may be stored differently. Some floating points have both a positive and negative zero. (You can get a negative zero by, say, multiplying -1.0E-25 * 1.0E-25. The answer has no significant digits, but is definitely negative.) Is this a problem? No, a point will have a negative zero y value if it IS below (but very close) to the x-axis. If the other side is above the x-axis then a crossing should be, and is, detected. -Chris Schoeneman Chris Schoeneman | I was neat, clean, shaved and sober, jindak@surfside.esd.sgi.com | and I didn't care who knew it. Silicon Graphics, Inc. | -Raymond Chandler Mountain View, CA | (The Big Sleep)
corkum@csri.toronto.edu (Brent Thomas Corkum) (08/24/90)
There's one question I have concerning this algorithm that I hope someone
can answer. It seems to me that a point exactly on an edge can evaluate
either inside or outside depending on which edge your on.
ie.
+----A----+
| |
| |
+----B----+
Using the algorithm I've seen, A would evaluate inside and B outside. Is
there a quick and dirty way to handle points exactly on edges? I realize
that a seperate calculation can be done to determine whether a point
is on a line but this is expensive computationally isn't it?
Brent
corkum@ecf.toronto.edurusmin@twinsun.com (Rusmin Dirgantoro) (08/24/90)
In article <1990Aug23.000651.10581@twinsun.com> rusmin@twinsun.com (Rusmin Dirgantoro) writes: >In article <2018@idunno.Princeton.EDU> markv@gauss.Princeton.EDU (Mark VandeWettering) writes: >>In article <1990Aug22.022743.4139@twinsun.com> rusmin@twinsun.com (Rusmin Dirgantoro) writes: >>>In article <5990@milton.u.washington.edu> benson@blake.acs.washington.edu (Dan Benson) writes: >>>> [ ask for points in polygon algorithms ] >>> [ responds with an "almost right" algorithm ] >But his solution has a little glitch. >... >Consider this polygon... > > +--------------------------------+ > | | > | | > | D E F | >y=0 - - o - - o----o----o - - o - - o----o-----o - - o - - > A B C | | G H I > | | > | | > +-----------+ Ooops... This is not right. I'm sorry about that. After rethinking about it, I found that Mark's algorithm works in this case as well. I was out of my head. Just ignore my previous posting. Thanks for your responses, though. Bye. -- Rusmin Dirgantoro rusmin@twinsun.com Twin Sun Inc. 1(213)640-6885 office El Segundo, California 1(213)837-3941 home
rusmin@twinsun.com (Rusmin Dirgantoro) (08/24/90)
In article <1990Aug23.181401.2258@jarvis.csri.toronto.edu> corkum@csri.toronto.edu (Brent Thomas Corkum) writes: > >There's one question I have concerning this algorithm that I hope someone >can answer. It seems to me that a point exactly on an edge can evaluate >either inside or outside depending on which edge your on. >ie. > +----A----+ > | | > | | > +----B----+ >Using the algorithm I've seen, A would evaluate inside and B outside. Is >there a quick and dirty way to handle points exactly on edges? I realize >that a seperate calculation can be done to determine whether a point >is on a line but this is expensive computationally isn't it? I'm not sure what you mean by 'handle'. But.. if I'm allowed to add something here, I believe this would lead to a discussion of PMC (Point Membership Classification), a classic issue in Geometric Modeling. Points on boundaries (here, edges) always have two neighborhoods, inside the object (here, polygon) and outside. For this reason people give points on boundaries *ON* classification, not *IN* nor *OUT*. Ex: In your problem.. A has *IN* neighborhood *BELOW* and *OUT* neighborhood *ABOVE*. B has *IN* neighborhood *ABOVE* and *OUT* neighborhood *BELOW*. So, both A & B have *ON* classification. Since we ought to pick either above or below on the ray (re: Eric/Mark's algorithm), the algorithm doesn't explore enough information for points with *ON* classification. It was meant to decide only either *IN* or *OUT*. That's why the algorithm could give different answers for *ON* points. What do you think Eric? Mark? Sorry, I am not answering your questions.. just sharing a little old idea. G'day.. -- Rusmin Dirgantoro rusmin@twinsun.com Twin Sun Inc. 1(213)640-6885 office El Segundo, California 1(213)837-3941 home
zap@lysator.liu.se ("Zap" Andersson) (08/28/90)
Forgive me if I am but a fool, but I have created a very simple algorithm and I have never had any problems with the Ambiguos Ray Intersections. But perhaps my polygon storage methods help me? I simply swap all endpoints so each line has its 'first' end upwards and its 'last' end downwards. Then I simply EXCLUDE the 'first' point from the intersection test, but INCLUDE the 'last' i.e. point is on the line if first < point <= last .... I havn't been thinking to hard if this is a loser method but it works ok in my implementations.....??? /Zap