Python：非凸网格的边界点

Question

我有以下情况，如下图所示：

我希望找出包围红点的网格点。红点是移动代理的轨迹。因此，在许多情况下，我们有一堆点，因此解决方案应该尽可能快速。

网格以点形式绘制。

第一步，我成功地减少了网格点的数量，如下所示（绘制为x）：

这是我的代码：

step = .5
gridX, gridY = np.meshgrid(np.arange(xmin-step, xmax+step, step), np.arange(ymin-step, ymax+step, step))
mask = False * np.empty_like(gridX, dtype=bool)
threshold = 0.5
for (x,y) in  zip(df_traj['X'], df_traj['Y']):
    pX = x * np.ones_like(gridX)
    pY = y * np.ones_like(gridY)
    distX = (pX - gridX)**2
    distY = (pY - gridY)**2
    dist = np.sqrt(distX + distY)
    condition = (dist < threshold)
    mask = mask | condition

gX = gridX*mask
gY = gridY*mask

• 第二步，这就是我需要一点帮助的地方：

如何高效地过滤掉网格内部的点，只保留“红区域”外部的“x点”？

编辑

在这种特殊情况下，我有92450个红点。

- Tengis

---

那么，你的输入是在变量df_traj中给出的？你大约有多少个点？“很多”？ - user202729

如果我理解你的问题正确，我认为你正在寻找你的点的非凸包。 - Cory Kramer

@gbtimmon，原帖明确指出边界不是凸的。 - user202729

听起来你需要一个凹壳算法。我一段时间前为Python实现了一个：hull（基于存储库中的一篇论文）。 - jsmolka

@jsmolka 看起来不错。我真的希望能找到一些简单的“掩码”，就像“x”点一样。谢谢，我会看看你的解决方案。 - Tengis

Answer 1

我认为如果你只是沿着边缘走，因为这是一个均匀间隔的网格，它应该可以工作。没有必要使用更复杂的非凸包来处理可以在任何地方的点。这不是针对你的代码进行调整的，我通过我的数据结构作弊，使代码易于编写，所以你需要处理一下，但我认为这个伪代码应该可以工作。

pnts = <<lists of points>>
edge_pnts = []
fpnt = pnt_with_min_x_then_min_y
cpnt = fpnt
npnt = None 

while npnt != fpnt:

    if   (cpnt[0] + 1, cpnt[1]    ) in pnts: npnt = (cpnt[0] + 1, cpnt[1]    )
    elif (cpnt[0] + 1, cpnt[1] + 1) in pnts: npnt = (cpnt[0] + 1, cpnt[1] + 1)
    elif (cpnt[0],     cpnt[1] + 1) in pnts: npnt = (cpnt[0]    , cpnt[1] + 1)
    elif (cpnt[0] - 1, cpnt[1] + 1) in pnts: npnt = (cpnt[0] - 1, cpnt[1] + 1)
    elif (cpnt[0] - 1, cpnt[1]    ) in pnts: npnt = (cpnt[0] - 1, cpnt[1]    )
    elif (cpnt[0] - 1, cpnt[1] - 1) in pnts: npnt = (cpnt[0] - 1, cpnt[1] - 1)
    elif (cpnt[0]    , cpnt[1] - 1) in pnts: npnt = (cpnt[0]    , cpnt[1] - 1)
    elif (cpnt[0] + 1, cpnt[1] - 1) in pnts: npnt = (cpnt[0] + 1, cpnt[1] - 1)
    else: raise ValueError("Oh no!")

    edge_pnts.append(npnt)
    cpnt = npnt

- gbtimmon

---

我认为你在这里假设了凸性。你应该考虑从哪个方向来，总是尝试从逆时针方向开始。 例如，在 OP 的例子中，你的算法将在点 (0.5,-0.5) 失败。 - shayelk

Answer 2

你只需要选择一个你知道在外壳上的点（我们可以选择最上面的最左边的点），并假设你是从上方到达的（因为我们知道它上面没有点）。 现在，当下一个点不在你的列表中时：

尝试从你来的方向逆时针移动。

代码看起来像这样：

matrix = [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
         [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
         [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

# Find the leftmost topmost point
first_point = None
for i in range(len(matrix)):
    if first_point:
        break
    for j in range(len(matrix[0])):
        if matrix[i][j]:
            first_point = [i, j]
            break

next_point = first_point
prev_direction = 'up'
next_direction_dict = {'up': 'left', 'left': 'down', 'down': 'right', 'right': 'up'}
opposite_direction = {'up': 'down', 'left': 'right', 'down': 'up', 'right': 'left'}
hull_points = []


def go_direction(point, direction):
    # Find the point to a given direction of a given point
    i = point[0]
    j = point[1]
    if direction == 'right':
        j += 1
    elif direction == 'up':
        i -= 1
    elif direction == 'left':
        j -= 1
    elif direction == 'down':
        i += 1
    else:
        raise ValueError
    return [i, j]


def find_next_point(matrix, point, prev_direction):
    next_direction = next_direction_dict[prev_direction]
    next_point = go_direction(point, next_direction)
    prev_direction = next_direction
    while not matrix[next_point[0]][next_point[1]]:
        next_direction = next_direction_dict[prev_direction]
        next_point = go_direction(point, next_direction)
        prev_direction = next_direction
    from_direction = opposite_direction[prev_direction]
    return next_point, from_direction

next_point, prev_direction = find_next_point(matrix, next_point, prev_direction)
while next_point != first_point:
    if next_point not in hull_points:
        hull_points.append(next_point)
    next_point, prev_direction = find_next_point(matrix, next_point, prev_direction)

编辑：现在还可以处理单点“触手”，通过迭代直到返回第一个点。

- shayelk

Answer 3

对于非凸多边形，例如您的示例，凸包不是解决方案。我的建议是，考虑到您已经有了离散网格，当样本出现在内部时，只需将值False赋给布尔网格单元格即可。像这样：

import numpy as np
import matplotlib.pyplot as plt

# Generic data production
X, Y = np.random.normal(0, 1, 100000), np.random.normal(0, 1, 100000)
ind = np.where((X > 0) & (Y > 0))
X[ind] = 0
Y[ind] = 0    

# Generic grid definition
step  = 0.5
xmin, xmax = X.min(), X.max()
ymin, ymax = Y.min(), Y.max()
firstx = xmin-step/2
firsty = ymin-step/2
lastx  = xmax+2*step/2
lasty  = ymax+2*step/2
gridX, gridY = np.meshgrid(np.arange(firstx, lastx, step), np.arange(firsty, lasty, step))

# This is the actual code that computes inside or outside
bool_grid = np.ones(gridX.shape, dtype="bool")
bool_grid[np.int_(0.5+(Y-firsty)/step), np.int_(0.5+(X-firstx)/step)] = False

# Plot code
plt.scatter(gridX.flatten(), gridY.flatten(), marker="+", color="black", alpha=0.3)
plt.scatter(gridX[bool_grid].flatten(), gridY[bool_grid].flatten(), marker="+", s=90, color="green")
plt.scatter(X, Y, s=10, color="red")
plt.show()

这将导致以下结果（绿色的叉子表示True值）：

注意：这很快，但有一些限制。 如果您的数据不紧凑，则会在形状内部有 True 值（因此可能存在孔）。可以处理图像以消除孔洞，例如使用泛洪填充或基于移动窗口的算法。另一个选择是调整网格的分辨率。

- armatita

---

如果不是分析形状，而是分析周围的空间并找到所有连接在外部的点，那么孔洞问题就迎刃而解了，你可以很快地找到边缘。你可能需要添加一些填充来确保有一个连续的外部区域，但这很简单。 - gbtimmon

等一下，我不明白你在这里是如何找到边缘的？ - gbtimmon

@gbtimmon 这是因为你在想形状（多边形），而算法是针对点的。把网格的单元格看作盒子。对于每个样本，计算其在网格上的位置，并将该“盒子”变为false。这是一种与广播兼容的操作。这就是为什么它如此快的原因。形状将是自然结果。至于孔洞，我的注释只是一个关于代码可能不完全明显的警告。但据我所知，“可能”的孔洞实际上可以成为OP的积极因素。 - armatita

根据我的理解，OP正在寻找与我示例图中的绿色十字标记相当的内容（因此是gridX [bool_grid]和gridY [bool_grid]）。至少这是我从他的问题中解释出来的：“如何有效地过滤掉网格内部的点并仅保留“红区域”外部的“x点”？” - armatita

太简单了，他已经有了内部点。那么就是所有点-内部点。我认为他想过滤多边形的内部点，得到一组边界点。 - gbtimmon

Answer 4

这是我脑海中浮现的另一个想法：

对空间进行泛洪填充：

pnts = <<lists of points>>
seen = set()
edges = []
stack = (0,0)

while stack:
    ele = stack.pop()
    if ele in pnts:
        edges.append(ele)
    else:
        seen.add(ele)
        if (ele[0] + 1, ele[1]) not in seen: 
            stack.append(ele[0] + 1, ele[1])
        if (ele[0] - 1, ele[1]) not in seen: 
            stack.append(ele[0] - 1, ele[1])
        if (ele[0], ele[1] + 1) not in seen: 
            stack.append(ele[0], ele[1] + 1)
        if (ele[0], ele[1] - 1) not in seen: 
            stack.append(ele[0], ele[1] - 1)

然后你需要对点进行排序，这应该不会太难。

- gbtimmon