Обнаружение нечеткого импульсного шума в Matlab

Question

Я пытаюсь реализовать метод обнаружения нечеткого импульсного шума в Matlab. Я определяю окно 3 * 3, которое для каждого неграничного пикселя изображения с серой шкалой будет вычислять разные градиенты во всех 8 возможных соседях центрального пикселя, проверять нечеткие правила и определять, является ли этот пиксель шумным или нет. Но он просто проходит первый пиксель и вычисляет его правильно; для второго пикселя я получаю следующую ошибку. Может ли кто-нибудь помочь? Кроме того, я пытаюсь определить функцию для вычисления градиентов, можно ли определить такую ​​функцию для всех направлений? ошибка:

Индекс превышает размеры матрицы.

Ошибка в главном (строка 29) g2 = double (img_temp (r, c + 1) - img_temp (r, c));

и вот мой код:

close all 
clc 

[file, path] = uigetfile('*.*' , 'Open an image'); 
filename = strcat(path, file); 
img = (imread(filename)); 

dim = ndims(img); 

if (dim==3) 
    img = rgb2gray(img); 

end 

figure, imshow(img); 

k = 1; 
[row , col] = size(img); 
for r=2:row-1 
    largeCount = 0; 

    for c=2:col-1 
     img_temp = img(r-1:r+1, c-1:c+1); 

    %% Gradient Calculation in Direction : N 
     g0 = double(img_temp(r-1,c) - img_temp(r,c)); 
     g1 = double(img_temp(r, c-1) - img_temp(r,c)); 
     g2 = double(img_temp(r, c+1) - img_temp(r,c)); 

     g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : NE 
    g0 = double(img_temp(r-1 , c+1) - img_temp(r,c)); 
    g1 = double(img_temp(r-1, c-1) - img_temp(r,c)); 
    g2 = double(img_temp(r+1 , c+1) - img_temp(r,c)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : E 
    g0 = double(img_temp(r,c+1) - img_temp(r,c)); 
    g1 = double(img_temp(r-1,c) - img_temp(r,c)); 
    g2 = double(img_temp(r+1 ,c) - img_temp(r,c)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : SE 
    g0 = double(img_temp(r+1, c+1) - img_temp(r,c)); 
    g1 = double(img_temp(r-1 , c+1) - img_temp(r,c)); 
    g2 = double(img_temp(r+1 , c-1) - img_temp(r,c)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : S 
    g0 = double(img_temp(r+1, c) - img_temp(r,c)); 
    g1 = double (img_temp(r , c+1) - img_temp(r,c)); 
    g2 = double(img_temp(r , c-1) - img_temp(r,c)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : SW 
    g0 = double(img_temp(r+1, c-1) - img_temp(r,c)); 
    g1 = double(img_temp(r+1 , c+1) - img_temp(r,c)); 
    g2 = double(img_temp(r-1, c-1) - img_temp(r,c)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : W 
    g0 = double(img_temp(r,c-1) - img_temp(r,c)); 
    g1 = double(img_temp(r+1, c) - img_temp(r,c)); 
    g2 = double(img_temp(r-1, c) - img_temp(r,c)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : NW 
    g0 = double(img_temp(r-1,c-1) - img_temp(r,c)); 
    g1 = double(img_temp(r+1 , c-1) - img_temp(r,c)); 
    g2 = double(img_temp(r-1 , c+1) - img_temp(r,c)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 
    k = k+1; 
    end 

    %% if largeCount > 4 then the pixel is noisy for sure 
    if largeCount> 4 
     %% Add the pixel value to histogram 
     out(r,c)= 0; 
    else 
     %% Don't change pixel value 

     output(r,c) = (temp(r,c)); 
    end 

end 

figure ; imshow(output); 

EDIT:

Я изменил код, однако, в то время как я бегу следующий код, он получает паузу (функция стека вызовов: ismemeber) и Дон Покажите мне img_out. Я получаю эту ошибку:

69 [sortuAuB, IndSortuAuB] = sort ([uA; uB]);

и вот мой отредактированный код:

close all 
clc 

[file, path] = uigetfile('*.*' , 'Open an image'); 
filename = strcat(path, file); 
img = (imread(filename)); 

dim = ndims(img); 

if (dim==3) 
    img = rgb2gray(img); 

end 

figure, imshow(img); 

k = 1; 
out = readfis('NoiseDetection.fis'); 

[row , col] = size(img); 

img_out = zeros(row , col , 'uint8'); 

for r=2:row-1 
    largeCount = 0; 

    for c=2:col-1 
     img_temp = img(r-1:r+1, c-1:c+1); 

    %% Gradient Calculation in Direction : N 
     g0 = double(img_temp(1,2) - img_temp(2,2)); 
     g1 = double(img_temp(2,1) - img_temp(2,2)); 
     g2 = double(img_temp(2,3) - img_temp(2,2)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : NE 
    g0 = double(img_temp(1,3) - img_temp(2,2)); 
    g1 = double(img_temp(1,1) - img_temp(2,2)); 
    g2 = double(img_temp(3,3) - img_temp(2,2)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : E 
    g0 = double(img_temp(2,3) - img_temp(2,2)); 
    g1 = double(img_temp(1,2) - img_temp(2,2)); 
    g2 = double(img_temp(3,2) - img_temp(2,2)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : SE 
    g0 = double(img_temp(3,3) - img_temp(2,2)); 
    g1 = double(img_temp(1,3) - img_temp(2,2)); 
    g2 = double(img_temp(3,1) - img_temp(2,2)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : S 
    g0 = double(img_temp(3, 2)- img_temp(2,2)); 
    g1 = double(img_temp(2, 3)- img_temp(2,2)); 
    g2 = double(img_temp(2 ,1)- img_temp(2,2)); 

    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : SW 
    g0 = double(img_temp(3,1) - img_temp(2,2)); 
    g1 = double(img_temp(3,3) - img_temp(2,2)); 
    g2 = double(img_temp(1,1) - img_temp(2,2)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : W 
    g0 = double(img_temp(2,1)- img_temp(2,2)); 
    g1 = double(img_temp(3,2)- img_temp(2,2)); 
    g2 = double(img_temp(1,2)- img_temp(2,2)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 

    %% Gradient Calculation in Direction : NW 
    g0 = double(img_temp(1,1)- img_temp(2,2)); 
    g1 = double(img_temp(3,1)- img_temp(2,2)); 
    g2 = double(img_temp(1,3)- img_temp(2,2)); 


    g_weight = max (max (evalfis([g0 g1 0],out) , evalfis([g0 0 g2],out)), evalfis([g0 g1 g2],out)); 
    if g_weight >128 
     largeCount = largeCount+1; 
    end 
    %% if largeCount > 4 then the pixel is noisy for sure 
     if largeCount> 4 

    %% Add the pixel value to histogram 
     img_out(r,c)=0; 

     else 
     %% Don't change pixel value 

     img_out(r,c) = img_temp(2,2); 
     end 
    end 
    k = k+1; 
end 

figure ; imshow(img_out); 

Answer 1

Вы извлекаете пиксельные окрестности в 3 х 3 для каждого центрального пиксела в изображении. Однако, когда вы оцениваете градиенты, вы используете систему координат относительно исходного изображения, а не выделенной окрестности пикселя.

Поэтому все ваши вычисления градиента неправильны. Это видно из-за того, что r и c могут выходить за пределы значений 3, которые являются наибольшим размером строк и столбцов вашего района.

Простейшим решением было бы исправить все операции индексирования в окрестности пикселя 3 x 3, чтобы они были относительно системы координат извлеченного изображения, а не исходной системы координат изображения.

Заменить все r и c с 2 и 2, все r-1 и c-1 с 1 и 1 и, наконец, все r+1 и c+1 с 3 и 3. Что касается точности вашего алгоритма, это должны быть единственные исправления, которые вам нужны, и если я правильно понимаю ваше описание, это должно более или менее дать вам то, что вы хотите. Однако есть несколько подозрительных переменных, к которым обращаются, но они не определены в вашем скрипте. Такими переменными являются out, output и temp. Я беру их по вере и предполагаю, что вы заявили об этом перед этим кодом. Есть более векторизованные подходы к тому, что вы хотите, но я оставлю это как упражнение.

В частности, изменения:

g0 = double(img_temp(r-1,c) - img_temp(r,c)); 
    g1 = double(img_temp(r, c-1) - img_temp(r,c)); 
    g2 = double(img_temp(r, c+1) - img_temp(r,c)); 

To:

g0 = double(img_temp(1,2) - img_temp(2,2)); 
    g1 = double(img_temp(2, 1) - img_temp(2,2)); 
    g2 = double(img_temp(2, 3) - img_temp(2,2)); 

Изменение:

g0 = double(img_temp(r,c+1) - img_temp(r,c)); 
g1 = double(img_temp(r-1,c) - img_temp(r,c)); 
g2 = double(img_temp(r+1 ,c) - img_temp(r,c)); 

To:

g0 = double(img_temp(2,3) - img_temp(2,2)); 
g1 = double(img_temp(1,2) - img_temp(2,2)); 
g2 = double(img_temp(3 ,2) - img_temp(2,2)); 

... и так далее и так далее. Их слишком много для изменения для вас, поэтому я предполагаю, что вы можете сделать это для остальных затронутых блоков кода.