% This batch procedure creates four figures illustrating process of
% fitting a spline surface through a set of points


load testpoints
x = x0; y = y0; z = z0;

vaz = 45; vel = 40;
axislim = [-2000 2000 -2000 2000 -2000 2000];

% Plotting 3D points as pins pinned into a sphere

figure(1); whitebg('white');

[sx,sy,sz] = sphere(15);
[az, el, r] = cart2sph(x,y,z);
R = min(min(r));
figure(1);
hold on;
plot3(x,y,z,'bo');
H = line([x R*cos(az).*cos(el)]',[y R*sin(az).*cos(el)]',[z R*sin(el)]');
for i = 1:length(H), set(H(i), 'Color', [0 0 0]); end; % black color
H = mesh(R*sx, R*sy, R*sz, 'b');
set(H,'EdgeColor',[0 0 0])
axis(axislim); view(vaz,vel); 
axis('square'); axis('off'); 
title('Points in spherical coordinate system');
drawnow; 
hold off;
disp('Press any key to continue'); pause;

% Plotting spherical coordinates in cartesian system, pin plot

figure(2); whitebg('white');

sPlotPin(az/pi*180, el/pi*180, r);
xtick = -180:45:180;
H = line([xtick;xtick], [-90;90]*ones(size(xtick)), zeros(2,size(xtick,2)) ...
   , 'LineStyle','-' );
for i = 1:length(H), set(H(i), 'Color', [0 0 0]); end; % black color
ytick = -90:30:90;
H = line([-180;180]*ones(size(ytick)), [ytick;ytick], zeros(2,size(ytick,2)) ...
   , 'LineStyle','-' );
for i = 1:length(H), set(H(i), 'Color', [0 0 0]); end; % black color
axis('off');
title('Spherical coordinate system unfolded');
xlabel('azimuth'); ylabel('elevation');
disp('Press any key to continue'); pause;

%%%%%%%%%%%%%%
if 0==1,  %just commenting out this portion of code
xlabel('azimuth [deg]');
ylabel('elevation [deg]');
zlabel('radius');
H = gca; 
axis([-180 180 -90 90 get(H,'ZLim') ]);
set(H, 'XTick', -180:90:180);
set(H, 'YTick', -90:45:90);
grid('on');

end; 


% Plotting spline fit of spherical coordinates in cartesian system

figure(4); whitebg('white');
figure(3); whitebg('white'); 

% sM = sFit( az, el, r, 32, 1/100*mean(mean(r)), 'plot' );
sM = sFit( az, el, r, 32, 1/100*mean(mean(r)) );
surf(sM);
axis('off'); axis('square');
title('Surface fit');
drawnow;


% Plotting spline fit 

figure(4); whitebg('white');

% hold on;
% [m,n] = size(sM);
% rfit = sValue(sM, (az+pi)/(2*pi), (el+pi/2)/pi);
% [x1,y1,z1] = sph2cart( az, el, rfit );
% plot3(x1,y1,z1,'o');  % plot points superimposed on the surface
sPlotSurf(sM);
axis(axislim); view(vaz,vel); 
axis('off'); axis('square');
title('Folded surface');
drawnow;
% hold ('off');