UFDC Home  myUFDC Home  Help 



Full Text  
SIMULTANEOUS RECONSTRUCTION AND 3D MOTION ESTIMATION FOR GATED MYOCARDIAL EMISSION TOMOGRAPHY By ZIXIONG CAO A DISSERTATION 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 2003 This dissertation is dedicated to my parents and my wife ACKNOWLEDGMENTS It is impossible to adequately thank my longtime advisor and mentor, David Gilland for his tremendous guidance, considerable patience and encouragement needed for me to proceed through the graduate study and complete this dissertation. His rich experience in the medicalimaging area, his thoughtful insight, his brilliant and creative ideas, and his great sense of humor have made my graduate study very rewarding. Bernard Mair deserves many thanks. His lecture on numerical analysis introduced me to this area and led to all of the computation work described in this dissertation. His constant guidance and supervision laid a smooth way for my studies and research work. I also thank Chris Batich and David Hintenlang for their discussions, suggestions, and encouragement during the development of this dissertation. It is a great honor to have them serve on my committee. I would like to thank Mu Chen and Karen Gilland for their helpful suggestion on my numerical observer study. I must express my sincere gratitude particularly to William Ditto for his support. I thank the Biomedical Engineering Department (especially AprilLane Derfinyak and Laura Studstill) for their help, kindness, and patience. I thank Rami, Uday, and PingFang for making the lab such a great place to work. I also thank Shu and Jiangyong for their comments and suggestions. Finally, I thank my dear parents, my sister Xia, and my wife Xueshuang for their endless love, support, encouragement, and understanding in dealing with all of the challenges that I have faced. I owe a great debt to them. TABLE OF CONTENTS Page ACKNOWLEDGMENT S ................. ................. iii........ .... LIST OF TABLES ................ ..............vii .......... .... LI ST OF FIGURE S ................. ................. viii............ AB STRAC T ................ .............. xi CHAPTER 1 INT RODUC TION ................. ...............1............ .... 1.1 Introduction ............ ...... ._ .............. 1... 1.2 Si gnificance ............ ..... ._ ...............3.... 2 BACKGROUND ............ ..... ._ ...............5.... 2. 1 Myocardial Emission Tomography ................. ...............5............ ... 2.2 Proj section and Image Reconstruction ................. ...............6............... 2.2. 1 Radon Transform .................. ........... .......... ... .. ......6 2.2.2 Maximum Likelihood and Image Reconstruction ................ ................9 2.3 Myocardium Motion Estimation ................ ...............11................ 2.4 Image Observer Evaluation ................. ...............13................ 2.4.1 Human Observer Shtdy................ ...............13. 2.4.2 Numerical Observer: Background .................. ............ ..................15 2.4.3 Hotelling Observer and Channelized Hotelling Observer...................16 3 FORMULATION OF SIMULTANEOUS ESTIMATION METHOD..................20 3.1 Objective Functions............... ...............2 3.2 EulerLagrange Equations ............ ......_ .. ...............21 3.3 Optimization Algorithm: RM ................. ........... ............... 26..... 3.3.1 Scheme of Optimization Algorithm............... ...............2 3.3.2 Computation of R Step .............. ...............27.... 3.3.3 Computation of M Step .............. ...............29.... 4 IMPLEMENTATION AND EXPERIMENTAL RESULTS: GEOMETRIC PHANTOM ................. ...............3.. 1.............. 4. 1 Simulated Source Obj ect and Proj section Data ................. ........................3 1 4. 1.1 Steps to Generate Geometric Phantom ................ ................ ...._.31 4. 1.2 Dimensions of Geometric Phantom ................. ......... ................32 4. 1.3 Intensity and Defect ................. ...............33.............. 4. 1.4 W all M otion Simulation .............. ............. ... ............3 4. 1.5 Forward Proj section Simulation and Detector Response ................... .....34 4.2 BiValue Lame Constants Model ................. ...............35............... 4.3 Convergence Properties of RM Algorithm .............. ...............38.... 4.4 Motion Estimation Results ................. ...............39........... ... 4.4. 1 Figures of M erit ................. ........... ...............39..... 4.4.2 Lame Constants Study Results .............. ...............41.... 4.4.3 Global Motion Evaluation .............. ...............42.... 4.4.4 Regional Motion Evaluation............... ...............4 4.5 Image Reconstruction Results ................. ...............45................ 4.5.1 Figures of M erit ................. ...............45.............. 4.5.2 Global Image Evaluation ................ ...............46............... 4.5.3 Defect Contrast Evaluation ................. ...............49............... 4.5.4 Regional Image Evaluation............... ...............4 5 IMPLEMENTATION AND EXPERIMENTAL RESULTS: NCAT PHANTOM.53 5.1 Simulated NCAT Phantom............... ...............53 5.2 Motion Estimation Results ................. ...............54........... ... 5.2.1 Global Motion Evaluation .............. ...............54.... 5.2.2 Regional Motion Evaluation............... ...............5 5.3 Image Reconstruction Results ................. ...............58........... ... 5.3.1 Global Image Evaluation ................ ...............58.............. 5.3.2 Defect Contrast Evaluation ................. ...............59........... .. 5.3.3 Regional Image Evaluation............... ...............6 6 OB SERVER EVALUATION OF IMAGE QUALITY ................. ............... .....64 6.1 M ethods ........... ......... ............ ..........6 6. 1.1 SingleSlice CHO Model ............... ... .............. ....................6 6. 1.2 MultiFrame MultiSlice CHOHO Model ................. .....................66 6.2 Results of Geometric Phantom Study .............. ...............69.... 6.2.1 SingleSlice CHO Results.................. ...... ...........6 6.2.2 MultiFrame MultiSlice CHOHO Results .............. .....................6 6.3 Results of NCAT Phantom Study .............. ...............72.... 6.3.1 SingleSlice CHO Results.................. ...... ...........7 6.3.2 MultiFrame MultiSlice CHOHO Results .............. .....................7 6.4 MultiFrame MultiSlice MultiView CHOHO Study .............. ...................74 7 CONCLUSIONS AND FUTURE WORK .............. ...............79.... APPENDIX NUMERICAL ANALYSIS FOR MOTION ESTIMATION .................. 82 LIST OF REFERENCES ................. ...............92................ BIOGRAPHICAL SKETCH .............. ...............98.... LIST OF TABLES Table pg 21. Conditional probabilities of a defect detection task .............. .....................14 41. Ellipsoid semiaxes lengths for source obj ect ................. ...............33............. 42. Poisson ratio and Lame constants ................. ...............37............... 43. Global motion evaluation results of the geometric phantom ................. ................ .44 44. Regional motion evaluation results of the geometric phantom ........._..... ...............45 51. Global motion evaluation results of the NCAT phantom ................. ............... .....56 52. Regional motion evaluation results of the NCAT phantom ................ ................. .58 53. The contrast and noise level of the NCAT phantom ................. ................ ...._.62 54. Five sets of hyperparameters for RM used for regional image evaluation ...............62 61. SCHO detectability index results of the geometric phantom .............. ...................71 62. MMCHO detectability index results of the geometric phantom .............. .... .........._.71 63. SCHO detectability index results of the NCAT phantom ................. .............. ....74 64. MMCHO detectability index results of the NCAT phantom ................ ................. 74 65. MMMCHO detectability index results of the NCAT phantom ................. ...............78 LIST OF FIGURES Figure pg 21. Anatomic orientation of the heart. ............ .....___ ...............6 22. Principle of data acquisition and geometric considerations for SPECT.............._._.....7 23. Coordinate systems transformation. ........._._. ......_ ....__ ......._.........8 24. Images of MLEM and PSMLEM ................. ...............11............... 25. Motion vector representing the displacement of voxels ................. ............. .......12 26. Klein's iterative motionestimation optimization method............. .. ........._ ....13 31. Scheme of the RM al gorithm. ................ ........... ........ ......... ...........27 41. Source objects of the geometric phantom............... ...............32 42. Simulation steps of the geometric phantom .............. ...............33.... 43. Ideal motion fields of the geometric phantom in 3 shortaxis slices. ................... ......3 5 44. Collapse highresolution proj section to lowresolution proj section ................... ...........36 45. Segmentation for bivalue Lame constants model. ............. .....................3 46. Convergence of the RM algorithm ................ ...............39............... 47. Convergence of R step............... ...............40.. 48. Convergence of M step............... ...............40.. 49. Global motion error vs. fl for several Poisson ratios. ............. .....................4 410. Ideal motion vector fields superimposed on true images. ............. .....................42 411. Estimated motion by M step using true images ................. ......... ................43 412. Individually reconstructed images by PSMLEM and estimated motion by M step.43 413. Estimated motion and reconstructed images by RM...................... ...............4 414. Chosen region to calculate regional motion error for the geometric phantom.........45 415. RMS errors of the RM and PSMLEM images. .......... ................ ................47 416. Hyperparameter effects on the RM images ................. ............... ......... ...48 417. RMS error of the RM images as a function of a and B............_.._. .....................48 418. Defect contrast of the PSMLEM and RM images ................. .......................49 419. The tradeoff between the contrast and noise level of the geometric phantom........5 1 420. Subregions for regional image evaluation for the geometric phantom. ..................51 421. The tradeoff between regional intensity bias and variance of the geometric phantom ..........._._ ....._. ._ ...............52..... 51. Long and short axis views of the NCAT cardiac phantom. ................ ................ .54 52. Ideal motion vector fields superimposed on true images .............. .....................5 53. Estimated motion by M step using true images............... ...............55. 54. Individually reconstructed images by PSMLEM and estimated motion by M step..56 55. Estimated motion and reconstructed images by RM .................... ...............5 56. Global motion errors of three methods............... ...............57 57. Chosen region to calculate regional motion error for the NCAT phantom ................58 58. RMS errors of the RM and PSMLEM images. ................ ................. ..........59 59. RMS error vs. iteration number of RM outer loop. .................. ................5 510. Normal and defect regions in the NCAT phantom ................. ................ ...._.60 51 1. Contrast vs. cutoff frequency of PSMLEM ................ ...............61............. 512. The tradeoff between the contrast and noise level of the NCAT phantom ............61 513. Subregions for regional image evaluation for the NCAT phantom. .......................63 514. The tradeoff between regional intensity bias and variance of the NCAT phantom..63 61. The RSCmodel channels used for channelization in singleslice CHO model.........65 62. Defect locations in the phantoms. ...._._._.. .... ..__.. ...............65.. 63. The RSCmodel channels represented as spatialdomain templates.. ......................67 64. Twostep procedure of the multiframe multislice CHOHO model. .................. .....68 65. SP images of the geometric phantom reconstructed by PSMLEM. .............................70 66. SP images of the geometric phantom reconstructed by RM. ............. ...................70 67. MMCHO detectability index results of the geometric phantom. ............. ..... ........._.71 68. SP true images of the NCAT phantom. ................ ...............72.............. 69. SP images of the NCAT phantom reconstructed by PSMLEM ........._..... ..............73 610. SP images of the NCAT phantom reconstructed by RM ................. ................ ...73 611i. MMCHO detectability index results of the NCAT phantom. ................ ...............75 612. Multiframe multislice multiview SP images reconstructed by PSMLEM.._......76 613. Multiframe multislice multiview SP images reconstructed by RM .....................77 614. Multiframe multislice multiview SP images reconstructed by RM .....................77 615. MMMCHO detectability index results of the NCAT phantom. .............. ..... ........._.78 A1. Trilinear interpolation coordinates .............. ...............90.... Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SIMULTANEOUS RECONSTRUCTION AND 3D MOTION ESTIMATION FOR GATED MYOCARDIAL EMIS SION TOMOGRAPHY By Zixiong Cao December 2003 Chair: David Gilland Maj or Department: Biomedical Engineering A new method for simultaneous image reconstruction and wallmotion estimation in gated myocardial emission tomography (ET) is presented. We implemented the method using simulated phantoms and evaluated the performance of the method. To reduce the blur in the images induced by the heart motion, gated ET uses the patient' s electrocardiogram (ECG) to trigger and acquire a time sequence of image frames, each capturing a particular stage of the cardiac cycle. However, the individual gated image frames reconstructed by conventional imagereconstruction algorithms are very noisy due to the reduced counts of detected gamma rays. Nonrigid deformation of myocardium can be characterized by means of a vector field called the motion field, which describes the relative displacement of each voxel, and thus establishes a voxelintensity correspondence between two image frames. Previous research incorporated a deformable elastic material model into a motionestimation algorithm to more accurately describe the deformation of the myocardium. The estimated motion field was used to warp each frame to a common reference volume. All warped frames were then recombined to generate an improved composite image. Instead of using the motion to warp and recombine frames, our simultaneous estimation method uses the motion field to regularize the individual reconstructed image frames, whereas the proj section data are also a constraint for motion estimation. Thus, the simultaneous method is expected to improve the motionestimation accuracy and the image quality of individual frames, relative to the independent motionestimation and imagereconstruction methods. We presented the mathematical formulation of the simultaneous estimation method, followed by an optimization algorithm. The algorithm was implemented on a simulated geometric phantom and a realistic nonuniformrationalBsplinesbased cardiac torso (NCAT) phantom. We evaluated the simultaneous method in terms of the motion accuracy and the image quality. A numerical observer model was used to assess the detectability provided by the simultaneous method for the defect detection task. Results showed that the simultaneous method produces better motion accuracy and image quality, relative to the independent motion estimation and an alternative imagereconstruction algorithm. CHAPTER 1 INTTRODUCTION 1.1 Introduction Nuclear medicine provides images of diseases by administering small amounts of radioactive materials into the patient's body and detecting the emitted gamma rays using special types of cameras. By processing the counts of gamma rays detected from various angles, a computer reconstructs images, which can provide information about the function of the body area being imaged. Nuclear medicine is unique in its ability to create functional images of blood flow or metabolic processes rather than the more conventional structural or anatomic images produced by Xray examination, computed tomography (CT), and magnetic resonance imaging (MRI) [1]. Tomography achieves detailed organ studies by separating a threedimensional (3D) obj ect into a stack of twodimensional (2D) cut sections or slices. Unlike the Xray CT that detects transmitted Xray photons, tomography in nuclear medicine detects emitted gamma rays and therefore is referred to as emission tomography, or ET in short [2]. Single Photon Emission Computed Tomography (SPECT) and Positron Emission Tomography (PET) are two widely used modalities in the clinical nuclear medicine. By measuring both the heart blood flow and the metabolic rate of the patients, physicians can find areas of decreased blood flow (such as that caused by blockages), and differentiate diseased from healthy muscle. This information is particularly important in diagnosing coronary artery diseases (CAD) [3, 4]. Myocardial ET is one way to measure blood flow and metabolism in the myocardium. SPECT imaging assesses the severity and extent of perfusion defects in patients with coronary stenosis. It accounts for over 90% of all myocardial perfusion imaging performed in the US today [5]. The motion of the patient' s heart causes blurring in myocardial ET images. This motion blurring can be reduced by gated myocardial ET, which uses the patient' s electrocardiogram (ECG) signal to trigger and acquire a time sequence (typically 8 or 16 frames) of acquisition over the cardiac cycle [69]. Data are acquired over many cardiac cycles to produce the final set of image frames. Each gated image frame records a particular stage of cardiac cycle. In addition to reducing motion blurring in each image frame, gated imaging has the advantage of allowing an estimation of cardiac wall motion and ej section fraction. However, the statistical quality of each gated image suffers as the total number of detected gamma rays is distributed over a series of time frames. Because of the low gamma ray counts, the images produced by conventional imagereconstruction algorithms usually are very noisy [4, 6]. Previous researchers [1013] presented several reconstruction methods to improve the quality of the gated image frames. Lalush et al. [10] considered the time sequence of gated frames as a fourdimensional (4D) image and reconstructed it using a penalized maximumlikelihood (ML) technique with spacetime Gibbs priors to ensure smoothness in individual image frames and between frames. Narayanan et al. [1l], Wernick et al. [12] and Brankov et al. [13] used an alternative technique based on principle component analysis, in which the time sequence of gated datasets were KarhunenLoeve (KL) transformed and then rapidly constructed framebyframe. The gated image frames were finally obtained by applying the inverse KL transform. Other previous studies [1416] presented methods to obtain a vector field describing the motion of the myocardium between successive frames. Klein and Huesman's method [15, 16] views the myocardium as an isotropic elastic material whose deformations are governed by the equations of continuum mechanics. It was shown that deforming and combining the gated image frames using the estimated motion resulted in improved composite images. Based on these previous studies, we generated a new idea to simultaneously reconstruct two gated image frames and estimate the motion vector field between them. Our hypotheses are: constraints from the material property will reduce the noise level in individual gated image frames, and lessnoise images will in turn produce more accurate motion estimation. 1.2 Significance Based on the 3D motionestimation method developed by Klein et al. [15, 16] and the conventional maximumlikelihood expectationmaximization (MLEM) algorithm [17, 18], we proposed a new estimation method for gated myocardial ET: the simultaneous imagereconstruction and 3D motionestimation method. From an image reconstruction viewpoint, the simultaneous estimation method uses the deformation of elastic material that enforces the consistency between frames as a penalty to perform penalized ML reconstruction. Unlike the usual penalized ML algorithm [10, 19], our penalty uses the deformation to regularize the reconstructed images, without any spatial smoothing. In their initial attempt, Gilland and Mair [20] used the PolakRibiere conjugate gradient optimization algorithm to solve the problem. Because of the complexity of the obj ective function, it is not clear that this algorithm converges and may produce negative image intensities. Our study developed a new twostep iterative algorithm (RM) that alternately updated the gated images and the motion vector field, and guaranteed nonnegativity in the images. This RM algorithm and some results were partly reported by Cao et al. [21]. We evaluated the motionestimation accuracy and image quality provided by the simultaneous method using a simulated geometric phantom and a realistic NonuniformrationalBsplinesbased Cardiac Torso (NCAT) cardiac phantom. We also evaluated the image quality using a numerical observer model for both phantoms. The numerical observer results have been reported by Cao et al. in "An observer model evaluation of simultaneous reconstruction and motion estimation for emission tomography" submitted to 2004 IEEE International Symposium on Biomedical Imaging, Arlington, Virginia, US. Myocardial ET provides perfusion metabolic and wallmotion information in a single imaging procedure. Our simultaneous estimation method promises to improve motionestimation accuracy and individual gated image quality simultaneously, which will improve the accuracy in diagnosing cardiac diseases. CHAPTER 2 BACKGROUND 2.1 Myocardial Emission Tomography While other cardiac imaging modalities (angiography, echocardiography, etc.) provide purely anatomical information, myocardial ET provides physiological information, such as that used to identify areas of relatively or absolutely reduced myocardial blood flow associated with ischemia or scar [1]. SPECT and PET are two main myocardial ET protocols used today. Although we presented and implemented our method on the basis of SPECT model because of its fundamental concept and simplicity, it can be expanded to estimate the image and motion for PET imaging. In SPECT myocardial perfusion imaging procedure, the patient is intravenously given a pharmaceutical bound with a radioisotope label (such as Tc99m). The radioactive pharmaceutical is then taken up in the patient's myocardium in proportion to perfusion or tissue metabolism. Single photons are emitted as the Tc99m radioisotope decays, and a fraction of these photons are able to penetrate the surrounding body tissue and reach the SPECT gamma camera [22]. The gamma camera records the spatial distribution of the photons and forms a 2D proj section of the 3D activity distribution. A 3D image of the organ is then reconstructed from the proj section data collected from various angles. Proj section and reconstruction procedures are discussed later. Patients with significant coronary artery stenosis usually have diminished radioactivepharmaceutical concentration in the area of decreased perfusion, which will cause a defective area in the perfusion image. By comparing the images taken in two cases, when the tracer is administered during stress and rest, physicians can distinguish ischemia from scar tissue. If the defect in the image taken in stress is worse than that taken at rest, it is most likely due to ischemia. If the defect remains unchanged, the lesion is most likely from scar [1]. The human heart usually points downward, approximately 450 to the left and anterior, but individual variations are considerable. Cuts perpendicular to the long axis are called coronal or short axis cuts. Cuts perpendicular to the short axis are called sagittal (vertical long axis) and transaxial (horizontal long axis) cuts [4]. Figure 21 shows anatomic orientation of the heart. Cozonal= Short Axis Figure 21. Anatomic orientation of the heart. 2.2 Projection and Image Reconstruction 2.2.1 Radon Transform Figure 22 shows the 2D data acquisition scheme in SPECT using a parallel collimator [22]. Rotating the detector allows one to observe the photon emission in the field of view from many angles. The number of scintillations detected at location s along the detector when the detector head is at an angle of 6 is defined as g(s, 6). We denote the position of a point in the 2D slice as (x, y) and the estimated number of photons emitted from this point as a function f(x, y). This estimated image function f(x, y) is assumed to be proportional to the tracer concentration of interest. Function g is the proj section of onto the detector as allowed by the parallel collimator. A sinogram is an image representative of g in which the horizontal axis represents the angular position of the detector and the vertical axis represents the count location on the detector. g(s,0) detector collimator 0 x f(x,y) object Figure 22. Principle of data acquisition and geometric considerations for SPECT. Mathematically, the proj section operator can be described by the Radon transform [23]. The Radon transform g(s, 6) of the function f(x, y) is the line integral of the values of f(x, y) along the line parallel to the collimator at a distance s from the origin: g(s, 6) = f (x(s,u), y(s, u))du (2 1) = f (s cos 6 u sin 6, s sin 6 + ucos 6)d~u where u denotes the location of the points along the integral line. Equation 21 involved a coordinate system transformation, which is illustrated in Figure 23. Assuming the two coordinate systems share the origin, the transformation between them is x = s cos 6 u sin 6 and y = s sin 6 + ucos 6 (2 2) From an operator' s point of view, the Radon transform can be represented as a procedure of an operator product: g(s, 6) = Hf(s, 6) = hs (r) f(r)dr (2 3 ) where r = (x, y) is the coordinate vector, H is called a forwardprojection operator, and h ()is a weight function that represents the probability that a photon is emitted from r and detected at (s, 6) Figure 23. Coordinate systems transformation. In numerical analysis by computers, (s, 6) and (x, y) are all represented as discrete variables, and g(s, 6) and f(x, y) are functions of these discrete variables. A discretization procedure for all the variables and the functions is inevitable. In our study, each element of the 3D obj ect image is called a voxel, and each measurement on the detector at a particular proj section angle is called a bin. In the discretized image and sinogram, the element f (j) is the jth VOxel in the image and g(i) is the ith bin in the sinogram. In the discretized domain, the sinogram and the image can both be considered as vectors. The elements in the sinogram vector are the counts of all the bins. The elements in the image vector are the intensities of all the voxels. In our study, for convenience, we use the same notation for function and vector. For example, when we use g(s, 6), g is a continuous function; when we use g(i), g is a vector with the element index of i. For the discretized image and sinogram, the forwardprojection operator H becomes a matrix. Vector g is then the matrix product of H and vector f : g = HJ; or in the form of components, g () = Hf (i) = h, f( j) i= 1, 2,.. .,p. (2 4 ) where H is now called a forwardproj section matrix with elements of h,,, q is the total number of voxels in the image slice, and p is the total number of bins in the sinogram. The entries of H can be carefully chosen to take into account the geometry of acquisition and, more precisely, the detector response, attenuation and scatter. 2.2.2 Maximum Likelihood and Image Reconstruction By comparing the proj section calculated from the current image estimate to the measured proj section data, an iterative imagereconstruction algorithm adjusts the image estimate iteratively, and eventually reconstructs an image [23]. The iterative algorithms differ in the comparison method and the correction applied to create the new estimate. Equation 24 is based on the assumption of the ideal noisefree situation. In experimental measurements, however, the counts are subj ect to randomly change under the Poisson statistics of the radioactive disintegrations. In emission tomography, the wellknown loglikelihood function can be written as L(gi, f ) = [gig(i)lgHf(i)) Hf(i)] (25) which measures the likelihood that the image produces the sinogram g. Maximizing the loglikelihood function leads to finding the most likely image fto produce the sinogram g. Maximumlikelihood expectationmaximization (MLEM) is an iterative algorithm that reconstructs this most likely image. Each MLEM iteration is divided into two steps [18, 24]: Expectation step (E step) in which a formula is formed to express the likelihood of any image given the measured data and current estimate f(k), and Maximization step (M step) in which a new image f(k+1) that has the greatest likelihood is found. The iterative imageupdating formula of the MLEM algorithm is given as [18, 24] f ~ fk (k+1 ( )=hg ,.,q MLEM shows its advantage of monotonicity of the likelihood function versus iteration [17]. The convergence, however, is very slow and the images appear unacceptably noisy before it converges [25]. To smooth the noisy image reconstructed by MLEM, a digital postprocessing filter is usually applied to the image reconstructed after many iterations of MLEM. This PostSmoothed MLEM (PSMLEM) algorithm has demonstrated excellent results in previous studies [26, 27]. The most important parameter in the PSMLEM algorithm is the cutoff frequency of the filter, which determines the smoothness of resulting images. In our study, we use 100 iterations MLEM and lowpass Hann filter because of its simplicity. Figure 24 (A) shows a sample image reconstructed by MLEM with 100 iterations. Figure 24 (B) demonstrates that the smoothness of the image reconstructed by PSMLEM is influenced significantly by the amount of the postsmoothing filter. Nyquist frequency (0.5 cycle/pixel) is usually used as a feasible cutoff. Because of this, when we compare our algorithm to PSMLEM algorithm, we usually use a range of cutoff frequencies (from 0. 1 to 1.0 cycle/pixel) for PSMLEM. Figure 24. Images of MLEM and PSMLEM. A) 100iteration MLEM image B) PSMLEM image using Hann filter with various cutoff frequencies. From 1 to 10, the cutoff frequency was 0. 1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 cycle/pixel. 2.3 Myocardium Motion Estimation A variety of methods have been proposed recently to estimate and represent the cardiac motion [1416, 2832]. Based on a deformable elastic model and the assumption that voxels corresponding to the same tissue in two frames conserve their intensity values, Klein and Huesman [15, 16, 32] presented a voxelbased method to estimate the cardiac motion vector Hield. This motion Hield is compatible for use by the simultaneous method since the images are also represented as voxelbased intensity scalar Hield. Klein's motionestimation method [32] involves minimizing a twocomponent obj ective function. The first component is an image matching obj ective function serves as a driving force to push deformed image intensities of two frames into correspondence. The second component is a regularization of strain energy function that prevents deformations that are physically unlikely or meaningless. We denote the frame 1 and 2 as f,(r) and f,(r) respectively, and the motion vector Hield as m(r) = (u(r), v(r), w(r)), where r = (x, y, z) are the voxel coordinates in the spatial domain. As shown in Figure 25, the motion vector m(r) presents the displacement from a voxel in Frame 1 to its correspondence that represents the same tissue in Frame 2. [32] Framel , Frame2: r r> mr Figure 25. Motion vector representing the displacement of voxels The image matching obj ective function is formulated as E,(J f;,m)= Cf(r) f(r +m(r))(dr (27) and the strain energy Es(m) is Es (m) = fIx +", Vy +W X dr+pus+ (28) where ii, pu are elastic weighting parameters called Lame constants that reflect the degree du of compressibility of the object, ux = is the derivative of the motion component u with respect to x. The entire obj ective function is grouped as E = EI + BEs, where a is a hyperparameter that control the balance between the two terms. In the 3D discrete numerical domain, E can be expressed as a function of all the variables (u,, v,, w, ) and their derivates, where j=1, 2, ... q denotes the voxel index in the image. E= [f (r~,)f( r +m ) +1; ap, it, ",3 l)i (29) ipp C u, + ", + .,,) u,, .L +v, i~+ ul +w, /",= l+ 3,N~ ]=1 L =1 Klein' s iterative optimization method [32] to minimize the entire obj ective function E consists of steps shown in Figure 26. Approximate E to a Find the minimum quadratic form (E) around I a point of E the current estimate m (k) denoted as m Use li as the new estimate m~kl Figure 26. Klein's iterative motionestimation optimization method 2.4 Image Observer Evaluation 2.4.1 Human Observer Study Besides the quantitative assessment such as the signaltonoise ratio and defect contrast, the image quality and the performance of a reconstruction algorithm can also be assessed by observer studies, in which the observers evaluate the images from the perspective of a clinically specific task. Defect detection in myocardial ET imaging is a specific twostate classification task that requires the observer to diagnose by viewing the image and classifying it as either normal or some particular disease state. The accuracy of the diagnosis can also be viewed as an observer measurement of the image quality. Psychophysical study with human observers (usually the physicians) is the standard means of observer evaluation of detection task performance in medical imaging [33]. For a defect detection task, there are two possible true states of the obj ect: signal present (SP) and signal absent (SA). There are also two possible responses: positive and negative, corresponding to the observer' s classification of SP and SA, respectively. Table 21 shows the conditional probabilities of a certain response given the occurrence of a certain stimulus [34]. Table 21. Conditional probabilities of a defect detection task Stimulus = SP Stimulus = SA Response True positive False positive = Positive P(TP)=P( Positive  SP ) P(FP)= P( Positive  SA) Response False negative True negative = Negative P(FN)=P( Negative  SP ) P(TN)= P( Negative  SA ) The probabilities in the Table 21 are not independent of each other. The relations between them are: P(TP)+P(FN)=1 and P(FP)+P(TN)=1. Thus we have only two independent probabilities in the defect detection problem. In an observer study to evaluate the performance of an imaging system, the human observers are shown a large number of images and asked to rate their confidence as to whether a defect is present or not [33]. The resulting P(TP) and P(FP) represent the probability of correctly distinguishing, or the detectability of the system. However, these two probabilities are obviously determined by the decision criteria the observer used. To remove the dependence on the decision criteria, a Receiver Operating Characteristics (ROC) curve [34] is plotted with P(TP) versus P(FP) as the decision criteria change. The area under the ROC curve is a common figure of merit to represent the detectability of the imaging system. The ROC curve can also be used to optimize the parameters to produce better images [3537]. 2.4.2 Numerical Observer: Background The human observer tests are very timeconsuming and difficult to perform in practice. Therefore, there is considerable interest in the use of numerical observers of which the performance indices can be calculated rather than measured. A numerical observer correlating well with the human observers would be extremely valuable for optimizing and evaluating imaging systems. The numerical observer procedure for the defect detection task is as follows: compute a scalar feature Z( f), where fis the image to be classified, then compare it to a threshold C. Choose positive if ii >C (or ii scalar feature ii is also called a test statistic or decision variable. A common figure of merit for the numerical observer is the detectability index da (also known as d) [38], defined by 22[E'(A  SP) (A  SA)]2 d = (210) Svar(il  SP) + var(il  SA) where E(AZ  k) and var(A2  k) are the conditional mean and variance of Z( f), respectively, given that fwas produced by an obj ect in state k (k = SP or SA). An ideal observer is defined as one who has full statistical knowledge of the task and who makes best use of the knowledge to minimize a suitably defined risk. The test statistic of the ideal observer is defined by [39] jldeal =IP( f  SP)P S (1 where P( f  k) is the conditional probability of given that fwas produced by an object in state k (k = SP or SA). The right side of Equation 211 is called loglikelihood ratio. (Note: loglikelihood ratio is unrelated to loglikelihood function discussed in 2.2.2). Ideal observer maximizes the area under the ROC curve and sets an upper limit to the performance obtainable by any other observers including the human, though it is only applicable for simple situations. In practice the loglikelihood ratio is rarely calculable, except for the detection of an exactly known signal superimposed on an exactly known background (SKEBKE) situation. In this case the likelihood ratio can be calculated by simple linear filtering process. If the noise is stationary, white and Gaussian, the likelihood ratio is the output of a matched filter. If the noise is stationary and Gaussian but not white, a prewhitening procedure is required before the matched filter [40]. As the performance of the human observer was found to be dramatically degraded by correlated noise [41, 42], and some spatialfrequency channels were found to exist in the human visual system [43, 44], Myers et al. [45] added a frequencyselective channel mechanism to the ideal observer. The channel model has demonstrated the ability to improve the correlation between numerical and human observers for detection tasks by other researchers [4650]. 2.4.3 Hotelling Observer and Channelized Hotelling Observer In practice, a linear observer model is usually more applicable than the ideal observer. A linear observer is defined by n( f ) = u f = Cu, f( j) (212) where u is called a discriminant vector, and T denotes the transpose. The fh COmponent value of u, denoted as u,, can be viewed as the importance assigned to f( j), the fh piXel value of the image vector f Before define Hotelling observer, we need to introduce two scatter matrices that are used to describe the first and second order statistics off[39]. The interclass scatter matrix S1 measures how far the state means deviate from the grand mean f and is defined by S, _C pf Pk( k k)T (213) where K is the number of states, Pk is the probability of occurrence of state k, fk is the state mean for the kth state. The grand mean is given by faf Pk k (2 14) k=1 The intraclass scatter matrix S2 is defined as the average covariance matrix across all states: S2 a PkKk (215) where the kth state covariance matrix Kk is given by Kk=(f kX/7k k,' (216) The angular brackets in Equation 216 have the same meaning as the overbar, that is, a full ensemble average over all obj ects fin state k. Based on the covariance matrices S1 and S2, a criterion measuring the separability of the states was proposed by Hotelling in his classic paper [5 1], and is now often referred to as the Hotelling trace: J = tr(S 'S,) (217) where tr denotes the trace (sum of the diagonal elements) of the matrix. The Hotelling trace is a common measure of classification performance in pattern recognition [39]. It increases if the variability in the image (due to noise of other factors) is decreased, since that corrnrnrespnd to reduci;ng, the covriance te~rmsn that go\ into S2 and hence increase S. It also increases if the system is modified in such a way that the two state means become more widely separated, since that increases the norm of the vectors that constitute S1. Previous research [40, 52, 53] presented that Hotelling Observer (HO) as a linear observer optimizing the Hotelling trace. The test statistic of HO is obtained by Arr = u of = (fe f)"S' f, (21 8) Implementing HO begins with creating a training set of realizations of obj ects in two states (SP and SA), and estimating the scatter matrix S2 %Om the training set. If there are N pixels in the image, then S2 is an NX Nmatrix with the rank of N The number of images in the training set should be larger than N2. Otherwise, the inverse of the estimated S2 Will HOt exist [54, 55]. The Channelized Hotelling Observer (CHO) incorporates the channel mechanism into HO [46]. Using a bank of2~bandpass frequency filters, CHO transform each image to an M~X 1 vector. The values in the resulting vector represent the frequency information of the image, so this channelization procedure can be viewed as a feature extraction procedure. Mathematically, it can be described by v =Uf (219) where U is an NX M matrix, which represents the channelization procedure and is determined by the channel model. For an SKE task, the mth row vector in the matrix U describes the mth channel's impulse response centered on the defect location. CHO calculates the test statistic of the resulting vector v by ACHO HO~ SP SA) TS V (220) where vSP and vSA are the state means of v for SP and SA, respectively, and Sv is the intraclass scatter matrix of v. CHO model has demonstrated to predict human observer performance for a variety of SKEBKE detection tasks. Gifford et al. [47] investigated the application of CHO to detect hepatic lesions using various channel models, eye filters and internal noise models. Burgess et al. [56] showed the correlation between the human observer and CHO for welldefined twocomponent noise fields (Poisson noise correlated with a Gaussian function) with various channel models. CHAPTER 3 FORMULATION OF SIMULTANEOUS ESTIMATION METHOD 3.1 Objective Functions Given two gated projection datasets g "' and g, 2, our method simultaneously reconstructs the image frames f, (r), f, (r) and estimates the motion vector field m(r) = (r(r), v(r), w(r)) that describes the motion from f, (r) to f, (r) .Here f,(r, fr) nd m(r) are presented as functions of a continuous variable in the spatial domain, r = (x, y, z), where x, y_, z are coordinates in three dimensions. We define an energy function as Elf,, f,, m)= ai(f,, f,)+ Eil is/, fm)+ PEs(m) (3 1 ) The first term is a sum of the negative loglikelihood functions of two frames L (f f, ) = H f ,, (i) J"' (i) lo gHf,, (i) (3 2 ) where m = 1, 2 denotes the frame index, i=1, 2, ..., p denotes the element index in the datasets, and H is the forwardproj section matrix that gives Hf ,c i)= [hes J, j ( 33 ) The second term in Equation 31 is the image matching term that measures the agreement of two frames under the assumed motion m: El 71,72,m)= Cf;(r)fl(r+m(r))]dr (34) The third term is the strain energy term of deforming elastic material: Es(m)= IAlrux y zI 2d+ lyu2 2 z2 (35) + pl(r) lu, v where the Lame constants 31 and p determine the material compressibility. Due to the different compressibility of the blood pool and the myocardium, we allow the values of 31 and p to change between regions by representing them as the functions 31(r), p(r). The values of the hyperparameters a and fl reflect the influence that the proj section data and the strain energy have on the estimates obtained by this method. In our study, we chose these values by experience. The reconstructed image frames f, (r), f2 (r) and estimated motion vector field m are obtained by minimizing the function E(fy, f2,m), that is, (J; f2 ,m )= argminE(;f f,m). (36) 3.2 EulerLagrange Equations We use the Frechet derivatives of the function E to obtain the EulerLagrange equations. As the first step, we compute the Gateaux differentials [57] of E with increment t. Let DmEV;fif, m; t), DuEV;fif, m; t), DEfifi,f2 m; t), DwEV;fi, fm; t) denote the Gateaux differentials of E with increment t relative to fm, u, v, w, respectively. The Gateaux differentials are computed by using the wellknown formula as DE(.fif2,m;t)= ~E~f, +srt, f2 0~s (37) Applying similar formulas to the functions EI, Es, L, which determine E, for differentials relative to all the variables, we obtain the Gateaux differentials D,E, if2,, ,m; t)= 25f [(r) ~f2(r + m(r))]t(r)dr (3 8) D2E,(f,, f2,m; t)= 25 [f, (r) f2(r + m(r))]t(r + m(r))dr (3 9) DmL~,, 2;t= [H~i) ge)(i)t~i/Hf~i)] m=1,2(310) Dz,:(I,72m~t= 2 [,(r f 2(r + m(r))t(r)dr (311) DE(II72,~t) 2[/() 2 2(r + m(r))t(r)dr (312) DzE,E(m; t) = [(jlV m + 2puxl)tr + pl(u~ + vx)tv + pr(uz, x )tz ]dr (314) DE,E(m; t) = [(lV m +2pny,)ty + p(ul + vr)tz*l 1Vz + y z]dr (315) D_,E (m; t) = [(AV m + 2pez)tz + pr(uz + x)tr #r(' y y~t ]dr (316) Now, let Cc (R) denote the set of all infinitely differentiable functions with compact support on s Then the EulerLagrange equations are DmE( fi, f2, m; t) = 0 , Dz,E( f,, f2, m; t) = 0, D,E( f,, f2, m; t) = 0, and D,E( f,, f2, m; t) = 0 for m = 1, 2 and all t E Cfm All these equations involve integrals containing the function t. To remove the integral and the dependence on t from these equations, note that if Sg(r)t(r)dr = 0 for all t E Cfm then g = 0. To use this fact, each of the Gateaux differentials of E needs to be expressed in the form Sg(r)t(r)dr for some g. Such a function g is called a Frechet derivative. The Frechet derivatives of E with respect to f,,, u, v, w are denoted by D~,,E Dz,E, DE, DwE respectively, and are weighted sums of the Frechet derivatives of EI, Es, and L. Some Frechet derivatives of E are immediately obtained from Equations 38, 311, 312, and 313 by easy calculation: DE,(fi,, 2,m) = 2[f,(r) f2(r + m(r))] (317) D E, (f,,f2, m) = 2[ f, (r) f2( 2 ~r) (r + m(r)) (318) DE E(II,2, fm) = 2[J;(r) f2( 2 ~r) (r + m(r)) (319) DE E(II,2, fm) = 2[f, (r) f2( 2 ~r) (r + m(r)) (320) To obtain D2Ef;(, f2, m), We need to express the integral in Equation 39 as g (r)t(r)dr . Consider the coordinate transformation r = r + m(r) = (I3 + m)(r) Where I3 is the identity 3 X 3 matrix. Then, the Jacobian is given by dr JWxr) [det(I 3 (321) where m'(r) =~ vx vy (2 So, if a is the image of R under the operator (I3+m), fTOm Equation 39 we obtain D2E:(fi72, ; t = 2 [ (IS F) (F)t(F dr(323) Now, the motion changes the location of pixels within the ROI R, but it leaves the entire ROI unaltered, hence R = 0 So, from Equations 321, 323 we obtain D2E, (ft,, 2,m) = 2[f, (r) f2(r)]/I det(3 + m' (r))1 (324) where r = (I; 1)(r 24 The Frechet derivatives of L are given by Hf, (i) (325) where h(r) [ ha (r) So, from Equations 317, 324, 325 and the fact that E does not =1 depend onf,,, we obtain the Frechet derivative of E with respect to fi and f2 g~'Hf, (i)(ih r (326) D2E( ,f,, fm ) We now determine the Frechet derivatives of Es with respect to u, v, w. The Gateaux differentials in Equations 314, 315, 316 are all of the form SF Vtdlr for some vector field F. By using the divergence theorem, for all t e Cf (R) SF Vtdr = (V o (tF) (V F)t)dr = 5 (V F)tdr (3 so the Frechet derivatives are V F. From Equations 314, 315, 316, we see that the vector field F for De,Es(m; t) D,Es (m; t), and D,Es (m; t) are (AZ(ux + v, + w,) + 2puux, p(u, + vx), pu(u, *x)) (3 (pU(u, + vx), A(ux + vy+ wZ) + 2pUy ,p(vZ y )) (3 (pU(u, *), pU(v W), *y x y z,)+ 2puw,) (3 respectively. In each case we obtain V F as 28) , 29) 30) 3 1) 2[ f (r 2 ()]/detI3 m'(r) + ~ r) ()H (i)(ih, (r)1 (327) (AZ + 2pu)uxx + pu(u, + uz)+ (/ + Pu)(xy xzw) (332) (AZ + pu)(uxy+ yz)+ xxv, z) + (AZ + 2pu)v,, (333) (AZ+ pu)(uxz + vyz)+ Pu(w + yy)+ (AZ+ 2pu)waz (334) Hence we obtain the Frechet derivatives of Es Dz,Es (m) = (AZ + 2pu)uxx pu(u, + uzz) (/ + Pu)(xy xzw) (335) DV,Es(m) = (AZ + pu)(uxy + yz) uvxx zz) (A + 2pu)v,, (336) Dw,Es(m) = (AZ + pu)(uxz + yz) u wx yyw) (AZ + 2pu)waz (3 3 7) Since L is independent of m, we obtain the Frechet derivatives ofE relative to u, v, w as Dz,E(f;, f2,m) = P[(il+ 2p)uxx pu(u, + uz) (i #~)(xy xz, (3f (33 8) 2[ f (r) f2( 2 ~)l ( ~) DE(f,, f2,m [i u(xy yz) uvxx zz) (il+ 2U)vw] afi (339) 2[J; (r) f2( 2 ~)l ( ~) DwE(f,, f2,Il [i ~ u(xz vyz uxx yyW) (jl+ 2p)waz i dfi (340) 2[f, (r) f2( 2 ~)l ( ~) Finally, the EulerLagrange equations for E, satisfied by (f, f2 ,m ), are 2[ zr)f2(~m~r)]+ h~) h(r~ *()/Hf(i)= 0(341) 2[ z(r) f2(r)]/ det(I3 + m'(r) + hj(r) h(r)rj()i/H2i (342) 2[ (r) f2 2ml~f Dx (343) B(A+2pU)ux + p(ul~ +uzz) (~ )(V xy xz1?) = 0 dy (344) Pi(A + 2p)ul, + p(v, + v ) + I~(A + p)(u +w = 0 2[f, (r) f,(r + m) r +m) 8: ( 3 45 ) B(A+2r)w_ +pI(w +Mw )+(Al+p)(ux= f's) =0 The Equations 343, 344, 345 are valid on regions with uniform Lame constants 32 and pu. So we assume throughout our study that the frames are segmented into two regions with uniform values for the Lame constants on each. 3.3 Optimization Algorithm: RM 3.3.1 Scheme of Optimization Algorithm We estimate the minimizer ( f, f m ) by a twostep iterative procedure, which updates the estimates ~f,(k),>2(k) m(k) iteratively to ~f,(k+1)> 2(,k+1) 'k+1) in the kth iteration. Figure 31 shows the scheme of this optimization algorithm. In R step (R for reconstruction), we assume the motion is reasonable and fixed, and minimize the sum of the image matching and likelihood functions to generate the image updates, that is, (fi (k+1), f(k+1)) = arg min Elf,, f,,m (k)) (346) = arg minaol(f,, f2) +E~ (i i /2,m(k) ) In M step (M for motion), we assume the images are reasonable and fixed, and minimize the sum of the image matching and strain energy functions to generate a new motion estimate, that is, ml(k+1) = arg min E( f,(k+1), 2(k+1),m = arg min E (fi"k+1) 2(k+1),m) Ss(m)] (347) Hence this iterative optimization algorithm is called "RM" algorithm. The image frames and motion vector field are simultaneously updated after one repetition of these two steps, called an "outer loop" of the RM algorithm. The minimization problems in R and M steps are also optimized using iterative algorithms, and therefore, they are referred to as "inner loops". J;(k), 2(k) R Step: Minimize f,'k+1) 2(k+1) M Step: i(k+1) f2(k+1) (k) Outer Loop Figure 31. Scheme of the RM algorithm. 3.3.2 Computation of R Step Equation 346 can be regarded as a penalized maximumlikelihood approach of image reconstruction in which the penalty term enforces natural constraints due to the limitedness of material deformations, rather than forcing (possibly unnatural) prior smoothing constraints. Green's onesteplate (OSL) algorithm [58] was proposed for the penalized maximumlikelihood estimation of a single frame using a general class of priors. It was used [10] for reconstructing gated frames, and convergence was obtained for a modified version [59]. We used a modified OSL for our Rstep optimization. Since EI and L are convex infi, f2, the KuhnTucker optimality conditions are necessary and sufficient for the minimizer. That is Dm E(fi, f2," (k) ) > 0 (348) and fmDmE( fi,f2, m(k)) = 0 for m =1,2. (349) From Equations 326, 327, we define af ha(r) g (i)Hf; (i) Az (f,, f27M m> =1 (350) ah7(r) + 2[ (r) f2( ~ afha(r)gl?(2) Hf2 A2 17, 2, _> 1=1 (351) ah(r) 2[ (r) f2 (r)]  det(I3+m'(^) Hence the necessary and sufficient conditions for Equation 346 are Am1 2mk) :, m =1,2 (352) fAm(11,2,m~) mr''>=, m =1,2 (353) So, it seems natural to estimate f,(ki+1) by the fixedpoint iterations given by my I 1+ = fml,t~l $1(f; 2f,17 (k)), l= 1,2,... (354) and f,( k+) p ~kl) If the image matching term EI is not present in the obj ective function in Equation 346, then Equation 354 reduces to MLEM algorithm. If E is replaced by a Gibbs prior, this approach yields Green's OSL algorithm. However, the multiplier Asm is not guaranteed to be nonnegative. We use the following modification of the iteration in Eqluation 354 to generate nonnegative approximations / 0~, J~ ~fk) ... off kll+1) R Step Alg~orithm 1. Initialization: f, o)_ = (k) 2. Given Jt compute the values of At (r) =' As '( f ,f fmk)FO 3. Determine constant tl > 0 satisfying t, mamaxax~lA' zr) ,0):rf e ,, j=; 1,2)< 1 (355) 4. he ff',= f '+ tz (A' ) is the update of ~f l. 5. Go back to 2. To apply this algorithm we need to compute the variable r = (I, + m) (r) for each re E which occurs in the quantities A','t(r) f u.. Hwee, the mathemaical form of the transformation m(k)(r) is completely unknown. Only numerical values of mk()are available, so we choose to use a fast, direct method of obtaining an approximation of i' by using a linear approximation of mk(k). By definition, r = r; + mtr) = r + m(r) + m' (r)(r r), hence, m(r) = (r, + m'(r))(r r), which gives the approximation as r rI,+m(r))m~r)(356) This approximation is valid when 1 is not an eigenvalue of m'(r) . 3.3.3 Computation of M Step M Step is to compute the motion update m(k+1) after the image frames have been updated to ~f,(k+1) 2(k+1)1 This has been reduced to the motionestimation problem considered by Klein and Huesman [15]. In the iterates of M step, the update mk) of m k) is obtained by minimizing a simpler, quadratic approximation of the original obj ective function, using Taylor approximation (357) Hence, approximately we need only to compute m'k+1) arg min Qk,/,m) + PEs (m)] (358) where Qk,l ( k)=5 i')(r _i 2(k (k _~ic, _mr (kx)cr, (kI)c mi(k)~ (359) From the Frechet derivatives of Qk,l and Es, the EulerLagrange equations for m = mI are 2 i(k) 2i(k) (k) _II (k) )7 2 (k)( (k 2I ~ nd(k)(k Dx (360) + P(t/+ 2pU)ux + iuO,.+> uzz (;3 xyU),' xz,) = 0 2 (k,() 2() ( k m ") _m (k) 2(k (k) (k):(k Sy (361) + P + 2pl)v,, + pl(v, zzv)+(i + pl)(uxy ~ y)] = 0 2 () 2() k) (k 2(k) () ] 2(k) (k)) 8z (362) + P(2 + 29UM~ zz yyU' xz yz~+~ + ~(* ),r) = 0 Using finite difference approximations for the derivatives, these partial differential equations (PDE) can be changed to an algebraic linear system for the motion components at all voxels, that is, (u,, v,, w,, j= 1,2,... q). The algebraic linear equations are solved using the regular conjugate gradient algorithm. The estimates m k) (k) (k) then iteratively obtained by solving Equations 360, 361, and 362 with the initialization m k) (k) Appendix gives details of the numerical analysis for motion estimation. CHAPTER 4 IMPLEMENTATION AND EXPERIMENTAL RESULTS: GEOMETRIC PHANTOM 4.1 Simulated Source Object and Projection Data Since it is difficult to obtain a standard data set that presents the "truth" volume, using real medical imaging data to test an algorithm is usually not applicable. To demonstrate some of the basic characteristics of the RM algorithm, this chapter implements the algorithm with a simulated 3D geometric phantom. Based on a hollow semiellipsoid model, the source obj ect of this phantom is designed to simulate the contracting left ventricle in a gated myocardial ET study. The inner and outer boundaries of the myocardium are defined by 2 concentric semiellipsoid surfaces with dimensions of short and long axes chosen to mimic the shape and size of a regular left ventricle. Three orthogonal views of the source obj ect are shown in Figure 41. 4.1.1 Steps to Generate Geometric Phantom The simulated geometric phantom is represented as a 30 X 30 X 30 voxel density volume. In order to create a more realistic source distribution in which obj ect boundaries have some degree of smoothness, we generated the 30 X 30 X 30 source obj ect from a highresolution (120 X 120 X 120) obj ect with the steps illustrated in Figure 42. Two simulated data files are generated for each frame: a noiseadded proj section data file and a true source obj ect image file. The proj section data file will be used to reconstruct the image by the RM algorithm and the PSMLEM algorithm, and the true image will be used to evaluate the reconstructed images. A Figure 41. Source objects of the geometric phantom. A) Frame 1. B) Frame 2. For both frames, the 3 orthogonal planes are in 2 long axis slices (two left columns) and a short axis slice (right column). 4.1.2 Dimensions of Geometric Phantom We generated 2 frames of the gated ET image (Frame 1 and Frame 2) representing the endsystolic and midsystolic phases of the heart cycle, respectively. The source obj ects of the 2 frames had different lengths of semiaxes in a manner that achieved a constant myocardium volume across the 2 frames and an anatomically realistic myocardialtochamb er volume ratio. Assuming a nominal voxel linear dimension of 3.5 mm, the myocardial volume was 120 ml, and the chamber volumes of Frame 1 and Frame 2 were 43.3 ml and 74.4 ml, respectively. These volumes were chosen to represent the normal human heart at these contraction phases [60]. The lengths in mm of the ellipsoid semiaxes for the 2 frames are shown in Table 41. Two of the three ellipsoid semiaxes were equal in length and represented the circular, shortaxis slices. Nyquist cutoff response Figure 42. Simulation steps of the geometric phantom Table 41. Ellipsoid semiaxes lengths for source object. (Unit: mm) Frame 1 Frame 1 Frame 2 Frame 2 Axis Inner boundary Outer boundary Inner boundary Outer boundary Short axis 21.0 38.5 26.2 40.6 Long axis 70.0 78.8 77.0 84.0 4.1.3 Intensity and Defect A uniform intensity was assigned to the myocardium region between the half ellipsoids, and zero intensity was assigned to the region in the inner chamber and external to the myocardium. For the contrast study purpose, two myocardial defects were simulated with 50% decreased intensity and located at basal and midventricular points along the long axis, as indicated by the arrows in Figure 41. The total counts of the proj section for each frame was scaled to 99,000, which was chosen to match that of a clinical, gated myocardial perfusion study using Tc99m sestamibi (70 y.o. male, 10 mCi injected activity). 4.1.4 Wall Motion Simulation To mimic more realistic movement of heart, a "wringing" motion was simulated in addition to the contraction. This motion was found in contrast ventriculography [61] and tagged MR images of the human heart [28]. We characterized this wringing motion by a rotation of the basal and apical myocardium about the heart long axis in opposite directions. In our study, the extreme basal and apical points on the heart long axis rotated 5 degrees in opposite directions, and points between rotated according to a linear gradient connecting these extremes. We calculated the motion geometrically and referred to it as the "ideal motion". To illustrate the wringing motion visually, the ideal motion fields in 3 shortaxis slices (#5, #15 and #20) are shown in Figure 43. 4.1.5 Forward Projection Simulation and Detector Response Collimatordetector response was simulated assuming a lowenergy highresolution parallelhole collimator. The proj section data were filtered using a Gaussian kernel with FWHM=1.9 pixels. With the nominal 3.5 mm pixel size, this kernel has an equivalent FWHM equal to 6.7 mm. Forward proj section was performed at 60 angles over 180 degrees. Effects of attenuation, scatter, and randoms were not considered. Figure 43. Ideal motion fields of the geometric phantom in 3 shortaxis slices. A) Slice #5. B) Slice #15. C) Slice #20. The highresolution proj section data were then collapsed to a lowresolution proj section (30 X 60 X 30) by a procedure described in Figure 44 (only a 2D slice is shown). While the final projection dimensions (30 X 60 X 30) are small relative to those used by large fieldofview imaging systems, the sampling rate used here, relative to the size of the heart, is comparable to these systems. For example, the thickness of the myocardium in humans at endsystole is typically close to 14 mm; in our heart model this thickness is approximately 4 pixels, which represents a nominal pixel size of 3.5 mm. Finally, Poisson noise was added to this scaled, collapsed proj section data. 4.2 BiValue LamC Constants Model The compressibility degree of the myocardium, the blood pool in the chamber, and other outside adj acent tissue differs from each other. Thus the Lame constants 32 and pu should have different values for these regions. In this chapter, due to the assumption of constant myocardium volumes across the 2 frames of the geometric phantom, and the fact that the region of the blood pool is not conserved during contraction, we propose a bivalue Lame constants model: a fairly incompressible model for myocardium and a compressible model for the other regions. Bin 1 1 2 1 . Bm I I 1le 60 Projectionn Angle Highresolution 120 X 60 Proj section Figure 44. Collapse highresolution proj section to lowresolution proj section In practice, there are two more often used constants presenting the compressibility of material: E, called Young' s elasticity modulus and v called Poisson ratio, which are related to Lame constants as pu(3il+ 2p) ii E=~ and v=~ (41) (a + p)> 2(il + pu) Young' s elasticity modulus E relates the tension of the object to its stretch in the same direction. Poisson ratio v is the ratio between lateral contraction and axial extension [62]. The ii term in the constraint equation penalizes nonzero divergence and the pu term penalizes sharp discontinuities in the flow field [32]. For highly incompressible material, the Poisson ratio approaches 0.5, which yields a 3 approaching infinity. At the other extreme situation of highly compressible material, 3A and v approach 0. Table 42 gives the Lame constants we used in our study for some specific Poisson ratios. Table 42. Poisson ratio and Lame constants Poisson ratio (v) 0 0.1 0.2 0.3 0.4 ~0.5 Lame constants 3A=0; 3A=0.25; 3A=0.67; 3A=1.5; 3A=4; 3A=100; (n, p)> iu= ; pu= ; pu= ; pu= ; pu= ; u= 1; To implement the bivalue Lame constants model, a preprocessing segmentation of the left ventricle from the background is required. Since the accurate automatic segmentation of the myocardium boundaries may be a formidable task itself, in our study we used a simple automatic thresholdbased procedure. Firstly, we reconstructed the image from the noisy data using the PSMLEM algorithm with Nyquist cutoff frequency. Then we compared each voxel intensity of the resulting image to a prechosen threshold and generated a segmentation bitmap. A feasible threshold was reassured by examining the segmentation bitmap visually. Because the segmentation is only to model elastic characteristics of different tissue types, the RM algorithm is expected not to be sensitive to slight errors of a few voxels along the myocardium boundaries. This insensitivity has also been demonstrated by the insignificant difference of the images reconstructed from using various thresholds for segmentation. The segmentation procedure and results are illustrated in Figure 45. In our study, the threshold was chosen to be 1/8 of the mean intensity of myocardium voxels. The segmentation bitmaps generated from the images with and without defect were similar to each other. A Figure 45. Segmentation for bivalue Lame constants model. A) With defect. B) Without defect. (Note: Left column are the images to be segmented with threshold. Right column are resulting segmentation bitmap where white region represents the myocardium and black region represents the background.) 4.3 Convergence Properties of RM Algorithm Since the RM algorithm is to minimize the obj ective function iteratively, its convergence property can be illustrated by the plot of the obj ective function value versus iteration number. Similarly, the convergence of the obj ective functions considered in the individual R and M steps illustrate the individual convergence properties. The convergence results for the RM algorithm are shown in Figure 46. The total obj ective function E and the negative loglikelihood function L decrease with iteration and are approximately constant after 20 iterations. Due to the initial zero motion estimate and initial uniform image estimates, the EI and Es functions present low initial value and increase with iteration before flattening out. Figure 47 shows the convergence of the individual R step. R step using either the true or the estimated motion largely converges by 25 iterations. The M step also has converged by 25 iterations as indicated in Figure 48. The obj ective functional value in M step achieved with the PSMLEM images is larger than that with true images. This reflects the inconsistencies between the estimated frames by PSMLEM due to the statistical noise. 4.4.1~~ Fiue oft Meri Th moio esiato iseautdb oprn h siae oint da moin lbal adreinll.Bae o heiag achn fntin w efn glba moinerrta esrsteareeto w reiaefae ne h esiae mto.It rerset th vrl cuayo temto siain 40 2.4 x 105 2.5 2.6 Using ideal emotion 2.7 U~sing: estimated motion 2.8 32 10 20 30 40 50r Iteration number Figure 47. Convergence of R step 11000 S Using ideal images 10000 ." Using PS MlLEM images 9000  8000 7000 5000 4(000 3000 0 10 20 30 40 50 Iteration Number Figure 48. Convergence of M step [,(r)f2F IT))2(42) where f, and f2 are the true images, and in is the estimated motion vector field. Regional motion evaluation is performed at some selected subregions. We define magnitude error as the absolute value of the magnitude difference between the estimated motion and the ideal motion, and angle error as the angle between the two vectors. 4.4.2 LamC Constants Study Results As we discussed in 4.2, motion estimation involves the assignment of Lame constants for myocardium and background regions. To investigate how the Lame constants influences motion estimation, we performed M step using true images with various Poisson ratio v (equivalent to various Lame constants, see Table 42) and hyperparameter B. We used six values (0, 0.1i, 0.2, 0.3, 0.4, and 0.5) of v for myocardium, and v =0 for the background. The results are shown in Figure 49. The global motion error using v = 0.5 for myocardium was much more sensitive to a than those using other v 's. The minimum global motion error was achieved with v = 0.4 and B = 0.02. Our results agree with Klein's claim [32], that is, for the fairly incompressible deformation, the motionestimation algorithm performs best with a chosen Poisson ratio value near 0.4, instead of the value of 0.5 for the ideally incompressible material. Thus in our following geometric phantom study, we kept the Poisson ratio at 0.4 for the myocardium and 0 for the background. 80 7 0 60 v=0.4 O ' 30 20* * O 0.02 0.04 0.06 0.OB 0.1 Hyperpararneter in M step: P Figure 49. Global motion error vs. B for several Poisson ratios. 4.4.3 Global Motion Evaluation To compare the motionestimation performance of RM and the individual motionestimation algorithm, we implemented the following 3 methods: Method 1: M step using true images with change ofa. Method 2: Individual method. M step using PSMLEM images with change of cutoff frequency and B; Method 3: RM algorithm with change of hyperparameter a and B. Figure 410 shows the ideal motion superimposed on the true Frame 1 images of the geometric phantom, in both long and short axis views. Figure 411 shows the motion estimated by Method 1, superimposed on the true images. Figure 412 shows the images and motion individually estimated by Method 2. Figure 413 shows the images and motion simultaneously estimated by RM. The minimum global motion error each method achieved is given in Table 43. The optimum conditions to produce the minimum errors for each method are also given. The numbers in parentheses are the standard deviations of the global motion error calculated across a 20image ensemble. Method 1 and 2 used 40 iterations for M step, and Method 3 used 40 outerloop iterations, each containing 1 Mstep iteration and 1 Rstep iteration. Figure 410. Ideal motion vector fields superimposed on true images. Figure 411. Estimated motion by M step using true images. (B = 0.02). Figure 412. Individually reconstructed images by PSMLEM and estimated motion by M step. (cutoff = 0.3 cycle/pixel; B = 0.01). Figure 413. Estimated motion and reconstructed images by RM. (a=0.02; J 0.01). Table 43. Global motion evaluation results of the geometric phantom Method Optimum Conditions Minimum global motion error 1 B = 0.02. 21.4 2 cutoff = 0.3 cycle/pixel; 46.4(1.2) B = 0.01. 3 a=0.02; 37.2(1.4) B = 0.01. The results in Table 43 show the smallest motion error was achieved by M step using the true images, as we expected. Of the methods those were operating on noisy data, the RM algorithm outperformed the individual method in terms of the global motion error. This result reflects the increased noise control of the RM algorithm. 4.4.4 Regional Motion Evaluation Three subregions within the myocardium of Frame 1 of the geometric phantom were selected for the regional motion evaluation: (1) basal and in normal region, (2) basal and in defect region, (3) midventricle and in normal region. Each subregion was a 2 X 2 X 2 voxel volume. Their locations are indicated by the squares in a slice view in Figure 414. To obtain the magnitude and angle errors of motion estimate, the mean magnitude and angle of the estimated motion in the subregions were computed and compared with the mean magnitude and angle of the ideal motion in those subregions accordingly. These errors were computed for a 20image ensemble, and the mean and standard deviation of the errors across the ensemble were obtained. The regional motion evaluation was performed for both the RM algorithm and the individual method. The results of the regional motion evaluation are shown in Table 44. The mean magnitude errors are in units of pixels, and mean angle errors are in units of degrees. Standard deviations are given in parentheses. Overall the table indicates the improved motionestimation accuracy for RM compared with the individual method. Figure 414. Chosen region to calculate regional motion error for the geometric phantom. Table 44. Regional motion evaluation results of the geometric phantom Individual method RM Subregion Ideal mag. Mag. error Angle error Mag. error Angle error 1 1.52 0.18(0.17) 31.7(14.1) 0.18(0.12) 28.7(11.5) 2 1.52 0.49(0.28) 30.1(14.8) 0.45(0.22) 25.1(11.1) 3 0.66 0.24(0.37) 61.0(16.8) 0.02(0.22) 43.3(17.3) 4.5 Image Reconstruction Results 4.5.1 Figures of Merit The RM images were evaluated both globally and regionally by comparing to the PSMLEM images. The evaluation methods include global image evaluation, regional image evaluation, and defect contrast evaluation. Since the cutoff frequency of the postfiltering used in the PSMLEM algorithm can greatly affect the quality of the resulting images, a range of cutoff frequencies was investigated. The figureofmerit for the global image evaluation was the Root Mean Squared (RMS) image error defined as RMS = frfr43 where f (r) is the true image, f (r) is the reconstructed image, and Nis the total number of voxels in the image. To get an overall image evaluation of the 2frame system, the mean RMS error was computed across the 2 frames. The regional image evaluation also considered the 3 subregions used for the regional motion evaluation. For each subregion, the mean bias and variance of the intensity of the estimated image were computed across a 20image ensemble. For the PSMLEM algorithm, a set of bias/variance points was computed across a range of cutoff frequencies from 0. 1 to 1.0 cycle/pixel. The bias/variance curves were plotted for each subregion to show the tradeoff between the absolute value of the bias and the variance, as the cutoff frequency changed. Defect contrast is defined as I I Contrast = normal defect (4 4) normal where Inormal and Idefect are mean image intensities within the normal and defect regions, respectively. The normal and defect regions were defined from the known source obj ect (true image) and included the entire region volume except for boundary voxels. Contrast was measured for each image out of the 20image ensemble, and then the mean contrast was computed. Like the regional image evaluation, a range of cutoff frequencies was considered for the PSMLEM algorithm. 4.5.2 Global Image Evaluation The results of the global image evaluation for RM and PSMLEM are given in Figure 415. The RMS error of the RM algorithm decreases with iteration until approximately iteration 20, then increases very slightly. This demonstrates the stability of the RM images with iteration, unlike MLEM. The PSMLEM algorithm shows a minimum RMS error with a cutoff frequency of approximately 0.4 cycles/pixel. This 47 minimum, however, is still substantially larger than the stable RMS error of the RM algorithm (0.045 compared with 0.03). The hyperparameter effects on the RM images are illustrated in Figure 416. The figure shows the image of Frame 2 as a and a are varied. As a is reduced, there is less reliance on the likelihood function (i.e., the measured proj section data), which causes the image of the myocardium to broaden and approach a more uniform intensity image. As a increases, the algorithm becomes more like MLEM, and a moderate increase in noise is evident. The sensitivity of the RM image to change in B is much less. Further illustration of the hyperparameter effect is shown in Figure 417. These surface plots show the change in the global motion and image errors as a 2D function of a and B. The minimum RMS error is located at a=0.02 and 8=0.04. When a was smaller than 0.02, the RMS error increased dramatically and became irregularly sensitive to B, which corresponded the image shown in Figure 416 (B). It r t o n m e r r R 1 m r o o s A a t n s .I . .1 o l Figre41. MS rrrsof heRM ndPSLE iage. ) R. ) PMEM 48 A B C D Figure 416. Hyperparameter effects on the RM images. A) a = 0.03, B = 0.01. B)B unchanged, smaller a (0.015). C)B unchanged, larger a (0.06). D) a unchanged, smaller B (0.002). E). a unchanged, larger B (0.1i). 0.07 0.065 ~006 S0.055 2 005 .0l 0.045 3 0 04 0.035, 0.1 0.080. Figure 417. RMS error of the RM images as a function of a and B 4.5.3 Defect Contrast Evaluation The results of the contrast measurements are shown in Figure 418. With the PSMLEM algorithm for both defects, the defect contrast increases as the cutoff frequency increases and, at the highest cutoff frequency, is slightly less than that of the RM image. Although at this cutoff frequency the 2 methods have similar defect contrast, the quality of the images is substantially different, as shown in Figure 419. S0.36 RM: 0.41 0.2t RM: 0.29 0.4 0.6 0.8 1 1.2 0.4 0.6 0.8 1 1.2 PSMLEM cutoff frequency (cyclelpixel) A PSMLEM cutoff frecluency (cyclelpixel) B Figure 418. Defect contrast of the PSMLEM and RM images. A) Defect 1. B) Defect 2. Figure 419 (A) and (B) show the images of Frame 1 for RM (40 iteration) and PSMLEM (1.4 cycles/pixel cutoff). It is evident that while the images have similar contrast, the RM image has substantially lower noise level. Also included in Figure 419 (C) are the PSMLEM images with Nyquist cutoff frequency. While these images have a noise level closer that of the RM images, the defect contrast is smaller. This is particularly evident in the short axis images in the right column. 4.5.4 Regional Image Evaluation The bias and variance of the voxel intensity within the 3 subregions (indicated by the squares in Figure 420) were calculated in both 2 frames. Figure 421 shows the plots of the absolute bias versus variance for the PSMLEM and RM images. In each plot, the RM bias/variance point is downward and left to the PSMLEM bias/variance curve, indicating that RM produced smaller bias for equal variance, or smaller variance for equal bias. This demonstrates that the RM algorithm achieved better regional image estimation, relative to PSMLEM. In summary, in this chapter we implemented the RM algorithm using a simulated geometric phantom that mimics the human left ventricle. We evaluated the motion accuracy by comparing the motion estimated by 3 methods with the ideal motion. The motion estimated by RM using noisy proj section data was superior to that estimated individually by applying the motion estimation method to PSMLEM reconstructed images. Also we evaluated the image reconstruction by comparing the RM images to the PSMLEM images in terms of the global image error, the regional intensity bias/variance, and the defect contrast. All the results showed that the RM algorithm produces better image quality. In Chapter 6 we will assess the image quality by observer study. Figure 419. The tradeoff between the contrast and noise level of the geometric phantom. A) RM images. B) PSMLEM images at similar defect contrast but higher noise level. C) PSMLEM images at similar noise level but lower contrast. Figure 420. Subregions for regional image evaluation for the geometric phantom. A) Frame 1. B) Frame 2. SPSMLEM \+ RM O 0.05 0.1 Variance 0 0.05 0.1 Variance + PSMLEM  + RM 0.05~ U.055 0.05 CO 0 045 0.035 0.03 n ncE SPSMLEM + RM m 0.08 ." 0.06 n1 0.04 0.02 0.12 0.1 .m 0.OB S0.06 .0 0.04 0.02 0.05 Variance 0.1 D Variance + PSMLEM + RM S0.03 S0.02 0.01 U O 0.05 0.1 Variance 0.15 0.05 0.1 Variance Figure 421. The tradeoff between regional intensity bias and variance of the geometric phantom. A) Subregion in basal normal region of Frame 1. B) Subregion in basal normal region of Frame 2. C) Subregion in basal defect region of Frame 1. D) Subregion in basal defect region of Frame 2. E) Subregion in midventricle normal region of Frame 1. F) Subregion in midventricle normal region of Frame 2. CHAPTER 5 IMPLEMENTATION AND EXPERIMENTAL RESULTS: NCAT PHANTOM 5.1 Simulated NCAT Phantom Besides the geometric phantom, we also implemented the RM algorithm on a more realistic, complex physical phantom: 4D N~URBSBased Cardiac Torso (NCAT) cardiac phantom. The NCAT phantom simulation software was developed at the Johns Hopkins University and fully described by Segars et al. [63, 64]. It uses NonUniform Rational BSplines (NURB S) method to define the geometric conditions of a human heart and produce the source obj ect of heart at a series of phrases over the cardiac cycle. The NCAT software also generates the myocardial motion vector field for the gated image frames from the tagged MRI data of normal human patients. We used the NCAT software to generate 8 gated image frames of the cardiac phantom, each represented as a 32 X 32 X 19 volume with a voxel linear dimension of 4.0 mm. The volume of the left myocardium was approximately 85 ml. We selected the 2 frames with maximal deformation to implement the RM algorithm: enddiastole and endsystole, which we referred to as Frame 1 and Frame 2, respectively. A defect region with 50% decreased voxel intensity was initially imposed on one side of the left ventricle in Frame 1. By deforming the defect in Frame I with the ideal motion field, we could impose the defect region onto Frame 2. The volume of the defect was approximately 5 ml. Similarly to the geometric phantom simulation, we scaled the total counts of each frame to 99,000 and added Poisson noise. The detector response was simulated with a FWHM of 7.5 mm. Figure 51 shows the longaxis and shortaxis slice views of these 2 frames. T B Figure 51. Long and short axis views of the NCAT cardiac phantom. A) Frame 1: enddiastole. B) Frame 2: endsystole. 5.2 Motion Estimation Results 5.2.1 Global Motion Evaluation As indicated in Klein' s study of the physical phantom [32], the best motion estimation is achieved with a Poisson ratio near 0.45 for the myocardium. Thus in this chapter we use 0.45 for the Poisson ratio of myocardium and the Lame constants are ii = 9, pu = 1 for voxels within the myocardium; Ai = 0, pu = 1 for voxels of the background. Similarly to the geometric phantom study, motion estimation for the NCAT phantom was also evaluated in terms of the global motion error and the regional motion error. We also compared the global motion error of the following 3 methods: *Method 1: M step using true images with change ofa. Method 2: Individual method. M step using PSMLEM images with change of cutoff frequency and B; Method 3: RM algorithm with change of hyperparameter a and B. Figure 52 shows the ideal motion superimposed on the true frame 1 images of the NCAT phantom, in both long and short axis views. Figure 53 shows estimated motion by M step using true images, superimposed on the true images. Figure 54 shows the images and motion individually estimated by Method 2. Figure 55 shows the images and motion simultaneously estimated by RM. Figure 52. Ideal motion vector fields superimposed on true images Figure 53. Estimated motion by M step using true images. (B = 0.10). Figure 54. Individually reconstructed images by PSMLEM and estimated motion by M step. (cutoff = 0.4 cycle/pixel; B = 0.05). Figure 55. Estimated motion and reconstructed images by RM. (a=0.10; B = 0.05). The global motion evaluation results are shown in Figure 56. The minimum global motion errors of the 3 methods are shown in Table 51. The numbers in the parentheses are the standard deviations of global motion error calculated from a 20image ensemble. Each method implemented 40 iterations. Table 51. Global motion evaluation results of the NCAT phantom Method Optimum conditions Minimum global motion error 1 B = 0.10. 267.7 2 cutoff = 0.4 cycle/pixel; 304.1 (11.3) B = 0.05. 3 a=0.10; 288.6 (5.8) B = 0.05. 57 340 330 S320 O 310 S300 S290 280 O 02 O 04 0.06 O0 OBO1 0 12 0.14 O 16 P A S500 700 S500 r~~ : 300 a 03 200 0.12 1 1 P 0.04 0.4 OO0 0.02 0.2 cutoff B C Figure 56. Global motion errors of three methods. A) global motion error vs. a for Method 1. B) global motion error vs. cutoff frequency and B for Method 2. C) global motion error vs. (a, B) for Method 3. 5.2.2 Regional Motion Evaluation A 2 X 2 subregion within the defect region in slice 7 (indicated by the square in Figure 57) was chosen for the regional motion evaluation. The results are shown in Table 52. The mean magnitude errors are in units of pixels, and mean angle errors are in units of degrees. The numbers in the parentheses are standard deviations of the errors. Overall the table indicates the improved motionestimation accuracy for the RM algorithm compared with the individual method. Figure 57. Chosen region to calculate regional motion error for the NCAT phantom Table 52. Regional motion evaluation results of the NCAT phantom Individual method RM Ideal mag. Mag. error Angle error Mag. error Angle error 1.56 0.84 (0.28) 20.2 (8.5) 0.52 (0.21) 18.6 (5.6) 5.3 Image Reconstruction Results To evaluate the imagereconstruction performance of the RM algorithm for the NCAT phantom study, we also compare the global image error (RMS), the regional image error (bias/variance), and the defect contrast of the RM images to those of the PSMLEM images over a range of cutoff frequencies. 5.3.1 Global Image Evaluation For PSMLEM algorithm, the mean RMS error of 2 frames changed with the cutoff frequency, as indicated in Figure 58 (A). The minimum RMS error (0.284) was achieved at approximately 0.6 cycles/pixel cutoff. For RM algorithm, a surface representing the change of the mean RMS error with the hyperparameters (a, a) is shown in Figure 58 (B). The minimum RMS error of RM (0.268) was achieved when a=0.13 and 8=0.07. Figure 59 shows the mean RMS error of the 2 frames versus the iteration number of RM outer loop. The mean RMS error decreased with iteration until approximately 40, 59 then stayed flat up to iteration 60. This demonstrates again the stability of the RM image with the outerloop iteration. 0.3 0.26. 2 0.3 .1 0. e 01 .0. O 0. .015 . ~~f0. 0. .2 0.2 * 0 ~ ~ 20 400 6050 Itertio nurnberc of RMM oute loo Figure 59. RMS errorvs.o ithe rato anume ofRM outger. loop. rrr s ct 532DfrqetCntras EvaPSluation err s (,B)fr Tocmpt tedfetcotatw ne irtdeemnete oml n dfc regions.~ Th eetrgosof2fae eebthcoe uigte iuaino h phato. W gneate te ora reinb oprngtevxlitnst ftetu images to a selected threshold. Figure 510 shows the normal and defect regions for 2 frames. White pixels present normal and gray pixels present defect. Most voxels in left myocardium were included into the normal or defect region. Figure 510. Normal and defect regions in the NCAT phantom. A) Frame 1. B) Frame2. As indicated in Figure 511i, PSMLEM using a higher cutoff frequency produced a higher contrast. The highest contrast of PSMLEM images (0.53) was achieved with the cutoff of 1.4 cycle/pixel. RM images had a lower contrast (0.5) compared with the highest contrast achieved by PSMLEM. However, RM images presented better tradeoff between the contrast and noise than PSMLEM. That is, RM had lower noise level while their contrasts were similar, and higher contrast while their noise levels were similar. Table 53 illustrates this tradeoff between the contrast and noise of RM and PSMLEM. We used the standard deviation of the voxel intensity in the normal myocardium region to represent the noise level. For RM, this standard deviation was 0.15. The PSMLEM image had comparable noise level at the cutoff of 0.4 cycle/pixel, but the contrast was substantially reduced to 0.43. At a higher cutoff of 0.7 cycle/pixel in which the contrast was comparable, the standard deviation was substantially higher (0.25). The images of the 3 cases are shown in Figure 512. 0.5 1 cutoff frequency Figure 511i. Contrast vs. cutoff frequency of PSMLEM Figure 512. The tradeoff between the contrast and noise level of the NCAT phantom. A) RM images. B) PSMLEM images with similar contrasts but higher noise levels. C) PSMLEM images with similar noise levels but lower contrasts. Table 53. The contrast and noise level of the NCAT phantom Method Contrast Std. Dev. RM 0.50 0.15 PSMLEM 0.50 0.25 (cutoff=0.7 cycle/pixel) PSMLEM 0.43 0.15 (cutoff=0.4 cycle/pixel) 5.3.3 Regional Image Evaluation Figure 513 shows the 2 subregions of interest in Frame 1 and Frame 2: one in the normal region and the other in the defect region. Figure 514 shows the tradeoff between intensity bias and variance at these subregions in Frame 1 and Frame 2. We investigated 5 sets of hyperparameters (see Table 54) for the RM algorithm. Set #1 (a=0. 13; 8=0.07) is the one that produces minimum RMS error. Set #2 increases a and keeps B unchanged, to investigate the influence of a larger value of a. Set #3, #4, #5 are to investigate the influence of a smaller a, a smaller B, and a larger B, respectively. These 5 sets of hyperparameters presented different bias/variance property, as indicated in Figure 514. For the subregion in the defect of Frame 1, all of the Hyve RM images presented close or slightly worse tradeoff between intensity bias and variance than PSMLEM. For the subregion in the defect of Frame 2 or the normal of Frame 1 and 2, all of the Hyve RM images presented better tradeoff between intensity bias and variance than PSMLEM. The one out of these 5 sets that achieved the best bias/variance tradeoff was the set #5, which grouped a medium a (0. 13) and a larger B(0. 15). It seems that the regional bias/variance property was benefited by the larger a which imposed more material property constraints and produced smaller and smoother motion Hield. Table 54. Five sets of hyperparameters for RM used for regional image evaluation Set #1 Set #2 Set #3 Set #4 Set #5 a 0.13 0.25 0.03 0.13 0.13 fl 0.07 0.07 0.07 0.01 0.15 PSMLEM \+ RM: set 1 + RM: set 2 RM: set 3 \x RM: set 4 RM: set 5 + 63 Figure 513. Subregions for regional image evaluation for the NCAT phantom. A) Frame 1. B) Frame 2. \ PSMLEM \+ RM: set 1 \+ RM: set 2 \RM: set 3 x RM: set 4 i~RM: set 5 ii+ 0.05 0.1 0.15 Variance 0.3 0.2 0.05 0.1 0.16 Variance 0.2 B 1.5 r 0 O o~ O 0.1 0.2 0.3 0.4 Variance 0.05 0.1 0.15 0.2 Variance Figure 514. The tradeoff between regional intensity bias and variance of the NCAT phantom. A) Subregion in the defect region of Frame 1. B) Subregion in the defect region of Frame 2. C) Subregion in the normal region of Frame 1. D) Subregion in the normal region of Frame 2. CHAPTER 6 OBSERVER EVALUATION OF IMAGE QUALITY 6.1 Methods 6.1.1 SingleSlice CHO Model Our study of observer evaluation of image quality for the RM algorithm started with a singleslice CHO (SCHO) model. A number of channel models have been presented since Myers et al. [45] added a channel mechanism to the ideal observer. Previous research [45, 49] revealed that the performance of the CHO model was rather insensitive to the details of the channel model, such as the bandpass cutoff frequencies and the geometric properties of the channels. Since Radially Symmetric Channel (RSC) model demonstrated excellent performance to predict human performance [47, 49, 56], we used this relatively simple but efficient channel model in our study. Our RSC model consisted of rotationally symmetric, nonoverlapping channels with cutoff frequencies of (0.03 125, 0.0625), (0.0625, 0. 125), (0. 125, 0.25), (0.25, 0.5) cycles/pixel. Choice of cutoff frequencies was somewhat arbitrary. Figure 61 shows the RSCmodel channels in the frequency domain. The channelization procedure was described in Equation 219 as v = Uf To compute the vector v from the image ffor a locationknown defect, we Fourier transformed each frequencydomain channel to a spatialdomain template centered at the defect location. The ith COmponent of the vector v was then calculated by integration of the pixelbypixel product of the image and the ith spatialdomain template. This is to avoid excessive Fourier transformation for the large number of images. Figure 61. The RSCmodel channels used for channelization in singleslice CHO model. (Note: Images are shown in the frequency domain.) Figure 62 shows the 2 defects in the geometric phantom and the defect in the NCAT phantom. Figure 63 shows the spatialdomain templates of the 4 RSCmodel channels used in our study. Three rows indicate the different defect locations in the two phantoms. A B Figure 62. Defect locations in the phantoms. A) Geometric phantom. B) NCAT phantom . Before implementing CHO to compute the test statistic ACHO We need to create a training set of obj ects in two states (SP and SA), and estimate the scatter matrix S2 and the mean vectors v, and v The size of the training set is usually determined by the channel number. As previous studies [54, 55] presented, for an Nchannel model, the size of the training set for each state (SP or SA) should be larger than N2. Otherwise the inverse of S2 Will HOt exist. In our study, we used 4channel RSC model, so the size of the training set should at least be 16. We used 40 to ensure the existence of S2. To compare RM with PSMLEM, we generated 280 images (140 SP and 140 SA) for each algorithm. We used 80 out of them (40 SP and 40 SA) as the training set. After generated the scatter matrix S2 and the mean vector v, and va, we applied the CHO discriminant vector on the other 200 images (100 SP and 100 SA) to compute the SCHO detectability index (SCHOd,) To acquire the standard deviation of the detectability index, we divided the 280 images arbitrarily to 7 groups, each containing 40 images (20 SP and 20 SA). Each time we selected 2 groups of images to be used as the training set, then applied the resulting CHO discriminant vector to the remaining 5 groups to estimate an SCHOd,. By randomly selecting 2 out of the 7 groups as the training set, we estimated 21 different values of SCHO,' s. Thus we are able to compute the mean and standard deviation of SCHOd,. 6.1.2 MultiFrame MultiSlice CHOHO Model Since a defect in a 3D obj ect may persist from slice to slice, in practice the physicians generally make their diagnostic decisions by observing multiple slices of images, and even multiple frames in the case of gated emission tomography imaging. Chen et al. [65] presented a multislice multiview CHOHO model for ungated SPECT myocardial perfusion imaging. We modified it to a multiframe multislice CHOHO (MMCHO) model for gated myocardial perfusion imaging. The MMCHO model involves 2 steps. In the first CHO step we apply a CHO model to each slice of all frames to assess the probability that the defect is present in that slice. The result of the first step for each multiframe multislice image is a vector of test statistics of which the element is an individual test statistic ASHO, Where i=1, 2 ,..., I (I is the number of slices) is the slice index and j=1, 2, ..., J (J is the number of frames) is the frame index. In the second HO step, we apply an HO model to the ensemble of vectors of A O(SP and SA) to compute a final MMCHO detectability index (MMCHOd,). Figure 64 shows the scheme for this 2step procedure. Figure 63. The RSCmodel channels represented as spatialdomain templates. The channels are Fourier transformed and shifted to the location centers of the defects. A) Channels for Defect 1 in the geometric phantom. B) Channels for Defect 2 in the geometric phantom. C) Channels for the defect in the NCAT phantom . In our MMCHO study, we considered 2 frames, each containing 3 chosen slices. The first CHO step applied the SCHO model to all the 3 slices of the 2 frames, and generated vectors containing 6 test statistics, A HO (i=1, 2; j=1, 2, 3). The second HO step used 60 teststatistic vectors as the training set to determine the HO discriminant vector, and then applied the discriminant to the remaining teststatistic vectors to compute the MMCHOd,. The mean and standard deviation of MMCHOd, were computed by the same procedure used in the SCHO model. frames /SP ensemble SA \ ensemble 3 slices 3 slices / SP ensemble ensemble S~A CHO HO O SP~ CHO HO Multiframe multislice CHOHO detectability index (MMCHOder) Figure 64. Twostep procedure of the multiframe multislice CHOHO model. CHO channels Step 1: CHO Step 2: HO 6.2 Results of Geometric Phantom Study Since the 2 defects in the geometric phantom were both centered in slice #15 of the long axis, we picked that slice and its 2 adj acent long axis slices, #14 and #16, to implement the SCHO and MMCHO study. The 3 images slices of 2 frames reconstructed by PSMLEM and RM are shown in Figure 65 and Figure 66, respectively. In each figure, the image slices for each frame are grouped in the order of #14, #15, and #16. The cutoff frequency for PSMLEM was 0.3 cycle/pixel. The RM algorithm used the hyperparameters (a=0.02, 8=0.04), which produced minimum RMS image error. The PSMLEM images present more smoothness than the RM images, particularly at the myocardium boundaries. 6.2.1 SingleSlice CHO Results We performed SCHO model to the 2 defects in the 3 longaxis slices of the 2 frames. For PSMLEM, we used a range of cutoff frequency from 0. 1 to 1.0 cycle/pixel and found the minimum SCHOd,. The results of PSMLEM and RM are shown in Table 61. Standard deviations of the SCHOd, are given in parentheses. The RM images had significantly improved SCHOd, in all of the cases except for the defect 1 in slice #15 of Frame 2. The improvement of RM for one slice was evidently different from another, while the SCHO performance of the slices differed, too. This also motivated us to develop a multiframe multislice observer model that produces a single detectability index for a series of gated myocardial ET images. 6.2.2 MultiFrame MultiSlice CHOHO Results We applied the MMCHO model to the 2 frames and 3 slices shown in Figure 65 and Figure 66. The results are shown in Table 62. The numbers in parentheses are the Figure 65. SP images of the geometric phantom reconstructed by PSMLEM. A) Frame 1. B) Frame 2. Figure 66. SP images of the geometric phantom reconstructed by RM. A) Frame 1. B) Frame 2. ~PSMLEM standard deviations of MMCHOd,. The results show that the RM images presented significantly higher MMCHO detectability than the PSMLEM images. Figure 67 shows how the MMCHOd, of the PSMLEM images changed with the cutoff frequency. For both defects, the maximum detectability indices of PSMLEM were achieved with the cutoff of 0.4 cycle/pixel. The MMCHO,' s of the RM images are also plotted. Table 61. SCHO detectability index results of the geometric phantom SCHOd, for Defect 1 Frame 1 Frame 2 6.98 (0.22) 6.92 (0.21) 9.26 (0.28) 8.85 (0.42) 8.38 (0.25) 9.06 (0.23) 8.95 (0.32) 8.03 (0.16) 9.07 (0.35) 8.58 (0.36) 12.89 (0.51) 11.87 (0.23) 25 SCHOd, for Defect 2 Frame 1 Frame 2 7.53 (0.13) 8.89 (0.24) 11.52 (0.63) 11.93 (0.46) 8.27 (0.19) 9.46 (0.28) 10.69 (0.30) 9.89 (0.23) 9.28 (0.35) 8.85 (0.35) 13.24 (0.49) 12.62 (0.41) Slice Algorithm 14 PSMLEM RM 15 PSMLEM RM 16 PSMLEM RM 25 x a, ri 20 i 15 n rd ~ io a, 0 5 I y L. r 0 0. 2 0. 4 0. 6 0. 8 1 1. 2 cutoff frequency (cycle,'pixel) 0 0.2 0.4 0.6 0.8 1 1.2 cutoff frequency (cycle,'pixel) Figure 67. MMCHO detectability index results of the geometric phantom. A) Defect 1. B) Defect 2. (Note: For the purpose of comparison, RM result is also represented as a point with error bar in each plot disregarding the cutoff frequency.) Table 62. MMCHO detectability index results of the geometric phantom MMCHOd, for Defect 1 MMCHOd, for Defect 2 PSMLEM 14.23 (0.25) 15.71 (0.23) RM 19.25 (1.13) 20.81 (0.91) 6.3 Results of NCAT Phantom Study For the NCAT phantom, we selected long axis slice #7, #8 and #9. The defect is centered in slice #7 and #8. Figure 68 shows the true image slices of two frames. Figure 69 shows the images reconstructed by the PSMLEM algorithm using the cutoff frequency of 0.3 cycle/pixel. Figure 610 shows the images reconstructed by the RM algorithm using the hyperparameters of (a = 0.03, a = 0.07). In each figure, the image slices for each frame are grouped in the order of #7, #8, and #9. 6.3.1 SingleSlice CHO Results The results of singleslice CHO performance for the PSMLEM and RM images are shown in Table 63. Again, the minimum SCHOda for PSMLEM was found by using a range of cutoff frequencies from 0. 1 to 1.0 cycle/pixel. Figure 68. SP true images of the NCAT phantom. A) Frame 1. B) Frame 2. Figure 69. SP images of the NCAT phantom reconstructed by PSMLEM. A) Frame 1. B) Frame 2 Figure 610. SP images of the NCAT phantom reconstructed by RM. A) Frame 1. B) Frame 2 Table 63. SCHO detectability index results of the NCAT phantom Slice Algorithm SCHOd, for Frame 1 SCHOd, for Frame 2 7 PSMLEM 11.85(0.23) 12.81(0.64) RM 20.35(0.80) 19.80(1.11) 8 PSMLEM 10.71(0.30) 13.01(0.81) RM 18.94(0.69) 18.45(0.45) 9 PSMLEM 8.01(0.26) 9.05(0.46) RM 8.32(0.35) 16.66(0.49) 6.3.2 MultiFrame MultiSlice CHOHO Results To investigate the algorithm performance for a lower defect contrast, we generated another simulated NCAT phantom which had a defect with 25% decreased intensity. We performed the MMCHO observer to both of the two phantoms (50% defect and 25% defect). The results are all shown in Table 64. It is evident that for both 50% and 25% defect phantoms, the RM images presented significantly higher MMCHO detectability than the PSMLEM images. Figure 611 shows how the MMCHOd, of the PSMLEM images changed with the cutoff frequency. Table 64. MMCHO detectability index results of the NCAT phantom MMCHOd, for 50% defect MMCHOd, for 25% defect PSMLEM 16.98(0.77) 8.14 (0.36) RM 33.41(0.66) 14.59 (0.52) 6.4 MultiFrame MultiSlice MultiView CHOHO Study Since a 3D defect in practice may present different shape, size, and probably different detectability in different axis views, it is meaningful and feasible to expand the multiframe multislice CHOHO model to a multiframe multislice multiview CHOHO (MMMCHO) model. There are also 2 steps involved in the MMMCHO model. In the first CHO step, the SCHO model is applied to each slice of all views and all frames. The first step produces a teststatistic vector with elements as AfHO,, Which is from the ith slice in the kth VieW Of the jth frame. The second HO step takes the resulting teststatistic vectors to compute a Einal MMMCHO detectability index (MMMCHO da). 40 16 35 ~ PSMLEM __APSMLEM RM HR 3012 25 10 10  0 0.2 0.4 0.6 0.8 1 20 0. 2 0. 4 0. 6 0. 8 1 1. 2 cutoff frequency (cycle/pixel) Acutoff frequency (cycle/pixel) B Figure 61 1. MMCHO detectability index results of the NCAT phantom. A) 50% defect. B) 25% defect. In our MMMCHO study, we considered 2 frames, 2 views (a long axis and a short axis), and 4 slices. The 16 images of a SP realization are showed in Figure 612, 613 and 614, presenting the PSMLEM algorithm (cutoff = 0.4 cycle/pixel), and the RM algorithm (a = 0.03 and 0.13), respectively. In each Eigure, the 1st row is a long axis view of Frame 1, the 2nd TOW iS a short axis view of Frame 1, the 3rd TOW iS a, long axis view of Frame 2 and the 4th TOW iS a short axis view of Frame 2. Because our previous studies showed that the RM images were more sensitive to a than to a, we implemented the algorithm using 2 different values: 0.03 and 0. 13 for a, while B was kept at 0.07. Since the teststatistic vector produced in the first CHO step has 16 elements, the size of the training set for each state (SP or SA) in the second HO step should be larger than 162 = 256. In our study, we used 350 images of each state as the training set. Thus, we generated 800 images (400 SP and 400 SA) for each algorithm. The first CHO step used 80 images (40 SP and 40 SA) as the training set. The second HO step used 700 sets of teststatistic vectors (3 50 SP and 3 50 SA) produced in the first step as the training set for the HO model. The HO discriminant vector was then applied to the remaining 100 teststatistic vectors (50 SP and 50 SA) to estimate the MMMCHO detectability index (MMMCHOd,). The MMMCHOd, versus cutoff frequency of PSMLEM is shown in Figure 615. The highest MMMCHOd, of PSMLEM was achieved with the cutoff of 0.4 cycle/pixel. The MMMCHO,' s and their standard deviations of the RM images and the PSMLEM images are shown in Table 65. The detectability of RM algorithm was affected by the value of a. A smaller value of a produced smoother images and achieved higher defect detectability. Future studies will investigate the MMMCHOd, as a function of a and B. Figure 612. Multiframe multislice multiview SP images reconstructed by PSMLEM. Figure 613. Multiframe multislice multiview SP images reconstructed by RM (a= 0.03). Figure 614. Multiframe multislice multiview SP images reconstructed by RM (a= 0.13). 78 20 A PSMLEM x 18 .~16 r 4RM 2 0 0. 2 0. 4 0. 6 0. 8 1 1. 2 cutof ffrequency (cycle/pixel) Figure 615. MMMCHO detectability index results of the NCAT phantom. Table 65. MMMCHO detectability index results of the NCAT phantom. MMMCHOda PSMLEM 9.93 (0.78) RM (a= 0.03) 16.55 (0.76) RM (a= 0.13) 13.58 (0.83) In summary, we presented a singleslice CHO model and two multiimage CHOHO models for observer evaluation of the gated myocardial ET images. Results of all the 3 numerical observers demonstrated that the RM algorithm produced images with higher detectability, relative to the PSMLEM algorithm. CHAPTER 7 CONCLUSIONS AND FUTURE WORK The obj ective of this study was to implement a new method that simultaneously reconstructed gated images and estimated cardiac motion for myocardial ET, and evaluate the performance of this method on a simulated geometric phantom and an NCAT cardiac phantom. The theoretical derivation of the optimization RM algorithm dealing with the nonrigid motion vector field and the emission images in two time frames of gated myocardial ET has been described. The RM algorithm may be viewed as a penalized version of Klein' s motionestimation algorithm in which the data likelihood forms an additional constraint. It may also be viewed as penalized likelihood reconstruction in which image smoothing is achieved by assuming the motion is a smooth deformation between frames, and that the intensity of the material points changes little between frames. The reconstructed images were compared with an independent reconstruction method consisting of high iterations of the MLEM algorithm applied to each proj section dataset followed by smoothing with the Hann filter, an algorithm often called PSMLEM. Our results showed that the RM images presented lower global image error and regional image error, relative to PSMLEM. For defect contrast study, the geometric phantom and the NCAT phantom were simulated one or more defective regions with 50% decrease in voxel intensity. Our results demonstrated that the RM images achieved substantial improvement in the contrastnoise tradeoff over the PSMLEM images. The RM images had a substantially higher defect contrast compared to the PSMLEM images except for the highest cutoff frequency PSMLEM images. However, at this high cutoff frequency, the PSMLEM images had substantially higher noise than the RM images. The RM estimated motion vector field was compared with that obtained from applying Klein's method to the PSMLEM images. Our method outperformed that method both in how well the resulting deformed image matched the reference image and in regional comparisons between the resulting motion fields and the ideal motion field. For the observer evaluation of image quality on specific detection task, we used the CHO model because it had shown good correlation to human observer. Based on the fundamental singleslice CHO model, we proposed a multiimage CHOHO model for gated myocardial ET. Results with both the singleslice CHO model and the multiimage CHOHO model demonstrated that the RM images had significantly improved detectability than the PSMLEM images. Our future work includes 1. Multiframe RM algorithm: Current RM algorithm deals with two time frames due to the restriction of huge amounts of computation. With the help of more powerful computers, multiframe processing will expand the current twoframe algorithm and is expected to present even higher quality gated images. 2. Human observer study: Though numerical observer is an efficient path for observer assessment of image quality, currently it is not a substitute of human observer. Human observer study is still necessary for further observer assessment study. 3. Spatial resolution study of RM: The regularization of MLEM by combing the likelihood with a penalty often results in position and image dependent spatial resolution and bias. More study about the spatial resolution and noise property of the RM algorithm is required. 4. More concerns such as attenuation, scatter, and randoms: Additional study is required to evaluate how these issues affect the performance of the RM algorithm. It plays an important role to ensure the advantage of RM in clinical application. In conclusion, the new RM algorithm has demonstrated the potential to improve both image quality and motionestimation accuracy for gated myocardial ET, relative to the independent motionestimation and imagereconstruction algorithms. More studies are required to investigate the performance of the RM algorithm in clinical application. APPENDIX A NUMERICAL ANALYSIS FOR MOTION ESTIMATION Numerical analysis is the mathematical method that uses numerical approximations to obtain numerical answers to the problems that do not have analytical solutions. The foundation of differential and integral calculus is based upon limits of approximations. The definition of the firstorder derivative is a difference quotient: du .u(x + x) u(x) = hm (A1) Therefore, there are various approximations for the derivative. Each of them has its own benefits and disadvantages. Some examples of the approximations are du u(x + &x) u(x) Forward difference: 4 ( A 2) dx~ hx du u(x) u(x &x) Backward difference: 4 (A3) dx~ bx 1 1 Central difference: 22 or (A4) 2hx and the central difference approximation for the secondorder derivative: d2u u(x + &) 2(x)+ u(x h) (A5) dx~C ( )2 The central difference requires information in front of, and behind the location being approximated, whereas the forward and backward differences use the information on only one side of the location being approximated. In our study, we use the central difference because of its better accuracy. To obtain the numerical solution of the motionestimation problem, we rewrite the Equations 360, 361, 362 as p[(A + 2~)u_ + p~(uY + uzz) + (I C P)(xy xzW ) 2m Vf2(k)( k) 2(,k) k Dx (A6) = 2(fi(k) ( f2(k) k) k) 2(k1k) ) 2(k) k) P$(A + 2pu)vy, + pu(v, + v, ) +(Al + pu)(uxy+ yz) . 2m1 Vf2(k) k) k)1 '' (A7) = 2 fi*(k) 2(k) + ik)) Ilk)) V2(k) ( lik)) k)+II ay $(A+ 2~w,+ pw +w )+ ( + )(u yz m Vf2k) (k)) 2(k) +(k)) 8z (A8) = f(k 2() (k) k) 2(k (k)1 2(k) (k)I 8z where k denotes the kth outer loop, I denotes the lth Mstep inner loop, and m k) is current motion estimation. The solution m of Equations A6, A7, A8 is the motion estimation for the (Il+)th Mstep inner loop. Equations A6, A7, A8 are linear to m and can be written in the matrix form as A M = B . Assuming the 3D image has dimensions ofIX JX K, the total voxel number is 11, 127 1(3N) N=IX JXK. Then, A= a21, 227 2(3N) iS a 3NX 3Nmatrix, a(3N)1, (3N)2, > (3N)(3N) M = [u1, u,.,27"' N> 1 2 N>,11 1'7 W2,l' NT iS the 3NX 1 motion vector field to estimate, and B = b,, b2,..., bN b ,bN+2,.,bN 2+,bN2.. 3 Sa3X1vco determined by image frames and current motion estimation. In numerical approach, the voxel coordinate components of r, (x, y, z) are represented as the discrete indices (i, j, k), respectively. Using the central difference approximation, the first and second differential of the motion component u at the location (i, j, k) can be expressed by the linear combination of the motion at its adjacent voxel locations. Substituting u with v and w gives the first and second differential of v and w. u u+1,],k 1Ik, +1,k 1, ,],k+1 ,, A 9 x2 2 a2 uxx = (uz+,l,,k ],k ,, 21], l+,,,k 2ur], 21] (A10) uW = (u",J+1,k U,,,k (U,,,k 2,J1,k Ui,J+1,k 2uI,,,k + 2,J1,k (A11) uzz = (ul,],k+1 U, ,k (U,,,k 2,,,k1 Ul,],k+1 2uI,,,k + 2,,,k1 (A12) u (u+1,+1, 2,J+,k +,J,k 1,J1,k(A13) v 4 u ( rJ1k1 2J1k+ J1k1 ,, (A14) u,,~ (A15) For each voxel (i, j, k), Equations A6, A7, A8 generate 3 rows of entries in the matrix A. The coefficients of (uz, u27"'>, N) for voxel (i, j, k) from Equation A6 can be filled into a chart as follows, where i1, i, i+1, j1,j +, k1, k, k+1 are voxel location indices in 3 dimensions, and the values in the chart are the coefficients for the motion component u at the voxel located by the indices. For example, the coefficients for uz+1,],k is (il+2u)/7. i1 0 0 0 0 (a +2pu)P 0 0 0 (2l+ 8pu)P i pf 0 upf df2 (k) 2 pf 0 up 2( )2 8x i+1 0 O2) l+pu i+pu 0 i1 0 0 0 P 0 p 0 0 4 4 df(k) df(k) 0 i 0 0 00 2(2 2 )00 0 i+1 l+pu i+pu 0 0 00 P 0 p 0 4 4 . O p o 0 f2(i 2(k) 2(k i+1 l+pu i+pu 0 0 p 0 0 0 0 0 P 4 4 j1 j j+1 j1 J j+1 Similarly, the coefficients of (v,, v27"'>, N) in Equation A6 for (i, j, k) can be expressed using the following chart: j1 j j+1 j1 j j+1 The coefficients of (w,,w27"' N,) in Equation A6 for voxel (i, j, k) j1 j j+1 j1 j j+1 The coefficients of (uz, u27"'>, N) for Equation A7 for voxel (i, j, k) are 4 4 (3/ (k) (f(k) 0 0 i 0 0 0 0 2( 2 O 8x S i+ o p o~ 0 0 4 4 k1 k k+1 The coefficients of (v, ,v27"'>v N) for Equation A7 for voxel (i, j, k) are j1 j j+1 j1 j j+1 j1 j j+1 i1 0 0 0 0 ##p 0 0 0 0 (2l+ 8pu)P pup (a + 2 p) f(k (a + 2 p) f i 0 0 of 2o 0 # 2( )2 i+1 k1 k k+1 The coefficients of (w ,w27"' N,) for Equation A7 for voxel (i, j, k) are j1 j j+1 j1 j j+1 j1 j j+1 i1 0 0 0 0 0 0 0 0 il+u p A+ p af2 (k) d2(k) lu c az2( )(/ 4 i+1 0 j1 j j+1 j1 j j+1 The coefficients of (u, ,u27"'>u N) for Equation A8 for voxel (i, j, k) are i Op o o o p0 4 4 8f~:k) ~S(k) 0 i 0 0 0 0 ~ 2(2 2) 0 0 0 i+1 l+pu i+pu 0 4 4 k1 k k+1 The coefficients of (v,, v27"'>, N) for Equation A8 for voxel (i, j, k) are j1 j +1 j1 j +1 j1 j +1 i1 0 0 0 0 0 0 0 0 0 i +u p f, (k) df(k) lu 2+#2( 2 4 y 8 4 0 4 i+1 k1 k k+1 The coefficients of (w,,w27"' N,) for Equation A8 for voxel (i, j, k) are j1 j +1 j1 j +1 j1 j +1 i1 0 0 0 0 ##p 0 0 0 0 (2il + 8pu)P 2( )2 8z i+1 # j1 j A numerical image in practice has certain boundaries. Numerical analysis methods have to consider the physical processes in the boundaries. Different boundary conditions may cause very different results. Because of the ideally compressible model for the space outside myocardium, it is reasonable to assume the motion at the image j+1 j1 boundaries to be approximately zero. Thus in our study, we used "fixed boundary conditions", in which m(xb, yb, Zb) = 0 where (xb, yb, Zb) is a boundary voxel. The calculation of f 'k(k) (kIl)) and its derivatives is not straightforward since m k) iS a floating number while r is integer indices (i, j, k). Two approaches to approximate f,(k)( ( mk)) and its derivatives are described and compared here. For notation convenience, we write f,(k) (~mk)) as f (r +m) in following text. 1. Taylor approximation approach Let m = m +Sm where m =(u, v,w) and (u, v, w) are the closest integer approximation to (u, v, w); 6m = m m is the difference vector. Then from Taylor approximation we obtain (A17) (rrm)= (r )m)+& (r +m1)+ Sv (r +m1)+ &v (r +m1) (A18) (r m =(r + m1) + At (r + m1) + Sv (r + m1) + Av (r + m1) (A19) The differential of f (r + m1) can be approximated using central difference. d, f,(i + u+1, j +v, k +w) f,(i + u1, j +v, k +w) (r +m) =(A20) Dx 2 