function [dst,cwh,clm] = mwatom(pict,cfs,ps,xform,cut)
% MWATOM - 2D transformace WaveAtom s omezenim poctu koeficientu pro rekonstrukci
%   MWATOM ziska pro vstupni obraz PICT vystup DST s omezenym poctem
%   nejvyznamnejsich koeficientu udavanou hodnotou CFS. Pomoci PS lze
%   zvolit typ kmitoctove deleni. Je mozne vybrat ze dvou typu:
%       - 'p'   - Navrhnul Lexing Ying
%       - 'q'   - Navrhnul Laurent Demanet
%   Pro podrobnosti: freq_pat.m
%   Promenna XFORM urcuje typ transformace:
%       - 'ortho'       > ortonormalni bazove funkce
%       - 'directional' > smerove bazove funkce
%       - 'complex'     > komplexni bazove funkce
%   
%   Obecny tvar funkce:
%   [DST,CWH,CLM] = MWATOM(PICT,CFS,PS,XFORM,CUT)
% 
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%   !!! POZOR!!! Transformace pracuje s rozmery 2^n !!!
%   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
% 
%   V pripade zadani vstupu s rozmerem jinym nez 2^n, 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.
%   Pri zadani vstupu s rozmerem 2^n neni potreba promennou CUT definovat  
% 
%   Krome omezeneho vystupu se ziskaji koeficienty vstupniho a vystupniho 
%   obrazu CWH a CLM. Koeficienty lze zobrazit pomoci WATOM2IM.
% 
%   Funkce dale urci hodnoty PSNR, SNR a MSE, ktere vypise do konzole.
% 
%   Priklad zadavani funkce:
%   [dst,cwh,clm] = mwatom('lena.png',4901,'p','directional','n');
%   
%   WaveAtom toolbox byl stazen z:
%   http://www.waveatom.org/software/WaveAtom-1.1.1.tar.gz
% 
%   Clanek:
%   Laurent Demanet, Lexing Ying
%   Wave Atoms and Sparsity of Oscillatory Patterns
%   http://www.waveatom.org/papers/WaveatomsImage.pdf

%% Nacteni obrazku
im = imread(pict);
im = im(:,:,1); % omezeni na jednu hladinu v pripade barevneho obrazku
%% Omezeni nebo doplneni rozmeru na potrebnou velikost
u = size(im);
%================ Porovnani rozmeru s 2^n ==================%
for i = 1:2
    dim = 2;
    while dim < u(i)
        dim = dim*2;
        if u(i) ~= dim
            res(i) = 1;
        else res(i) = 2;
        end
    end
end
%========= Omezeni/Doplneni do potrebneho rozmeru ==========%
if res(1) == 1 || res(2) == 1
pict = im2pwsym(im,cut);
im = pict;
end
%===========================================================%
%% Viceurovnova dekompozice Ridgelet
cwh = fwa2sym(im,ps,xform);
%% Omezeni na CFS koeficientu
x = cwh(:); 
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 - VAL urcuje prah
%% Rekonstrukce z omezeneho poctu koeficientu
dst = iwa2sym(clm,ps,xform);
%=========== Oriznuti vystupu po doplneni nulami ===========%
if (res(1) == 1 && cut == 'n') || (res(2) == 1 && cut == 'n')
    dst = dst(1:u(1),1:u(2));
    im = im(1:u(1),1:u(2));
end
%===========================================================%
%% Vypocty SNR, PSNR, MSE 
xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (dst/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
ssm = ssim(im,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('WaveAtom 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(abs(PSNRdb)),'dB']);
disp(['SNR  = ',num2str(abs(snr)),'dB']);
disp(['MSE  = ',num2str(abs(error))]);
disp(['SSIM = ',num2str(abs(ssm))]);
disp(' ');
disp(' ');
dst = uint8(dst);