Lee filter

Let’s realize a Lee filter using MATLAB and C for despeckling of an image. For Lee filter using MATLAB explicitly check : https://www.imageeprocessing.com/2020/07/lee-filter-using-matlab_16.html

Since it’s a patch based processing, the computation cost will be high.

In order to reduce the same, a part of the code is realized in C language for improved performance.

MATLAB code:

%Read an Image

I = imread('coins.png');

%Add multiplicative noise to the image

J = imnoise(I,'speckle',0.01);

figure,imshow(J);

%Apply Lee filter

K = Lee_filter_C(I,J,[5 5]);

figure,imshow(uint8(K));

Create MATLAB function Lee_filter_C :

This function takes the reference image, speckled/noisy image and the window size as input and performs the following steps.

1.     The variance of the reference image is found. Variance can be found either by using MATLAB built-in function or user defined function. Here in this case, a user defined function is used to find the variance.

2.     Based on the size of the kernel, the noisy image is padded with zeros on all sides.

3.     The center index of the kernel is found

4.     The noisy image is processed patch by patch.

5.     The code written in C computelee.c to despeckle the image is used.

MATLAB code:

function Y = Lee_filter_C(R,E,sz)

%R is the Reference Image

%E is the Error or Noisy Image

%K is the Kernel or Window

%Y is the Output Image

% Y = mean(K)+W*(C-mean(K);

% W = variance(K)/(variance(K)+variance(R))

%Define the type

R = double(R);

E = double(E);

%Preallocate the Output Matrix

Y = zeros(size(R));

mn = round((sz-1)/2);

Tot = sz(1,1)*sz(1,2);

EImg = padarray(E,mn);

%Variance of the reference Image

Rvar = myvar(R);

%Rvar = var(R(:));

Indx = floor(median(1:Tot));

for i = 1:size(R,1)

    for j = 1:size(R,2)

        K = EImg(i:i+sz(1,1)-1,j:j+sz(1,2)-1);

        Y(i,j) = computelee(Rvar,K(Indx),K(:)');

       

    end

end

end

NOTE: save the above function as Lee_filter_C.m 

Function to find the variance:

MATLAB code:

function var_v = myvar(I)

I = I(:);

var_v = sum((I-mean(I)).^2)/numel(I);

end

NOTE: Save the above function as myvar.m

C program: computelee.c

Basics on MEX files in MATLAB:

If you are familiar with C programming then we just need to understand the gateway from MATLAB to C programming. The rest of the procedures will be same as we code in c.

The main gateway function for writing the C program is Mex function and it normally takes 4 parameters.

Nlhs – To find number of left hand side parameters

Plhs – Contains all the output parameters in an array

Nrhs – To find the number of right hand side parameters

Prhs – Contains the input parameters in an array

C code:

#include "mex.h"

/* Find variance and mean of the pixels in the window */

void arrayProduct(double v, double In, double *y, double *z, mwSize n)

{

    double VarW,MeanW,W;

    mwSize i;

   

    MeanW=0;

    for (i=0; i

        MeanW=MeanW+y[i];

    }

    MeanW=MeanW/n;

    VarW=0;

    for(i=0;i

    {

        VarW=VarW+((y[i]-MeanW)*(y[i]-MeanW));

    }

    VarW=VarW/n;

   

    W=VarW/(VarW+v);

    z[0]=MeanW+W*(In-MeanW);

}

void mexFunction( int nlhs, mxArray *plhs[],

                  int nrhs, const mxArray *prhs[])

{

    double Rvariance,CIndx,*inMatrix,*out;              

    size_t ncols;               

   

   

    Rvariance = mxGetScalar(prhs[0]);

    CIndx  = mxGetScalar(prhs[1]);

    inMatrix = mxGetPr(prhs[2]);

    ncols = mxGetN(prhs[2]);

   

    plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL) 

    out = mxGetPr(plhs[0]);

    /* Call the function arrayProduct */

    arrayProduct(Rvariance,CIndx,inMatrix,out,(mwSize)ncols);

}

Let’s compare the calling function in MATLAB and the C code for better understanding.

MATLAB code:

computelee(Rvar,K(Indx),K(:)');

Three parameters are passed to the function computeelee.c

1. Variance of the reference image

2. Centre pixel from the kernel

3. All the pixels from the kernel based on the window size

C code:

Rvariance = mxGetScalar(prhs[0]);

CIndx  = mxGetScalar(prhs[1]);

inMatrix = mxGetPr(prhs[2]);

Rvariance contains the variance of the reference image

CIndx contains the centre pixel of the corresponding kernel/patch

inMatrix contains all the pixels of the corresponding kernel/patch

ncols = mxGetN(prhs[2]);

ncols will obtain total number of elements in kernel.

plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); 

out = mxGetPr(plhs[0]);

The final result will be stored in the variable ‘outMatrix’.

The function ‘arrayProduct’ written in C language is called to perform the computation.

Steps to be performed:

Let’s see how Lee filter works on homogeneous area and edge based area.

Homogeneous Patch

Edges and lines

The result shows that the filter works well on homogeneous area rather than on edges and lines.

Points to remember when using Mex file:

1.     If you are using the C code for the first time, its better to find the appropriate compiler and configure to use it.

2.     Compile the C code if the compiler is already present and make sure the compilation is successful.

3.     Syntax to setup the compiler: mex  - setup
The above syntax will identify the available compilers in your local system and you can configure it manually.

4.     For compiling the above c code:  mex computelee.c

For successful compilation, the sample output will look like below:

Building with 'lcc-win32'.

MEX completed successfully.

5.     The result of successful compilation will generate computelee.mexw32 for 32 bit system and for 64 bit system, you can find computelee.mexw64 file in your local directory.

To sum up, the two functions Lee_filter_C.m and myvar.m should be placed in the same directory. Computelee.c should also be placed in the same directory and the successful compilation of the same is mandatory.

I hope the readers of my blog will be familiar with working on Mex files. If you need any tutorial or post on this topic, kindly post or mail your suggestions.

One Dimensional Low pass , High Pass and band pass filtering

Consider a one dimensional signal in time domain.

For instance, generate cosine waves of different amplitudes and different frequencies and combine them to form a complicated signal.

Consider four cosines waves with frequencies 23, 18, 10 and 5 and amplitudes 1.2, 0.8, 1.1 and 2 respectively and combine them to generate a signal.

Eg.

Cosine wave with frequency 23 and amplitude 1.2.

 

Figure.1

Cosine wave with frequency 18 and amplitude 0.8

 

Figure.2

Cosine wave with frequency 10 and amplitude 1.1

 

Figure.3

Cosine wave with frequency 5 and amplitude 2.

Figure.4

MATLAB code:

%To generate a signal

%Input function A cos(2pift)

freq =[23 18 10 5];

Amp = [1.2 0.8 1.1 2];

fs = 100;

ts = 1/fs;

t = 0:ts:5;

xt1 = 0;

for Ind = 1:numel(freq)

Freq = freq(Ind);

A = Amp(Ind);

xt1 = xt1+A*cos(2*pi*Freq*t);

end

figure,plot(t,xt1);

 

Figure.5

Let’s make a little understanding of this signal in frequency domain.

i.e in Fourier domain.

Fourier transform is one of the various mathematical transformations known which is used to transform signals from time domain to frequency domain.
The main advantage of this transformation is it makes life easier for many problems when we deal a signal in frequency domain rather than time domain.

After fourier transform, the signal will look like the below figure.

MATLAB code:

f1 = linspace(-fs/2,fs/2,numel(xt1));

xf = fftshift(fft(xt1));

str = 'Frequency domain';

figure,plot(f1,abs(xf)),title(str);

 

Figure.6

EXPLANATION:

In the figure.6, there are 4 different colors to denote four different signals.

Due to symmetric property, we see two peaks of the same signal.

fftshift is useful for visualizing the Fourier transform with the zero-frequency component in the middle of the spectrum.

The yellow color peaks denote the frequency with 5

The blue color peaks denote the frequency with 10

The green color peaks denote the frequency with 18

And the red color denotes the frequency with 23.

From this, we can understand that the low frequency components are close to the centre of the plot and the high frequency components are away from the centre.

Now, lets try to separate each of the original 4 signals from the combined ones.

We will keep only one signal at a time and remove/filter the rest.

MATLAB code:

       Kernel_BP = abs(f1) < 8;

   BPF = Kernel_BP .* xf;

  

   figure,plot(f1,abs(BPF));

Figure.7

Figure.8

A kernel is used to obtain the essential frequency component alone.

This can be done by making the values one for the frequency range less than 8 and rest zero.

We have obtained the lowest frequency component and the kernel we designed is for low pass filtering. The value which we use to find the limit is the cut off frequency. Here 8 is the cut of frequency.

Now , let’s go back to time domain.

MATLAB code:

BP_t = ifft(ifftshift(BPF));

plot(t,real(BP_t));

Compare the figure.4 and figure.9, both have same frequency right!!

Figure.9

Now let’s try to design the kernel for filtering high frequency component.

High Pass filter:MATLAB code :

    

   Kernel_BP = abs(f1) > 20;

   BPF = Kernel_BP .* xf;

  

   figure,plot(f1,abs(BPF));

 

Figure.10

   

Figure.11

  BP_t = ifft(ifftshift(BPF));

  plot(t,real(BP_t));

Figure.12

We have obtained the same signal as in figure.1

We will try to filter the signal of frequency 18.

This is band pass filter.

Band Pass Filter:
MATLAB code:

   Kernel_BP = abs(f1) > 15 & abs(f1) < 20;

   BPF = Kernel_BP .* xf;

  

   figure,plot(f1,abs(BPF));

         

Figure.13

   BP_t = ifft(ifftshift(BPF));

   plot(t,real(BP_t));

 

Figure.14

Figure.15

This same type of filtering can be applied on noisy signals to remove undesired signals.

MATLAB code:

Noise = 2 * randn(size(t));

xt1 = xt1 + Noise;

Using this signal as input, we can apply low, high and band pass to remove noise and obtain the original signal(to some extent).

Wait for my next post on two dimensional low and high pass filtering. i.e on images in space domain.

Gaussian Filter without using the MATLAB built_in function

Gaussian Filter

Gaussian Filter is used to blur the image. It is used to reduce the noise and the image details.

          The Gaussian kernel's center part  ( Here 0.4421 )  has the highest value and  intensity of other pixels  decrease as the distance from the center part increases.

               On convolution of the local region and the Gaussian kernel gives the highest intensity value to the center part of the local region(38.4624) and the remaining pixels have less intensity as the distance from the center increases.

               Sum up the result and store it in the current pixel location(Intensity = 94.9269) of the image.

MATLAB CODE:

%Gaussian filter using MATLAB built_in function

%Read an Image

Img = imread('coins.png');

A = imnoise(Img,'Gaussian',0.04,0.003);

%Image with noise

figure,imshow(A);

H = fspecial('Gaussian',[9 9],1.76);

GaussF = imfilter(A,H);

figure,imshow(GaussF);

MATLAB CODE for Gaussian blur WITHOUT built_in function:

%Read an Image

Img = imread('coins.png');

A = imnoise(Img,'Gaussian',0.04,0.003);

%Image with noise

figure,imshow(A);

 

Image with Noise

I = double(A);

%Design the Gaussian Kernel

%Standard Deviation

sigma = 1.76;

%Window size

sz = 4;

[x,y]=meshgrid(-sz:sz,-sz:sz);

M = size(x,1)-1;

N = size(y,1)-1;

Exp_comp = -(x.^2+y.^2)/(2*sigma*sigma);

Kernel= exp(Exp_comp)/(2*pi*sigma*sigma);

 

Gaussian Kernel 9x9 size with Standard Deviation =1.76

%Initialize

Output=zeros(size(I));

%Pad the vector with zeros

I = padarray(I,[sz sz]);

%Convolution

for i = 1:size(I,1)-M

    for j =1:size(I,2)-N

        Temp = I(i:i+M,j:j+M).*Kernel;

        Output(i,j)=sum(Temp(:));

    end

end

%Image without Noise after Gaussian blur

Output = uint8(Output);

figure,imshow(Output);

After Filtering