function [dst,cwh,clm]=mcontour(pict,cfs,nlevs)
% MCONTOUR - 2D transformace Contourlet s omezenim poctu koeficientu pro rekonstrukci
%   MCONTOUR ziska pro vstupni obraz PICT vystup DST s omezenym poctem
%   nejvyznamnejsich koeficietu udavanou hodnotou CFS do urovne dekompozice 
%   dane poctem hodnot v NLEVS. Do NLEVS se uklada pocet smeru na jednotlivych 
%   urovnich dekompozice (2^NLEVS). Hodnoty NLEVS jsou cela kladna cisla. 
%   Cisla jsou oddelena carkou!
%   Priklad zapisu NLEVS: [0,1,2,4,4,5]
% 
%   Jako vstup se zadava primo nazev souboru.
% 
%   Obecny tvar funkce:
%   [DST,CWH,CLM]=MCONTOUR(PICT,CFS,NLEVS)
% 
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%   !!! POZOR!!! Transformace pracuje s obrazky stejneho rozmeru!!!
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%   
%   Krome omezeneho vystupu se ziskaji koeficienty vstupniho a vystupniho 
%   obrazu CWH a CLM. Koeficienty lze zobrazit pomoci SHOWPDFB a CONT2IM.
% 
%   Funkce dale urci hodnoty PSNR, SNR a MSE, ktere vypise do konzole.
%
%   Priklad zadavani funkce:
%   [dst,cwh,clm] = mcontour('lena.png',3000,[0,1,2,4,4,5]);
% 
%   Contourlet toolbox byl stazen z:
%   http://www.ifp.illinois.edu/~minhdo/software/contourlet_toolbox.zip
%   
%   Clanek:
%   Minh N. Do, Martin Vetterli 
%   The contourlet transform: an efficient directional multiresolution image representation
%   http://www.ifp.illinois.edu/~minhdo/publications/contourlet_txform.pdf

%% Nacteni obrazku a nastaveni filtru
x = imread(pict);
x = x(:,:,1); % Omezeni na jednu hladina v pripade barevneho obrazku
pfil = '5-3'; % Nazev filtru pyramidove dekompozice
dfil = 'pkva';% Nazev filtru smerove dekompozice
% pkva6', 'pkva8', 'pkva12', 'pkva'
%% Viceurovnova dekompozice
cwh = pdfbdec(double(x),pfil,dfil, nlevs);
%% Omezeni na CFS keoficientu
[c, s] = pdfb2vec(cwh); % Vystup dekompozice do radkoveho vektoru
len = length(c);
CC = zeros(1,len); % Definice vystupni promenne
[~,b]=sort(abs(c),'descend'); % Serazeni koeficientu sestupne
%======== Zapsani zvoleneho poctu koeficientu do CC ========%
i = 1;
while i < len+1 && i<cfs+1
    CC(b(i)) =c(b(i));
    i = i+1;
end
%===========================================================%
clm = vec2pdfb(CC, s); % Radkovy vektor do tvaru povudni dekompozice
%% Rekonstrukce z omezeneho poctu koeficientu
dst = pdfbrec(clm,pfil,dfil);
%% Vypocty SNR, PSNR, MSE 
xx = double(x)/256; %% xx - Vstupni hodnoty
xhat1 = (double(dst)/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
ssm = ssim(x,dst);
% snr = -10*log10( sum((xx(:) - xhat1(:)).^2) / sum(xx(:).^2) ); %SNR
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR
%% Vypis hodnot do konzole
disp('Contourlet transform');
disp(' ');
disp(['Transformace pro ',num2str(cfs),' koeficientu']);
disp(' ');
if  cfs > len
    disp('!!!Zadany pocet koeficientu je vetsi nez celkovy pocet koeficientu!!!')
    disp(['!!!Vysledek je zobrazen pro ',num2str(i-1),' koeficientu!!!']);
    disp(' ');
end
disp(['PSNR = ',num2str(PSNRdb),'dB']);
disp(['SNR  = ',num2str(snr),'dB']);
disp(['MSE  = ',num2str(error)]);
disp(['SSIM = ',num2str(ssm)]);
disp(' ');
disp(' ');
dst = uint8(dst);
% imm=showpdfb(cwh);