% tento soubor vyuziva funkce knihovny WaveAtom!
close all;    % Volitelne zavreni vsech oken Figure(X)
clear all;    % Volitelne mazani vsech promennych
%% Nacteni a zasumeni obrazu
im = imread('text.png');   % nacteni obrazku
im = im(:,:,1);                 % omezeni na 1 hladinu
obb = imnoise(im,'gaussian',0.005,0.005); % "zasumeny obraz"
%% Nastaveni vstupnich parametru
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
trans = 'ortho';
pat = 'q';
tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
tt2(5) = 0; % |
d = mefwa2sym(double(obb),pat,trans); % Dekompozice WaveAtom
%% 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.
    if strcmp(trans,'directional') == 1 % Prahovani pro trans = directional
        Cr1 = d(:,:,1)/max(max(d(:,:,1))); 
        Cr2 = d(:,:,2)/max(max(d(:,:,2)));
        y1 = threshold(Cr1,t1(e),typ{r},t2);
        y2 = threshold(Cr2,t1(e),typ{r},t2);
        y1(isnan(y1)) = 0;
        y2(isnan(y2)) = 0;
        CSS1 = y1.*(max(max(d(:,:,1))));
        CSS2 = y2.*(max(max(d(:,:,2))));
        CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
    elseif strcmp(trans,'ortho') == 1 % Prahovani pro trans = ortho
        Cr1 = d./max(max(d));
        y1 = threshold(Cr1,t1(e),typ{r},t2);
        CSS = y1.*(max(max(d)));
    elseif strcmp(trans,'complex') == 1 % Prahovani pro trans = complex
        Cr1 = d(:,:,1)/max(max(d(:,:,1)));
        Cr2 = d(:,:,2)/max(max(d(:,:,2)));
        Cr3 = d(:,:,3)/max(max(d(:,:,3)));
        Cr4 = d(:,:,4)/max(max(d(:,:,4)));
        y1 = threshold(Cr1,t1(e),typ{r},t2);
        y2 = threshold(Cr2,t1(e),typ{r},t2);
        y3 = threshold(Cr3,t1(e),typ{r},t2);
        y4 = threshold(Cr4,t1(e),typ{r},t2);
        y1(isnan(y1)) = 0;
        y2(isnan(y2)) = 0;
        y3(isnan(y2)) = 0;
        y4(isnan(y2)) = 0;
        CSS1 = y1.*(max(max(d(:,:,1))));
        CSS2 = y2.*(max(max(d(:,:,2))));
        CSS3 = y3.*(max(max(d(:,:,3))));
        CSS4 = y4.*(max(max(d(:,:,4))));
        CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
        CSS(:,:,3) = CSS3; CSS(:,:,4) = CSS4;
    end
    X = meiwa2sym(CSS,pat,trans); % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    error = mean(mean((abs(xx)-abs(xhat1)).^2));  %MSE
    if error<ssm1(r)
        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.
    if strcmp(trans,'directional') == 1 % Prahovani pro trans = directiona
        Cr1 = d(:,:,1)/max(max(d(:,:,1)));
        Cr2 = d(:,:,2)/max(max(d(:,:,2)));
        y1 = threshold(Cr1,ee,typ{r},t2);
        y2 = threshold(Cr2,ee,typ{r},t2);
        y1(isnan(y1)) = 0;
        y2(isnan(y2)) = 0;
        CSS1 = y1.*(max(max(d(:,:,1))));
        CSS2 = y2.*(max(max(d(:,:,2))));
        CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
    elseif strcmp(trans,'ortho') == 1 % Prahovani pro trans = ortho
        Cr1 = d./max(max(d));
        y1 = threshold(Cr1,ee,typ{r},t2);
        CSS = y1.*(max(max(d)));
    elseif strcmp(trans,'complex') == 1 % Prahovani pro trans = complex
       Cr1 = d(:,:,1)/max(max(d(:,:,1)));
        Cr2 = d(:,:,2)/max(max(d(:,:,2)));
        Cr3 = d(:,:,3)/max(max(d(:,:,3)));
        Cr4 = d(:,:,4)/max(max(d(:,:,4)));
        y1 = threshold(Cr1,ee,typ{r},t2);
        y2 = threshold(Cr2,ee,typ{r},t2);
        y3 = threshold(Cr3,ee,typ{r},t2);
        y4 = threshold(Cr4,ee,typ{r},t2);
        y1(isnan(y1)) = 0;
        y2(isnan(y2)) = 0;
        y3(isnan(y2)) = 0;
        y4(isnan(y2)) = 0;
        CSS1 = y1.*(max(max(d(:,:,1))));
        CSS2 = y2.*(max(max(d(:,:,2))));
        CSS3 = y3.*(max(max(d(:,:,3))));
        CSS4 = y4.*(max(max(d(:,:,4))));
        CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
        CSS(:,:,3) = CSS3; CSS(:,:,4) = CSS4;
    end
    X = meiwa2sym(CSS,pat,trans); % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    error = mean(mean((abs(xx)-abs(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
%% Vypocty, zobrazeni...
[aaa,bbb] = sort(ssm1eee); % Serazeni hodnot MSE pro jednotlive prahovaci techniky
if strcmp(trans,'directional') == 1 % Prahovani pro trans = directiona
    Cr1 = d(:,:,1)/max(max(d(:,:,1)));
    Cr2 = d(:,:,2)/max(max(d(:,:,2)));
    y1 = threshold(Cr1,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    y2 = threshold(Cr2,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    y1(isnan(y1)) = 0;
    y2(isnan(y2)) = 0;
    CSS1 = y1.*(max(max(d(:,:,1))));
    CSS2 = y2.*(max(max(d(:,:,2))));
    CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
elseif strcmp(trans,'ortho') == 1 % Prahovani pro trans = ortho
    Cr1 = d./max(max(d));
    y1 = threshold(Cr1,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    CSS = y1.*(max(max(d)));
elseif strcmp(trans,'complex') == 1 % Prahovani pro trans = complex
    Cr1 = d(:,:,1)/max(max(d(:,:,1)));
    Cr2 = d(:,:,2)/max(max(d(:,:,2)));
    Cr3 = d(:,:,3)/max(max(d(:,:,3)));
    Cr4 = d(:,:,4)/max(max(d(:,:,4)));
    y1 = threshold(Cr1,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    y2 = threshold(Cr2,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    y3 = threshold(Cr3,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    y4 = threshold(Cr4,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1)));
    y1(isnan(y1)) = 0;
    y2(isnan(y2)) = 0;
    y3(isnan(y2)) = 0;
    y4(isnan(y2)) = 0;
    CSS1 = y1.*(max(max(d(:,:,1))));
    CSS2 = y2.*(max(max(d(:,:,2))));
    CSS3 = y3.*(max(max(d(:,:,3))));
    CSS4 = y4.*(max(max(d(:,:,4))));
    CSS(:,:,1) = CSS1; CSS(:,:,2) = CSS2;
    CSS(:,:,3) = CSS3; CSS(:,:,4) = CSS4;
end
X = meiwa2sym(CSS,pat,trans); % 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 = abs(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(abs(error))]);
disp(['SSIM = ',num2str(abs(ssm))]);
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
end
disp(' ');
disp(' ');