"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.