%IMAGE DECOMPOSITION AND RECONSTRUCTION USING SVD

clear all, close all, clc
% First consider decomposition and reconstruction  of images in
load gallery.mat
% the above command reads array b with nim=30 rows and (ny x nx)=24x16 
% columns, each row corresponds to one image
szb=size(b)
nim=szb(1);
ny=24;
nx=16;
colormap(gray(256));              %this assigns colors to the image
for i=1:nim
%   image(reshape(b(i,:),ny,nx)); % display all images
%   pause                         % hit any key to continue
end
bm=ones(nim,1)*mean(b);           % matrix with nim rows with the 
                                  % average image
image(reshape(bm(1,:),ny,nx));    % display the average image
pause 
b1=b-bm;             % subtract the average image from each original
[u,s,v]=svd(b1);                 % svd of the result
size(u), size(s), size(v)        % sizes of svd matrices
img=1                            % original image number img will be 
                                 % reconstructed
image(reshape(b(img,:),ny,nx))   % display the original image
figure(2)
colormap(gray(256));
figure(3)
colormap(gray(256));
for k=1:nim-1
   Neigenvectors=k               % the number of eigenvectors used in
                                 % reconstruction
   figure(2)
   b2=bm+u(:,1:k)*s(1:k,1:k)*v(:,1:k)'; % matrix of reconstructed images
   image(reshape(b2(img,:),ny,nx));     % display reconstructed image

% display the difference between the original and reconstructed 
% images: black-completely different, white-no difference

   figure(3)
   image(reshape(255*ones(1,nx*ny)-b2(img,:)+b(img,:),ny,nx));
% note how the average difference between the original and 
% reconstructed image decreases:
   disp('Average difference between the original and reconstructed')
   disp(['images is ' num2str(mean(abs(b2(img,:)-b(img,:))))])
   pause
end