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