function [mindist,route]=dtw(A,B,r,flag);

% Dynamic Time Warping
%  [mindist,route] = dtw(A,B,r,flag)
%
%  Compute distance btwn two sequences 
%  of feature vectors, A and B.
%  (based on Sakoe and Chibas' paper)
%   
%   A, B - sequences of feature vectors
%      r - adjustment window size
%   flag - warping function display flag
% 

%% Initialization

[m,n]=size(A);
I=n;
[m,n]=size(B);
J=n;
g=zeros(I,J);    % minimal cost
ind=zeros(I,J);  % array of trace
g(1,1)= 2*mydist(A,B,1,1);
LargeNumber=realmax;
w=sqrt(-1);

%%compute local distances

imax=min(J+r, 2*J-1);
imin=max(J-r, floor(J/2)+1);
if(I>imax | I<imin) 
    %%error('Endpoint does not fall in adjustement area!!!');
    route=[];
    mindist=10000;
else

  for j=2:J
    imax=min([j+r,2*j-1,I]);
    imin=max([j-r,floor(j/2)+1]);

    for i=imin:imax
       if(inrange(i-2,j-1,I,J,r))
          path(1)=g(i-2,j-1)+2*mydist(A,B,i-1,j)+mydist(A,B,i,j);
       else
          path(1)= LargeNumber;
       end
       if(inrange(i-1,j-1,I,J,r))
          path(2)= g(i-1,j-1)+2*mydist(A,B,i,j);
       else
          path(2) = LargeNumber;
       end
       if(inrange(i-1,j-2,I,J,r))
          path(3)=g(i-1,j-2)+2*mydist(A,B,i,j-1)+mydist(A,B,i,j);
       else
          path(3)=LargeNumber;
       end
       [g(i,j),index]=min(path);
       if index==1
          ind(i,j)=(i-2)+(j-1)*w;
       elseif index ==2
          ind(i,j)=(i-1)+(j-1)*w;
       elseif index ==3
          ind(i,j)=(i-1)+(j-2)*w;
       end
    end
  end

  %% minimal distance btwn two sequences
  mindist=g(I,J)/(I+J);

  %% determine the best route (warping function)
  route=[];
  route(1)=I+J*w;
  route(2)=ind(I,J);
  x=real(route(2));y=imag(route(2));
  index=2;
  while(route(index) ~= 1+w)
   index=index+1;
   route(index)=ind(x,y);
   x=real(route(index));y=imag(route(index));
  end
  route=fliplr(route);

  %% if display flag is set, display computed
  %% warping function
  clf;
  if (flag == 1)
   axis([0,I+2,0,J+2]);grid;
   hold
   %% plot adjustment window
   xx=0:I+2;
   plot(xx,xx-r);
   plot(xx,xx+r);
   %% plot allowable points
   pts=ind;
   plot(1,1,'+');
   for j=2:J
     for i=1:I
       if(pts(i,j) ~= 0)
           plot(i,j,'+');
       end
     end
   end
   %plot warping function
   plot(real(route),imag(route),'r');
  end
end
return


                                                    