%============================================================================
% Beeldverwerking cursus - Prof. W. Philips
% Practica
% Demo voor DCT kompressie
% Aleksandra Pizurica TELIN/RUG 02/2001
%============================================================================


[filename,directory]=uigetfile ('*.tif');
I=imread(filename,'tif');
I=double(I);

M=input('block size (4,8,16,32): ');
N=M;
percent=input('Percent of the coefficients to remain in each block (%): ');
v=round(M*N*(1-percent/100));

dctm = dctmtx(M);

% Old call:  blocksDCT = blkproc(I,[M N],'dct2');
blocksDCT = blkproc(I,[M N],'P1*x*P2',dctm,dctm.');

DCTvar = im2col(blocksDCT,[M N],'distinct').';
n = size(DCTvar,1);
%DCTvar = (sum(DCTvar.*DCTvar) - (sum(DCTvar)/n).^2)/n;
DCTvar = (sum(DCTvar.*DCTvar)/n - (sum(DCTvar)/n).^2);

[dum,order] = sort(DCTvar);

mask=ones(M,N);
mask(order(1:v))=0;
rem=length(find(mask));

%Old call: Icompr = blkproc(blocksDCT,[M N],'idct2(x.*P1)',mask(1:M,1:N));
Icompr = blkproc(blocksDCT,[M N],'P1*(x.*P2)*P3',dctm.',mask(1:M,1:N),dctm);
Idiff=I-Icompr;
PSNR=20*log10(max(max(I))/std2(Icompr-I));

figure
subplot(2,2,1),imshow(I,[0 255]);
title('Original image');
subplot(2,2,2),imshow(Icompr,[0 255]);
Title_rec=strcat('Reconstructed image',' PSNR=',num2str(PSNR),'dB');
title(Title_rec);
colormap(gray);
subplot(2,2,3),imshow(Idiff,[]);
title('Difference image');
subplot(2,2,4),imagesc(mask);
Title_mask=strcat('Selected coefficients:.', num2str(rem), ' .out of .', num2str(M*N));
title(Title_mask);
colormap(gray);