Math 584, Mathematics of Medical Imaging

Math 584, Mathematics of Medical Imaging. 1

Matlab Lessons. 1


Matlab Lesson 2: Images and the Radon Transform.. 11

Images and Phantoms. 11

The Radon and inverse Radon transforms. 12

Exercises (Lesson 2) 12

Matlab tips
. 12


Matlab Lessons


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. 

Matlab Lesson 2: Images and the Radon Transform

Images and Phantoms


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.


The Radon and inverse Radon transforms


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:

 ??? Error using ==> minus
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);


Exercises (Lesson 2)


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.


Matlab tips