
%Fourierova transformace
function [output,coef,ocoef] = FOURIER(img,cfs)
%   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]  = FOURIER(IM,CFS)
%
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output,coef,ocoef] = FOURIER('lena.png',1024);
% 
%   funkce fftshift a fft jsou soucasti MATLABU

%% Nacteni  vstupniho obrazku
im = imread(img);
im = im(:,:,1); % omezeni na jednu hladina tedy na sedou v pripade barevneho obrazku

%% Výpočet transformacnich koeficientu pomocí Fourierové transformace

coef = fftshift(fft2(im));

%% Rozklad koeficientu z matice na vektor a nasledne omezeni podle zvoleneho poctu cfs koeficientu

m = size(coef);
vel = m(1)*m(2);
V = reshape(coef,1,vel); % vytvoreni vektoru z matice koeficientu a dosazeni hodnot coef
XX = zeros(1,vel); % vytvoreni prazdne matice
[~,S]=sort(abs(V),'descend'); % serazeni koeficientu dosazenych do vektoru po porade a to sestupne od nejvetsich po nejmensi
%% vlozeni vybraneho poctu koeficientu ze serazeneho vektoru do prazdne
%% vytvorene matice

i=1;
while i < vel + 1 && i < cfs + 1 
    XX(S(i)) =xe(S(i));
    i = i+1;
end

%%
ocoef = reshape(xe,s(1),s(2)); % matice z radkoveho vektoru xe

%% Rekonstrukce z omezeneho počtu koeficientů
output = ifft2(ifftshift(ocoef));

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

%% zobrazení rekonstruovaného obrazu
imshow(uint8(output))

%% Vypocty, PSNR, SNR, MSE, SSIM
xx = double(im)/256; 
xhat1 = (output/256); 
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;                   
err = abs(err(:));
ssm = ssim(im,output);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR

snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2));

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

output=uint8(output);