% icomp.m 
% threshhold an image using global DFT


% parameter definition
% compression percentage = % retained coefficients 
cp = 5;


% get image from MATLAB image collection and show image
  load('clown');
% load('mandrill');

figure(1);
imagesc(real(X)); , colormap(gray), axis equal, axis tight 
title('Original')
 
% calculate and show fft 
Z=fft2(X);
maz =  max(abs(Z(:)));
figure(2);
imagesc(log(1+abs(fftshift(Z)))); , colormap(gray), axis equal, axis tight 
title('shifted log(1+|dft|)');
drawnow
disp('Press a key for thresholded signals'), pause

% threshold and reconstruct
sZ = sort(abs(Z(:)));
ci = round((1-cp/100)*length(sZ));
cutlev = sZ(ci);  
J = find(abs(Z) < cutlev);
cZ = Z;
cZ(J) = zeros(size(J));
cM = ones(size(Z));
cM(J) = zeros(size(J));

% reconstruct
cX = real(ifft2(cZ));
figure(3)
imagesc(cX);  colormap(gray), axis equal, axis tight 
title('Reconstruction');
drawnow;

figure(4)
imagesc(fftshift(cM)), colormap(gray), axis equal, axis tight 
title('Threshold Map');
drawnow;



% show compression distortion
sX = size(X,1)*size(X,2);
sJ = size(J,1)*size(J,2);
PercentRetainedCoefficients =100*(sX-sJ)/sX

DistortionPercent = 100*energy(cX-X)/energy(X)

% show cut off level (logarithmic)
  logcutlev =-log10(cutlev/maz)