DistanceTransform

DistanceTransform - Compute distance transform of a binary matrix.

  This function assigns to every point (x,y) in a binary matrix B the distance
  to the nearest point in B with value 1. It uses the Euclidean metric for
  computing distances. The code is based on the algorithm by Meijster et al.
  (2000).

  USAGE

    dt = DistanceTransform(b)

    b              binary matrix

0001 function dt = DistanceTransform(b)
0002 
0003 %DistanceTransform - Compute distance transform of a binary matrix.
0004 %
0005 %  This function assigns to every point (x,y) in a binary matrix B the distance
0006 %  to the nearest point in B with value 1. It uses the Euclidean metric for
0007 %  computing distances. The code is based on the algorithm by Meijster et al.
0008 %  (2000).
0009 %
0010 %  USAGE
0011 %
0012 %    dt = DistanceTransform(b)
0013 %
0014 %    b              binary matrix
0015 %
0016 
0017 % Copyright (C) 2014-2018 by Michaël Zugaro
0018 %
0019 % This program is free software; you can redistribute it and/or modify
0020 % it under the terms of the GNU General Public License as published by
0021 % the Free Software Foundation; either version 3 of the License, or
0022 % (at your option) any later version.
0023 
0024 % The code is not self explanatory, but is a direct transcription of the algorithms provided
0025 % in Meijster et al. (2000)
0026 
0027 [m,n] = size(b);
0028 g = nan(m,n);
0029 
0030 for x = 1:n,
0031     % Scan 1
0032     if b(1,x),
0033         g(1,x) = 0;
0034     else
0035         g(1,x) = m + n;
0036     end
0037     for y = 2:m,
0038         if b(y,x),
0039             g(y,x) = 0;
0040         else
0041             g(y,x) = 1 + g(y-1,x);
0042         end
0043     end
0044     % Scan 2
0045     for y = m-1:-1:1,
0046         if g(y+1,x) < g(y,x),
0047             g(y,x) = 1 + g(y+1,x);
0048         end
0049     end
0050 end
0051 
0052 for y = 1:m,
0053     q = 1;
0054     s(1) = 1;
0055     t(1) = 1;
0056     for u = 2:n,
0057         % Scan 3
0058         while true,
0059             x = t(q);
0060             i = s(q);
0061             f1 = (x-i)^2 + g(y,i)^2;
0062             i = u;
0063             f2 = (x-i)^2 + g(y,i)^2;
0064             if f1 <= f2, break; end
0065             q = q - 1;
0066             if q < 1, break; end
0067         end
0068         if q < 1,
0069             q = 1;
0070             s(1) = u;
0071         else
0072             i = s(q);
0073             a = (u^2-i^2+g(y,u)^2-g(y,i)^2);
0074             b = 2*(u-i);
0075             w = 1 + floor(a/b);
0076             if w <= n,
0077                 q = q + 1;
0078                 s(q) = u;
0079                 t(q) = w;
0080             end
0081         end
0082     end
0083     for u = n:-1:1,
0084         % Scan 4
0085         x = u;
0086         i = s(q);
0087         dt(y,u) = sqrt((x-i)^2 + g(y,i)^2);
0088         if u == t(q),
0089             q = q - 1;
0090         end
0091     end
0092 end