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

numpy: get the indices within a polygon without creating a mask

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

如何获取从这些点生成的多边形内的所有索引(或坐标)?

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

你可以使用 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]

使用口罩可能是您所能达到的最有效的方法。这是一些效率相当低的算法,但可能可以优化以接近掩码方法。这本质上做了一个掩码,但在线上。

方法是:

  1. 求所有边的直线方程

  2. 查找边界框

  3. 对于边界框内的每个 y(或 x,以较小者为准),计算在该 y 处与水平线 (y=yi) 相交的边,并找出它们与 x 相交的边.

  4. 对于边界框内的每个 x,找出 x 右侧与直线 y=yi 相交的边数。如果边数为奇数,则点 (x,y) 位于多边形内部。

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

import numpy as np

# taken from: 
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)