Blog reader Bogdan asked this question last week:

I am trying to find the length and location of the shortest path between objects. It seems like some application of bwdist should be able to do this. Any advice? Example:

1 0 1
0 1 0
0 0 0
0 1 0
1 0 1

Should yield

0 0 0
0 0 0
0 p 0
0 0 0
0 0 0

With p = path length = 1

I agree, some application of bwdist should do it. Actually, there are several other useful functions that can be applied to this task, including imregionalmin and bwdistgeodesic, a new function that just shipped in R2011b.

I've been tinkering with this problem off and on for a few days, and I think it might take a few posts to explore some of the interesting ideas and issues posed by this problem. Today I'll start with the basics of using bwdist and imregionalmin to identify a set of shortest paths from one object to another.

Let's use this simple test image to start our exploration.

bw = [ ...
    0 0 0 0 0 0 0 0 0 0 0 0
    0 1 1 0 0 0 0 0 0 0 0 0
    0 1 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 1 0
    0 0 0 0 0 0 0 0 0 1 1 0
    0 0 0 0 0 0 0 0 0 0 0 0   ]

bw =

     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     1     0     0     0     0     0     0     0     0     0
     0     1     1     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     1     1     0
     0     0     0     0     0     0     0     0     0     1     1     0
     0     0     0     0     0     0     0     0     0     0     0     0

We want to find the shortest path (or paths!) from the upper left block of foreground pixels to the lower right block. The key algorithmic idea is this:

(For reference, see P. Soille, Morphological Image Analysis: Principles and Applications, 2nd edition, Springer, 2003, section 7.3.4, "Application to minimal path detection," pp. 233-234.)

But we can't use the normal Euclidean distance transform because it doesn't correspond to a path traversal from one pixel center to another. Rather, we have to use one of the distance transform approximations supported by bwdist: 'chessboard', 'cityblock', and 'quasi-euclidean'.

Here's how to do it using 'quasi-euclidean'.

bw1 = [ ...
    0 0 0 0 0 0 0 0 0 0 0 0
    0 1 1 0 0 0 0 0 0 0 0 0
    0 1 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0   ];

bw2 = [ ...
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 1 0
    0 0 0 0 0 0 0 0 0 1 1 0
    0 0 0 0 0 0 0 0 0 0 0 0   ];

D1 = bwdist(bw1, 'quasi-euclidean');
D2 = bwdist(bw2, 'quasi-euclidean');

D_sum = D1 + D2;
D_sum(3:5, 3:10)

ans =

    7.8284    7.8284    7.8284    7.8284    7.8284    7.8284    8.4142    9.0000
    8.4142    7.8284    7.8284    7.8284    7.8284    7.8284    7.8284    8.4142
    9.0000    8.4142    7.8284    7.8284    7.8284    7.8284    7.8284    7.8284

You can see that the smallest value of D_sum is about 7.8284, and there's a set of pixels in-between the two objects that all have that value. The function imregionalmin identifies them all:

path_pixels = imregionalmin(D_sum)

path_pixels =

     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     1     1     1     1     1     0     0     0     0
     0     0     0     1     1     1     1     1     1     0     0     0
     0     0     0     0     1     1     1     1     1     1     0     0
     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0

In the quasi-euclidean approximation to the distance transform, the distance from a pixel to one of its horizontal or vertical neighbors is 1, and the distance from a pixel to one of its diagonal neighbors is sqrt(2). The result above says that the shortest such path between the two objects is 7.8284. Furthermore, it says that there is more than one shortest path.

I plan to explore this topic more fully this month. When using the quasi-euclidean distance transform, for example, there are floating-point round-off issues that have to be dealt with when you try to apply the technique above to larger images. Another issue is the effect of using different path-based distance transform approximations, such as 'cityblock' or 'chessboard'. Yet another issue is how to identify specific paths from the family of paths identified by imregionalmin. Along the way I'll experiment with different ways to visualize the results.

Finally, I'll explore how to use a new function, bwdistgeodesic, to find shortest paths when there are constraints on where the paths can go.

All the posts in this series

• the basic idea of finding shortest paths by adding two distance transforms together (part 1)
• the nonuniqueness of the shortest paths (part 2)
• handling floating-point round-off effects (part 3)
• using thinning to pick out a single path (part 4)
• using bwdistgeodesic to find shortest paths subject to constraint (part 5)

Get the MATLAB code

Published with MATLAB® 7.13