Median Filtering using MATLAB

In this blog we will learn how we can do a time domain filtering for an image. We will consider median filtering for this example. We will see two methods - first one is the iterative methhod which is time consumingg. Second methods uses matlab's vectorization and performs well.

Example Image
For this blog, we will take a very short image to understand how filtering is being done. Our sample image is shown in Fig1 (7x7).

Fig1 : Image A (magnified)

A =
8 5 8 0 6 8 7
9 9 1 8 1 6 7
1 9 4 9 7 3 1
9 1 9 6 0 9 4
6 9 7 7 2 0 4
0 9 9 7 0 4 6
2 4 6 3 0 3 7

You can create this by the following command

A=[8 5 8 0 6 8 7;9 9 1 8 1 6 7;1 9 4 9 7 3 1;9 1 9 6 0... 9 4; 6 9 7 7 2 0 4;0 9 9 7 0 4 6;2 4 6 3 0 3 7];

In actual cases you will read a image by using imread

A=imread('C:\users\Myname\codes\image1.jpg');
A=rgb2gray(A);

[l m]=size(A);

Kernel
We would be filtering with 3x3 kernel. It means that for each pixel of the output image, we would take a 3x3 neighborhood around the respective pixel in the original image, and assign the median color value in this neighborhood to the pixel of output image.

Output
Output would be a image of 7x7 too. Let us preallocate it with all zeros.

B= zeros(l,m);

Iterative method:
For each i and j, we first extract a 3x3 matrix around this point as shown in figure.

Fig2: Image Kernel Elaboration

K=A(i-1:i+1,j-1:j+1);

This will select a 3x3 matrix containing elements from i-1,i,i+1 row and j-1,j,j+1 column and assign it to K.

Now calculate median of this K and assign it to B(i,j);

B(i,j)=median(K);

We run a loop on i,j = 2 to 6. We cannot start with 1 because 3x3 kernel should fit inside the image to avoid any matrix index out of bound error. In the same way we cannot have i or j = 7. Note that we can take only those i where i-1, i,i+1 all three should come in the range 1 to 7. So with the j.

for i=2:6

  for j=2:6

    K= A(i-1:i+1,j-1:j+1);

    B(i,j)=median(K);

  end

end

 So for example for i=3,j=4, we take the kernel consisting of 2-3-4 row and 4-5-6 column as shown in figure 2. Elements are 1,4,6,7,7,9,9,9,9 and the median is 7 , so the output image's B(3,4)=7.

 This is done for all the indexes.

Padding:
We will first pad the original image to take care of the corner cases. To pad we will create a new matrix with l+2, m+2 dimension and assign the middle region from i=2 to l+1, j=2 to m+1 from the original image.

Y=zeros(l+2,m+2);
Y(2:l+1,2:m+1)=A;

Now remember rth,sth element of Y is equal to ith,jth element of A where i=r-1,j=s-1 for r,s>2.
Now you can go again i=2:l+1, 2:m+1, but remember to assign it as following

B(i-1,j-1)=median(K);

Vectorized Method:
The iterative method is little slow. MATLAB is optimized for vectorized operation. Let us try the same with vectorized way but it is little tricky. As you know, to calculate median, you need to pass MATLAB a matrix.

First we will create 9 shifted version of A (let us say A1, A2,... A9) so that for any i,j,A1(i,j) to A9(i,j) represent the 9 values in 3x3 neighborhood of A(i,j).
since we are considering i=2 to 6, these 9 matrixes will be 5x5 size;

A1=A(1:5,1:5);
A2=A(1:5,2:6);
A3=A(1:5,3:7);


A4=A(2:6,1:5);
A5=A(2:6,2:6);
A6=A(2:6,3:7);


A7=A(3:7,1:5);
A8=A(3:7,2:6);
A9=A(3:7,3:7);

Now we need to find median of 9 values A1(i,j) to A9(i,j) for all i,j. We first combine them into a 3D matrix and calculate the median in its first dimension.

C(1,:,:)=A1;
C(2,:,:)=A2;
C(3,:,:)=A3;
C(4,:,:)=A4;
C(5,:,:)=A5;
C(6,:,:)=A6;
C(7,:,:)=A7;
C(8,:,:)=A8;
C(9,:,:)=A9;
D=median(C,1);
B=squeeze(D);

And Done! Here is you answer B. This is pretty fast specially for large matrixes. There are only 9 operations (equal to kernel size) are happening, while in iterative one, m*n operations are happening, so if m=n=1000, there are 1000000 operations !!