function varargout = mctlib_CalcVel(varargin)
% function varargout = mctlib_CalcVel(varargin)
%
%
% Calculates the velocity for the smooth regularized velocity field using
% in MCT.
%
% Inputs:
% 'x0' - the resting points
% 'y0' - the points after the velocity field has been applied
% 'sigma' - controls the high pass filtering (large values make for a
%           smoother field (default = 1))
% 'beta' - the parameters calculated from mctlib_CalcParam
% 
% See also: mctlib_CalcParam, mctlib_GaussGram, mctlib_normalise
%
% Dr. A. I. Hanna (2007)
error(nargchk(0,inf,nargin));
x0 = [];
y0 = [];
beta = [];
sigma = 1;
lambda = 1;
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename, ' must come param/value in pairs.'])
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'x0'
            x0 = varargin{i+1};
        case 'y0'
            y0 = varargin{i+1};
        case 'sigma'
            sigma = varargin{i+1};  
        case 'beta'
            beta = varargin{i+1};
        case 'lambda'
            lambda = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(x0)
    error('Parameter X0 must be supplied');
end
if isempty(y0)
    error('Parameter Y0 must be supplied');
end
if isempty(beta)
    error('Parameter BETA must be supplied');
end
[N_x0, dim_x0] = size(x0);
[N_y0, dim_y0] = size(y0);
if (dim_x0 ~= dim_y0)
    error('Dimensions of input vectors must be the same');
end
%function V = mct_reg_interp(X, x0, beta, sigma)
%x0 = x0 - ones(N_x0, 1)*mean(x0);
%y0 = y0 - ones(N_y0, 1)*mean(y0);

G = mctlib_GaussGram('V1', x0, 'V2', y0, 'sigma', sigma);
%V = (beta*(G + lambda*eye(size(G))))';
V = (beta*G)';
varargout{1} = V;
return