numpy：在不创建掩码的情况下获取多边形内的索引答案

【问题标题】：numpy: get the indices within a polygon without creating a masknumpy：在不创建掩码的情况下获取多边形内的索引
【发布时间】：2018-02-07 12:37:33
【问题描述】：

我有一个多边形矢量，格式如下(x1,y1,x2,y2, ...., xn,yn)。例如，考虑这个多边形数组： polyPoints = [3,5,7,8,9,5]

如何获取从这些点生成的多边形内的所有索引（或坐标）？

到目前为止，我看到的答案要求您先创建 2D 蒙版，然后才能获取多边形内的索引。

【问题讨论】：

所有索引（x,y坐标）都是整数吗？
是的，(x,y) 索引
直截了当的蛮力方法可能是迭代从 min(x)、min(y) 到 max(x)、max(y) 的所有整数点，并为每个点测试是否在内部或外部多边形。要对每个点进行测试，请使用光线投射或缠绕数算法。有关这些测试，请参见此处：en.wikipedia.org/wiki/Point_in_polygon
谢谢！但这听起来效率很低。

标签： python numpy 2d polygon

【解决方案1】：

你可以使用 scikit-image:

import numpy as np
from skimage.draw import polygon
points = [3,5,7,8,9,5]
r, c = polygon(points[1::2], points[::2])
print(r, c)

输出是：

[5 5 5 5 5 5 6 6 6 6 7 7] [3 4 5 6 7 8 5 6 7 8 6 7]

【讨论】：

【解决方案2】：

使用口罩可能会尽可能高效。这是一些效率相当低的算法，但可能可以优化为接近掩码方法。这本质上是做一个面具，但在线条上。

方法是：

求所有边线的方程
寻找边界框
对于边界框内的每个 y（或 x，以较小者为准），计算在该 y 处与水平线 (y=yi) 相交的边，并找出它们在哪个 x 处相交。
对于边界框内的每个 x，找出 x 右侧与 y=yi 相交的边数。如果边数是奇数，则点 (x,y) 在多边形内。

它确实适用于简单的方形几何体。

import numpy as np

# taken from: https://stackoverflow.com/questions/20677795/how-do-i-compute-the-intersection-point-of-two-lines-in-python
def line(p1, p2):
    A = (p1[1] - p2[1])
    B = (p2[0] - p1[0])
    C = (p1[0]*p2[1] - p2[0]*p1[1])
    return A, B, -C


def intersection(L1, L2):
    D  = L1[0] * L2[1] - L1[1] * L2[0]
    Dx = L1[2] * L2[1] - L1[1] * L2[2]
    Dy = L1[0] * L2[2] - L1[2] * L2[0]
    if D != 0:
        x = Dx / D
        y = Dy / D
        return x,y
    else:
        return False


# polyPoints = np.array([0, 0, 4, 0,4, 4, 0, 4])
polyPoints = np.array([[3,5,7,8,9,5]])
polyPoints = polyPoints.reshape(-1, 2)
npoints = polyPoints.shape[0]
polyEgdes = []

for i in range(npoints):
    point1, point2 = polyPoints[i, :], polyPoints[(i+1) % npoints, :]
    polyEgdes.append(line(point1, point2))

# bounding box
boundingBox = np.vstack((polyPoints.min(axis=0), polyPoints.max(axis=0)))


inside_points = []

for y in range(boundingBox[0, 1], boundingBox[1, 1]):
    x_intersect = []
    for l in polyEgdes:
        # y_ins should be same as y
        insect_point = intersection(l, [0, y, 0])
        if insect_point:
            x_intersect.append(insect_point[0])

    x_intersect = np.array(x_intersect)
    for x in range(boundingBox[0, 0]+1, boundingBox[1, 0]-1):
        x_int_points = x_intersect[(x_intersect - x) >= 0]

        if len(x_int_points) % 2 == 1:
            inside_points.append((x, y))

print(inside_points)

【讨论】：