如何确定一个二维点是否在多边形内？（涉及IT技术）

这是@nirg答案（Philipp Lenssen的JavaScript版本）的Rust版本。 我提供这个答案是因为我从这个网站得到了很多帮助，我将JavaScript版本翻译成Rust作为练习，并希望能帮助一些人，最后一个原因是在我的工作中，我将把这段代码翻译成wasm以提高我的画布性能，这是一个开始。我的英语很差，请原谅。

pub struct Point {
    x: f32,
    y: f32,
}
pub fn point_is_in_poly(pt: Point, polygon: &Vec<Point>) -> bool {
    let mut is_inside = false;

    let max_x = polygon.iter().map(|pt| pt.x).reduce(f32::max).unwrap();
    let min_x = polygon.iter().map(|pt| pt.x).reduce(f32::min).unwrap();
    let max_y = polygon.iter().map(|pt| pt.y).reduce(f32::max).unwrap();
    let min_y = polygon.iter().map(|pt| pt.y).reduce(f32::min).unwrap();

    if pt.x < min_x || pt.x > max_x || pt.y < min_y || pt.y > max_y {
        return is_inside;
    }

    let len = polygon.len();
    let mut j = len - 1;

    for i in 0..len {
        let y_i_value = polygon[i].y > pt.y;
        let y_j_value = polygon[j].y > pt.y;
        let last_check = (polygon[j].x - polygon[i].x) * (pt.y - polygon[i].y)
            / (polygon[j].y - polygon[i].y)
            + polygon[i].x;
        if y_i_value != y_j_value && pt.x < last_check {
            is_inside = !is_inside;
        }
        j = i;
    }
    is_inside
}


let pt = Point {
    x: 1266.753,
    y: 97.655,
};
let polygon = vec![
    Point {
        x: 725.278,
        y: 203.586,
    },
    Point {
        x: 486.831,
        y: 441.931,
    },
    Point {
        x: 905.77,
        y: 445.241,
    },
    Point {
        x: 1026.649,
        y: 201.931,
    },
];
let pt1 = Point {
    x: 725.278,
    y: 203.586,
};
let pt2 = Point {
    x: 872.652,
    y: 321.103,
};
println!("{}", point_is_in_poly(pt, &polygon));// false
println!("{}", point_is_in_poly(pt1, &polygon)); // true
println!("{}", point_is_in_poly(pt2, &polygon));// true

`

当使用qt（Qt 4.3+）时，可以使用QPolygon的函数containsPoint

要检测多边形的命中，我们需要测试两件事：

1. 点是否在多边形区域内（可以通过射线投射算法实现）。
2. 点是否在多边形边界上（可以通过与用于检测折线（线）上的点的相同算法实现）。

在射线投射算法中处理以下特殊情况：

1. 射线与多边形的一条边重合。
2. 点位于多边形内部，且射线穿过多边形的一个顶点。
3. 点位于多边形外部，且射线刚好触及多边形的一个角。

请查看确定复杂多边形内部点的方法。该文章提供了一种简单的解决方案，因此上述情况不需要特殊处理。

您可以通过检查连接所需点与多边形顶点形成的区域是否与多边形本身的面积相匹配来实现此操作。

或者，您可以检查从您的点到每对相邻多边形顶点到您的检查点的内角之和是否为360，但我觉得第一种选项更快，因为它不涉及除法或反三角函数的计算。

如果您的多边形内部有一个洞，我不知道会发生什么，但我觉得主要思想可以适应这种情况

您也可以在数学社区中发布问题。我敢打赌他们有一百万种方法来做到这一点

答案取决于您是否拥有简单或复杂的多边形。简单多边形不能有任何线段相交。因此，它们可以有孔洞，但是线条不能相互交叉。复杂区域可以有线段交点-因此它们可以具有重叠区域或仅通过单个点接触彼此的区域。

对于简单多边形，最佳算法是射线投射（穿越数）算法。对于复杂多边形，该算法无法检测位于重叠区域内部的点。因此，对于复杂多边形，您必须使用绕数算法。

这里有一篇优秀的文章，其中包含两种算法的C实现。我尝试过它们，它们运行良好。

http://geomalgorithms.com/a03-_inclusion.html

from typing import Iterable

def pnpoly(verts, x, y):
    #check if x and/or y is iterable
    xit, yit = isinstance(x, Iterable), isinstance(y, Iterable)
    #if not iterable, make an iterable of length 1
    X = x if xit else (x, )
    Y = y if yit else (y, )
    #store verts length as a range to juggle j
    r = range(len(verts))
    #final results if x or y is iterable
    results = []
    #traverse x and y coordinates
    for xp in X:
        for yp in Y:
            c = 0 #reset c at every new position
            for i in r:
                j = r[i-1] #set j to position before i
                #store a few arguments to shorten the if statement
                yneq       = (verts[i][1] > yp) != (verts[j][1] > yp)
                xofs, yofs = (verts[j][0] - verts[i][0]), (verts[j][1] - verts[i][1])
                #if we have crossed a line, increment c
                if (yneq and (xp < xofs * (yp - verts[i][1]) / yofs + verts[i][0])):
                    c += 1
            #if c is odd store the coordinates        
            if c%2:
                results.append((xp, yp))
    #return either coordinates or a bool, depending if x or y was an iterable
    return results if (xit or yit) else bool(c%2)

这个 Python 版本非常灵活。您可以输入单个 x 和单个 y 值以获得 True/False 结果，或者您可以使用 range 来遍历整个点网格的 x 和 y。如果使用范围，则返回所有 True 点的 x/y 对的 list。 vertices 参数需要一个二维的 Iterable 的 x/y 对，例如：[(x1,y1), (x2,y2), ...]

示例用法：

vertices = [(25,25), (75,25), (75,75), (25,75)]
pnpoly(vertices, 50, 50) #True
pnpoly(vertices, range(100), range(100)) #[(25,25), (25,26), (25,27), ...]

实际上，即使这些也可以工作。

pnpoly(vertices, 50, range(100)) #check 0 to 99 y at x of 50
pnpoly(vertices, range(100), 50) #check 0 to 99 x at y of 50

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

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

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

GitHub Gist的链接。

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

它是新颖而优雅的。我们花费 O(n * log(n)) 的时间来构建一个表，使我们能够在 O(log(n) + k) 的时间内测试点是否在多边形内。

与射线追踪或角度不同，使用扫描光束表可以显著加快对同一多边形进行多次检查的结果。我们必须预先构建一个扫描光束活动边缘表，这就是大部分代码所做的工作。

我们计算在y方向上的扫描线和活动边。我们按照它们的y分量排序，生成一个点列表，并且对于两个事件（Start-Y和End-Y），我们遍历这个列表并跟踪处理列表时的活动边。我们按顺序处理每个事件，对于每个扫描线，我们记录事件的y值以及每个事件（即start-y和end-y）处的活动边，但仅当我们的事件y与上次不同时才记录这些内容（因此在我们将其标记为表中的内容之前，事件点上的所有内容都已被处理）。

因此，我们得到了以下表格：

1. []
2. p6p5，p6p7
3. p6p5，p6p7，p2p3，p2p1
4. p6p7，p5p4，p2p3，p3p1
5. p7p8，p5p4，p2p3，p2p1
6. p7p8，p5p4，p3p4，p2p1
7. p7p8，p2p1
8. p7p8，p1p0
9. p8p0，p1p0
10. []

实际执行工作的代码在构建此表格后仅有几行。

• 注意：此处的代码使用复数值作为点。因此，.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) 时间。每次检查一个点时，它会迭代整个多边形。因此，即使这种明显不太优化的实现方式，仍然完胜其他所有方法。

有一些性能技巧可以使用，例如，我们可以将时间复杂度降低到O(log(n) + log(k))。为了做到这一点，我们会在扫描线中实现 Bentley-Ottmann，并且不再将交点处理为不同的事件，而是在交点处分割线段。然后，我们还按照它们的x截距对活动边缘进行排序。我们知道，在扫描线期间不会发生任何交点，因为我们已经对线段进行了排序（即使它们从相同的初始点开始，也要正确地对它们进行排序（需要查看斜率或仅比较线段的中点）。然后，我们就有了一个排序好的无交点活动列表扫描线表，这意味着我们也可以在活动边缘列表中进行二进制搜索。虽然这听起来很费力，但对于通常为2或4的k值来说，这是值得的。

此外，由于这基本上变成了一个查找表和一些最小的x截距计算，所以更适合使用GPU。您不再需要循环遍历多边形。因此，您实际上可以使用像numpy这样的工具批量计算这些点，从而获得同时执行所有计算的性能提升。

这似乎在R中可以工作（对于丑陋的外观表示歉意，希望能看到更好的版本！）。

pnpoly <- function(nvert,vertx,verty,testx,testy){
          c <- FALSE
          j <- nvert 
          for (i in 1:nvert){
              if( ((verty[i]>testy) != (verty[j]>testy)) && 
   (testx < (vertx[j]-vertx[i])*(testy-verty[i])/(verty[j]-verty[i])+vertx[i]))
            {c <- !c}
             j <- i}
   return(c)}