使用 numpy 进行循环测试的快速点

Question

我有大量具有整数坐标的(x,y) 网格点，我想测试它们是否位于由半径和中心给出的少量圆中。这些点是图像的一些标记部分，这意味着有少量不规则形状的块，其中包含这些点。在那里我想检查碰撞并计算圆内的点数。我目前的方法相当慢（使用 python 和 numpy）。

现在我有两个任务：

1. 测试，如果集合 A 的任何点在任何圆内
2. 计算集合 B 中的点数，这些点在一个圆圈中

我当前的实现如下所示（setA 和 setB 是 Nx2 numpy 数组，center 是 1x2 数组。）：

1) 为每个圆创建一个point - center 的数组，将其按元素平方并求和，然后检查它是否小于radius**2

for circle in circles:
    if (((setA - circle.center)**2).sum(axis=1) < circle.radius**2).any():
        return "collision"
return "no collision"

这可以通过使用 python 循环并在第一次碰撞时中断来优化，但通常 numpy 循环比 python 循环快很多，实际上两个版本都比预期慢。

2) 为每个圆创建一个距离数组并进行元素小于半径测试。将所有数组相加并计算结果的非零元素。

pixels = sp.zeros(len(setB))
for circle in circles:
    pixels += (((setB - circle.center)**2).sum(axis=1) < circle.radius**2)
return np.count_nonzero(pixels)

有没有一种简单的方法可以加快速度？

我不想过度优化（并使程序变得更复杂），只想以最有效的方式使用 numpy，尽可能使用 numpy 矢量化。

所以构建最完美的空间树或类似的空间树不是我的目标，但我认为在普通桌面上应该可以以最快的方式处理几千个点和 10-20 个圆的 O(n^2) 算法今天的电脑。

Answer 1

利用坐标是整数：

创建查找图像

radius = max([circle.radius for circle in circles])
mask = np.zeros((image.shape[0] + 2*radius, image.shape[1] + 2*radius), dtype=int)
for circle in circles:
    center = circle.center + radius
    mask[center[0]-circle.radius:center[0]+circle.radius + 1,
         center[1]-circle.radius:center[1]+circle.radius + 1] += circle.mask

circle.mask 是一个小正方形补丁，包含内部点圆盘的掩码

计算碰撞次数现在就像

mask[radius:-radius, radius:-radius][setB[:,0], setB[:,1]].sum()

快速创建光盘（无乘法，无平方根）：

r = circle.radius
h2 = np.r_[0, np.add.accumulate(np.arange(1, 2*r+1, 2))]
w = np.searchsorted(h2[-1] - h2[::-1], h2)
q = np.zeros(((r+1), (r+1)), dtype=int)
q[np.arange(r+1), w[::-1]] = 1
q[1:, 0] -= 1
q = np.add.accumulate(q.ravel()).reshape(r+1, r+1)
h = np.c_[q, q[:, -2::-1]]
circle.mask = np.r_[h, h[-2::-1]]