// Exemple de compression d'image par SVD
// M. Kern, d'apres un exemple de J. Demmel

// Lire l'image, calculer la SVD
m=320; n=200;
A=read('clown', m, n);
x=1:m; y=1:n;
[U, S, V]=svd(A);
xset("window",0); xbasc(); 
plot2d('nl', (1:n)', diag(S), style=-2);

gr=1;
// Image originale
xset("window",gr);xset('colormap',graycolormap(64)); gr=gr+1;
 xbasc(); grayplot(x,y,A);

chapeau=["SVD et compression d''image";
    'k : nombre de valeurs singulieres a retenir';
    ];
labels=["k"];
dims=list("vec",1);
defaults=["10"];

// Divers taux de compressions

ok=%t;
while ok do
  [ok, k] = getvalue(chapeau,labels,dims,defaults);
  if ok then
    if k>(n-1) then
      printf('Valeur de k trop grande %3d, k doit etre inferieur a %3d\n', k, n);
    else
      Ak=U(:, 1:k)*S(1:k,1:k)*V(:,1:k)';   
      // statistiques
      compress=S(1,1)/S(k+1,k+1); stock =(1-(m+n)*k/(m*n))*100;
      printf('k= %3d, Sig(1)/Sig(k): %6.3f, Economie de stockage: %6.3f %%', ...
	  k, compress, stock) 
      
      //Voir l'image compressee 
      xset("window",gr); xset('colormap',graycolormap(64)); gr=gr+1;
       xbasc(); grayplot(x,y,Ak);   
    end
  end
end