geodesic distance transform in python

Question

In python there is the distance_transform_edt function in the scipy.ndimage.morphology module. I applied it to a simple case, to compute the distance from a single cell in a masked numpy array.

However the function remove the mask of the array and compute, as expected, the Euclidean distance for each cell, with non null value, from the reference cell, with the null value.

Below is an example I gave in 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')

However I want to compute the geodesic distance transform that take into account the masked elements of the array. I do not want to compute the Euclidean distance along a straight line that go through masked elements.

I used The Dijkstra algorithm to obtain the result I want. Below is the implementation I proposed:

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()

The function implemented above works well but is too slow for the application I developed which needs to compute the geodesic distance transform several times.

Below is the time benchmark of the euclidean distance transform and the geodesic distance transform:

%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

How can I obtained a faster geodesic distance transform?

Answer 1

First of all, thumbs up for a very clear and well written question.

There is a very good and fast implementation of a Fast Marching method called scikit-fmm to solve this kind of problem. You can find the details here: http://pythonhosted.org//scikit-fmm/

Installing it might be the hardest part, but on Windows with Conda its easy, since there is 64bit Conda package for Py27: https://binstar.org/jmargeta/scikit-fmm

From there on, just pass your masked array to it, as you do with your own function. Like:

distance = skfmm.distance(m)

The results looks similar, and i think even slightly better. Your approach searches (apparently) in eight distinct directions resulting in a bit of a 'octagonal-shaped` distance.

On my machine the scikit-fmm implementation is over 200x faster then your function.