function [dst,cwh,clm] = mfourier(pict,cfs)
% MFOURIER - 2D Fourierova transformace s omezenim poctu koeficientu pro rekonstrukci
%   MFOURIER ziska pro vstupni obraz PICT vystup DST s omezenym poctem
%   nejvyznamnejsich koeficientu udavanou hodnotou CFS. Krome omezeneho vystupu 
%   se ziska vstupni a vystupni obraz v kmitoctove oblasti FWH a FLM v logaritmickem meritku.
% 
%   Obecny tvar funkce:
%   [DST,CWH,CLM] = MFOURIER(PICT,CFS)
% 
%   Jako vstup se zadava nazev souboru.
% 
%   Funkce dale urci hodnoty PSNR, SNR a MSE, ktere vypise do konzole.
%   
%   Priklad zadavani funkce:
%   [dst,cwh,clm] = mfourier('lena.png',1200);

%% Nacteni obrazku
im = imread(pict);
im = im(:,:,1); % omezeni na jednu hladina v pripade barevneho obrazku
%% Fourierova transformace obrazku
d = fftshift(fft2(im));
%% Fourierova oblast ze vsech koeficientu
cwh = log(1+abs(d));
cwh = im2uint8(cwh/max(max(cwh)));
%% Rozklad a omezeni na CFS koeficientu
s = size(d);
vel = s(1)*s(2);
x = reshape(d,1,vel); % Radkovy vektor z matice d
xe = zeros(1,vel); % Definice vystupni promenne
[~,Y]=sort(abs(x),'descend'); % Serazeni koeficientu sestupne
%======== Zapsani zvoleneho poctu koeficientu do xe ========%
i=1;
while i < vel+1 && i<cfs+1 
    xe(Y(i)) =x(Y(i));
    i = i+1;
end
%===========================================================%
dd = reshape(xe,s(1),s(2)); % matice z radkoveho vektoru xe
%% Fourierova oblast ze zvoleneho poctu koeficientu
clm = log(1+abs(dd));
clm = im2uint8(clm/max(max(clm)));
%% Rekonstrukce z omezeneho poctu koeficientu
dst = ifft2(ifftshift(dd));
%% Vypocty, SNR, PSNR, MSE, SSIM
xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (dst/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;                   
err = abs(err(:));
ssm = ssim(im,dst);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
% ssm = ssim(f,dst);
% snr = -10*log10( sum((xx(:) - xhat1(:)).^2) / sum(xx(:).^2) ); %SNR
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2));
%% Vypis hodnot do konzole
disp('Discrete Fourier transform');
disp(' ');
disp(['Transformace pro ',num2str(cfs),' koeficientu']);
disp(' ');
if cfs > vel
    disp('!!!Zadany pocet koeficientu je vetsi nez celkovy pocet koeficientu!!!')
    disp(['!!!Vysledek je zobrazen pro ',num2str(i-1),' koeficientu!!!']);
    disp(' ');
end
disp(['PSNR = ',num2str(abs(PSNRdb)),'dB']);
disp(['SNR  = ',num2str(abs(snr)),'dB']);
disp(['MSE  = ',num2str(abs(error))]);
disp(['SSIM = ',num2str(abs(ssm))]);
disp(' ');
dst = real(dst);
dst = uint8(dst);