
close all;
clear all;

img = double(imread('Lena.jpg'));
n = size(img,1);
sigma = 20;        
is_real = 1;
cfs=8192;  
noisy_img = img + sigma*randn(n);

disp('Compute all thresholds');
X = randn(n);
tic; C = fdct_usfft(X,is_real); toc; 

% Compute norm of curvelets (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

% Take curvelet transform
disp(' ');
disp('Take curvelet transform: fdct_usfft');
tic; C = fdct_usfft(noisy_img,is_real); toc;


%Poz
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


maxValueLength = 10;


bigVector = [];
borders = [1];

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);
% x = reshape(d,1,vel); % Radkovy vektor z matice d
xe = zeros(1,vel); % Definice vystupni promenne
[~,Y]=sort(abs(bigVector),'descend'); % Serazeni koeficientu sestupne
%======== Zapsani zvoleneho poctu koeficientu do xe ========%
i=1;
while i < vel+1 && i<cfs+1
    xe(Y(i)) = bigVector(Y(i));
    i = i+1;
end

actualBorder = 1
% creation of new cell
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,:) = xe(borders(actualBorder):borders(actualBorder+1)-1);
            actualBorder = actualBorder + 1;
        end
        C2{1,i}{1,ii} = tmp2;
    end 
end
restored_img = ifdct_usfft(C2,is_real);


 imagesc(restored_img); colormap gray; ...
  

% disp(' ');
% disp('This method is of course a little naive.');
% R=input('Do you want to see the outcome of a more sophisticated experiment? [Y/N] ','s');
% disp(' ');
% 
% if strcmp(R,'') + strcmp(R,'y') + strcmp(R,'Y'), 
%   combined = double(imread('LenaCombined.jpg'));
%   figure; imagesc(combined); colormap gray; axis('image');
%   
% disp('This image (courtesy of Jean-Luc Starck) is the result of a more sophisticated');
% disp('strategy which involves both curvelets and wavelets. For a reference, please see');
% disp('J.L. Starck, D.L. Donoho and E. Candes, Very High Quality Image Restoration,')
% disp('in SPIE conference on Signal and Image Processing: Wavelet Applications');
% disp('in Signal and Image Processing IX, A. Laine, M. A. Unser and A. Aldroubi Eds,'); 
% disp('Vol 4478, 2001.');
% disp(' ');
% disp('For other experiments, please check the "enhanced denoise demo" in the wrapping');
% disp('directory.');
% disp(' ')
% end
