计算以每个数据点为中心的固定半径球内数据点数的有效方法

Question

我有一个包含许多数据点的数据库，每个数据点都有一个 x,y,z 坐标。我想计算与相邻点一定距离内的点数。有些点会有一对在半径 R 内，有些则没有。我只是想计算一定距离内的对数。我可以轻松地编写一个算法来执行此操作，但效率不够（因为我会遍历每个数据点）。

这似乎已经存在于astropy、scipy等中，但我似乎找不到我要找的东西。有什么东西可以做到这一点吗？

Answer 1

我没有这方面的直接经验，但 scipy.spatial.distance.pdist 可能就是您要找的东西。

这个也可能有帮助。它给出了一个我理解的如何解决问题的例子。

Answer 2

正如@Davis Herring 在评论中提到的，一个有效的选择是 k-d 树。

k-d 树是一种算法，它避免了 brute-force 方法并允许进行有效的距离计算*（背景见答案底部）。

有几种 Python 实现方式，其中一种是通过 SciPy:

SciPy k-d tree in Cython（更快，因为它使用 C/Cython）

SciPy k-d tree in pure Python

您可以通过首先为您的 xyz 数据构建一个 k-d 树来使用它：

import numpy as np  #for later code
from scipy.spatial import cKDTree

kdtree = cKDTree(xyzData)

然后，您必须查询具有点 point 的 k-d 树，以计算 point 与其最近邻居之间的距离。此查询的输出是 point 与其最近邻居之间的距离 NN_dist 以及该邻居的索引 NN_idx。要为你的所有点计算这个，我们需要一个 for 循环，但是考虑到 k-d 树算法，这比 brute-force 计算快得多：

NN_dists = np.zeros(numPoints)  #pre-allocate an array to store distances
for i in range(numPoints):
    point = xyzData[i]

    NN_dist, NN_idx = kdtree.query(point,k=[1])

    #Note: 'k' specifies the kth neighbor distance to compute, 
    #so set k=2 if you end up finding the point as its own "neighbor":
    if NN_dist == 0:
        NN_dist, NN_idx = targetTree.query(curCoord,k=[2])
    
    NN_dists[i] = NN_dist

（有关详细信息，请参阅 k-d tree query）。

然后，要找到低于某个阈值的距离，您可以在使用比较运算符（如 <）时使用 NumPy 数组的 built-in 实用程序：

distanceThres = 10
goodIdx = NN_dists < distanceThres
goodPoints = xyzData[goodIdx]

这将为您提供指定距离阈值 distanceThres 内的索引 goodIdx 和点 goodPoints（尽管您可能需要根据 shape/format 你的 xyz 坐标数据）。

*k-d 树上的浅色背景（掩盖了细节——更多信息请参阅参考资料）：k-d 树方法以这样一种方式对数据集进行分区，避免计算它们之间的距离每个点（即蛮力法）。它通过将数据集划分为二进制 space 分区以构建 k-d 树来实现此目的。这些分区使得距离计算（例如 nearest-neighbor 搜索）可以忽略远处分区中的数据点。此外，同样的 k-d 树被重复用于每个点。

一般来说，网上有很多关于 k-d 树的资源。当我学习这个算法时，我发现这些参考资料最有帮助：Stanford k-d trees or Princeton k-d trees.

如果您有任何问题，请告诉我 -- 我自己在天文学项目中遇到了这个问题，所以我可以提供更多帮助。

计算以每个数据点为中心的固定半径球内数据点数的有效方法

Efficient method for counting number of data points inside sphere of fixed radius centered on each data point

python

astronomy

astropy