一种快速计算三维矩阵中一个点到其余点距离的变换矩阵的方法

---

需求：

要快速获取一个三维矩阵中某个点到其余点的距离矩阵。

---

分析：

在二维矩阵中我们可以通过分解坐标，分别对每个坐标进行计算距离，之后用numpy进行对应维度的扩充和向量加法运算，我们可以得到二维矩阵中任意点到其余点的距离矩阵。附上相关文章链接：

一种利用numpy快速获取“二维矩阵中一个点到其余点距离矩阵”的计算方法 - 简书 (jianshu.com)

在立体几何学中，我们知道空间中一点到另一点的距离可以用以下公式表示：

可以看到公式中是把三个分量分开计算然后加在一起开方的。我们是不是也可以把一个高维矩阵拆成相应维度的向量，先对向量操作，最后再利用 numpy 的广播机制将若干向量融合在一起形成一个蕴含多点到一点欧式距离信息的高维矩阵。

---

实验：

代码：

import numpy as np
import time
from matplotlib import pyplot as plt
def getDisMat_H(c,h,w):
    strat = time.time()
    tmp = [np.arange(c), np.arange(h), np.arange(w)] 
    center = -0.5 + np.array([c, h, w])/2
    dis_mat = (np.expand_dims((tmp[0] - center[0])**2, axis=(1,2)) 
             + np.expand_dims((tmp[1] - center[1])**2, axis=(0,2))
             + np.expand_dims((tmp[2] - center[2])**2, axis=(0,1)))**0.5
    end = time.time()
    return end-strat

def getDisMat_L(c,h,w):
    strat = time.time()
    center = -0.5 + np.array([c, h, w])/2
    dis_mat = np.zeros((c,h,w))
    for i in range(c):
        for j in range(h):
            for k in range(w):
                dis_mat[i,j,k] = ((center[0]-i)**2 
                                + (center[1]-j)**2 
                                + (center[2]-k)**2)**0.5
    end = time.time()
    return end-strat

if __name__ == '__main__':
    X, Ltime, Htime = [], [], []
    for i in range(2,31):
        x = i**2 
        h,w,c = x, x, 10  # src mat shape
        X.append(x)
        Ltime.append(getDisMat_L(c,h,w))
        Htime.append(getDisMat_H(c,h,w))
    plt.figure()
    plt.subplot(211)
    plt.plot(X, Ltime, color="r", linestyle="-", marker="^", linewidth=1, label="for")
    plt.plot(X, Htime, color="g", linestyle="-", marker="s", linewidth=1, label="np")
    plt.legend(loc='upper left', bbox_to_anchor=(0.05, 0.95))
    plt.ylabel("time(s)")
    plt.title("Method compare")

    plt.subplot(212)
    plt.plot(X, Htime, color="g", linestyle="-", marker="s", linewidth=1,label="np")
    plt.legend(loc='upper left', bbox_to_anchor=(0.05, 0.95))
    plt.xlabel("size(px)")
    plt.ylabel("time(s)")

    plt.show()

结果：

Compare.png

结论：

由图可见，随着输入矩阵尺寸的增大，传统for循环的方法耗时逐渐离谱，而改进方法并未增减太多时间消耗（对人感觉而言）。

---

拓展：

同理可以把此方法推广到 N 维矩阵中应用 。

---

注意：

文章为本人在简书平台原创，盗用必究！