% Tento soubor vykonava vsechny metody pro odstranovani sumu.
% Po ukonceni vsech metod se zobrazi jejich vysledky.
% Funkce uklada obrazky ve formatu png.
% Na konci simulace se zobrazi doba potrebna na vykoneni cele funkce.
% Jednotlive metody lze najit v souborech mThreshContour.m,
% mThreshContourSD.m, mThreshFourier.m ...
clc; 
close all;    % Volitelne zavreni vsech oken Figure(X)
clear all;    % Volitelne mazani vsech promennych
tic
%% Nacteni a zasumeni obrazu
img = 'zoneplate.png';
im = imread(img);   % nacteni obrazku
im = im(:,:,1);                 % omezeni na 1 hladinu 
obb = imnoise(im,'gaussian',0.1,0.1); % "zasumeny obraz"
% obb = im2prim(obb,'y'); % omezeni velikosti obrazu pro FRAT a FRIT
% im = im2prim(im,'y');
%%%%%% Contourlet
%% Nastaveni vstupnich parametru
pfil = '5-3';
dfil = 'pkva';
nlev = [0,2,4,4];   % Rozklad na jednotlivych urovnich dekompozice
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
tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
tt2(5) = 0; % |
%% Dekompozice
yy = pdfbdec(double(obb),pfil,dfil, nlev); % dekompozice pomoci transformace Contourlet
[d, S] = pdfb2vec(yy);  % Prevedeni matic koeficientu do sloupcoveho vektoru
%% 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
    yyy = vec2pdfb(CSS, S);             % Prevedeni sloupcoveho vektoru do matic koeficientu
    X = pdfbrec(yyy,pfil,dfil);         % Rekonstrukce prahovanych koeficientu

    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)    % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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
    yyy = vec2pdfb(CSS, S);         % Prevedeni sloupcoveho vektoru do matic koeficientu
    X = pdfbrec(yyy,pfil,dfil);     % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1eee(r) = errr;
        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
y = threshold(Cr,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1))); % Provedeni prahovani pouze pro nejvhodnejsi typ prahovani
CSS = y.*(max(max(d)));     % Zpetne nastaveni puvodnich hodnot koeficientu
yyy = vec2pdfb(CSS, S);     % Prevedeni sloupcoveho vektoru do matic koeficientu
X = pdfbrec(yyy,pfil,dfil); % Rekonstrukce prahovanych koeficientu
figure; imshow(obb);        % Zobrazeni zasumeneho obrazu
imwrite(uint8(obb),[img, 'sum.png']);
xx = double(im)/256;        % xx - Vstupni hodnoty
xhat1 = (double(obb)/256);    % xhat1 - Vystupni hodnoty
errr = 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,obb);   % SSIM
wwsum{1} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwsum{2} = ['SNR  = ',num2str(snr),'dB'];
wwsum{3} = ['MSE  = ',num2str(errr)];
wwsum{4} = ['SSIM = ',num2str(ssm)];
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('Contourlet')
imwrite(uint8(X),[img, 'contourlet.png']);

xx = double(im)/256;        % xx - Vstupni hodnoty
xhat1 = (double(X)/256);    % xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
wwcont{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwcont{2} = ['type = ',typ{bbb(1)}];
wwcont{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwcont{4} = ['SNR  = ',num2str(snr),'dB'];
wwcont{5} = ['MSE  = ',num2str(errr)];
wwcont{6} = ['SSIM = ',num2str(ssm)];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwcont{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');
%%%%% ContourletSD

%% Nastaveni vstupnich parametru
Pyr_mode = 2;
dfilt = 'pkva';
smooth_func = @rcos;
nlev_SD = [1 2 4 4];
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
tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
tt2(5) = 0; % |
%% Dekompozice
yy = ContourletSDDec(double(obb), nlev_SD, Pyr_mode, smooth_func, dfilt); % dekompozice pomoci transformace ContourletSD
[d, S] = pdfb2vec(yy);  % Prevedeni matic koeficientu do sloupcoveho vektoru
%% 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
    yyy = vec2pdfb(CSS, S);             % Prevedeni sloupcoveho vektoru do matic koeficientu
    X = ContourletSDRec(yyy, Pyr_mode, smooth_func, dfilt); % Rekonstrukce prahovanych koeficientu
    
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256);     %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)    % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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  
    yyy = vec2pdfb(CSS, S);         % Prevedeni sloupcoveho vektoru do matic koeficientu
    X = ContourletSDRec(yyy, Pyr_mode, smooth_func, dfilt); % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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
y = threshold(Cr,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1))); % Provedeni prahovani pouze pro nejvhodnejsi typ prahovani
CSS = y.*(max(max(d)));     % Zpetne nastaveni puvodnich hodnot koeficientu
yyy = vec2pdfb(CSS, S);     % Prevedeni sloupcoveho vektoru do matic koeficientu
X = ContourletSDRec(yyy, Pyr_mode, smooth_func, dfilt); % Rekonstrukce prahovanych koeficientu
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('ContourletSD')
imwrite(uint8(X),[img, 'contourletSD.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu
xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (double(X)/256); %% xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
wwcontSD{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwcontSD{2} = ['type = ',typ{bbb(1)}];
wwcontSD{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwcontSD{4} = ['SNR  = ',num2str(snr),'dB'];
wwcontSD{5} = ['MSE  = ',num2str(errr)];
wwcontSD{6} = ['SSIM = ',num2str(ssm)];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwcontSD{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');

%%%% Fourier

%% 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
tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
tt2(5) = 0; % |
%% FFT2
d = fftshift(fft2(obb));
%% 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
    CSS = y.*(max(max(d)));             % Zpetne nastaveni puvodnich hodnot koeficientu   
    X = ifft2(ifftshift(CSS));          % Rekonstrukce prahovanych koeficientu


    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)    % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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)+tt(r)
    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
    CSS = y.*(max(max(d)));         % Zpetne nastaveni puvodnich hodnot koeficientu
    X = ifft2(ifftshift(CSS));      % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1eee(r) = errr;
        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
y = threshold(Cr,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1))); % Provedeni prahovani pouze pro nejvhodnejsi typ prahovani
CSS = y.*(max(max(d)));     % Zpetne nastaveni puvodnich hodnot koeficientu
X = ifft2(ifftshift(CSS));  % Rekonstrukce prahovanych koeficientu
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('Fourier');
imwrite(uint8(X),[img, 'fourier.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu
xx = double(im)/256;        %% xx - Vstupni hodnoty
xhat1 = (double(X)/256);  %% xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
wwfour{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwfour{2} = ['type = ',typ{bbb(1)}];
wwfour{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwfour{4} = ['SNR  = ',num2str(snr),'dB'];
wwfour{5} = ['MSE  = ',num2str(errr)];
wwfour{6} = ['SSIM = ',num2str(ssm)];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwfour{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');

%%%% WaveAtom

%% 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 = fwa2sym(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 = iwa2sym(CSS,pat,trans); % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean((abs(xx)-abs(xhat1)).^2));  %MSE
    if errr<ssm1(r)
        ssm1(r) = errr;
        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 = iwa2sym(CSS,pat,trans); % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean((abs(xx)-abs(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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 = iwa2sym(CSS,pat,trans); % Rekonstrukce prahovanych koeficientu
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('WaveAtom');
imwrite(uint8(X),[img, 'watom.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu

xx = double(im)/256;        %% xx - Vstupni hodnoty
xhat1 = (double(X)/256);  %% xhat1 - Vystupni hodnoty
errr = 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(abs(errr))]);
disp(['SSIM = ',num2str(abs(ssm))]);
wwwato{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwwato{2} = ['type = ',typ{bbb(1)}];
wwwato{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwwato{4} = ['SNR  = ',num2str(snr),'dB'];
wwwato{5} = ['MSE  = ',num2str(abs(errr))];
wwwato{6} = ['SSIM = ',num2str(abs(ssm))];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwwato{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');

%%%% WaveAtomME

%% 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
    errr = mean(mean((abs(xx)-abs(xhat1)).^2));  %MSE
    if errr<ssm1(r)
        ssm1(r) = errr;
        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
    errr = mean(mean((abs(xx)-abs(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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
title('WaveAtomME');
imwrite(uint8(X),[img, 'watomME.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu

xx = double(im)/256;        %% xx - Vstupni hodnoty
xhat1 = (double(X)/256);  %% xhat1 - Vystupni hodnoty
errr = 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(errr))]);
disp(['SSIM = ',num2str(abs(ssm))]);
wwwatome{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwwatome{2} = ['type = ',typ{bbb(1)}];
wwwatome{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwwatome{4} = ['SNR  = ',num2str(snr),'dB'];
wwwatome{5} = ['MSE  = ',num2str(abs(errr))];
wwwatome{6} = ['SSIM = ',num2str(abs(ssm))];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwwatome{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');

%%%%% Wavelet

%% 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 = 1:7
    [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
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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
errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        tteee(r) = ee;
        tt2(r) = t2;
        disp('_');
        disp(num2str(tteee));
        disp(num2str(ssm1eee));
        disp(wavename{q});
    end
end
end
    if mss>errr
        wname = wavename{q}; % ulozeni nazvu nejlepsi vlnky
        mss = errr;
    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
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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
errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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
title('Wavelet');
imwrite(uint8(X),[img, 'wavelet.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu

xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (double(X)/256); %% xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
disp(['wname= ',wname]);
wwwave{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwwave{2} = ['type = ',typ{bbb(1)}];
wwwave{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwwave{4} = ['SNR  = ',num2str(snr),'dB'];
wwwave{5} = ['MSE  = ',num2str(abs(errr))];
wwwave{6} = ['SSIM = ',num2str(abs(ssm))];
wwwave{8} = ['wname= ',wname];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwwave{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');

%%%% Wpak
%% Nastaveni dwtmode
dt = dwtmode;
k = 'per';
chng = 0;
if isequal(dt,k) ~= 1
    dwtmode('per');
    chng = 1;
end
%% 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
for q = 1:7
    d = wpakdec(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.
    y = wpakthresh(d,t1(e),typ{r},t2);  % Prahovani koeficientu
    X = wpakrec(y,wavename{q});
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)
        ssm1(r) = errr;
        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.
    y = wpakthresh(d,ee,typ{r},t2); % Prahovani koeficientu
    X = wpakrec(y,wavename{q});     % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        tteee(r) = ee;
        tt2(r) = t2;
        disp('_');
        disp(num2str(tteee));
        disp(num2str(ssm1eee));
        disp(wavename{q});
    end
end
end
    if mss>errr
        wname = wavename{q}; % ulozeni nazvu nejlepsi vlnky
        mss = errr;
    end
end
d = wpakdec(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.
    y = wpakthresh(d,t1(e),typ{r},t2);  % Prahovani koeficientu
    X = wpakrec(y,wname);
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)
        ssm1(r) = errr;
        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.
    y = wpakthresh(d,ee,typ{r},t2); % Prahovani koeficientu
    X = wpakrec(y,wname);     % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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 = wpakthresh(d,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1))); % Provedeni prahovani pouze pro nejvhodnejsi typ prahovani
X = wpakrec(y,wname);      % Rekonstrukce prahovanych koeficientu
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('Wpak');
imwrite(uint8(X),[img, 'wpak.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu

xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (double(X)/256); %% xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
disp(['wname= ',wname]);
wwwpak{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwwpak{2} = ['type = ',typ{bbb(1)}];
wwwpak{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwwpak{4} = ['SNR  = ',num2str(snr),'dB'];
wwwpak{5} = ['MSE  = ',num2str(abs(errr))];
wwwpak{6} = ['SSIM = ',num2str(abs(ssm))];
wwwpak{8} = ['wname= ',wname];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwwpak{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');
%% Nastaveni puvodniho dwtmode
if chng == 1;
    dwtmode(dt);
end


%%%% FRAT

obb = im2prim(obb,'y'); % Nastaveni rozmeru obrazku na prvocislo
im = im2prim(im,'y');
%% 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
tt(5)=0;    % |=> Inicializace hodnot pro nejvhodnejsi hodnotu jednotlivych prahovacich technik
tt2(5) = 0; % |
%% Dekompozice
[d,S] = frat(double(obb)); % dekompozice pomoci transformace FRAT
%% 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
    CSS = y.*(max(max(d)));             % Zpetne nastaveni puvodnich hodnot koeficientu
    X = ifrat(CSS,S);                   % Rekonstrukce prahovanych koeficientu

    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)    % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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
    CSS = y.*(max(max(d)));         % Zpetne nastaveni puvodnich hodnot koeficientu
    X = ifrat(CSS,S);               % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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
y = threshold(Cr,tteee(bbb(1)),typ{bbb(1)},tt2(bbb(1))); % Provedeni prahovani pouze pro nejvhodnejsi typ prahovani
CSS = y.*(max(max(d)));     % Zpetne nastaveni puvodnich hodnot koeficientu
X = ifrat(CSS,S);           % Rekonstrukce prahovanych koeficientu
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('FRAT');
imwrite(uint8(X),[img, 'frat.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu
xx = double(im)/256;        % xx - Vstupni hodnoty
xhat1 = (double(X)/256);    % xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
wwfrat{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwfrat{2} = ['type = ',typ{bbb(1)}];
wwfrat{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwfrat{4} = ['SNR  = ',num2str(snr),'dB'];
wwfrat{5} = ['MSE  = ',num2str(abs(errr))];
wwfrat{6} = ['SSIM = ',num2str(abs(ssm))];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwfrat{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');

%%%% FRIT

%% 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 = 1:7 % Nastaveni typu vlnek, pro ktere se provede simulace
    [d, l, m] = frit(double(obb), level, wavename{q}); % Dekompozice FRIT
    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 = ifrit(CSS, l, m, wavename{q});  % Rekonstrukce prahovanych koeficientu

    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1(r) = errr;
        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 = ifrit(CSS, l, m, wavename{q});  % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r) % Ulozeni a zobrazeni hodnot pri splneni podminky
        ssm1eee(r) = errr;
        tteee(r) = ee;
        tt2(r) = t2;
        disp('_');
        disp(num2str(tteee));
        disp(num2str(ssm1eee));
        disp(wavename{q});
    end
end
end
if errr<mss
    mss = errr;
    wname = wavename{q}; % ulozeni nazvu nejlepsi vlnky
end
end
[d, l, m] = frit(double(obb), level, wname); % dekompozice FRIT
ssm1(1:5) = 1;
%% Hruby odhad prahu pro nejlepsi vlnku
for r = 1:5
    disp(typ);
for e = 1:11
    t2 = t1(e)*2;
    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 = ifrit(CSS, l, m, wname);        % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1(r)
        ssm1(r) = errr;
        tt(r) = t1(e);
        disp(num2str(ssm1));
        disp(num2str(tt));
    end
end
end
%% Jemnyy odhad 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;
    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 = ifrit(CSS, l, m, wname);    % Rekonstrukce prahovanych koeficientu
    xx = double(im)/256; %% xx - Vstupni hodnoty
    xhat1 = (X/256); %% xhat1 - Vystupni hodnoty
    errr = mean(mean(((xx)-(xhat1)).^2));  %MSE
    if errr<ssm1eee(r)
        ssm1eee(r) = errr;
        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 = ifrit(CSS, l, m, wname);    % Rekonstrukce prahovanych koeficientu
figure; imshow(obb);        % Zobrazeni zasumeneho obrazu
title('frxt sum');
imwrite(uint8(obb),[img, 'sumfrxt.png']);
xx = double(im)/256;        % xx - Vstupni hodnoty
xhat1 = (double(obb)/256);    % xhat1 - Vystupni hodnoty
errr = 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,obb);   % SSIM
wwsumfrxt{1} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwsumfrxt{2} = ['SNR  = ',num2str(snr),'dB'];
wwsumfrxt{3} = ['MSE  = ',num2str(errr)];
wwsumfrxt{4} = ['SSIM = ',num2str(ssm)];
    
figure; imshow(uint8(X));   % Zobrazeni "odsumeneho" obrazu
title('FRIT');
imwrite(uint8(X),[img, 'frit.png']);
% figure; imshow(obb);        % Zobrazeni zasumeneho obrazu

xx = double(im)/256; %% xx - Vstupni hodnoty
xhat1 = (double(X)/256); %% xhat1 - Vystupni hodnoty
errr = 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(errr)]);
disp(['SSIM = ',num2str(ssm)]);
disp(['wname= ',wname]);
wwfrit{1} =['t1   = ',num2str(tteee(bbb(1)))];  
wwfrit{2} = ['type = ',typ{bbb(1)}];
wwfrit{3} = ['PSNR = ',num2str(PSNRdb),'dB'];
wwfrit{4} = ['SNR  = ',num2str(snr),'dB'];
wwfrit{5} = ['MSE  = ',num2str(abs(errr))];
wwfrit{6} = ['SSIM = ',num2str(abs(ssm))];
wwfrit{8} = ['wname= ',wname];
if strcmp(typ{bbb(1)},'firm')==1
    disp(['t2   = ',num2str(tt2(bbb(1)))]);
    wwfrit{7} = ['t2   = ',num2str(tt2(bbb(1)))];
end
disp(' ');
disp(' ');
disp('Contourlet');
wwcont'
disp('ContourletSD');
wwcontSD'
disp('Fourier');
wwfour'
disp('FRAT');
wwfrat'
disp('FRIT');
wwfrit'
disp('WaveAtom');
wwwato'
disp('WaveAtomME');
wwwatome'
disp('Wavelet');
wwwave'
disp('Wpak');
wwwpak'
disp('sum');
wwsum'
disp('sumfrxt');
wwsumfrxt'
disp('________________________________________');
disp(' ');
disp('Celkova doba vypoctu"');
disp(' ');
disp([num2str(floor(toc/60)),' min. ', num2str(mod(toc,60)),' sec.']);
error('KONEC!!!');