Exploring shortest paths – part 2
Earlier this month I started exploring the question of computing the shortest path from one object in a binary image to another. I described the following basic algorithm:
- Compute the distance transform for just the upper left block of foreground pixels.
- Compute the distance transform for just the lower right block of foreground pixels.
- Add the two distance transforms together.
- The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.
This algorithm only works for path-based approximations to the Euclidean distance transform; it doesn't work for the Euclidean distance transform itself.
The Image Processing Toolbox function supports three path-based approximations to the Euclidean distance transform: 'cityblock', 'chessboard', and 'quasi-euclidean'.
Today I want to compare these three and discuss ambiguities associated with the idea of "shortest path" on an image.
Let's work with another small test image:
bw = logical([ ...
0 0 0 0 0 0 0
0 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 0 0 1 0
0 0 0 0 0 0 0 ]);One way to get two binary images, each containing one object from the original, is to use a label matrix:
L = bwlabel(bw); bw1 = (L == 1)
bw1 =
0 0 0 0 0 0 0
0 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 0 0 0 0
0 0 0 0 0 0 0
bw2 = (L == 2)
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 0 0 0 0 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 0
Step 1 of the algorithm is to compute the distance transform for the first object. I'll use the 'cityblock' distance transform approximation. In this approximation, only horizontal and vertical path segments are allowed. Diagonal segments are not allowed. The distance between a pixel and any of its N, E, S, or W neighbors is 1.
D1 = bwdist(bw1, 'cityblock')
D1 =
2 1 2 3 4 5 6
1 0 1 2 3 4 5
2 1 2 3 4 5 6
3 2 3 4 5 6 7
4 3 4 5 6 7 8
5 4 5 6 7 8 9
6 5 6 7 8 9 10
7 6 7 8 9 10 11
8 7 8 9 10 11 12
9 8 9 10 11 12 13
10 9 10 11 12 13 14
Step 2 is to compute the distance transform for the second object.
D2 = bwdist(bw2, 'cityblock')
D2 =
14 13 12 11 10 9 10
13 12 11 10 9 8 9
12 11 10 9 8 7 8
11 10 9 8 7 6 7
10 9 8 7 6 5 6
9 8 7 6 5 4 5
8 7 6 5 4 3 4
7 6 5 4 3 2 3
6 5 4 3 2 1 2
5 4 3 2 1 0 1
6 5 4 3 2 1 2
Step 3 is to add the distance transforms together.
D = D1 + D2
D =
16 14 14 14 14 14 16
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
14 12 12 12 12 12 14
16 14 14 14 14 14 16
The smallest value of D, 12, is the minimum path length (for paths consisting only of horizontal and vertical segments) from one object to the other.
Step 4 is to find the set of pixels in the regional minimum of D. This set represents the shortest path.
paths = imregionalmin(D);
To help visualize the path I'll use imoverlay from the MATLAB Central File Exchange. The code below will overlay the pixels on the shortest path in gray.
P = false(size(bw)); P = imoverlay(P, paths, [.5 .5 .5]); P = imoverlay(P, bw, [1 1 1]); imshow(P, 'InitialMagnification', 'fit')
Uh oh, why do we have a big rectangular block of gray instead of a single path? The answer is that there isn't a unique shortest path. There are many paths you could travel to get from one object to the other that all have the same path length, 12. Here are just a view of them:
subplot(2,2,1) imshow(P, 'InitialMagnification', 'fit') x = [2 6 6]; y = [2 2 10]; hold on plot(x, y, 'g', 'LineWidth', 2) hold off subplot(2,2,2) imshow(P, 'InitialMagnification', 'fit') x = [2 2 6]; y = [2 10 10]; hold on plot(x, y, 'g', 'LineWidth', 2); hold off subplot(2,2,3) imshow(P, 'InitialMagnification', 'fit') x = [2 3 3 4 4 5 5 6 6 6 6 6 6]; y = [2 2 3 3 4 4 5 5 6 7 8 9 10]; hold on plot(x, y, 'g', 'LineWidth', 2) hold off subplot(2,2,4) imshow(P, 'InitialMagnification', 'fit') x = [2 2 2 2 2 2 3 3 4 4 5 5 6]; y = [2 3 4 5 6 7 7 8 8 9 9 10 10]; hold on plot(x, y, 'g', 'LineWidth', 2) hold off
Next, let's look at the 'chessboard' distance. In this distance transform approximation, a path can consist of horizontal, vertical, and diagonal segments, and the path length between a pixel and any of its neighbors (N, NE, E, SE, S, SW, W, and NW) is 1.0.
D1 = bwdist(bw1, 'chessboard') D2 = bwdist(bw2, 'chessboard') D = D1 + D2
D1 =
1 1 1 2 3 4 5
1 0 1 2 3 4 5
1 1 1 2 3 4 5
2 2 2 2 3 4 5
3 3 3 3 3 4 5
4 4 4 4 4 4 5
5 5 5 5 5 5 5
6 6 6 6 6 6 6
7 7 7 7 7 7 7
8 8 8 8 8 8 8
9 9 9 9 9 9 9
D2 =
9 9 9 9 9 9 9
8 8 8 8 8 8 8
7 7 7 7 7 7 7
6 6 6 6 6 6 6
5 5 5 5 5 5 5
5 4 4 4 4 4 4
5 4 3 3 3 3 3
5 4 3 2 2 2 2
5 4 3 2 1 1 1
5 4 3 2 1 0 1
5 4 3 2 1 1 1
D =
10 10 10 11 12 13 14
9 8 9 10 11 12 13
8 8 8 9 10 11 12
8 8 8 8 9 10 11
8 8 8 8 8 9 10
9 8 8 8 8 8 9
10 9 8 8 8 8 8
11 10 9 8 8 8 8
12 11 10 9 8 8 8
13 12 11 10 9 8 9
14 13 12 11 10 10 10
You can see that the shortest path length for the 'chessboard' distance is 8 instead of 12. Let's look at the pixels on the various shortest paths.
paths = imregionalmin(D) P = false(size(bw)); P = imoverlay(P, paths, [.5 .5 .5]); P = imoverlay(P, bw, [1 1 1]); clf imshow(P, 'InitialMagnification', 'fit')
paths =
0 0 0 0 0 0 0
0 1 0 0 0 0 0
1 1 1 0 0 0 0
1 1 1 1 0 0 0
1 1 1 1 1 0 0
0 1 1 1 1 1 0
0 0 1 1 1 1 1
0 0 0 1 1 1 1
0 0 0 0 1 1 1
0 0 0 0 0 1 0
0 0 0 0 0 0 0
Now the set of gray pixels is not rectangular, but it still encompasses multiple possible shortest paths. And, surprisingly, it seems to include some pixels on which the path actually moves temporarily away from its destination.
Here are several possible paths, all of length 8 (according to the 'chessboard' distance transform).
subplot(2,2,1) imshow(P, 'InitialMagnification', 'fit') x = [2 3 4 5 6 6 6 6 6]; y = [2 3 4 5 6 7 8 9 10]; hold on plot(x, y, 'g', 'LineWidth', 2) hold off subplot(2,2,2) imshow(P, 'InitialMagnification', 'fit') x = [2 2 2 2 2 3 4 5 6]; y = [2 3 4 5 6 7 8 9 10]; hold on plot(x, y, 'g', 'LineWidth', 2); hold off subplot(2,2,3) imshow(P, 'InitialMagnification', 'fit') x = [2 1 2 1 2 3 4 5 6]; y = [2 3 4 5 6 7 8 9 10]; hold on plot(x, y, 'g', 'LineWidth', 2) hold off subplot(2,2,4) imshow(P, 'InitialMagnification', 'fit') x = [2 3 4 3 2 3 4 5 6]; y = [2 3 4 5 6 7 8 9 10]; hold on plot(x, y, 'g', 'LineWidth', 2) hold off
Finally, let's look at the 'quasi-euclidean' distance transform. In this variation, paths from one pixel to another can consist of horizontal, vertical, and diagonal segments. The distance from a pixel to one of its horizontal or vertical neighbors is 1, while the distance to one of its diagonal neighbors is sqrt(2).
Here's the computation again, this time specifying 'quasi-euclidean':
D1 = bwdist(bw1, 'quasi-euclidean'); D2 = bwdist(bw2, 'quasi-euclidean'); D = D1 + D2
D = 12.4853 11.6569 11.6569 12.2426 12.8284 13.4142 14.8284 11.0711 9.6569 10.2426 10.8284 11.4142 12.0000 13.4142 10.4853 9.6569 9.6569 10.2426 10.8284 11.4142 12.8284 10.4853 9.6569 9.6569 9.6569 10.2426 10.8284 12.2426 10.4853 9.6569 9.6569 9.6569 9.6569 10.2426 11.6569 11.0711 9.6569 9.6569 9.6569 9.6569 9.6569 11.0711 11.6569 10.2426 9.6569 9.6569 9.6569 9.6569 10.4853 12.2426 10.8284 10.2426 9.6569 9.6569 9.6569 10.4853 12.8284 11.4142 10.8284 10.2426 9.6569 9.6569 10.4853 13.4142 12.0000 11.4142 10.8284 10.2426 9.6569 11.0711 14.8284 13.4142 12.8284 12.2426 11.6569 11.6569 12.4853
You can see that the shortest path length is approximately 9.6569.
min(D(:))
ans =
9.6569
As before, let's find the set of pixels that belong to a shortest path.
paths = imregionalmin(D)
paths =
0 0 0 0 0 0 0
0 1 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 1 1 0
0 0 0 0 0 1 0
0 0 0 0 0 0 0
That looks like trouble! The shortest-path pixels don't actually connect the two objects, as you can see below:
P = false(size(bw)); P = imoverlay(P, paths, [.5 .5 .5]); P = imoverlay(P, bw, [1 1 1]); clf imshow(P, 'InitialMagnification', 'fit')
Where has our algorithm gone astray? That's the question I'll tackle next time.
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)
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.