function svd_demo_with_image ()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SVD_DEMO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 451-551
%% this demo shows an application of SVD (Singular Value Decomposition) 
%% to image processing,
%% references: [J. Daniel, B. Noble, "Applied Linear Algebra"]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% Jason Kyle, 2004/12 (based on svd_demo.m by M. Peszynska)
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% INSTRUCTIONS: 
%% calling sequence: 
%%        svd_demo_with_image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
[a,map] = imread('arcdbk.bmp');   %% Import Picture and Color Map
a = double(a);             %% Convert to Correct Data Type
n = size(a);
%spy(a);pause;                           %% Display Matrix

colormap(map);                          %% Apply Colormap
figure (1);
h = pcolor(a);                            %% Display Picture
set(h,'EdgeColor','interp');            %% Data Props
rotate(h,[0,0,1],180);
pause;
%%%%%%%%% find SVD of a
[u,s,v] = svd(a);
diag(s),                %% show the singular values of a
m=rank(a);              %% what is the rank of a ?

%%%%%%%%% find successive low rank approximations of a in matrix <appr>
appr = zeros(n(1),n(2));
i=1;
p=1;
figure(2);
while (i<m) 
    if (i+p > m) p = m-i; end;
    for j=i:(i+p)
      appr=appr+u(:,j)*s(j,j)*v(:,j)'; %% p Steps at a Time
      i = i+1;
    end;
    %% report how good is the current approximation
    fprintf('approximation of order %d error =%g\n',i,norm(a-appr,2)/norm(a));
    colormap(map);          %% Apply Correct ColorMap
    h=pcolor(appr);         %% Plot
    rotate(h,[0,0,1],180);  %% Matlab Imports Images Upside Down
    set(h,'EdgeColor','interp');
    pause;                  %% and display it

end;