%this is a simple example of using a thin plate spline radial basis
%function to fit a surface to a number of points
%
%Copyright (c) 2009, Travis Wiens
%All rights reserved.
%
%Redistribution and use in source and binary forms, with or without
%modification, are permitted provided that the following conditions are
%met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
%THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%POSSIBILITY OF SUCH DAMAGE.
%
% If you would like to request that this software be licensed under a less
% restrictive license (i.e. for commercial closed-source use) please
% contact Travis at travis.mlfx@nutaksas.com
%
N_t=20;%number of training points
N_dim=2;%number of dimensions in input
X_t=rand(N_dim,N_t);%training input
Y_t=rand(1,N_t);%training output
[a, Xc]=train_thin_plate_spline(X_t,Y_t);%train thin plate spline network
X_range=[0 1;0 1];%limits of X to plot at
N_grid=[256 256];%number of points to plot at in form [x y]
[x_mesh y_mesh]=meshgrid(linspace(X_range(1,1),X_range(1,2),N_grid(1)),linspace(X_range(2,1),X_range(2,2),N_grid(2)));
X_mesh=[reshape(x_mesh,1,[]);reshape(y_mesh,1,[])];%reshape into X=[x;y] matrix
Y_mesh=sim_thin_plate_spline(X_mesh,Xc,a);%perform tps calculation on mesh points
Y_mesh_r=reshape(Y_mesh,N_grid(1),N_grid(2));%reshape into rectangular matrix
figure(1)
surf(x_mesh,y_mesh,Y_mesh_r)
colormap(jet)
shading interp
hold on
plot3(X_t(1,:),X_t(2,:),Y_t,'*k')
hold off
