在Python或C中实现Matlab / Octave的bwdist()

Question

有没有人知道用Python替代Matlab或Octave中的bwdist()函数的方法？这个函数的作用是计算给定矩阵中每个单元格到最近的非零单元格的欧几里得距离。我看到过Octave的C语言实现和纯Matlab的实现，我在想有没有人用ANSI C（不包含任何Matlab或Octave的头文件，这样我可以方便地从Python中集成）或者纯Python实现过这个功能。

我提到的两个链接如下：

C++

Matlab M-File

作为测试，Matlab的代码和输出大概是这样的：

bw= [0   1   0   0   0;
     1   0   0   0   0;
     0   0   0   0   1;
     0   0   0   0   0;
     0   0   1   0   0]

D = bwdist(bw)

D =

   1.00000   0.00000   1.00000   2.00000   2.00000
   0.00000   1.00000   1.41421   1.41421   1.00000
   1.00000   1.41421   2.00000   1.00000   0.00000
   2.00000   1.41421   1.00000   1.41421   1.00000
   2.00000   1.00000   0.00000   1.00000   2.00000

我在Python中测试了一个推荐的distance_transform_edt调用，结果是这样的：

import numpy as np
from scipy import ndimage

a = np.array(([0,1,0,0,0],
              [1,0,0,0,0],
              [0,0,0,0,1],
              [0,0,0,0,0],
              [0,0,1,0,0]))

res = ndimage.distance_transform_edt(a)
print res

[[ 0.  1.  0.  0.  0.]
 [ 1.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.]]

这个结果似乎和Octave或Matlab的输出不一致。

Answer 1

不需要做1-a的部分

>>> distance_transform_edt(a==0)
    array([[ 1.        ,  0.        ,  1.        ,  2.        ,  2.        ],
           [ 0.        ,  1.        ,  1.41421356,  1.41421356,  1.        ],
           [ 1.        ,  1.41421356,  2.        ,  1.        ,  0.        ],
           [ 2.        ,  1.41421356,  1.        ,  1.41421356,  1.        ],
           [ 2.        ,  1.        ,  0.        ,  1.        ,  2.        ]])
    

Answer 2

你觉得 scipy.ndimage.morphology.distance_transform_edt 这个工具能满足你的需求吗？

Answer 3

在Matlab中，bwdist这个函数会返回到最近的非零单元格的距离，而在Python中，distance_transform_edt这个函数返回的是“到最近的背景元素的距离”。不过，SciPy的文档没有明确说明什么是“背景”，实际上它背后有一些类型转换的机制；简单来说，0被视为背景，而非零的值则是前景。

假设我们有一个矩阵 a：

>>> a = np.array(([0,1,0,0,0],
              [1,0,0,0,0],
              [0,0,0,0,1],
              [0,0,0,0,0],
              [0,0,1,0,0]))

如果我们想得到相同的结果，就需要把矩阵中的1换成0，把0换成1，比如考虑矩阵 1-a：

>>> a
array([[0, 1, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0],
       [0, 0, 1, 0, 0]])
>>> 1 - a
array([[1, 0, 1, 1, 1],
       [0, 1, 1, 1, 1],
       [1, 1, 1, 1, 0],
       [1, 1, 1, 1, 1],
       [1, 1, 0, 1, 1]])

在这种情况下，scipy.ndimage.morphology.distance_transform_edt会给出我们预期的结果：

>>> distance_transform_edt(1-a)
array([[ 1.        ,  0.        ,  1.        ,  2.        ,  2.        ],
       [ 0.        ,  1.        ,  1.41421356,  1.41421356,  1.        ],
       [ 1.        ,  1.41421356,  2.        ,  1.        ,  0.        ],
       [ 2.        ,  1.41421356,  1.        ,  1.41421356,  1.        ],
       [ 2.        ,  1.        ,  0.        ,  1.        ,  2.        ]])