Help with this algorithm

Question

I have an algorithm that can find if a point is inside a polygon.

 int CGlEngineFunctions::PointInPoly(int npts, float *xp, float *yp, float x, float y)
 {
     int i, j, c = 0;
     for (i = 0, j = npts-1; i < npts; j = i++) {
         if ((((yp[i] <= y) && (y < yp[j])) ||
             ((yp[j] <= y) && (y < yp[i]))) &&
             (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
             c = !c;
     }
     return c;
 }

My only issue with it is it assumes an odd winding rule. What I mean by this is that if the polygon is self intersecting, certain parts that it would considered to be 'empty' will return as false. What I'd need in even if it self intersects, anything inside the polygon will return true.

Thanks

Answer 1

Beware: this answer is wrong. I have no time to fix it right now, but see the comments.

This casts a ray from the point to infinity, and checks for intersections with each of the polygon's edges. Each time an intersection is found, the flag c is toggled:

c = !c;

So an even number of intersections means an even number of toggles, so c will be 0 at the end. An odd number of intersections means an odd number of toggles, so c will be 1.

What you want instead is to set the c flag if any intersection occurs:

c = 1;

And for good measure, you can then eliminate c entirely, and terminate early:

 int CGlEngineFunctions::PointInPoly(int npts, float *xp, float *yp, float x, float y)
 {
     int i, j;
     for (i = 0, j = npts-1; i < npts; j = i++) {
         if ((((yp[i] <= y) && (y < yp[j])) ||
             ((yp[j] <= y) && (y < yp[i]))) &&
             (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
             return 1;
     }
     return 0;
 }

Answer 2

To translate your original algorithm to english: You're determining if the number of polygon segments to the right of your point are even or odd. If it's even (including zero) your point is outside, if it's odd your point is inside. This means if there are two segments to the right and also two segments to the left, the point is not considered inside the polygon.

What you need to do is change the algorithm so that it checks for segments on both sides; if there's a segment on both sides of the point, then the point is within the polygon.

 int CGlEngineFunctions::PointInPoly(int npts, float *xp, float *yp, float x, float y) 
 { 
     int i, j;
     bool hasSegmentLeft = false;
     bool hasSegmentRight = false;
     for (i = 0, j = npts-1; i < npts; j = i++) { 
         if ((((yp[i] <= y) && (y < yp[j])) || 
             ((yp[j] <= y) && (y < yp[i]))))
         {
             if (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])
             {
                 hasSegmentRight = true;
                 if (hasSegmentLeft)  // short circuit early return
                     return true;
             }
             else
             {
                 hasSegmentLeft = true;
                 if (hasSegmentRight)  // short circuit early return
                     return true;
             }
     } 
     return hasSegmentLeft && hasSegmentRight; 
 }

P.S. that for statement construct is a very clever way of dealing with a circular list that wraps around to the beginning; I'd never seen it before.