Search code examples
c++algorithmpoint-in-polygon

What would make the point in polygon algorithm only work correctly when run in both forward and reverse?

The most popular answer I found for "point in polygon algorithm for c++" seems to be this C algorithm: https://wrf.ecse.rpi.edu/Research/Short_Notes/pnpoly.html

However, when I implemented this in my program, it was not able to work correctly for certain shapes. Like when I draw a circle but then one side is a straight line at slight angle. After hours of guessing and trying other algorithms, I came up with a solution where I run the algorithm once and then I run the algorithm again but with the polygon points reversed. This fixed the problem for all shapes, but I don't understand why. How would I go about eliminating this extra loop? Is it possible the input data needs to have the points sorted or processed somehow to work correctly?

My application of the algorithm is that I am letting the user draw any polygon shape on a google map with the javascript api, then the points of the shape are sent to the server where C++ then runs the algorithm to exclude locations from the search results. Because the user is drawing the polygon, it could be drawn clockwise or counterclockwise and be any complex or even overlapping shape so I can't really force a predictable drawing order or a simple polygon in my application.

My data structure is a little different, but I think the logic is the same as the original algorithm. I also have used >= in 2 of the comparisons on purpose which I think covers some cases where the point is on the edge of the polygon, but I think the rest is the same.

The arrays are storing the points like x,y,x,y. My array always has the same point at the beginning and the end of the array. I'm also using int on purpose, because I don't need double/float in my use case. I am displaying a small portion of the world and have converted degrees to feet to use int.

bool pointInPolygon = false;
int lastPointIndex=search->polygonSize-2;
int lastX=search->polygon[lastPointIndex];
int lastY=search->polygon[lastPointIndex+1];
for (int i=0; i < search->polygonSize; i+=2) {
    int x=search->polygon[i];
    int y=search->polygon[i+1];
    if ( ((y>=listing->longitude) != (lastY>=listing->longitude)) && 
(listing->latitude <= (lastX-x) * (listing->longitude-y) / (lastY-y) + x) ){
        pointInPolygon = !pointInPolygon;
    }
    lastX=x;
    lastY=y;
    lastPointIndex=i;
}
if(!pointInPolygon){
    return false;
}
pointInPolygon=false;
lastPointIndex=search->polygonSize-2;
lastX=search->polygon[lastPointIndex];
lastY=search->polygon[lastPointIndex+1];
for (int i=0; i < search->polygonSize; i+=2) {
    int x=search->polygonReverse[i];
    int y=search->polygonReverse[i+1];
    if ( ((y>=listing->longitude) != (lastY>=listing->longitude)) && 
(listing->latitude <= (lastX-x) * (listing->longitude-y) / (lastY-y) + x) ){
        pointInPolygon = !pointInPolygon;
    }
    lastX=x;
    lastY=y;
    lastPointIndex=i;
}
if(!pointInPolygon){
    return false;
}

This is an example of the points with the problem: 8245533,28415352,8236498,28392000,8220158,28373794,8205897,28364691,8195806,28361525,8176659,28360734,8170041,28362317,8161681,28366275,8152969,28374190,8143906,28390417,8140769,28403082,8137979,28435932,8137979,28451368,8139374,28459679,8148438,28482239,8153318,28489759,8159590,28495695,8165861,28499258,8174918,28500445,8245533,28415352

Here is a point that is not filtered correctly when the algorithm is only run once. This problem does not occur when it is run twice like my code example. 8230098,28456924

This image shows the same polygon and point on the map. Note that I am offsetting my longitude coordinates 180 degrees to have positive numbers.

Update: It looks like aschepler spotted my ints are overflowing in this case of a long line in the polygon and I needed to use double. I will continue to use int but convert to double in a polygon search to correct this. I will mark them as the answer if they choose to post it as the answer. Thank you!

I thought it would be clever to use int instead of double for coordinates because I only needed to search the USA, and I have combined 8 of my search criteria into AVX2 SIMD instructions that are stored as int in a 256 bit data type in C++ to make a very fast in memory search system. After the initial rectangle filtering of the coordinates, I now have to cast to double to do the polygon search on a smaller number of records. It is still pretty efficient. I'll have to remember that multiplying my coordinates needs to be in the higher precisions if I need this again someday. Searching for just a rectangle only needed > or < operations so there was plenty of precision for that.

map example

Solution

  • @aschepler answered my problem in the comments above. The int was overflowing when the points were a long straight line. I had to use higher precision casting for that multiply operation like double or long long to fix it. Thanks!