% tento soubor vyuziva funkce knihovny WaveAtom!
close all;    % Volitelne zavreni vsech oken Figure(X)
clear all;    % Volitelne mazani vsech promennych
%% Kontrola a nastaveni dwtmode
dt = dwtmode;
k = 'sym';
chng = 0;
if isequal(dt,k) ~= 1
    dwtmode('sym');
    chng = 1;
end
%% Nastaveni vstupnich promennych
PSNRdb2 = 1;
X = zeros(512,512);     % Inicializace
vstup = 'barbara128.png'; % 'lena64.png' , 'lena128.png', 
                          % 'barbara64.png' , 'barbara128.png' 
%%%%%%%%%%%%%%%%%%%% Nastaveni Wavelet %%%%%%%%%%%%%%%%%%%%
level = 2;  % 2 pro 128, 3 pro 64
wavename = 'db7';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%% Nastaveni WaveAtom %%%%%%%%%%%%%%%%%%%
trans = 'ortho';
pat = 'p';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Nastaveni Threshold %%%%%%%%%%%%%%%%%%%
% t1 =    0.004;  % nastaveni prahu
% sig =   0.0001; % 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 t22 = 0.02:0.02:0.02 % lze nastavit vice hodnot t2 jednim spustenim
    t2 = t22;            % nacteni aktualniho nastaveni t22 do t2 pri simulaci s vice t2
for t11 = 0.0002:0.0002:0.0002 % lze nastavit vice hodnot t1 jednim spustenim
    sig = t11/10;        % Volba kroku snizeni prahu/ lze nastavit napevno
    t1 = t11;            % nacteni aktualniho nastaveni t11 do t1 pri simulaci s vice t1
for i = 1:iter
    disp([num2str(i),'. iterace']);
    disp(['Threshold = ',num2str(t1)]);
    clear C;
    clear S;
%% wavedec2
    [C,S] = wavedec2(double(X),level,wavename);
%%  nacteni obrazku pro interpolaci
    imm = imread(vstup);
    %%%% kompenzace vlivu funkce wavedec2 na obrazek %%%%
    immm = wextend('2D','sym',imm,5);                   %
    immm = circshift(immm,[4 4]);                       %
    immm = immm(:,:,1);                                 %
    immm = double(immm);                                %
    immm = immm/max(max(immm));                         %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    siz = size(immm);
    CCut = immm(2:138,2:138); % Hodnoty platne jenom pro db7 !!
%% nahrani novych hodnot na patricne misto v promenne C
    CCCut = reshape(CCut,1,S(1)*S(2));
    if i == 1
        C = zeros(1,length(C)); % vymazani promenne
    else
        C = C/max(max(C));      % "normalizace" C
    end
    C(1:S(1)*S(2)) = CCCut;     % nahrani novych hodnot
%% waverec2
% C*(2^level)*256)  ==> hodnota koeficientu se po kazde urovni dekompozice
%                       zvysuje priblizne timto vztahem
%                       Tedy:   1. uroven - cca.  512
%                               2. uroven - cca. 1024
%                               2. uroven - cca. 2048
    temp = imread(vstup);
    temp = temp(:,:,1);
    temp = max(max(temp));
    temp = double(temp);
    kesz = waverec2(C.*((2^level)*temp),S,wavename);
%% Dekompozice pomoci WaveAtom
C = mefwa2sym(double(kesz),pat,trans);
%% Prahovani
    if strcmp(trans,'directional') == 1
        Cr1 = C(:,:,1)/max(max(C(:,:,1))); % "normalizace" hodnot
        Cr2 = C(:,:,2)/max(max(C(:,:,2)));
        y1 = threshold(Cr1,t1,typ,t2);
        y2 = threshold(Cr2,t1,typ,t2);
        CSS1 = y1.*(max(max(C(:,:,1))));
        CSS2 = y2.*(max(max(C(:,:,2))));
        CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
    elseif strcmp(trans,'ortho') == 1
        Cr1 = C./max(max(C));
        y1 = threshold(Cr1,t1,typ,t2);
        CSS = y1.*(max(max(C)));
    end

%% Rekonstrukce pomoci WaveAtom
    X = meiwa2sym(CSS,pat,trans);
    t1 = t1-sig; % snizeni prahu hodnotou sig
end
%% SNR, PSNR, MSE a SSIM interpolace opakovanym prahovanim
    img = imread('barbara.png');
    ff = double(img(2:512,2:512));
    xx = ff/256;
    figure;imshow(uint8(ff));
    xhat1 = (double(X(1:511,1:511))/256);
    figure;
    imshow(uint8(X(1:511,1:511)));
    error = mean(mean((abs(xx)-abs(xhat1)).^2));
    err = abs(xx) - abs(xhat1);
    err = err(:);
    PSNRdb1 = 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 WaveAtom koeficientu ');
    disp(' ');
    disp(['PSNR = ',num2str(PSNRdb1),'dB']);
    disp(['SNR = ',num2str(snr),'dB']);
    disp(['MSE = ',num2str(error)]);
    disp(['SSIM = ',num2str(abs(ssim(img(2:512,2:512),X(1:511,1:511))))]);
    disp(' ');
    disp(' ');
%% SNR, PSNR, MSE a SSIM bikubicke interpolace
    ff = double(img(2:512,2:512));
    xx = ff/256;
    clear xhat1;
    xhat1 = (double(imresize(imm, [512 512], 'bicubic'))/256);
    XX = xhat1*256;
    xhat1 = xhat1(2:512,2:512);
    figure;
    imshow(xhat1);
    figure;
    imshow(uint8(ff));
    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('Bikubicka interpolace');
    disp(' ');
    disp(['PSNR = ',num2str(PSNRdb),'dB']);
    disp(['SNR = ',num2str(snr),'dB']);
    disp(['MSE = ',num2str(error)]);
    disp(['SSIM = ',num2str(ssim(img,XX))]);
    disp(' ');
    disp(' ');
%% Uchova informaci o nejlepsich hodnotach v pripade spusteni vice simulaci
%  pomoci cyklu pro t11 a t22
    if PSNRdb2 < PSNRdb1 
            PSNRdb2 = PSNRdb1;
            tt1 = t11;
            tt2 = t2;
            sigg = sig;
%             disp(['PSNR = ',num2str(PSNRdb2),'dB']);
%             disp(['t1 = ',num2str(tt1)]);
%             disp(['t2 = ',num2str(tt2)]);
    end
    if rem(t1,0.02) ==0
            disp(['t1 = ',num2str(t1)]);
            disp(['t2 = ',num2str(t2)]);
    end
        disp(' ');
        disp('*****************************');
        disp(['aktualni t2 = ',num2str(t2),' t1 = ',num2str(t11)]);
        disp(['PSNR = ',num2str(abs(PSNRdb2)),'dB']);
        disp(['t1 = ',num2str(tt1)]);
        disp(['t2 = ',num2str(tt2)]);
        disp(['sig = ',num2str(sigg)]);
        disp('*****************************');
        disp(' ');
end
end
%% Nastavi povodni dwtmode
if chng == 1;
    dwtmode(dt);
end