%RBFFIT Fits radial basis function (RBF) to the data.
%    RBFFIT(x,f,RBF) finds the coefficients of the thin-plate spline RBF that 
%    interpolates the given data.  The thin-plate spline is defined as
%                          phi(r) = (r^2)*log(r)
%
%   Inputs: 
%          x : the nodes arranged in a n-by-d matrix, where n is the number of
%              nodes and d is dimension the interpolation is supposed to be
%              perfomed in.  For example for 200 nodes in 2-D, x would be a 
%              200-by-2 matrix where each row corresponds to a node point.
%
%          f : the data values at the node points.  If f contains more than one
%              column, then an RBF will be fit to each column of f.  Thus, 
%              f should be of size n-by-k, where n is the number of samples and
%              k is the number of different functions.  For example, if you
%              want to fit four different functions sampled at 200 points in 
%              2-D, then f would be of size 200-by-4.
%
%  See also: rbfval.
function lam = rbffit(x,f)

[n,d] = size(x);
[m,k] = size(f);

%
% Compute the pairwise distances between all the nodes.
%

% Prepare the distance squared matrix.
A = zeros(n);
for j=1:d
   [xd1,xd2] = meshgrid(x(:,j));
   A = A + (xd1-xd2).^2;
end
clear xd1 xd2;

% Need to apply the TPS radial function to A.  However, the TPS is r^2*log(r)
% and will have numerical trouble when r=0.  It should be 0 when r=0, so we
% explicity make this the case.
id = 1:(n+1):n^2;    % indicies for the diagonal of A (i.e. where r=0).
A(id) = 1;           % make the distance 1 since log(1) = 0.
% Compute the TPS interpolation matrix.  Note that we are dealing with the
% distances squared.  The TPS kernel with distances squared reduces to 
% 0.5*A*log(A)
A = 0.5*A.*log(A);

% Add on the extra polynomial conditions.
ev = ones(n,1);
A = [[A ev x];[[ev x]' zeros(d+1)]];

% Solve for the RBF expansion coefficients
lam = A\[f;zeros(d+1,k)];