function [output,C,C2] = CURVELET(im,cfs)

%transformační metoda Curvelet

%   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,C,C2]  = CURVELET(IM,CFS)
%
%   kde C - vypoctene transformacni koeficienty
%       C2 - omezene nejvyznamejsi koeficienty podle hodnoty cfs
%
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output,C,C2] = CURVELET('lena.png',256);

% pro správnou funkcnost toolboxu je potreba vyuzit i jednotlive cesty k
% adresarim v hlavni slozce


%% načtení vstupního obrazu

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

%% Použití curvelet transformace pro získání výstupních koeficientů obrazu
is_real=1;
C = fdct_usfft(img,is_real);

%% funkce pro omezení počtu koeficientů


bigVector = [];         %vytvoreni prazdneho vektoru pro prevod koficientu
borders = [1];          %hranice pro vytvoření buňky C2


for i = 1:size(C,2);
    for ii = 1:size(C{1,i},2)
        tmp = C{1,i}{1,ii};
        for iii = 1:size(tmp,1)
           bigVector(end+1:(end+size(tmp,2))) = [tmp(iii,:)]; 
           borders(end+1) = borders(end) + size(tmp,2); 
        end
    end    
end

vel = length(bigVector);
XX = zeros(1,vel); % Vytvoreni prazdneho vektoru
[~,S]=sort(abs(bigVector),'descend'); % Serazeni koeficientu sestupne

%% Zapis zvoleneho poctu koeficientu do vystupni promenne XX
i = 1;
while i < vel + 1 && i < cfs+1
    XX(S(i)) = bigVector(S(i));
    i = i + 1;
end

actualBorder = 1;

%% vytvoření nových buněk pro zapsání omezených koeficientů

for i = 1:size(C,2)
    C2{1,i} = cell(size(C{1,i}));
    for ii = 1:size(C{1,i},2)
        tmp = C{1,i}{1,ii};
        tmp2 = [];
        for iii = 1:size(tmp,1)
            tmp2(iii,:) = XX(borders(actualBorder):borders(actualBorder+1)-1);
            actualBorder = actualBorder + 1;
        end
        C2{1,i}{1,ii} = tmp2;
    end 
end

%% zpětná rekonstrukce obrazu z omezeného počtu koeficientů

output = ifdct_usfft(C2,is_real);

%% vypocty parametrů výstupního obrazu pro metodu Curevelt PSNR, SNR, MSE a
%% SSIM 

vst = double(img)/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(img, output);               
SNR = 10 * log10(var(vst(:), 1) / mean(chyba(:).^2));

%% Zobrazení vypoctenych vystupnich parametru na displej

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

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

%% zobrazení rekonstruovaného obrazu

imagesc(output); colormap gray; axis image off;
  
output=uint8(output);

