function [dst,cwh,clm]  = mwavelet(pict,cfs,level,wavename)
% MWAVELET - 2D waveletova transformace s omezenim poctu koeficientu pro rekonstrukci
%   MWAVELET ziska pro vstupni obraz PICT vystup DST s omezenym poctem
%   nejvyznamnejsich koeficientu udavanou hodnotou CFS do urovne dekompozice 
%   LEVEL pomoci vlnky WAVENAME. Funkce vyuziva vlnky implementovane 
%   v programu MATLAB.
%   Seznam dostupny pomoci HELP WAVEINFO 
%   Pomoci CWH lze zobrazit dekompozici ze vsech koeficientu.
%   CLM obsahuje obraz dekompozice ze zvoleneho poctu keoficientu.
%   CWH a CLM lze zobrazit pomoci IMSHOW.
% 
%   Obecny tvar funkce:
%   [DST,CWH,CLM]  = MWAVELET(PICT,CFS,LEVEL,WAVENAME)
%   
%   Jako vstup se zadava primo nazev souboru.
% 
%   Funkce dale urci hodnoty PSNR, SNR a MSE, ktere vypise do konzole.
%   
%   Priklad zadavani funkce:
%   [dst,cwh,clm] = mwavelet('lena.png',4000,3,'db7');

%% Nacteni obrazku
f = imread(pict);
f = f(:,:,1); % omezeni na jednu hladina v pripade barevneho obrazku
%% Viceurovnova waveletova dekompozice
[C,S] = wavedec2(f,level,wavename);
%% Omezeni na CFS koeficientu
len = length(C);
CC = zeros(1,len); % Definice vystupni promenne
cwh = wave2im(C,S,level,wavename);
[~,y]=sort(abs(C),'descend');  % Serazeni koeficientu sestupne
%======== Zapsani zvoleneho poctu koeficientu do CC ========%
i = 1;
while i < len+1  && i<cfs+1
    CC(y(i)) =C(y(i));
    i = i+1;
end
%===========================================================%
%% Rekonstrukce z omezeneho poctu koeficientu
dst = waverec2(CC,S,wavename);
clm = wave2im(CC,S,level,wavename);
%% Vypocty SNR, PSNR, MSE 
xx = double(f)/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(f,dst);
% snr = -10*log10( sum((xx(:) - xhat1(:)).^2) / sum(xx(:).^2) ); %SNR
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2));
%% Vypis hodnot do konzole
disp('Discrete Wavelet 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);