%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   matlab script to calculate the SVD of an image
%   This script first reads in an image (a detail of an etching by Durer of a
%   magic square) and plots it.
%
%      It then takes the SVD of the  image and calculates the
%     spectrum of singular values.  Note that only the first 50 or so of the
%     singular values are large (and really only the first 4 or so).  Because the
%     high singular values are negligible, we can reconstruct much of the image just using
%    the first 50 "Empirical orthogonal Functions" (which are just the first 50 Eigenvectors
%    in V.  This script builds up the reconstruction one EOF at a time.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
%%%%%%
%  load and display the original image
%%%%%%
load detail
[m,n]=size(X);
imagesc(X); colormap(map); axis image; axis off;
set(gca,'fontweight','bold','fontsize',[14]);
title(sprintf('%s, %dx%d image',caption(1,:),m,n));
print -depsc2 Original.eps
%%%%%
%  now calculate the SVD of the image
%     and set VT=V transpose
%%%%%

[U,S,V]=svd(X,0);
VT=V';

%%%%%%
%  plot the spectrum of the Singular values sigma_1-sigma_m
%%%%%%

figure
semilogy(diag(S),'b-o');
set(gca,'fontsize',[16],'fontweight','bold');
title('Singular Values');
grid;
print -depsc2 Spectrum.eps
 
figure
i=1;
Xi=S(i,i)*U(:,i)*VT(i,:);
imagesc(Xi); colormap(map); axis image; axis off;
set(gca,'fontsize',[16],'fontweight','bold');
title(sprintf('EOF reconstruction with %d modes',i))
eval(sprintf('print -depsc2 EOF_%3.3d_modes.eps',i));
disp('Hit return for more modes')
pause;
for i=[10 20 30 40 50 100 200]
  Xi=U(:,1:i)*S(1:i,1:i)*VT(1:i,:);
  imagesc(Xi); colormap(map); axis image; axis off;
  set(gca,'fontsize',[16],'fontweight','bold')
  title(sprintf('EOF reconstruction with %d modes',i))
  eval(sprintf('print -depsc2 EOF_%3.3d_modes.eps',i));
  drawnow;
end