点在多边形内算法为何仅正反双向运行时结果才准确？

问题背景 #

我搜索"C++ 点在多边形内算法"得到的最热门答案是C语言实现的pnpoly算法。我在程序中实现该算法后，部分形状的判定结果存在错误，比如某条边为微斜直线的类圆形。
我尝试了数小时、测试了其他算法后找到了临时解决方案：先正常运行一次算法，再将多边形顶点顺序反转后再运行一次，该方案解决了所有形状的判定问题，但我不清楚背后的原因。我想要去掉这额外的遍历步骤，不确定是否需要对输入顶点进行排序或其他预处理才能让算法正常工作。

应用场景 #

允许用户通过JavaScript API在谷歌地图上绘制任意多边形，形状顶点会发送到服务端，由C++运行该算法过滤不符合要求的搜索结果。由于多边形是用户手动绘制的，可能是顺时针、逆时针、任意复杂形状甚至自相交形状，无法强制规定绘制顺序或要求只能绘制简单多边形。

实现细节 #

• 算法逻辑和原pnpoly算法一致，特意在两处比较中使用>=运算符，用于覆盖点落在多边形边上的场景
• 顶点数组按x,y,x,y的格式存储，数组首尾的顶点始终相同
• 坐标使用int类型存储：仅需要展示世界的一小部分区域，已经将经纬度转换为英尺单位，不需要使用double/float类型

问题代码 #

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

复现示例 #

• 出现问题的多边形顶点：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
• 仅运行一次算法时判定错误的点：8230098,28456924

下图展示了地图上的对应多边形和点，经度坐标已偏移180度以获得正数：

根因与优化方案 #

问题根源为处理多边形长边时int类型发生了运算溢出，不需要保留正反两次遍历的逻辑，也不需要修改顶点预处理规则，只需要保留int存储坐标的方案，在多边形搜索时将坐标转换为double类型运算即可修复问题。
原本用int替代double存储坐标是为了配合C++的AVX2 SIMD指令，以256位int类型存储实现极快的内存搜索，坐标初始矩形过滤阶段仅需要大小比较，int精度完全足够，只有过滤后的少量数据需要转换为double进行多边形搜索，整体效率不受影响。

内容的提问来源于stack exchange，提问作者Bruce Kirkpatrick