%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% In this example, a set of gridable data (we call it a "matroid") from a
% paper by Farouki is turned into a surface. From 18 points we construct a
% smooth surface that <b>interpolates</b> the data. 
% 
% The matroid in this case is given by the dimensions (in this case 6x3) and
% then three columns of data (each 6x3=18 elements long). The first column
% represents the x coordinates, the second the y, and the third the z
% coordinates. Each is then turned into an ordinary matrix (one for the x
% coordinates, one for the y, etc.). 
% 
% The SVD is performed on the three matrices separately, and then each svd
% is fitted using the "flesh" procedure of Long and Long (Long, A. and
% C. Long. Surface Approximation and Interpolation Via Matrix SVD. The
% College Mathematics Journal, Vol. 32, #1, January, 2001:20-25.)
% 
% The interpolators are each functions, so we've got a continuous function
% of each coordinate, and they can be recombined to produce a continuous
% function in three space that interpolates the original Farouki data.
% 
% http://www.nku.edu/~longa/papers/CMJ020-0.PDF
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Farouki data
% Original: 6-20-96  svd3dfars.m
% 
% Modification for the addition of this documentation, to improve
% readability of the code, and to improve the stereo viewing
% 	(in loving memory of Clifford Long, by his son Andy): 
% 	Sat Oct 9 10:11:50 EDT 2010 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
far=[-3.0 -1.2  0.0
    -3.0  0.0  0.1
    -2.7  1.3 -0.1
    -2.0 -1.0 -0.2
    -1.7  0.2  0.3
    -1.5  1.6  0.0
    -1.1 -1.5  1.0
    -0.6 -0.1  1.1
    -0.8  1.2  0.3
    -0.4 -0.8 -0.3
    0.0  0.1  0.1
    0.2  1.0  0.5
    0.8 -1.0  0.3
    1.5 -0.3  0.4
    2.2  1.5  0.6
    2.3 -1.7  0.2
    3.0 -0.1  0.3
    3.3  1.2  0.5 ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% number of data points in each direction: rows, columns
xdata = 3		
ydata = 6
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% x-coordinates are ax
% y-coordinates are ay
% z-coordinates are az
ax = far(:,1);
ay = far(:,2);
az = far(:,3);
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Now turn them into matrices, with the proper shape:
% 
axx=reshape(ax,xdata,ydata)';
ayy=reshape(ay,xdata,ydata)';
azz=reshape(az,xdata,ydata)';
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% number of rows & columns in surface representation:
% 
k=10;
xcoords=k*(xdata-1)+1;
ycoords=k*(ydata-1)+1;
xx=linspace(1,xdata,xdata)';
yy=linspace(1,ydata,ydata)';
xxx=linspace(1,xdata,xcoords)';
yyy=linspace(1,ydata,ycoords)';
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
[u,s,v]=svd(axx);
s = diag(s);
z = zeros(ycoords,xcoords);
for j=1:3
  spx = spline(xx,v(:,j), xxx);
  spy = spline(yy,u(:,j), yyy);
  z = z + s(j)*spy*spx';
end
splx=z;
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
[u,s,v]=svd(ayy);
s = diag(s);
z = zeros(ycoords,xcoords);
for j=1:3
  spx = spline(xx,v(:,j), xxx);
  spy = spline(yy,u(:,j), yyy);
  z = z + s(j)*spy*spx';
end
sply=z;
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
[u,s,v]=svd(azz);
s = diag(s);
z = zeros(ycoords,xcoords);
for j=1:3
  spx = spline(xx,v(:,j), xxx);
  spy = spline(yy,u(:,j), yyy);
  z = z + s(j)*spy*spx';
end
splz=z;
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This is a plot of the raw data:
% 
subplot(1,1,1)
title("Plot of the Raw Farouki Rata");
plot3(axx,ayy,azz,'o');
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Plot in stereo (two different views):
% 
title("Farouki Data (stereo -- wall-eyed -- viewing)");
subplot(2,2,1)
view(150,30)
subplot(2,2,2)
view(145,30)
subplot(2,2,3)
view(50,20)
subplot(2,2,4)
view(45,20)
% clf
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
subplot(1,1,1)
title("The \"Farouki Surface\"");
surf(splx,sply,splz,splz);
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Plot in stereo (two different views):
% 
subplot(2,2,1)
title("Stereo Views (wall-eyed) of a \"Farouki Surface\"");
view(30,30)
subplot(2,2,2)
view(25,30)
subplot(2,2,3)
view(135,20)
subplot(2,2,4)
view(130,20)
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


