计算点p与每个多边形顶点之间的有向角和。如果总倾斜角是360度，那么这个点在里面。如果总数为0，则点在外面。

我更喜欢这种方法，因为它更健壮，对数值精度的依赖更小。

计算交集数量的均匀性的方法是有限的，因为你可以在计算交集数量的过程中“击中”一个顶点。

编辑:顺便说一下，这种方法适用于凹凸多边形。

编辑:我最近在维基百科上找到了一篇关于这个话题的完整文章。

为了完整性，这里是nirg提供的算法的lua实现，由Mecki讨论:

function pnpoly(area, test)
    local inside = false
    local tx, ty = table.unpack(test)
    local j = #area
    for i=1, #area do
        local vxi, vyi = table.unpack(area[i])
        local vxj, vyj = table.unpack(area[j])
        if (vyi > ty) ~= (vyj > ty)
        and tx < (vxj - vxi)*(ty - vyi)/(vyj - vyi) + vxi
        then
            inside = not inside
        end
        j = i
    end
    return inside
end

变量区域是一个点的表，这些点依次存储为2D表。例子:

> A = {{2, 1}, {1, 2}, {15, 3}, {3, 4}, {5, 3}, {4, 1.5}}
> T = {2, 1.1}
> pnpoly(A, T)
true

GitHub Gist的链接。

当我还是Michael Stonebraker手下的一名研究员时，我做了一些关于这方面的工作——你知道，就是那位提出了Ingres、PostgreSQL等的教授。

我们意识到最快的方法是首先做一个边界框，因为它非常快。如果它在边界框之外，它就在外面。否则，你就得做更辛苦的工作……

如果你想要一个伟大的算法，看看开源项目PostgreSQL的源代码的地理工作…

我想指出的是，我们从来没有深入了解过左撇子和右撇子(也可以表达为“内”和“外”的问题……

更新

BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.

此外，我的工作使用单独的函数进行“在线”测试。

.．.因为有人问:我们发现当垂直的数量超过某个数字时，边界盒测试是最好的——如果有必要，在做更长的测试之前做一个非常快速的测试……边界框是通过简单地将最大的x，最小的x，最大的y和最小的y放在一起，组成一个框的四个点来创建的……

另一个提示是:我们在网格空间中进行了所有更复杂的“调光”计算，都是在平面上的正点上进行的，然后重新投影到“真实”的经度/纬度上，从而避免了在经度180线交叉时和处理极地时可能出现的环绕错误。工作好了!

令人惊讶的是之前没有人提出这个问题，但是对于需要数据库的实用主义者来说:MongoDB对Geo查询提供了出色的支持，包括这个查询。

你需要的是:

db.neighborhoods。findOne({geometry: {$geoIntersects: {$geometry: { type: "Point"，坐标:["经度"，"纬度"]}}} })

communities是存储一个或多个标准GeoJson格式多边形的集合。如果查询返回null，则表示不相交，否则为。

这里有详细的记录: https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/

在330个不规则多边形网格中，超过6000个点分类的性能不到一分钟，没有任何优化，包括用各自的多边形更新文档的时间。

真的很喜欢Nirg发布的解决方案，由bobobobo编辑。我只是让它javascript友好，更容易读懂我的使用:

function insidePoly(poly, pointx, pointy) {
    var i, j;
    var inside = false;
    for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
        if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
    }
    return inside;
}

在大多数情况下，这是一个比其他算法都快的算法。

它又新又雅致。我们花费O(n * log(n))时间构建一个表，允许我们在O(log(n) + k)时间内测试多边形中的点。

与光线跟踪或角度不同，使用扫描光束表可以更快地对同一多边形进行多次检查。我们必须预先构建一个扫描束活动边表，这是大多数代码正在做的事情。

We calculate the scanbeam and the active edges for that position in the y-direction. We make a list of points sorted by their y-component and we iterate through this list, for two events. Start-Y and End-Y, we track the active edges as we process the list. We process the events in order and for each scanbeam we record the y-value of the event and the active edges at each event (events being start-y and end-y) but we only record these when our event-y is different than last time (so everything at the event point is processed before we mark it in our table).

我们得到我们的表格:

[] p6p5、p6p7 p6p5, p6p7, p2p3, p2p1 p6p7, p5p4, p2p3, p3p1 p7p8, p5p4, p2p3, p2p1 p7p8, p5p4, p3p4, p2p1 p7p8 p2p1、 p7p8、p1p0 p8p0、p1p0 []

在构建该表之后，实际执行工作的代码只有几行。

注意:这里的代码使用复数值作为点。所以。real是。x。imag是。y。

def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

我们对事件进行二进制搜索，以找到特定值的actives_at_y。对于在y点的所有活动，我们计算在我们点的特定y点的x段值。每次x截距大于点的x分量时加1。然后对总数乘以2。(这是偶数-奇数填充规则，你可以很容易地适应任何其他填充规则)。

完整的代码:

from bisect import bisect

def build_edge_list(polygon):
    edge_list = []
    for i in range(1, len(polygon)):
        if (polygon[i].imag, polygon[i].real) < (polygon[i - 1].imag, polygon[i - 1].real):
            edge_list.append((polygon[i], i))
            edge_list.append((polygon[i - 1], ~i))
        else:
            edge_list.append((polygon[i], ~i))
            edge_list.append((polygon[i - 1], i))

    def sort_key(e):
        return e[0].imag, e[0].real, ~e[1]

    edge_list.sort(key=sort_key)
    return edge_list


def build_scanbeam(edge_list):
    actives_table = []
    events = []
    y = -float("inf")
    actives = []
    for pt, index in edge_list:
        if y != pt.imag:
            actives_table.append(list(actives))
            events.append(y)
        if index >= 0:
            actives.append(index)
        else:
            actives.remove(~index)
        y = pt.imag
    return actives_table, events

def point_in_polygon(polygon, point):
    def x_intercept(e, y):
        pt0 = polygon[e-1]
        pt1 = polygon[e]
        if pt1.real - pt0.real == 0:
            return pt0.real
        m = (pt1.imag - pt0.imag) / (pt1.real - pt0.real)
        b = pt0.imag - (m * pt0.real)
        return (y - b) / m

    edge_list = build_edge_list(polygon)
    actives_table, events = build_scanbeam(edge_list)
    try:
        if len(point):
            return [point_in_scantable(actives_table, events, x_intercept, p) for p in point]
    except TypeError:
        return point_in_scantable(actives_table, events, x_intercept, point)

def point_in_scantable(actives_table, events, xi, point):
    beam = bisect(events, point.imag) - 1  # Binary search in sorted array.
    actives_at_y = actives_table[beam]
    total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
    return bool(total % 2)

如果忽略，则扫描表的构建时间为O(n*log(n))。我们实际上是在O(log(n) + k)时间内查到的。其中n是多边形中段数的大小，k是该多边形中典型的活动边数。其他的光线追踪方法实际上需要O(n)时间。每次我们检查一个点，它迭代整个多边形。所以即使有这个明显的次优实现，它也轻而易举地打败了其他所有的。

There's a few performance tricks that could be done, for example, we can lower the time complexity to O(log(n) + log(k)) time. To do this we would implement Bentley-Ottmann into the sweep line, and rather than processing the intersections as different events, we split the lines at the intersections. We then also sort the active edges by their x-intercepts. We then know that no intersections occur during a scanbeam and since we sorted our segments (taking care to order them correctly within the scanbeam even if they start at the same initial point (you need to look at the slopes, or just compare midpoints of the segments). We then have a sorted intersection-less actives lists scanbeam table which means we can binary search into active edge list as well. Though that sounds like a lot of work for a value of k which is going to be typically 2 or maybe 4.

此外，由于这基本上变成了一个查找表和一些x截距的最小计算，它更能用GPU完成。你不再需要在多边形上循环了。所以你可以用numpy这样的东西来大量计算这些点，这样你就可以一次做所有的计算，从而提高性能。