{"id":396,"date":"2011-11-26T16:25:29","date_gmt":"2011-11-26T21:25:29","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2011\/11\/26\/exploring-shortest-paths-part-2\/"},"modified":"2019-10-29T17:00:11","modified_gmt":"2019-10-29T21:00:11","slug":"exploring-shortest-paths-part-2","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2011\/11\/26\/exploring-shortest-paths-part-2\/","title":{"rendered":"Exploring shortest paths – part 2"},"content":{"rendered":"
\r\n

Earlier this month<\/a> I started exploring the question of computing the shortest path from one object in a binary image to another. I described\r\n the following basic algorithm:\r\n <\/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. 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

    This algorithm only works for path-based approximations to the Euclidean distance transform; it doesn't work for the Euclidean\r\n distance transform itself.\r\n <\/p>\r\n

    The Image Processing Toolbox function supports three path-based approximations to the Euclidean distance transform: 'cityblock',\r\n 'chessboard', and 'quasi-euclidean'.\r\n <\/p>\r\n

    Today I want to compare these three and discuss ambiguities associated with the idea of \"shortest path\" on an image.<\/p>\r\n

    Let's work with another small test image:<\/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  ]);<\/pre>
    One way to get two binary images, each containing one object from the original, is to use a label matrix:<\/p>
    L = bwlabel(bw);\r\nbw1 = (L == 1)<\/pre>
    \r\nbw1 =\r\n\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     0     0\r\n     0     0     0     0     0     0     0\r\n\r\n<\/pre>
    bw2 = (L == 2)<\/pre>
    \r\nbw2 =\r\n\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     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\n\r\n<\/pre>
    Step 1 of the algorithm is to compute the distance transform for the first object. I'll use the 'cityblock' distance transform\r\n      approximation. In this approximation, only horizontal and vertical path segments are allowed. Diagonal segments are not allowed.\r\n      The distance between a pixel and any of its N, E, S, or W neighbors is 1.\r\n   <\/p>
    D1 = bwdist(bw1, 'cityblock'<\/span>)<\/pre>
    \r\nD1 =\r\n\r\n     2     1     2     3     4     5     6\r\n     1     0     1     2     3     4     5\r\n     2     1     2     3     4     5     6\r\n     3     2     3     4     5     6     7\r\n     4     3     4     5     6     7     8\r\n     5     4     5     6     7     8     9\r\n     6     5     6     7     8     9    10\r\n     7     6     7     8     9    10    11\r\n     8     7     8     9    10    11    12\r\n     9     8     9    10    11    12    13\r\n    10     9    10    11    12    13    14\r\n\r\n<\/pre>
    Step 2 is to compute the distance transform for the second object.<\/p>
    D2 = bwdist(bw2, 'cityblock'<\/span>)<\/pre>
    \r\nD2 =\r\n\r\n    14    13    12    11    10     9    10\r\n    13    12    11    10     9     8     9\r\n    12    11    10     9     8     7     8\r\n    11    10     9     8     7     6     7\r\n    10     9     8     7     6     5     6\r\n     9     8     7     6     5     4     5\r\n     8     7     6     5     4     3     4\r\n     7     6     5     4     3     2     3\r\n     6     5     4     3     2     1     2\r\n     5     4     3     2     1     0     1\r\n     6     5     4     3     2     1     2\r\n\r\n<\/pre>
    Step 3 is to add the distance transforms together.<\/p>
    D = D1 + D2<\/pre>
    \r\nD =\r\n\r\n    16    14    14    14    14    14    16\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    14    12    12    12    12    12    14\r\n    16    14    14    14    14    14    16\r\n\r\n<\/pre>
    The smallest value of D<\/tt>, 12, is the minimum path length (for paths consisting only of horizontal and vertical segments) from one object to the other.\r\n   <\/p>\r\n   
    Step 4 is to find the set of pixels in the regional minimum of D<\/tt>. This set represents the shortest path.\r\n   <\/p>
    paths = imregionalmin(D);<\/pre>
    To help visualize the path I'll use imoverlay<\/a> from the MATLAB Central File Exchange. The code below will overlay the pixels on the shortest path in gray.\r\n   <\/p>
    P = 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> 
    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\r\n      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\r\n      are just a view of them:\r\n   <\/p>
    subplot(2,2,1)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 6 6];\r\ny = [2 2 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2)\r\nhold off<\/span>\r\n\r\nsubplot(2,2,2)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 2 6];\r\ny = [2 10 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2);\r\nhold off<\/span>\r\n\r\nsubplot(2,2,3)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 3 3 4 4 5 5 6 6 6 6 6 6];\r\ny = [2 2 3 3 4 4 5 5 6 7 8 9 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2)\r\nhold off<\/span>\r\n\r\nsubplot(2,2,4)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 2 2 2 2 2 3 3 4 4 5 5 6];\r\ny = [2 3 4 5 6 7 7 8 8 9 9 10 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2)\r\nhold off<\/span><\/pre> 
    Next, let's look at the 'chessboard' distance. In this distance transform approximation, a path can consist of horizontal,\r\n      vertical, and diagonal segments, and the path length between a pixel and any of its neighbors (N, NE, E, SE, S, SW, W, and\r\n      NW) is 1.0.\r\n   <\/p>
    D1 = bwdist(bw1, 'chessboard'<\/span>)\r\nD2 = bwdist(bw2, 'chessboard'<\/span>)\r\nD = D1 + D2<\/pre>
    \r\nD1 =\r\n\r\n     1     1     1     2     3     4     5\r\n     1     0     1     2     3     4     5\r\n     1     1     1     2     3     4     5\r\n     2     2     2     2     3     4     5\r\n     3     3     3     3     3     4     5\r\n     4     4     4     4     4     4     5\r\n     5     5     5     5     5     5     5\r\n     6     6     6     6     6     6     6\r\n     7     7     7     7     7     7     7\r\n     8     8     8     8     8     8     8\r\n     9     9     9     9     9     9     9\r\n\r\n\r\nD2 =\r\n\r\n     9     9     9     9     9     9     9\r\n     8     8     8     8     8     8     8\r\n     7     7     7     7     7     7     7\r\n     6     6     6     6     6     6     6\r\n     5     5     5     5     5     5     5\r\n     5     4     4     4     4     4     4\r\n     5     4     3     3     3     3     3\r\n     5     4     3     2     2     2     2\r\n     5     4     3     2     1     1     1\r\n     5     4     3     2     1     0     1\r\n     5     4     3     2     1     1     1\r\n\r\n\r\nD =\r\n\r\n    10    10    10    11    12    13    14\r\n     9     8     9    10    11    12    13\r\n     8     8     8     9    10    11    12\r\n     8     8     8     8     9    10    11\r\n     8     8     8     8     8     9    10\r\n     9     8     8     8     8     8     9\r\n    10     9     8     8     8     8     8\r\n    11    10     9     8     8     8     8\r\n    12    11    10     9     8     8     8\r\n    13    12    11    10     9     8     9\r\n    14    13    12    11    10    10    10\r\n\r\n<\/pre>
    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\r\n      various shortest paths.\r\n   <\/p>
    paths = imregionalmin(D)\r\n\r\nP = false(size(bw));\r\nP = imoverlay(P, paths, [.5 .5 .5]);\r\nP = imoverlay(P, bw, [1 1 1]);\r\nclf\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre>
    \r\npaths =\r\n\r\n     0     0     0     0     0     0     0\r\n     0     1     0     0     0     0     0\r\n     1     1     1     0     0     0     0\r\n     1     1     1     1     0     0     0\r\n     1     1     1     1     1     0     0\r\n     0     1     1     1     1     1     0\r\n     0     0     1     1     1     1     1\r\n     0     0     0     1     1     1     1\r\n     0     0     0     0     1     1     1\r\n     0     0     0     0     0     1     0\r\n     0     0     0     0     0     0     0\r\n\r\n<\/pre> 
    Now the set of gray pixels is not rectangular, but it still encompasses multiple possible shortest paths. And, surprisingly,\r\n      it seems to include some pixels on which the path actually moves temporarily away from its destination.\r\n   <\/p>\r\n   
    Here are several possible paths, all of length 8 (according to the 'chessboard' distance transform).<\/p>
    subplot(2,2,1)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 3 4 5 6 6 6 6 6];\r\ny = [2 3 4 5 6 7 8 9 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2)\r\nhold off<\/span>\r\n\r\nsubplot(2,2,2)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 2 2 2 2 3 4 5 6];\r\ny = [2 3 4 5 6 7 8 9 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2);\r\nhold off<\/span>\r\n\r\nsubplot(2,2,3)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 1 2 1 2 3 4 5 6];\r\ny = [2 3 4 5 6 7 8 9 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2)\r\nhold off<\/span>\r\n\r\nsubplot(2,2,4)\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)\r\nx = [2 3 4 3 2 3 4 5 6];\r\ny = [2 3 4 5 6 7 8 9 10];\r\nhold on<\/span>\r\nplot(x, y, 'g'<\/span>, 'LineWidth'<\/span>, 2)\r\nhold off<\/span><\/pre> 
    Finally, let's look at the 'quasi-euclidean' distance transform. In this variation, paths from one pixel to another can consist\r\n      of horizontal, vertical, and diagonal segments. The distance from a pixel to one of its horizontal or vertical neighbors is\r\n      1, while the distance to one of its diagonal neighbors is sqrt(2).\r\n   <\/p>\r\n   
    Here's the computation again, this time specifying 'quasi-euclidean':<\/p>
    D1 = bwdist(bw1, 'quasi-euclidean'<\/span>);\r\nD2 = bwdist(bw2, 'quasi-euclidean'<\/span>);\r\nD = D1 + D2<\/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>
    You can see that the shortest path length is approximately 9.6569.<\/p>
    min(D(:))<\/pre>
    \r\nans =\r\n\r\n    9.6569\r\n\r\n<\/pre>
    As before, let's find the set of pixels that belong to a shortest path.<\/p>
    paths = imregionalmin(D)<\/pre>
    \r\npaths =\r\n\r\n     0     0     0     0     0     0     0\r\n     0     1     0     0     0     0     0\r\n     0     1     1     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     1     1     0\r\n     0     0     0     0     0     1     0\r\n     0     0     0     0     0     0     0\r\n\r\n<\/pre>
    That looks like trouble! The shortest-path pixels don't actually connect the two objects, as you can see below:<\/p>
    P = false(size(bw));\r\nP = imoverlay(P, paths, [.5 .5 .5]);\r\nP = imoverlay(P, bw, [1 1 1]);\r\nclf\r\nimshow(P, 'InitialMagnification'<\/span>, 'fit'<\/span>)<\/pre> 
    Where has our algorithm gone astray? That's the question I'll tackle next time.<\/p>\r\n\r\n
    All the posts in this series<\/h4>\r\n