

function [output,coef,ocoef] = SURELET(img,cfs)
%transformační metoda SURELET

%   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, vypoctene koeficienty ze vstupniho obrazu a
%   omezene koeficienty.
%
%   Funkce v zaveru i urci hodnoty kvality obrazu jako jsou SNR,PSNR,MSE a
%   SSIM.
%
%   Obecny tvar funkce:
%   [OUTPUT,COEF,OCOEF]  = SURELET(IM,CFS)
%
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output,coef,ocoef] = SURELET('lena.png',4000);


%%
% Načtení vstupního obrazu

im = double(imread(img));
im = im(:,:,1); % omezeni na jednu hladinu v pripade barevneho obrazku
%% vstupní parametry
sigma = 20;             % hodnota standardu zasumeni
wtype = 'sym8';         % typ orthonormální filtru
PLOT  = 'on';           % volba vykresleni
parent=0;
%%
%%upravení velikosti vstupu
if(nargin<2)
    wtype = 'sym8';
end
[nx,ny] = size(im);
 
%%
%Výpočet počtu iterací J

J = aux_num_of_iters([nx,ny]);
if(max(J)==0)
    disp('The size of the signal is too small to perform a reliable denoising based on statistics.');
    output = im;
    return;
end 
%% Viceurovnova dekompozice
%výpočet koeficientu 
coef = fft_wavedec(im,wtype,J);
%% převod matice koeficientů ná řádkový vektor a jejich následné seřazení

s=size(coef);
vel = m(1)*m(2);

X = reshape(coef,1,vel); % Radkovy vektor z matice d
XX = zeros(1,vel); % Definice vystupni promenne
[~,S]=sort(abs(X),'descend'); % Serazeni koeficientu sestupne
%%
%======== Zapsani zvoleneho poctu koeficientu do xe ========%
i=1;
while i < vel + 1 && i < cfs + 1 
    (S(i)) =X(S(i));
    i = i+1;
end
%%
%omezená matice koeficientů dd
dd = reshape(xe,s(1),s(2));
ocoef=dd;
%%
%Zpetná rekonstrukce output pomocí funkce fft_waverec 
output = fft_waverec(dd,wtype,J);


%%
%zobrazení rekonstruovaného obrazu
imagesc(output,[0 255]);axis image off;colormap gray(256);


%%
%%výpočty parametrů výstupního obrazu
xx = double(im)/256;
xhat1 = (output/256);
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
ssm = ssim(im, output);
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR

%% zobrazení parametru na displej
disp(['Transformace pro ',num2str(cfs),' koeficientu']);
disp(' ');

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

output=uint8(output);

