In part 2 of Exploring shortest paths, I noted a problem with using the 'quasi-euclidean' distance transform to find the shortest paths between two objects in a binary image. Specifically, our algorithm resulted in a big, unfilled gap between the two objects.

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  ]);
L = bwlabel(bw);
bw1 = (L == 1);
bw2 = (L == 2);

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

paths = imregionalmin(D);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

The pixels marked in gray should be the set of pixels that lie along a shortest path from the two objects in the original image. You can plainly, see, however, that there's a large gap.

To see why, let's examine more closely the values of D.

D

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

There appears to be an unbroken set of pixels between the two objects with value 9.6569. Actually, it turns out that the pixels do not all have the same value. For example, D(3,3) appear to have the same value D(4,4) but are actually slightly different.

D(3,3)

ans =

    9.6569

D(4,4)

ans =

    9.6569

D(3,3) - D(4,4)

ans =

 -9.5367e-007

The difference between the two is the corresponding single-precision relative floating-point precision:

eps(D(3,3))

ans =

  9.5367e-007

This floating-point round-off error comes into play because of the way that multiples of sqrt(2) are being added in different orders.

So we need to adjust our algorithm to account for floating-point round-off error. One way to do that is to round the distance transform values to some lower precision. For example, the code below rounds the distance transform to be a multiple of (1/32).

D = round(D*32)/32;

Now D(3,3) and D(4,4) are equal:

D(3,3) - D(4,4)

ans =

     0

And imregionalmin works as expected to extract the set of pixels belonging to shortest paths. (I'm using imoverlay from the MATLAB Central File Exchange.)

paths = imregionalmin(D);
P = false(size(bw));
P = imoverlay(P, paths, [.5 .5 .5]);
P = imoverlay(P, bw, [1 1 1]);
imshow(P, 'InitialMagnification', 'fit')

Here's our revised algorithm, including the new rounding step:

Next time I'll look into how to pick a particular path among the many shortest-path choices.

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