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)
Get
the MATLAB code
Published with MATLAB® 7.13
This technique will work most of the time, but will occasionally fail. The basic issue is that a == b is only true if a is identical to b down to the last bit. round(a*32)/32 == round(b*32)/32 will be true for most values of a and b that are close. It fails however, for a=1/64 – 0.0000001 and b = 1/64; round(a*32)/32 = 0 whereas round(b*32)/32 = 0.03125 (i.e. 1/32).
Hence, you have replaced a bug which will happen fairly often with a bug that will happen only when you least expect it. If you are using single precision, this will happen more than one time in a million. In double precision the odds are much longer and perhaps acceptable to most people.
The following untested code should be more robust:
D = ( D – min(D(:)) ) / max(D(:))
D ( D < .00001 ) = 0
Ken—I don’t agree. Relative tolerance tests (like your suggestions) and absolute tolerances (like I used in this post) both have their places, and this is an appropriate use of an absolute tolerance test. I need to distinguish between a distance value corresponding to the shortest path and a distance corresponding to a path that is one step longer. The difference in distances could be 2.0 – sqrt(2). So really I could have gotten away with a much higher tolerance than I used above. However, it is not appropriate to scale the tolerance even higher according to the magnitude of D. Besides, the geodesic distance algorithm works by incremental additions of 1.0 or sqrt(2.0), so if the distances were huge enough to cause the comparison failure you are worried about, the distance transform step itself would have failed. But since single precision is good enough to handle distance transform computations for something like a 5 million x 5 million image (200 terabytes or so), I don’t think it’s worth worrying about.