

%     program to illustrate the SVD in data compression .  The 
% original image is given by one of the inline functions . 
n = input('enter the meshsize n  ')
h = 2/n;
x = linspace(-2+h,2-h, n); y  = x;
[X,Y] = meshgrid(x,y);
disp(' The first three choices of images are built in. The fourth ')
disp(' choice allows you to enter your own image.  You must construct')
disp(' a matrix that you must name myZ. It should be nxn. Start with')
disp(' myZ = zeros(n,n) and fill in numbers between zero and 1 ')
disp('  ')
k = input('enter the choice of data 1,2, 3, or 4   ');

if k == 1
      f = inline('rand(size(x))', 'x', 'y')
      Z = f(X,Y);
 elseif k == 2
      f = inline('exp(-5*(sqrt(x.^2 +y.^2) - 1.5).^2)', 'x', 'y');
      Z = f(X,Y);
 elseif k == 3 
      f = inline(' 6*exp(-3*x.^2 - y.^2) + .5*x +y', 'x', 'y');
      z = f(X,Y);
      Z = 0*z;
      for j = 1:8
      Z = Z + (z< (j-3)+.25)  - (z < (j-3)-.25);
      end
 else 
      Z = myZ   % put in your own image matrix here.
 end

 % make singular value decomp. of the image matrix Z.
[U,S,V] = svd(Z);
 s = diag(S); % the vector of singular values.
clf
   subplot(1,2,1)
   image(64*Z); colormap(gray)

disp('Hit return to continue ')

% begin adding up the terms of the expansion.
Z1 = zeros(n, n);
 for j = 1:n
   pause
   u = U(:,j); v = V(:,j);
   Z1 = Z1 + u*s(j)*v';
   subplot(1,2,2)
   image(64*Z1); colormap(gray)
   disp(' ')
   disp('       k    singular value   ')
   format compact
   [j, s(j)]
end
   
 
