Point in a poly


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
Discuss this article in the forums
Date this article was posted to GameDev.net: 7/16/1999
(Note that this date does not necessarily correspond to the date the article was written)
See Also:
Polygons