# Image reconstruction in Nuclear Medicine

### Pages to are hidden for

"Image reconstruction in Nuclear Medicine"

```					Image reconstruction in Nuclear Medicine

Project 2 : ENEE 739R
Aniruddha Kembhavi
I. Introduction

In this project I have implemented Hudson and Larkin s reconstruction technique for PET (Positron
Emission Tomography) transmission images using the                Ordered Subsets         Expectation
Maximization (OS-EM) algorithm. I also compute a reconstructed image using the simple filtered
back-projection technique with the help of the inbuilt MATLAB functions. This project goes to
demonstrate the effectiveness of using Ordered Subsets as opposed to a basic EM algorithm, with
regards to image quality and reconstruction speed.

II. Using the filtered back-projection technique

The image obtained by the filtered back-projection technique is not of the highest quality. It is
shown below in Fig.1. However, this quality improves when radial filtering is performed prior to the
back-projection. The figure also shows that the choice of the filter is quite crucial to obtain a good
reconstructed image.

Fig.1 Reconstructed images in the absence and presence of radial filtering
III. OS-EM approach                Calculation of the weight matrix A

The outline of my implementation of the OS-EM algorithm is as follows:

1. Generate starting image
2. Loop over : iteration no (till convergence)
3.         Loop over all subsets
4.                 Based on subset no., choose the rays from all the projections
5.                 Generate new image based on old one (hence update old one)
6.         End loop for subsets
7. End loop for iterations

Steps 4 and 5 are used to create a new image based on projection data and the old image. This is the
only place where the A-matrix is used in our calculations. The equation used for this step is given
by:

The matrix A has (o1 * o2 * p1 * p2) number of values, where o1 and o2 are the dimensions of the
output image and p1 and p2 are the dimensions of the projection data. Given that these are 4-byte
floating point values, for our case, this data would require a space upwards of 2 Gb, which is quite
impractical; more so the case if we expect to have all this data in the RAM at once. Hence I took an
approach to calculate the values aij on-the-fly keeping the number of calculations to a minimum.

Now for every projection angle (total of 192 for this case), the slope of the rays is the same. The only
changing factor between the rays is the y-intercept. Hence for every projection angle and every ray,
the A-matrix changes and has to be separately calculated. For this purpose, I have a separate
function:
A-matrix = function (parameters a, b, c of ray where ax + by + c = 0)
The matrix calculated here consists of aij for a particular i and all j. (Here i is the variable for the
projection rays and j is the variable for the image pixels.) The elements of this matrix are then
summed up (after pre-multiplying with the old image pixels) to give us:

Now in the same loop, for the same i the following can be calculated as we have qi the projection
data and the matrix A (just calculated) to give us:

This matrix is then put in an accumulator matrix so that it adds to itself for every i to generate:

Note here, that the memory requirements here are bare-minimum as no matrix has been stored
except 1 accumulator matrix having (o1 x o2) floating point values. At the same time, within this
loop, we add the bare A-matrix to another accumulator matrix so that after we go through all i
values, we arrive at the value:

The last step now becomes trivial where we update each pixel to get:

Thus it can be seen that we do NOT store the A-matrix at any point, and at the same time do not
recalculate its values, so as to reduce the total number of computations.
Fig.2 Operation of the accumulator matrix

IV. OS-EM approach Choice of subsets

The choice of subsets can sometimes be quite crucial to the reconstructed image quality. The
approach while choosing subsets must be to ensure that the data in each subset (the projection data)
is well spread, so as to resemble a quantized version of the original projection data.

Consider the projection data given in Fig.3. We need to choose a few projection values from each
projection angle to form one subset and so on. A good way of doing this (which is also the method I
have implemented), is shown in Fig.4 for choosing 2 subsets. This ensures that:
a. Each subset has information from each projection angle
b. Each subset has information spread over the entire region covered for each projection angle.
The method used to implement this was shown in the previous section. After each subset, the image
obtained was used as the base image for the new subset. At the end of all subsets, the image obtained
was used as the base image for the starting subset of the next iteration. Changing the number of
subsets just involved changing the selection of rays for each subset.
Fig.3 Transmission scan image

Fig.4 Transmission scan divided into 2 subsets

Another way of choosing the subsets is shown in Fig.5 where as a poor method for choosing them is
shown in Fig.6 Obviously this is a poor method as the distribution over the rays is not uniform.

Fig.5 Another method to choose subsets
Fig.6 A poor choice for the subsets

V. Results

Following are some results for varied number of subsets and varied number of iterations. I start by
showing the results for the OS(1)-EM case, which basically boils down to a simple EM algorithm. I
then display results when using 2, 4 and 8 subsets. For each case, I have shown the image after 1
iteration, after convergence and a few results in between to maintain continuity.

Fig.7 Reconstructed images     1 Subset. Iterations    1,5,10,15,20,23
(Top row first, then bottom row)
Fig.8 Images for 2 Subsets. Iterations    1,5,10,15

Fig.9 Images for 4 Subsets. Iterations    1,4,8,11

Fig.10 Images for 8 Subsets. Iterations    1,3,6,8

Watching the images after every iteration is even more pleasing to the eye, when the images are
grouped together to form a video. This shows us how the image changes drastically for the first few
iterations and then it slowly stabilizes and converges to form the final output image. Frame rate is 5
fps. The videos can be viewed at:

http://www.umiacs.umd.edu/~ani/movie_os1.avi
http://www.umiacs.umd.edu/~ani/movie_os2.avi
http://www.umiacs.umd.edu/~ani/movie_os4.avi
http://www.umiacs.umd.edu/~ani/movie_os8.avi
No. of subsets             No. of iterations run        Iterations for convergence

1                            48                               23
2                            24                               15
4                            24                               11
8                            24                                8
16                            24                                5

Table1. Number of iterations required for convergence with the %-change metric

VI. Convergence and the quality of the image

As the number of iterations kept increasing, the current and previous image kept getting closer and
closer. One method I used to quantify this closeness was the following ratio:
Average % change = 100 * (absolute(new_image        old_image)) / old_image
The values for various iterations, for the different number of subsets are shown below.

Fig.11 % change in image Vs Iteration number
However it is not just this measure that helps us gauge convergence. As the number of iterations
increases, the % change reduces, but ALSO, the rate of this change becomes constant. Hence the
slope of the above graphs also played a role in deciding the convergence. The slopes of the above
graphs are shown below.

Fig.12 Slope of the % change in an image Vs Iteration no

Hence it can be seen that by choosing appropriate thresholds for the two corresponding graphs for a
fixed number of subsets can give us a measure of convergence of the image sequence. This is how I
obtained the image number for convergence in the previous table.

Another metric I used to quantify image quality was Chi-square metric. As given by Hudson and
Larkin, the Chi-square of an image after an iteration is given by:

where yi represent the detector values, and ui represent the expected values for cumulated counts y
for each detector. As was expected, the Chi-square values decrease with an increase in the number of
iterations, and slowly converge upon to a minimum. Chi-square values for varied number of subsets
and iterations are shown.

Fig. 13 Chi-square Vs No of iterations

Just like the previous metric that I defined and used, the slope of the Chi-square plot can give us a
good estimate as to where to stop the number of iterations. Hence setting an appropriate threshold
for the slope of the Chi-square plot can give us a measure of the image quality that can be obtained
and hence the convergence.
Fig. 14 Slope of the Chi-square metric Vs No. of iterations

VII. Computational time

Regardless of the number of subsets, the time required to run a complete iteration is always the
same. This is because regardless of the number of subsets, the total number of transmission rays that
have to be processed remains the same. 1 iteration of my program takes between 1 min and 1min,
5sec to run. The machine I am using has a 3 GHz processor. Though each iteration takes the same
amount of time irrespective of the number of subsets, convergence is obtained faster with more
subsets. Hence the overall computational time is reduced as we choose a larger number of subsets.

VIII. Conclusion

In this project a method for PET reconstruction using ordered subsets, proposed by Hudson and
Larkin was implemented the effect of using these ordered subsets was studied. 2 metrics to measure
the image quality with respect to convergence of the iterative procedure were proposed.
This document was created with Win2PDF available at http://www.daneprairie.com.
The unregistered version of Win2PDF is for evaluation or non-commercial use only.

```
DOCUMENT INFO
Shared By:
Categories:
Stats:
 views: 86 posted: 1/12/2010 language: English pages: 13