%This is an exmple of using a Thin Plate Spline (TPS) Radial Bais Function
%(RBF) to interpolate non-uniformly spaced data. In this case, the x,y,z,c
%point cloud data of a scan of the head of Radiohead's Thom Yorke, which
%was graciously released by the band. The non-uniform data is interpolated
%on a uniform grid and plotted using the surf function.
%In the interests of speed and memory usage, the training scheme only uses a
%portion of the data and the mesh grid is relatively sparse. To increase
%the accuracy, increase N_t (up a maximum of N_pix) and the increase the
%resolution of the plot grid, increase N_grid. The code is not currently
%optimized for large matrices so it may be very slow on your computer if
%you increase these values. The default values run in about a minute on a
%Core2Duo @ 2.16 GHz
load frame100XZC.mat
%this mat file has data from frame 100 of Radiohead's "House of Cards" data
%set. This data set is an irregularly spaced cloud of data points in form
%X=[x;y],z,c where c is a colour intensity. The data was culled from the
%noisy data by skipping any points with C<30.
%The full data set is available from
%http://code.google.com/creative/radiohead/
%and is licensed under the Creative Commons Attribution-Noncommercial-Share Alike
%3.0 License
N_pix=size(X,2);%number of points in data set
N_t=1000;%approx number of points to use for training
p=N_t/N_pix;%probability of any point going into training set
idx_t=rand(1,N_pix)<p;%logical indices of training set data
X_range=[-75 275;-15 310];%limits of X to plot at
%we're going to set an outer boundary and "pin" the z-values down at a
%number of points in a circle.
N_pin=8;%number of points to pin down at
X_cent=mean(X_range,2);%centre of circle
r_pin=min(diff(X_range,[],2))/2;%radius of circle
z_pin=-250*ones(1,N_pin);%z-value to pin at
c_pin=30*ones(1,N_pin);%colur value to pin at
theta_pin=linspace(0,2*pi,N_pin+1);%angles to set pins at
theta_pin(end)=[];
X_pin=X_cent*ones(1,N_pin)+r_pin*[cos(theta_pin); sin(theta_pin)];%[x;y] pin positions
X_t=[X(:,idx_t) X_pin];%X training points plus pinned points
z_t=[z(:,idx_t) z_pin];%z training points plus pinned points
c_t=[c(:,idx_t) c_pin];%colour training points plus pinned points
tic;%start timer
[a_z, Xc]=train_thin_plate_spline(X_t,z_t);%train thin plate spline network
[a_c, Xc]=train_thin_plate_spline(X_t,c_t);%note that the radial basis
%function centres (Xc) are the same in both cases because X_t are the same
t1=toc;%elapsed time
fprintf('Training time: %f s\n',t1)
N_grid=[100 100];%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
tic%start timer
z_mesh=sim_thin_plate_spline(X_mesh,Xc,a_z);%perform tps calculation on mesh points
z_mesh=reshape(z_mesh,N_grid(1),N_grid(2));%reshape into rectangular matrix
z_mesh(z_mesh<min(z_pin))=min(z_pin);%set a floor to data
c_mesh=sim_thin_plate_spline(X_mesh,Xc,a_c);%repeat for colour values
c_mesh=reshape(c_mesh,N_grid(1),N_grid(2));
c_mesh(c_mesh<min(c_pin))=min(c_pin);
t2=toc;%elapsed time
fprintf('Mesh Calculation Time: %f s\n',t2)
figure(1)
surf(x_mesh,y_mesh,z_mesh,c_mesh)
colormap(gray)
shading interp
axis equal
%
%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
