Matlab Programming

MATLAB is an integrated environment that is used for solving many problems in scientific domain. Matlab is an abbreviation of the word – MATrix LABoratory. It is designed to perform matrix operations. Since images are 2D and 3D in nature, Matlab is suitable for performing image manipulations. It is a powerful tool and can be used to implement imaging projects effectively. MATLAB package comes with several functions that facilitate the image processing. The power of Matlab comes with a set of functions called toolboxes. A toolbox is a collection of functions that are designed to do image processing. Using these functions it is easier to load, save, and perform custom functions on images. The image processing toolbox allows us to do various tasks such as 

1.      Reading/ Writing of images

2.      Colour space and format conversions

3.      Visualization of images

4.      Image filtering

5.      Image Transforms

6.      Image arithmetic

7.      Morphological operations

Matlab also comes with toolboxes for statistics, Wavelet, Neural networks. Apart from the official toolboxes, many public domain toolboxes are available to perform various tasks.

The subsequent sections discuss some of the basic programming capability of Matlab.

Basics of Matlab

Based on the versions of Matlab, the visual appearance of Matlab varies. But essentially the matlab environment has three basic windows.

1. Command Windows – In this main window, the matlab command prompt is >>. Matlab is an interactive environment. The commands can be given one by one and can see the execution interactively. While this acceptable for simple scripts, for length programs this is not acceptable. To facilitate the project development, matlab provides an interactive environment window called Edit Window.

2. Edit Window – In this window, the programs can be written. The programs can be entered, edited and executed in this window. One can create programs with .m extension, called ‘M-files’. One can compare this with C program or Java program. A M-file is a collection of programming commands just like C programs to accomplish a task.

3. Graphics Window – This is the window used to visualize the graphics that are created in the matlab programs. This window also provides help to manipulate the figures such as changing the labels for axes, titles and for exporting the graphics.

Matlab comes with lot of on-line help – Commands such as helpdesk, help, lookfor and helpwin to locate the context oriented help.

Data Types

The Fundamental data type of Matlab is an array. Since images are matrices, Matlab is a perfect fit for manipulating images. Apart from array, the matlab also supports many data types such as integers, double (including floating point numbers), character strings, structures and cells.

File Types

Matlab provides three types of files.

M-Files – These are flat Ascii files that can be created using Matlab environment or any other text editors. M-files are collection of matlab commands. There are two types of M-files.

            1. Script Files

            2. Function Files

Mat-Files – These are native data files with .mat extension. These are binary data that is created with the command save in the command prompt. In the next session, the data can be loaded with the command load. The advantage is that this is a special format that can be read only with Matlab.

Mex-File – These are programs that can call other programming files that are written in programming languages like C or Fortran.

Matlab can be launched by double-clicking the Matlab icon or by navigating the matlab program.

Some simple commands like this can be tested

                        >> 2 + 4

Matlab promptly returns the result. Instead of typing the Matlab commands one by one, a set of valid Matlab commands can be put into a file with .m extension. This is called script file. The script file can be executed in the command prompt by typing the file name. Caution should be taken that the scripts that are created by us should have a name that is different from the key words of matlab.

Function files are similar to script files. The main difference is all the variables defined in a function file are local. One can compare the function file with the function or subroutine of a programming language.

The structure of    a function file is given as

function [outvariables] = functionname (Inputparameters)

% First comment – This appears when the on-line help is sought called h1 line

% This symbol is an indication of comments that are used for documentation

Statement(s)

As the structure is given above, functions are created with the keyword function   followed by the output variables. The definition of a function includes function name also with the necessary input arguments. The first line is called h1 line which is used to provide the comments that are useful to understand the reason for creating that function. This can be followed by many comments. Comments are useful to document the problem so that the third person can understand it. The core part of the function is a set of commands that are necessary to execute it.

Once a function is defined, the function can be executed by invoking it either by calling the function with or without the output variables. If the objective of the function is to display a graph, then storing the resultant may not be necessary. In that case, the output variables can be ignored.      

The functions can be called within another function also as

            t1 = @samplefunction

Here a handle t1 is created for a sample function. This can be used inside for another function as

            t2 = samplefunction2(t1,inputvariables)

                                    or

            t2 = samplefunction2(@samplefunction,inputvariables)

This command is very useful in image processing for block processing.

To increase the speed of execution, the functions can be compiled readily. So the function is interpreted and translated to an intermediate form and stored for later uses. This approach increased the fastness. Also parsing creates a file that is in protected form. This is immensely useful to protect the identity of the programming code so that no one else can see the code and alter it.

Matlab Programming language

Like any programming language like C or C++, Matlab uses a set of programming constructs. This covers all the programming aspects like sequence, branching and repetition. Like other programming languages the global variables can be created with the command global.

                        For example, the statement  

global  m,n

creates two global variables.

Matlab provides branching using if-elseif-else statements as

                                    if condition1

                                       statement(s)

                                   else if condition2

                                       statement(s)

                                   else

   statement(s)

                                   end

Like C language, the conditions are executed and a logical branching is done. Based on the truth of the conditions the statements are executed. Else- Part is executed when both the conditions are false. When conditions become more, writing a recursive if statements with more conditions makes program difficult to understand. In that case the ‘switch’ statement is very helpful. The syntax for switch is given as

                                    switch indicator

                                    case value1

                                            Statement(s)

                                    case value2

                                            Statement(s)

                                    otherwise

                                           Defaultstatement(s)

                                    end

Similarly matlab provides two repetitive statements for and while. The syntax for the statement for is given as

                                    for counter:=start:Increment:End

                                                statement(s)

                                    End

The increment can be positive or negative also. Initially the counter value is the start. then based on the increment , the counter value is incremented. Once it reaches the end value, the program control comes out of the loop. The statements are thus executed as per the start, end and increment value.

Similarly the statement ‘while’ can be used to execute a set of statements till the condition specified by the while is no longer satisfied. The syntax for ‘while’ statement is given as

                        while Condition

                                    Statement(s)

                        end

Other statements that are useful are

            1. Break – This statement is useful break out of a loop in a loop. This is applicable even if the condition of the execution is true.

            2. Return – This statement simply returns the program control to the function that invokes it.

            3. error – When there is a problem inside a function or script, this statement is useful for returning the control to the keyboard.

            4. pause – This statement is useful to temporarily pause or halt the current process and waits for the user response. pause(n) cause the ‘pause’ command  to last n seconds and resume after that.

Input Commands

The syntax of the input command is

                        output_variable = input(‘String’)

For example, the command   n = input(‘Enter the value of n’) display the string in the command prompt expecting the user response. The user response is collected and stored in the variable ‘n’.

When the choices are more, the command ‘menu’ can be used. The syntax of te command ‘menu’ is given as follows.

menu(‘Menu name’,’choice 1’, ‘choice 2’,….’choice n’)

This command creates an onscreen menu. Based on the response given by the user either by mouse or keyboard, the respective commands can be used. Matlab provides advanced versions of GUI. One can use ‘menu’ command for some primitive level usage.

When the data becomes more, the files can be used. One can compare this with C programming language and some of its commands are shown below

1. fpen – Open a file or create a new file.

2. fclose – Close the file

3. fread – read the value from the file

4. fwrite – write the value into the file

5. fscanf – read the formatted data

6. fwritef – write the formatted data  

Some of the advanced matlab usages are structures and cells that can be referred in Matlab manuals.

                                                            ************

Designing and Implementing Linear Filters in the Spatial Domain

Overview

Filtering is a technique for modifying or enhancing an image. For example, you can filter an image to emphasize certain features or remove other features. Image processing operations implemented with filtering include smoothing, sharpening, and edge enhancement. Filtering is a neighborhood operation, in which the value of any given pixel in the output image is determined by applying some algorithm to the values of the pixels in the neighborhood of the corresponding input pixel. A pixel's neighborhood is some set of pixels, defined by their locations relative to that pixel.  Linear filtering is filtering in which the value of an output pixel is a linear combination of the values of the pixels in the input pixel's neighborhood.

Convolution

Linear filtering of an image is accomplished through an operation called convolution. Convolution is a neighborhood operation in which each output pixel is the weighted sum of neighboring input pixels. The matrix of weights is called the convolution kernel, also known as the filter. A convolution kernel is a correlation kernel that has been rotated 180 degrees.

For example, suppose the image is

A = [17  24   1   8  15

     23   5   7  14  16

      4   6  13  20  22

     10  12  19  21   3

     11  18  25   2   9]

and the convolution kernel is

h = [8   1   6

     3   5   7

     4   9   2]

The following figure shows how to compute the (2,4) output pixel using these steps:

1. Rotate the convolution kernel 180 degrees about its center element.
2. Slide the center element of the convolution kernel so that it lies on top of the (2,4) element of A.
3. Multiply each weight in the rotated convolution kernel by the pixel of A underneath.
4. Sum the individual products from step 3.

Performing Linear Filtering of Images Using imfilter

Filtering of images, either by correlation or convolution, can be performed using the toolbox function imfilter. This example filters an image with a 5-by-5 filter containing equal weights. Such a filter is often called an averaging filter.

I = imread('coins.png');

h = ones(5,5) / 25;

I2 = imfilter(I,h);

imshow(I), title('Original Image');

figure, imshow(I2), title('Filtered Image')

Data Types

The imfilter function handles data types similarly to the way the image arithmetic functions do The output image has the same data type, or numeric class, as the input image. The imfilter function computes the value of each output pixel using double-precision, floating-point arithmetic. If the result exceeds the range of the data type, the imfilter function truncates the result to that data type's allowed range. If it is an integer data type, imfilter rounds fractional values.

Because of the truncation behavior, you might sometimes want to consider converting your image to a different data type before calling imfilter. In this example, the output of imfilter has negative values when the input is of class double.

A = magic(5)

A =

    17    24     1     8    15

    23     5     7    14    16

     4     6    13    20    22

    10    12    19    21     3

    11    18    25     2     9

h = [-1 0 1]

h =

    -1     0     1

imfilter(A,h)

ans =

    24   -16   -16    14    -8

     5   -16     9     9   -14

     6     9    14     9   -20

    12     9     9   -16   -21

    18    14   -16   -16    -2

Notice that the result has negative values. Now suppose A is of class uint8, instead of double.

A = uint8(magic(5));

imfilter(A,h)

ans =

    24     0     0    14     0

     5     0     9     9     0

     6     9    14     9     0

    12     9     9     0     0

    18    14     0     0     0

Since the input to imfilter is of class uint8, the output also is of class uint8, and so the negative values are truncated to 0. In such cases, it might be appropriate to convert the image to another type, such as a signed integer type, single, or double, before calling imfilter.

However, if you want to perform filtering using convolution instead, you can pass the string 'conv' as an optional input argument to imfilter. For example:

A = magic(5);

h = [-1 0 1]

imfilter(A,h)   % filter using correlation

ans =

    24   -16   -16    14    -8

     5   -16     9     9   -14

     6     9    14     9   -20

    12     9     9   -16   -21

    18    14   -16   -16    -2

imfilter(A,h,'conv')   % filter using convolution

ans =

   -24    16    16   -14     8

    -5    16    -9    -9    14

    -6    -9   -14    -9    20

   -12    -9    -9    16    21

   -18   -14    16    16     2

Boundary Padding Options

When computing an output pixel at the boundary of an image, a portion of the convolution or correlation kernel is usually off the edge of the image

When the Values of the Kernel Fall Outside the Image

The imfilter function normally fills in these off-the-edge image pixels by assuming that they are 0. This is called zero padding and is illustrated in the following figure.

When you filter an image, zero padding can result in a dark band around the edge of the image, as shown in this example.

I = imread('eight.tif');

h = ones(5,5) / 25;

I2 = imfilter(I,h);

imshow(I), title('Original Image');

figure, imshow(I2), title('Filtered Image with Black Border')

To eliminate the zero-padding artifacts around the edge of the image, imfilter offers an alternative boundary padding method called border replication. In border replication, the value of any pixel outside the image is determined by replicating the value from the nearest border pixel. 

Replicated Boundary Pixels

To filter using border replication, pass the additional optional argument 'replicate' to imfilter.

I3 = imfilter(I,h,'replicate');

figure, imshow(I3);

title('Filtered Image with Border Replication')

The imfilter function supports other boundary padding options, such as 'circular' and 'symmetric'. See the reference page for imfilter for details.

Multidimensional Filtering

The imfilter function can handle both multidimensional images and multidimensional filters. A convenient property of filtering is that filtering a three-dimensional image with a two-dimensional filter is equivalent to filtering each plane of the three-dimensional image individually with the same two-dimensional filter. This example shows how easy it is to filter each color plane of a truecolor image with the same filter:

1. Read in an RGB image and display it.

2.  rgb = imread('peppers.png');

imshow(rgb);

3. Filter the image and display it.

4.  h = ones(5,5)/25;

5.  rgb2 = imfilter(rgb,h);

figure, imshow(rgb2)

Relationship to Other Filtering Functions

MATLAB^® has several two-dimensional and multidimensional filtering functions. The function filter2 performs two-dimensional correlation, conv2 performs two-dimensional convolution, and convn performs multidimensional convolution. Each of these filtering functions always converts the input to double, and the output is always double. These other filtering functions always assume the input is zero padded, and they do not support other padding options.

Filtering an Image with Predefined Filter Types

The fspecial function produces several kinds of predefined filters, in the form of correlation kernels. After creating a filter with fspecial, you can apply it directly to your image data using imfilter. This example illustrates applying an unsharp masking filter to a grayscale image. The unsharp masking filter has the effect of making edges and fine detail in the image more crisp.

I = imread('moon.tif');

h = fspecial('unsharp');

I2 = imfilter(I,h);

imshow(I), title('Original Image')

figure, imshow(I2), title('Filtered Image')

Designing Linear Filters in the Frequency Domain

FIR Filters

The Image Processing Toolbox™ software supports one class of linear filter: the two-dimensional finite impulse response (FIR) filter. FIR filters have a finite extent to a single point, or impulse. All the Image Processing Toolbox filter design functions return FIR filters.

FIR filters have several characteristics that make them ideal for image processing in the MATLAB^® environment:

• FIR filters are easy to represent as matrices of coefficients.
• Two-dimensional FIR filters are natural extensions of one-dimensional FIR filters.
• There are several well-known, reliable methods for FIR filter design.
• FIR filters are easy to implement.
• FIR filters can be designed to have linear phase, which helps prevent distortion.

Another class of filter, the infinite impulse response (IIR) filter, is not as suitable for image processing applications. It lacks the inherent stability and ease of design and implementation of the FIR filter. Therefore, this toolbox does not provide IIR filter support.

Frequency Transformation Method

The frequency transformation method transforms a one-dimensional FIR filter into a two-dimensional FIR filter. The frequency transformation method preserves most of the characteristics of the one-dimensional filter, particularly the transition bandwidth and ripple characteristics. This method uses a transformation matrix, a set of elements that defines the frequency transformation.

The toolbox function ftrans2 implements the frequency transformation method. This function's default transformation matrix produces filters with nearly circular symmetry. By defining your own transformation matrix, you can obtain different symmetries. (See Jae S. Lim, Two-Dimensional Signal and Image Processing, 1990, for details.)

The frequency transformation method generally produces very good results, as it is easier to design a one-dimensional filter with particular characteristics than a corresponding two-dimensional filter. For instance, the next example designs an optimal equiripple one-dimensional FIR filter and uses it to create a two-dimensional filter with similar characteristics. The shape of the one-dimensional frequency response is clearly evident in the two-dimensional response.

b = remez(10,[0 0.4 0.6 1],[1 1 0 0]);

h = ftrans2(b);

[H,w] = freqz(b,1,64,'whole');

colormap(jet(64))

plot(w/pi-1,fftshift(abs(H)))

figure, freqz2(h,[32 32])

Frequency Sampling Method

The frequency sampling method creates a filter based on a desired frequency response. Given a matrix of points that define the shape of the frequency response, this method creates a filter whose frequency response passes through those points. Frequency sampling places no constraints on the behavior of the frequency response between the given points; usually, the response ripples in these areas. (Ripples are oscillations around a constant value. The frequency response of a practical filter often has ripples where the frequency response of an ideal filter is flat.)

The toolbox function fsamp2 implements frequency sampling design for two-dimensional FIR filters. fsamp2 returns a filter h with a frequency response that passes through the points in the input matrix Hd. The example below creates an 11-by-11 filter using fsamp2 and plots the frequency response of the resulting filter. (The freqz2 function in this example calculates the two-dimensional frequency response of a filter.

Hd = zeros(11,11); Hd(4:8,4:8) = 1;

[f1,f2] = freqspace(11,'meshgrid');

mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64))

h = fsamp2(Hd);

figure, freqz2(h,[32 32]), axis([-1 1 -1 1 0 1.2])

Notice the ripples in the actual frequency response, compared to the desired frequency response. These ripples are a fundamental problem with the frequency sampling design method. They occur wherever there are sharp transitions in the desired response.

You can reduce the spatial extent of the ripples by using a larger filter. However, a larger filter does not reduce the height of the ripples, and requires more computation time for filtering. To achieve a smoother approximation to the desired frequency response, consider using the frequency transformation method or the windowing method.

Windowing Method

The windowing method involves multiplying the ideal impulse response with a window function to generate a corresponding filter, which tapers the ideal impulse response. Like the frequency sampling method, the windowing method produces a filter whose frequency response approximates a desired frequency response. The windowing method, however, tends to produce better results than the frequency sampling method.

The toolbox provides two functions for window-based filter design, fwind1 and fwind2. fwind1 designs a two-dimensional filter by using a two-dimensional window that it creates from one or two one-dimensional windows that you specify. fwind2 designs a two-dimensional filter by using a specified two-dimensional window directly.

fwind1 supports two different methods for making the two-dimensional windows it uses:

• Transforming a single one-dimensional window to create a two-dimensional window that is nearly circularly symmetric, by using a process similar to rotation
• Creating a rectangular, separable window from two one-dimensional windows, by computing their outer product

The example below uses fwind1 to create an 11-by-11 filter from the desired frequency response Hd. The example uses the Signal Processing Toolbox hamming function to create a one-dimensional window, which fwind1 then extends to a two-dimensional window.

Hd = zeros(11,11); Hd(4:8,4:8) = 1;

[f1,f2] = freqspace(11,'meshgrid');

mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64))

h = fwind1(Hd,hamming(11));

figure, freqz2(h,[32 32]), axis([-1 1 -1 1 0 1.2])

Creating the Desired Frequency Response Matrix

The filter design functions fsamp2, fwind2, and fwind2 all create filters based on a desired frequency response magnitude matrix. Frequency response is a mathematical function describing the gain of a filter in response to different input frequencies.

You can create an appropriate desired frequency response matrix using the freqspace function. freqspace returns correct, evenly spaced frequency values for any size response. If you create a desired frequency response matrix using frequency points other than those returned by freqspace, you might get unexpected results, such as nonlinear phase.

For example, to create a circular ideal lowpass frequency response with cutoff at 0.5, use

[f1,f2] = freqspace(25,'meshgrid');

Hd = zeros(25,25); d = sqrt(f1.^2 + f2.^2) < 0.5;

Hd(d) = 1;

mesh(f1,f2,Hd)

Ideal Circular Lowpass Frequency Response

Note that for this frequency response, the filters produced by fsamp2, fwind1, and fwind2 are real. This result is desirable for most image processing applications. To achieve this in general, the desired frequency response should be symmetric about the frequency origin (f1 = 0, f2 = 0).

Computing the Frequency Response of a Filter

The freqz2 function computes the frequency response for a two-dimensional filter. With no output arguments, freqz2 creates a mesh plot of the frequency response. For example, consider this FIR filter,

h =[0.1667    0.6667    0.1667

    0.6667   -3.3333    0.6667

    0.1667    0.6667    0.1667];

This command computes and displays the 64-by-64 point frequency response of h.

freqz2(h)

Frequency Response of a Two-Dimensional Filter

To obtain the frequency response matrix H and the frequency point vectors f1 and f2, use output arguments

[H,f1,f2] = freqz2(h);

freqz2 normalizes the frequencies f1 and f2 so that the value 1.0 corresponds to half the sampling frequency, or π radians.

For a simple m-by-n response, as shown above, freqz2 uses the two-dimensional fast Fourier transform function fft2. You can also specify vectors of arbitrary frequency points, but in this case freqz2 uses a slower algorithm.

What Is the Image Processing Toolbox?
The Image Processing Toolbox is a collection of functions that extend the capability of the MATLAB^® numeric computing environment. The toolbox supports a wide range of image processing operations, including:

• Spatial image transformations
• Morphological operations
• Neighborhood and block operations
• Linear filtering and filter design
• Transforms
• Image analysis and enhancement
• Image registration
• Deblurring
• Region of interest operations

Many of the toolbox functions are MATLAB M-files, a series of MATLAB statements that implement specialized image processing algorithms. You can view the MATLAB code for these functions using the statement

·         type function_name

•  

You can extend the capabilities of the Image Processing Toolbox by writing your own M-files, or by using the toolbox in combination with other toolboxes, such as the Signal Processing Toolbox and the Wavelet Toolbox.

Read and Display an Image
Clear the MATLAB workspace of any variables and close open figure windows.

·         clear, close all

•  

To read an image, use the imread command. Let's read in a TIFF image named pout.tif (which is one of the sample images that is supplied with the Image Processing Toolbox), and store it in an array named I.

·         I=imread('pout.tif');

•  

Now call imshow to display I.

·         imshow(I)

Check the Image in Memory
Enter the whos command to see how I is stored in memory.

·         whos

MATLAB responds with

·         Name      Size         Bytes  Class

·          

·           I       291x240        69840  uint8 array

·          

·         Grand total is 69840 elements using 69840 bytes

Perform Histogram Equalization
As you can see, pout.tif is a somewhat low contrast image. To see the distribution of intensities in pout.tif in its current state, you can create a histogram by calling the imhist function. (Precede the call to imhist with the figure command so that the histogram does not overwrite the display of the image I in the current figure window.)

·         figure, imhist(I) % Display a histogram of I in a new figure.

·          

·         

·          

Notice how the intensity range is rather narrow. It does not cover the potential range of [0, 255], and is missing the high and low values that would result in good contrast.
Now call histeq to spread the intensity values over the full range, thereby improving the contrast of I. Return the modified image in the variable I2.

·         I2 = histeq(I);   % Equalize I and output in new array I2.

•  

Display the new equalized image, I2, in a new figure window.

·         figure, imshow(I2)  % Display the new equalized image I2.

·          

·         

·          

Call imhist again, this time for I2.

·         figure, imhist(I2) % Show the histogram for the new image I2.

·          

·         

·          

See how the pixel values now extend across the full range of possible values.

Write the Image
Write the newly adjusted image I2 back to disk. Let's say you'd like to save it as a PNG file. Use imwrite and specify a filename that includes the extension 'png'.

·         imwrite (I2, 'pout2.png');

Check the Contents of the Newly Written File
Now, use the imfinfo function to see what was written to disk. Be sure not to end the line with a semicolon so that MATLAB displays the results. Also, be sure to use the same path (if any) as you did for the call to imwrite, above.

·         imfinfo('pout2.png')

•  

MATLAB responds with

·         ans =

•  

Filename:	'pout2.png'
FileModDate:	'03-Jun-1999 15:50:25'
FileSize:	36938
Format:	'png'
FormatVersion:	[]
Width:	240
Height:	291
BitDepth:	8
ColorType:	'grayscale'

Read and Display an Image
Clear the MATLAB workspace of any variables and close open figure windows. Read and display the intensity image rice.tif.

·         clear, close all

·         I = imread('rice.tif');

·         imshow(I)

·          

·         

·          

Use Morphological Opening to Estimate the Background
Notice that the background illumination is brighter in the center of the image than at the bottom. Use the imopen function to estimate the background illumination.

·         background = imopen(I,strel('disk',15));

•  

To see the estimated background image, type

·         imshow(background)

•  

Display the Background Approximation as a Surface
Use the surf command to create a surface display of the background approximation, background. The surf function requires data of class double, however, so you first need to convert background using the double command.

·         figure, surf(double(background(1:8:end,1:8:end))),zlim([0 255]);

·         set(gca,'ydir','reverse');

•  

The example uses MATLAB indexing syntax to view only 1 out of 8 pixels in each direction; otherwise the surface plot would be too dense. The example also sets the scale of the plot to better match the range of the uint8 data and reverses the y-axis of the display to provide a better view of the data (the pixels at the bottom of the image appear at the front of the surface plot).

Subtract the Background Image from the Original Image
Now subtract the background image, background, from the original image, I, to create a more uniform background.

·         I2 = imsubtract(I,background);

•  

Now display the image with its more uniform background.

·         figure, imshow(I2)

·          

·         

·          

Adjust the Image Contrast
The image is now a bit too dark. Use imadjust to adjust the contrast.

·         I3 = imadjust(I2, stretchlim(I2), [0 1]);

•  

Display the newly adjusted image.

·         figure, imshow(I3);

·          

·         

·          

Apply Thresholding to the Image
Create a new binary thresholded image, bw, by using the functions graythresh and im2bw.

·         level = graythresh(I3);

·         bw = im2bw(I3,level);

·         figure, imshow(bw)

·          

·         

·          

Now call the whos command to see what type of array the thresholded image bw is.

·         whos

•  

MATLAB responds with

·         Name                Size         Bytes  Class

·          

·         I                 256x256        65536  uint8 array

·         I2                256x256        65536  uint8 array

·         I3                256x256        65536  uint8 array

·         background        256x256        65536  uint8 array

·         bw                256x256        65536  logical array

·         level               1x1              8  double array

·          

·         Grand total is 327681 elements using 327688 bytes

•  

Determine the Number of Objects in the Image
To determine the number of grains of rice in the image, use the bwlabel function. This function labels all of the connected components in the binary image bw and returns the number of objects it finds in the image in the output value, numobjects.

·         [labeled,numObjects] = bwlabel(bw,4);% Label components.

·         numObjects =

·          

·             80

•  

The accuracy of your results depends on a number of factors, including:

• The size of the objects
• The accuracy of your approximated background
• Whether you set the connectivity parameter to 4 or 8
• Whether or not any objects are touching (in which case they may be labeled as one object) In the example, some grains of rice are touching, so bwlabel treats them as one object.

Examine the Label Matrix
You may find it helpful to take a closer look at labeled to see what bwlabel has created. Use the imcrop command to select and display pixels in a region of labeled that includes an object and some background.
To ensure that the output is displayed in the MATLAB window, do not end the line with a semicolon. In addition, choose a small rectangle for this exercise, so that the displayed pixel values don't wrap in the MATLAB command window.
The syntax shown below makes imcrop work interactively. Your mouse cursor becomes a cross-hair when placed over the image. Click at a position in labeled where you would like to select the upper left corner of a region. Drag the mouse to create the selection rectangle, and release the button when you are done.

·         grain = imcrop(labeled) % Crop a portion of labeled.

•  

We chose the left edge of a grain and got the following results.

·         grain =

·          

·              0     0     0     0     0     0     0    60    60

·              0     0     0     0     0    60    60    60    60

·              0     0     0    60    60    60    60    60    60

·              0     0     0    60    60    60    60    60    60

·              0     0     0    60    60    60    60    60    60

·              0     0     0    60    60    60    60    60    60

·              0     0     0    60    60    60    60    60    60

·              0     0     0     0     0    60    60    60    60

·              0     0     0     0     0     0     0     0     0

·              0     0     0     0     0     0     0     0     0

•  

A good way to view a label matrix is to display it as a pseudo-color indexed image. In the pseudo-color image, the number that identifies each object in the label matrix maps to a different color in the associated colormap matrix. When you view a label matrix as an RGB image, the objects in the image are easier to distinguish.
To view a label matrix in this way, use the label2rgb function. Using this function, you can specify the colormap, the background color, and how objects in the label matrix map to colors in the colormap.

·         RGB_label = label2rgb(labeled, @spring, 'c', 'shuffle');

·         imshow(RGB_label);

·          

·         

·          

9. Measure Object Properties in the Image
The regionprops command measures object or region properties in an image and returns them in a structure array. When applied to an image with labeled components, it creates one structure element for each component. Use regionprops to create a structure array containing some basic properties for labeled.

·         graindata = regionprops(labeled,'basic')

•  

MATLAB responds with

·         graindata =

·          

·         80x1 struct array with fields:

·             Area

·             Centroid

·             BoundingBox

•  

To find the area of the component labeled with 51's, use dot notation to access the Area field in the 51st element in the graindata structure array. Note that structure field names are case sensitive, so you need to capitalize the name as shown.

·         graindata(51).Area

•  

returns the following results

·         ans =

·          

·            296

•  

To find the smallest possible bounding box and the centroid (center of mass) for the same component, use this code:

·         graindata(51).BoundingBox, graindata(51).Centroid

·         ans =

·          

·           142.5000   89.5000   24.0000   26.0000

·         ans =

·          

·           155.3953  102.1791

•  

To create a new vector, allgrains, which holds just the area measurement for each grain, use this code:

·         allgrains = [graindata.Area];

·         whos allgrains

•  

Call the whos command to see how MATLAB allocated the allgrains variable.

·           Name            Size         Bytes  Class

·          

·           allgrains       1x80           640  double array

·          

·         Grand total is 80 elements using 640 bytes

•  

allgrains is a one-row array of 80 elements, where each element contains the area measurement of a grain. Check the area of the 51st element of allgrains.

·         allgrains(51)

•  

returns

·         ans =

·          

·            296

•  

Compute Statistical Properties of Objects in the Image
Now use MATLAB functions to calculate some statistical properties of the thresholded objects. First use max to find the size of the largest grain. (If you have followed all of the steps in this exercise, the "largest grain" is actually two grains that are touching and have been labeled as one object).

·         max(allgrains)

•  

returns

·         ans =

·          

·            695

•  

Use the find command to return the component label of this large-sized grain.

·         biggrain = find(allgrains==695)

•  

returns

·         biggrain =

·          

·             68

•  

Find the mean grain size.

·         mean(allgrains)

•  

returns

·         ans =

·          

·            249

•  

Make a histogram containing 20 bins that show the distribution of rice grain sizes.

·         hist(allgrains,20)

·          

·         

References:
www.mathworks.com/products/image/
www.mathtools.net/MATLAB/Image_Processing/
www.amath.colorado.edu/courses/4720/2000Spr/Labs/Worksheets/Matlab_tutorial/matlabimpr.html
www.imageprocessingplace.com/DIPUM/dipum_book_description/book_description.htm