planting@speedy.cs.wisc.edu (W. Harry Plantinga) (02/27/88)
In article <7626@pur-ee.UUCP> beatyr@pur-ee.UUCP (Robert Beaty) writes: >In article <20533@amdcad.AMD.COM> msprick@amdcad.UUCP (Rick Dorris) writes: >>In article <971@ut-emx.UUCP> jezebel@ut-emx.UUCP (Jim @ MAC) writes: >>>Have a nice one here: >>>Have a boundary defined on the screen. Boundary is composed of points >>>joined by lines... Now some random points are generated and I want to check >>>if a given point is within or outside the existing boundary.. Any algorithm for >> <rick suggests the infinite ray/intersection counting algorithm> > >I have seen this type of algorithm before and the one I thought up >in an afternoon is FAR surperior and vastly simpler to code up. > ><robert beaty suggests an algorithm that measures sums the angles of >the lines connecting the poitns and the test point> Haven't I seen this discussion of point-in-polygon algorithms about 15 times before? :-) Robert, your algorithm is simpler only in the kind of problems it handles. It will only work for convex polygons, whereas the other works for any polygons. It's not even much simpler; measuring angles and determining whether line segments intersect are of comparable difficulty. Maybe you should have said that althogh your algorithm is far inferior, it's not any easier to code?
kenm@sci.UUCP (Ken McElvain) (02/28/88)
In article <5305@spool.cs.wisc.edu>, planting@speedy.cs.wisc.edu (W. Harry Plantinga) writes: > In article <7626@pur-ee.UUCP> beatyr@pur-ee.UUCP (Robert Beaty) writes: > >In article <20533@amdcad.AMD.COM> msprick@amdcad.UUCP (Rick Dorris) writes: > >>In article <971@ut-emx.UUCP> jezebel@ut-emx.UUCP (Jim @ MAC) writes: > >>>Have a nice one here: > >>>Have a boundary defined on the screen. Boundary is composed of points > >>>joined by lines... Now some random points are generated and I want to check > >>>if a given point is within or outside the existing boundary.. Any algorithm for > >> <rick suggests the infinite ray/intersection counting algorithm> > > > >I have seen this type of algorithm before and the one I thought up > >in an afternoon is FAR surperior and vastly simpler to code up. > > > ><robert beaty suggests an algorithm that measures sums the angles of > >the lines connecting the poitns and the test point> > > Haven't I seen this discussion of point-in-polygon algorithms about 15 > times before? :-) > > Robert, your algorithm is simpler only in the kind of problems it > handles. It will only work for convex polygons, whereas the other > works for any polygons. It's not even much simpler; measuring angles > and determining whether line segments intersect are of comparable > difficulty. > > Maybe you should have said that althogh your algorithm is far > inferior, it's not any easier to code? Robert's suggestion is a good one. The sum of the angles about the test point is known as the winding number. For non intersecting polygons if the winding number is non-zero then the test point is inside the polygon. It works just fine for convex and concave polygon's. Intersecting polygon's give reasonable results too. A figure 8 will give a negitive winding number for a test point in one of the loops and a positive winding number for the other loop, with all points outside having a winding number of 0. Other advantages of the winding number calculation are that it does not suffer from the confusion of the infinite ray algorithm when the ray crosses a vertex or is colinear with an edge. Here is my version of a point in poly routine using a quadrant granularity winding number. The basic idea is to total the angle changes for a wiper arm with its origin at the point and whos end follows the polygon points. If the angle change is 0 then you are outside, otherwise you are in some sense inside. It is not necessary to compute exact angles, resolution to a quadrant is sufficient, greatly simplifying the calculations. I pulled this out of some other code and hopefully didn't break it in doing so. There is some ambiguity in this version as to whether a point lying on the polygon is inside or out. This can be fairly easily detected though, so you can do what you want in that case. ----------------------------------------------------------------- /* * Quadrants: * 1 | 0 * ----- * 2 | 3 */ typedef struct { int x,y; } point; pointinpoly(pt,cnt,polypts) point pt; /* point to check */ int cnt; /* number of points in poly */ point *polypts; /* array of points, */ /* last edge from polypts[cnt-1] to polypts[0] */ { int oldquad,newquad; point thispt,lastpt; int a,b; int wind; /* current winding number */ wind = 0; lastpt = polypts[cnt-1]; oldquad = whichquad(lastpt,pt); /* get starting angle */ for(i=0;i<cnt;i++) { /* for each point in the polygon */ thispt = polypts[i]; newquad = whichquad(thispt,pt); if(oldquad!=newquad) { /* adjust wind */ /* * use mod 4 comparsions to see if we have * advanced or backed up one quadrant */ if(((oldquad+1)&3)==newquad) wind++; else if((((newquad+1)&3)==oldquad) wind--; else { /* * upper left to lower right, or * upper right to lower left. Determine * direction of winding by intersection * with x==0. */ a = lastpt.y - thispt.y; a *= (pt.x - lastpt.x); b = lastpt.x - thispt.x; a += lastpt.y * b; b *= pt.y; if(a > b) wind += 2; else wind -= 2; } } lastpt = thispt; } return(wind); /* non zero means point in poly */ } /* * Figure out which quadrent pt is in with respect to orig */ int whichquad(pt,orig) point pt; point orig; { if(pt.x < orig.x) { if(pt.y < orig.y) quad = 2; else quad = 1; } else { if(pt.y < orig.y) quad = 3; else quad = 0; } } Ken McElvain decwrl!sci!kenm