%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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));

%%%%%
%  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;
 
 
%
% okay now do it as sums of rank 1 matrices
%
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))
disp('Hit return for more modes')
pause;
for i=2:50
  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))
  pause;
  drawnow;
end
%
% and show for 100 and 200 modes
%
figure
i=100;
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))
figure
i=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))