Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UFE0015871/00001
## Material Information- Title:
- Penalized Maximum Likelihood Methods in Emission Tomography
- Creator:
- ZAHNEN, JEFFREY A. (
*Author, Primary*) - Copyright Date:
- 2008
## Subjects- Subjects / Keywords:
- Image reconstruction ( jstor )
Mathematical constants ( jstor ) Mathematical vectors ( jstor ) Maximum likelihood estimations ( jstor ) Objective functions ( jstor ) Penalty function ( jstor ) Perceptron convergence procedure ( jstor ) Pets ( jstor ) Positron emission tomography ( jstor ) Tomography ( jstor )
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Jeffrey A. Zahnen. Permission granted to University of Florida to digitize and display this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Embargo Date:
- 3/1/2007
- Resource Identifier:
- 659560804 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads:
zahnen_j ( .pdf )
zahnen_j_Page_54.txt zahnen_j_Page_60.txt zahnen_j_Page_59.txt zahnen_j_Page_23.txt zahnen_j_Page_13.txt zahnen_j_Page_03.txt zahnen_j_Page_37.txt zahnen_j_Page_35.txt zahnen_j_Page_56.txt zahnen_j_Page_31.txt zahnen_j_Page_16.txt zahnen_j_Page_20.txt zahnen_j_Page_33.txt zahnen_j_Page_51.txt zahnen_j_Page_39.txt zahnen_j_Page_19.txt zahnen_j_Page_02.txt zahnen_j_Page_45.txt zahnen_j_Page_10.txt zahnen_j_Page_43.txt zahnen_j_Page_53.txt zahnen_j_Page_07.txt zahnen_j_Page_44.txt zahnen_j_Page_17.txt zahnen_j_Page_41.txt zahnen_j_Page_28.txt zahnen_j_Page_62.txt zahnen_j_Page_11.txt zahnen_j_Page_52.txt zahnen_j_Page_06.txt zahnen_j_Page_30.txt zahnen_j_Page_15.txt zahnen_j_Page_57.txt zahnen_j_Page_50.txt zahnen_j_Page_55.txt zahnen_j_Page_38.txt zahnen_j_Page_09.txt zahnen_j_Page_34.txt zahnen_j_Page_04.txt zahnen_j_Page_42.txt zahnen_j_Page_27.txt zahnen_j_Page_12.txt zahnen_j_Page_26.txt zahnen_j_Page_58.txt zahnen_j_Page_05.txt zahnen_j_Page_01.txt zahnen_j_Page_29.txt zahnen_j_Page_24.txt zahnen_j_Page_49.txt zahnen_j_Page_08.txt zahnen_j_Page_63.txt zahnen_j_Page_25.txt zahnen_j_Page_22.txt zahnen_j_Page_48.txt zahnen_j_Page_21.txt zahnen_j_Page_18.txt zahnen_j_Page_36.txt zahnen_j_Page_32.txt zahnen_j_Page_14.txt zahnen_j_Page_61.txt zahnen_j_Page_40.txt zahnen_j_Page_47.txt zahnen_j_Page_46.txt zahnen_j_pdf.txt |

Full Text |

PENALIZED MAXIMUM LIKELIHOOD RECONSTRUCTION METHODS FOR EMISSION TOMOGRAPHY By JEFFREY A. ZAHNEN A DISSERATATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006 Copyright 2006 by Jeffrey A. Zahnen This dissertation is dedicated to my grandmother, Lena Zahnen, for ah--iv believing in me. ACKNOWLEDGMENTS For his patience and advice through the years, I offer my greatest thanks to my advisor, Dr. Bernard Mair. In addition, I would like to thank Dr. David Gilland, Dr. William Hager, Dr. James Keesling, Dr. Scott McCullough and Dr. David Wilson for taking the time to serve on my committee and providing sl.-.-. -1 i..o:, to improve the work. Finally, I would like to thank my parents for their unwavering support of my studies. TABLE OF CONTENTS page ACKNOWLEDGMENTS ................... ...... iv LIST OF TABLES ...................... ......... vii LIST OF FIGURES ................... ......... viii ABSTRACT ....................... ........... ix CHAPTER 1 MATHEMATICS OF POSITRON EMISSION TOMOGRAPHY (PET) 1 1.1 Introduction To PET ......... ................. 1 1.2 Mathematical Model for PET ....... ............. 2 1.3 PET Probability Function ..... .................. 3 1.4 Reconstruction Techniques ........ .............. 4 1.4.1 Maximum Likelihood ...... ... ......... .. 4 1.4.2 Penalized Maximum Likelihood (PML) ...... ... ... 6 1.5 Minimizing Function ........ .................. 9 2 PML RECONSTRUCTION METHODS .................. 10 2.1 Green's One Step Late Algorithm ...... ......... 10 2.1.1 Lange's Improvement to Green's OSL ..... ... .... 11 2.1.2 Mair-Gilland Method ....... ......... .... 12 2.2 Accelerated Penalized Maximum Likelihood Method (APML) .12 3 GENERAL PML ALGORITHM .................. ..... .15 3.1 General Algorithm Statement .................. .. 15 3.1.1 General Algorithm Update .................. .. 15 3.1.2 Penalty function criteria ............. .. .. .16 3.2 Proof of Convergence .................. ....... .. 17 3.3 New PML Algorithm .................. ....... .. 21 3.3.1 Non-Uniform Step-Size Method ............... .. 22 3.3.2 Line Search Methods .................. ..... 24 3.3.3 Stopping Criteria .................. ... .. 26 4 MULTIPLE IMAGE-MOTION RECONSTRUCTION . ... 27 4.1 Algorithm Desription .................. ....... .. 27 4.1.1 RM Algorithm ...... 4.1.2 Image-Matching Penalty 4.1.3 Modified RM Iteration 4.2 Convergence Results ...... 5 RECONSTRUCTIONS. .................. ........ 36 5.1 PML Reconstructions .................. ....... 36 5.2 Modified RM Reconstructions .................. .. 42 5.3 Conclusions .................. ............ 45 APPENDIX: SAMPLE ALGORITHM FORTRAN CODE . .47 REFERENCES .............. ................... 53 BIOGRAPHICAL SKETCH .................. ......... 54 .................... .................... .................... .................... LIST OF TABLES Table page 4-1 RM voxel neighborhood .................. ........ .. 30 5-1 PML reconstruction data .................. ........ .. 38 5-2 RM data a = 0.02, 0 = 0.005 .................. ... .. 44 LIST OF FIGURES Figure page 1-1 A 2-D PET tube registering an emission ............. 1 1-2 Tube regions ............... ............ .. 3 1-3 Sample neighborhood system with associated weights . . .. 7 1-4 Derivatives of the penalty terms ................ ..... 8 3-1 Armijo rule .................. ............... .. 25 5-1 Phantom chest source image .................. ..... .. 36 5-2 Comparing the decrease in the objective function . ..... 38 5-3 Comparing the root-mean-square (RMS) error . . 39 5-4 Phantom image reconstruction using bisection rule . . 40 5-5 Phantom image reconstruction using Armijo rule . 41 5-6 Line plot of heart region ............... ....... .. 41 5-7 Comparing reconstruction algorithms ............. .. .. 42 5-8 Modified RM function plots ................ .... .... 43 5-9 Quiver plot of Slice 20 ............... ........ .. 44 5-10 RM reconstructions ............... .......... .. 45 Abstract of Disseratation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy PENALIZED MAXIMUM LIKELIHOOD RECONSTRUCTION METHODS FOR EMISSION TOMOGRAPHY By Jeffrey A. Zahnen December 2006 C'!h In: Bernard Mair Major Department: Mathematics Positron emission tomography (or PET) has been utilized for years in medical diagnosis. The PET reconstruction problem consists of recovering the image from the data gathered by the detectors in the machine. Penalized maximum likelihood techniques allow for the recovery of the image with a smoothing term to improve the image. Previous methods using this technique forced the use of a weak penalty term in order to gain convergence. Using a new method, this restriction is removed while retaining convergence. In addition, the method is extended to the problem of recovering multiple PET images along with a vector describing the motion between the images. CHAPTER 1 MATHEMATICS OF POSITRON EMISSION TOMOGRAPHY (PET) 1.1 Introduction To PET Emission tomography has been used for years to study metabolic processes in the body. In this procedure, a patient is injected with a short-lived radioactive isotope. While the substance is being absorbed by the tissues in the body, the patient is placed in a scanner. As the isotope decays, it begins to emit positrons. These emitted positrons travel only a short distance before interacting with a nearby electron. In positron emission tomography (PET), the subsequent annihilation produces a pair of photons traveling in opposite directions on a line oriented in a uniformly random direction. A PET scanner consists of individual detector crystals arranged around the patient in a cylindrical manner. Each pair of detectors forms a virtual detection tube. Figure 1-1 di- -1ph,- a tube formed by detectors a and 3. A near-simultaneous arrival of two photons at opposite ends of a o j/ Figure 1-1: A 2-D PET tube registering an emission the tube is then registered as a single detection. The tube shown in Figure 1-1 will register the di-pl i '1, emission at point B but not the emission at point A. In single photon emission computed tomography (SPECT), a single photon is emitted. In this method (more widely applied than PET due to the relative costs of the scanning machines), the detectors are located only on one side of the patient and a single detection results in a count. This scan is performed for different angles between the patient and detector to provide a full image. In both procedures, the aim is to determine the number of emissions from the various locations in the patient. The resulting intensity images are used to make medical diagnoses. This work will consider the problem of PET image reconstruction using the maximum likelihood (\I b) method with the addition of a penalty term. Although developed for PET reconstruction, the method can be extended for SPECT reconstruction. Unlike previous methods, the revised algorithm presented uses a more general penalty term, allowing for a wider range of settings. Proof of convergence and reconstruction results are obtained using this method. 1.2 Mathematical Model for PET The 2-dimensional PET scanner consists of a circle of detectors surrounding the region of interest, Q, which can be taken to be rectangular. There are I virtual detector tubes, each one formed by two of the detectors composing the scanner ring. f is partitioned into J smaller-sized rectangles, called voxels. An emission in voxel j is detected by tube i if the emitted positrons are detected simultaneously in the two detectors forming the tube. Shepp and Vardi [16] introduced a mathematical model for PET in 1982. In the model, the random emissions from voxel j are modeled by a Poisson process with mean fj. Thus the detector count in tube i, Yi, will also be governed by a Poisson process with mean Ai = pj fj (11) where pij is the geometric probability that an emission in voxel j is registered by detector tube i. Using the properties of the Poisson distribution, the probability that there are yi counts registered in detector i is given by e-Cx i y y'! To recover the emission intensity, it is necessary to estimate f = [f, f2, ..., fJ] based on a random sample from Y = [Y, Y2, ..., Y]. 1.3 PET Probability Function The probability function pij for the 2-dimensional PET scanner can be defined using the geometry of the scanner. Given any tube i formed by two detectors, it can be divided into six regions as shown in Figure 1-2. The voxels B A SR22, F3 \R /AA Figure 1-2: Tube regions are assumed to be so small that the probability can be assumed approximately constant on the entire voxel. If z is the center of voxel j, then pij is defined by the probability that the path of a pair of photons emitted from the center of the voxel, z, intersects both detectors defining the tube. This probability is proportional to the angle-of-view of the tube from z. Thus, if CzD is defined to be the acute angle formed by the lines passing through C and z and D and z, the following definition results [12]. CzD if z e R1 7 -BzC if z R2 7TPi, AzB if z e R3 (1-2) 7 DzA if z e R4 0 if z E R22 or 44 1.4 Reconstruction Techniques The reconstruction problem is to recover the emission intensity f given the scanner data y, assumed to be Poisson-distributed. Several methods have been used to accomplish this reconstruction. If the detectors are treated as points, the model reduces to line integrals along detector lines. Filtered-back projection (FBP) is a popular method utilized to recover the emission intensity, as it provides a quick way to invert the Radon transform. However, it has been seen to produce images with streaking artifacts [16]. In 1998, Carroll [3] decomposed the data, probability matrix and intensity functions into basis functions, transforming the problem into an infinite set of linear equations to solve. Our algorithm though is based on a third method, Maximum likelihood estimation (! I). 1.4.1 Maximum Likelihood In the ML method, a search is made for the emission intensity f that maximizes the probability of obtaining the Poisson-distributed data. Each detector count yi is assumed to be an independent sample from a Poisson distribution with a mean corresponding to i = Pf i = j pijfj. Thus the probability that there are 5 yi counts in detector i is given by P(Y,= Uyt=A)- (13) The detector counts are independent random variables, so these probabilities are multiplied together to yield the probability of obtaining the data from an emission intensity map of f. An ML estimate is an intensity f that maximizes the following probability. --i Y In addition, the intensity map, f is constrained to be non-negative ({f lf > 0}). An equivalent problem is to apply the negative logarithm and perform a minimization of the term -(Pf, y log Pf, log yi!). As yi! does not depend on the unknown intensity f, it can be safety dropped from the term resulting in the ML objective function. L(f) = Pf, y, logPf, (1 4) The Expectation-Maximization (or EM) algorithm developed by Shepp and Vardi [16] can then applied to find the minimizer of this term, resulting in the EMML algorithm. f7+ fni If Z iIf (1 5) Z Pi,j where f" refers to the nth step in the iteration. Alternately, the EMML iteration can be expressed as a descent algorithm. fn+1 fn fin6 f )/VL,. JY i ] *3Qf ' Restated in matrix form yields the equation f"+l f" aD(f")VL(f"). where a = 1 is the step-size and D is the diagonal matrix whose diagonal entries consist of dj,j = fj/ ipj. While the EM algorithm has been shown to monotonically decrease the objective function and be globally convergent, a well-known deficiency is a tendency to produce noisy images if the iteration is run too long [8]. 1.4.2 Penalized Maximum Likelihood (PML) To help reduce the noise, a smoothing term can be added to the likelihood term resulting in the penalized maximum likelihood (PML) method. Equivalently, this method can be considered as an application of MAP (maximum a-posteriori), in that the recovered image is assumed to have a known local Gibbs probability distribution. This method, used by Geman and McClure [7], is employ, 1l by first choosing a penalty function Q and then forming the term U(f) E '(f- fJ) (1-6) i j where w is a weighting factor that associates pixels. This can be used to link neighboring pixels, forcing them to have similar intensities. For example, Ni, can be defined on each voxel i to be the set of pixels in a neighborhood lattice of that voxel. Thus the weight :.' i 0 if j N(i). If j N(i), the weight is defined to be a nonzero constant. One possibility is to let Ni consist of the eight neighboring voxels as shown in Figure 1-3 with associated weights. In addition, a constant value, 7 is multiplied to the penalty term to control the influence of the penalty in the reconstruction. As 7 0, the reconstructions approach those of the EM algorithm. However, if 7 is chosen to be too large, the images will be oversmoothed. Figure 1-3: Sample neighborhood system with associated weights There exist many choices for the penalty function 0(r). A frequently used choice is a penalty of the form 0(r) = log cosh [8],where 6 is a previously chosen constant that allows for some fine-tuning of the penalty. This constant determines how the prior is scaled with respect to r, as the log cosh penalty will smooth all differences above a certain cutoff value, determined by 6. As 6 -+ o, the method will approach results of the EM Algorithm (as no smoothing will occur). This penalty is useful as it approximates quadratic smoothing for small values of r, but its derivative remains bounded. In 1993, Lalush and Tsui [9] developed a new penalty function that allowed for even greater control over the penalized range of intensities. Noting that current PML methods did not require the penalty function to be calculated, they instead built the function from properties of its derivative. Thus, they produced S r [1 + 2( )2]-2 (17) Or Irl 6 r where a and 6 are constants that affect the spread and location respectively of peak of the derivative. As a increases, the peak becomes more narrow and as 6 1/sqrt(2) 1 1/sqrt(2) 1 0 1 1/sqrt(2) 1 1/sqrt(2) 8 increases, the peak away from the origin. The plot of the derivative shows the in i i iiiy of the smoothing occurring in the region under the peak of the derivative [9]. Thus altering the values of a and 6 will affect the range of smoothing in the image. By comparing this function to the derivative of the logcosh function, '(r) = (1/6) tanh r it can be seen that while the logcosh function does not have a peaked derivative 002 09 S0018 08 0016 0016 07 0014 06 0012 05 001- 04 0008 03 0006 02 0004 01 0002 0 0 0 20 40 60 80 100 20 2 40 60 80 100 (a) Derivative of the Lalush-Tsui penalty (b) Derivative of the logcosh penalty Figure 1-4: Derivatives of the penalty terms (it increases to an .i-vmptotic limit), it can be considered as an infinitely-wide peak. Thus, Green's logcosh penalty can be approximated as a Lalush-Tsui penalty by taking a close to zero (producing a near infinitely-wide peak). Integrating the Lalush-Tsui penalty term results in their penalty function: 6 T2 2a4 ( = log[1 2a2(1- )+ a2(6 + r4 2r2(1 -2a2)]. 2a 62 62 While Green's original OSL algorithm does not require the calculation of this term, other algorithms do require the frequent calculation of the objective function. In that case, the Lalush-Tsui penalty will have a larger computation time than other penalties. In addition, the general Lalush-Tsui penalty does not have the strict convexity condition needed for convergence to a unique minimizer. 1.5 Minimizing Function Given a choice of penalty term U and constant 7, the PML objective function thus becomes E(f) = Pf, y logPf, + TU(f) (1 8) In Maximum Likelihood, this objective function must be minimized over the set of positive intensity values {f fj > 0}. This minimum must satisfy the Kuhn-Tucker optimality conditions for the problem. In this situation, those two conditions are that for all j OE (f) > 0 (1-9) afj aE fj- (f)= 0 t9afj It will be shown that the revised algorithm presented in C'i lpter 3 will meet these conditions at convergence. CHAPTER 2 PML RECONSTRUCTION METHODS 2.1 Green's One Step Late Algorithm Applying the Maximum Likelihood method, a search is made to minimize the penalized objective function: E(f) (Pf,- ylogPfi) + TU(f) i Differentiating the objective function and replacing it in the 2nd Kuhn-Tucker optimality condition (1-9) yields the condition that for all j: hfy [ (pi,- ) + 7 a (f)] 0. Regrouping terms in this sum produces i 9U f ii Pfi Now define rj(f) (EPj + 7 (f))- a i 'afj With this definition, it can be seen that the second Kuhn-Tucker condition is equivalent to fj fjri(f) Y YiPi,j (2-1) i Pfi This form for the optimality condition is useful as it reformulates the equation into a fixed point condition. In 1990, Green developed the One-Step-Late (OSL) algorithm [8] to apply the EM algorithm to the Penalized Maximum Likelihood problem. The One-Step-Late name comes from the choice to evaluate the penalty at the current iterate rather than the (unknown) update, as it allowed the EM algorithm to be applied to the PML problem. In Green's OSL, the update to iterate f" is given by f~+ = f-rj (ff) 7 ) iPij i i Using the fact that 1f(f) E pj y U (f), the update can be expressed in the alternative form: f"+l fn + v(fn) (2-2) where aE vj(f) --fjrhrj (f). (2-3) Thus if Green's OSL algorithm converges, it will satisfy the same fixed point condition as the second Kuhn-Tucker condition. Note that if 7 = 0, Green's OSL algorithm reduces to the EMML Algorithm (1-5). However, Green's OSL algorithm does not guarantee positivity of the iterates or convergence. 2.1.1 Lange's Improvement to Green's OSL In 1990, Lange showed that altering Green's algorithm to include a line search could enforce convergence and positivity [10]. In Lange's algorithm, rj(f") is required to be be positive for all j. To meet this condition, the penalty parameter 7 be must be chosen small enough to keep rj(fn) > 0 for all possible iterates f" and voxels j. In the method, Green's OSL update is altered to include a step-size parameter s (taking s=1 reduces the method to the original OSL method). Thus the update becomes f"+l f" + s"v(f") (2-4) where v is defined as before. The step-size s" is chosen to be a positive value that gives both positivity of the new iterate (taking s" < 1 enforces this condition) and a decrease in the objective function E. To decrease the objective function, define P(s) = E(f + sv(f")) (2-5) If P'(0) < 0, the algorithm will be decreasing for small values of the step-size. With this condition that rj(f") > 0, it follows that P'(O) =- j frj (f")(2)2 (f ) < 0 if fL (ff") / 0 for some j. Since P'(0) < 0, the objective function will be decreasing for small values of the step-size s. A line-search is then performed for an appropriate step-size that will give a positive update and a decrease in the objective function. In Lange's algorithm, the step-size that gives the greatest decrease in the objective function is denoted s and a search is performed to find a value s < s within c of s. This condition requires an exact line search, eliminating the possibility of faster inexact search algorithms. 2.1.2 Mair-Gilland Method One way to remove the restriction of positivity of the r-term is to allow the constant s to be a negative value. In this method (developed for use in the RM Algorithm [14]), the step-size s" is allowed to be negative. The sign of s" depends on the sign of P'(0). If P'(0) < 0, the objective function will be decreasing for small positive values of the step-size. In this case, the step-size is chosen as in Lange's method. A search is made for a value that ensures positivity of the iterate, but is also within c of the value that gives the greatest decrease. If P'(0) > 0, then E will be increasing for small positive values of s. In this case, s" is allowed to be negative. The search is made in the negative direction, again for a value that ensures positivity and is within c of the value that gives the greatest decrease in the objective function. Although this method will decrease the objective function and convergence, the limit may not satisfy the second Kuhn-Tucker optimality condition. 2.2 Accelerated Penalized Maximum Likelihood Method (APML) The accelerated penalized maximum likelihood method (APML) was developed in 2004 by ('!ih ig et al [5]. In the APML method, surrogate functions to decrease the objective function and provide convergence. In the method, a penalty function 0(r) and associated neighborhood system is chosen as described in Chapter 1. In addition, the related function 7(r) = r(r)/r is required to be non-increasing. In the method, the trial iterate is formed by f'1 -Gj(f) + /Gj(f)2 8Ej(f)F()/4F ) (2 6) with the defined functions given by YiPi,j fj Z j Pij f4, F, 2/3 ,7(f-, f) kENj G p 2 i(f fk)(f + fk) i kENj As this was seen to converge slowly, an acceleration factor (the APML Method) was added to the algorithm [4]. As in Lange's method, a step-size factor is added to the algorithm. In APML, the accelerated update is given by f"+i f" + sv(f) (2-7) with v(f") = f" f". In Lange's method, a search method is employ, l1 to find a suitable choice for the step-size. Here, surrogate functions are used to find the best possible step-size needed to decrease the objective function. Surrogate functions were first applied to PET reconstruction by DePierro [6]. In order to minimize a function f over a set S using the method of surrogate functions, another function i must be found satisfying f(x") = (x") f(x) > Q(x) for all x e S. The difficulty in using this method lies in finding an appropriate function Q. In APML, creating the surrogate function requires the calculation of four separate functions related to the penalty function, possibly increasing the computation 14 time. Reconstructions gathered using the APML algorithm will be compared to reconstructions using the revised algorithms presented in this paper. CHAPTER 3 GENERAL PML ALGORITHM 3.1 General Algorithm Statement In the maximum likelihood method, a search is made for an emission intensity f that minimizes the objective function E(f) = L(f) + 7U(f) under the restriction that fj > 0 for all j. In this chapter, an iterative algorithm to find this minimum will be presented. Under specific conditions, convergence to a minimizer will be proven. 3.1.1 General Algorithm Update The new algorithm is formed by taking Lange's modification to OSL, but treating the step-size as a vector instead of a constant. To apply the algorithm, initialize first by taking fo = Yi yi/J. Given a current iterate f", the next iterate, f"+l can be constructed by the following update: fn+l= f" + sn v(f") (3-1) where s" is the step-size vector rj(f) (i pi,j + 7l-(f)) 1 vj(f) fjrj(f) (f) (x 0 y)j = xjyj (the Hadamard product). Note that v(f) and r(f) are the same defined functions as in the OSL algorithm. Given a current iterate f", the step-size vector s" = [s,, ..., s"] is chosen to produce the next step in the algorithm. This vector must satisfy the conditions: 0 < |sj| < K for some large constant K (3-2) QE l < ( rj(fn) E(fJ)) 1 (33) E(f + s'+ v(f)) < E(fn)if lv(fn)l / 0 (3-4) (s' r(f'))j > 0. (3-5) The first three conditions on the step-size vector enforce respectively boundedness of the step-size, positivity of the iterate and decrease in the objective function. The final condition requires that for each voxel j, the sign of sj matches the sign of rj(f). Note that a unique s is not guaranteed. The revised algorithm presented after the proof provides an example of one way to choose a suitable step-size vector. Once a penalty function and method of constructing a step-size vector (both matching the required criteria) have been chosen, the algorithm can proceed. Continue iterating until a suitable stopping criteria has been reached. Convergence of this algorithm will now be investigated. 3.1.2 Penalty function criteria The general penalty term U(f) is of the form U(f) E iy(f, fy) (3-6) i j where 0(r) is required to satisfy the following properties: S(0) 0 limroo (r) -+ o (-r) = (r) (0 is an even function) c C1 (0 has a differentiable first derivative) b is a strictly convex function (0"(r) > 0). The first three conditions are basic conditions for a penalty function, while the latter two conditions are required for convergence. Functions that match these conditions include the logcosh and quadratic penalties. Note that the general Lalush-Tsui penalty does not meet these conditions due its lack of convexity. 3.2 Proof of Convergence Theorem 1. Let < be a p ,'.l'//l; function .-,/->;bi;/ the conditions stated in Section (3.1.2). Let E be the PML objective function E(f) = L(f) + 7U(f). DA i, fo = E yi/J for all j (the il1.,>rithm is initialized with a flat emission ';. iilu:). Let the ilj,.rithm update be given by (3-1) with the step-size vector chosen such that conditions (3-2 35) are -.,i/r.f ,. The il'ji .':thm will then converge to a unique minimizer. The proof (extended from Lange's PML proof [10]) will demonstrate that the algorithm converges to a unique vector that satisfies both Kuhn-Tucker optimality conditions, hence is the minimizer of E. It consists of two parts. In the first, Zangwill's convergence theorem is applied to prove that the objective function converges and that all limit points will fulfil the second Kuhn-Tucker condition. In the second part of the proof, it will be shown that there is a single limit point and this point will also fulfil the first Kuhn-Tucker condition. Definition 1. A stat.':'r,:, r; point is 1. I;,,, l to be ,i ':', i' ,'Ilq function f where for allj, f VE(f) 0. By this definition, the 2nd Kuhn-Tucker condition will be met at a stationary point. Zangwill's Convergence Theorem [11] is now applied. Theorem 2. Z.i: ;,.jI :'s Convergence Theorem. Let T be a point-to-set function /, /7;,., on a metric space X. Generate a sequence {x"} such that x"'+ E T(x"). Let a solution set, S C X, be given il. ,:i with a continuous function E. Suppose that 1. the sequence {x"} is contained in a compact subset of X 2. if x i S, E(y) < E(x) for all y c T(x) 3. if x e S, E(y) < E(x) for all y e T(x) 4. the i'nij'l:',j T is closed at points outside of S. Then the limit of i,,;, convergent subsequence lies in S and {E(x")} converges monotor,::. allii to Es = E(y) for some y E S. To apply Zangwill, let X = R", S be the set of stationary points and let E the PML objective function. For any emission intensity f, define the set Sf to be the set of all possible step-size vectors s fulfilling condition (3-2 3-5). Then define the point-to-set function T by the equation: T(f) {f+s v(f) : s Sf}. (3-7) Note that if f e S, v(f) = 0, thus T(f) {f}. In the PML algorithm, the next iterate is created by taking f"+l = f" + s 0 v(f") for a particular choice of s E Sfn, hence f"+l e T(f'). Lemma 1. The point-set i,,. j'.,:I' T is closed at all nonstationary points 0. The closure definition requires that if 0" 0 S, ', E T(0O), and, , then E T(O). Hence it must be shown that there exists a step-size vector s E So such that =- 0 + s ( v(O) holds for all 0 S. By the definition of the update, for each n there is a step-size vector s" E Son such that 0 + s" O v(08). Since 2- is assumed to be continuous, rj(f) is continuous and hence vj(f) will also be continuous. Let W = {j : vj(O) / 0}. Since 0 is not a stationary point, W is nonempty. For each j E W, there exists a constant Mj > 0 such that for all mj > My, vj(0") / 0. Therefore for all j e W and m > maxj(Mj), m j 0 i Vj (0m) which converges to bj 0j Sj(0) For j e W, define sj to be ^(. Therefore for all j W, j = + sjvj(0). If j ( W, vj(0) = 0 and hence j = Oj. Thus for all j, j = Oj + sjvj(0) where sj = 1 for j W. By the continuity of E, v, 0 and b, s E So and hence b E T(O). Lemma 2. The iterates belong to a bounded and compact set. It will be shown that if IIf"l -- oc, E(f") -+ o as well. The log-likelihood term L is a function of the form Y ci yilog(ci), where ci = Pfi and yi are both positive terms. If IIf"ll oo, then for some i, f7 oo as well. For each voxel j, there exists at least one tube i with pij / 0. Hence if the intensity at voxel j becomes unbounded, c = Pfi -+ o forcing ci yilog(ci) -- o0. Thus if the iterates become unbounded, the log-likelihood function will also become unbounded. Hence if IIf"ll o0, E(f") = L(f) + 7U(f") o0. By construction, E(f") < E(fO) for all n. Thus, it cannot be true that E(f") -- oo. Therefore, the iterates must belong to a bounded set. By construction, E(f"+l) < E(f"). Thus E(f") < E(f0) for all n. Therefore the iterates belong to the set {OIE(O) < E(fo)}. Due to the the continuity of E, the set is closed. Since the set is closed and bounded, it is a compact set. Theorem 3. The limit of i,:1, convergent subsequence is a stationary point. Zangwill's theorem is now applied. Let X = R" and S =the set of stationary points. If f" e S, v(f') = 0 and hence T(f") = {f"}. If f" P S, by construction, E(f" + s" ( v(f")) < E(f") for any step-size vector s" satisfying conditions (3-2 3-5). By Lemma 1, the point-to-set mapping T is closed at non-stationary points. Thus Zangwill's theorem applies and all limit points are stationary points. In addition, E(f") Es = E(y) for some y E S. It must be shown that Es is well-defined. If x E S, then x is a limit point of the sequence. Hence there is a subsequence {ff } converging to x. Since E(f") -+ E(y), it must be true that E(x) = E(y). Lemma 3. If E(f) = Es for some n, then f" e S. In the proof of Zangwill's conditions, E(f-+1) < E(f') if f i S and E(f-+1) < E(f') if f E S. Thus if f" S, E(f'+1) < E(f") = Es. But E(f") converges monotonically to Es, giving a contradiction. Thus f" e S. Lemma 4. The set of limit points is finite. The results of Zangwill's Theorem state that all limit points of the sequence are stationary points. A stationary point has been defined to be one where for each index j, f 2(f) = 0. With the requirement that U is a strictly convex (penalty) term, E will be a strictly convex function [10]. Therefore, given any subset W of [1, J], E will remain a strictly convex function on the plane {fj 0 j E W}. Hence there is a unique point on this plane that is a minimizer for E. Denote this point fw. Thus for any subset W, there is only one point fw with -(fw) = 0 for j ( W and fj = 0 for j E W. Since there are a finite number of possible subsets W, the set of stationary points is finite. Lemma 5. The distance between successive iterates tends to zero. I||f+l fnl| II||n Qs v(ff)||. Since Is'| < K for all n, j, it suffices to prove that vj(f') goes to zero. The function v is continuous and is zero for any stationary point. Since the set of limit points is finite, given e > 0, a single 6 > 0 can be found such that if If" yll < 6 for any y E S, then lv(f")l = Ilv(f") v(y) | < Let Z be the set of all sequence points within 6 of a limit point and the W be the remaining sequence points. If W is infinite, it contains a convergent subsequence (the iterates are contained in a compact set). Thus there are points in W arbitrarily close to a limit point, a contradiction. Hence W must be finite. Let N be the largest sequence index in W. Thus if n > N, I f" y|| <6 for some y E S. Hence Ilv(f") I < c. Therefore v(f") 0. Lemma 6. The set of limit points consists of a single station point. By Lemma 5, the distance between iterates goes to zero. Ostrowski proved that the set of limit points of such a sequence is both compact and connected [15]. Since this set is finite, it must consist of a single point. Denote this intensity by f. Theorem 4. The iterates converge to the minimizer of E. There is a single limit point for the sequence, hence the sequence will converge to f. Assume this function is not the minimizer of E. Since the limit is a stationary point, the second Kuhn-Tucker criterion fj (f) = 0 is satisfied, so the first optimality condition must be violated. Therefore, an index j can be found such that fj = 0 and f(f) < 0 and hence there is an N such that (f") < 0 for all n > N. By construction, rj(fn) and sj(fn) have the same sign and fj > 0. Thus for n > N, fj+ fj + s vj(f) ff sj fj r(f)- (ff) > fl and fj does not converge to zero, contrary to the assumption. Therefore, both Kuhn-Tucker conditions hold at the convergence point and hence it is the minimizer of E. If the penalty function lacks convexity (as in the generalized Lalush-Tsui penalty (1-7), the proof of Theorem 4 will not apply as it requires convexity. Thus, as with the modified RM algorithm discussed in the following chapter, the only result that can be proved is that all of the limit points generated from the sequence will be minimizers. 3.3 New PML Algorithm The general statement does not guarantee the existence of a method of choosing the step-size vector with the needed conditions. However, there do exist several v-wi of meeting the necessary conditions. In Lange's algorithm [10], all entries in the step-size vector are the same positive value. In this method, conditions (3-2 3-5) are satisfied by assuming that rj(f) is positive for all emission intensities f and voxels j. With this restriction, the objective function will be decreasing for small positive values of the step-size s. However, in order to ensure the positivity of r(f), the strength of the penalty term 7 may need to be decreased, which alters the objective function needed to be minimized. The revised algorithm presented below has been designed to remove this strong condition on the penalty term by using a different method to choose the step-size, while still providing for the convergence of the iterates. 3.3.1 Non-Uniform Step-Size Method In this algorithm, the step-size vector consists of two values, one positive and one negative. To find these values at the nth step in the iteration, form the two sets, D+(f") := {jlrj(fj) > 0} and D_(f") := {jrj(fj) < 0}. Since Spi,j + 724 (fn) will be finite for each j, r(f") can never be zero. Therefore, D+(f") and D_(f") will partition the set of voxels. These sets will determine which value of the step-size will be applied to the appropriate voxel. Define the equations g"(tl, t2 { vj(fn) if j D(fP) 8 gilt1, = t2 (3-8) fj + t2 (fn) if j E D_(fF) (t1, 2) = E(g"(t1, t2)). (3-9) The current value of the objective function at step n is given by b(0, 0) = E(f"). Since the goal is to decrease the value of E, a search is made for a new update in the direction of the negative gradient. This will give the direction of the steepest descent of E. If a single value is used for our constant (as in Lange's method [10] or the Mair-Gilland method [14]), it corresponds to using (t,t) as the search direction in (3-8). This may not be the direction of the greatest decrease. In order to achieve the greatest decrease in the objective function, the vector (tl, 2) is assumed to be in the direction of the negative gradient. Therefore, define (t, 2) = -sVQ(0, 0) for some s > 0 to be chosen later. Calculating the two partial derivatives yields: S(0,0) E0 (f>)v(f") E frj (fn)( )2(fn) at, aF aft9 j6D" j6D" (0,0) (f")v (f") f rj(fn)( )2(fn). jCD" jCD" Therefore by defining f E (f) fjrj(f)( -)2(f) (3 10) jED (f) '(f) E jrj(f)(j)2(f) jED (f) it can be seen that (t1,t2) s(r(f"+), T-r(f)). Since (T, r_) provides a search direction only, it can be normalized at this stage to produce a unit vector (7-(f), 7_(f)) (7r and 7_ will now refer to the normalized vector). The new iterate is thus formed by: f+l = f" + Tr(f") v(f") (3 11) where Tj(f) = +(f) if j E D+(f) and Tj(f) r(f) if j c D_(f). The problem has thus been reduced from a 2-dimensional search for (ti, t2) to a 1-dimensional search for s. Defining P(s) = E(f + sT(f") e v(f)), yields P'(0) = VE(f" + s(ff) O v(f")) d(f" + sr(f") ) v(f")) I o ds E -(T)2(fn)j(f1)frj(f')_< 0. 8E The final inequality comes from the fact that rj has the same sign as rj (and so (s 7jrj > 0 for all j, proving the fourth criteria for the step-size vector). Thus for small positive values of the step-size coefficient s, E will be decreasing. It remains to be shown that the method of selecting the step-size coefficient s will result in a step-size vector, st, fulfilling the conditions stated in (3-2 3-5). 3.3.2 Line Search Methods In the non-uniform step-size method discussed above, the step-size vector consists of sT. Now, a search must be made for an appropriate step-size coefficient to ensure the decrease in the objective function and positivity of the iterate. First, a bounding value, smax must be found to ensure positivity. Using the definition of the update, f+' = fj + sn v(f)vj(f) aE sf- l(f/)fi(f)l---f(ff) f/( s"(f')rj(f") (f )). 9fj Therefore, by defining 8E s = minj(|rj(fn)Tj(fn)- (f )|)- (3 12) then fI 0.If (3-t2) it can be seen that if s" < s~,, then n+l > 0. If K < sax, let sa, = K. Now a search is made for a step-size coefficient within the bounding value that decreases the algorithm. If P(s) is defined to be E(fn+l(s)) = E(f+sr(f") v(f")), then the algorithm will decrease the objective function if the constant s" satisfies P'(s") < 0. Using the definition of P(s), df"+l P'(s) = VE(f+(s)) d ds OE El 0E (f"+ (s))rj(ff) fr(ff") (ff). In order to find a constant s' to produce a negative value, a search algorithm is emploi-. 1 One possible choice is the bisection rule, which produces a result within a specified precision. The algorithm proceeded with this choice, but with slow computation speeds as it requires multiple calculations of P'(s). In order to decrease the iteration time, the Armijo rule [1] was emploi-. 1 The Armijo rule is an inexact search algorithm, here used to approximate where P'(s) = 0. The value returned by the Armijo rule will be a constant s that will decrease the objective function as required. To apply the Armijo rule, the search begins with an initial guess (here taken to be the minimum of sma and K). The initial value is decremented this by a constant factor 3 (taken to be 1/3) until the following relationship holds (6 is a positive constant chosen to be close to 0): P(0) + s6P'(0) > P(s). The following figure illustrates one application of the Armijo rule and the value that would be returned by the method. The use of the Armijo rule produced much decrement acceptable value start point x P(s) Figure 3-1: Armijo rule smaller computation times than the bisection method. The former requires the calculation of the derivative only at the origin and thereafter calculates values of P(s) to give a rough approximation to the solution. Even though it does not give the best value for the step-size coefficient, it was seen that the number of iterations needed to reach a stopping point did not significantly increase when the Armijo Rule was applied. Hence, the savings in computation time caused it to be used in the reconstructions. Since the step-size coefficient s' satisfies P'(s") < 0, the objective function decreases as needed. The update is then formed by fn+1 fn + szn_(fr) o v(fr) and the iteration proceeds. It remains to be shown that this choice of step-size vector will satisfy the required conditions for convergence. By construction, (sTj(f)) < K. By specifying that (sTj(f)) < smax, (3-3) is met. (sr) decreases the objective function, satisfying (3-4). Finally, by construction, Tj(f) and rj(f) have the same sign, proving (3-5). Since conditions (3-2 3-5) have been satisfied, the iteration will converge as proven in Theorem 1. 3.3.3 Stopping Criteria The algorithm should terminate when the Kuhn-Tucker optimality conditions (1-9) are realized. The projected gradient function 1 can be used to accomplish this. This function is defined by aE pgd(f) = max(fj (f)),0) fJ o. (3-13) 9fj Theorem 5. pgd(f) = 0 if and only if the Kuhn-Tucker oi'nl,h,/,:li, conditions are met. Assume the Kuhn-Tucker conditions are not met. If the first condition is not true, then (f) < 0, resulting in pgd(f) > 0. If the second condition is not met, fj and j-(f) are nonzero for some j. If fj > | (f) pgd(f) > |-(f)| > 0. If the other inequality holds, then pgd(f) > fj > 0. Thus if the Kuhn-Tucker conditions do not hold, the projected gradient function is nonzero. Now assume the Kuhn-Tucker conditions are met. For each j, (f) > 0 and at least one of fj and E(f) are zero. If fj > 1-(f), then E(f) 0 and |max(fj E(f),0) fj f) 0. If E(f) > f,, then fj 0 and max(fj (f),0) 0. -(f aO) fjf fj O. Given some choice of small c > 0, the algorithm can be terminated when pgd(f) < e. 1 as I-'-. -1 I by Prof. Hager CHAPTER 4 MULTIPLE IMAGE-MOTION RECONSTRUCTION 4.1 Algorithm Desription 4.1.1 RM Algorithm The Penalized Maximum Likelihood method discussed in Chapter 3 can be extended for use in the reconstruction of time-d, 1 -1, images. In this problem, a series of CT-generated images of an object are reconstructed along with any motion of the object between the scans. The RM algorithm, as developed by Gilland and Mair [13] [2] [14], can be used to reconstruct this series. However, its convergence to a minimizer has not been proven. A modification of their work allows for convergence of an altered objective function. As in their work, the algorithm will be described initially in the continuous case, then discretized to produce the iteration. In this reconstruction, a total of K time-lapsed scans (or frames) are analyzed. In the continuous case, the frame number will be denoted in the subscript, fk(z). In the discrete case, the frame number will be the first subscript, fk,j referring to the jth voxel in frame k. Along with the unknown emission intensities in each frame k, fk(w), the motion vector, mk(z) from frame k to frame k + 1 is also unknown. This motion vector describes the movement of each point z in frame k, as z moves to z + mk(z) in frame k + 1. In the continuous case, the negative log-likelihood of the reconstruction is given by K L(f) [Pfk(i) -k(i)logPf(i)] (4 1) k=0 i where Pfk(i) is the continuous projection operator given by Pfk(i) pi(z)fk(z)dz and pi(z) is the probability that an emission at w will be registered in detector i (compare with the discrete projection operator given by Pfi = iPijfj). The penalty term U now measures how well the motion term matches up points between frames. Ideally, fk(z) = fk+l(z + mk(z)) but it is to be expected that this will not occur with real data. Hence the image-matching term between frames k and k + 1 is given by Uk(fk, fk+, mk) fk() fk+(z + mk(z))]2dz (4-2) and the penalty term becomes the sum of these individual terms. K U(f, m) Uk(fk, fk+1,m) (4 3) k= 1 In addition, the objective function includes a third term Es = 1= Esk modeling the elastic strain of the material. Esk(mk) = 1/2 A(z)(uk,x +' + k, )2dz + (z) + ,z)dz +1/2 (z)[(uk,y + Vk,)2 + (Uk,x + Wk,x)2 + (Vk,z+ "' ,)2 d where A and p are material parameters,uk, 1, i, w are the components of mk and Uk,x is the partial derivative of Uk with respect to x. Thus after adding constants to balance the contribution of each term, the continuous objective function becomes E(f, m) a cL(f) + U(f, m) + 3Es(m). (4-4) The reconstruction problem thus becomes a minimization of E(f, m) with the restriction that the intensities in each frame remain positive. At each iteration of the RM Algorithm (the outer loop), two processes are run in succession the R-Step and the M-Step. In the R-Step (an inner loop), the previous motion iterate is used to improve the reconstruction of the frames. Since Es does not depend on f, the R-Step is applied only to aL(f) + U(f, m) with the motion vector held constant. A PML reconstruction method can be applied (such as the non-uniform step-size method described C'!i pter 3) to decrease this term and produce the next iterate. In the M-Step (an inner loop), the motion vector is updated using the frames generated from the R-Step. In this step, the conjugate gradient algorithm is applied to decrease U(f, m) + 3Es(m) with f now held constant. 4.1.2 Image-Matching Penalty In order to run the algorithm, the image space is partitioned into J equally sized voxels. However, after discretizing the image space, it cannot be assumed that the center of a voxel is mapped to the center of another voxel by the motion vector. Thus, it is necessary to distribute the motion over adjoining voxels. To accomplish this, let (uk,i, Vk,i, Wk,i) be the components of the motion vector mapping voxel i in frame k into frame k + 1. Now define di, dj and dk to be the decimal portions of u, v and w respectively. If (ai, bi, ci) are the coordinates of the center of voxel i, then the motion vector will map this voxel to the "voxel" in frame k + 1 centered at (ai + Uk,i, bi + vk,i, Ci + Wk,i). The target voxel to be approximated will lie in at most eight real voxels. Define dooo to be the voxel located at ([ai + Uk,i], [bi + Vk,i], [ci + i,' ]). Then the target voxel will lie in a cube of dimension 2-voxels with dooo in the lower left corner. The emission intensity in the target voxel can then be found by a trilinear approximation using the voxels forming the cube, with the weights proportional to the distance from center of each voxel to the center of the target voxel. Thus, the discrete version of the image-matching term (4-3) is given by: Ukcfk, fk+l,mk) fk,i- Ck,i,jfk+l,j)2 (4-5) i jENk,i where Nk,i is the 8-voxel cube formed by mk,i and the weights, k,ij, are given by Table 4-1. This form for the image-matching term is currently implemented in the voxel voxel center weight dooo ([ai + uk,i], [b + vk,i], [ci + ) (1 )( dj) dk) dloo ([ai + Uk,i] + 1, [bi + Vk,i], [Ci + ) ( dj)( dk) dolo ([ai + Uk,i], [bi + Vk,i] + 1, [ci + (1 di) d( dk) dllo ([ai + Uk,i] + 1, [bi + vk,i] + 1, [ci + ,, ) didj( dk) dool ([ai + Uk,i], [bi + Vk,i], [ci + wk,i] + 1) (1- d)(- dj)dk dlol ([ai + Uk,i] + 1, [bi + vk,i], [Ci + Wk,i] + 1) d dj)dk doll ([ai + Uk,i], [bi + k,i] + 1, [ci + wk,i] + 1) ( di)djdk d111 ([ai + Uk,i] + 1, [bi + vk,i] + 1, [ci + wk,i] + 1) didjdk Table 4-1: RM voxel neighborhood objective function used in the RM Algorithm. However, it is difficult to perform the conjugate gradient (I -Step) on this term. it has proven difficult to implement the M-Step on the this term since the neighborhood system Nk,i depends on the motion vector. To solve this problem, the M-Step has instead been performed on a linearization of the image-matching term [14] given by Qk(m) fk(w)- fk+(W)- (mk(w) m(w)). Vfk+(w + imk(w)dw with the previous motion estimate given by m. However, it has not been proven that a decrease in the linearized term Qk will result in a decrease in the original image-matching term Uk as required. One possible way to solve this problem is to perform the linearization at the objective function level rather than at each inner loop. Hence, the original objective function is modified to replace Uk with a linearized version Uk. With this change, the algorithm will now be referred to as the modified RM algorithm. Since the objective functions are different, a minimizer for the modified RM algorithm is not expected to be a minimizer for the original RM algorithm. However, a decrease in the modified objective function can be proven, allowing for convergence results. To describe the modified objective function, define the term K U = YUk(fkfk+1,mk) k= 1 Uk (fk()- fk+1(z)- mk(z)) Vf-k+l(Z))2dz (4-6) and update the objective function to become E(f, m) = aL(f) + U(f, m) + 3Es(m). (4-7) In addition, linearizing the image-matching term gains convexity in the objective function (the RM image-matching term is not convex in m). However, this linearization is valid only for small displacements. If the motion between the frames is large, the reconstruction may not be valid. Using this modified objective function, the two inner loop sub-processes can now be described in detail. 4.1.3 Modified RM Iteration The R-Step consists of minimizing aL(f) + U(f, m) over the set of positive intensities f. As this in the same form as the objective function in ('!i Ilter 3, the non-uniform step-size algorithm (3-1) can be utilized. Thus, the continuous R-Step update for each frame k is given by f+ = f + sv(f") (4-8) f s frk(fn)DkE(f, m) with rk(f)(z) (= EiPi(z)+DkU(f,m))-'. As derived in the Appendix of [14], the derivatives DkE (taken with respect to fk) are given by DkE(f, m) DkU(f, m) DkUk(fk, fk+I, mk) Dk k- I(fk-1, fk, mk-1) a (pi~z -k (ipi ()z+ Dki(fZ, M) N y ((z) Pf(i)" ) + DnU(f, m) Pfk(i) Dkk (fk, fk+, mk) + DkUk- (f-1, fk, k m1) 2(fk fk+l mk Vfk+l) -2(fkl fk mki Vfk) +2(Vfk- 1 V7 V(mk Vfk)) mk-_ +2(f-_1 f i m_ Vfk)(V mkl). The last derivative can be seen (in the 2-frame case without loss of generality) by using the formula d D2 U (fl, f2, ml) E(fI(z), f2(z) + th(z), mi(z))It=-o dt (4-9) to obtain D2U (fl, f2, mi) -2 (fl(z) f2(z) mi(z) Vf2(z))h(z)dz -2 (fi(z) f2(z) mi(z) Vf2(z))(mi(z) Vh(w))dz. D2U(fl f, m) = if equation (4-9) can be written as f tb(z)h(z)dz. As the first term is already in this form, only the second needs to be simplified. This can be done by applying the identity Sgm Vh(z) S(V gm)h(z) S(Vg m)h(z) + (gV. m)h(z) to yield the equality for the second term: 2 J((V(fi(z) f2(z) mi(z) Vf2(z)) mi)h(z)dz +2 (fi(z) f2(z) mi(z) Vf2(z))(V mi)h(z)dz. Distributing the V in the first integral and adding the results to the first term yields the desired derivative. After discretizing the voxels (and approximating derivatives with central differences), (4-8) reduces to the discrete iteration: fJ = ft + s",jk,j (4-10) where fkjrk,(f) (f) and rk,(f) (- Ca k,i,j +DkU(f, m))-1. As before, the step-size s" vector must be chosen so that (3-2 3-5) are satisfied. The method described in C'! Ilpter 3 can be used to find such a vector. Once the emission intensity has been updated by the R-Step, the M-Step is run to update the motion vector. This step consists of minimizing U(f, m) + 3Es(m) over all possible motion vectors m. The conjugate gradient algorithm (using the Polak-Ribiere variant) is then applied to the discretized voxels to update the motion vector [14]. In practice, at each iteration of the RM algorithm (the outer loop), one iteration each of the R-Step and M-Step occurs (the inner loop). Continue running the outer loop under a stopping criteria has been reached, such as the projected gradient (3-13) falling under a chosen critical value. 4.2 Convergence Results Although the new objective function E aL + U + 3Es is convex in f, Lemma 4 does not apply in this situation (unless a single motion vector is fixed). Hence, convergence to a single minimizer is not guaranteed. Theorem 6. If {ff} is a sequence of iterates produced by (4-10), then each limit point of the sequence is a minimizer for the objective function E and E(f") converges. As before, Zangwill's Convergence Theorem is applied. The definition of the point-to-set function T is expanded to the following: Ti(f,m) {(g, mcG)|g E T(f)} (4-11) where T is the previous point-to-set function defined in (3-7) (with a range of possible step-size vectors) and mce is the motion vector outputted from applying the conjugate gradient algorithm to U(g, m) + 3Es(m). Let E be the linearized objective function and define the set of stationary points S to be the set of all points (f, m) with v(f) = 0. If f e S, then T(f) = {f} and the R-Step will not change the value of L + yU. The M-Step (which does not affect the reconstructed image f) has been defined to decrease yU + Es Thus E(Ti(f, m)) < E(f, m) if f e S. If f S, then the R-Step will decrease L + yU as desired. Following this, the M-Step will decrease the value of 7U + Es. Hence E(TI(f, m)) < E(f, m) if f S. The closure condition is proved by noting that T1 can be split into its two parts, thus allowing the critical closure condition on f to be proved as in Lemma 1. Thus by Zangwill, all limit points of the sequence {(f", m"n)} lie in S and E(f", m") will monotonically converge to E(y, m) for some (y, m) e S. Now to prove all limit points are minimizers, let y be a limit point for the algorithm. Consider the subsequence f" where for each j, II f" y|| < 1/j (the definition of a limit point gives the existence of such points). Now as in Theorem 4, assume y does not meet the first K-T condition for a minimizer (By Zangwill, y e S, hence the second K-T condition minimizer of E is met). Therefore, an index k can be found such that yk = 0 and OJ(y) < 0 and hence there is an J such that E(f j) < 0 for all j > J. By construction, rk(fnd) and Sk have the same sign and fj > 0. Thus for j > J, f' = f + "' (f) f fsf3kk f fkk9 35 and fk, does not converge to yk = 0, contrary to the assumption. Therefore, both K-T conditions hold at the limit point. If the modified M-Step is terminated after iteration K (locking in a specific value for the motion vector) and the modified R-Step is used exclusively afterwards, the convergence proof for penalized maximum likelihood presented in Chapter 3 will suffice to give convergence for the modified RM algorithm. CHAPTER 5 RECONSTRUCTIONS 5.1 PML Reconstructions To test the non-uniform step-size algorithm presented in C'!h pter 3, a series of reconstructions was performed on a phantom source image 1 A = [A, ..., Aj], of dimension J 128x128. As di-',11l -I below, the phantom image is that of a typical cardiac scan, with lungs, heart, spine and other tissues. Higher intensity values correspond to regions of greater metabolic activity (in particular the heart and lung regions). The detector data, Y = [y, ..., yi], was created by calculating the 3 Figure 5-1: Phantom chest source image probability matrix P using (1-2). The chest emission was then projected into the detector space using yi = Z pijAj. Finally, random Poisson noise was added to the 1 created by Dr.Anderson data. The PML objective function is formed by the sum E(f) L(f) + U(f). The penalty term, U(f) was formed using the logcosh function with 6 taken to be 50. U(f) E log cosh f (5 1) i jENi The corresponding neighborhood system, N, is formed by the 8 closest voxels, with the diagonal elements having a weight of and the 4 axial neighbors a weight of 1, as seen in Figure 1-3. The strength of the penalty term, 7, ranged over a set of values to provide a comparison. The reconstructions were formed by the non-uniform step-size method (3-1). In separate trials, the search for the step-size coefficient, s, was performed with the Armijo Rule and the bisection method. For each value of 7, the algorithm was run until the projected gradient (3-13) fell under 10-2. The resulting reconstructions were then compared to the source image with a root-mean-square calculation Zj(fj- Aj)2 128 128 to determine optimal values to use for the penalty strength. As a comparison, the same phantom scanner data was used in the APML algorithm (2-7). The reconstructions outputted from APML were visually identical to the reconstructions gained using the non-uniform step-size algorithm. In all trials, the objective function decreased rapidly, then leveled off. It was seen that the image quality did not improve significantly after this leveling off of the objective function. As can be seen in Figure 5-2, initially the objective function decreases faster with the non-uniform step method than with APML, but eventually both fall to the same level. The computation times for APML were approximately the same as the computation time using the non-uniform S- APML Non Unif L- -3 5 Iteration Figure 5-2: Comparing the decrease in the objective function step-size method with the Armijo Rule. The APML algorithm was able to reduce the projected gradient factor in less iterations than with the non-uniform step-size method. Table 5-1 gives both the RMS error and the iterations required for algorithm termination, occurring when the projected gradient fell below 10-2. In Figure 5-3, 7 Method Iterations CPU time RMS error 0.015 Bisect 219 475 49.5 0.020 Bisect 227 492 46.5 0.025 Bisect 235 631 45.5 0.030 Bisect 183 530 45.4 0.015 Armijo 220 118 49.6 0.020 Armijo 228 125 46.4 0.025 Armijo 268 129 45.5 0.030 Armijo 229 100 45.4 Table 5-1: PML reconstruction data the RMS error is plotted against the penalty paramter 7. Reconstruction using the non-uniform step-size algorithm (applying either bisection or Armijo in the line search produced an identical RMS error) produced a smaller RMS error than reconstruction of the image than using the APML reconstruction method. Thus, these methods were able to produce a reconstruction closer to the source phantom image. Due to its inexact search, the Armijo rule reduced the projected gradient 50------------------ .. APML 49 Non Unit S485 0 0 48 S475 S47 465 46 45 5 451 3 4 5 Gamma Figure 5-3: Comparing the root-mean-square (RMS) error criteria in a significantly shorter cputime than the bisection method (although the iterations were approximately the same). The RMS plot in Figure 5-3 shows that a penalty parameter taken from [0.03,0.04] resulted in the closest reconstruction to the original phantom. The reconstructions seen in Figure 5-4 were formed by applying the bisection method. Although the iteration time was longer than using the Armijo method, the bisection method gave a visually better image with higher values of the penalty parameter, 7. The next set of reconstructions found in Figure 5-5 were formed by applying the Armijo Rule. In this case, as 7 increased, the Armijo rule did not perform as well in high-intensity regions. This can be observed in the heart region of the later reconstructions. Figure 5-6 di- l1' a plot formed by a vertical line thru the heart region. In this region, the exact bisection rule was able to produce a smoother image in the high-intensity heart region than the non-exact Armijo rule. Note that both methods produced similar intensities except in the heart region. With further iterations, it was seen that the non-smooth intensity in the heart region was reduced (using a smaller value for the projected gradient cutoff would accomplish this). A tradeoff has thus been made between reconstruction speed (the Armijo rule) and accuracy in the step-size choice (bisection). One possible solution (a) 7 = 0.020 (b) 7 = 0.025 (c) 7 = 0.030 (d) 7 = 0.035 Figure 5-4: Phantom image reconstruction using bisection rule (a) 7 0.020 (c) 7 = 0.030 (d) 7 = 0.035 Figure 5-5: Phantom image reconstruction using Armijo rule -Armijo \, Bisect Figure 5-6: Line plot of heart region (b) 7 = 0.025 is to combine these methods begin with the Armijo method to approximate the solution, then use bisection for precision. Originally, the new algorithm was designed to remove the restrictions on the penalty function imposed by Lange's algorithm (2-4). Figure 5-7 shows two reconstructions of the chest, the first using the non-uniform step-size method and the second using Lange's algorithm (and the Armijo Rule). In this case, a quadratic penalty was used in place of the logcosh function. With a penalty parameter of 7 = 0.03, the r(f) function becomes negative (a situation not covered in Lange's algorithm). As can be seen from the reconstructions, the non-uniform step-size method is able to reconstruct the image correctly, but Lange's algorithm (together with the Armijo Rule) produces artifacts. (a) Non-uniform step-size method (b) Lange's method Figure 5-7: Comparing reconstruction algorithms 5.2 Modified RM Reconstructions The modified RM algorithm presented in C'! lpter 4 was tested on a 2-frame NCAT phantom image of a heart. The algorithm reconstructs both images with the assistance of the motion vector between the frames. The objective function used in the modified RM algorithm is given by E(f, m) = aL(f) + U(f, m) + 3Es(m) with the constants taken to be a = 0.02 and 3 = 0.005 as in previous reconstructions. The linearized penalty term U was used in the reconstruction instead of the virtual voxel method used in the RM algorithm of Mair and Gilland [14]. Although the use of the virtual voxels was seen to produce useful reconstructions, the convergence of the RM objective function was unproven (but was observed numerically). In addition, the use of the linearized penalty produced a decreasing objective function as seen in Figure 5-8. The second plot in Figure 5-8 shows the functions involved in the M-Step, U and sEs. As can be seen, both U and 3Es are 15 2200 2000 u 0 1800 10 E > != C 3 1600 0 LL. 1400 - 1200 1000 0 0 10 20 30 40 50 60 0 10 20 30 40 50 60 Iteration Iteration (a) Modified RM objective function (b) U and 3Es Figure 5-8: Modified RM function plots increasing for each outer loop. Despite this, the M-Step is reducing their sum for each M-Step as can be seen in the sample RM data found in Table 5-2. Typically, the R-step increases the value of U while reducing aL. Then the M-Step decreases U and increases 3Es to give a decrease in the objective function. The reconstructions formed using the modified RM algorithm (with the linearized penalty, U), were comparable to those outputted using the original RM algorithm (using the approximated voxels penalty method (4-1)). Figure 5-9 di-pli-v' a quiver plot of the reconstructed Slice 20 combined with the motion vector. Figure 5-10 di-pl i-1 a comparison of the methods. The top reconstructions were Iteration 1 R-Step 1 M-Step 2 R-Step 2 M-Step 3 R-Step 3 M-Step 4 R-Step 4 M-Step 5 R-Step 5 M-Step aL U 2866.4 1.326 1.262 2728.4 1.678 1.673 2569.1 2.341 2.232 2374.0 3.367 3.315 2116.4 5.638 5.330 3Es Objective 2867.76 0.0203 2867.72 2730.11 0.0228 2730.11 2571.48 0.0750 2571.42 2377.44 0.104 2377.42 2122.11 0.277 2121.97 Table 5-2: RM data a 0.02, 3 = 0.005 :..M 10 20 30- 40 - 50- 60- * U 10 20 30 40 50 60 Figure 5-9: Quiver plot of Slice 20 performed with the modified RM algorithm. The lower plots were reconstructed using the original RM algorithm, both using the parameters a 0.02 and 3 = 0.005. Although the original RM method gives a better estimate of the motion penalty, (a) Frame 1 Slice 20 Modified RM (c) Frame 1 Slice 20 Original RM Figure 5-10: RM reconstructions (b) Frame 2 Slice 20 Modified RM (d) Frame 2 Slice 20 Original RM a decrease in the objective function (and hence convergence) remains unproven. The use of the linearized penalty in the modified RM algorithm produces a similar reconstruction, while in addition giving the needed decrease in the objective function. 5.3 Conclusions Several methods exist to solve the problem of image recovery in Emission Tomography. Green's OSL algorithm as extended by Lange provides for convergence in a reduced set of penalty functions (which do not include the image-matching penalty used in the RM algorithm). Although Lange's algorithm can be employ, 1 in larger settings, neither its convergence to a minimizer nor a decrease in the objective function is guaranteed. If the step-size is allowed to become negative [2], a decrease in the objective function and algorithm convergence is once again guaranteed, but this minimizer may not satisfy the first Kuhn-Tucker optimality condition. However, a further extrapolation to a step-size vector will both guarantee convergence to a minimizer and allow for a larger set of penalty functions. Unlike previous methods, the convergence proof of the generalized method allows for the use of inexact faster line searches, such as the Armijo Rule. The non-uniform step-size algorithm provides one way of constructing a step-size vector that meets the conditions required for convergence. This generalized method can be extended to other image recovery problems, such as multiple image-motion recovery. In this case, the generalized algorithm is incorporated into the modified RM algorithm for reconstruction. Future work in this area will concentrate on refining the algorithm to increase speed without losing image quality. APPENDIX SAMPLE ALGORITHM FORTRAN CODE This Fortran code can be utilized to run the non-uniform step method as described in ('!Ch pter 3. Code to calculate the penalty term (and its derivative) is required. The probability matrix pi,j is stored in a sparce matrix format consisting of hrow, hcol and hv preserving the non-zero entries of p. The iteration loop consists of the following code: do j=l,nhdim temp = hcol(j) hsum(temp) = hsum(temp) + hv(j) enddo do j=l,nhdim tempi = hrow(j) temp2 = hcol(j) AX(templ) = AX(templ) + hv(j)*f(temp2) enddo do j=l,nhdim tempi = hrow(j) temp2 = hcol(j) if(AX(templ).ne.0)then AY(temp2) = AY(temp2) + hv(j)*Y(templ)/AX(templ) endif enddo call findderivpenalty(f,dU) do j=l,NVOXX dE(j) = hsum(j) AY(j) + gamma*dU(j) enddo tau(1) = 0 tau(2) = 0 do j=l,NVOXX r(j) = sA(j) + gamma*dU(j) fl(j) = f(j)*AY(j)/r(j) v(j) = fl(j) f(j) enddo tempi = 0 do j=l,NVOXX if(r(j).le.0)templ = 1 enddo if (templ.eq. )then do j=l,NVOXX if(r(j) .ge.0.)then tau(l)=tau(l)+(f(j)*dE(j)**2)/r(j) else tau(2)=tau(2)+(f(j)*dE(j)**2)/r(j) endif enddo else tau(1) = 1 endif tau(3) = sqrt(tau(1)**2 + tau(2)**2) tau(1) = tau(1) / tau(3) tau(2) = tau(2) / tau(3) do j=l,NVOXX if (r(j).ge.0.)then v(j) = v(j)*tau(1) else v(j) = v(j)*tau(2) endif enddo smax=0 do j=l,NVOXX if(r(j) .ge.0)then smax = max(smax, tau(1)*dE(j)/r(j)) else smax = max(smax, tau(2)*dE(j)/r(j)) endif enddo if (smax.eq.0.)then smax = 1 else smax = 1/smax endif sopt =find-s(f,smax,v,hcol,hrow, hv,Y,dE,gamma) do j=l,NVOXX f(j) = f(j) + sopt*v(j) enddo Utilizing the Armijo rule in the line-search can be accomplished by the code (called by find-s in the main loop): dEdotv = 0 do j=l,NVOXX dEdotv = dEdotv + dE(j)*v(j) enddo tol = 0.1 k= 0 rho = 3 30 continue alpha = smax / (rho**k) do j=l,NVOXX f1(j) = f(j) + alpha*v(j) enddo do j=l,nhdim temple = hrow(j) temp2 = hcol(j) AX(templ) = AX(templ) + hv(j)*fl(temp2) enddo call findpenalty(fl,U) E= 0 do j=l,ndet if(AX(j).ne.0)then E = E + AX(j) Y(j)*log(AX(j)) endif enddo 51 do j=l,NVOXX E = E + gamma*U(j) enddo flag = EO + alpha tol dEdotv k = k+1 if(E.gt.flag)goto 30 finds = alpha REFERENCES [1] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty, 1993. Nonlinear r'ii'i"'.r .: ':r theory and l..', .:thms. New York : J. Wiley & Sons, 2nd edition. [2] Z. Cao, D. Gilland, B. Mair, and R. Jaszczak. Three-dimensional motion estimation with image reconstruction for gated cadiac ECT. IEEE Trans. Nuclear Science, pages 384-388. [3] R. B. Carroll. An or'/,. '.,, ,';l series approach to positron emission '.i,,,.'/,,/',; PhD thesis, University of Florida, Gainesville, FL, October 1998. [4] J. C'I ,i:;, J.M.M. Anderson, and B.Mair. An accelerated penalized maximum likelihood algorithm for positron emission tomography. IEEE Trans. Image Processing, forthcoming, 2006. [5] J. C'!I in,; J.M.M. Anderson, and J. Votaw. Regularized image reconstruction algorithms for positron emission tomography. IEEE Trans. Med. Imag., 23(9):1165-1175, 2004. [6] A.R. DePierro. A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography. IEEE Trans. Med. Imag., 14(1):132-137, 1995. [7] S. Geman and D. McClure. B -i ,'i image analysis: an application to single photon emission tomography. Proc. American Statistical S .. I:, pages 12-18, 1985. [8] P. Green. On use of the EM algorithm for penalized likelihood estimation. J.R. Statist. Soc, 52(3):443-452, 1990. [9] D. Lalush and B. Tsui. A generalized Gibbs prior for maximum a posterioro reconstruction in SPECT. Phys. Med. Biol., 38:729-741, 1993. [10] K. Lange. Convergence of EM image reconstruction algorithms with Gibbs smoothing. IEEE Trans. Med. Imag., 9(4):439-446, 1990. [11] D. Luenberger, 1973. Introduction to linear and nonlinear "'"-','i'rn.':', Reading : Addison-Wesley. [12] B. Mair. A mathematical model incorporating the effects of detector width in 2D PET. Inverse Problems, 16:223 224, 2000. 53 [13] B. Mair, D. Gilland, and Z. Cao. Simultaneous motion estimation and image reconstruction from gated data. Proc. 2002 IEEE Internat. Symp. Biomed. Imaging, pages 661-664, 2002. [14] B. Mair, D. Gilland, and J.Sun. Estimation of images and nonrigid deformations in gated emission CT. IEEE Trans. Med. Imag., 25(9):1130-1144, 2006. [15] A. M. Ostrowski, 1973. Solutions of equations in euclidean and banach spaces. New York : Academic Press, 3rd edition. [16] L. A. Shepp and Y. Vardi. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imag., 1:113-122, 1982. BIOGRAPHICAL SKETCH Jeff Zahnen was born on August 7, 1974, in Chicago, Illinois. He received his Bachelor of Science degree in mathematics with a minor in history from the University of Florida in 1996. Following this, he received his master's degree in 1998. In 2000, he began to study medical imaging techniques under the supervision of Dr. Bernard Mair. In 2006, he completed his Ph.D. work on emission tomography methods. During this time, he also served as the coach of the University of Florida College Bowl Team. |

Full Text |

PAGE 4 Forhispatienceandadvicethroughtheyears,Ioermygreatestthankstomyadvisor,Dr.BernardMair.Inaddition,IwouldliketothankDr.DavidGilland,Dr.WilliamHager,Dr.JamesKeesling,Dr.ScottMcCulloughandDr.DavidWilsonfortakingthetimetoserveonmycommitteeandprovidingsuggestionstoimprovethework.Finally,Iwouldliketothankmyparentsfortheirunwaveringsupportofmystudies. iv PAGE 5 page ACKNOWLEDGMENTS ............................. iv LISTOFTABLES ................................. vii LISTOFFIGURES ................................ viii ABSTRACT .................................... ix CHAPTER 1MATHEMATICSOFPOSITRONEMISSIONTOMOGRAPHY(PET) 1 1.1IntroductionToPET .......................... 1 1.2MathematicalModelforPET ..................... 2 1.3PETProbabilityFunction ....................... 3 1.4ReconstructionTechniques ....................... 4 1.4.1MaximumLikelihood ...................... 4 1.4.2PenalizedMaximumLikelihood(PML) ............ 6 1.5MinimizingFunction .......................... 9 2PMLRECONSTRUCTIONMETHODS .................. 10 2.1Green'sOneStepLateAlgorithm ................... 10 2.1.1Lange'sImprovementtoGreen'sOSL ............. 11 2.1.2Mair-GillandMethod ...................... 12 2.2AcceleratedPenalizedMaximumLikelihoodMethod(APML) ... 12 3GENERALPMLALGORITHM ...................... 15 3.1GeneralAlgorithmStatement ..................... 15 3.1.1GeneralAlgorithmUpdate ................... 15 3.1.2Penaltyfunctioncriteria .................... 16 3.2ProofofConvergence .......................... 17 3.3NewPMLAlgorithm .......................... 21 3.3.1Non-UniformStep-SizeMethod ................ 22 3.3.2LineSearchMethods ...................... 24 3.3.3StoppingCriteria ........................ 26 4MULTIPLEIMAGE-MOTIONRECONSTRUCTION .......... 27 4.1AlgorithmDesription .......................... 27 v PAGE 6 .......................... 27 4.1.2Image-MatchingPenalty .................... 29 4.1.3ModiedRMIteration ..................... 31 4.2ConvergenceResults .......................... 33 5RECONSTRUCTIONS ............................ 36 5.1PMLReconstructions .......................... 36 5.2ModiedRMReconstructions ..................... 42 5.3Conclusions ............................... 45 APPENDIX:SAMPLEALGORITHMFORTRANCODE ........... 47 REFERENCES ................................... 53 BIOGRAPHICALSKETCH ............................ 54 vi PAGE 7 Table page 4{1RMvoxelneighborhood ........................... 30 5{1PMLreconstructiondata ........................... 38 5{2RMdata=0:02,=0:005 ........................ 44 vii PAGE 8 Figure page 1{1A2-DPETtuberegisteringanemission .................. 1 1{2Tuberegions ................................. 3 1{3Sampleneighborhoodsystemwithassociatedweights ........... 7 1{4Derivativesofthepenaltyterms ....................... 8 3{1Armijorule .................................. 25 5{1Phantomchestsourceimage ......................... 36 5{2Comparingthedecreaseintheobjectivefunction ............. 38 5{3Comparingtheroot-mean-square(RMS)error ............... 39 5{4Phantomimagereconstructionusingbisectionrule ............ 40 5{5PhantomimagereconstructionusingArmijorule ............. 41 5{6Lineplotofheartregion ........................... 41 5{7Comparingreconstructionalgorithms .................... 42 5{8ModiedRMfunctionplots ......................... 43 5{9QuiverplotofSlice20 ............................ 44 5{10RMreconstructions .............................. 45 viii PAGE 9 Positronemissiontomography(orPET)hasbeenutilizedforyearsinmedicaldiagnosis.ThePETreconstructionproblemconsistsofrecoveringtheimagefromthedatagatheredbythedetectorsinthemachine.Penalizedmaximumlikelihoodtechniquesallowfortherecoveryoftheimagewithasmoothingtermtoimprovetheimage.Previousmethodsusingthistechniqueforcedtheuseofaweakpenaltyterminordertogainconvergence.Usinganewmethod,thisrestrictionisremovedwhileretainingconvergence.Inaddition,themethodisextendedtotheproblemofrecoveringmultiplePETimagesalongwithavectordescribingthemotionbetweentheimages. ix PAGE 10 1{1 displaysatubeformedbydetectorsand.Anear-simultaneousarrivaloftwophotonsatoppositeendsof Figure1{1: A2-DPETtuberegisteringanemission 1 PAGE 11 thetubeisthenregisteredasasingledetection.ThetubeshowninFigure 1{1 willregisterthedisplayedemissionatpointBbutnottheemissionatpointA. Insinglephotonemissioncomputedtomography(SPECT),asinglephotonisemitted.Inthismethod(morewidelyappliedthanPETduetotherelativecostsofthescanningmachines),thedetectorsarelocatedonlyononesideofthepatientandasingledetectionresultsinacount.Thisscanisperformedfordierentanglesbetweenthepatientanddetectortoprovideafullimage. Inbothprocedures,theaimistodeterminethenumberofemissionsfromthevariouslocationsinthepatient.Theresultingintensityimagesareusedtomakemedicaldiagnoses.ThisworkwillconsidertheproblemofPETimagereconstructionusingthemaximumlikelihood(ML)methodwiththeadditionofapenaltyterm.AlthoughdevelopedforPETreconstruction,themethodcanbeextendedforSPECTreconstruction.Unlikepreviousmethods,therevisedalgorithmpresentedusesamoregeneralpenaltyterm,allowingforawiderrangeofsettings.Proofofconvergenceandreconstructionresultsareobtainedusingthismethod. 16 ]introducedamathematicalmodelforPETin1982.Inthemodel,therandomemissionsfromvoxeljaremodeledbyaPoissonprocesswithmeanfj.Thusthe PAGE 12 detectorcountintubei,Yi,willalsobegovernedbyaPoissonprocesswithmean wherepi;jisthegeometricprobabilitythatanemissioninvoxeljisregisteredbydetectortubei.UsingthepropertiesofthePoissondistribution,theprobabilitythatthereareyicountsregisteredindetectoriisgivenbyeiyii 1{2 .Thevoxels Figure1{2: Tuberegions areassumedtobesosmallthattheprobabilitycanbeassumedapproximatelyconstantontheentirevoxel.Ifzisthecenterofvoxelj,thenpi;jisdenedbythe PAGE 13 probabilitythatthepathofapairofphotonsemittedfromthecenterofthevoxel,z,intersectsbothdetectorsdeningthetube.Thisprobabilityisproportionaltotheangle-of-viewofthetubefromz.Thus,ifC^zDisdenedtobetheacuteangleformedbythelinespassingthroughCandzandDandz,thefollowingdenitionresults[ 12 ]. 16 ].In1998,Carroll[ 3 ]decomposedthedata,probabilitymatrixandintensityfunctionsintobasisfunctions,transformingtheproblemintoaninnitesetoflinearequationstosolve.Ouralgorithmthoughisbasedonathirdmethod,Maximumlikelihoodestimation(ML). PAGE 14 Thedetectorcountsareindependentrandomvariables,sotheseprobabilitiesaremultipliedtogethertoyieldtheprobabilityofobtainingthedatafromanemissionintensitymapoff.AnMLestimateisanintensityfthatmaximizesthefollowingprobability.ieiyii Inaddition,theintensitymap,fisconstrainedtobenon-negative(ffjfj0g). AnequivalentproblemistoapplythenegativelogarithmandperformaminimizationofthetermXi(PfiyilogPfilogyi!): TheExpectation-Maximization(orEM)algorithmdevelopedbySheppandVardi[ 16 ]canthenappliedtondtheminimizerofthisterm,resultingintheEMMLalgorithm. wherefnreferstothenthstepintheiteration. Alternately,theEMMLiterationcanbeexpressedasadescentalgorithm.fn+1j=fnjfnj@L @fj(fn)=Xipi;j PAGE 15 Restatedinmatrixformyieldstheequationfn+1=fnD(fn)rL(fn): 8 ]. 7 ],isemployedbyrstchoosingapenaltyfunctionandthenformingtheterm where!isaweightingfactorthatassociatespixels.Thiscanbeusedtolinkneighboringpixels,forcingthemtohavesimilarintensities.Forexample,Ni,canbedenedoneachvoxelitobethesetofpixelsinaneighborhoodlatticeofthatvoxel.Thustheweight!i;j=0ifj=2N(i).Ifj2N(i),theweightisdenedtobeanonzeroconstant.OnepossibilityistoletNiconsistoftheeightneighboringvoxelsasshowninFigure 1{3 withassociatedweights.Inaddition,aconstantvalue,ismultipliedtothepenaltytermtocontroltheinuenceofthepenaltyinthereconstruction.As!0,thereconstructionsapproachthoseoftheEMalgorithm.However,ifischosentobetoolarge,theimageswillbeoversmoothed. PAGE 16 Figure1{3: Sampleneighborhoodsystemwithassociatedweights Thereexistmanychoicesforthepenaltyfunction(r).Afrequentlyusedchoiceisapenaltyoftheform(r)=logcoshr [ 8 ],whereisapreviouslychosenconstantthatallowsforsomene-tuningofthepenalty.Thisconstantdetermineshowthepriorisscaledwithrespecttor,asthelogcoshpenaltywillsmoothalldierencesaboveacertaincutovalue,determinedby.As!1,themethodwillapproachresultsoftheEMAlgorithm(asnosmoothingwilloccur).Thispenaltyisusefulasitapproximatesquadraticsmoothingforsmallvaluesofr,butitsderivativeremainsbounded. In1993,LalushandTsui[ 9 ]developedanewpenaltyfunctionthatallowedforevengreatercontroloverthepenalizedrangeofintensities.NotingthatcurrentPMLmethodsdidnotrequirethepenaltyfunctiontobecalculated,theyinsteadbuiltthefunctionfrompropertiesofitsderivative.Thus,theyproduced @r=r r)2]1=2(1{7) whereandareconstantsthataectthespreadandlocationrespectivelyofpeakofthederivative.Asincreases,thepeakbecomesmorenarrowandas PAGE 17 increases,thepeakawayfromtheorigin.Theplotofthederivativeshowsthemajorityofthesmoothingoccurringintheregionunderthepeakofthederivative[ 9 ].Thusalteringthevaluesofandwillaecttherangeofsmoothingintheimage.Bycomparingthisfunctiontothederivativeofthelogcoshfunction,0(r)=(1=)tanhr itcanbeseenthatwhilethelogcoshfunctiondoesnothaveapeakedderivative (b)Derivativeofthelogcoshpenalty Derivativesofthepenaltyterms (itincreasestoanasymptoticlimit),itcanbeconsideredasaninnitely-widepeak.Thus,Green'slogcoshpenaltycanbeapproximatedasaLalush-Tsuipenaltybytakingclosetozero(producinganearinnitely-widepeak).IntegratingtheLalush-Tsuipenaltytermresultsintheirpenaltyfunction:(r)= 2p PAGE 18 InMaximumLikelihood,thisobjectivefunctionmustbeminimizedoverthesetofpositiveintensityvaluesffjfj0g.ThisminimummustsatisfytheKuhn-Tuckeroptimalityconditionsfortheproblem.Inthissituation,thosetwoconditionsarethatforallj @fj(f)0(1{9)fj@E @fj(f)=0 ItwillbeshownthattherevisedalgorithmpresentedinChapter3willmeettheseconditionsatconvergence. PAGE 19 Dierentiatingtheobjectivefunctionandreplacingitinthe2ndKuhn-Tuckeroptimalitycondition( 1{9 )yieldstheconditionthatforallj:fj[Xi(pi;jyipi;j @fj(f)]=0: @fj(f)]=fjXiyipi;j @fj(f))1: Thisformfortheoptimalityconditionisusefulasitreformulatestheequationintoaxedpointcondition. In1990,GreendevelopedtheOne-Step-Late(OSL)algorithm[ 8 ]toapplytheEMalgorithmtothePenalizedMaximumLikelihoodproblem.TheOne-Step-Latenamecomesfromthechoicetoevaluatethepenaltyatthecurrentiteraterather 10 PAGE 20 thanthe(unknown)update,asitallowedtheEMalgorithmtobeappliedtothePMLproblem.InGreen'sOSL,theupdatetoiteratefnisgivenbyfn+1j=fnjrj(fn)Xiyipi;j @fj(f)=Pipi;jyipi;j @fj(f),theupdatecanbeexpressedinthealternativeform: where @fj(f):(2{3) ThusifGreen'sOSLalgorithmconverges,itwillsatisfythesamexedpointconditionasthesecondKuhn-Tuckercondition.Notethatif=0,Green'sOSLalgorithmreducestotheEMMLAlgorithm( 1{5 ).However,Green'sOSLalgorithmdoesnotguaranteepositivityoftheiteratesorconvergence. 10 ].InLange'salgorithm,rj(fn)isrequiredtobebepositiveforallj.Tomeetthiscondition,thepenaltyparameterbemustbechosensmallenoughtokeeprj(fn)>0forallpossibleiteratesfnandvoxelsj.Inthemethod,Green'sOSLupdateisalteredtoincludeastep-sizeparameters(takings=1reducesthemethodtotheoriginalOSLmethod).Thustheupdatebecomes wherevisdenedasbefore.Thestep-sizesnischosentobeapositivevaluethatgivesbothpositivityofthenewiterate(takingsn<1enforcesthiscondition)andadecreaseintheobjectivefunctionE.Todecreasetheobjectivefunction,dene PAGE 21 IfP0(0)<0,thealgorithmwillbedecreasingforsmallvaluesofthestep-size.Withthisconditionthatrj(fn)>0,itfollowsthatP0(0)=Pjfnjrj(fn)(@E @fj)2(fn)<0iffnj@E @fj(fn)6=0forsomej.SinceP0(0)<0,theobjectivefunctionwillbedecreasingforsmallvaluesofthestep-sizes.Aline-searchisthenperformedforanappropriatestep-sizethatwillgiveapositiveupdateandadecreaseintheobjectivefunction.InLange'salgorithm,thestep-sizethatgivesthegreatestdecreaseintheobjectivefunctionisdenotedsandasearchisperformedtondavalues PAGE 22 ~fn+1j=Gj(fn)+q withthedenedfunctionsgivenbyEj=yipi;jfnj Asthiswasseentoconvergeslowly,anaccelerationfactor(theAPMLMethod)wasaddedtothealgorithm[ 4 ].AsinLange'smethod,astep-sizefactorisaddedtothealgorithm.InAPML,theacceleratedupdateisgivenby withv(fn)=~fn+1fn.InLange'smethod,asearchmethodisemployedtondasuitablechoiceforthestep-size.Here,surrogatefunctionsareusedtondthebestpossiblestep-sizeneededtodecreasetheobjectivefunction.SurrogatefunctionswererstappliedtoPETreconstructionbyDePierro[ 6 ].InordertominimizeafunctionfoverasetSusingthemethodofsurrogatefunctions,anotherfunctionmustbefoundsatisfying Thedicultyinusingthismethodliesinndinganappropriatefunction.InAPML,creatingthesurrogatefunctionrequiresthecalculationoffourseparatefunctionsrelatedtothepenaltyfunction,possiblyincreasingthecomputation PAGE 23 time.ReconstructionsgatheredusingtheAPMLalgorithmwillbecomparedtoreconstructionsusingtherevisedalgorithmspresentedinthispaper. PAGE 24 undertherestrictionthatfj0forallj.Inthischapter,aniterativealgorithmtondthisminimumwillbepresented.Underspecicconditions,convergencetoaminimizerwillbeproven. where @fj(f))1 @fj(f) Notethatv(f)andr(f)arethesamedenedfunctionsasintheOSLalgorithm. 15 PAGE 25 Givenacurrentiteratefn,thestep-sizevectorsn=[sn1;:::;snJ]ischosentoproducethenextstepinthealgorithm.Thisvectormustsatisfytheconditions: 0 PAGE 26 Therstthreeconditionsarebasicconditionsforapenaltyfunction,whilethelattertwoconditionsarerequiredforconvergence.Functionsthatmatchtheseconditionsincludethelogcoshandquadraticpenalties.NotethatthegeneralLalush-Tsuipenaltydoesnotmeettheseconditionsdueitslackofconvexity. Theorem1. 3.1.2 ).LetEbethePMLobjectivefunctionE(f)=L(f)+U(f).Denef0j=Piyi=Jforallj(thealgorithmisinitializedwithaatemissionintensity).Letthealgorithmupdatebegivenby( 3{1 )withthestep-sizevectorchosensuchthatconditions( 3{2 3{5 )aresatised.Thealgorithmwillthenconvergetoauniqueminimizer. 10 ])willdemonstratethatthealgorithmconvergestoauniquevectorthatsatisesbothKuhn-Tuckeroptimalityconditions,henceistheminimizerofE.Itconsistsoftwoparts.Intherst,Zangwill'sconvergencetheoremisappliedtoprovethattheobjectivefunctionconvergesandthatalllimitpointswillfullthesecondKuhn-Tuckercondition.Inthesecondpartoftheproof,itwillbeshownthatthereisasinglelimitpointandthispointwillalsofulltherstKuhn-Tuckercondition. 11 ]isnowapplied. 1. thesequencefxngiscontainedinacompactsubsetofX 2. ifx=2S,E(y) PAGE 27 ifx2S,E(y)E(x)forally2T(x) themappingTisclosedatpointsoutsideofS. ThenthelimitofanyconvergentsubsequenceliesinSandfE(xn)gconvergesmonotonicallytoES=E(y)forsomey2S. 3{2 3{5 ).Thendenethepoint-to-setfunctionTbytheequation: Notethatiff2S,v(f)=0,thusT(f)=ffg.InthePMLalgorithm,thenextiterateiscreatedbytakingfn+1=fn+sv(fn)foraparticularchoiceofs2Sfn,hencefn+12T(fn). @fjisassumedtobecontinuous,rj(f)iscontinuousandhencevj(f)willalsobecontinuous.LetW=fj:vj()6=0g.Sinceisnotastationarypoint,Wisnonempty.Foreachj2W,thereexistsaconstantMj>0suchthatforallmj>Mj,vj(m)6=0.Thereforeforallj2Wandm>maxj(Mj),smj=mjmj PAGE 28 whichconvergestojj 3{2 3{5 ).ByLemma 1 ,thepoint-to-setmappingTisclosedatnon-stationarypoints.ThusZangwill'stheoremappliesandalllimitpointsarestationarypoints.Inaddition,E(fn)!ES=E(y)forsomey2S.ItmustbeshownthatESiswell-dened.Ifx2S,thenxisalimitpointofthesequence.Hencethereis PAGE 29 asubsequenceffnjgconvergingtox.SinceE(fn)!E(y),itmustbetruethatE(x)=E(y). @fj(f)=0.WiththerequirementthatUisastrictlyconvex(penalty)term,Ewillbeastrictlyconvexfunction[ 10 ].Therefore,givenanysubsetWof[1;J],Ewillremainastrictlyconvexfunctionontheplaneffj=0jj2Wg.HencethereisauniquepointonthisplanethatisaminimizerforE.DenotethispointfW.ThusforanysubsetW,thereisonlyonepointfWwith@E @fj(fW)=0forj=2Wandfj=0forj2W.SincethereareanitenumberofpossiblesubsetsW,thesetofstationarypointsisnite. PAGE 30 5 ,thedistancebetweeniteratesgoestozero.Ostrowskiprovedthatthesetoflimitpointsofsuchasequenceisbothcompactandconnected[ 15 ].Sincethissetisnite,itmustconsistofasinglepoint.Denotethisintensitybyf. @fj(f)=0issatised,sotherstoptimalityconditionmustbeviolated.Therefore,anindexjcanbefoundsuchthatfj=0and@E @fj(f)<0andhencethereisanNsuchthat@E @fj(fn)<0foralln>N.Byconstruction,rj(fn)andsj(fn)havethesamesignandfnj>0.Thusforn>N,fn+1j=fnj+snjvj(fn)=fnjsnjfnjrj(fn)@E @fj(fn)>fnj Ifthepenaltyfunctionlacksconvexity(asinthegeneralizedLalush-Tsuipenalty( 1{7 ),theproofofTheorem 4 willnotapplyasitrequiresconvexity.Thus,aswiththemodiedRMalgorithmdiscussedinthefollowingchapter,theonlyresultthatcanbeprovedisthatallofthelimitpointsgeneratedfromthesequencewillbeminimizers. 10 ],allentriesinthestep-sizevectorarethesamepositivevalue.Inthismethod,conditions( 3{2 3{5 )aresatisedbyassumingthatrj(f)ispositiveforall PAGE 31 emissionintensitiesfandvoxelsj.Withthisrestriction,theobjectivefunctionwillbedecreasingforsmallpositivevaluesofthestep-sizes.However,inordertoensurethepositivityofr(f),thestrengthofthepenaltytermmayneedtobedecreased,whichalterstheobjectivefunctionneededtobeminimized.Therevisedalgorithmpresentedbelowhasbeendesignedtoremovethisstrongconditiononthepenaltytermbyusingadierentmethodtochoosethestep-size,whilestillprovidingfortheconvergenceoftheiterates. @fj(fn)willbeniteforeachj,r(fn)canneverbezero.Therefore,D+(fn)andD(fn)willpartitionthesetofvoxels.Thesesetswilldeterminewhichvalueofthestep-sizewillbeappliedtotheappropriatevoxel.Denetheequations Thecurrentvalueoftheobjectivefunctionatstepnisgivenby(0;0)=E(fn).SincethegoalistodecreasethevalueofE,asearchismadeforanewupdateinthedirectionofthenegativegradient.ThiswillgivethedirectionofthesteepestdescentofE.Ifasinglevalueisusedforourconstant(asinLange'smethod[ 10 ]ortheMair-Gillandmethod[ 14 ]),itcorrespondstousing(t,t)asthesearchdirectionin( 3{8 ).Thismaynotbethedirectionofthegreatestdecrease.Inordertoachievethegreatestdecreaseintheobjectivefunction,thevector(t1;t2)isassumedtobeinthedirectionofthenegativegradient.Therefore,dene(t1;t2)=sr(0;0)for PAGE 32 somes0tobechosenlater.Calculatingthetwopartialderivativesyields:@ @t1(0;0)=Xj2Dn+@E @fj(fn)vj(fn)=Xj2Dn+fnjrj(fn)(@E @fj)2(fn)@ @t2(0;0)=Xj2Dn@E @fj(fn)vj(fn)=Xj2Dnfnjrj(fn)(@E @fj)2(fn): @fj)2(f) (3{10) @fj)2(f) itcanbeseenthat(t1;t2)=s(+(fn);(fn)).Since(+;)providesasearchdirectiononly,itcanbenormalizedatthisstagetoproduceaunitvector(+(f);(f))(+andwillnowrefertothenormalizedvector).Thenewiterateisthusformedby: wherej(f)=+(f)ifj2D+(f)andj(f)=(f)ifj2D(f).Theproblemhasthusbeenreducedfroma2-dimensionalsearchfor(t1;t2)toa1-dimensionalsearchfors.DeningP(s)=E(fn+s(fn)v(fn)),yieldsP0(0)=rE(fn+s(fn)v(fn))d ds(fn+s(fn)v(fn))js=0=Xj(@E @fj)2(fn)j(fn)fnjrj(fn)0: 3{2 3{5 ). PAGE 33 @fj(fn)=fnj(1snj(fn)rj(fn)@E @fj(fn)): @fj(fn)j)1;(3{12) itcanbeseenthatifsn PAGE 34 ruleisaninexactsearchalgorithm,hereusedtoapproximatewhereP0(s)=0.ThevaluereturnedbytheArmijorulewillbeaconstantsthatwilldecreasetheobjectivefunctionasrequired.ToapplytheArmijorule,thesearchbeginswithaninitialguess(heretakentobetheminimumofsmaxandK).Theinitialvalueisdecrementedthisbyaconstantfactor(takentobe1/3)untilthefollowingrelationshipholds(isapositiveconstantchosentobecloseto0):P(0)+sP0(0)P(s): Figure3{1: Armijorule smallercomputationtimesthanthebisectionmethod.TheformerrequiresthecalculationofthederivativeonlyattheoriginandthereaftercalculatesvaluesofP(s)togivearoughapproximationtothesolution.Eventhoughitdoesnotgivethebestvalueforthestep-sizecoecient,itwasseenthatthenumberofiterationsneededtoreachastoppingpointdidnotsignicantlyincreasewhentheArmijoRulewasapplied.Hence,thesavingsincomputationtimecausedittobeusedinthereconstructions. Sincethestep-sizecoecientsnsatisesP0(sn)<0,theobjectivefunctiondecreasesasneeded.Theupdateisthenformedbyfn+1=fn+sn(fn)v(fn) PAGE 35 andtheiterationproceeds.Itremainstobeshownthatthischoiceofstep-sizevectorwillsatisfytherequiredconditionsforconvergence.Byconstruction,(sj(f)) PAGE 36 4.1.1RMAlgorithm 13 ][ 2 ][ 14 ],canbeusedtoreconstructthisseries.However,itsconvergencetoaminimizerhasnotbeenproven.Amodicationoftheirworkallowsforconvergenceofanalteredobjectivefunction.Asintheirwork,thealgorithmwillbedescribedinitiallyinthecontinuouscase,thendiscretizedtoproducetheiteration.Inthisreconstruction,atotalofKtime-lapsedscans(orframes)areanalyzed.Inthecontinuouscase,theframenumberwillbedenotedinthesubscript,fk(z).Inthediscretecase,theframenumberwillbetherstsubscript,fk;jreferringtothejthvoxelinframek.Alongwiththeunknownemissionintensitiesineachframek,fk(w),themotionvector,mk(z)fromframektoframek+1isalsounknown.Thismotionvectordescribesthemovementofeachpointzinframek,aszmovestoz+mk(z)inframek+1.Inthecontinuouscase,thenegativelog-likelihoodofthereconstructionisgivenby wherePfk(i)isthecontinuousprojectionoperatorgivenbyPfk(i)=Rpi(z)fk(z)dzandpi(z)istheprobabilitythatanemissionatwwillberegisteredindetectori PAGE 37 (comparewiththediscreteprojectionoperatorgivenbyPfi=Pipi;jfj).ThepenaltytermUnowmeasureshowwellthemotiontermmatchesuppointsbetweenframes.Ideally,fk(z)=fk+1(z+mk(z))butitistobeexpectedthatthiswillnotoccurwithrealdata.Hencetheimage-matchingtermbetweenframeskandk+1isgivenby andthepenaltytermbecomesthesumoftheseindividualterms. Inaddition,theobjectivefunctionincludesathirdtermES=PKk=1ESkmodelingtheelasticstrainofthematerial.ESk(mk)=1=2Z(z)(uk;x+vk;y+wk;z)2dz+Z(z)(u2k;x+v2k;y+w2k;z)dz+1=2Z(z)[(uk;y+vk;x)2+(uk;x+wk;x)2+(vk;z+wk;y)2]dz Thusafteraddingconstantstobalancethecontributionofeachterm,thecontinuousobjectivefunctionbecomes ThereconstructionproblemthusbecomesaminimizationofE(f;m)withtherestrictionthattheintensitiesineachframeremainpositive. AteachiterationoftheRMAlgorithm(theouterloop),twoprocessesareruninsuccession-theR-StepandtheM-Step.IntheR-Step(aninnerloop),the PAGE 38 previousmotioniterateisusedtoimprovethereconstructionoftheframes.SinceESdoesnotdependonf,theR-StepisappliedonlytoL(f)+U(f;m)withthemotionvectorheldconstant.APMLreconstructionmethodcanbeapplied(suchasthenon-uniformstep-sizemethoddescribedChapter3)todecreasethistermandproducethenextiterate.IntheM-Step(aninnerloop),themotionvectorisupdatedusingtheframesgeneratedfromtheR-Step.Inthisstep,theconjugategradientalgorithmisappliedtodecreaseU(f;m)+ES(m)withfnowheldconstant. 4{3 )isgivenby: PAGE 39 whereNk;iisthe8-voxelcubeformedbymk;iandtheweights,ck;i;j,aregivenbyTable 4{1 .Thisformfortheimage-matchingtermiscurrentlyimplementedinthe voxelvoxelcenterweight RMvoxelneighborhood objectivefunctionusedintheRMAlgorithm.However,itisdiculttoperformtheconjugategradient(M-Step)onthisterm.ithasprovendiculttoimplementtheM-SteponthethistermsincetheneighborhoodsystemNk;idependsonthemotionvector.Tosolvethisproblem,theM-Stephasinsteadbeenperformedonalinearizationoftheimage-matchingterm[ 14 ]givenbyQk(m)=Zfk(w)fk+1(w)(mk(w)~mk(w))rfk+1(w+~mk(w)dw PAGE 40 function,denetheterm ~U=KXk=1~Uk(fk;fk+1;mk)~Uk=Z(fk(z)fk+1(z)mk(z))rfk+1(z))2dz andupdatetheobjectivefunctiontobecome ~E(f;m)=L(f)+~U(f;m)+ES(m):(4{7) Inaddition,linearizingtheimage-matchingtermgainsconvexityintheobjectivefunction(theRMimage-matchingtermisnotconvexinm).However,thislinearizationisvalidonlyforsmalldisplacements.Ifthemotionbetweentheframesislarge,thereconstructionmaynotbevalid.Usingthismodiedobjectivefunction,thetwoinnerloopsub-processescannowbedescribedindetail. 3{1 )canbeutilized.Thus,thecontinuousR-Stepupdateforeachframekisgivenby (4{8) =fnksfnkrk(fn)DkE(f;m) PAGE 41 withrk(f)(z)=(Pipi(z)+DkU(f;m))1.AsderivedintheAppendixof[ 14 ],thederivativesDkE(takenwithrespecttofk)aregivenbyDkE(f;m)=Xi(pi(z)yk(i)pi(z) dtE(f1(z);f2(z)+th(z);m1(z))jt=0(4{9) toobtainD2~U1(f1;f2;m1)=2Z(f1(z)f2(z)m1(z)rf2(z))h(z)dz2Z(f1(z)f2(z)m1(z)rf2(z))(m1(z)rh(w))dz:D2~U1(f1;f2;m1)=ifequation( 4{9 )canbewrittenasR(z)h(z)dz.Asthersttermisalreadyinthisform,onlythesecondneedstobesimplied.ThiscanbedonebyapplyingtheidentityZgmrh(z)=Z(rgm)h(z)=Z(rgm)h(z)+(grm)h(z) toyieldtheequalityforthesecondterm:=2Z((r(f1(z)f2(z)m1(z)rf2(z))m1)h(z)dz+2Z(f1(z)f2(z)m1(z)rf2(z))(rm1)h(z)dz: PAGE 42 Distributingtherintherstintegralandaddingtheresultstothersttermyieldsthedesiredderivative. Afterdiscretizingthevoxels(andapproximatingderivativeswithcentraldierences),( 4{8 )reducestothediscreteiteration: wherevk;j=fk;jrk;j(f)@~E @fk;j(f)andrk;j(f)=(Pipk;i;j+Dk~U(f;m))1.Asbefore,thestep-sizesnvectormustbechosensothat( 3{2 3{5 )aresatised.ThemethoddescribedinChapter3canbeusedtondsuchavector. OncetheemissionintensityhasbeenupdatedbytheR-Step,theM-Stepisruntoupdatethemotionvector.Thisstepconsistsofminimizing~U(f;m)+ES(m)overallpossiblemotionvectorsm.Theconjugategradientalgorithm(usingthePolak-Ribierevariant)isthenappliedtothediscretizedvoxelstoupdatethemotionvector[ 14 ]. Inpractice,ateachiterationoftheRMalgorithm(theouterloop),oneiterationeachoftheR-StepandM-Stepoccurs(theinnerloop).Continuerunningtheouterloopunderastoppingcriteriahasbeenreached,suchastheprojectedgradient( 3{13 )fallingunderachosencriticalvalue. 4 doesnotapplyinthissituation(unlessasinglemotionvectorisxed).Hence,convergencetoasingleminimizerisnotguaranteed. 4{10 ),theneachlimitpointofthesequenceisaminimizerfortheobjectivefunction~Eand~E(fn)converges. PAGE 43 Asbefore,Zangwill'sConvergenceTheoremisapplied.Thedenitionofthepoint-to-setfunctionTisexpandedtothefollowing: whereTisthepreviouspoint-to-setfunctiondenedin( 3{7 )(witharangeofpossiblestep-sizevectors)andmCGisthemotionvectoroutputtedfromapplyingtheconjugategradientalgorithmto~U(g;m)+ES(m).Let~EbethelinearizedobjectivefunctionanddenethesetofstationarypointsStobethesetofallpoints(f;m)withv(f)=0.Iff2S,thenT(f)=ffgandtheR-StepwillnotchangethevalueofL+~U.TheM-Step(whichdoesnotaectthereconstructedimagef)hasbeendenedtodecrease~U+ES.Thus~E(T1(f;m))~E(f;m)iff2S.Iff=2S,thentheR-StepwilldecreaseL+~Uasdesired.Followingthis,theM-Stepwilldecreasethevalueof~U+ES.Hence~E(T1(f;m))<~E(f;m)iff=2S.TheclosureconditionisprovedbynotingthatT1canbesplitintoitstwoparts,thusallowingthecriticalclosureconditiononftobeprovedasinLemma 1 ThusbyZangwill,alllimitpointsofthesequencef(fn;mn)glieinSand~E(fn;mn)willmonotonicallyconvergeto~E(y;m)forsome(y;m)2S. Nowtoprovealllimitpointsareminimizers,letybealimitpointforthealgorithm.Considerthesubsequencefnjwhereforeachj,jjfnjyjj<1=j(thedenitionofalimitpointgivestheexistenceofsuchpoints).NowasinTheorem 4 ,assumeydoesnotmeettherstK-Tconditionforaminimizer(ByZangwill,y2S,hencethesecondK-TconditionminimizerofEismet).Therefore,anindexkcanbefoundsuchthatyk=0and@E @fk(y)<0andhencethereisanJsuchthat@E @fj(fnj)<0forallj>J.Byconstruction,rk(fnj)andskhavethesamesignandfnjk>0.Thusforj>J,fnj+1k=fnjk+snjkvk(fnj)=fnjksnjkfnjkrk(fnj)@E @fk(fnj)>fnjk PAGE 44 andfnjkdoesnotconvergetoyk=0,contrarytotheassumption.Therefore,bothK-Tconditionsholdatthelimitpoint. IfthemodiedM-StepisterminatedafteriterationK(lockinginaspecicvalueforthemotionvector)andthemodiedR-Stepisusedexclusivelyafterwards,theconvergenceproofforpenalizedmaximumlikelihoodpresentedinChapter3willsucetogiveconvergenceforthemodiedRMalgorithm. PAGE 45 Figure5{1: Phantomchestsourceimage probabilitymatrixPusing( 1{2 ).Thechestemissionwasthenprojectedintothedetectorspaceusingyi=Pjpi;jj.Finally,randomPoissonnoisewasaddedtothe 36 PAGE 46 data.ThePMLobjectivefunctionisformedbythesumE(f)=L(f)+U(f): Thecorrespondingneighborhoodsystem,Niisformedbythe8closestvoxels,withthediagonalelementshavingaweightof1 1{3 .Thestrengthofthepenaltyterm,,rangedoverasetofvaluestoprovideacomparison.Thereconstructionswereformedbythenon-uniformstep-sizemethod( 3{1 ).Inseparatetrials,thesearchforthestep-sizecoecient,s,wasperformedwiththeArmijoRuleandthebisectionmethod.Foreachvalueof,thealgorithmwasrununtiltheprojectedgradient( 3{13 )fellunder102.Theresultingreconstructionswerethencomparedtothesourceimagewitharoot-mean-squarecalculations Pj(fjj)2 todetermineoptimalvaluestouseforthepenaltystrength. Asacomparison,thesamephantomscannerdatawasusedintheAPMLalgorithm( 2{7 ).ThereconstructionsoutputtedfromAPMLwerevisuallyidenticaltothereconstructionsgainedusingthenon-uniformstep-sizealgorithm.Inalltrials,theobjectivefunctiondecreasedrapidly,thenleveledo.Itwasseenthattheimagequalitydidnotimprovesignicantlyafterthislevelingooftheobjectivefunction.AscanbeseeninFigure 5{2 ,initiallytheobjectivefunctiondecreasesfasterwiththenon-uniformstepmethodthanwithAPML,buteventuallybothfalltothesamelevel.ThecomputationtimesforAPMLwereapproximatelythesameasthecomputationtimeusingthenon-uniform PAGE 47 Figure5{2: Comparingthedecreaseintheobjectivefunction step-sizemethodwiththeArmijoRule.TheAPMLalgorithmwasabletoreducetheprojectedgradientfactorinlessiterationsthanwiththenon-uniformstep-sizemethod. Table 5{1 givesboththeRMSerrorandtheiterationsrequiredforalgorithmtermination,occuringwhentheprojectedgradientfellbelow102.InFigure 5{3 0.015Bisect21947549.50.020Bisect22749246.50.025Bisect23563145.50.030Bisect18353045.40.015Armijo22011849.60.020Armijo22812546.40.025Armijo26812945.50.030Armijo22910045.4 Table5{1: PMLreconstructiondata theRMSerrorisplottedagainstthepenaltyparamter.Reconstructionusingthenon-uniformstep-sizealgorithm(applyingeitherbisectionorArmijointhelinesearchproducedanidenticalRMSerror)producedasmallerRMSerrorthanreconstructionoftheimagethanusingtheAPMLreconstructionmethod.Thus,thesemethodswereabletoproduceareconstructionclosertothesourcephantomimage.Duetoitsinexactsearch,theArmijorulereducedtheprojectedgradient PAGE 48 Figure5{3: Comparingtheroot-mean-square(RMS)error criteriainasignicantlyshortercputimethanthebisectionmethod(althoughtheiterationswereapproximatelythesame).TheRMSplotinFigure 5{3 showsthatapenaltyparametertakenfrom[0.03,0.04]resultedintheclosestreconstructiontotheoriginalphantom. ThereconstructionsseeninFigure 5{4 wereformedbyapplyingthebisectionmethod.AlthoughtheiterationtimewaslongerthanusingtheArmijomethod,thebisectionmethodgaveavisuallybetterimagewithhighervaluesofthepenaltyparameter,.ThenextsetofreconstructionsfoundinFigure 5{5 wereformedbyapplyingtheArmijoRule.Inthiscase,asincreased,theArmijoruledidnotperformaswellinhigh-intensityregions.Thiscanbeobservedintheheartregionofthelaterreconstructions.Figure 5{6 displaysaplotformedbyaverticallinethrutheheartregion.Inthisregion,theexactbisectionrulewasabletoproduceasmootherimageinthehigh-intensityheartregionthanthenon-exactArmijorule.Notethatbothmethodsproducedsimilarintensitiesexceptintheheartregion.Withfurtheriterations,itwasseenthatthenon-smoothintensityintheheartregionwasreduced(usingasmallervaluefortheprojectedgradientcutowouldaccomplishthis).Atradeohasthusbeenmadebetweenreconstructionspeed(theArmijorule)andaccuracyinthestep-sizechoice(bisection).Onepossiblesolution PAGE 49 (b)=0:025 (c)=0:030 (d)=0:035 Phantomimagereconstructionusingbisectionrule PAGE 50 (b)=0:025 (c)=0:030 (d)=0:035 PhantomimagereconstructionusingArmijorule Figure5{6: Lineplotofheartregion PAGE 51 istocombinethesemethods-beginwiththeArmijomethodtoapproximatethesolution,thenusebisectionforprecision. Originally,thenewalgorithmwasdesignedtoremovetherestrictionsonthepenaltyfunctionimposedbyLange'salgorithm( 2{4 ).Figure 5{7 showstworeconstructionsofthechest,therstusingthenon-uniformstep-sizemethodandthesecondusingLange'salgorithm(andtheArmijoRule).Inthiscase,aquadraticpenaltywasusedinplaceofthelogcoshfunction.Withapenaltyparameterof=0:03,ther(f)functionbecomesnegative(asituationnotcoveredinLange'salgorithm).Ascanbeseenfromthereconstructions,thenon-uniformstep-sizemethodisabletoreconstructtheimagecorrectly,butLange'salgorithm(togetherwiththeArmijoRule)producesartifacts. (b)Lange'smethod Comparingreconstructionalgorithms PAGE 52 withtheconstantstakentobe=0:02and=0:005asinpreviousreconstructions.Thelinearizedpenaltyterm~UwasusedinthereconstructioninsteadofthevirtualvoxelmethodusedintheRMalgorithmofMairandGilland[ 14 ].Althoughtheuseofthevirtualvoxelswasseentoproduceusefulreconstructions,theconvergenceoftheRMobjectivefunctionwasunproven(butwasobservednumerically). Inaddition,theuseofthelinearizedpenaltyproducedadecreasingobjectivefunctionasseeninFigure 5{8 .ThesecondplotinFigure 5{8 showsthefunctionsinvolvedintheM-Step,~UandES.Ascanbeseen,both~UandESare (b)~UandES ModiedRMfunctionplots increasingforeachouterloop.Despitethis,theM-StepisreducingtheirsumforeachM-StepascanbeseeninthesampleRMdatafoundinTable 5{2 .Typically,theR-stepincreasesthevalueof~UwhilereducingL.ThentheM-Stepdecreases~UandincreasesEStogiveadecreaseintheobjectivefunction.ThereconstructionsformedusingthemodiedRMalgorithm(withthelinearizedpenalty,~U),werecomparabletothoseoutputtedusingtheoriginalRMalgorithm(usingtheapproximatedvoxelspenaltymethod( 4{1 )).Figure 5{9 displaysaquiverplotofthereconstructedSlice20combinedwiththemotionvector.Figure 5{10 displaysacomparisonofthemethods.Thetopreconstructionswere PAGE 53 IterationL~UESObjective 1-R-Step2866.41.3262867.761-M-Step1.2620.02032867.72 2-R-Step2728.41.6782730.112-M-Step1.6730.02282730.11 3-R-Step2569.12.3412571.483-M-Step2.2320.07502571.42 4-R-Step2374.03.3672377.444-M-Step3.3150.1042377.42 5-R-Step2116.45.6382122.115-M-Step5.3300.2772121.97 Table5{2: RMdata=0:02,=0:005 Figure5{9: QuiverplotofSlice20 PAGE 54 performedwiththemodiedRMalgorithm.ThelowerplotswerereconstructedusingtheoriginalRMalgorithm,bothusingtheparameters0:02and=0:005.AlthoughtheoriginalRMmethodgivesabetterestimateofthemotionpenalty, (b)Frame2Slice20ModiedRM (c)Frame1Slice20OriginalRM (d)Frame2Slice20OriginalRM RMreconstructions adecreaseintheobjectivefunction(andhenceconvergence)remainsunproven.TheuseofthelinearizedpenaltyinthemodiedRMalgorithmproducesasimilarreconstruction,whileinadditiongivingtheneededdecreaseintheobjectivefunction. PAGE 55 penaltyusedintheRMalgorithm).AlthoughLange'salgorithmcanbeemployedinlargersettings,neitheritsconvergencetoaminimizernoradecreaseintheobjectivefunctionisguaranteed.Ifthestep-sizeisallowedtobecomenegative[ 2 ],adecreaseintheobjectivefunctionandalgorithmconvergenceisonceagainguaranteed,butthisminimizermaynotsatisfytherstKuhn-Tuckeroptimalitycondition.However,afurtherextrapolationtoastep-sizevectorwillbothguaranteeconvergencetoaminimizerandallowforalargersetofpenaltyfunctions.Unlikepreviousmethods,theconvergenceproofofthegeneralizedmethodallowsfortheuseofinexactfasterlinesearches,suchastheArmijoRule.Thenon-uniformstep-sizealgorithmprovidesonewayofconstructingastep-sizevectorthatmeetstheconditionsrequiredforconvergence.Thisgeneralizedmethodcanbeextendedtootherimagerecoveryproblems,suchasmultipleimage-motionrecovery.Inthiscase,thegeneralizedalgorithmisincorporatedintothemodiedRMalgorithmforreconstruction.Futureworkinthisareawillconcentrateonreningthealgorithmtoincreasespeedwithoutlosingimagequality. PAGE 56 ThisFortrancodecanbeutilizedtorunthenon-uniformstepmethodasdescribedinChapter3.Codetocalculatethepenaltyterm(anditsderivative)isrequired.Theprobabilitymatrixpi;jisstoredinasparcematrixformatconsistingofhrow,hcolandhvpreservingthenon-zeroentriesofp. Theiterationloopconsistsofthefollowingcode: PAGE 59 UtilizingtheArmijoruleintheline-searchcanbeaccomplishedbythecode(calledbynd-sinthemainloop): PAGE 61 [1] M.S.Bazaraa,H.D.Sherali,andC.M.Shetty,1993.Nonlinearprogramming:theoryandalgorithms.NewYork:J.Wiley&Sons,2ndedition. [2] Z.Cao,D.Gilland,B.Mair,andR.Jaszczak.Three-dimensionalmotionestimationwithimagereconstructionforgatedcadiacECT.IEEETrans.NuclearScience,pages384{388. [3] R.B.Carroll.Anorthogonalseriesapproachtopositronemissiontomography.PhDthesis,UniversityofFlorida,Gainesville,FL,October1998. [4] J.Chang,J.M.M.Anderson,andB.Mair.Anacceleratedpenalizedmaximumlikelihoodalgorithmforpositronemissiontomography.IEEETrans.ImageProcessing,forthcoming,2006. [5] J.Chang,J.M.M.Anderson,andJ.Votaw.Regularizedimagereconstructionalgorithmsforpositronemissiontomography.IEEETrans.Med.Imag.,23(9):1165{1175,2004. [6] A.R.DePierro.Amodiedexpectationmaximizationalgorithmforpenalizedlikelihoodestimationinemissiontomography.IEEETrans.Med.Imag.,14(1):132{137,1995. [7] S.GemanandD.McClure.Bayesianimageanalysis:anapplicationtosinglephotonemissiontomography.Proc.AmericanStatisticalSociety,pages12{18,1985. [8] P.Green.OnuseoftheEMalgorithmforpenalizedlikelihoodestimation.J.R.Statist.Soc,52(3):443{452,1990. [9] D.LalushandB.Tsui.AgeneralizedGibbspriorformaximumaposteriororeconstructioninSPECT.Phys.Med.Biol.,38:729{741,1993. [10] K.Lange.ConvergenceofEMimagereconstructionalgorithmswithGibbssmoothing.IEEETrans.Med.Imag.,9(4):439{446,1990. [11] D.Luenberger,1973.Introductiontolinearandnonlinearprogramming.Reading:Addison-Wesley. [12] B.Mair.Amathematicalmodelincorporatingtheeectsofdetectorwidthin2DPET.InverseProblems,16:223{224,2000. 52 PAGE 62 [13] B.Mair,D.Gilland,andZ.Cao.Simultaneousmotionestimationandimagereconstructionfromgateddata.Proc.2002IEEEInternat.Symp.Biomed.Imaging,pages661{664,2002. [14] B.Mair,D.Gilland,andJ.Sun.EstimationofimagesandnonrigiddeformationsingatedemissionCT.IEEETrans.Med.Imag.,25(9):1130{1144,2006. [15] A.M.Ostrowski,1973.Solutionsofequationsineuclideanandbanachspaces.NewYork:AcademicPress,3rdedition. [16] L.A.SheppandY.Vardi.Maximumlikelihoodreconstructionforemissiontomography.IEEETrans.Med.Imag.,1:113{122,1982. PAGE 63 JeZahnenwasbornonAugust7,1974,inChicago,Illinois.HereceivedhisBachelorofSciencedegreeinmathematicswithaminorinhistoryfromtheUniversityofFloridain1996.Followingthis,hereceivedhismaster'sdegreein1998.In2000,hebegantostudymedicalimagingtechniquesunderthesupervisionofDr.BernardMair.In2006,hecompletedhisPh.D.workonemissiontomographymethods.Duringthistime,healsoservedasthecoachoftheUniversityofFloridaCollegeBowlTeam. 54 |