• # Point in a poly

Math and Physics

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 thispt = polypts; 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

Report Article

## User Feedback

You need to be a member in order to leave a review