% tento soubor vyuziva funkce knihovny ContourletSD!
close all;    % Volitelne zavreni vsech oken Figure(X)
clear all;    % Volitelne mazani vsech promennych
%% Kontrola a nastaveni dwtmode
dt = dwtmode;
k = 'per';
chng = 0;
if isequal(dt,k) ~= 1
    dwtmode('per');
    chng = 1;
end
%% Nastaveni vstupnich promennych
X = zeros(512,512);     % Inicializace
%%%%%%%%%%%%%%%%%%%% Nastaveni Wavelet %%%%%%%%%%%%%%%%%%%%
level = 2;  % 2 pro 128, 3 pro 64
wavename = 'db7';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Nastaveni ContourletSD %%%%%%%%%%%%%%%%%
nlev_SD = [4  2 2 ];
smooth_func = @rcos;
Pyr_mode = 2;
dfilt = 'pkva';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Nastaveni Threshold %%%%%%%%%%%%%%%%%%%
t1 =    0.009;  % nastaveni prahu
sig =   0.0009; % hodnota zmensujici prah v kazde iteraci
t2 =    0.01;   % pro garrote, musi platit: t2>t1
                % V nektere literature se uvadi 3-krat vetsi threshold
                % pro hard nez pro soft. zatim jsem nechal 1. Lze menit 
                % poslednim parametrem prikazu z 1 na 3.
iter =  9;      % pocet opakovani
typ =   'garrote'; % 'hard','soft','firm','garrote','hyperbolic'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(['Interpolace pro pocet iteraci: ',num2str(iter)]);
disp(['Prahovani typu: ', typ]);

for i = 1:iter
    disp([num2str(i),'. iterace']);
    disp(['Threshold = ',num2str(t1)]);
    clear C; % kdyz jsem to nevynuloval, nekdy to delalo
    clear S; % divne vypocty   
%% wavedec2, v prvnim kroku ze vstupu X, v dalsich krocich z obrazku s jiz prahovanymi koeficienty
    [C,S] = wavedec2(double(X),level,wavename);
%%  nacteni obrazku pro interpolaci
    imm = wavedec2(imread('barbara.png'),level,wavename);
    immm = imm(1:128*128);
    immm = immm(:,:,1);
    immm = double(immm);
    siz = size(imm);
%% nahrani novych hodnot na patricne misto v promenne C
    if i == 1
        C = zeros(1,length(C)); % vymazani promenne
    end
    C(1:S(1)*S(2)) = double(immm); 
%% waverec2
    kesz = waverec2(C,S,wavename);
        if i == 1
            figure;
            imshow(uint8(kesz));
            title('Prvni zpetna vlnkova transformace')
        end    
%% Dekompozice pomoci ContourletSD
    y = ContourletSDDec(double(kesz), nlev_SD, Pyr_mode, smooth_func, dfilt);
    [CCont, SCont] = pdfb2vec(y);   % rozklad na koeficienty a rozmery matic
%% Prahovani
    Cr = CCont/max(max(CCont));     % "normalizace" koeficientu
    y = threshold(Cr,t1,typ,t2);    % Prahovani
    CSS = y.*(max(max(CCont)));     % Nastaveni keoficientu na jejich puvodni hodnoty
%% Rekonstrukce pomoci ContourletSD
    yy = vec2pdfb(CSS, SCont);
    X = ContourletSDRec(yy, Pyr_mode, smooth_func, dfilt);
    if i == 1;
        figure;
        X = uint8(X);
        imshow(X);
        title('Prvni iterace rekonstrukce');
    end
t1 = t1-sig; % snizeni prahu hodnotou sig
end
imm = reshape(immm,128,128);
imm = imm/((2^level)*256);
imm = circshift(imm,[-4 -4]);
figure;
imshow(uint8(X));  
title([num2str(i),'. iterace interpolace']);
%% SNR, PSNR, MSE a SSIM interpolace opakovanym prahovanim
xhat1 = (double(X(1:512,1:512))/256);
ime = imread('barbara.png');
xx =  (double(ime(1:512,1:512))/256);
error = mean(mean((abs(xx)-abs(xhat1)).^2));
err = abs(xx) - abs(xhat1);
err = err(:);
PSNRdb = 10 * log10(255/sqrt(mean(err .^2)));
snr = -10*log10( sum((abs(xx(:)) - abs(xhat1(:))).^2) / sum(xx(:).^2) );
disp(' ');
disp(' ');
disp('Interpolace pomoci prahovani ContourletSD koeficientu ');
disp(' ');
disp(['PSNR = ',num2str(PSNRdb),'dB']);
disp(['SNR = ',num2str(snr),'dB']);
disp(['MSE = ',num2str(error)]);
disp(['SSIM = ',num2str(ssim(ime,X))]);
disp(' ');
disp(' ');
%% Bikubicka interpolace
figure;
imc = imresize(ime,[128 128]);
imre = imresize(imc, [512 512], 'cubic');
imshow(imre);
title('Interpolace obrazku klasickou interpolaci');
%% Nastaveni puvodni dwtmode
if chng == 1;
    dwtmode(dt);
end
figure;
imshow(imread('barbara.png'));
title('Puvodni obrazek')