{"id":397,"date":"2011-12-02T13:00:55","date_gmt":"2011-12-02T18:00:55","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=397"},"modified":"2019-10-29T17:01:01","modified_gmt":"2019-10-29T21:01:01","slug":"exploring-shortest-paths-part-3","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2011\/12\/02\/exploring-shortest-paths-part-3\/","title":{"rendered":"Exploring shortest paths – part 3"},"content":{"rendered":"
\r\n

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

bw = logical([ ...<\/span>\r\n    0 0 0 0 0 0 0\r\n    0 1 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 0 0\r\n    0 0 0 0 0 1 0\r\n    0 0 0 0 0 0 0  ]);\r\nL = bwlabel(bw);\r\nbw1 = (L == 1);\r\nbw2 = (L == 2);\r\n\r\nD1 = bwdist(bw1, 'quasi-euclidean'<\/span>);\r\nD2 = bwdist(bw2, 'quasi-euclidean'<\/span>);\r\nD = D1 + D2;\r\n\r\npaths = imregionalmin(D);\r\nP = false(size(bw));\r\nP = imoverlay(P, paths, [.5 .5 .5]);\r\nP = imoverlay(P, bw, [1 1 1]);\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
The pixels marked in gray should be the set of pixels that lie along a shortest path from the two objects in the original\r\n      image. You can plainly, see, however, that there's a large gap.\r\n   <\/p>\r\n   
To see why, let's examine more closely the values of D<\/tt>.\r\n   <\/p>
D<\/pre>
\r\nD =\r\n\r\n   12.4853   11.6569   11.6569   12.2426   12.8284   13.4142   14.8284\r\n   11.0711    9.6569   10.2426   10.8284   11.4142   12.0000   13.4142\r\n   10.4853    9.6569    9.6569   10.2426   10.8284   11.4142   12.8284\r\n   10.4853    9.6569    9.6569    9.6569   10.2426   10.8284   12.2426\r\n   10.4853    9.6569    9.6569    9.6569    9.6569   10.2426   11.6569\r\n   11.0711    9.6569    9.6569    9.6569    9.6569    9.6569   11.0711\r\n   11.6569   10.2426    9.6569    9.6569    9.6569    9.6569   10.4853\r\n   12.2426   10.8284   10.2426    9.6569    9.6569    9.6569   10.4853\r\n   12.8284   11.4142   10.8284   10.2426    9.6569    9.6569   10.4853\r\n   13.4142   12.0000   11.4142   10.8284   10.2426    9.6569   11.0711\r\n   14.8284   13.4142   12.8284   12.2426   11.6569   11.6569   12.4853\r\n\r\n<\/pre>
There appears to be an unbroken set of pixels between the two objects with value 9.6569. Actually, it turns out that the pixels\r\n      do not all have the same value. For example, D(3,3)<\/tt> appear to have the same value D(4,4)<\/tt> but are actually slightly different.\r\n   <\/p>
D(3,3)<\/pre>
\r\nans =\r\n\r\n    9.6569\r\n\r\n<\/pre>
D(4,4)<\/pre>
\r\nans =\r\n\r\n    9.6569\r\n\r\n<\/pre>
D(3,3) - D(4,4)<\/pre>
\r\nans =\r\n\r\n -9.5367e-007\r\n\r\n<\/pre>
The difference between the two is the corresponding single-precision relative floating-point precision:<\/p>
eps(D(3,3))<\/pre>
\r\nans =\r\n\r\n  9.5367e-007\r\n\r\n<\/pre>
This floating-point round-off error comes into play because of the way that multiples of sqrt(2) are being added in different\r\n      orders.\r\n   <\/p>\r\n   
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\r\n      transform values to some lower precision. For example, the code below rounds the distance transform to be a multiple of (1\/32).\r\n   <\/p>
D = round(D*32)\/32;<\/pre>
Now D(3,3)<\/tt> and D(4,4)<\/tt> are equal:\r\n   <\/p>
D(3,3) - D(4,4)<\/pre>
\r\nans =\r\n\r\n     0\r\n\r\n<\/pre>
And imregionalmin<\/tt> works as expected to extract the set of pixels belonging to shortest paths. (I'm using imoverlay<\/a> from the MATLAB Central File Exchange.)\r\n   <\/p>
paths = imregionalmin(D);\r\nP = false(size(bw));\r\nP = imoverlay(P, paths, [.5 .5 .5]);\r\nP = imoverlay(P, bw, [1 1 1]);\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
Here's our revised algorithm, including the new rounding step:<\/p>\r\n   
\r\n      
    \r\n         
  1. Compute the distance transform for just the upper left block of foreground pixels.<\/li>\r\n         
  2. Compute the distance transform for just the lower right block of foreground pixels.<\/li>\r\n         
  3. Add the two distance transforms together.<\/li>\r\n         
  4. Round to lower precision.<\/li>\r\n         
  5. The pixels in the regional minimum of the result lie along one or more of the shortest paths from one block to the other.<\/li>\r\n      <\/ol>\r\n   <\/div>\r\n   
    Next time I'll look into how to pick a particular path among the many shortest-path choices.<\/p>\r\n
    All the posts in this series<\/h4>\r\n