% function [output,coef,ocoef] = WAVELET(im,cfs,level,wavename)

%%Transformační metoda schearlet
%% přidání transformační metody wavelet pro porovnání

%   ziska pro vstupni obraz im vystup output, ktery je s omezenym poctem
%   nejvyznamnejsich koeficientu. Ty jsou zadany pod hodnotou cfs.
%   Z možných vystupnich hodnot je na vyber výsledná rekonstrukce obrazu z
%   daneho poctu koeficientu a vypoctene koeficienty dekompozici vstupniho
%   obrazu
%
%   Funkce v zaveru i urci hodnoty kvality obrazu jako jsou SNR,PSNR,MSE a
%   SSIM.
%
%   Obecny tvar funkce:
%   [OUTPUT2,COEF,OCOEF]  = WAVELET(IM,CFS, LEVEL, WAVENAME)
%
%   level - uroven dekompozice  
%   Wavename - typ vlnky
% 
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output2,coef,ocoef] = WAVELET('lena.png',4000,3,'db7');


cfs=8192;
level=3;
wavename='db7';



%% Nacteni obrazku
f = zeros(256.256);

%% Viceurovnova waveletova dekompozice
[R,S] = wavedec2(f,level,wavename);

%% Omezeni na CFS koeficientu
len = length(R);
CC = zeros(1,len); % Definice vystupni promenne
coef = wave2im(R,S,level,wavename);
[~,y]=sort(abs(R),'descend');  % Serazeni koeficientu sestupne
%% ======== Zapsani zvoleneho poctu koeficientu do CC ========%
i = 1;
while i < len+1  && i<cfs+1
    CC(y(i)) =R(y(i));
    i = i+1;
end



%% Rekonstrukce z omezeneho poctu koeficientu
output = waverec2(CC,S,wavename);
ocoef = wave2im(CC,S,level,wavename);



%% Vypocty SNR, PSNR, MSE 
xx = double(f)/256; %% xx - Vstupni hodnoty
xhat1 = (double(output)/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,output);
% 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 pro metodu Wavelet']);
disp(' ');
disp(['PSNR = ',num2str(PSNRdb),'dB']);
disp(['SNR  = ',num2str(snr),'dB']);
disp(['MSE  = ',num2str(error)]);
disp(['SSIM = ',num2str(ssm)]);
disp(' ');
disp(' ');
%%
output = uint8(output);