load clown
figure(1);image(X);axis image;
colormap('gray');
%colormap(map); % for color
title('Original image');
pause(2)

%% gray scale
[U,S,V] = svd(X,0);
for k=[1:9, 10:2:18, 21:3:30],
  Xk=U(:,1:k)*S(1:k,1:k)*V(:,1:k)'; 
  figure(2);image(Xk);axis image;
  colormap('gray')
  title(sprintf('k=%d',k));
  pause(1)
end

%% color : 
%T = [map(X,1) map(X,2) map(X,3)];
%X = reshape(T,[size(X) 3]);
%[m,n,p] = size(X);
%X = reshape(X,m,p*n);
%[V,S,U] = svd(X',0);
%k = 100; % or whatever 
%kv = 1:k; 
%Xend = U(:,kv)*S(kv,kv)*V(:,kv)';
%Xend = reshape(Xend,m,n,p);
%Xend(Xend<0)=0;
%Xend(Xend>1)=1;
%figure(2);image(Xend);axis image;
%title(sprintf('k=%d',k));
%colormap(map);

figure(3);plot(diag(S)./max(diag(S)),'o')
hold on;plot([1 length(diag(S))], [S(k,k) S(k,k)]./max(diag(S)));
energy=sum(diag(S(1:k,1:k)))/sum(diag(S));
title(sprintf('k=%d, ||A-A_k||_2/||A||_2=%g, energy=%g%%',...
	      k,S(k+1,k+1)/S(1,1),100*energy));

