matthiasm@8: function [p,q,D,sc] = dpfast(M,C,T,G)
matthiasm@8: % [p,q,D,sc] = dpfast(M,C,T,G) 
matthiasm@8: %    Use dynamic programming to find a min-cost path through matrix M.
matthiasm@8: %    Return state sequence in p,q; full min cost matrix as D and 
matthiasm@8: %    local costs along best path in sc.
matthiasm@8: %    This version gives the same results as dp.m, but uses dpcore.mex
matthiasm@8: %    to run ~200x faster.
matthiasm@8: %    C is a step matrix, with rows (i step, j step, cost factor)
matthiasm@8: %    Default is [1 1 1.0;0 1 1.0;1 0 1.0];
matthiasm@8: %    Another good one is [1 1 1;1 0 1;0 1 1;1 2 2;2 1 2]
matthiasm@8: %    T selects traceback origin: 0 is to any edge; 1 is top right (default);
matthiasm@8: %    T > 1 finds path to min of anti-diagonal T points away from top-right.
matthiasm@8: %    Optional G defines length of 'gulleys' for T=0 mode; default 0.5
matthiasm@8: %    (i.e. accept path to only 50% of edge nearest top-right)
matthiasm@8: % 2003-04-04,2005-04-04 dpwe@ee.columbia.edu $Header: /Users/dpwe/projects/dtw/RCS/dpfast.m,v 1.6 2008/03/14 14:40:50 dpwe Exp $
matthiasm@8: 
matthiasm@8: % Copyright (c) 2003 Dan Ellis <dpwe@ee.columbia.edu>
matthiasm@8: % released under GPL - see file COPYRIGHT
matthiasm@8: 
matthiasm@8: if nargin < 2
matthiasm@8:   % Default step / cost matrix
matthiasm@8:   C = [1 1 1.0;0 1 1.0;1 0 1.0];
matthiasm@8: end
matthiasm@8: 
matthiasm@8: if nargin < 3
matthiasm@8:   % Default: path to top-right
matthiasm@8:   T = 1;
matthiasm@8: end
matthiasm@8: 
matthiasm@8: if nargin < 4
matthiasm@8:   % how big are gulleys?
matthiasm@8:   G = 0.5;  % half the extent
matthiasm@8: end
matthiasm@8: 
matthiasm@8: if sum(isnan(M(:)))>0
matthiasm@8:   error('dpwe:dpfast:NAN','Error: Cost matrix includes NaNs');
matthiasm@8: end
matthiasm@8: 
matthiasm@8: if min(M(:)) < 0
matthiasm@8:   disp('Warning: cost matrix includes negative values; results may not be what you expect');
matthiasm@8: end
matthiasm@8: 
matthiasm@8: [r,c] = size(M);
matthiasm@8: 
matthiasm@8: % Core cumulative cost calculation coded as mex
matthiasm@8: [D,phi] = dpcore(M,C);
matthiasm@8: 
matthiasm@8: p = [];
matthiasm@8: q = [];
matthiasm@8: 
matthiasm@8: %% Traceback from top left?
matthiasm@8: %i = r; 
matthiasm@8: %j = c;
matthiasm@8: 
matthiasm@8: if T == 0
matthiasm@8:   % Traceback from lowest cost "to edge" (gulleys)
matthiasm@8:   TE = D(r,:);
matthiasm@8:   RE = D(:,c);
matthiasm@8:   % eliminate points not in gulleys
matthiasm@8:   TE(1:round((1-G)*c)) = max(max(D));
matthiasm@8:   RE(1:round((1-G)*r)) = max(max(D));
matthiasm@8:   if (min(TE) < min(RE))
matthiasm@8:     i = r;
matthiasm@8:     j = max(find(TE==min(TE)));
matthiasm@8:   else
matthiasm@8:     i = max(find(RE==min(RE)));
matthiasm@8:     j = c;
matthiasm@8:   end
matthiasm@8: else
matthiasm@8:   % Traceback from min of antidiagonal
matthiasm@8:   %stepback = floor(0.1*c);
matthiasm@8:   stepback = T;
matthiasm@8:   slice = diag(fliplr(D),-(r-stepback));
matthiasm@8:   [mm,ii] = min(slice);
matthiasm@8:   i = r - stepback + ii;
matthiasm@8:   j = c + 1 - ii;
matthiasm@8: end
matthiasm@8: 
matthiasm@8: p=i;
matthiasm@8: q=j;
matthiasm@8: 
matthiasm@8: sc = M(p,q);
matthiasm@8: 
matthiasm@8: while i > 1 & j > 1
matthiasm@8: %  disp(['i=',num2str(i),' j=',num2str(j)]);
matthiasm@8:   tb = phi(i,j);
matthiasm@8:   i = i - C(tb,1);
matthiasm@8:   j = j - C(tb,2);
matthiasm@8:   p = [i,p];
matthiasm@8:   q = [j,q];
matthiasm@8:   sc = [M(i,j),sc];
matthiasm@8: end