python - 计算周期性二维晶格中的成对空间距离

Question

我一直在徒劳地寻找一种解决方案来计算周期性/包裹边格子上离散的等距正方形站点之间的空间距离。例如在 9x9 的格子上设置相邻站点为 2 个单位:

m = 9
lattice = np.zeros((m,m))
for i in arange(0,6+1,3):
    for j in arange(0,6+1,3):
        lattice[i,j]=1

In  []: lattice
Out []: array([[ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
               [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

In  []: plt.imshow(lattice, origin='origin', cmap='binary', interpolation='none')

其中站点用索引注释，红色箭头表示与站点 0(2 个单位)的等距离，我如何创建一个数组 sdist 给我站点之间的最短成对距离，即 sdist [7,1]=2，sdist[2,6] = 2.8284。常用工具 scipy.spatial.distance.cdist 在这里似乎不适用。

Answer 1

由于 distance.cdist 接受任意度量，作为可调用提供，问题在于为您的度量编写一个函数。

如果它是点 p 和 q 之间的包裹距离，那就是

def wrapped_euclidean_points(p, q):
    diff = np.abs(p - q)
    return np.linalg.norm(np.minimum(diff, m - diff))

其中 m 是晶格的大小(不幸的是 cdist 不支持将额外的参数传递给距离函数，所以 m 必须取自更大的范围。)

但在您的情况下，您希望边长为 1 的正方形之间的距离最小。这意味着在计算包裹差异向量 np.minimum(diff, m - diff) 之后，我们可以减少它的每个分量都减 1，因为比方说，两个正方形点之间的最小 x 差比这些正方形中心之间的 x 差小 1。当然，这种减法不应使差异负。所以函数就变成了

def wrapped_euclidean_squares(p, q):
    diff = np.abs(p - q)
    return np.linalg.norm(np.clip(np.minimum(diff, m - diff) - 1, 0, None))

其中 clip 处理两个中心具有相同 x 坐标或相同 y 坐标的情况。

剩下的只有两行(不算from scipy.spatial import distance):

coords = np.vstack(np.nonzero(lattice)).T
dist = distance.cdist(coords, coords, metric=wrapped_euclidean_squares)

示例的输出:

[[ 0.          2.          2.          2.          2.82842712  2.82842712       2.          2.82842712  2.82842712]
 [ 2.          0.          2.          2.82842712  2.          2.82842712       2.82842712  2.          2.82842712]
 [ 2.          2.          0.          2.82842712  2.82842712  2.       2.82842712  2.82842712  2.        ]
 [ 2.          2.82842712  2.82842712  0.          2.          2.          2.       2.82842712  2.82842712]
 [ 2.82842712  2.          2.82842712  2.          0.          2.       2.82842712  2.          2.82842712]
 [ 2.82842712  2.82842712  2.          2.          2.          0.       2.82842712  2.82842712  2.        ]
 [ 2.          2.82842712  2.82842712  2.          2.82842712  2.82842712       0.          2.          2.        ]
 [ 2.82842712  2.          2.82842712  2.82842712  2.          2.82842712       2.          0.          2.        ]
 [ 2.82842712  2.82842712  2.          2.82842712  2.82842712  2.          2.       2.          0.        ]]