列出到线的距离内的所有点

Question

我有一个坐标列表a。我想指定半径r，然后列出a 中任意点的指定半径内的二维网格上的所有点，以及每个网格点到a 中任意点的最小距离。由于a 很大（~1000-2000 点），我希望尽可能高效。

到目前为止，我的实现使用this code 列出给定点半径内的所有坐标；然后迭代a中的所有坐标；然后展平输出，取一个集合（会有很多重复） - 这是在线上任意点的指定半径内的一组坐标 - 然后使用 scipy.spatial.distance.cdist 计算该集合到 a 的最小距离：

import numpy as np
np.random.seed(123)
a = np.random.randint(50, size=(100, 2))

def collect(coord, sigma: float =3.0):
    """ create a small collection of points in a neighborhood of some point 
    """
    neighborhood = []
    
    x=coord[0]
    y=coord[1]

    X = int(sigma)
    for i in range(-X, X + 1):
        Y = int(pow(sigma * sigma - i * i, 1/2))
        for j in range(-Y, Y + 1):
            neighborhood.append((x + i, y + j))

    return neighborhood

rad = 10
coord_list = [collect(coord, rad) for coord in a]
coord_set = np.array(list(set([val for sublist in coord_list for val in sublist])))

from scipy.spatial import distance

dists = distance.cdist(coord_set, a).min(axis=1)

结果：

coord_set
> array([[26, 21],
       [18, 17],
       [50,  6],
       ...,
       [14, 47],
       [15, 12],
       [ 7,  8]])

dists
> array([2.23606798, 5.65685425, 1.41421356, ..., 3.16227766, 3.        ,
       1.41421356])

有没有人有什么办法可以改善这一点并更有效地做到这一点？

Answer 1

让我们仔细检查您的实施。请注意，您已经准备好并部分计算 collect 中的距离度量。

1. 如果您将neighborhood 设为dict 而不是list，并由具有最小距离值的网格点键入会怎样。您可以完全消除对set 和cdist 的调用。
2. 如果a 可以包含浮点值，则应将范围从int(coord[0] - rad) 到int(coord[0] + rad) + 1。 int(0.5 - 10) 是 -9，而 int(0.5) - 10 是 -10。
3. 您可以与平方半径进行比较，因此您不需要取平方根，只需要一次即可获得最终结果。

第 2 点和第 3 点是相对较小的改进。

这是一个例子：

rad = 10
rad2 = rad**2

points = {}

for coord in a:
    for x in range(int(np.ceil(coord[0] - rad)), int(coord[0] + rad) + 1):
        dx2 = (x - coord[0])**2
        chord = np.sqrt(rad2 - dx2)
        for y in range(int(np.ceil(coord[1] - chord)), int(coord[1] + chord) + 1):
            d2 = dx2 + (y - coord[1])**2
            key = x, y
            points[key] = min(d2, points.get(key, rad2))

要将其转换为 numpy 数组：

grids = np.array(list(points.keys()))
nearest = np.sqrt(np.fromiter(points.values(), count=len(points), dtype=float))

Answer 2

您可以修改 the answer 我以非常直接的方式回答您的其他链接问题。结果也非常快（a 中的 10K 点约为 425 毫秒）。

编辑：对于稀疏的情况（实际过滤的点数只是整个网格的一小部分），我还在下面添加了sparse_grid_within_radius() 版本。

（重要：使用scipy >= 1.6.0，其中KDTree的Python实现已被cKDTree替换。参见release notes）。

# please use scipy >= 1.6.0
from scipy.spatial import KDTree


def grid_within_radius(a, radius=10, resolution=(1,1), origin=(0,0)):
    pmin = (a - origin).min(axis=0) - radius
    pmax = (a - origin).max(axis=0) + radius
    imin = np.floor(pmin / resolution)
    imax = np.ceil(pmax / resolution) + 1
    xv, yv = np.meshgrid(np.arange(imin[0], imax[0]), np.arange(imin[1], imax[1]))
    grid = np.stack((xv, yv), axis=-1) * resolution + origin
    dist, _ = KDTree(a).query(grid, k=1, distance_upper_bound=radius + 0.01)
    idx = dist <= radius
    return grid[idx], dist[idx]

用法

首先，OP的例子：

np.random.seed(123)
a = np.random.randint(10, size=(100, 2))
g, d = grid_within_radius(a)

为了与 OP 的结果进行比较，我们需要对他们的解决方案进行排序 (coord_set, dists)：

def sort2d(a, other=None):
    other = a if other is None else other
    return other[np.lexsort((a[:, 0], a[:, 1]))]

这样，我们可以检查我们的解决方案是否相同：

>>> np.allclose(g, sort2d(coord_set))
True

>>> np.allclose(d, sort2d(coord_set, dists))
True

还有另一个例子（使用不同的网格分辨率和半径）：

g, d = grid_within_radius(a, radius=0.6, resolution=(.11, .47))
plt.scatter(*a.T, s=10, c='r')
plt.scatter(*g.T, s=1)

速度

a = np.random.randint(1000, size=(10_000, 2))
%timeit grid_within_radius(a)
# 425 ms ± 528 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

稀疏版

当返回的点数占总网格的很大一部分（例如 30% 或更多）时，上述方法效果很好。但是对于非常稀疏的情况（例如，a、radius 和网格resolution 的边界框的组合导致生成一个巨大的网格，然后消除其中的大部分），那么它就很慢。为了说明，以下是a = np.random.randint(0, 200, (10, 2)) 的稀疏案例：

下面的版本通过在量化位置而不是整个网格周围生成网格“补丁”来解决这个问题。

import numpy as np
import pandas as pd
from scipy.spatial import KDTree


def unique_rows(a):
    # np.unique(a, axis=0) is slow, in part because it sorts;
    # using pandas.drop_duplicates() instead.
    # see https://github.com/numpy/numpy/issues/11136#issuecomment-822000680
    return pd.DataFrame(a).drop_duplicates().values

def sparse_grid_within_radius(a, radius=10, resolution=(1,1), origin=(0,0), rbox=1):
    resolution = np.array(resolution)
    box_size = radius * rbox
    nxy0 = np.floor(radius / resolution).astype(int)
    nxy1 = np.ceil((box_size + radius) / resolution).astype(int) + 1
    patch = np.stack(list(map(np.ravel, np.indices(nxy0 + nxy1))), axis=1) - nxy0
    
    ar = np.floor((a - origin) / box_size) * box_size
    ia = unique_rows(np.floor(ar / resolution).astype(int))
    grid = unique_rows((ia[:, None] + patch).reshape(-1, 2)) * resolution + origin

    dist, _ = KDTree(a).query(grid, k=1, distance_upper_bound=radius * 1.01)
    idx = dist <= radius
    return grid[idx], dist[idx]

这样，即使是非常稀疏的结果也很快。

示例

a = np.random.randint(0, 4000, (100, 2))

%timeit g, d = grid_within_radius(a)
# 3.88 s ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit gs, ds = sparse_grid_within_radius(a)
# 29.6 ms ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

速度对比

import perfplot

perfplot.show(
    setup=lambda n: np.random.randint(0, n, (100, 2)),
    kernels=[grid_within_radius, sparse_grid_within_radius],
    n_range=[2 ** k for k in range(14)],
    equality_check=lambda a, b: np.allclose(sort2d(a[0]), sort2d(b[0])),
)