mathematics - Adding tolerance to a point in polygon test

I've been using this method which was taken from Game Coding Complete to detect whether a point is inside of a polygon. It works in almost every case, but is failing on a few edge cases, and I can't figure out the reason.

For example, given a polygon with vertices at (0,0) (0,100) and (100,100), the algorithm is returning:

• True for any point strictly inside the polygon


• False for any of the vertices


• False for (0, 50) which lies on one of the edges of the polygon


• True (?) for (50,50) which is also on one of the edges of the polygon

I'd actually like to relax the algorithm so that it returns true in all of these cases. In other words, it should return true for points that are strictly inside, for the vertices themselves, and for points on the edges of the polygon.

If possible I'd also like to give it enough tolerance so that it always tend towards "true" in face of floating point fluctuations. For example, I have another method, that given a line segment and a point, returns the closest location on the line segment to the given point.

Currently, given any point outside the polygon and one of its edges, there are cases where the result is categorized as being inside by the method above, while other points are considered outside. I'd like to give it enough tolerance so that it always returns true in this situation.

The way I've currently solved the problem is an hack, which consists of using an external library to inflate the polygon by a few pixels, and performing the tests on the inflated polygon, but I'd really like to replace this with a proper solution.

Answer

Since the algorithm already parses all edges of the polygon to see how many times a ray cast from test crosses them, I think it's reasonable to add a check to see whether test lies on the edge exactly (or within an epsilon).

In order to avoid too much additional complexity (the usual point - segment distance computation is really awful), I suggest approximating the edge to an extremely thin ellipse and see whether test is on that ellipse. This is the resulting code (only the inner loop shown):

oldPoint = polygon[points-1];
float oldSqDist = sqlength(oldPoint - test);

for (unsigned int i=0 ; i < points; i++)
{
    newPoint = polygon[i];

    float newSqDist = sqlength(newPoint - test);

    if (oldSqDist + newSqDist + 2.0f * std::sqrt(oldSqDist * newSqDist) - sqlength(newPoint - oldPoint) < EPSILON)
        return true;

    if (newPoint.x > oldPoint.x)
    {
        left = oldPoint;
        right = newPoint;
    }

    else
    {
        left = newPoint;
        right = oldPoint;
    }

    if ((newPoint.x < test.x) == (test.x <= oldPoint.x)
       && (test.y-left.y) * (right.x-left.x)
        < (right.y-left.y) * (test.x-left.x) )
    {

        inside=!inside;
    }

    oldPoint = newPoint;
    oldSqDist = newSqDist;
}

Where sqlength() is your favourite way of doing x * x + y * y. I tested the code with EPSILON = 1e-10f and got good results. In real life you should have EPSILON vary relatively to the size of your polygon.