BASIC IMAGE PROCESSING TUTORIAL

BASIC IMAGE PROCESSING TUTORIAL

This tutorial is intended for a basic understanding of some topics on image processing. They are developed using the  MATLAB package and it only contains the commands to be issued. Results should be "experienced" in MATLAB.

Tutorials on standard use of MATLAB can be found at here while the specific use for image processing is here. 
 

Fourier Transform

MATLAB tutorial on Fourier Transform is quite complete. Anyway, in this tutorial the main commands are summarized and some extra experiments done.

To compute the Fourier Transform of an image do the following

     % Prepare image 
     f = zeros(30,30); 
     f(5:24,13:17) = 1; 
     imshow(f,'notruesize')

     % Compute Fourier Transform 
     F = fft2(f,256,256); 
     F = fftshift(F); % Center FFT

     % Measure the minimum and maximum value of the transform amplitude 
     min(min(abs(F))) %   0 
     max(max(abs(F))) % 100 
     imshow(abs(F),[0,100]); colormap(jet); colorbar 
     imshow(log(1+abs(F)),[0,3]); colormap(jet); colorbar 
            * What is the main difference between representing the amplitude and its logarithm?

     % Look at the phases 
     imshow(angle(F),[-pi,pi]); colormap(jet); colorbar

            * Try with other images 
     f = imread('saturn.tif'); 
     f = ind2gray(f,gray(256));

            * Or more patterned from Yamitani (at Caltech Univ.) homepage. 
            * Checkered images 
     function I=checkered(N,s) 
     % This function returns a checkered image of size N 
     % and with a check size of s

     I=zeros(N,N); 
     for i=1:N, 
         for j=1:N, 
             if (mod(floor(i/s)+floor(j/s),2)==0) I(i,j)=1; end; 
         end; 
     end;

     f=checkered(256,32); imshow(f); 
 

Sampling rate (Downsampling)

The sampling rate determines the pixel width. Once you have sampled an image you may downsample it by a factor of 'm'. See  this page for a theoretical explanation of sampling and  this one for a comprenhension of aliasing.

Downsampling can be done in two fashions: by interpolation and by sampling theory.

By interpolation (explore different possibilities of interpolation schemes used by imresize)

     f=imread('blood1.tif'); 
     f=ind2gray(f,gray(256)); 
     f=f(1:256,1:256); 
     f1=imresize(f,0.2);      % Downsample to a 1/5th of the size 
     f2=imresize(f1,5);       % Go back to original size 
     figure; imshow(f); 
     figure; imshow(f1); 
     figure; imshow(f2);

Via sampling theory

     f=imread('blood1.tif'); 
     f=ind2gray(f,gray(256)); 
     f=f(1:256,1:256);

     function Idown=downsampling(I,m) 
     % Downsample the square image I by a factor of m

     [N,M]=size(I);

     % Apply ideal filter 
     w=1/m; 
     F=fftshift(fft2(I)); 
     for i=1:N 
         for j=1:N 
             r2=(i-round(N/2))^2+(j-round(N/2))^2; 
             if (r2>round((N/2*w)^2)) F(i,j)=0; end; 
         end; 
     end; 
     Idown=real(ifft2(fftshift(F)));

     % Now downsample 
     Idown=imresize(Idown,[N/m,N/m],'nearest');

     % ==================================================================== 
     function Iup=upsampling(I,m) 
     % Upsample the square image I by a factor of m

     [N,M]=size(I); 
     Iup=zeros(m*N,m*N);

     % Expand input image 
     for i=1:N 
         for j=1:N 
            Iup(m*(i-1)+1,m*(j-1)+1)=I(i,j); 
         end; 
     end;

     % Ideal filter 
     [N,M]=size(Iup); 
     w=1/m; 
     F=fftshift(fft2(Iup)); 
     for i=1:N 
         for j=1:N 
             r2=(i-round(N/2))^2+(j-round(N/2))^2; 
             if (r2>round((N/2*w)^2)) F(i,j)=0; end; 
         end; 
     end; 
     Iup=(m*m)*abs(ifft2(fftshift(F)));

     f1=downsampling(f,2); 
     f2=upsampling(f1,2); 
     figure; imshow(f); 
     figure; imshow(f1); 
     figure; imshow(f2); 
 

Image Quantization

Image quality strongly depends on the number of bits used for coding grey levels. This is called image quantization. With the following example these concepts should become clear.

     f=load('saturn.tif'); 
     imshow(f); colormap(gray); 
     imshow(f); colormap(gray(32)); 
     imshow(f); colormap(gray(4)); 
     imshow(f); colormap(jet(16)); % False color 
 

Fourier space filtering

We have already made use of sharp frequency defined filters when upsampling/downsampling. The basic idea is to multiply the Fourier transform of the image by the Fourier transform of the filter. For instance, a hard lowpass filter at frequency w is implemented as follows

     function Ifilter=ideal_lowpass(I,w) 
     % ideally filter image I up to a frequency w

     [N,M]=size(I); 
     F=fftshift(fft2(I)); 
     for i=1:N 
         for j=1:N 
             r2=(i-round(N/2))^2+(j-round(N/2))^2; 
             if (r2>round((N*w)^2)) F(i,j)=0; end; 
         end; 
     end; 
     Ifilter=real(ifft2(fftshift(F)));

     f=imread('blood1.tif'); 
     f=ind2gray(f,gray(256)); 
     f=f(1:256,1:256); 
     figure; imshow(ideal_lowpass(f,0.25)); 
     figure; imshow(ideal_lowpass(f,0.10));

Real space filtering

Real space filtration is done via  convolution and in MATLAB the command is  filter2. Many filters are given by a special filter design (see next section) or via  fspecial.

Next example shows you how to perform a low pass filter via convolution with a small kernel

     f=imread('blood1.tif'); 
     f=ind2gray(f,gray(256)); 
     H=0.25*[0.25 0.5 0.25; 0.5 1 0.5; 0.25 0.5 0.25]; 
     f1=filter2(H,f); 
     figure; imshow(f); 
     figure; imshow(f1);

The frequency respose of the filter can be shown with

     freqz2(H); 
 

Filter design

Matlab tutorial on image filtering is quite representative of how to design specific filters. 
 

Histogram operations

Histogram is a measure of the distribution of density within pixels. Image histogram is computed as

     imhist(f);

while imadjust and histeq changes the histogram shape or limits

     f1=imadjust(f,[0.4 0.6],[0 1]); 
     f2=histeq(f); 
     figure; imhist(f); 
     figure; imhist(f1); 
     figure; imhist(f2); 
     figure; imshow(f); 
     figure; imshow(f1); 
     figure; imshow(f2); 
 

Image enhancement

Histogram equalization and adjustment together with linear filtering are common image enhancement operations. However, there exist other operations such us

• Thresholding

  
  f=imread('blood1.tif'); 
  f=ind2gray(f,gray(256)); 
  imhist(f); 
  f2=max(f1,0.4); % remove all values below 0.4 
  imshow(f); 
  figure; imshow(f2); 
   

• Median filtering for salt & pepper noise

  
  f1=imnoise(f,'salt & pepper',0.05); 
  f2=medfilt2(f1,[3 3]); 
  figure; imshow(f); 
  figure; imshow(f1); 
  figure; imshow(f2); 
   

• Adaptative filtering for speckle noise
   

   

  f1=imnoise(f,'speckle',0.04); 
  f2=wiener2(f1,[7 7]); 
  figure; imshow(f); 
  figure; imshow(f1); 
  figure; imshow(f2); 
   

• Edge reinforcement
   

   

  f=imread('saturn.tif'); 
  f=in2gray(f,gray(256)); 
  f1=filter2(fspecial('unsharp'),f); 
  figure; imshow(f); 
  figure; imshow(f1);

Edge detection

Click  here to see an example on edge detection using Canny filters.

Segmentation

Segmentation can be done using different techniques. Here are ones of the most simple

1. Thresholding
   All values above some threshold are set to 1 and those below the threshold are set to 0 or viceversa.

  
  
  
  
  
  
  
 

Morphological operations

 

Image processing practices by other people

    Dong Zhao image processing homework 
    Chong Sze Tong filtering tutorial 
    Manoj Varghese image filtering homework