%RBFVAL Evaluates the thin-plate spline (TPS) radial basis function interpolant 
%       specified by the RBF expansion coefficients lambda at the points t.
%
%   Inputs: 
%
%        lam : The rbf expansion coefficients.  Can be a column vector or
%              matrix.  If it is a matrix then separate evaluations are done for
%              each column.
%
%          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. 
%
%          t : The points at which to evaluate the RBF.  t should be arranged
%              as a m-by-d matrix, where m is the number of evaluation points 
%              and d is dimension of the interpolation.
%
%  See also: rbffit
function p = rbfval(lam,x,t)

% Get the sizes of the nodes and evaluation points.
[n,d] = size(x);
[m,d] = size(t);

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

% Need to apply the TPS radial function to B.  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.
locz = find(abs(B)<4*eps);
B(locz) = 1;
% Compute the TPS interpolation matrix.  Note that we are dealing with the
% distances squared.  The TPS kernel with distances squared reduces to 
% 0.5*B*log(B)
B = 0.5*B.*log(B);

% Add on the extra polynomial conditions.
ev = ones(m,1);
B = [B ev t];

% Evaluate the RBF
p = B*lam;