% Дослідження методів фільтраії випадкових викривлень
% 15-01-11

clear all;
% Параметри зображень
Dir_Test='D:\My_doc 02-11-10\MM_3\Test\'; % Папка для зберіганя тестових зображень
Dir_Im='D:\My_doc 02-11-10\MM_3\CFA_Natural_Stone\'; % Папка для зберігання реальних зображень
FileMask='.bmp'; Type_Real_Im=1; % Тип реального зображення: 1 - граніт Покостівського родовища ЦФА; 2 - погранична крива ЦФА; 3 - лаб. мікроскоп; 4 - тестовий об'єкт АВК
Ideal_Im_name='D:\My_doc 02-11-10\MM_3\Ideal_Image\Ideal_1982_8s_12s_RGB.bmp'; % Ідеальне зображення без випадкових викривлень (при дослідженні реальних зображень)
M=1024; N=1024;               % Розмір зображення (по верт. та гор.) тестового об'єкта вимірювань, д.т.
Y_Background=27/255;          % яскравість фону на зображенні тестового об'єкта вимірювань, д.р.
Y_Figure=227/255;             % яскравість тестового об'єкта вимірювань, д.р.
% для реального зображення задати на основі його дослідження
Threshold_Segm=(Y_Background+Y_Figure)/2; % поріг яскравості для сегментаціїї, д.р.
N_Im=3;                       % кількість зображень для моделювання та розрахунку похибок 
N_Frame_Filt=10;              % кількість нерухомих  зображень для методу фільтрації усередненням
Type_Figure=2;                % Тип зображення: 1 - тест прямокутник; 2 - тест еліпс; 3 - реальне зображеня
Boun=10;                      % ширина смуги (д.т.) на краю зображення, де відсутні тестові об'єкти вимірювань
% Параметри випадкових та динамічних викривлень
Type_Noise=1;      % Тип випадкових викривлень: 0 - некорельовані; 1 - корельовані 
AlfaX_OV=0.700;    % Параметр кореляційної функції випадкових викривлень, (д.т.)^-1 
NG=50;             % Кількість випадкових реалізацій, що додаються для отримання випадкових викривлень 
L_PK=12;           % ширина пограничної кривої, д.т. 
PSF_Size=L_PK+2;   % розмір маски фільтра для моделювання динамічних викривлень, д.т.
SNR=40;            % співвідношення сигнал/шум, дБ
PDov=0.95;         % довірча ймовірність для розрахунку похибок вимірювань
% Метод фільтрації зображень
Filtered_Method=3; % 1 - прямокутний двовимірний фільтр в частотній області
                   % 2 - прямокутний двовимірний фільтр в часовій області
                   % 3 - оптимальний двовимірний фільтр в частотній області
                   % 4 - оптимальный двовимірний фільтр в часовій області
                   % 5 - прямокутний одновимірний фільтр в реальному часі
                   % 6 - оптимальний одновимірний фільтр в реальному часі
                   % 7 - усереднення послідовності нерухомих зображень 
Filt_M_Rect=0.25;  % Нормована частота зрізу для прямокутного фільтра 
Filt_N_Rect=0.25;  % по двом просторовим координатам: миним. 0; макс. 0.5
Filt_M_Rect=round(Filt_M_Rect*M); Filt_N_Rect=round(Filt_N_Rect*N); 
Filt_N_Mask=3;     % Розмір маски фільтра в часовій області, д.т.
F_Order_1D=1;      % порядок одновимірного фільтра в реальному часі
% Параметри виведення результатів розрахунку
Im_Show=0;         % виведення зображень
Delta_Show=1;      % розрахунок та виведення величини похибок вимірювань
Im_Write=0;        % запис зображень на диск
% Параметри тестових об'єктів вимірювань (прямокутника та еліпса)
X_Center=round(N/2); Y_Center=round(M/2);         % центр об'єкта вимірювань
RectXMax=round(M/sqrt(2)-2*Boun-2*L_PK);          % Велика сторона прямокутника
RectYMin=round((M/sqrt(2)-2*Boun-2*L_PK)/10);     % Мала сторона прямокутника
EllipsAMax=round((M/sqrt(2)-2*Boun-2*L_PK)/2);    % Велика вісь еліпса
EllipsBMin=round((M/sqrt(2)-2*Boun-2*L_PK)/2/3);  % Мала вісь еліпса

% Масиви для запису результатів розрахунку
A_SKO_Im_M=zeros(N_Im,1); A_SNR_dB_Im_M=zeros(N_Im,1); A_PSNR_dB_Im_M=zeros(N_Im,1);
A_SKO_Im_R=zeros(N_Im,1); A_SNR_dB_Im_R=zeros(N_Im,1); A_PSNR_dB_Im_R=zeros(N_Im,1);
B_M_XOR=zeros(N_Im,1); B_R_XOR=zeros(N_Im,1); 
S_Figure=zeros(N_Im,1); S_Figure_M=zeros(N_Im,1); S_Figure_R=zeros(N_Im,1); 
Centroid_X=zeros(N_Im,1); Centroid_X_M=zeros(N_Im,1); Centroid_X_R=zeros(N_Im,1); 
Centroid_Y=zeros(N_Im,1); Centroid_Y_M=zeros(N_Im,1); Centroid_Y_R=zeros(N_Im,1); 
W_Figure=zeros(N_Im,1); W_Figure_M=zeros(N_Im,1); W_Figure_R=zeros(N_Im,1); 
H_Figure=zeros(N_Im,1); H_Figure_M=zeros(N_Im,1); H_Figure_R=zeros(N_Im,1); 
ED_Figure=zeros(N_Im,1); ED_Figure_M=zeros(N_Im,1); ED_Figure_R=zeros(N_Im,1);
% rand('state',sum(100*clock)); randn('state',sum(100*clock));
rand('state',0); randn('state',0);

% ФОРМУВАННЯ ПЕРЕЛІКУ ФАЙЛІВ В ПАПЦІ З ДВОВИІМРНОЮ ІНФОРМАЦІЄЮ
ListFileInDir=dir(Dir_Im); 
NumberOfFileInDir=length(ListFileInDir(:,1)); 
NumberOfFileToMask=0;
ListFileToMask=[];
for i=1:NumberOfFileInDir,
    if strfind(ListFileInDir(i).name, upper(FileMask))>=1
        NumberOfFileToMask=NumberOfFileToMask+1;
        ListFileToMask(NumberOfFileToMask).name=ListFileInDir(i).name;
    end;
    if strfind(ListFileInDir(i).name, lower(FileMask))>=1
        NumberOfFileToMask=NumberOfFileToMask+1;
        ListFileToMask(NumberOfFileToMask).name=ListFileInDir(i).name;
    end;
end;
% Результат пошуку файлів:
% ListFileToMask - перелік файлів, що відповідають масці для пошуку
% NumberOfFileToMask - кількість файлів, що відповідають масці для пошуку

% ----- Цикл моделювання та обробки зображень, а також визначення похибок вимырювань -----
h_waitbar = waitbar(0,'Obrobka zobrazhen');
for i_Im=1:N_Im
waitbar(i_Im/N_Im,h_waitbar);

% Масиви для початкових даних та запису результатів розрахунку
Im_Frame=zeros(M,N); % Початкове зображення 
Im_Frame_M=zeros(M,N);     % Зображення з випадковими викривленнями
Im_Frame_R=zeros(M,N);     % Зображення після фільтрації випадкових виркивлень

% Формування початкового зображення
switch Type_Figure
    case 1 % формування тестового прямокутника з випадковими геометричними параметрами
        RectX=round(RectXMax-(RectXMax-RectYMin)*rand);
        RectY=round(RectXMax-(RectXMax-RectYMin)*rand);
        RectTheta=180*rand;
        Im_Frame_J=zeros(RectY, RectX)+Y_Figure;
        Im_Frame_J = imrotate(Im_Frame_J,RectTheta,'bicubic','loose');
        [M_J N_J]=size(Im_Frame_J);
        Im_Frame(round(Y_Center-M_J/2+1):round(Y_Center+M_J/2), round(X_Center-N_J/2+1):round(X_Center+N_J/2))=Im_Frame_J;
        Im_Frame(find(Im_Frame==0))=Y_Background;
        Im_Frame=mat2gray(Im_Frame, [0 1]);
    case 2 % формування тестового еліпса з випадковими геометричними параметрами
        EllipsA=round(EllipsAMax-(EllipsAMax-EllipsBMin)*rand);
        EllipsB=round(EllipsAMax-(EllipsAMax-EllipsBMin)*rand);
        EllipsTheta=180*rand;
        Im_Frame_J=zeros(EllipsB*2, EllipsA*2);
        for j=1:EllipsB*2
            for i=1:EllipsA*2
                if abs(i-EllipsA)<=EllipsA*sqrt(1-(EllipsB-j)^2/EllipsB^2)
                    Im_Frame_J(j,i)=Y_Figure;
                end;
            end;
        end;
        Im_Frame_J = imrotate(Im_Frame_J,EllipsTheta,'bicubic','loose');
        [M_J N_J]=size(Im_Frame_J);
        Im_Frame(round(Y_Center-M_J/2+1):round(Y_Center+M_J/2), round(X_Center-N_J/2+1):round(X_Center+N_J/2))=Im_Frame_J;
        Im_Frame(find(Im_Frame<Y_Background))=Y_Background;
        Im_Frame=mat2gray(Im_Frame, [0 1]);
    case 3 % завантаження реального зображення
        if i_Im<=NumberOfFileToMask
            Im_Frame=rgb2gray(im2double(imread(Ideal_Im_name)));
            Im_Frame_M=rgb2gray(im2double(imread(strcat(Dir_Im,ListFileToMask(i_Im).name))));
            M=size(Im_Frame,1); N=size(Im_Frame,2); 
            Im_Frame_R=zeros(M,N); 
        else
            N_Im=i_Im-1; break;
        end;
end; % switch Type_Figure

% Моделювання випадкових та динамічних викривлень тестових зображень
if Type_Figure<3
    PSF_DVC=fspecial('gaussian',PSF_Size,L_PK/6);
    Im_Frame=imfilter(Im_Frame, PSF_DVC, 'circular');
    switch Type_Noise 
        case 0
            Disp_Gaussian=1/10^(2*SNR/20);
            Im_Frame_M=imnoise(Im_Frame,'gaussian',0,Disp_Gaussian);
        case 1
            gamma1=rand(NG,1); gamma2=rand(NG,1); gamma4=rand(NG,1)*2*pi; 
            psbp=rand(NG,1); 
            for lll=1:NG 
                if psbp(lll)>=0.5 psbp(lll)=1; else psbp(lll)=-1; end;
            end; 
            ksi_k=psbp.*sqrt(-log(gamma1));
            e1=cos(gamma4); e2=sin(gamma4); % орт x орт y
            eta_mod=AlfaX_OV.*sqrt(1./gamma2.^2-1);
            eta1_k=eta_mod.*e1; eta2_k=eta_mod.*e2;     
            h2 = waitbar(0,strcat('Wait ... generated Im Noise ',num2str(i_Im)));
            Im_Noise=zeros(M,N);
            for lll=1:NG
                waitbar(lll/NG,h2);
                for jjj=1:M for iii=1:N
                    Im_Noise(jjj,iii)=Im_Noise(jjj,iii)+sqrt(2)*ksi_k(lll)*sin(eta1_k(lll)*iii+eta2_k(lll)*jjj+pi/4);
                end; end;
            end;
            close(h2);       
            Im_Noise=1/10^(SNR/20)/sqrt(NG)*Im_Noise; 
            Im_Frame_M=Im_Frame+Im_Noise;
    end; % switch Type_Noise 
end; % if Type_Figure<3

% фільтрация
switch Filtered_Method
    case 1 % 1 - прямокутний двовимірний фільтр в частотній області
        Snoise2=abs(fft2(Im_Frame_M-Im_Frame, M*2, N*2)).^2;
        Simage2=abs(fft2(Im_Frame, M*2, N*2)).^2;
        HF=zeros(M*2,N*2);
        HF=fftshift(Simage2./(Simage2+Snoise2));
        h17=fspecial('average',17);
        HF=imfilter(HF,h17);
        HF=zeros(M*2,N*2);
        HF((M-2*Filt_M_Rect+1):(M+2*Filt_M_Rect),(N-2*Filt_N_Rect+1):(N+2*Filt_N_Rect))=1;
        Fourie_Im=fftshift(fft2(Im_Frame_M, M*2, N*2));
        iFourie_Im=real(ifft2(fftshift(HF.*Fourie_Im)));
        Im_Frame_R=iFourie_Im(1:M, 1:N);
    case 2 % 2 - прямокутний двовимірний фільтр в часовій області
        HF=zeros(M,N);
        HF((round(M/2)-Filt_M_Rect+1):(round(M/2)+Filt_M_Rect),...
           (round(N/2)-Filt_N_Rect+1):(round(N/2)+Filt_N_Rect))=1;
        [f1 f2]=freqspace([M N], 'meshgrid');
        hf=real(fsamp2(f1,f2,HF,[Filt_N_Mask Filt_N_Mask])); hf=hf/sum(sum(hf)); 
        Im_Frame_R=imfilter(Im_Frame_M,hf);
    case 3 % 3 - оптимальний двовимірний фільтр в частотній області
        % варіант з розширенням массиву
        Snoise2=abs(fft2(Im_Frame_M-Im_Frame, M*2, N*2)).^2;
        Simage2=abs(fft2(Im_Frame, M*2, N*2)).^2;
        HF=zeros(M*2,N*2);
        HF=fftshift(Simage2./(Simage2+Snoise2));
        h17=fspecial('average',17);
        HF=imfilter(HF,h17);
        Fourie_Im=fftshift(fft2(Im_Frame_M, M*2, N*2));
        iFourie_Im=real(ifft2(fftshift(HF.*Fourie_Im)));
        Im_Frame_R=iFourie_Im(1:M, 1:N);
    case 4 % 4 - оптимальний двовимірний фільтр в часовій області
        Snoise=abs(fft2(Im_Frame_M-Im_Frame)).^2;
        Simage=abs(fft2(Im_Frame)).^2;
        HF=zeros(M,N);
        HF=fftshift(Simage./(Simage+Snoise));
        h17=fspecial('average',17);
        HF=imfilter(HF,h17);
        [f1 f2]=freqspace([M N], 'meshgrid');
        hf=real(fsamp2(f1,f2,HF,[Filt_N_Mask Filt_N_Mask])); hf=hf/sum(sum(hf)); 
        Im_Frame_R=imfilter(Im_Frame_M,hf);
    case 5 % 5 - прямокутний одновимірний фільтр в реальному часі
        F_1D=linspace(0,1,round(N/2));
        HF_1D=zeros(1,round(N/2));
        HF_1D(1:Filt_N_Rect)=1;
        [B_1D,A_1D]=yulewalk(F_Order_1D,F_1D,HF_1D);
        for jjj=1:M
            Im_Frame_R(jjj,:) = filter(B_1D, A_1D, Im_Frame_M(jjj,:));
        end;
    case 6 % 6 - оптимальний одновимірний фільтр в реальному часі
        F_1D=linspace(0,1,round(N/2));
        Snoise1=abs(fft(Im_Frame_M(round(M/2),:)-Im_Frame(round(M/2),:), N)).^2;
        Simage1=abs(fft(Im_Frame(round(M/2),:), N)).^2;
        HF_1D=Simage1./(Simage1+Snoise1);
        h17=fspecial('average',17);
        HF_1D=imfilter(HF_1D(1:round(N/2)),h17);
        [B_1D,A_1D]=yulewalk(F_Order_1D,F_1D,HF_1D);
        B_1D=B_1D*sum(A_1D)/sum(B_1D);
        for jjj=1:M
            Im_Frame_R(jjj,:) = filter(B_1D, A_1D, Im_Frame_M(jjj,:));
        end;
    case 7 % 7 - усереднення послідовності нерухомих зображень 
        Im_Frame_R=zeros(M,N);
        h1_waitbar = waitbar(0,'Obrobka neruchomych zobrazhen');
        if Type_Figure<3
            for iii=1:N_Frame_Filt
                waitbar(iii/N_Frame_Filt,h1_waitbar);
                switch Type_Noise 
                    case 0
                        Disp_Gaussian=1/10^(2*SNR/20);
                        Im_Frame_M=imnoise(Im_Frame,'gaussian',0,Disp_Gaussian);
                    case 1
                        gamma1=rand(NG,1); gamma2=rand(NG,1); gamma4=rand(NG,1)*2*pi; 
                        psbp=rand(NG,1); 
                        for lll=1:NG 
                            if psbp(lll)>=0.5 psbp(lll)=1; else psbp(lll)=-1; end;
                        end; 
                        ksi_k=psbp.*sqrt(-log(gamma1));
                        e1=cos(gamma4); e2=sin(gamma4); % орт x орт y
                        eta_mod=AlfaX_OV.*sqrt(1./gamma2.^2-1);
                        eta1_k=eta_mod.*e1; eta2_k=eta_mod.*e2;     
                        h2 = waitbar(0,strcat('Wait ... generated Im Noise ',num2str(iii)));
                        Im_Noise=zeros(M,N);
                        for lll=1:NG
                            waitbar(lll/NG,h2);
                            for jjj=1:M for iii=1:N
                                Im_Noise(jjj,iii)=Im_Noise(jjj,iii)+sqrt(2)*ksi_k(lll)*sin(eta1_k(lll)*iii+eta2_k(lll)*jjj+pi/4);
                            end; end;
                        end;
                        close(h2);       
                        Im_Noise=1/10^(SNR/20)/sqrt(NG)*Im_Noise; 
                        Im_Frame_M=Im_Frame+Im_Noise;
                end; % switch Type_Noise 
                Im_Frame_R=Im_Frame_R+Im_Frame_M;
            end;
        else
            for iii=1:N_Frame_Filt
                waitbar(iii/N_Frame_Filt,h1_waitbar);
                if iii<=NumberOfFileToMask
                    Im_Frame_M=rgb2gray(im2double(imread(strcat(Dir_Im,ListFileToMask(iii).name))));
                    Im_Frame_R=Im_Frame_R+Im_Frame_M;
                else
                    N_Frame_Filt=iii-1; break;
                end;
            end;
        end;
        Im_Frame_R=Im_Frame_R/N_Frame_Filt;
        close(h1_waitbar);
end; % switch Filtered_Method

% Сегментація
Pixel_Value=zeros(M,N); Pixel_Value_M=zeros(M,N); Pixel_Value_R=zeros(M,N);
if Y_Figure<Threshold_Segm 
    Pixel_Value(1+Boun:M-Boun,1+Boun:N-Boun)=roicolor(Im_Frame(1+Boun:M-Boun,1+Boun:N-Boun), 0, Threshold_Segm); 
    Pixel_Value_M(1+Boun:M-Boun,1+Boun:N-Boun)=roicolor(Im_Frame_M(1+Boun:M-Boun,1+Boun:N-Boun), 0, Threshold_Segm); 
    Pixel_Value_R(1+Boun:M-Boun,1+Boun:N-Boun)=roicolor(Im_Frame_R(1+Boun:M-Boun,1+Boun:N-Boun), 0, Threshold_Segm); 
else
    Pixel_Value(1+Boun:M-Boun,1+Boun:N-Boun)=roicolor(Im_Frame(1+Boun:M-Boun,1+Boun:N-Boun), Threshold_Segm, 1); 
    Pixel_Value_M(1+Boun:M-Boun,1+Boun:N-Boun)=roicolor(Im_Frame_M(1+Boun:M-Boun,1+Boun:N-Boun), Threshold_Segm, 1); 
    Pixel_Value_R(1+Boun:M-Boun,1+Boun:N-Boun)=roicolor(Im_Frame_R(1+Boun:M-Boun,1+Boun:N-Boun), Threshold_Segm, 1); 
end;

% Пошук об'єктів 
[L num]=bwlabel(Pixel_Value, 8); 
[L_M num_M]=bwlabel(Pixel_Value_M, 8); 
[L_R num_R]=bwlabel(Pixel_Value_R, 8); 
% Обчислення геометричних параметрів об'єктів: площадей и координат центрів мас 
feats=imfeature(L, 'Area', 'Centroid', 8); 
feats_M=imfeature(L_M, 'Area', 'Centroid', 8); 
feats_R=imfeature(L_R, 'Area', 'Centroid', 8); 
% Перезапис значень параметрів з масива структур feats в окремі масиви 
Areas=zeros(num); CentX=zeros(num); CentY=zeros(num); 
Areas_M=zeros(num_M); CentX_M=zeros(num_M); CentY_M=zeros(num_M); 
Areas_R=zeros(num_R); CentX_R=zeros(num_R); CentY_R=zeros(num_R); 
for i=1:1:num 
    Areas(i)=feats(i).Area; CentX(i)=feats(i).Centroid(1); CentY(i)=feats(i).Centroid(2); 
end; 
for i=1:1:num_M 
    Areas_M(i)=feats_M(i).Area; CentX_M(i)=feats_M(i).Centroid(1); CentY_M(i)=feats_M(i).Centroid(2); 
end; 
for i=1:1:num_R 
    Areas_R(i)=feats_R(i).Area; CentX_R(i)=feats_R(i).Centroid(1); CentY_R(i)=feats_R(i).Centroid(2); 
end; 
% Формування бінарного зображення з об'єктами, площа яких більш ніж 1000 (д.т.)^2, 
% а центр мас знаходиться на відстані більш ніж 50 д.т. від краю зображення 
[rows cols]=size(Pixel_Value); 
idx=find(Areas>1000&CentX>50&CentY>50&CentX<(cols-50)&CentY<(rows-50)); 
idx_M=find(Areas_M>1000&CentX_M>50&CentY_M>50&CentX_M<(cols-50)&CentY_M<(rows-50)); 
idx_R=find(Areas_R>1000&CentX_R>50&CentY_R>50&CentX_R<(cols-50)&CentY_R<(rows-50)); 
Pixel_Value=ismember(L, idx); 
Pixel_Value_M=ismember(L_M, idx_M); 
Pixel_Value_R=ismember(L_R, idx_R); 
% Знову пошук об'єктів 
[L num]=bwlabel(Pixel_Value, 8); 
[L_M num_M]=bwlabel(Pixel_Value_M, 8); 
[L_R num_R]=bwlabel(Pixel_Value_R, 8); 
% Обчислення геометричних параметрів об'єктів
feats=imfeature(L, 'Area', 'Centroid', 'BoundingBox', 'EquivDiameter', 8); 
feats_M=imfeature(L_M, 'Area', 'Centroid', 'BoundingBox', 'EquivDiameter', 8); 
feats_R=imfeature(L_R, 'Area', 'Centroid', 'BoundingBox', 'EquivDiameter', 8); 

% Обчислення похибок згідно визначених критеріїв
% Амплітудні критерії
% СКЗ амплітудної похибки
A_SKO_Im_M(i_Im)=255.0*sqrt(sum(sum((Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun)-Im_Frame_M(Boun+1:M-Boun,Boun+1:N-Boun)).^2))/((M-2*Boun)*(N-2*Boun)));
A_SKO_Im_R(i_Im)=255.0*sqrt(sum(sum((Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun)-Im_Frame_R(Boun+1:M-Boun,Boun+1:N-Boun)).^2))/((M-2*Boun)*(N-2*Boun)));
% Співвідношення сигнал/шум
A_SNR_dB_Im_M(i_Im)=10*log10(sum(sum(Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun).^2))/sum(sum((Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun)-Im_Frame_M(Boun+1:M-Boun,Boun+1:N-Boun)).^2)));
A_SNR_dB_Im_R(i_Im)=10*log10(sum(sum(Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun).^2))/sum(sum((Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun)-Im_Frame_R(Boun+1:M-Boun,Boun+1:N-Boun)).^2)));
% Пікове співвідношення сигнал/шум
A_PSNR_dB_Im_M(i_Im)=10*log10(1.0*(M-2*Boun)*(N-2*Boun)/sum(sum((Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun)-Im_Frame_M(Boun+1:M-Boun,Boun+1:N-Boun)).^2)));
A_PSNR_dB_Im_R(i_Im)=10*log10(1.0*(M-2*Boun)*(N-2*Boun)/sum(sum((Im_Frame(Boun+1:M-Boun,Boun+1:N-Boun)-Im_Frame_R(Boun+1:M-Boun,Boun+1:N-Boun)).^2)));

% Геометричні параметри: площа, центр масс, шириина и высота
S_Figure(i_Im)=feats.Area;
Centroid_X(i_Im)=feats.Centroid(1); Centroid_Y(i_Im)=feats.Centroid(2); 
W_Figure(i_Im)=feats.BoundingBox(3); H_Figure(i_Im)=feats.BoundingBox(4); 
S_Figure_M(i_Im)=feats_M.Area;
Centroid_X_M(i_Im)=feats_M.Centroid(1); Centroid_Y_M(i_Im)=feats_M.Centroid(2); 
W_Figure_M(i_Im)=feats_M.BoundingBox(3); H_Figure_M(i_Im)=feats_M.BoundingBox(4); 
S_Figure_R(i_Im)=feats_R.Area;
Centroid_X_R(i_Im)=feats_R.Centroid(1); Centroid_Y_R(i_Im)=feats_R.Centroid(2); 
W_Figure_R(i_Im)=feats_R.BoundingBox(3); H_Figure_R(i_Im)=feats_R.BoundingBox(4); 
ED_Figure(i_Im)=feats.EquivDiameter; ED_Figure_M(i_Im)=feats_M.EquivDiameter; ED_Figure_R(i_Im)=feats_R.EquivDiameter; 

% порівняння контурів об'єктів вимірювань 
B_M=xor(Pixel_Value, Pixel_Value_M);
B_R=xor(Pixel_Value, Pixel_Value_R);
if Im_Show==1
    figure; subplot(2, 3, 1); imshow(Pixel_Value); title('Segm. pochatkovogo zobrazhennya');
    subplot(2, 3, 2); imshow(Pixel_Value_M); title('Segm. zobrazhennya z vykryvlennyamy');
    subplot(2, 3, 5); imshow(Pixel_Value_R); title('Segm. filtr. zobrazhennya');
    subplot(2, 3, 3); imshow(B_M); title('Vidminnist pochatkovogo i vykr.');
    subplot(2, 3, 6); imshow(B_R); title('Vidminnist pochatkovogo i filtr.');
end;
B_M_XOR(i_Im)=sum(sum(B_M(Boun+1:M-Boun,Boun+1:N-Boun)));
B_R_XOR(i_Im)=sum(sum(B_R(Boun+1:M-Boun,Boun+1:N-Boun)));

% запис зображень на диск
if Im_Write==1
    Name_Frame=strcat(Dir_Test,'Test_',num2str(i_Im,'%d'),'.bmp');
    imwrite(Im_Frame,Name_Frame);
    Name_Frame=strcat(Dir_Test,'Test_Noise_',num2str(i_Im,'%d'),'.bmp');
    imwrite(Im_Frame_M,Name_Frame);
    Name_Frame=strcat(Dir_Test,'Filtered_',num2str(i_Im,'%d'),'.bmp');
    imwrite(Im_Frame_R,Name_Frame);
    Name_Frame=strcat(Dir_Test,'Segm_Test_',num2str(i_Im,'%d'),'.bmp');
    imwrite(Pixel_Value,Name_Frame);
    Name_Frame=strcat(Dir_Test,'Segm_Test_Noise_',num2str(i_Im,'%d'),'.bmp');
    imwrite(Pixel_Value_M,Name_Frame);
    Name_Frame=strcat(Dir_Test,'Segm_Filtered_',num2str(i_Im,'%d'),'.bmp');
    imwrite(Pixel_Value_R,Name_Frame);
end;
% Виведення зображень
if Im_Show==1
    figure; subplot(2, 2, 1); imshow(Im_Frame); title('Pochatkove zobrazhennya');
    subplot(2, 2, 2); imshow(Im_Frame_M); title('Zobrazhennya z vykryvlennyamy');
    subplot(2, 2, 3); imshow(Im_Frame_R); title('Filtrovane zobrazhennya');
    subplot(2, 2, 4); imshow(Pixel_Value_R); title('Segm. filtr. zobrazhennya');
end;

end; % for i_Im=1:N_Im
close(h_waitbar);

% Обчислення похибок згідно визначених критеріїв
if Delta_Show==1
    
    % Амплітудні критерії
    Max_A_SKO_Im_M=max(A_SKO_Im_M); Max_A_SKO_Im_R=max(A_SKO_Im_R); 
    Mean_A_SKO_Im_M=mean(A_SKO_Im_M); Mean_A_SKO_Im_R=mean(A_SKO_Im_R); 
    Min_A_SKO_Im_M=min(A_SKO_Im_M); Min_A_SKO_Im_R=min(A_SKO_Im_R); 
    Delta_Dov_M=sort(A_SKO_Im_M, 'descend'); Delta_Dov_R=sort(A_SKO_Im_R, 'descend'); 
    MaxDov_A_SKO_Im_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_A_SKO_Im_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    
    Max_A_SNR_dB_Im_M=max(A_SNR_dB_Im_M); Max_A_SNR_dB_Im_R=max(A_SNR_dB_Im_R); 
    Mean_A_SNR_dB_Im_M=mean(A_SNR_dB_Im_M); Mean_A_SNR_dB_Im_R=mean(A_SNR_dB_Im_R); 
    Min_A_SNR_dB_Im_M=min(A_SNR_dB_Im_M); Min_A_SNR_dB_Im_R=min(A_SNR_dB_Im_R); 
    Delta_Dov_M=sort(A_SNR_dB_Im_M, 'descend'); Delta_Dov_R=sort(A_SNR_dB_Im_R, 'descend'); 
    MaxDov_A_SNR_dB_Im_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_A_SNR_dB_Im_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    
    Max_A_PSNR_dB_Im_M=max(A_PSNR_dB_Im_M); Max_A_PSNR_dB_Im_R=max(A_PSNR_dB_Im_R); 
    Mean_A_PSNR_dB_Im_M=mean(A_PSNR_dB_Im_M); Mean_A_PSNR_dB_Im_R=mean(A_PSNR_dB_Im_R); 
    Min_A_PSNR_dB_Im_M=min(A_PSNR_dB_Im_M); Min_A_PSNR_dB_Im_R=min(A_PSNR_dB_Im_R); 
    Delta_Dov_M=sort(A_PSNR_dB_Im_M, 'descend'); Delta_Dov_R=sort(A_PSNR_dB_Im_R, 'descend'); 
    MaxDov_A_PSNR_dB_Im_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_A_PSNR_dB_Im_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
   
    % площа і центр мас
    Delta_S_Figure_M=(S_Figure_M-S_Figure); D1_S_Figure_M=Delta_S_Figure_M.*100./S_Figure;
    Delta_S_Figure_R=(S_Figure_R-S_Figure); D1_S_Figure_R=Delta_S_Figure_R.*100./S_Figure;
    Delta_Centroid_X_M=Centroid_X_M-Centroid_X; D1_Centroid_X_M=Delta_Centroid_X_M.*100./Centroid_X;
    Delta_Centroid_X_R=Centroid_X_R-Centroid_X; D1_Centroid_X_R=Delta_Centroid_X_R.*100./Centroid_X;
    Delta_Centroid_Y_M=Centroid_Y_M-Centroid_Y; D1_Centroid_Y_M=Delta_Centroid_Y_M.*100./Centroid_Y;
    Delta_Centroid_Y_R=Centroid_Y_R-Centroid_Y; D1_Centroid_Y_R=Delta_Centroid_Y_R.*100./Centroid_Y;

    Max_Delta_S_Figure_M=max(abs(Delta_S_Figure_M)); Max_Delta_S_Figure_R=max(abs(Delta_S_Figure_R)); 
    Delta_Dov_M=sort(abs(Delta_S_Figure_M), 'descend'); Delta_Dov_R=sort(abs(Delta_S_Figure_R), 'descend'); 
    MaxDov_Delta_S_Figure_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_Delta_S_Figure_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_Delta_S_Figure_M=mean(Delta_S_Figure_M); Mean_Delta_S_Figure_R=mean(Delta_S_Figure_R); 
    Std_Delta_S_Figure_M=std(Delta_S_Figure_M,1); Std_Delta_S_Figure_R=std(Delta_S_Figure_R,1); 

    Max_Delta_Centroid_M=max(abs([Delta_Centroid_X_M; Delta_Centroid_Y_M])); 
    Max_Delta_Centroid_R=max(abs([Delta_Centroid_X_R; Delta_Centroid_Y_R])); 
    Delta_Dov_M=sort(abs([Delta_Centroid_X_M; Delta_Centroid_Y_M]), 'descend'); 
    Delta_Dov_R=sort(abs([Delta_Centroid_X_R; Delta_Centroid_Y_R]), 'descend'); 
    MaxDov_Delta_Centroid_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_Delta_Centroid_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_Delta_Centroid_M=mean([Delta_Centroid_X_M; Delta_Centroid_Y_M]); 
    Mean_Delta_Centroid_R=mean([Delta_Centroid_X_R; Delta_Centroid_Y_R]); 
    Std_Delta_Centroid_M=std([Delta_Centroid_X_M; Delta_Centroid_Y_M],1); 
    Std_Delta_Centroid_R=std([Delta_Centroid_X_R; Delta_Centroid_Y_R]); 

    Max_D1_S_Figure_M=max(abs(D1_S_Figure_M)); Max_D1_S_Figure_R=max(abs(D1_S_Figure_R)); 
    Delta_Dov_M=sort(abs(D1_S_Figure_M), 'descend'); Delta_Dov_R=sort(abs(D1_S_Figure_R), 'descend'); 
    MaxDov_D1_S_Figure_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_D1_S_Figure_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_D1_S_Figure_M=mean(D1_S_Figure_M); Mean_D1_S_Figure_R=mean(D1_S_Figure_R); 
    Std_D1_S_Figure_M=std(D1_S_Figure_M,1); Std_D1_S_Figure_R=std(D1_S_Figure_R,1); 

    Max_D1_Centroid_M=max(abs([D1_Centroid_X_M; D1_Centroid_Y_M])); 
    Max_D1_Centroid_R=max(abs([D1_Centroid_X_R; D1_Centroid_Y_R])); 
    Delta_Dov_M=sort(abs([D1_Centroid_X_M; D1_Centroid_Y_M]), 'descend'); 
    Delta_Dov_R=sort(abs([D1_Centroid_X_R; D1_Centroid_Y_R]), 'descend'); 
    MaxDov_D1_Centroid_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_D1_Centroid_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_D1_Centroid_M=mean([D1_Centroid_X_M; D1_Centroid_Y_M]); 
    Mean_D1_Centroid_R=mean([D1_Centroid_X_R; D1_Centroid_Y_R]); 
    Std_D1_Centroid_M=std([D1_Centroid_X_M; D1_Centroid_Y_M],1); 
    Std_D1_Centroid_R=std([D1_Centroid_X_R; D1_Centroid_Y_R]); 

    % ширина і висота
    Delta_W_Figure_M=(W_Figure_M-W_Figure); D1_W_Figure_M=Delta_W_Figure_M.*100./W_Figure;
    Delta_W_Figure_R=(W_Figure_R-W_Figure); D1_W_Figure_R=Delta_W_Figure_R.*100./W_Figure;
    Delta_H_Figure_M=(H_Figure_M-H_Figure); D1_H_Figure_M=Delta_H_Figure_M.*100./H_Figure;
    Delta_H_Figure_R=(H_Figure_R-H_Figure); D1_H_Figure_R=Delta_H_Figure_R.*100./H_Figure;

    Max_Delta_L_Figure_M=max(abs([Delta_W_Figure_M; Delta_H_Figure_M])); 
    Max_Delta_L_Figure_R=max(abs([Delta_W_Figure_R; Delta_H_Figure_R])); 
    Delta_Dov_M=sort(abs([Delta_W_Figure_M; Delta_H_Figure_M]), 'descend'); 
    Delta_Dov_R=sort(abs([Delta_W_Figure_R; Delta_H_Figure_R]), 'descend'); 
    MaxDov_Delta_L_Figure_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_Delta_L_Figure_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_Delta_L_Figure_M=mean([Delta_W_Figure_M; Delta_H_Figure_M]); 
    Mean_Delta_L_Figure_R=mean([Delta_W_Figure_R; Delta_H_Figure_R]); 
    Std_Delta_L_Figure_M=std([Delta_W_Figure_M; Delta_H_Figure_M],1); 
    Std_Delta_L_Figure_R=std([Delta_W_Figure_R; Delta_H_Figure_R],1); 

    Max_D1_L_Figure_M=max(abs([D1_W_Figure_M; D1_H_Figure_M])); 
    Max_D1_L_Figure_R=max(abs([D1_W_Figure_R; D1_H_Figure_R])); 
    Delta_Dov_M=sort(abs([D1_W_Figure_M; D1_H_Figure_M]), 'descend'); 
    Delta_Dov_R=sort(abs([D1_W_Figure_R; D1_H_Figure_R]), 'descend'); 
    MaxDov_D1_L_Figure_M=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_D1_L_Figure_R=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_D1_L_Figure_M=mean([D1_W_Figure_M; D1_H_Figure_M]); 
    Mean_D1_L_Figure_R=mean([D1_W_Figure_R; D1_H_Figure_R]); 
    Std_D1_L_Figure_M=std([D1_W_Figure_M; D1_H_Figure_M],1); 
    Std_D1_L_Figure_R=std([D1_W_Figure_R; D1_H_Figure_R],1); 

    % похибки контуров
    B_M_XOR=B_M_XOR./(ED_Figure.*pi); B_R_XOR=B_R_XOR./(ED_Figure.*pi); 
    Max_B_M_XOR=max(abs(B_M_XOR)); Max_B_R_XOR=max(abs(B_R_XOR)); 
    Delta_Dov_M=sort(abs(B_M_XOR), 'descend'); Delta_Dov_R=sort(abs(B_R_XOR), 'descend'); 
    MaxDov_B_M_XOR=Delta_Dov_M(1+round(N_Im*(1-PDov))); MaxDov_B_R_XOR=Delta_Dov_R(1+round(N_Im*(1-PDov)));
    Mean_B_M_XOR=mean(B_M_XOR); Mean_B_R_XOR=mean(B_R_XOR); 
    Min_B_M_XOR=min(B_M_XOR); Min_B_R_XOR=min(B_R_XOR); 

end; % {Delta_Show}

fprintf(1,'----- Дослідження точності методів фільтрації двовимірної інформації -----\n');
fprintf(1,'Розмір зображення, д.т. %d x %d\n', M, N); 
fprintf(1,'Об"єкт вимірювань: 1 - прямокутник; 2 - еліпс; 3 - реальний %d\n', Type_Figure);
fprintf(1,'Кількість зображень для дослідження %d\n', N_Im);
fprintf(1,'Ширина пограничної кривої, д.т. %d\n', L_PK);
fprintf(1,'Розмір маски фільтра, д.т. %d\n', Filt_N_Mask);
fprintf(1,'Співвідношення сигнал/шум, дБ %d\n', SNR);
fprintf(1,'Метод фільтрації %d\n', Filtered_Method);
if Delta_Show==1
    fprintf(1,'\nАмпл похибка                  без фільт. з фільт.\n');
    fprintf(1,'Середньокв. зн., д.р.     %10.3f %10.3f\n', Mean_A_SKO_Im_M, Mean_A_SKO_Im_R); 
    fprintf(1,'Співвід. с/ш, дБ          %10.3f %10.3f\n', Mean_A_SNR_dB_Im_M, Mean_A_SNR_dB_Im_R); 
    fprintf(1,'Пікове співвід. с/ш, дБ   %10.3f %10.3f\n', Mean_A_PSNR_dB_Im_M, Mean_A_PSNR_dB_Im_R); 

    fprintf(1,'Абсолютні похибки геометричних параметрів, д.т.\n');
    fprintf(1,'Без фільт       площа   центр мас   лінійні розміри\n');
    fprintf(1,'Максимальна  %10.3f %7.3f %7.3f\n', Max_Delta_S_Figure_M, Max_Delta_Centroid_M, Max_Delta_L_Figure_M);
    fprintf(1,'Макс. дов.   %10.3f %7.3f %7.3f\n', MaxDov_Delta_S_Figure_M, MaxDov_Delta_Centroid_M, MaxDov_Delta_L_Figure_M);
    fprintf(1,'Серед. зн.   %10.3f %7.3f %7.3f\n', Mean_Delta_S_Figure_M, Mean_Delta_Centroid_M, Mean_Delta_L_Figure_M);
    fprintf(1,'СКЗ          %10.3f %7.3f %7.3f\n', Std_Delta_S_Figure_M, Std_Delta_Centroid_M, Std_Delta_L_Figure_M);
    fprintf(1,'З фільт         площа   центр мас   лінійні розміри\n');
    fprintf(1,'Максимальна  %10.3f %7.3f %7.3f\n', Max_Delta_S_Figure_R, Max_Delta_Centroid_R, Max_Delta_L_Figure_R);
    fprintf(1,'Макс. дов.   %10.3f %7.3f %7.3f\n', MaxDov_Delta_S_Figure_R, MaxDov_Delta_Centroid_R, MaxDov_Delta_L_Figure_R);
    fprintf(1,'Серед. зн.   %10.3f %7.3f %7.3f\n', Mean_Delta_S_Figure_R, Mean_Delta_Centroid_R, Mean_Delta_L_Figure_R);
    fprintf(1,'СКЗ          %10.3f %7.3f %7.3f\n', Std_Delta_S_Figure_R, Std_Delta_Centroid_R, Std_Delta_L_Figure_R);

    fprintf(1,'Відносні похибки геометричних параметрів, %%\n');
    fprintf(1,'Без фільт       площа   центр мас   лінійні розміри\n');
    fprintf(1,'Максимальна  %10.3f %7.3f %7.3f\n', Max_D1_S_Figure_M, Max_D1_Centroid_M, Max_D1_L_Figure_M);
    fprintf(1,'Макс. дов.   %10.3f %7.3f %7.3f\n', MaxDov_D1_S_Figure_M, MaxDov_D1_Centroid_M, MaxDov_D1_L_Figure_M);
    fprintf(1,'Серед. зн.   %10.3f %7.3f %7.3f\n', Mean_D1_S_Figure_M, Mean_D1_Centroid_M, Mean_D1_L_Figure_M);
    fprintf(1,'СКЗ          %10.3f %7.3f %7.3f\n', Std_D1_S_Figure_M, Std_D1_Centroid_M, Std_D1_L_Figure_M);
    fprintf(1,'З фільт         площа   центр мас   лінійні розміри\n');
    fprintf(1,'Максимальна  %10.3f %7.3f %7.3f\n', Max_D1_S_Figure_R, Max_D1_Centroid_R, Max_D1_L_Figure_R);
    fprintf(1,'Макс. дов.   %10.3f %7.3f %7.3f\n', MaxDov_D1_S_Figure_R, MaxDov_D1_Centroid_R, MaxDov_D1_L_Figure_R);
    fprintf(1,'Серед. зн.   %10.3f %7.3f %7.3f\n', Mean_D1_S_Figure_R, Mean_D1_Centroid_R, Mean_D1_L_Figure_R);
    fprintf(1,'СКЗ          %10.3f %7.3f %7.3f\n', Std_D1_S_Figure_R, Std_D1_Centroid_R, Std_D1_L_Figure_R);

    fprintf(1,'\nПохибки відтворення контурів\n                        без фільт   з фільт\n');
    fprintf(1,'макс., д.т.        %10.3f %10.3f\n', Max_B_M_XOR, Max_B_R_XOR); 
    fprintf(1,'макс. дов., д.т.   %10.3f %10.3f\n', MaxDov_B_M_XOR, MaxDov_B_R_XOR); 
    fprintf(1,'серед., д.т.       %10.3f %10.3f\n', Mean_B_M_XOR, Mean_B_R_XOR); 
    fprintf(1,'мінім.,  д.т.      %10.3f %10.3f\n', Min_B_M_XOR, Min_B_R_XOR); 
end; % {Delta_Show}