Math
584,
Mathematics of Medical Imaging

Matlab Lesson 2:
Images and the Radon Transform

The Radon and
inverse Radon transforms

This is a simple
tutorial for using
Matlab for Math 584, Mathematics of Medical Imaging.

The black text in bold is what you should type at the Matlab command prompt. The non-bold black text which follows is the expected Matlab output. The red text will explain what you should be learning. Please experiment by changing the commands so you can get a feel for the language of Matlab.

In Matlab,
images are real-valued matrices. The
matrix entries represent the intensities or densities at the pixels of
the
image. One way to work with images is to load one from a file using the
imread function.
Download the file circle.jpg
into your Matlab working directory and
then type

>> **X=imread('circle.jpg','jpeg');**

Matlab can read several different image file formats. Type **>>
help imread** for more information. To display the image, type

>> **figure****; imagesc(X); colormap(gray);**

The image X is a
real-valued
matrix. You can view the entries just
like any other matrix:

>> **size(****X), X(100,100)**

ans =

256 256 This means that the image is 256x256 pixels.

ans =

255 The value of the
pixel at location (100,100) is 255.

Experiment by defining
your own
images. For example, try:

>> **figure****; imagesc(eye(150)); colormap(gray);**

A phantom is an image
which can be
used to test the numerical accuracy of a 2-D image reconstruction
algorithm. Matlab
has several built-in phantoms. The
default phantom is a ‘modified’ Shepp-Logan
phantom:

>> **P=phantom(256);
figure; imagesc(P); colormap(gray);**

You can also easily
generate your
own ellipse-based phantoms. Type
>> **help phantom** for more
information.

For more information on
these
topics, see the following help topics: image, imagesc,
imread, imwrite,
phantom.

Begin with a phantom:

>> **P=phantom(256);**

To compute the radon
transform of
this matrix, type

>> **R=radon(P,0:1:359);**

The vector 0:1:359
represents the
discrete set of angles (in degrees) for the Radon transform. Plot the resulting matrix, and spend some time
trying to understand why the image appears as it does.
What are the horizontal and vertical
axes? What symmetries or approximate
symmetries exist? What does the dark
region of the image represent?

>> **figure****; imagesc(R); colormap(gray);**

Now reconstruct the
image using
filtered back projection (this may take several minutes if you have a
slower
computer) :

>> **P2=iradon(R,0:1:359); figure; imagesc(P2);
colormap(gray);**

You can experiment with
different
options and filters for iradon. See >> **help iradon**.

We can compare the reconstructed image to the original by displaying
the difference:

>> **figure; imagesc(P2-P);
colormap(gray);
**

This command produces the following error:

Matrix dimensions must agree.

We check

>> size(P2)

ans =

258 258

The
reconstruction has produced a slightly larger image. We can select a
256X256 part to compare to the original image.

>>figure; imagesc(P2(2:257,2:257)-P); colormap(gray);

1.
Define an
image (matrix) consisting of all zeros except for a 1 at the (m,n)^{th}
pixel. What do you expect the Radon
transform of this data to look like? Now
test your conjecture using Matlab. Try this for a few different pairs (m,n).

2.
Let **R** be the Radon transform of the 'modified' Shepp-Logan phantom as above. What is
the effect of the Matlab statement
>> **R(:,****10:20****)=0**.
How does this change affect the reconstructed
image?

3. In the final step above we
selected P2(2:257,2:257) to compare to P. There are two other
reasonable possibilities. Find them, compare them to P and state which
of the three possibilities is most like the original image.

- Always end Matlab statements with a semicolon, unless you want Matlab to display the resulting value of the statement (or expression). You can get the value of any variable by simply typing its name (press enter without a semicolon).
- You can use the up and down arrows on your keyboard to retrieve previously executed Matlab statements. For more efficient retrieval, type the first few characters of the statement before using the up arrow.
- From the Matlab prompt, type
**help [function name]**to get help on any built-in Matlab function. - It is sometimes tempting to use "for loops" to define and manipulate matrices and arrays. However this should be avoided whenever possible because Matlab is very inefficient when implementing such loops. There are certain Matlab operators and functions which can often be used in place of "for loops". See the following help topics: times, divide, colon, find.
- You can combine a sequence of Matlab statements into a single Matlab-executable file called an m-file. See help function.
- You can save large data objects in files (for future use). See the help topics: save, load.
- You can have more than one statement on a single line: just separate the statements by semicolons.
- Have you discovered other useful tips and tricks? Send them to cle@math.upenn.edu .