%OPTICFLOW   Computes the optic flow of a 2D/3D time series
%  Image_in is a time series of 2D/3D dip_images, ordering x,y,t/ x,y,z,t
%  stored in one spatiotemporal dip_image 3D/4D.
%  The flow is computed via; dI/dt =0
%  v = -(\nabla I \nabla I^t)^{-1} \nabla I I_time
%
% SYNOPSIS:
%  out = opticflow(image_in, SpatialSig, TimeSig, TensorSig, Frame, Method)
%
% PARAMETERS:
%  SpatialSig   sigma of the spatial derivatives
%  TimeSig      sigma for the time derivatives
%  TensorSig    smoothing of the matrix of derivatives
%  Frame        frame number at which to compute the optic flow, otherwise
%               for all frames the flow is computed at once; a 3*TensorSig
%               window to both sides of the current time slice is used to
%               compute the flow
%  Method       'direct', fast but high memory usage
%               'sequential' slow but lower memory usage
%               the first and last 3*TensorSig frames are not computed
%
% DEFAULTS:
%  SpatialSig = 1
%  TimeSig    = 1
%  TensorSig  = 4
%  Frame      = -1 (all frames)
%  Method     = 'direct'
%
%  A 2D/3D vector images is returned with the x-coordinate
%  of the velocity field in the first component out{1}, etc.
%
% LITERATURE:
% Lucas, B.D. and Kanade, T.,
% An Iterative Image Registration Technique with an Application to Stereo Vision
% Proc. Image Understanding Workshop, p. 121-130, 1981

% (C) Copyright 1999-2007               Pattern Recognition Group
%     All rights reserved               Faculty of Applied Physics
%                                       Delft University of Technology
%                                       Lorentzweg 1
%                                       2628 CJ Delft
%                                       The Netherlands
%
% Bernd Rieger, Apr 2001.
% Bernd, Aug 2001, 6sigma window for time slice computation
%                  option to compute flow sequential/direct
% Bernd, March 2002, changed computation of slices to use less time
% Bernd, Jan 2004, removed 4D dimension limit, KvW implemented higher dim inverse
% MvG, Jun 2005, added undocumented feature: for (x,y,t) [3D] input only,
%                a confidence measure is optionally returned.
% 10 October 2007: Calling DERIVATIVE instead of DIP_GAUSS directly. (CL)

function [out,varargout] = opticflow(varargin)
d = struct('menu','Analysis',...
   'display','Optic flow',...
   'inparams',struct('name',       {'in','sigma','stime','stensor', 't0','method'},...
            'description',{'Input image','Spatial Sigma','Temporal Sigma',...
                         'Tensor Sigma','Time frame','Method'},...
            'type',       {'image','array','array','array','array','option'},...
            'dim_check',  {0, 0,   0,   -1,0,0},...
            'range_check',{[],'R+','R+','R+','Z',{'direct','sequential'}},...
            'required',   {1, 0,   0,   0,0,0},...
            'default',    {'a',1,   1,   4,-1,'direct'}...
             ),...
    'outparams',struct('name',{'out'},...
                      'description',{'Output image'},...
                      'type',{'image'}...
                      )...
          );
if nargin == 1
   s = varargin{1};
   if ischar(s) & strcmp(s,'DIP_GetParamList')
      out = d;
      return
   end
end
try
   [in, sigma, stime, stensor, t0, method] = getparams(d,varargin{:});
catch
   if ~isempty(paramerror)
      error(paramerror)
   else
      error(firsterr)
   end
end
di = ndims(in); %dimesions of the input image

%code is quite unreadable due to ND support and slice computation

if t0 <=0
   s = whos('in');
   me = s.bytes*(di-1);
   fprintf('Memory for output vector image +- %d MB\n',round(me/2^20));
   if round(me/2^20)>500
      ant=input('Are you sure the 4D image has reasonable size?[y/n]','s');
      if strcmp(ant,'n')
         return
      end
   end
   in2=in;
end
if strcmp(method,'sequential') & t0>0
   method = 'direct';
end
if strcmp(method,'sequential')
   fprintf('Sequential computation, get some coffee...\n');
   out = opt_step(in, sigma, stime, stensor);
   return;
end
if t0>=0
   sT = stensor(end);
   if t0 < 3*sT | t0 > size(in,di)-3*sT
      error(['You request a time slice that cannot be calculated' ...
      ' as it is within  3 sigmar tensor form the image borders, 99% confidence.']);
   end
   % compute 6sigma window indices and use this as input
   s1 = [num2str(t0) '-3*sT:3*sT+' num2str(t0)];
   s = ['in2 = in(' repmat(':,',1,di-1) s1 ');'];
   eval(s);
end


sigma = [repmat(sigma,1,di-1) stime];

% spatial derivative, only di-1 -> cannot use gradient
Ix = dip_image('array',di-1);

for ii =1:di-1
   order=zeros(1,di);
   order(ii)=1;
   Ix{ii} = derivative(in2,sigma,order);
end

% time derivative
It = derivative(in2,sigma,[repmat(0,1,di-1) 1]);

%mem_usage/2^20
A = smooth(Ix*Ix',stensor);
b = smooth(Ix.* It, stensor);

out = -1*inv(A)*b;

if (nargout>1) & (di==3)
   It2 = gaussf( It*It, stensor );
   varargout{1}=dip_symmetriceigensystem3(...
      A{1,1},A{1,2},b{1},A{2,2},b{2},It2,{'cylindrical'});
end

if t0 >= 0
   s2 = ['out(' repmat(':,',1,di-1) num2str(3*sT) ')'];
   out = squeeze(eval(s2));
end

function out = opt_step(in,sg,stime,sT)
sz = size(in);
di = ndims(in); %spatial + time

out = dip_image('array',di-1);
for jj = 1:di-1
   out{jj} = dip_image('zeros',sz);
end
s1 = repmat(':,',1,di-1);
s2 = ['out{jj}(' s1 'ii) = tmp{jj};'];

%for ii=0:sz(end)-1
sTt =sT(end);
fprintf('Computation from time slice %d to %d\n',3*sTt,sz(end)-1-3*sTt);
for ii=3*sTt:sz(end)-1-3*sTt
   tmp = opticflow(in,sg,stime,sT,ii,'direct');
   for jj = 1:di-1
      eval(s2)
   end
end
