
close all;
clear all;
%% načtení vstupního obrazu
img = double(imread('flinstones.png'));
n = size(img,1);
sigma = 20;        
is_real = 1;
%% stanovení počtu koeficientů
cfs=8192;
%% výpočty prahů
X = randn(n);
tic; C = fdct_usfft(X,is_real); toc; 

%% vypočet norem cuvrelet (Monte Carlo)
E = cell(size(C));
 for s=1:length(C)
   E{s} = cell(size(C{s}));
   for w=1:length(C{s})
     A = C{s}{w};
     E{s}{w} = median(abs(A(:) - median(A(:))))/.6745; % Estimate noise level with robust estimator
   end
 end

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


%% Použití prahování

for s = 2:length(C)
  thresh = 3*sigma + sigma*(s == length(C));
  for w = 1:length(C{s})
    Ct{s}{w} = C{s}{w}.* (abs(C{s}{w}) > thresh*E{s}{w});
  end
end

%% 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); % Definice vystupni promenne
[~,S]=sort(abs(bigVector),'descend'); % Serazeni koeficientu sestupne do vytvoreneho prazdneho vektoru
%% 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
%%  vytvoreni nove bunky 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);

%% zobrazení rekonstruovaného obrazu
imagesc(output); colormap gray; axis image off;
  
%% vypocty parametrů výstupního obrazu pro metodu Curevelt PSNR, SNR, MSE a
%% SSIM 

vst = double(noisy_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(noisy_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))]);
output=uint8(output);