Alternative tests for finding whether a point is inside a 2D polygon as standard tests fail

Question

I wrote the following code to find whether a point is inside a polygon or not -- the algorithm using crossing number and winding number tests. However, I find that these tests fail when my polygon is: (10,10),(10,20),(20,20) and (20,10) and the point which I want to find whether it is in the polygon is: (20,15). My points here represent the (latitude, longitude) of locations.

(20,15) lies on the boundary, therefore should be inside the polygon

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using namespace std;

typedef struct {int x, y;} Pt;

inline int
isLeft( Pt P0, Pt P1, Pt P2 )
{
    return ( (P1.x - P0.x) * (P2.y - P0.y)
            - (P2.x -  P0.x) * (P1.y - P0.y) );
}
int crossingNumTest( Pt P, Pt* V, int n )
{
    int    cn = 0;    

    for (int i=0; i<n; i++) {    
       if (((V[i].y <= P.y) && (V[i+1].y > P.y)) 
        || ((V[i].y > P.y) && (V[i+1].y <=  P.y))) { 
            float vt = (float)(P.y  - V[i].y) / (V[i+1].y - V[i].y);
            if (P.x <  V[i].x + vt * (V[i+1].x - V[i].x)) 
                 ++cn;   
        }
    }
    return (cn&1);    
}
int windingNumTest( Pt P, Pt* V, int n )
{
    int    wn = 0;    

    for (int i=0; i<n; i++) {   
        if (V[i].y <= P.y) {    
            if (V[i+1].y  > P.y)  
                 if (isLeft( V[i], V[i+1], P) > 0)  
                     ++wn;           
        }
        else {                       
            if (V[i+1].y  <= P.y)    
                 if (isLeft( V[i], V[i+1], P) < 0)  
                     --wn;          
        }
    }
    return wn;
}
//===================================================================
int main(int argc, char** argv) {

    Pt arr[11];
    arr[0].x=10;
    arr[0].y=10;
    arr[1].x=10;
    arr[1].y=20;
    arr[2].x=20;
    arr[2].y=20;
    arr[3].x=20;
    arr[3].y=10;

    Pt test;
    test.x=20; test.y=15;

    cout<<"1.polygon inside="<<crossingNumTest(test,arr,4)<<"\n";
    cout<<"2.polygon inside="<<windingNumTest(test,arr,4)<<"\n";

    return (EXIT_SUCCESS);
}

If you want the points _on_ the boundary, you can just add something like `0.01` to each point in your polygon. — Mads Marquart, Nov 20 '15 at 15:04
what is your input - exact numbers or floating point numbers? — maxim1000, Nov 20 '15 at 19:37
the posted code is mixing C headers with C++ headers, So, at best this is C++ code. Please remove the 'c' tag — user3629249, Nov 21 '15 at 19:04
In the calculations, allow for a 0 difference between a line and a point, so the point can be considered within the poly, if the difference is 0. However, note that floating numbers should not be compared to exact (==) values, however, relative values (< and >) are acceptable. — user3629249, Nov 21 '15 at 19:07
amongst other problems in the code, this (and some other lines in the psoted code) line: `if (((V[i].y <= P.y) && (V[i+1].y > P.y))` will, on the last iteration of the loop, access beyond the bounds of the `V[]` array. This is undefined behaviour and can lead to a seg fault event — user3629249, Nov 21 '15 at 19:21
the code will fail if any point is defined in other than ascending magnitude order, (example). if arr[0] values were greater than arr[1] values. This is because the sub functions do not poperly handle the cases where the calculated values are negative — user3629249, Nov 21 '15 at 19:27
@MadsMarquart, adding .01 to the boundary markers would not work as it would just (effectively) move the whole polygon up/right by the .01 offset. — user3629249, Nov 21 '15 at 19:37
@user3629249 Sorry, but i think you understand my point, to make the square a little bit bigger. — Mads Marquart, Nov 22 '15 at 21:33

score 1 · Answer 1 · edited May 23 '17 at 12:15

1

This is most likely a floating point accuracy issue.

The point you are checking is exactly on the line of the polygon, so it's quite possible that floating point is not able to represent that correctly and shows it as being outside when it is not.

You can either accept this behavior or add a small error tolerance on your checking. See What is the most effective way for float and double comparison? for a bigger picture. Effectively you need to identify the tolerance you will accept and expand your polygon a tiny amount when checking.

edited May 23 '17 at 12:15

Community

1
1

answered Nov 20 '15 at 15:03

Tim B

40,716
16
83
128

@JannatArora See expanded answer. – Tim B Nov 20 '15 at 15:07

score 1 · Answer 2 · edited May 23 '17 at 11:44

No need for float and incur potential rounding issues.
Stay with exact (likely faster) integer math. Reform test to mathematically equivalent.
(Note: no divide by 0 possible with this new approach as no division is done)

// float vt = (float)(P.y  - V[i].y) / (V[i+1].y - V[i].y);
// if (P.x <  V[i].x + vt * (V[i+1].x - V[i].x)) ++cn;   

//  Note: Use long long math if needed
int dy0 = P.y  - V[i].y;
int dy1 = V[i+1].y - V[i].y;
int dx0 = P.x  - V[i].x;
int dx1 = (V[i+1].x - V[i].x;
if (dy1*dx0 < dy0*dx1) ++cn;

[Edit] @user3629249 comment points out an array access violation. Suggest:

int crossingNumTest( Pt P, Pt* V, int n ) {
    // for (int i=0; i<n; i++) {    
    for (int i=1; i<n; i++) {    
       // if (((V[i].y <= P.y) && (V[i+1].y > P.y)) 
       if (((V[i-1].y <= P.y) && (V[i].y > P.y)) 
    ...
}

Alternative tests for finding whether a point is inside a 2D polygon as standard tests fail

2 Answers2