Lets Learn together... Happy Reading

" Two roads diverged in a wood, and I,
I took the one less traveled by,
And that has made all the difference "-Robert Frost

Upsampling in Frequency Domain




1.1  Upsampling using MATLAB built-in function
1.2  Upsampling in 1D
1.3  Upsampling in 2D or image matrix
2.1  Upsampling a 1D signal
2.2  Upsampling a image matrix

In Frequency domain, upsampling means nothing but the padding of zeros at the end of high frequency components on both sides of the signal.

STEPS TO PERFORM:

1.      Read an image
2.      Obtain the ratio to upsample
3.      Perform Fast Fourier Transform
4.      Shift the Low frequency components to the centre and High frequency components outside.
5.      Add zeros on both the sides of the image
6.      Shift the High frequency components to the centre and Low frequency components to the exterior (Inverse of fftshift)
7.      Perform Inverse Fast Fourier Transform
8.      Display the Upsampled Image

MATLAB CODE:

%UPSAMPLING IN FREQUENCY DOMAIN
%1D UPSAMPLING

FS = 100;
t  = 0:(1/FS):1;

A = 10*sin(2*pi*5*t);
figure,plot(t,A);
 



%FOURIER DOMAIN

FT = fft(A);
fq =linspace(-1/FS,1/FS,101);
figure,plot(fq,abs(FT));title('Before FFTSHIFT');


FT_s = fftshift(FT);
figure,plot(fq,abs(FT_s));title('After FFTSHIFT');


pad_zero = padarray(FT_s,[0 50]);

fq =linspace(-1/FS,1/FS,201);
figure,plot(fq,abs(pad_zero));title('After PADDING WITH ZEROS');


%INVERSE FOURIER TRANSFORM
IFT = ifft(ifftshift(pad_zero));

%UPSAMPLED SIGNAL
t1 = linspace(0,1,201);
figure,plot(t1,(IFT*2),'r',t,A,'g');




EXPLANATION:

Amplitude of the input original signal is 10 and the frequency is 5.
Similarly, the amplitude of the upsampled signal is 10 and the frequency is 5. The number of samples used to plot the signal is increased in the later case.



IMAGE UPSAMPLING IN FOURIER DOMAIN

MATLAB CODE:

%READ AN INPUT IMAGE
A = imread('cameraman.tif');

%RATIO
RatioM = 3;
RatioN = 2;

%UPSAMPLING OVER EACH ROW

mnrow = round(size(A,2)*(RatioM-1)/2);
% 1D FFT ON EACH ROW
row_fft = fft(A,[],2);

%PAD WITH ZEROS ON BOTH SIDES OF EACH ROW
pad_row = padarray(fftshift(row_fft,2),[0 mnrow]);

 
Logarthmic scale
%UPSAMPLING OVER EACH COLUMN
mncol = round(size(A,1)*(RatioN-1)/2);

% 1D FFT ON EACH COLUMN
col_fft = fft(pad_row,[],1);

%PAD WITH ZEROS ON BOTH SIDES OF EACH COLUMN
pad_col = padarray(fftshift(col_fft,1),[mncol 0]);
 
Logarthmic scale
%PERFORM 1D IFFT ON EACH COLUMN
ifft_col = ifft(ifftshift(pad_col,1),[],1);

%PERFORM 1D IFFT ON EACH ROW
ifft_col_row = ifft(ifftshift(ifft_col,2),[],2);

%DISPLAY THE IMAGE AFTER UPSAMPLING
res = abs(ifft_col_row);
res = uint8(res*(numel(res)/numel(A)));


figure,imagesc(res);




SIMPLE VERSION :

A = imread('cameraman.tif');

%RATIO
RatioM = 3;
RatioN = 2;



mnrow = round(size(A,2)*(RatioM-1)/2);
mncol = round(size(A,1)*(RatioN-1)/2);

%FFT ON 2D MATRIX
FT = fftshift(fft2(A));


%PADDING WITH ZEROS
pad_rc = padarray(FT,[mncol mnrow]);


%INVERSE FOURIER TRANSFORM
IFT = ifft2(ifftshift(pad_rc));

Img = uint8(abs(IFT)*(numel(IFT)/numel(A)));
figure,imagesc(Img);




like button Like "IMAGE PROCESSING" page

Convolution in MATLAB

                 Let us try to understand convolution by performing spatial averaging on a matrix without using MATLAB built in function ‘conv2()’.

%CONVOLUTION IN MATLAB

%INPUT  MATRIX
A = zeros(5);
A(:) = 1:25;

%KERNEL
avg3 = ones(3)/9;

%CONVOLUTION
Result = conv2(A,avg3,'same');
display(Result);

Steps to be performed:





NOTE :
To define a kernel for spatial averaging, fill the kernel with ones and divide it by the number of elements in it.
For instance, consider kernel of size 4x4 , fill the  matrix with ones and divide it by 16. i.e the total number of elements in the matrix.

MATLAB CODE:
%INPUT  MATRIX
A = zeros(5);
A(:) = 1:25;

%KERNEL
avg3 = ones(3)/9;

%PAD THE MATRIX WITH ZEROS
B = padarray(A,[1 1]);
% PRE-ALLOCATE THE MATRIX
Output = zeros([size(A,1) size(A,2)]);

%PERFORM COONVOLUTION
for i = 1:size(B,1)-2
    for j = 1:size(B,2)-2
        Temp = B(i:i+2,j:j+2).*avg3;
        Output(i,j) = sum(Temp(:));
    end
 end

display(Output);

like button Like "IMAGE PROCESSING" page

Multilook Technique for speckle reduction




Consider a stack of images  affected by speckle noise of a same scene. 

In order to reduce the speckle, the availability of the stack of images can be utilized to obtain speckle reduced image.
In order to understand the multilook technique, let’s first generate noisy images and then apply speckle reduction.




Generate stack of noisy images:
1.      Read a noise free image
2.      Generate exponential noise with mean 1 and multiply it with the noise free image

3.      Repeat the process of generating exponential noise and multiplying  it with the noise free image to create stack of ‘n’ images

MATLAB code:
%Multilook - Speckle reduction

%Read a noisy free Image
I = imread('cameraman.tif');
I = double(I);
figure,imagesc(I);colormap(gray);title('Original Image');

[m, n] = size(I);
numofimg = 9;
Stack = zeros([m n numofimg]);

for i = 1:numofimg
    %Generate exponential noise
    noise=random('exp',1,m,n);
   
    %Multiply Noise free Image with noise
    Stack(:,:,i)=I.*noise;
   
end

Explanation:

Random noise with exponential distribution is multiplied with the noise free image to generate image affected by speckle.





Fig: Probability Density Function (PDF) of the noise generated for ‘n’ images that follows exponential distribution.

Multilook(ML) Technique:
Multilook image is the average of the stack of images. This technique is used widely in the field of Synthetic Aperture Radar(SAR) to reduce speckle on the same scene but different acquisition periods.
The stack of data is created for the same scene but with different time of acquisition and ML is applied to reduce the speckle considerably since the noise on the single image will be quite high.

MATLAB code:
ML_image = mean(Stack,3);
figure,imagesc(Stack(:,:,1));colormap(gray);title('Single Noisy Image');
figure,imagesc(ML_image);colormap(gray);title('Multilook Image');



Explanation:


For instance, in the above shown image, consider the red dot as the pixel value of the first pixel position of each image in the stack. Add all the pixel values represented in red dot and divide by the number of images in the stack. Similarly, perform the average for the second pixel position (green dot) and so on and so forth for the whole image.

Matlab code ‘mean(Stack,3) ‘ finds the mean of the image in 3rd dimension and the final result is the Multilooked Image with reduced speckle.


Moving Average Filter:
In order to reduce speckle further, a moving average filter can be used.
Moving average filter of window size 3x3 and 5x5 is applied on the multilook image and window of size 3x3 on a single speckled image in order to appreciate the multilook technique performance.

MATLAB code:

box = ones(3)/9;
ML_avg3 = conv2(ML_image,box,'same');
figure,imagesc(ML_avg3);colormap(gray);title('Moving Average - Window 3 x 3');
box = ones(5)/25;
ML_avg5 = conv2(ML_image,box,'same');
figure,imagesc(ML_avg5);colormap(gray);title('Moving Average - Window 5 x 5');
box = ones(3)/9;
avg_flt = conv2(Stack(:,:,1),box,'same');
figure,imagesc(avg_flt);colormap(gray);title('Moving Average - Single Image');







like button Like "IMAGE PROCESSING" page

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.








like button Like "IMAGE PROCESSING" page
Previous Post Next Post Home
Google ping Hypersmash.com