function [dst,cwh,clm] = mwpak(pict,cfs,level,wname)
% MWPAK - 2D transformace Wavelet packet s omezenim poctu koeficientu rekonstrukce
%   MWPAK ziska pro vstupni obraz PICT vstup DST s omezenym poctem
%   nejvyznamnejsich koeficientu udavanou hodnotou CFS do urovne
%   dekompozice LEVEL pomoci vlnky WNAME. Funkce vyuziva vlnky
%   implementovane v programu MATLAB.
%   Seznam dostupny pomoci HELP WAVEINFO 
% 
%   MWPAK urci nejlepsi baze ohodnocenim podstromu pomoci Shannonovy 
%   entropie.
%   
%   Obecny tvar funkce:
%   [DST,CWH,CLM] = MWPAK(PICT,CFS,LEVEL,WNAME)
% 
%   Momentalne neplati:
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
%   !!! POZOR!!! Transformace pracuje s rozmery !!!
%   !!!       delitelne 2 do urovne LEVEL       !!!
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% 
%   V pripade zadani vstupu s jinym rozmerem bude vstup zmenen na potrebnou
%   velikost. Ve vysledku se doplnene nuly neobjevi, mohou svak ovlivnit
%   vypocet koeficientu.
%   ________________________________________________________________
% 
%   Transformace pracuje presne pro DWTMODE: 'per'
%   Pri pouziti 'per' se transformace provede pomoci DWT2/IDWT2.
%   Pro ostatni mody se dekompoice provede pomoci WAVEDEC2/WAVEREC2
%   WAVEDEC2/WAVEREC2 je pro delsi vlnky redundantni. V zavislosti 
%   na zvolenem modu je ovlivnen vyber nejvyznamnejsich koeficientu.
% 
%   Krome omezeneho vystupu se ziskaji koeficienty vstupniho a vystupniho 
%   obrazu CWH a CLM. Koeficienty lze zobrazit pomoci WPAK2IM. Kvarterni
%   strom lze zobrazit z koeficientu CWH pomoci WPAKTREE.
% 
%   Priklad zadavani funkce:
%   [dst,cwh,clm] = mwpak('lena.png',4000,3,'db2');

%% Nacteni obrazku
im = imread(pict);
im = im(:,:,1); % omezeni na jednu hladina v pripade barevneho obrazku
s = size(im);
%% Nastaveni potrebneho rozmeru obrazku
% if mod(s,2^(level+1)) > 0
%     obrazek = im2sc(im,level);
% else obrazek = im;
% end
obrazek = im;
%% DWTMODE= 'per'
if strcmp(dwtmode,'per') == 1
%% Viceurovnova dekompozice wavelet packet
    cwh = wpakdec(obrazek,level,wname);
%% Omezeni na CFS koeficientu  
    clm = wpaklim2(cwh,cfs);
%% Rekonstrukce z omezeneho poctu koeficientu
    dst = wpakrec(clm,wname);
    dst = dst(1:s(1),1:s(2)); % oriznuti vystupu na potrebnou velikost
else
%% Viceurovnova dekompozice wavelet packet pomoci wavedec2 pro DWTMODE~='per'
    [cwh,S] = wpakdec2(obrazek,level,wname);
%% Omezeni na CFS koeficientu
    clm = wpaklim(cwh,cfs);
%% Rekonstrukce z omezeneho poctu koeficientu
    dst = wpakrec2(clm,S,wname);
    dst = dst(1:s(1),1:s(2)); % oriznuti vystupu na potrebnou velikost
end
%% Vypocty SNR, PSNR, MSE 
xx = double(im)/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(im,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('Wavelet Packet transform');
disp(' ');
disp(['Transformace pro ',num2str(cfs),' koeficientu']);
disp(' ');
if cfs > wpaknnz(cwh)
    disp('!!!Zadany pocet koeficientu je vetsi nez celkovy pocet koeficientu!!!')
    disp(['!!!Vysledek je zobrazen pro ',num2str(wpaknnz(cwh)),' 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);