How does the scipy distance_transform_edt function work?

I'm having trouble understanding how the Euclidean distance transform function works in Scipy. From what I understand, it is different than the Matlab function (bwdist). As an example, for the input:

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

The scipy.ndimage.distance_transform_edt function returns the same array:

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

But the matlab function returns this:

1.4142 1.0000 1.4142 2.2361 3.1623
1.0000 0 1.0000 2.0000 2.2361
1.4142 1.0000 1.4142 1.0000 1.4142
2.2361 2.0000 1.0000 0 1.0000
3.1623 2.2361 1.4142 1.0000 1.4142

which makes more sense, as it is returning the "distance" to the nearest one.

2 Answers

It is not clear from the docstring, but distance_transform_edt computes the distance from non-zero (i.e. non-background) points to the nearest zero (i.e. background) point.

For example:

In [42]: x
Out[42]:
array([[0, 0, 0, 0, 0, 1, 1, 1], [0, 1, 1, 1, 0, 1, 1, 1], [0, 1, 1, 1, 0, 1, 1, 1], [0, 0, 1, 1, 0, 0, 0, 1]])
In [43]: np.set_printoptions(precision=3) # Easier to read the result with fewer digits.
In [44]: distance_transform_edt(x)
Out[44]:
array([[ 0. , 0. , 0. , 0. , 0. , 1. , 2. , 3. ], [ 0. , 1. , 1. , 1. , 0. , 1. , 2. , 2.236], [ 0. , 1. , 1.414, 1. , 0. , 1. , 1. , 1.414], [ 0. , 0. , 1. , 1. , 0. , 0. , 0. , 1. ]])

You can get the equivalent of Matlab's bwdist(a) by applying distance_transform_edt() to np.logical_not(a) (i.e. invert the foreground and background):

In [71]: a
Out[71]:
array([[ 0., 0., 0., 0., 0.], [ 0., 1., 0., 0., 0.], [ 0., 0., 0., 0., 0.], [ 0., 0., 0., 1., 0.], [ 0., 0., 0., 0., 0.]])
In [72]: distance_transform_edt(np.logical_not(a))
Out[72]:
array([[ 1.414, 1. , 1.414, 2.236, 3.162], [ 1. , 0. , 1. , 2. , 2.236], [ 1.414, 1. , 1.414, 1. , 1.414], [ 2.236, 2. , 1. , 0. , 1. ], [ 3.162, 2.236, 1.414, 1. , 1.414]])

Warren has already explained how distance_transform_edt works. In your case,you could change sampling units along x and y

ndimage.distance_transform_edt(a)
array([[ 0., 0., 0., 0., 0.], [ 0., 1., 0., 0., 0.], [ 0., 0., 0., 0., 0.], [ 0., 0., 0., 1., 0.], [ 0., 0., 0., 0., 0.]])

But

>>> ndimage.distance_transform_edt(a, sampling=[2,2])
array([[ 0., 0., 0., 0., 0.], [ 0., 2., 0., 0., 0.], [ 0., 0., 0., 0., 0.], [ 0., 0., 0., 2., 0.], [ 0., 0., 0., 0., 0.]])

Or

ndimage.distance_transform_edt(a, sampling=[3,3])
array([[ 0., 0., 0., 0., 0.], [ 0., 3., 0., 0., 0.], [ 0., 0., 0., 0., 0.], [ 0., 0., 0., 3., 0.], [ 0., 0., 0., 0., 0.]])