function varargout = mctlib_GaussGram(varargin)
% function varargout = mctlib_GaussGram(varargin)
%
% Calculates the Gaussian Gram matrix for the two vectors V1 and V2 where
% V1 is a NxD matrix and V2 is a MxD matrix and D is the dimension.
% The output G is a NxM matrix where each element G(i,j) = dot(V1(i,:), V2(j,:))
%
% Inputs:
% 'v1' - the first set of points
% 'v2' - the second set of points
% 'sigma' - controls the high pass filtering (large values make for a
%           smoother field (default = 1))
% 
% See also: mctlib_CalcParam, mctlibCalcVel, mctlib_normalise
%
% Dr. A. I. Hanna (2007)
error(nargchk(0,inf,nargin));
v1 = [];
v2 = [];
sigma = 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 'v1'
            v1 = varargin{i+1};
        case 'v2'
            v2 = varargin{i+1};
        case 'sigma'
            sigma = varargin{i+1};  
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(v1)
    error('Parameter V1 must be supplied');
end
if isempty(v2)
    error('Parameter V2 must be supplied');
end
[N_v1, dim_v1] = size(v1);
[N_v2, dim_v2] = size(v2);
if (dim_v1 ~= dim_v2)
    error('Dimensions of input vectors must be the same');
end
k=-2*sigma^2;
G=repmat(v1,[1 1 N_v2])-permute(repmat(v2,[1 1 N_v1]),[3 2 1]);
G=squeeze(sum(G.^2,2));
G=G/k;
G=exp(G);
G = G./(2*pi*sigma^2);
varargout{1} = G;
return