3

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

Improve this question
3
  • Perhaps you could just reorder the points such that the order you draw them in will no self intersect? Commented Aug 9, 2010 at 18:14
  • No, my application is a vector drawing one which allows different kinds of winding rules Commented Aug 9, 2010 at 18:24
  • 2
    The algorithm is based on the Jordon curve theorem. There are essentially 2 spaces associated with a bounded object - either in or out. Casting a ray out from the point in question will intersect the boundary either an odd or even number of times. If it is odd then the point is in the boundary, if its even then the point is outside of the boundary. Also this particular version of the algorithm assume a simple optimisation based on the fact that the ray being cast is horizontal. Commented Aug 10, 2010 at 9:58

2 Answers 2

Reset to default
5

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;
 }
Improve this answer
Sign up to request clarification or add additional context in comments.

7 Comments

And fix the design by returning true/false and change the return type to bool. Also, npts (and therefore i and j) should be unsigned, since a signed value doesn't make sense. Change it to size_t. Lastly, a more C++-type interface would be a template function that takes iterators to the first and last of the points (or even better, just a collection of points wrapped in a class). But I suppose that's a larger re-design.
Thanks for the tips GMan, I have a more templatized and C++ version of it, I just wanted to show this to make it easier to understand
Oh, from the style of coding, I just assumed that this was C... I didn't look to hard at the first line, I guess.
@Jex, I believe this algorithm fails if the point is to the left of the entire polygon.
@Mark: Better is to add/subtract 1, depending whether the line between (x,y) and infinity crosses the edge from left to right (looking from start to end) or right to left. Then the Odd rule is calculated by (c % 2), while the Zero Winding rule is (c != 0). (Almost the) same calculation for both rules!
|
1

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.

Improve this answer

1 Comment

What about a U-shaped polygon where the point is in the bowl of the U (outside the polygon, but inside its convex hull)? You'll find intersections on both sides, but the point is not inside the polygon.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.