% Simple distance transform function using Chamfer metrics
% DT(x,y) with integer or fractional metric (Borgefors/Butt and Maragos) and 5x5 template, 101x101 cells
% Programd yields a symmetric (hexadecagon) surface about the sampled point (50,50)
% M J de Smith - 18/5/03 for MATLAB (as a template prog)

% clear command window and set program timer
home;tic;

% define arrays and dimensions
DT=ones(101,101)*999; DT(50,50)=0; % input image
LDM=zeros(9); DX=zeros(9); DY=zeros(9); % Local Distance Metric (LDM) and mask arrays
ydim = 99; xdim = 99;

% define chamfer values
a1=2.2062; a2=1.4141; a3=0.9866;

% forward scan
LDM(1)=a1; LDM(2)=a1; LDM(3)=a1; LDM(4)=a2; LDM(5)=a3; LDM(6)=a2; LDM(7)=a1; LDM(8)=a3; LDM(9)=0;
DX(1)=-2; DX(2)=-2; DX(3)=-1; DX(4)=-1; DX(5)=-1; DX(6)=-1; DX(7)=-1; DX(8)=0;  DX(9)=0;
DY(1)=-1; DY(2)=1;  DY(3)=-2; DY(4)=-1; DY(5)=0;  DY(6)=1;  DY(7)=2;  DY(8)=-1; DY(9)=0;

for i = 3:xdim
for j = 3:ydim
    d0=DT(i,j);
    for k = 1:9
        r=i+DX(k); c=j+DY(k); d=DT(r,c); d0=min(d+LDM(k),d0);
    end % k
    DT(i,j)=d0;
end % j
end % i

% backwards scan
LDM(1)=0; LDM(2)=a3; LDM(3)=a1; LDM(4)=a2; LDM(5)=a3; LDM(6)=a2; LDM(7)=a1; LDM(8)=a1; LDM(9)=a1;
DX(1)=0; DX(2)=0; DX(3)=1;  DX(4)=1;  DX(5)=1; DX(6)=1; DX(7)=1; DX(8)=2;  DX(9)=2; 
DY(1)=0; DY(2)=1; DY(3)=-2; DY(4)=-1; DY(5)=0; DY(6)=1; DY(7)=2; DY(8)=-1; DY(9)=1;

for i = xdim:-1:3
for j = ydim:-1:3
    d0=DT(i,j);
    for k = 1:9
        r=i+DX(k); c=j+DY(k); d=DT(r,c); d0=min(d+LDM(k),d0);
    end % k
    DT(i,j)=d0;
end % j
end % i
DT(100:101,:)=[];DT(:,100:101)=[];DT(1:2,:)=[];DT(:,1:2)=[];

% shaded contour plot
whitebg([0 .4 .6]); rect = [100, 100, 650, 650];
figure(1); set(1,'Position',rect); str = datestr(now,0); plotrangeX=xdim-2; plotrangeY=ydim-2;
contourf (DT,8); colormap jet; alpha (.3); title([' Generated on: ',str]); xlabel('x-position'); ylabel ('y-position'); axis ([1 plotrangeX 1 plotrangeY]);

% display elapsed time
toc