python - python中的测地距离变换

Question

在 python 中，scipy.ndimage.morphology 模块中有一个 distance_transform_edt 函数。我将它应用于一个简单的案例，以计算与屏蔽 numpy 数组中单个单元格的距离。

然而，该函数移除了数组的掩码，并按预期计算了每个具有非空值的单元格与具有空值的引用单元格的欧几里得距离。

下面是我在 my blog post 中给出的例子:

%pylab
from scipy.ndimage.morphology import distance_transform_edt
l = 100
x, y = np.indices((l, l))
center1 = (50, 20)
center2 = (28, 24)
center3 = (30, 50)
center4 = (60,48)
radius1, radius2, radius3, radius4 = 15, 12, 19, 12
circle1 = (x - center1[0])**2 + (y - center1[1])**2 < radius1**2
circle2 = (x - center2[0])**2 + (y - center2[1])**2 < radius2**2
circle3 = (x - center3[0])**2 + (y - center3[1])**2 < radius3**2
circle4 = (x - center4[0])**2 + (y - center4[1])**2 < radius4**2
# 3 circles
img = circle1 + circle2 + circle3 + circle4
mask = ~img.astype(bool)
img = img.astype(float)
m = ones_like(img)
m[center1] = 0
#imshow(distance_transform_edt(m), interpolation='nearest')
m = ma.masked_array(distance_transform_edt(m), mask)
imshow(m, interpolation='nearest')

但是我想计算考虑到数组的屏蔽元素的测地线距离变换。我不想沿着穿过屏蔽元素的直线计算欧几里得距离。

我使用 The Dijkstra 算法得到了我想要的结果。以下是我提出的实现方案:

def geodesic_distance_transform(m):
    mask = m.mask
    visit_mask = mask.copy() # mask visited cells
    m = m.filled(numpy.inf)
    m[m!=0] = numpy.inf
    distance_increments = numpy.asarray([sqrt(2), 1., sqrt(2), 1., 1., sqrt(2), 1., sqrt(2)])
    connectivity = [(i,j) for i in [-1, 0, 1] for j in [-1, 0, 1] if (not (i == j == 0))]
    cc = unravel_index(m.argmin(), m.shape) # current_cell
    while (~visit_mask).sum() > 0:
        neighbors = [tuple(e) for e in asarray(cc) - connectivity 
                     if not visit_mask[tuple(e)]]
        tentative_distance = [distance_increments[i] for i,e in enumerate(asarray(cc) - connectivity) 
                              if not visit_mask[tuple(e)]]
        for i,e in enumerate(neighbors):
            d = tentative_distance[i] + m[cc]
            if d < m[e]:
                m[e] = d
        visit_mask[cc] = True
        m_mask = ma.masked_array(m, visit_mask)
        cc = unravel_index(m_mask.argmin(), m.shape)
    return m

gdt = geodesic_distance_transform(m)
imshow(gdt, interpolation='nearest')
colorbar()

上面实现的函数运行良好，但对于我开发的需要多次计算测地线距离变换的应用程序来说太慢了。

下面是欧氏距离变换和测地线距离变换的时间基准:

%timeit distance_transform_edt(m)
1000 loops, best of 3: 1.07 ms per loop

%timeit geodesic_distance_transform(m)
1 loops, best of 3: 702 ms per loop

如何获得更快的测地距离变换？

Answer 1

首先，对一个非常清晰且写得很好的问题竖起大拇指。

有一个名为 scikit-fmm 的 Fast Marching 方法可以很好且快速地实现来解决此类问题。您可以在此处找到详细信息: http://pythonhosted.org//scikit-fmm/

安装它可能是最难的部分，但在带有 Conda 的 Windows 上它很容易，因为有用于 Py27 的 64 位 Conda 包: https://binstar.org/jmargeta/scikit-fmm

从那里开始，只需将掩码数组传递给它，就像您对自己的函数所做的那样。喜欢:

distance = skfmm.distance(m)

结果看起来相似，我认为甚至更好。您的方法(显然)在八个不同的方向上搜索，从而产生了一点“八角形”距离。

在我的机器上，scikit-fmm 实现比您的函数快 200 多倍。