%function [output,coef,ocoef] = ODSUMCONTOUR(im,cfs,NLEVS)

%%Transformační metoda ODSUMCONTOURLET

%   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 a vypoctene koeficienty dekompozici vstupniho
%   obrazu
%
%   Funkce v zaveru i urci hodnoty kvality obrazu jako jsou SNR,PSNR,MSE a
%   SSIM.
%
%   Obecny tvar funkce:
%   [OUTPUT2,COEF,OCOEF]  = CONTOUR(IM,CFS,NLEVS)
%
%    nlevs = pocet smeru na jednotlivých urovních dekompozice
 
% 
%   Vstupem se zadava primo nazev souboru.
%
%   Priklad zadavani funkce:
%   [output,coef,ocoef] = ODSUMCONTOUR('lena.png',4000,[0, 0, 4, 4, 5]);
% 
%   Vystupy:    output  -rekonstuovaný obraz z omezeného počtu koeficientů
%               coef    -transformacni koeficienty
%               ocoef   -omezene transformacni koeficienty
%
%   Contourlet toolbox byl stazen z:
%   http://www.ifp.illinois.edu/~minhdo/so

%% Parametry
cfs=1;                   % zvoleny pocet koeficientu pro omezeni
pfilt = '5-3';
dfilt = 'pkva';
%% nacteni vstupniho obrazu

im = imread('lena.png');
im = double(im);

%% Vytvoreni zasumeneho obrazu

n = size(im,1);
sigma = 20;        
is_real = 1;
nim = im + sigma*randn(n);

%% Contourlet transform
coef = pdfbdec(nim, pfilt, dfilt, nlevs);

%% Omezeni transformovanych koeficientu na cfs keoficientu

[x, y] = pdfb2vec(coef); % rozlozeni matice vznikle z dekompozice do radkoveho vektoru
len = length(x);
XX = zeros(1,len); % definica promene pro zapsani omezenych koeficientu
[~,a]=sort(abs(x),'descend'); % Serazeni sestupne v radkovem vektoru

%% Vlozeni omezenych koeficientu do vytvorene promenne XX

i = 1;
while i < len+1 && i < cfs+1
    XX(a(i)) =x(a(i));
    i = i+1;
end

%% Vytvoreni matice s omezenymi nejvyznamnejsimi koeficienty

ocoef = vec2pdfb(XX, y); % Prevod radkoveho vektoru do tvaru povudni matice dekompozice

%% Rekonstrukce obrazu z omezenych nejvyznamnejsich koeficientu

output = pdfbrec(ocoef, pfilt, dfilt);

%% Vypocty PSNR, SNR, MSE a SSIM 

vst = double(nim)/256; %%  vstupni hodnoty
vys = (double(output)/256); %%  ystupni hodnoty
chyb = mean(mean(((vst)-(vys)).^2));  % vypocet stredni kvadraticke chyby MSE
chyba = vst - vys;
chyba = chyba(:);
PSNR = 10 * log10(1/mean(chyba .^2));   % vypocet spockoveho pomeru signal-sum: PSNR
ssimm = ssim(nim, output);               
SNR = 10 * log10(var(vst(:), 1) / mean(chyba(:).^2));

%% Zobrazení vypoctenych vystupnich parametru na displej
disp(['Transformace pro ',num2str(cfs),' koeficientu Contourlet metody']);

disp(' ');

disp(['PSNR = ',num2str(PSNR),'dB']);
disp(['SNR  = ',num2str(SNR),'dB']);
disp(['MSE  = ',num2str(chyb)]);
disp(['SSIM = ',num2str(abs(ssimm))]);

%%vykresleni vysptuniho odsumeneho obrazu

imagesc(output,[0 255]);axis image off;colormap gray(256);

output=uint8(output);