function [p,q,D] = dp_dtw(M)
% [p,q] = dp(M)
% Use dynamic programming to find a min-cost path through matrix M.
% Return state sequence in p,q
% 2003-03-15 dpwe@ee.columbia.edu
% Copyright (c) 2003 Dan Ellis <dpwe@ee.columbia.edu>
% released under GPL - see file COPYRIGHT
[r,c] = size(M);
% costs
D = zeros(r+1, c+1);
D(1,:) = NaN;
D(:,1) = NaN;
D(1,1) = 0;
D(2:(r+1), 2:(c+1)) = M;
% traceback
phi = zeros(r,c);
for i = 1:r;
for j = 1:c;
[dmax, tb] = min([D(i, j), D(i, j+1), D(i+1, j)]);
D(i+1,j+1) = D(i+1,j+1)+dmax;
phi(i,j) = tb;
end
end
% Traceback from top left
i = r;
j = c;
p = i;
q = j;
while i > 1 & j > 1
tb = phi(i,j);
if (tb == 1)
i = i-1;
j = j-1;
elseif (tb == 2)
i = i-1;
elseif (tb == 3)
j = j-1;
else
error;
end
p = [i,p];
q = [j,q];
end
p=[1,p];
q=[1,q];
% Strip off the edges of the D matrix before returning
D = D(2:(r+1),2:(c+1));