function [output,coef] = zasumWavelet(img,cfs)

%%Transformační metoda wavelet pro odsumovani
%   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.
%   
%   Funkce MakeONFilter slouzi k vyberu orthonormalniho filtru. Je mozna
%   volba 'Haar', 'Beylkin', 'Coiflet', 'Daubechies', 'Symmlet','Vaidyanathan','Battle'
%            
%   Funkce Ti_Denoise_Hard slouzi k odsumeni obrazu
%       Vstupni parametry:    Zasumeny obraz
%                               L - hloubka rozkladu
%                               qmf - quadraturni zrcadlovy filtr
%                               thresh - prahovani pouzivane na koeficienty
% 
% 


                            
%   Obecny tvar funkce:
%   [OUTPUT,COEF]  = zasumWavelet(IM,CFS)
%
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output,coef] = zasumWhearlet('lena.png',1);
%
%   output -  vystup pro metodu Shearlet
%% vstupní obraz
im = imread('lena.png');
im = double(im);
%% vytviření zašuměný
zasumeny = GWNoisy2(im,20);

%% imagesc(zasumeny),colormap(gray), axis image off

%% volání 1D zrcadlového kvadraturního filtru pro schearlet a wavelet
qmf2 = MakeONFilter('Daubechies',4);
qmf = MakeONFilter('Daubechies',4);

%% Použití wavelet transformace s ostrým prahem a využití rekonstrukce
im_wrec = TIDenoiseHard2(zasumeny,5,qmf,3*30);

%% Zobrazení pouze části obrazu
ind1 = 1:length(im);
ind2 = 1:length(im);

original =   im(ind1, ind2);
output =  ( im_wrec(ind1, ind2));

%% Vypocty PSNR, SNR, MSE a SSIM 

vst = double(zasumeny)/256; %%  vstupni hodnoty
vys = (double(output)/256); %%  ystupni hodnoty
chyb = mean(mean(((vst)-(vys)).^2));  % vypocet stredni kvadraticke chyby MSE
chyba = vst - vys;
chyba = chyba(:);
PSNR = 10 * log10(1/mean(chyba .^2));   % vypocet spockoveho pomeru signal-sum: PSNR
ssimm = ssim(zasumeny, output);               
SNR = 10 * log10(var(vst(:), 1) / mean(chyba(:).^2)); %vypocet pomer signal-sum: SNR

%% vypsani vystupnich parametru na displej 
disp(['Transformace pro ',num2str(cfs),' koeficientu wavelet metody']);

disp(' ');

disp(['PSNR = ',num2str(PSNR),'dB']);
disp(['SNR  = ',num2str(SNR),'dB']);
disp(['MSE  = ',num2str(chyb)]);
disp(['SSIM = ',num2str(abs(ssimm))]);


output=uint8(output);
