%function [output,coef,ocoef] = CONTOUR(img,cfs)

%%Transformační metoda CONTOURLET

%   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] = CONTOUR('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/software/contourlet_toolbox.zip


cfs=4096;

% Parametry
pfilt = '9-7';
dfilt = 'pkva';
nlevs = [0, 0, 4, 4, 5];    % pocet smeru na jednotlivých urovních dekompozice

%% nacteni vstupniho obrazu
im = zeros(256.256); %Vytvoreni 

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

%% Omezeni na CFS keoficientu
[c, s] = pdfb2vec(coef); % Vystup dekompozice do radkoveho vektoru
len = length(c);
CC = zeros(1,len); % Definice vystupni promenne
[~,b]=sort(abs(c),'descend'); % Serazeni koeficientu sestupne
%% ======== Zapsani zvoleneho poctu koeficientu do CC ========%
i = 1;
while i < len+1 && i<cfs+1
    CC(b(i)) =c(b(i));
    i = i+1;
end
%% ===========================================================%
ocoef = vec2pdfb(CC, s); % Radkovy vektor do tvaru povudni dekompozice


%% Rekonstrukce
ocoef = vec2pdfb(CC, s);
ocoef{3}{2}(12,10) = 1;


output = pdfbrec(ocoef, pfilt, dfilt);

%% Vypocty SNR, PSNR, MSE 
xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (double(output)/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;
err = err(:);
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
ssm2 = ssim(im, output);
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2));


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


%imagesc(output,[0 255]);axis image off;colormap gray(256);
surf(output,'edgecolor','none');
figure;mesh(output);
o = fftshift(fft2(output));
mesh(abs(o));



