close all;    % Volitelne zavreni vsech oken Figure(X)
clear all;    % Volitelne mazani vsech promennych
%% Nacteni a zasumeni obrazu
im = imread('zoneplate.png');   % nacteni obrazku
im = im(:,:,1);                 % omezeni na 1 hladinu
obb = imnoise(im,'gaussian',0.005,0.005); % "zasumeny obraz"
%% Nastaveni vstupnich parametru
wavename{1} = 'db1';        % \
wavename{2} = 'db2';        %  \
wavename{3} = 'db3';        %   |
wavename{4} = 'db4';        %   |
wavename{5} = 'db6';        %   |
wavename{6} = 'db7';        %   | =>  Nejsou vypsany vsechny vlnky.
wavename{7} = 'db10';       %   | =>  Lze tuto skupinu libovolne doplnit.
wavename{8} = 'bior1.1';    %   |
wavename{9} = 'bior1.3';    %   |
wavename{10} = 'bior2.2';   %   | 
wavename{11} = 'bior3.3';   %  /
wavename{12} = 'bior4.4';   % /
level = 3;  % uroven dekompozive waveletove transformace
typ{1} = 'hard';        % \
typ{2} = 'soft';        %  |
typ{3} = 'firm';        %  | => vsechny volitelne typy prahovani
typ{4} = 'garrote';     %  |
typ{5} = 'hyperbolic';  % /
ssm1(1:5) = 1;  % Inicializacni hodnota MSE pro jednotlive prahovaci techniky
t1 = [0.00000000001,0.00000001,0.0000001, 0.000001,0.00001,0.0001, 0.001, 0.01, 0.1, 0.5, 0.8]; % Inicializace hodnot pro t1
mss = 1; % inicializace hodnoty pro vyber nejvhodnejsi vlnky
%% Odhad pro jednotlive typy vlnek
for q = 4:6
    [d,S] = wavedec2(double(obb),level,wavename{q});
    tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
    tt2(5) = 0; % |
%% Hruby odhad hodnoty prahu
for r = 1:5
    disp(typ);
for e = 1:11
	t2 = t1(e)*2;   % Pro "firm" je treba t2 nastavit na vetsi hodnotu nez t1.
                    % Timto se zaruci, ze t2 bude vzdy vetsi nez t1.
                    % Cislo nasobku je mozne zmenit podle potreb.
    Cr = d/max(max(d));                 % "normalizace" koeficientu
    y = threshold(Cr,t1(e),typ{r},t2);  % Prahovani koeficientu
    y(isnan(y)) = 0;                    % Osetreni "prazdnych" koeficientu
    CSS = y.*(max(max(d)));             % Zpetne nastaveni puvodnich hodnot koeficientu
    X = waverec2(CSS,S,wavename{q});    % Rekonstrukce prahovanych koeficientu

    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    error = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if error<ssm1(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = error;
        tt(r) = t1(e);
        disp('_');
        disp(num2str(tt));
        disp(num2str(ssm1));
    end
end
end
%% Jemny odhad hodnoty prahu
tteee = 0;
ssm1eee(1:5) = 1;
for r = 1:5
    disp(typ);
for ee = tt(r)/20:tt(r)/20:tt(r)*4
    t2 = ee*5; % Cislo nasobku pro prah t2 je mozne opet zmenit podle potreb.
    Cr = d/max(max(d));                 % "normalizace" koeficientu
    y = threshold(Cr,ee,typ{r},t2);     % Prahovani koeficientu
    y(isnan(y)) = 0;                    % Osetreni "prazdnych" koeficientu
    CSS = y.*(max(max(d)));             % Zpetne nastaveni puvodnich hodnot koeficientu
    X = waverec2(CSS,S,wavename{q});    % Rekonstrukce prahovanych koeficientu

nnz(CSS);
xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if error<ssm1eee(r)
        ssm1eee(r) = error;
        tteee(r) = ee;
        tt2(r) = t2;
        disp('_');
        disp(num2str(tteee));
        disp(num2str(ssm1eee));
    end
end
end
    if mss>error
        wname = wavename{q}; % ulozeni nazvu nejlepsi vlnky
        mss = error;
    end
end
[d,S] = wavedec2(double(obb),level,wname);
ssm1(1:5) = 1;
tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
tt2(5) = 0; % |
%% Hruby odhad prahu pro nejlepsi vlnku
for r = 1:5
    disp(typ);
for e = 1:11
	t2 = t1(e)*2;   % Pro "firm" je treba t2 nastavit na vetsi hodnotu nez t1.
                    % Timto se zaruci, ze t2 bude vzdy vetsi nez t1.
                    % Cislo nasobku je mozne zmenit podle potreb.
    Cr = d/max(max(d));                 % "normalizace" koeficientu
    y = threshold(Cr,t1(e),typ{r},t2);  % Prahovani koeficientu
    y(isnan(y)) = 0;                    % Osetreni "prazdnych" koeficientu
    CSS = y.*(max(max(d)));             % Zpetne nastaveni puvodnich hodnot koeficientu
    X = waverec2(CSS,S,wname);          % Rekonstrukce prahovanych koeficientu

    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    error = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if error<ssm1(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = error;
        tt(r) = t1(e);
        disp('_');
        disp(num2str(tt));
        disp(num2str(ssm1));
    end
end
end
%% Jemny odhad hodnoty prahu pro nejlepsi vlnku
tteee = 0;
ssm1eee(1:5) = 1;
for r = 1:5
    disp(typ);
for ee = tt(r)/20:tt(r)/20:tt(r)*4
    t2 = ee*5; % Cislo nasobku pro prah t2 je mozne opet zmenit podle potreb.
    Cr = d/max(max(d));                 % "normalizace" koeficientu
    y = threshold(Cr,ee,typ{r},t2);     % Prahovani koeficientu
    y(isnan(y)) = 0;                    % Osetreni "prazdnych" koeficientu
    CSS = y.*(max(max(d)));             % Zpetne nastaveni puvodnich hodnot koeficientu
    X = waverec2(CSS,S,wname);          % Rekonstrukce prahovanych koeficientu

nnz(CSS);
xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if error<ssm1eee(r)
        ssm1eee(r) = error;
        tteee(r) = ee;
        tt2(r) = t2;
        disp('_');
        disp(num2str(tteee));
        disp(num2str(ssm1eee));
    end
end
end
[aaa,bbb] = sort(ssm1eee);  % Serazeni hodnot MSE pro jednotlive prahovaci techniky
y = threshold(Cr,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1))); % Provedeni prahovani pouze pro nejvhodnejsi typ prahovani
y(isnan(y)) = 0;            % Osetreni "prazdnych" koeficientu
CSS = y.*(max(max(d)));     % Zpetne nastaveni puvodnich hodnot koeficientu        
X = waverec2(CSS,S,wname);  % Rekonstrukce prahovanych koeficientu
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
figure; imshow(obb);        % Zobrazeni zasumeneho obrazu

xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (double(X)/256); %% xhat1 - Vystupni hodnoty
error = mean(mean(((xx)-(xhat1)).^2));  %MSE
err = xx - xhat1;                   
err = abs(err(:));
PSNRdb = 10 * log10(1/mean(err .^2));   %PSNR
snr = 10 * log10(var(xx(:), 1) / mean(err(:).^2)); %SNR
ssm = ssim(im,X);   % SSIM
disp('________________________________');
disp(['t1   = ',num2str(tteee(bbb(1)))]);
disp(['type = ',typ{bbb(1)}]);
disp(['PSNR = ',num2str(PSNRdb),'dB']);
disp(['SNR  = ',num2str(snr),'dB']);
disp(['MSE  = ',num2str(error)]);
disp(['SSIM = ',num2str(ssm)]);
disp(['wname= ',wname]);
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
end
disp(' ');
disp(' ');