Python中曲线的距离矩阵

Question

我有一组由2D数组（点数，坐标数）定义的曲线。我正在使用豪斯多夫距离为它们计算距离矩阵。我的当前代码如下。不幸的是，当每个曲线都有50-100个3D点时，500-600个曲线的速度太慢了。有没有更快的方法？

def distanceBetweenCurves(C1, C2):
    D = scipy.spatial.distance.cdist(C1, C2, 'euclidean')

    #none symmetric Hausdorff distances
    H1 = np.max(np.min(D, axis=1))
    H2 = np.max(np.min(D, axis=0))

    return (H1 + H2) / 2.

def distanceMatrixOfCurves(Curves):
    numC = len(Curves)

    D = np.zeros((numC, numC))
    for i in range(0, numC-1):
        for j in range(i+1, numC):
            D[i, j] = D[j, i] = distanceBetweenCurves(Curves[i], Curves[j])

    return D

- ahmethungari

---

scipy.spatial.distance.cdist 是慢的部分还是 distanceMatrixOfCurves 里面的双重循环？如果这些曲线是凸的，那么可能可以优化第一个潜在的缓慢部分。这些曲线是否相交或包含在其他曲线中？我觉得你可以重用早先发现的距离来加速新的计算。当然，这只是胡言乱语，我自己也遇到了类似于 min(min(..)) 测量的问题，并且难以将这些考虑到一般情况中。你有尝试过或思考过代码之外的内容吗？ - mmgp

@ahmethungari，你应该对你的代码进行性能分析，以确保（使用cProfile + runsnakerun来解释结果非常好）找到确切的瓶颈所在。我对这些事情没有很好的感觉，但是如果你可以添加一些生成一些小例子数据的代码，那么帮助你就会更容易了。 - YXD

1. 你真的需要完整的矩阵D吗？还是只需要上三角或下三角矩阵就可以了？这种形式的D[i,j] = D[j,i] =...对于数据局部性来说绝对不好。
2. 你尝试过使用列表推导式或map代替双重循环吗？

- ev-br

这并没有以任何方式减少计算复杂度。就目前而言，我认为原帖作者只会得到一些实现上的小改进，例如使用库/语言/包X，因为它可以比原来运行Y（其中Y可能是相同的方法，也许有一些小修改）快！等等。我希望我在这里被证明是错误的。 - mmgp

@Zhenya 我不需要整个矩阵，你是对的。你认为列表推导比循环快吗？我可以试一下... - ahmethungari

Answer 1

你的问题可能与这个问题有关。

这是一个比较困难的问题。一种可能的方法是自己实现欧几里得距离，完全放弃scipy并利用pypy的JIT编译器。但最有可能的情况是这样做不会让你获得太多好处。

个人而言，我建议您使用C语言编写程序。

问题不在于实现，而在于解决问题的方式。您选择了通过计算度量空间子集中每对不同点的欧几里得距离来采用暴力方法。这是计算密集型的：

假设您有500条曲线，每条曲线都有75个点。使用暴力方法，您将需要计算500 * 499 * 75 * 75 = 1,403,437,500次欧几里得距离。因此，这种方法需要很长时间才能运行。

我不是这方面的专家，但我知道Hausdorff距离在图像处理中被广泛使用。我建议您浏览文献以获取速度优化算法。一个起点可能是这篇文章，或者这篇文章。此外，在Hausdorff距离中经常与之提到的是Voroni图。

希望这些链接可以帮助您解决问题。

- j-i-l

---

第一篇链接的论文很有趣，解决了问题的一部分。现在另一个问题可以或多或少地被表述为：给定三个集合A、B、C，如果我们知道从A到B的豪斯多夫距离以及从该计算得出的内部细节，是否可能以比计算A和B更小的计算复杂度找到从A到C的豪斯多夫距离，如果存在从A到C的路径（考虑某种度量），跨越B吗？好吧，这比我开始写的要长得多。希望有人能理解。 - mmgp

感谢您提供的所有建议。我明白了：“使用更智能的算法”。希望在使用当前版本完成论文截止日期后，我能够实现它 :) - ahmethungari

@ahmethungari 很抱歉我不能为您提供更多帮助，但是这个问题没有简单的解决方案。我有一个加速计算的想法：不要计算Curve1中给定点与Curve2中所有点之间的euclid.dist.，而是随机选择Curve2中的一个点，计算euclid.dist.（称其为r），然后对于Curve2中的每个其他点，可以计算曼哈顿距离（计算速度更快），并检查每个维度的曼哈顿距离是否大于r，如果是... - j-i-l

@ahmethungari 任意选择的点更接近，因此您无需计算欧几里得距离（如果您不确定原因：三角形不等式）。如果不是这样，请计算欧几里得距离。如果它小于“r”，请更新“r”。您仍将通过Curve2中的所有点，但通常只需计算曼哈顿距离，这更快。...好吧，这只是一个想法，肯定有更复杂的方法来加速这个过程。 - j-i-l

如果您对此答案感到满意，请考虑接受它。 :) - j-i-l

Answer 2

我最近在类似问题的这里回复了： 3D网格之间的Hausdorff距离 希望这可以帮到你，我需要对25 x 25,000个点进行成对比较（总共是25 x 25 x 25,000个点），我的代码运行时间从1分钟到3-4小时不等（取决于点数）。我没有看到数学上提高速度的太多选择。

另一种选择是使用不同的编程语言（C / C ++）或将此计算引入GPU（CUDA）。我现在正在尝试CUDA方法。

2015年3月12日编辑：

我通过进行并行CPU计算来加速此比较。那是最快的方法。我使用了“pp”软件包（“parallel python”）的很好的示例，并在三台不同的计算机和python组合上运行。 不幸的是，我一直遇到python 2.7 32位的内存错误，因此我安装了WinPython 2.7 64位和一些实验性的numpy 64位软件包。

所以对我来说，这个工作相当有帮助，而且没有CUDA那么复杂...祝好运

- Akos Gulyban

Answer 3

你可以尝试以下几种方法：

1. 使用numpy-MKL，它使用英特尔高性能数学核心库而不是numpy；
2. 使用Bootleneck进行数组函数处理；
3. 使用Cpython进行计算。

- foool