Exploring shortest paths – part 16

Posted by Steve Eddins,

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:

1. Compute the distance transform for just the upper left block of foreground pixels.
2. Compute the distance transform for just the lower right block of foreground pixels.
3. Add the two distance transforms together.
4. The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.

(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

6 CommentsOldest to Newest

Sepideh replied on : 1 of 6

Hi,
Thanks for the example and explanation.

I guess bwdist will be helpful for my problem too.
I have a binary image which contains one pixel width skeleton of an image. I have got the skeleton using the bwmorph function.
As you may know basically in this binary image, I have ones on the skeleton and zeros anywhere else. how can I find the shortest ‘quasi-euclidean’ distance between two nonzero elements which are on the skeleton?
The path should be on the skeleton itself too.

I will really appreciate any help.

Thanks,

Mark Hayworth replied on : 2 of 6

It will be interesting to see your explanation of the new bwdistgeodesic() function. In the past I’ve solved this problem with dynamic programming (http://en.wikipedia.org/wiki/Dynamic_programming#Checkerboard) but it’s tricky and if you had a one-liner function call, that would be great.

Dynamic programming offers the ability to define a cost function so that the path cost depends on the gray levels of the image under the path, so for example you could use it to find the bright path or dark path along a blood vessel in a radiograph. It’s not simply and strictly distance-based but weighted by the underlying pixel values.

Steve replied on : 3 of 6

Mark—Sounds like the new function graydist may be of interest to you.

Steve replied on : 4 of 6

Sepideh—I will explore this idea in a later post. If you want to try to figure it out in the meantime, take a look at bwdistgeodesic.

You could try converting the 2D image to a mesh, represented as a sparse matrix, and then use matrix-vector multiply to do a breadth-first search. I use that to compute the Bacon-number for the IMDB graph: http://www.cise.ufl.edu/research/sparse/matrices/Pajek/IMDB.html , which computes the shortest path from Kevin Bacon to all other actors in the graph. I’m not sure how the pixel values themselves would be accounted for, though. Jeremy Kepner and John Gilbert have a new book out on linear algebra used for graph computations that would be well worth a look: http://www.ec-securehost.com/SIAM/SE22.html . The book’s at my office and I’m at home, but I think they look at the shortest path problem, if I remember.

Steve replied on : 6 of 6

Tim—Thanks for the book suggestion.

These postings are the author's and don't necessarily represent the opinions of MathWorks.