% SVD demo - Bringing Lena in focus  
% This is ported to Octave -
%   corresponding Matlab commands are indicated as comments
% Requires mat2gray.m from the "image" package from the octave-forge repository

close all;
clear all;

A=imread("lenna.png"); % [A,map]=imread("lena.gif");
B=A;                  % B=im2double(A,'indexed');

imshow(B), axis("image") % imshow(B,map)

mn=size(B); m1=mn(1); m2=mn(2); depth=mn(3);
origsize=m1*m2*depth
disp([' Image dimensions in pixels: ',num2str(mn)])
disp([' Bytes in image = ',num2str(m1*m2*depth)])


% Compute the SVD
for i=1:depth
  %[u(:,:,i),s(:,:,i),v(:,:,i)]=svd(fft( B(:,:,i)));
  [u(:,:,i),s(:,:,i),v(:,:,i)]=svd( ( B(:,:,i)));
endfor

figure(2)
for pct=[5, 10, 15, 20, 25, 30];
  k= floor( pct * min( m1, m2 ) / 100 );

  for i=1:depth
    %C(:,:,i)=real( ifft( u(:,1:k,i)*s(1:k,1:k,i)*(v(:,1:k,i)') ) ) ;
    C(:,:,i)=u(:,1:k,i)*s(1:k,1:k,i)*(v(:,1:k,i)');
  endfor
  C = uint8( floor( C  ) );
  %C=mat2gray(C); %C=floor(C);

  imshow(C), axis("image")  % imshow(C,map)
  imwrite( C, strcat( 'leena-svd-', num2str(pct), '-pct.png') );
  newsize=depth * (m1*k+m2*k+k)
  disp(['Rank k=', num2str(k),'    Bytes needed = ', num2str(newsize), ...
    '	Compression = ', num2str(newsize/origsize * 100), '%' ])
  pause
end
