• 07/16/99 05:58 PM
    Sign in to follow this  
    Followers 0

    Point in a poly

    Math and Physics

    Myopic Rhino
    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
    0


    Sign in to follow this  
    Followers 0


    User Feedback

    Create an account or sign in to leave a review

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

    Create an account

    Sign up for a new account in our community. It's easy!


    Register a new account

    Sign in

    Already have an account? Sign in here.


    Sign In Now

    There are no reviews to display.