【问题标题】:Calculate euclidean distance with numpy用numpy计算欧几里得距离
【发布时间】:2015-09-23 22:26:15
【问题描述】:

我有一个点集,我将它的坐标存储在三个不同的数组(xa、ya、za)中。现在,我想计算这个点集(xa[0]、ya[0]、za[0] 等)的每个点与另一个点集(xb、yb、zb)的所有点之间的欧几里得距离) 并且每次都将最小距离存储在一个新数组中。

假设 xa.shape = (11,), ya.shape = (11,), za.shape= (11,)。分别为xb.shape = (13,), yb.shape = (13,), zb.shape = (13,)。我想做的是每次取一个xa[],ya[],za[],并计算它与xb,yb,zb的所有元素的距离,最后将最小值存储到xfinal中。形状 = (11,) 数组。

你认为 numpy 可以做到这一点吗?

【问题讨论】:

  • 换句话说,对于每个xa/ya/za,您想计算到xb/yb/zb中最近点的距离?
  • 是的,如果它会更容易一些......

标签: python numpy


【解决方案1】:

另一种解决方案是使用 scipy 的空间模块,尤其是 KDTree。

这个类从一组数据中学习,并且可以在给定一个新数据集的情况下进行询问:

from scipy.spatial import KDTree
# create some fake data
x = arange(20)
y = rand(20)
z = x**2
# put them togheter, should have a form [n_points, n_dimension]
data = np.vstack([x, y, z]).T
# create the KDTree
kd = KDTree(data)

现在,如果您有一个点,您可以通过以下方式询问最近点(或 N 个最近点)的距离和索引:

kd.query([1, 2, 3])
# (1.8650720813822905, 2)
# your may differs

或者,给定一个位置数组:

#bogus position
x2 = rand(20)*20
y2 = rand(20)*20
z2 = rand(20)*20
# join them togheter as the input
data2 = np.vstack([x2, y2, z2]).T
#query them
kd.query(data2)

#(array([ 14.96118553,   9.15924813,  16.08269197,  21.50037074,
#    18.14665096,  13.81840533,  17.464429  ,  13.29368755,
#    20.22427196,   9.95286671,   5.326888  ,  17.00112683,
#     3.66931946,  20.370496  ,  13.4808055 ,  11.92078034,
#     5.58668204,  20.20004206,   5.41354322,   4.25145521]),
#array([4, 3, 2, 4, 2, 2, 4, 2, 3, 3, 2, 3, 4, 4, 3, 3, 3, 4, 4, 4]))

【讨论】:

  • 远远优于其他答案,您还可以使用更快的 cKDTree(使用新的 scipy 可能相同)。
  • 好的,显然,这个选项效果更好。但是,最后我不知道这是否可以在我运行这些脚本的 abaqus 中工作。您知道如何从 kd.query 中提取具有最小值的数组吗?提前非常感谢。
  • 好的,我做到了。这不是那么困难。再次非常感谢。
【解决方案2】:

您可以使用np.subtract.outer(xa, xb) 计算每个 xa 到每个 xb 的差异。到最近的 xb 的距离由下式给出

np.min(np.abs(np.subtract.outer(xa, xb)), axis=1)

要将其扩展到 3D,

distances = np.sqrt(np.subtract.outer(xa, xb)**2 + \
    np.subtract.outer(ya, yb)**2 + np.subtract.outer(za, zb)**2)
distance_to_nearest = np.min(distances, axis=1)

如果您真的想知道 哪个 b 点是最近的,请使用 argmin 代替 min

index_of_nearest = np.argmin(distances, axis=1)

【讨论】:

  • 是的,我认为它有效。它似乎非常快。我的数组很大!因此,我认为我必须遍历 xa,ya,za 长度,获取每一行并计算它与整个 xb,yb,zb 的距离。我认为它有效。我会交叉检查它的有效性,我会告诉你的。无论如何,非常感谢!
  • 如果您可以使用 scipy,请参阅@EnricoGiampieri 的解决方案。
  • 不幸的是,它说数组太大了。它似乎不适用于我的数组。
【解决方案3】:

有不止一种方法可以做到这一点。最重要的是,内存使用和速度之间存在权衡。这是浪费的方法:

s = (1, -1)
d = min((xa.reshape(s)-xb.reshape(s).T)**2
     + (ya.reshape(s)-yb.reshape(s).T)**2
     + (za.reshape(s)-zb.reshape(s).T)**2), axis=0)

另一种方法是遍历b 中设置的点,以避免扩展到完整的矩阵。

【讨论】: