Xref: utzoo comp.graphics:1837 comp.sys.ibm.pc:12567
Path: utzoo!mnetor!uunet!seismo!sundc!pitstop!sun!decwrl!pyramid!voder!tolerant!sci!kenm
From: kenm@sci.UUCP (Ken McElvain)
Newsgroups: comp.graphics,comp.sys.ibm.pc
Subject: Re: point-in-polygon again
Message-ID: <16460@sci.UUCP>
Date: 28 Feb 88 00:21:05 GMT
References: <971@ut-emx.UUCP> <20533@amdcad.AMD.COM> <7626@pur-ee.UUCP> <5305@spool.cs.wisc.edu>
Distribution: na
Organization: Silicon Compilers Systems Corp. San Jose, Ca
Lines: 132
Keywords: READ THIS ONLY IF YOUR IQ <<125
Summary: point in poly by winding number

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