## Description

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