Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Matlab / Octave bwdist() in Python or C

Does anyone know of a Python replacement for Matlab / Octave bwdist() function? This function returns Euclidian distance of each cell to the closest non-zero cell for a given matrix. I saw an Octave C implementation, a pure Matlab implementation, and I was wondering if anyone had to implement this in ANSI C (which doesn't include any Matlab / Octave headers, so I can integrate from Python easily) or in pure Python.

Both links I mentioned are below:

C++

Matlab M-File

As a test, a Matlab code / output looks something like this:

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

I tested a recommended distance_transform_edt call in Python, which gave this result:

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.]]

This result does not seem to match the Octave / Matlab output.

like image 259
BBSysDyn Avatar asked Mar 10 '11 13:03

BBSysDyn


3 Answers

While Matlab bwdist returns distances to the closest non-zero cell, Python distance_transform_edt returns distances “to the closest background element”. SciPy documentation is not clear about what it considers to be the “background”, there is some type conversion machinery behind it; in practice 0 is the background, non-zero is the foreground.

So if we have matrix 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]))

then to calculate the same result we need to replaces ones with zeros and zeros with ones, e.g. consider matrix 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]])

In this case scipy.ndimage.morphology.distance_transform_edt gives the expected results:

>>> 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.        ]])
like image 61
sastanin Avatar answered Oct 20 '22 06:10

sastanin


Does scipy.ndimage.morphology.distance_transform_edt meet your needs?

like image 22
asthasr Avatar answered Oct 20 '22 07:10

asthasr


No need to do the 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.        ]])
    
like image 2
luoenzhen Avatar answered Oct 20 '22 08:10

luoenzhen