function [dst,cwh,clm] = mfrat(pict,cfs,cut)
% MFRAT - 2D transformace Radon s omezenim poctu koeficientu pro rekonstrucki
%   MFRAT ziska pro vstupni obraz PICT vystup DST s omezenym poctem
%   nejvyznamnejsich koeficientu udavanou hodnotou CFS. 
%   MFRAT provadi dekompozici na jedne urovni. 
% 
%   Jako vstup se zadava primo nazev souboru.
% 
%   Obecny tvar funkce:
%   [DST,CWH,CLM] = MFRAT(PICT,CFS,CUT)
% 
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
%   !!! POZOR!!! Transformace pracuje s prvocisly !!!
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%   
%   V pripade zadani vtupu s rozmerem jinym nez prvocislo, bude vstup
%   zmenen na potrebnou velikost. V promenne CUT lze zvolit 'y' pro omezeni
%   rozmeru na potrebnou velikost, nebo 'n' pro doplneni nulami na
%   potrebnou velikost. Ve vysledku se doplnene nuly neobjevi, mohou vsak
%   ovlivnit vypocet koeficientu.
%   Jestlize jsou rozmery vstupu stejna prvocisla, CUT se nemusi zadavat.
%
%   Krome omezeneho vystupu se ziskaji koeficienty vstupniho a vystupniho 
%   obrazu CWH a CLM. Koeficienty lze zobrazit pomoci FRXT2IM.
% 
%   Funkce dale urci hodnoty PSNR, SNR a MSE, ktere vypise do konzole.
% 
%   Priklad zadavani funkce:
%   [dst,cwh,clm] = mfrat('lena.png',3400,'y');
% 
%   FRIT Toolbox obsahujici FRAT byl stazen z:
%   http://www.ifp.illinois.edu/~minhdo/software/frit_toolbox.tar.gz
%   
%   Clanek:
%   Minh N. Do, Martin Vetterli
%   The finite ridgelet transform for image representation
%   http://www.ifp.illinois.edu/~minhdo/publications/FRIT_trans.pdf

%% Nacteni obrazku
im = imread(pict);
im = im(:,:,1); % omezeni na jednu hladina v pripade barevneho obrazku
%% Omezeni nebo doplneni rozmeru na potrebnou velikost
u = size(im);
g = 0;
%============= Nastaveni na potrebny rozmer ================%
if isprime(u(1))==0 || isprime(u(2)) == 0
    out = im2prim (im, cut);
    im = out;
    if cut == 'n'
        g = 2;
    end
else g=1;
end
im = double(im); % FRAT pracuje s double
%% Dekompozice Ridgelet
[cwh,m] = frat(im);
%% Omezeni na CFS koeficientu
s = size(cwh);
x = reshape(cwh,1,s(1)*s(2)); % Rozklad na radkovy vektor
len = length(x);
cffs = cfs; % Ulozeni hodnoty kvuli pozdejsimu pouziti
if cffs > len       % Pri volbe CFS vetsi nez pocet keoficientu,
    cffs = len;     % se pocet omezi na celkovy pocet koeficientu
end
xx = sort(abs(x),'descend'); % Serazeni koeficientu sestupne
val = xx(cffs); % cffs-ta nejvetsi hodnota (nejmensi hodnota pro omezeni)
clm = cwh .* (abs(cwh)>=val); % Prahovani hodnot

%% Rekonstrukce z omezeneho poctu koeficientu
imrec = ifrat(clm,m);
%=========== Oriznuti vystupu po doplneni nulami ===========%
if g == 2;
    im = imread(pict);
    im = im(:,:,1);
    imsiz = size(im);
    imrec = imrec(1:imsiz(1),1:imsiz(2));
end
%===========================================================%
xx = double(im)/256;
xhat1 = (imrec/256);
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
ssm = ssim(im, imrec);
% snr = -10*log10( sum(abs((xx(:)) - abs(xhat1(:))).^2) / sum(xx(:).^2) );
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR
%% Vypis hodnot do konzole
disp('Finite Radon transform');
disp(' ');
disp(['Transformace pro ',num2str(cfs),' koeficientu']);
disp(' ');
if  cfs > len
    disp('Zadany pocet koeficientu je vetsi nez celkovy pocet koeficientu!');
    disp(['Vysledek je zobrazen pro ',num2str(len),' koeficientu.']);
    disp(' ');
end
disp(['PSNR = ',num2str(PSNRdb),'dB']);
disp(['SNR  = ',num2str(snr),'dB']);
disp(['MSE  = ',num2str(error)]);
disp(['SSIM = ',num2str(ssm)]);
disp(' ');
disp(' ');
dst = uint8(imrec);