两个不同 Numpy 数组中的点之间的最小欧几里得距离，不在

Posted 2023-02-23

【中文标题】两个不同 Numpy 数组中的点之间的最小欧几里得距离，不在

【英文标题】：Minimum Euclidean distance between points in two different Numpy arrays, not within

【发布时间】：2010-12-24 16:47:21

【问题描述】：

我有两个 x-y 坐标数组，我想在一个数组中找到 每个 点之间的最小欧几里得距离all 是另一个数组中的点。数组不一定大小相同。例如：

xy1=numpy.array(
[[  243,  3173],
[  525,  2997]])

xy2=numpy.array(
[[ 682, 2644],
[ 277, 2651],
[ 396, 2640]])

我当前的方法循环遍历xy1 中的每个坐标xy，并计算该坐标与其他坐标之间的距离。

mindist=numpy.zeros(len(xy1))
minid=numpy.zeros(len(xy1))

for i,xy in enumerate(xy1):
    dists=numpy.sqrt(numpy.sum((xy-xy2)**2,axis=1))
    mindist[i],minid[i]=dists.min(),dists.argmin()

有没有办法消除 for 循环并以某种方式在两个数组之间进行逐个元素的计算？我设想生成一个距离矩阵，我可以为其找到每行或每列中的最小元素。

另一种看待问题的方式。假设我将xy1（长度m）和xy2（长度p）连接成xy（长度n），然后我存储原始数组的长度。从理论上讲，我应该能够从那些我可以从中获取 m x p 子矩阵的坐标生成一个 n x n 距离矩阵。有没有办法有效地生成这个子矩阵？

【问题讨论】：

如果你需要加速你的代码，你应该删除不必要的numpy.sqrt（并且只有在你找到它时才取最小平方距离的平方根）。

【参考方案1】：

对于你想要做的事情：

dists = numpy.sqrt((xy1[:, 0, numpy.newaxis] - xy2[:, 0])**2 + (xy1[:, 1, numpy.newaxis - xy2[:, 1])**2)
mindist = numpy.min(dists, axis=1)
minid = numpy.argmin(dists, axis=1)

编辑：你可以使用numpy.hypot，而不是调用sqrt，做方块等

dists = numpy.hypot(xy1[:, 0, numpy.newaxis]-xy2[:, 0], xy1[:, 1, numpy.newaxis]-xy2[:, 1])

【讨论】：

天哪，太棒了。我没有意识到逐个元素也可以这样工作。所以xy1[:,0,numpy.newaxis] 通过列向量有效地替换了我的 for 循环，从中减去 xy2 的所有 x 值。非常棒，谢谢。

是的。有关更通用和优雅的方法，请参阅 Alex 的答案。

@fideli: help(numpy.subtract.outer) 告诉你 Alok 的 numpy.newaxis 技巧在 Alex 的回答中也起作用。

【参考方案2】：

要计算 m x p 的距离矩阵，应该可以：

>>> def distances(xy1, xy2):
...   d0 = numpy.subtract.outer(xy1[:,0], xy2[:,0])
...   d1 = numpy.subtract.outer(xy1[:,1], xy2[:,1])
...   return numpy.hypot(d0, d1)

.outer 调用生成两个这样的矩阵（沿两个轴的标量差异），.hypot 调用将它们转换为相同形状的矩阵（标量欧几里得距离）。

【讨论】：

在这种情况下我会选择 cdist，但是 +1 并且我从这个解决方案中学到了很酷的东西

@Alex：有没有一种有效的方法可以将其推广到两列以上？

@Mac，查看接受的答案和scipy.spatial.distance.cdist。

@Alex：谢谢，我确实看到了。我宁愿希望有一种方法可以特别概括你的代码，这样我就不必仅仅为那个函数添加对 SciPy 的依赖。还是谢谢。

@Mac，不确定我是否理解你的问题，但为什么不只是做 np.hypot(np.hypot(d0, d1), d2) 等（这里 d2 硬编码，只是为了举个例子）更多维度？

【参考方案3】：

（几个月后） scipy.spatial.distance.cdist( X, Y ) 给出所有成对的距离， X 和 Y 2 暗，3 暗 ... 它还做了 22 种不同的规范，详细 here.

# cdist example: (nx,dim) (ny,dim) -> (nx,ny)

from __future__ import division
import sys
import numpy as np
from scipy.spatial.distance import cdist

#...............................................................................
dim = 10
nx = 1000
ny = 100
metric = "euclidean"
seed = 1

    # change these params in sh or ipython: run this.py dim=3 ...
for arg in sys.argv[1:]:
    exec( arg )
np.random.seed(seed)
np.set_printoptions( 2, threshold=100, edgeitems=10, suppress=True )

title = "%s  dim %d  nx %d  ny %d  metric %s" % (
        __file__, dim, nx, ny, metric )
print "\n", title

#...............................................................................
X = np.random.uniform( 0, 1, size=(nx,dim) )
Y = np.random.uniform( 0, 1, size=(ny,dim) )
dist = cdist( X, Y, metric=metric )  # -> (nx, ny) distances
#...............................................................................

print "scipy.spatial.distance.cdist: X %s Y %s -> %s" % (
        X.shape, Y.shape, dist.shape )
print "dist average %.3g +- %.2g" % (dist.mean(), dist.std())
print "check: dist[0,3] %.3g == cdist( [X[0]], [Y[3]] ) %.3g" % (
        dist[0,3], cdist( [X[0]], [Y[3]] ))


# (trivia: how do pairwise distances between uniform-random points in the unit cube
# depend on the metric ? With the right scaling, not much at all:
# L1 / dim      ~ .33 +- .2/sqrt dim
# L2 / sqrt dim ~ .4 +- .2/sqrt dim
# Lmax / 2      ~ .4 +- .2/sqrt dim

【讨论】：

@denis cdist 计算所有对之间的距离。假设X 和Y 的长度相同N，如何仅在对应元素之间保持距离，例如[ dist(X[0],Y[0]), dist(X[1],Y[1]), ... dist(X[N],Y[N]) ]？

@LWZ，正是你所拥有的——np.array([ dist( x, y ) for x, y in zip( X, Y )])

这行得通！而且相当快。请注意，要计算的元素的长度必须为 2，否则 python 将引发错误。对于opencv2检测到的轮廓中的点列表，我需要先使用numpy的reshape函数对其进行reshape ...

@denis：好吧，例如cv2.findContours 的结果是这样的：[ [[x1 y1]] [[x2 y2]].....[[xn yn]] ]。我不知道他们为什么将坐标放在 2 方括号中，但是如果您将 cdist 直接应用于列表，每个元素的长度将只有 1（一个列表，其中包含 1 个列表），我必须对其进行整形所以轮廓成为长度为 2 元素的列表（这意味着平掉双括号）

@Jim Raynor，对 - cdist 期望 ndim == 2 的数组，应该检查

【参考方案4】：

import numpy as np
P = np.add.outer(np.sum(xy1**2, axis=1), np.sum(xy2**2, axis=1))
N = np.dot(xy1, xy2.T)
dists = np.sqrt(P - 2*N)

【讨论】：

这非常有效！这是唯一可以推广到任意维数的解决方案。

【参考方案5】：

接受的答案没有完全解决问题，它要求找到两组点之间的最小距离，而不是两组点之间每个点之间的距离套。

虽然对原始问题的直接解决方案确实包括计算每个对之间的距离并随后找到最小的，但如果一个人只对最小 距离。对于后一个问题，存在一个更快的解决方案。

所有提议的解决方案的运行时间都可以缩放为m*p = len(xy1)*len(xy2)。这对于小型数据集是可以的，但可以编写一个最佳解决方案，将其缩放为 m*log(p)，从而为大型 xy2 数据集节省大量成本。

可以使用scipy.spatial.cKDTree 实现这种最佳执行时间缩放，如下所示

import numpy as np
from scipy import spatial

xy1 = np.array(
    [[243,  3173],
     [525,  2997]])

xy2 = np.array(
    [[682, 2644],
     [277, 2651],
     [396, 2640]])

# This solution is optimal when xy2 is very large
tree = spatial.cKDTree(xy2)
mindist, minid = tree.query(xy1)
print(mindist)

# This solution by @denis is OK for small xy2
mindist = np.min(spatial.distance.cdist(xy1, xy2), axis=1)
print(mindist)

其中mindist 是xy1 中的每个点与xy2 中的点集之间的最小距离

【参考方案6】：

我认为下面的函数也可以。

import numpy as np
from typing import Optional
def pairwise_dist(X: np.ndarray, Y: Optional[np.ndarray] = None) -> np.ndarray:
    Y = X if Y is None else Y
    xx = (X ** 2).sum(axis = 1)[:, None]
    yy = (Y ** 2).sum(axis = 1)[:, None]
    return xx + yy.T - 2 * (X @ Y.T)

说明 假设X和Y的每一行是两组点的坐标。 让它们的大小分别为m X p 和p X n。 结果将产生一个大小为m X n 的numpy 数组，其中(i, j)-th 条目是i-th 行与X 和Y 的j-th 行之间的距离分别。

【参考方案7】：

我强烈建议使用一些内置方法来计算平方，并为它们定制根以优化计算方式并且非常安全地防止溢出。

@alex 下面的答案在溢出方面是最安全的，而且应该非常快。同样对于单点，您可以使用现在支持超过 2 维的 math.hypot。

>>> def distances(xy1, xy2):
...   d0 = numpy.subtract.outer(xy1[:,0], xy2[:,0])
...   d1 = numpy.subtract.outer(xy1[:,1], xy2[:,1])
...   return numpy.hypot(d0, d1)

安全问题

i, j, k = 1e+200, 1e+200, 1e+200
math.hypot(i, j, k)
# np.hypot for 2d points
# 1.7320508075688773e+200

np.sqrt(np.sum((np.array([i, j, k])) ** 2))
# RuntimeWarning: overflow encountered in square

overflow/underflow/speeds