python中的测地距离变换

Posted 2023-02-23

【中文标题】python中的测地距离变换

【英文标题】：geodesic distance transform in python

【发布时间】：2015-03-27 02:15:30

【问题描述】：

在 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')

但是，我想计算考虑到数组的掩码元素的测地距离变换。我不想计算穿过蒙版元素的直线上的欧几里得距离。

我使用了 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

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

【问题讨论】：

distance_transform_edt 没有给出你想要的，但你的 m 在输入函数之前也永远不会被屏蔽。

【参考方案1】：

实现与geodesic_distance_transform 相同的结果的稍快（大约 10 倍）的实现：

def getMissingMask(slab):

    nan_mask=numpy.where(numpy.isnan(slab),1,0)
    if not hasattr(slab,'mask'):
        mask_mask=numpy.zeros(slab.shape)
    else:
        if slab.mask.size==1 and slab.mask==False:
            mask_mask=numpy.zeros(slab.shape)
        else:
            mask_mask=numpy.where(slab.mask,1,0)
    mask=numpy.where(mask_mask+nan_mask>0,1,0)

    return mask

def geodesic(img,seed):

    seedy,seedx=seed
    mask=getMissingMask(img)

    #----Call distance_transform_edt if no missing----
    if mask.sum()==0:
        slab=numpy.ones(img.shape)
        slab[seedy,seedx]=0
        return distance_transform_edt(slab)

    target=(1-mask).sum()
    dist=numpy.ones(img.shape)*numpy.inf
    dist[seedy,seedx]=0

    def expandDir(img,direction):
        if direction=='n':
            l1=img[0,:]
            img=numpy.roll(img,1,axis=0)
            img[0,:]==l1
        elif direction=='s':
            l1=img[-1,:]
            img=numpy.roll(img,-1,axis=0)
            img[-1,:]==l1
        elif direction=='e':
            l1=img[:,0]
            img=numpy.roll(img,1,axis=1)
            img[:,0]=l1
        elif direction=='w':
            l1=img[:,-1]
            img=numpy.roll(img,-1,axis=1)
            img[:,-1]==l1
        elif direction=='ne':
            img=expandDir(img,'n')
            img=expandDir(img,'e')
        elif direction=='nw':
            img=expandDir(img,'n')
            img=expandDir(img,'w')
        elif direction=='sw':
            img=expandDir(img,'s')
            img=expandDir(img,'w')
        elif direction=='se':
            img=expandDir(img,'s')
            img=expandDir(img,'e')

        return img

    def expandIter(img):
        sqrt2=numpy.sqrt(2)
        tmps=[]
        for dirii,dd in zip(['n','s','e','w','ne','nw','sw','se'],\
                [1,]*4+[sqrt2,]*4):
            tmpii=expandDir(img,dirii)+dd
            tmpii=numpy.minimum(tmpii,img)
            tmps.append(tmpii)
        img=reduce(lambda x,y:numpy.minimum(x,y),tmps)

        return img

    #----------------Iteratively expand----------------
    dist_old=dist
    while True:
        expand=expandIter(dist)
        dist=numpy.where(mask,dist,expand)
        nc=dist.size-len(numpy.where(dist==numpy.inf)[0])

        if nc>=target or numpy.all(dist_old==dist):
            break
        dist_old=dist

    return dist

另请注意，如果蒙版形成超过 1 个连接区域（例如，添加另一个不接触其他圆圈的圆圈），您的函数将陷入无限循环。

更新：

我发现了一个 Cython 实现的快速扫描方法in this notebook，它可以用来实现与scikit-fmm 相同的结果，并且速度可能相当。只需将二进制标志矩阵（以 1 作为可行点，inf 否则）作为 GDT() 函数的成本。

【参考方案2】：

首先，为一个非常清晰且写得很好的问题点赞。

有一个非常好的和快速的实现 Fast Marching 方法称为 scikit-fmm 来解决这类问题。您可以在这里找到详细信息： http://pythonhosted.org//scikit-fmm/

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

从那里开始，只需将您的掩码数组传递给它，就像您使用自己的函数一样。喜欢：

distance = skfmm.distance(m)

结果看起来很相似，而且我认为稍微好一点。您的方法（显然）在八个不同的方向上搜索，从而产生了一点“八角形”距离。

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

【讨论】：

谢谢你的方法。结果看起来非常好而且很快。这就是我想要的！但是，唯一的难点是安装。目前我没有编译成功，我更喜欢使用更多标准python库的解决方案来方便代码的分发。

我也使用了这种方法，但是距离计算与 scipy edt 函数相比存在差异。例如，如果我的网格空间为 0.5，并且我想计算一个点到所有周围点的距离，则对角线单元格在 scipy edt 版本中的距离为 0.70710678，在 scikit fmm 版本中为 0.85355339。为什么会这样？

【参考方案3】：

现在可从 Christoph Gohlke 获得 scikit-fmm 的 64 位 Windows 二进制文件。

http://www.lfd.uci.edu/~gohlke/pythonlibs/#scikit-fmm