CODEWITHSUNDEEP
c++cgeolocationgeometrygeo
Jannat Arora
Jannat Arora

Reputation: 2989

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

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);
}

Upvotes: 1

Views: 87

Answers (2)

chux
chux

Reputation: 153602

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)) 
    ...
}

Upvotes: 1

Tim B
Tim B

Reputation: 41188

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.

Upvotes: 1

Related Questions