Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!mnetor!seismo!ut-sally!husc6!sri-unix!amdahl!dlb!auspyr!sci!kenm
From: kenm@sci.UUCP (Ken McElvain)
Newsgroups: comp.graphics
Subject: Re: Point inside of polygon revisited...
Message-ID: <6169@sci.UUCP>
Date: Thu, 18-Jun-87 18:37:27 EDT
Article-I.D.: sci.6169
Posted: Thu Jun 18 18:37:27 1987
Date-Received: Mon, 22-Jun-87 01:48:17 EDT
References: <948@elrond.CalComp.COM>
Organization: Silicon Compilers Systems Corp. San Jose, Ca
Lines: 108
Keywords: ray tracing, graphics, polygon intersection
Summary: winding numbers

In article <948@elrond.CalComp.COM>, amamaral@elrond.CalComp.COM (Alan Amaral) writes:
> 
> A while ago someone asked for a method for determining whether or not
> a point was inside or outside of a polygon.  I have the replies, and
> am in the process of implementing just such a function, but have come
> up against what I view as a significant problem.
> 
> To recap, the theorum put forth was that given a bounded area on a plane
> and a point on the plane, one could determine whether the point was inside
> the area, or outside the area by determining how many edges an arbitrary
> ray starting at the given point crossed. (and a partridge in a pear tree...)
> If the ray crossed an even number of edges, then the point is outside of the
> area, and it's inside if the number of crossings is odd.
> 
> This all works well and good, EXCEPT for one particular case. What if the
> ray happens to cross 2 edges at exactly the same point (i.e. at a corner)?
> 
> Has anyone solved this particular problem, or do you just ignore it?
> I have tried all sorts of schemes to work around this particular problem,
> but have always been able to come up with some kind of pathological case
> where a concave polygon will fail.
> 
> HELP!


Winding numbers work a bit better.  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.

typedef struct point_ {
	float 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;
	float a,b;
	int wind;	/* current winding number */

	wind = 0;
	lastpt = polypts[cnt-1];
	oldquad = whichquad(lastpt,pt);
	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
 */
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;
	}
}

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.

I tried to mail this but it bounced, maybe others will be interested.

Ken McElvain
decwrl!sci!kenm