# Exploring shortest paths – part 3

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:

- 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.
- Round to lower precision.
- The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.

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)

## 댓글

댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.