Fuzzy impulse noise detection in Matlab

Question

I am trying to implement fuzzy impulse noise detection method in Matlab. I define a 3*3 window that for each non-border pixel of a gray scale image will calculate different gradients in all 8 possible neighbors of central pixel, check fuzzy rules and detect if that pixel is noisy or not. But it just goes through the first pixel and calculate it correctly; for the second pixel I get following error. Can anyone help please? In addition, I am trying to define a function for calculating gradients , is this possible to define such a function for all directions? error :

Index exceeds matrix dimensions.

Error in Main (line 29) g2 = double(img_temp(r, c+1) - img_temp(r,c));

and here is my code :

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 :

I have changed my code however, while I am running the following code, it gets paused (Function call stack : ismemeber) and don't show me img_out. I get this error :

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

and here is my edited code:

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

You are extracting a 3 x 3 pixel neighbourhood for each central pixel in the image. However, when you are evaluating the gradients, you are using the coordinate system with respect to the original image and not the extracted pixel neighbourhood itself.

Therefore all of your gradient calculations are wrong. This is apparent due to the fact that r and c may go beyond the values of 3 which is the largest dimension of the rows and columns of your neighbourhood.

The simplest fix would be to correct all indexing operations into your 3 x 3 pixel neighbourhood so that they are with respect to the the coordinate system of the extracted image and not the original image coordinate system.

Replace all r and c with 2 and 2, all r-1 and c-1 with 1 and 1 and finally all r+1 and c+1 with 3 and 3. As to the accuracy of your algorithm, this should be the only fixes you need and if I understand your description correctly, this should more or less give you what you want. There are however a few suspect variables that are being accessed but aren't defined in your script. Such variables are out, output and temp. I'm taking these by faith and assuming you declared these before this code. There are more vectorized approaches to doing what you want but I'll leave that as an exercise.

Specifically, change:

   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));

Change:

 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));

... and so on and so forth. There are too many of them to change for you so I'm assuming you can do that for the rest of the affected code blocks.