%
%   This script reads a JPEG image, converts it to double
%   and computes the SVD.
%   
%   Then a succession of rank k reconstructions are
%   visualized for  k = 5:5:r
%
%       where  r is the smallest integer s.t.  s(r) < .01*s(1)
%
%       with s(j) the j-th singular value of the image
%
%--------------------------------------------------------
%
%   D.C. Sorensen
%   10 Sep. 04
%   Revised 21 Sep 09
%

clear

%
%  input the name of a file containing the JPEG image
%
%     e.g.  'Boffo.jpg'
%
   imstring = input('image file name = ')
   Z = imread(imstring);
%
% Gets representation as 64 bit arith. (double) words
%
   Z1 = im2double(Z);
%  Z1 = Z;
%
%  Z1 will contain an RGB representation of the image
%  This happend even when I put in a B and W image
%
%  So, Z1 is  m by n by 3   
%
%  Z1(m,n,:) = [R,G,B]    the Red, Green, and Blue
%                         values for the (m,n) pixel
%
%  For color you need an SVD for each
%  For B and W you only need one
%
   [U1,S1,V1] = svd(Z1(:,:,1));

%
%  Determine the cut off index r
%
%           s(r - 1) > .01*s(1))    but   s(r) .le. .01*s(1)
%
   s = diag(S1);
   r = 1;
   while (s(r) > .01*s(1)),
       r = r+1;
   end
%
%  Display original image
%
  iptsetpref('ImshowBorder', 'loose')
  mag = 2.7;
  [mz,nz] = size(Z1); 

   figure(1)
   imshow(Z1,64)
%  truesize(1)
%  truesize(1,[mz*mag,nz*mag])
disp('strike key')
pause
%
%  Display Singular Value Plot 
%

   figure(2)
   plot(s,'o')
disp('strike key')
pause
   semilogy(s,'o')
   hold
      semilogy(s(1:r),'r*')
   hold off
%
%  Display sequence of rank k approximations 
%

disp('strike key')
pause
   for k = 5:5:r, 
       G = U1(:,1:k)*S1(1:k,1:k)*V1(:,1:k)';
       imshow(G,64); 
%      truesize(1,[mz*mag,nz*mag])
       rankG = k, 
disp('strike key')
pause; 
   end