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

%%Transformační metoda schearlet
%   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 test_dst slouzi k vypoctu a omezeni kouficientu. Vstupnimi
%   hodnotami  je originalni obraz, velikost shearletu na kazde urovni,
%   pocet koeficientu a pouzity tyo filtru pro dekompozici
%
%   Obecny tvar funkce:
%   [OUTPUT,COEF]  = SHEARLET(IM,CFS)
%
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output,output1,ocoef] = SHEARLET('lena.png',1);
%
%   output -  vystup pro metodu Shearlet

%% Vstupní obraz:
im = imread(img);
im = double(im);
im = im(:,:,1); % omezeni na jednu hladina v pripade barevneho obrazku


%% volání 1D kvadraturního zrcadlového filtru pro shearlet
qmf2 = MakeONFilter('Daubechies',4);  %generováni orthonormalniho filtru pro transformaci 

%% Použití schearlet transformace s ostrým prahem a využití rekonstrukce ve
% dvou směrech (vertikální a horizontální)

[coef num1] = test_dst(im,[2 3 3 4 4],cfs,qmf2,qmf2);

%% výpočet poměru sígnál šum (psnr)
psnr1 = 20*log10(255/(1/512*norm(double(im(:))-coef(:))));

%% volání 9-7 kvadraturního filtru pro wavelet transformaci
[qmf,dqmf] = MakeBSFilter('CDF',[4 4]);

%% Použití wavelet transformace s ostrým prahem a využití rekonstrukce
[im_wrec num2 psnr2] = wave_comp(double(im),dqmf,qmf,4,cfs,1);

%% Zobrazení pouze části obrazu
 
 ind1 = 1:512;
 ind2 = 1:512;
 original =   im(ind1, ind2);
 
 output1 =  ( im_wrec(ind1, ind2));
 output2 =  ( coef(ind1, ind2));
 
 %% výpočty parametrů výstupního obrazu pro metofy wavelet
xx = double(im)/256;
xhat1 = (output1/256);
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
ssm = ssim(im, output1);
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR


disp(['Transformace pro ',num2str(cfs),' koeficientu metody wavelet']);
disp(' ');

disp(['PSNR = ',num2str(abs(PSNRdb)),'dB']);
disp(['SNR  = ',num2str(abs(snr)),'dB']);
disp(['MSE  = ',num2str(abs(error))]);
disp(['SSIM = ',num2str(abs(ssm))]);

 %% výpočty parametrů výstupního obrazu
xx = double(im)/256;
xhat1 = (output2/256);
error2 = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb2 = 10 * log10(1/mean(err .^2));   %PSNR
ssm2 = ssim(im, output2);
snr2 = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR


disp(['Transformace pro ',num2str(cfs),' koeficientu metody shearlet']);
disp(' ');

disp(['PSNR = ',num2str(abs(PSNRdb2)),'dB']);
disp(['SNR  = ',num2str(abs(snr2)),'dB']);
disp(['MSE  = ',num2str(abs(error2))]);
disp(['SSIM = ',num2str(abs(ssm2))]);
 





 
%  imagesc(original),colormap(gray), axis image off
  %imagesc(output1),colormap(gray), axis image off
%   imagesc(output2),colormap(gray), axis image off
% 


output1=uint8(output1);

output2=uint8(output2);

