%metoda pro zobrazeni baze pro jeden zvoleny koeficient transformační metody Curvelet

close all;
clear all;

%načtení vstupního obrazu
img = zeros(256,256);

is_real = 1;

%stanovení počtu koeficientů
cfs=256;  
%vypocet transformacnich koeficientu
C = fdct_usfft(img,is_real);

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

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

C2{3}{2}(10,10) = 1;
%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;
  
 %%výpočty parametrů výstupního obrazu pro metofy wavelet
xx = double(img)/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(img, output);
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR


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

disp(['PSNR = ',num2str(abs(PSNRdb)),'dB']);
disp(['SNR  = ',num2str(abs(snr)),'dB']);
disp(['MSE  = ',num2str(abs(error))]);
disp(['SSIM = ',num2str(abs(ssm))]);
surf(output,'edgecolor','none');
figure;mesh(output);
o = fftshift(fft2(output));
mesh(abs(o));

