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 !!

Multi-variate Model Fitting using Taylor Series Method

In the previous example, we have learnt how to calculate a singlevariate model using a taylor series method when no knowledge about form of relation is known between input and output. This example will extend the same using iterative application of singlevariate model fitting.

We will consider a bivariate model and same can be easily generalized for higher dimensions.

Problem Statement:
We consider a system where there are two inputs/parameters (x and y) which we can change and observe the output (z). Our aim is to find a relation between this output and inputs but again we don't know what relation form(whether it is exponential, logarithmic, linear or trigonometric) exists.

Taylor Series Method
see this article for Taylor Series Method http://matlabbyexamples.blogspot.com/2011/12/singlevariate-model-fitting-for.html.

Iterative Taylor Series Method
We will first make one parameter (y) constant (let us say at y0) and change the other one (x) only. Then we try to fit a taylor series model in to it and we get coefficient b_n|at y=y0 such that  \[ f(x) = \sum_n{b_n x^n} \].
Now we will change the parameter y to y1 and do the same fit again to get a different set of coefficients  b_n|at y=y1. Do it for all value of y and each time you get a 1xN vector b containing b_n's. Now stack them such that  we get a matrix B each row represent b for one value of y. Now observe that ith column (let us call it B_i) represents b_i's for all values of y.

Now it is clear that these columns are only dependent on value of y. So again we can fit a taylor series to each of these columns such that  \[ b_n = \sum_m{c_m y^m} \].  This can be done by solving \[B_n = Y*C\]

So we have the final model as
\[ f(x) = \sum_n{\sum_m{c_{m,n} y^m} x^n} \].

MATLAB Implementation:

function C=taylorseriesfit2D(X,Y,Z);
%assumes X and Y are meshgrid format matrix
for i=1:length(y)
z=Z(i,:)';
B(i,:)=taylorseriesfit1D(X(1,:),z);
end

for i=1:size(B,2)
C(:,i)=taylorseriesfit1D(Y(:,1),B(:,i));
end

where taylorseriesfit1D is given as

function B= taylorseriesfit1D(x,z)
Xmat=[ones(size(x)) x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7];
B=Xmat\z;

Simulation Data generation:

x=[0:0.1:1]';
y=[0:.1:pi/2]';
[X Y ]=meshgrid(x,y);
Z=exp(X).*exp(Y);
C=taylorseriesfit2D(X,Y,Z)

C will come like for N=7.

1.0000 1.0000 0.5000 0.1667 0.0416 0.0085 0.0012 0.0003
1.0000 1.0000 0.5000 0.1667 0.0416 0.0085 0.0012 0.0003
0.5000 0.5000 0.2500 0.0833 0.0208 0.0043 0.0006 0.0002
0.1669 0.1669 0.0835 0.0278 0.0069 0.0014 0.0002 0.0001
0.0410 0.0410 0.0205 0.0068 0.0017 0.0004 0.0000 0.0000
0.0092 0.0092 0.0046 0.0015 0.0004 0.0001 0.0000 0.0000
0.0008 0.0008 0.0004 0.0001 0.0000 0.0000 0.0000 0.0000
0.0004 0.0004 0.0002 0.0001 0.0000 0.0000 0.0000 0.0000

Final Guesses 
Again we can guess the model by observing these coefficients and plotting them with respect to m,n (Fig1). First thing we can observer that it is symmetric, therefor exihibit the same relations for both the inputs to ouput.

Fig1: Plot C vrs m,n

Again, we cannot guess, even then these matrixes are sufficient for all numerical purposes.

Singlevariate Model fitting for a Simulation Data without any prior i/o relation knowledge using Taylor Series Method

Singlevariate Model fitting for a Simulation Data without any prior i/o relation knowledge using Taylor Series Method

In this example we will talk about how we can fit a model/function in a single variate simulation data. In many practical cases, we get some observation data after changing an input variable and our aim is to find the relation between input and output. If we know the relation form atleast (like y=mx+c or y=d*sin(x)^2), then we can find the values of these coeefficients y,c or d by linear/non linear regression. But what if we dont know the relation form, we dont have any information about whether they follow exponential form or logarthim relation. Here is a trick which we are going to discuss in this example.

Taylor Series fit
Let i/o relation be represented by a function z=f(x). We assume here that f(x) is continuous and differentiable, which is generally the case with practical/real time systems. Now each such function can be represented using a taylor's series
\[ f(x) = \sum_n{b_n x^n} \] we will limit the summation to n=N
Since we know the data x and z=f(x), we can calculate the coefficients bi .
Let \[B=\left[\begin{array}{ccccc} b_0 & b_1 & b_2  & ... & b_N \end{array} \right]^T \] and \[ Xmat =\left[ \begin{array}{ccccc} 1 & x^1 & x^2 & ... & x^N \end{array} \right] \]
so \[ z=f(x)= Xmat B \]

This is also true when x is a column vector mx1 and contains m data points and z contains values of f at those m point.
So B can be calculated as \[B=Xmat^{-1} Z\]

function: taylorseriesfit()
Writing the function in matlab is very easy. we just need to follow the steps as explained above.
Since our function need to matrixs x and z, these will be input to the function and B will be the output of the function. We will take N=4 here. Increasing the value may result in accuracy of your fit.

function B= taylorseriesfit(x,z)
Xmat=[ones(size(x)) x x.^2 x.^3 x^4];
B=Xmat\z; 

Here x and z are column vectors representing data points.

Generation of Simulation Data 
We will generate a sample Simulation data for our example. In real example, you will get this data from some experiment after changing an input parameter(x) and observing the output(z).

x=[0:0.1:1]';

z=exp(x);

Calling the defined function

B=taylorseriesfit(x,z);
disp(B);

You will get B as

B =
1.0000
0.9988
0.5099
0.1400
0.0696

Now your model is
\[ z= x*B = B(1)+B(2)x+B(3)x^2+B(4)x^3+B(5)x^4 \] \[ z= 1+0.99x+0.50 x^2+0.14x^3+.07x^4 \]

Error
we can calculate the modelling error (MMSE) by comparing the actual output and the simulation output.

z_o = Xmat*B;
e=sqrt(sum((z-z_o).^2));
disp(e);

you will get something like

e= 6.9400e-005

which is pretty good.

Final Observation
Now if we observe the coefficients and plot them with respect with n, we may guess the actual function. Fig1 represents such a graph for N=7 case.

Fig1 : b_i plotted versus i

like for this example, it looks like 1/n!, so this is an exponential relations. Well, if we can guess, it is good. But if we cannot, still we have a relation which is sufficient for all numerical purposes.