Lab #6: Image Processing: FIR Filtering of Images

ECES-352

Lab    #6:        Image Processing: FIR Filtering of Images
        (Lab Report Due at beginning of next lab)
This    is    the    official    Lab    #6    description.        Formal    Lab    Report:        You    must    write    a    formal    
lab    report    that    describes    your    work    in    Section    3.    You    should    read    the    Pre-Lab    
section    of    the    lab    and    do    all    the    exercises    in    the    Pre-Lab    section    before
your    assigned    lab    time.    
Important:        When    it    instructs    you    to    get    an    updated    matlab    file    (like    
specgram.m),        please refer to the zipfile on the following site
http://dspfirst.gatech.edu/matlab/toolbox/
Introduction
The    objective    in    this    lab    is    to    introduce    digital    images    as    a    second    useful    signal    type.    
We    will    show    how    the    A-to-D    sampling    and    the    D-to-A    reconstruction    processes    are    
carried    out    for    digital    images.    In    particular,    we    will    show    a    commonly    used    method    
of    image    zooming    (reconstruction)    that    gives    “poor”    results—a    later    lab    will    revisit    
this    issue    and    do    a    better    job.
1 Pre-lab
1.1 Digital    Images    
In    this    lab,    we    introduce    digital    images    as    a    signal    type    for    studying    the    effect    of    
sampling,    aliasing    and    reconstruction.        An    image    can    be    represented    as    a    function
x(t1,    t2)    of    two    continuous    variables    representing
the horizontal (t2) and vertical (t1) coordinates of a point in space.1 For monochrome images, the signal
x(t1,t2) would be a scalar function of the two spatial variables, but for color images the function x(·, ·)
would have to be a vector-valued function of the two variables.2 Moving images (such as TV) would add a
time variable to the two spatial variables.
Monochrome images are displayed using black and white and shades of gray, so they are called grayscale images. In this lab we will consider only sampled gray-scale still images. A sampled gray-scale still
image would be represented as a two-dimensional array of numbers of the form
x[m, n] = x(mT1, nT2) 1 ⌅ m ⌅ M, and 1 ⌅ n ⌅ N
where T1 and T2 are the sample spacings in the horizontal and vertical directions. Typical values of M and
N are 256 or 512; e.g., a 512 ⇤ 512 image which has nearly the same resolution as a standard TV image. In
MATLAB we can represent an image as a matrix, so it would consist of M rows and N columns. The matrix
entry at (m, n) is the sample value x[m, n]—called a pixel (short for picture element).
An important property of light images such as photographs and TV pictures is that their values are
always non-negative and finite in magnitude; i.e.,
0 ⌅ x[m, n] ⌅ Xmax
This is because light images are formed by measuring the intensity of reflected or emitted light which must
always be a positive finite quantity. When stored in a computer or displayed on a monitor, the values of
x[m, n] have to be scaled relative to a maximum value Xmax. Usually an eight-bit integer representation is
used. With 8-bit integers, the maximum value (in the computer) would be Xmax = 28  1 = 255, and there
would be 28 = 256 gray levels for the display, from 0 to 255.
1.2 Displaying Images
As you will discover, the correct display of an image on a computer monitor can be tricky, especially after
some processing has been performed on the image. We have provided the function show img.m in the
DSP First Toolbox to handle most of these problems,3 but it will be helpful if the following points are
⇤⇥⌅
CD-ROM
show img.m noted:
1. All image values must be non-negative for the purposes of display. Filtering may introduce negative
values, especially if differencing is used (e.g., a high-pass filter).
2. The default format for most gray-scale displays is eight bits, so the pixel values x[m, n] in the image
must be converted to integers in the range 0 ⌅ x[m, n] ⌅ 255 = 28  1.
3. Displaying images on the monitor should be done using the show img function.4 The show img
function will handle the color map and the “true” size of the image. The appearance of the image
can be altered by running the pixel values through a “color map.” In our case, we want “grayscale
display” where all three primary colors (red, green and blue, or RGB) are used equally, creating what
is called a “gray map.” In MATLAB the gray color map is set up via
colormap(gray(256))
1
The variables t1 and t2 do not denote time, they represent spatial dimensions. Thus, their units would be inches or some other
unit of length. 2
For example, an RGB color system needs three values at each spatial location: one for red, one for green and one for blue. 3
If you have the MATLAB Image Processing Toolbox, then the function imshow.m can be used instead. 4
If the MATLAB function imagesc.m is used to display the image, two features will be missing: (1) the color map may be
incorrect because it will not default to gray, and (2) the size of the image will not be a true pixel-for-pixel rendition of the image on
the computer screen.
2
which gives a 256⇤3 matrix where all 3 columns are equal. The function colormap(gray(256))
creates a linear mapping, so that each input pixel amplitude is rendered with a screen intensity proportional to its value (assuming the monitor is calibrated). For our lab experiments, non-linear color
mappings would introduce an extra level of complication, so we won’t use them.
4. When the image values lie outside the range [0,255], or when the image is scaled so that it only
occupies a small portion of the range [0,255], the display may have poor quality. In this lab, we will
use show img.m to automatically rescale the image: This requires a linear mapping of the pixel
values:5
xs[m, n] = µx[m, n] +
The scaling constants µ and  can be derived from the min and max values of the image, so that all
pixel values are recomputed via:
xs[m, n] =
⇤
255.999 x[m, n]  xmin
xmax  xmin
⇥⌅
where ⇧x⌃ is the floor function, i.e., the greatest integer less than or equal to x.
1.3 MATLAB Function to Display Images
You can load the images needed for this lab from *.mat files. Any file with the extension *.mat is in
⇤⇥⌅
CD-ROM
IMAGE DATA FILES
MATLAB format and can be loaded via the load command. To find some of these files, look for *.mat
in the DSP First toolbox or in the MATLAB directory called toolbox/matlab/demos. Some of the
image files are named lenna.mat, echart.mat and zone.mat, but there are others within MATLAB’s
demos. The default size is 256 ⇤ 256, but alternate versions are available as 512 ⇤ 512 images under names
such as lenna512.mat and zone512.mat. After loading, use the command whos to determine the
name of the variable that holds the image and its size.
Although MATLAB has several functions for displaying images on the CRT of the computer, we have
written a special function show img() for this lab. It is the visual equivalent of soundsc(), which
⇤⇥⌅
CD-ROM
show im.m we used when listening to speech and tones; i.e., show img() is the “D-to-C” converter for images. This
function handles the scaling of the image values and allows you to open up multiple image display windows.
Here is the help on show img:
function [ph] = show_img(img, figno, scaled, map)
%SHOW_IMG display an image with possible scaling
% usage: ph = show_img(img, figno, scaled, map)
% img = input image
% figno = figure number to use for the plot
% if 0, re-use the same figure
% if omitted a new figure will be opened
% optional args:
% scaled = 1 (TRUE) to do auto-scale (DEFAULT)
% not equal to 1 (FALSE) to inhibit scaling
% map = user-specified color map
% ph = figure handle returned to caller
%----
Notice that unless the input parameter figno is specified, a new figure window will be opened.
5
The MATLAB function show img has an option to perform this scaling while making the image display.
3
1.4 Get Test Images
In order to test your understanding of image display, do the following simple exercises:
(a) Load and display the 326 ⇤ 426 “lighthouse” image from lighthouse.mat. This image can
be downloaded from Web-CT. The MATLAB command load lighthouse will put the sampled
image into the array ww. Use whos to check the size of ww after loading. When you display the image
it might be necessary to set the colormap via colormap(gray(256)).
(b) Use the colon operator to extract the 200th row of the “lighthouse” image, and make a plot of that row
as a one-dimensional discrete-time signal.
ww200 = ww(200,:);
Observe that the range of signal values is between 0 and 255. Which values represent white and which
ones black? Can you identify the region where the 200th row crosses the fence?
2 Warm-up
The instructor verification sheet may be found at the end of this lab.
2.1 Synthesize a Test Image
In order to probe your understanding of the relationship between MATLAB matrices and image display, you
can generate a synthetic image from a mathematical formula.
(a) Generate a simple test image in which all of the columns are identical by using the following outer
product:
xpix = ones(256,1)*cos(2*pi*(0:255)/16);
Display the image and explain the gray-scale pattern that you see. How wide are the bands in number
of pixels? How can you predict that width from the formula for xpix?
(b) In the previous part, which data value in xpix is represented by white? which one by black?
(c) Explain how you would produce an image with bands that are horizontal. Give the formula that would
create a 480⇤480 image with 6 horizontal black bands separated by white bands. Write the MATLAB
code to make this image and display it.
Instructor Verification (separate page)
(d) Generate another test image using the following MATLAB command:
xpix = cos(2*pi*(0:255)/16)’*cos(2*pi*(0:255)/16);
Explain what you see.
4
2.2 Printing Multiple Images on One Page
The phrase “what you see is what you get” can be elusive when dealing with images. It is VERY TRICKY
to print images so that the hard copy matches exactly what is on the screen, because there is usually some
interpolation being done by the printer or by the program that is handling the images. One way to think
about this in signal processing terms is that the screen is one kind of D-to-A and the printer is another kind,
but they use different (D-to-A) reconstruction methods to get the continuous-domain (analog) output image
that you see.
Furthermore, if you try to put two images of different sizes into subplots of the same MATLAB figure, it
won’t work because MATLAB wants to force them to be the same size. Therefore, you should display your
images in separate MATLAB Figure windows. In order to get a printout with MULTIPLE IMAGES ON
THE SAME PAGE, use the following procedure:
1. In MATLAB, use show img and trusize to put your images into separate figure windows at the
correct pixel resolution.
2. Use the Windows program called PAINT to assemble the different images onto one page. This program can be found under Accessories.
3. For each MATLAB figure window, do ALT-PRINT-SCREEN which will copy the active window
contents to the clipboard.
4. After each “window capture” in step 3, paste the clipboard contents into PAINT.6
5. Arrange the images so that you can make a comparison for your lab report.
6. Print the assembled images from PAINT to a printer.
2.3 Sampling of Images
Images that are stored in digital form on a computer are sampled images because they are stored in an
M ⇤ N array (i.e., a matrix). The sampling rate in the two spatial dimensions is chosen at the time the
image is digitized (in units of samples per inch if the original was a photograph). For example, an image
may have been “sampled” by a scanner where the resolution was chosen to be 300 dpi (dots per inch).7 If
we want a different sampling rate, we can simulate a lower sampling rate by simply throwing away samples
in a periodic way. For example, if every other sample is removed, the sampling rate will be halved (in
our example, the 300 dpi image would become a 150 dpi image). Usually this is called sub-sampling or
down-sampling.
8
Down-sampling throws away samples and, therefore, will reduce the size of an image. Downsampling by a factor of p may be accomplished with the following command:
wp = ww(1:p:end,1:p:end);
(a) Apotential problem with down-sampling is that aliasing may occur. This can be illustrated in a dramatic fashion with the lighthouse image.
6
An alternative is to use the free program called IRFANVIEW, which can do image editing and also has screen capture capability.
It can be obtained from http://stud4.tuwien.ac.at/˜e9227474/english.htm. 7
For this example, the sampling periods would be T1 = T2 = 1/300 inches. 8
The Sampling Theorem applies to digital images, so there is a Nyquist Rate that depends on the maximum spatial frequency in
the image.
5
Load the lighthouse.mat file which has the image stored in a variable called ww. When you
check the size of the image, you’ll find that it is not square. Now down-sample the lighthouse
image by a factor of 2. What is the size of the down-sampled image? Notice the aliasing in the downsampled image, which is surprising since no new values are being created by the down-sampling
process. Describe how the aliasing appears visually.
9 Which parts of the image show the aliasing
effects most dramatically?
Instructor Verification (separate page)
3 Lab Exercises: Sampling, Aliasing and Reconstruction
3.1 Down-Sampling
For the lighthouse picture, downsampled by two in the warm-up section:
(a) Describe how the aliasing appears visually. Compare the original to the downsampled image. Which
parts of the image show the aliasing effects most dramatically?
(b) This part is challenging: explain why the aliasing happens in the lighthouse image by using a
“frequency domain” explanation. In other words, estimate the frequency of the features that are being
aliased. Give this frequency as a number in cycles per pixel. (Note that the fence provides a sort of
“spatial chirp” where the spatial frequency increases from left to right.) Can you relate your frequency
estimate to the Sampling Theorem?
You might try zooming in on a very small region of both the original and downsampled images.
3.2 Reconstruction of Images
When an image has been sampled, we can fill in the missing samples using a procedure known as interpolation. For images, this would be analogous to the examples shown in Chapter 4 for sine-wave interpolation
which is part of the reconstruction process in a D-to-A converter. We could use a “square pulse” or a
“triangular pulse” or other pulse shapes for the reconstruction.
Repeat Down
 the Columns
EXPANDED
x[m,n] y[m,n] Repeat Along
 the Rows
REPEAT
SAMPLES
Figure 1: 2-D Interpolation by a factor of three, broken down into row and column operations: the gray dots
indicate repeated data values created by a zero-order hold; or, in the case of linear interpolation, they are the
interpolated values.
9
One difficulty with showing aliasing is that we must display the pixels of the image exactly. This almost never happens because
most monitors and printers will perform some sort of interpolation to adjust the size of the image to match the resolution of the
device. In MATLAB we can override these size changes by using the function truesize which is part of the Image Processing
Toolbox. In the DSP First Toolbox, an equivalent function called trusize.m is provided.
6
For these reconstruction experiments, use the stinger image, down-sampled by a factor of 4 (similar
to what you did in Section 2.3). You will have to generate this by loading in the image from stinger.mat
that may be extracted from the zipped file (lab05images.zip) to get the image that will be in the array called
xx. A down-sampled stinger image should be created and stored in the variable xx4. The objective will
be to reconstruct an approximation to the original stinger image, which is 410 ⇤ 544, from the smaller
down-sampled image.
(a) The simplest interpolation would be reconstruction with a square pulse which produces a “zero-order
hold.” Here is a method that works for a one-dimensional signal (i.e., one row or one column of the
image), assuming that we start with a row vector xr1, and the result is the row vector xr1hold.
xr1 = cos(2*pi*(0:7)/8);
L = length(xr1);
nn = ceil((0.999:1:5*L)/5); %<-- Round up to the integer part
xr1hold = xr1(nn);
Plot the vector xr1hold to verify that it is a zero-order hold version derived from xr1. Explain
what values are contained in the indexing vector nn. If xr1hold is treated as an interpolated version
of xr1, then what is the interpolation factor? Your lab report should include an explanation for this
part, but plots are optional—use them if they simplify the explanation.
(b) Now return to the down-sampled stinger image, and process all the rows of xx4 to fill in the
missing points. Use the zero-order hold idea from part (a), but do it for an interpolation factor of 4.
Call the result xholdrows. Display xholdrows as an image, and compare it to the downsampled
image xx4; compare the size of the images as well as their content.
(c) Now process all the columns of xholdrows to fill in the missing points in each column and and call
the result xhold. Compare the result (xhold) to the original image stinger. Include your code
for parts (b) and (c) in the lab report.
(d) Linear interpolation can be done in MATLAB using the interp1 function (that’s “interp-one”).
When unsure about a command, use help.
Its default mode is linear interpolation, which is equivalent to using the ’*linear’ option, but
interp1 can also do other types of polynomial interpolation. Here is an example on a 1-D signal:
n1 = 0:7;
xr1 = cos(2*pi*n1/8);
tti = 0:0.1:7; %-- locations between the n1 indices
xr1linear = interp1(n1,xr1,tti); %-- function is INTERP-ONE
stem(tti,xr1linear)
For the example above, what is the interpolation factor when converting xr1 to xr1linear?
(e) In the case of the stinger image, you need to carry out a linear interpolation operation on both the
rows and columns of the down-sampled image xx4. This requires two calls to the interp1 function,
because one call will only process all the columns of a matrix.10 Name the interpolated output image
xxlinear. Include your code for this part in the lab report.
10Use a matrix transpose in between the interpolation calls. The transpose will turn rows into columns.
7
(f) Compare xxlinear to the original image stinger. Comment on the visual appearance of the “reconstructed” image versus the original; point out differences and similarities. Can the reconstruction
(i.e., zooming) process remove the aliasing effects from the down-sampled stinger image?
(g) Compare the quality of the linear interpolation result to the zero-order hold result. Point out regions
where they differ and try to justify this difference by estimating the local frequency content. In other
words, look for regions of “low-frequency” content and “high-frequency” content and see how the
interpolation quality is dependent on this factor.
(h) Display the stinger image and be sure that it is “truesize” by using the following:
show_img(xx)
truesize
Now expand the display window to full size. What type of interpolation do you think that MATLAB
does when it expands a window holding an image?
Here are a couple of questions to think about: Are edges low frequency or high frequency features? Are
the feathers low frequency or high frequency features? Is the background a low frequency or high frequency
feature?
Comment: You might use MATLAB’s zooming feature to show details in small patches of the output
image. However, be careful because zooming does its own interpolation.
8
3.3 More about Images in MATLAB (Optional)
This section11 is included for those students who might want to relate these MATLAB operations to previous
experience with software such as Photoshop. There are many image processing functions in MATLAB. For
example, try the help command:
help images
for more information, but keep in mind that the Image Processing Toolbox, which is available in the ECE
labs, may not be on your computer.
3.3.1 Zooming in Software
If you have used an image editing program such as Adobe’s Photoshop, you might have observed how well
or how poorly image zooming (i.e., interpolation) is done. For example, if you try to blow up a JPEG file
that you’ve downloaded from the web, the result is usually disappointing. Since MATLAB has the capability
to read lots of different formats, you can apply the image zooming via interpolation to any photograph that
you can acquire. The MATLAB function for reading JPEG images is imread( ) which would be invoked
as follows:
xx = imread(’foo.jpg’,’jpeg’);
Since imread( ) is part of the image processing toolbox, this test can be done in the ECE computer labs,
but may not be possible on your home computer.
3.3.2 Warnings
Images obtained from JPEG files might come in many different formats. Two precautions are necessary:
1. If MATLAB loads the image and stores it as 8-bit integers, then MATLAB will use an internal data
type called uint8. The function show img( ) cannot handle this format, but there is a conversion
function called double( ) that will convert the 8-bit integers to double-precision floating-point for
use with filtering and processing programs.
yy = double(xx);
You can convert back to 8-bit values with the function uint8().
2. If the image is a color photograph, then it is actually composed of three “image planes” and MATLAB
will store it as a 3-D array. For example, the result of whos for a 545 ⇤ 668 color image would give:
Name Size Bytes Class
xx 545x668x3 1092180 uint8 array
In this case, you should use MATLAB’s image display functions such as imshow( ) to see the color
image Or you can convert the color image to gray-scale with the function rgb2gray( ). For more
information on the image processing functions in MATLAB, try help:
help images
11Optional mean that you don’t have to include this in a lab report. This section provided in case you are curious and want to
learn more on your own.
9