|UFDC Home||myUFDC Home | Help|
This item has the following downloads:
SIMULTANEOUS RECONSTRUCTION AND 3D MOTION
ESTIMATION FOR GATED MYOCARDIAL EMISSION TOMOGRAPHY
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
This dissertation is dedicated to my parents and my wife
It is impossible to adequately thank my long-time 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 medical-imaging 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
April-Lane Derfinyak and Laura Studstill) for their help, kindness, and patience. I thank
Rami, Uday, and Ping-Fang 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
ACKNOWLEDGMENT S ................. ................. iii........ ....
LIST OF TABLES ................ ..............vii .......... ....
LI ST OF FIGURE S ................. ................. viii............
AB STRAC T ................ .............. xi
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 Euler-Lagrange 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 Bi-Value 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 Single-Slice CHO Model ............... ... .............. ....................6
6. 1.2 Multi-Frame Multi-Slice CHO-HO Model ................. .....................66
6.2 Results of Geometric Phantom Study .............. ...............69....
6.2.1 Single-Slice CHO Results.................. ...... ...........6
6.2.2 Multi-Frame Multi-Slice CHO-HO Results .............. .....................6
6.3 Results of NCAT Phantom Study .............. ...............72....
6.3.1 Single-Slice CHO Results.................. ...... ...........7
6.3.2 Multi-Frame Multi-Slice CHO-HO Results .............. .....................7
6.4 Multi-Frame Multi-Slice Multi-View CHO-HO 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
2-1. Conditional probabilities of a defect detection task .............. .....................14
4-1. Ellipsoid semi-axes lengths for source obj ect ................. ...............33.............
4-2. Poisson ratio and Lame constants ................. ...............37...............
4-3. Global motion evaluation results of the geometric phantom ................. ................ .44
4-4. Regional motion evaluation results of the geometric phantom ........._..... ...............45
5-1. Global motion evaluation results of the NCAT phantom ................. ............... .....56
5-2. Regional motion evaluation results of the NCAT phantom ................ ................. .58
5-3. The contrast and noise level of the NCAT phantom ................. ................ ...._.62
5-4. Five sets of hyper-parameters for RM used for regional image evaluation ...............62
6-1. SCHO detectability index results of the geometric phantom .............. ...................71
6-2. MMCHO detectability index results of the geometric phantom .............. .... .........._.71
6-3. SCHO detectability index results of the NCAT phantom ................. .............. ....74
6-4. MMCHO detectability index results of the NCAT phantom ................ ................. 74
6-5. MMMCHO detectability index results of the NCAT phantom ................. ...............78
LIST OF FIGURES
2-1. Anatomic orientation of the heart. ............ .....___ ...............6
2-2. Principle of data acquisition and geometric considerations for SPECT.............._._.....7
2-3. Coordinate systems transformation. ........._._. ......_ ....__ ......._.........8
2-4. Images of MLEM and PS-MLEM ................. ...............11...............
2-5. Motion vector representing the displacement of voxels ................. ............. .......12
2-6. Klein's iterative motion-estimation optimization method............. .. ........._ ....13
3-1. Scheme of the RM al gorithm. ................ ........... ........ ......... ...........27
4-1. Source objects of the geometric phantom............... ...............32
4-2. Simulation steps of the geometric phantom .............. ...............33....
4-3. Ideal motion fields of the geometric phantom in 3 short-axis slices. ................... ......3 5
4-4. Collapse high-resolution proj section to low-resolution proj section ................... ...........36
4-5. Segmentation for bi-value Lame constants model. ............. .....................3
4-6. Convergence of the RM algorithm ................ ...............39...............
4-7. Convergence of R step............... ...............40..
4-8. Convergence of M step............... ...............40..
4-9. Global motion error vs. fl for several Poisson ratios. ............. .....................4
4-10. Ideal motion vector fields superimposed on true images. ............. .....................42
4-11. Estimated motion by M step using true images ................. ......... ................43
4-12. Individually reconstructed images by PS-MLEM and estimated motion by M step.43
4-13. Estimated motion and reconstructed images by RM...................... ...............4
4-14. Chosen region to calculate regional motion error for the geometric phantom.........45
4-15. RMS errors of the RM and PS-MLEM images. .......... ................ ................47
4-16. Hyper-parameter effects on the RM images ................. ............... ......... ...48
4-17. RMS error of the RM images as a function of a and B............_.._. .....................48
4-18. Defect contrast of the PS-MLEM and RM images ................. .......................49
4-19. The trade-off between the contrast and noise level of the geometric phantom........5 1
4-20. Sub-regions for regional image evaluation for the geometric phantom. ..................51
4-21. The trade-off between regional intensity bias and variance of the geometric
phantom ..........._._ ....._. ._ ...............52.....
5-1. Long- and short- axis views of the NCAT cardiac phantom. ................ ................ .54
5-2. Ideal motion vector fields superimposed on true images .............. .....................5
5-3. Estimated motion by M step using true images............... ...............55.
5-4. Individually reconstructed images by PS-MLEM and estimated motion by M step..56
5-5. Estimated motion and reconstructed images by RM .................... ...............5
5-6. Global motion errors of three methods............... ...............57
5-7. Chosen region to calculate regional motion error for the NCAT phantom ................58
5-8. RMS errors of the RM and PS-MLEM images. ................ ................. ..........59
5-9. RMS error vs. iteration number of RM outer loop. .................. ................5
5-10. Normal and defect regions in the NCAT phantom ................. ................ ...._.60
5-1 1. Contrast vs. cut-off frequency of PS-MLEM ................ ...............61.............
5-12. The trade-off between the contrast and noise level of the NCAT phantom ............61
5-13. Sub-regions for regional image evaluation for the NCAT phantom. .......................63
5-14. The trade-off between regional intensity bias and variance of the NCAT phantom..63
6-1. The RSC-model channels used for channelization in single-slice CHO model.........65
6-2. Defect locations in the phantoms. ...._._._.. .... ..__.. ...............65..
6-3. The RSC-model channels represented as spatial-domain templates.. ......................67
6-4. Two-step procedure of the multi-frame multi-slice CHO-HO model. .................. .....68
6-5. SP images of the geometric phantom reconstructed by PS-MLEM. .............................70
6-6. SP images of the geometric phantom reconstructed by RM. ............. ...................70
6-7. MMCHO detectability index results of the geometric phantom. ............. ..... ........._.71
6-8. SP true images of the NCAT phantom. ................ ...............72..............
6-9. SP images of the NCAT phantom reconstructed by PS-MLEM ........._..... ..............73
6-10. SP images of the NCAT phantom reconstructed by RM ................. ................ ...73
6-11i. MMCHO detectability index results of the NCAT phantom. ................ ...............75
6-12. Multi-frame multi-slice multi-view SP images reconstructed by PS-MLEM.._......76
6-13. Multi-frame multi-slice multi-view SP images reconstructed by RM .....................77
6-14. Multi-frame multi-slice multi-view SP images reconstructed by RM .....................77
6-15. MMMCHO detectability index results of the NCAT phantom. .............. ..... ........._.78
A-1. 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
Chair: David Gilland
Maj or Department: Biomedical Engineering
A new method for simultaneous image reconstruction and wall-motion 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 image-reconstruction algorithms are
very noisy due to the reduced counts of detected gamma rays.
Non-rigid 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 voxel-intensity correspondence between two image frames. Previous
research incorporated a deformable elastic material model into a motion-estimation
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 motion-estimation accuracy and the
image quality of individual frames, relative to the independent motion-estimation and
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 non-uniform-rational-B-splines-based 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 image-reconstruction
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 X-ray examination, computed tomography
(CT), and magnetic resonance imaging (MRI) .
Tomography achieves detailed organ studies by separating a three-dimensional (3D)
obj ect into a stack of two-dimensional (2D) cut sections or slices. Unlike the X-ray CT
that detects transmitted X-ray photons, tomography in nuclear medicine detects emitted
gamma rays and therefore is referred to as emission tomography, or ET in short .
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 .
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 [6-9]. 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 image-reconstruction
algorithms usually are very noisy [4, 6].
Previous researchers [10-13] presented several reconstruction methods to improve
the quality of the gated image frames. Lalush et al.  considered the time sequence of
gated frames as a four-dimensional (4D) image and reconstructed it using a penalized
maximum-likelihood (ML) technique with space-time Gibbs priors to ensure smoothness
in individual image frames and between frames. Narayanan et al. [1l], Wernick et al. 
and Brankov et al.  used an alternative technique based on principle component
analysis, in which the time sequence of gated datasets were Karhunen-Loeve (KL)
transformed and then rapidly constructed frame-by-frame. The gated image frames were
finally obtained by applying the inverse KL transform.
Other previous studies [14-16] 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 less-noise images will in turn produce more accurate
Based on the 3D motion-estimation method developed by Klein et al. [15, 16] and
the conventional maximum-likelihood expectation-maximization (MLEM) algorithm [17,
18], we proposed a new estimation method for gated myocardial ET: the simultaneous
image-reconstruction and 3D motion-estimation 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  used the Polak-Ribiere 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 two-step iterative algorithm (RM) that
alternately updated the gated images and the motion vector field, and guaranteed
non-negativity in the images. This RM algorithm and some results were partly reported
by Cao et al. .
We evaluated the motion-estimation accuracy and image quality provided by the
simultaneous method using a simulated geometric phantom and a realistic
Non-uniform-rational-B-splines-based 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 wall-motion information in a
single imaging procedure. Our simultaneous estimation method promises to improve
motion-estimation accuracy and individual gated image quality simultaneously, which
will improve the accuracy in diagnosing cardiac diseases.
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 . 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 Tc-99m). 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 Tc-99m radioisotope
decays, and a fraction of these photons are able to penetrate the surrounding body tissue
and reach the SPECT gamma camera . 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
radioactive-pharmaceutical 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 .
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 . Figure 2-1
shows anatomic orientation of the heart.
Figure 2-1. Anatomic orientation of the heart.
2.2 Projection and Image Reconstruction
2.2.1 Radon Transform
Figure 2-2 shows the 2D data acquisition scheme in SPECT using a parallel
collimator . 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.
Figure 2-2. Principle of data acquisition and geometric considerations for SPECT.
Mathematically, the proj section operator can be described by the Radon transform
. 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
= 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 2-1 involved a
coordinate system transformation, which is illustrated in Figure 2-3. 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 forward-projection operator, and
h ()is a weight function that represents the probability that a photon is emitted from r
and detected at (s, 6)
Figure 2-3. 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 forward-projection 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 forward-proj 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 image-reconstruction algorithm adjusts the image
estimate iteratively, and eventually reconstructs an image . The iterative algorithms
differ in the comparison method and the correction applied to create the new estimate.
Equation 2-4 is based on the assumption of the ideal noise-free 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
well-known log-likelihood function can be written as
L(gi, f ) = [gig(i)lgHf(i)) Hf(i)] (2-5)
which measures the likelihood that the image produces the sinogram g.
Maximizing the log-likelihood function leads to finding the most likely image fto
produce the sinogram g. Maximum-likelihood expectation-maximization (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 image-updating 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 . The convergence, however, is very slow and the images appear
unacceptably noisy before it converges . To smooth the noisy image reconstructed by
MLEM, a digital post-processing filter is usually applied to the image reconstructed after
many iterations of MLEM. This Post-Smoothed MLEM (PS-MLEM) algorithm has
demonstrated excellent results in previous studies [26, 27]. The most important parameter
in the PS-MLEM algorithm is the cut-off frequency of the filter, which determines the
smoothness of resulting images. In our study, we use 100 iterations MLEM and
low-pass Hann filter because of its simplicity. Figure 2-4 (A) shows a sample image
reconstructed by MLEM with 100 iterations. Figure 2-4 (B) demonstrates that the
smoothness of the image reconstructed by PS-MLEM is influenced significantly by the
amount of the post-smoothing filter. Nyquist frequency (0.5 cycle/pixel) is usually used
as a feasible cut-off. Because of this, when we compare our algorithm to PS-MLEM
algorithm, we usually use a range of cut-off frequencies (from 0. 1 to 1.0 cycle/pixel) for
Figure 2-4. Images of MLEM and PS-MLEM. A) 100-iteration MLEM image B)
PS-MLEM image using Hann filter with various cut-off frequencies. From 1
to 10, the cut-off frequency was 0. 1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0
2.3 Myocardium Motion Estimation
A variety of methods have been proposed recently to estimate and represent the
cardiac motion [14-16, 28-32]. 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 voxel-based 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 voxel-based intensity scalar Hield.
Klein's motion-estimation method  involves minimizing a two-component
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 2-5, 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. 
Framel ,- Frame2:
r r> mr
Figure 2-5. 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 (2-7)
and the strain energy Es(m) is
Es (m) = fIx +", Vy +W X dr+pus+
where ii, pu are elastic weighting parameters called Lame constants that reflect the degree
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
hyper-parameter 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
ipp C u, + ", + .,,) u,, .L +v, i~+ ul +w, /",= l+ 3,-N~
]=1 L =1
Klein' s iterative optimization method  to minimize the entire obj ective function
E consists of steps shown in Figure 2-6.
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 2-6. Klein's iterative motion-estimation optimization method
2.4 Image Observer Evaluation
2.4.1 Human Observer Study
Besides the quantitative assessment such as the signal-to-noise 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 two-state 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 .
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
2-1 shows the conditional probabilities of a certain response given the occurrence of a
certain stimulus .
Table 2-1. 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 2-1 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 . 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  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 [35-37].
2.4.2 Numerical Observer: Background
The human observer tests are very time-consuming 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) , defined by
22[E'(A | SP) (A | SA)]2
d = (2-10)
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 
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 2-11 is called log-likelihood ratio.
(Note: log-likelihood ratio is unrelated to log-likelihood 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 log-likelihood ratio is rarely calculable,
except for the detection of an exactly known signal superimposed on an exactly known
background (SKE-BKE) 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 .
As the performance of the human observer was found to be dramatically degraded
by correlated noise [41, 42], and some spatial-frequency channels were found to exist in
the human visual system [43, 44], Myers et al.  added a frequency-selective 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 [46-50].
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) (2-12)
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. The interclass scatter
matrix S1 measures how far the state means deviate from the grand mean f and is
S, _C pf Pk( k k)T (2-13)
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)
The intraclass scatter matrix S2 is defined as the average covariance matrix across all
S2 a PkKk (2-15)
where the kth state covariance matrix Kk is given by
Kk=(f kX/-7k k,' (2-16)
The angular brackets in Equation 2-16 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,) (2-17)
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 . 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, (2-1 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 . Using a bank of2~band-pass 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 (2-19)
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 (2-20)
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 SKE-BKE detection tasks. Gifford et al.  investigated the application of CHO to
detect hepatic lesions using various channel models, eye filters and internal noise models.
Burgess et al.  showed the correlation between the human observer and CHO for
well-defined two-component noise fields (Poisson noise correlated with a Gaussian
function) with various channel models.
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 log-likelihood 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 forward-proj section matrix that gives
Hf ,c i)= [hes J, j ( 3-3 )
The second term in Equation 3-1 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 (3-4)
The third term is the strain energy term of deforming elastic material:
Es(m)= IAlrux y zI 2d+ lyu2 2 z2
+ 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 hyper-parameters 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). (3-6)
3.2 Euler-Lagrange Equations
We use the Frechet derivatives of the function E to obtain the Euler-Lagrange
equations. As the first step, we compute the Gateaux differentials  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 well-known formula as
DE(.fif2,m;t)= ~E~f, +srt, f2 0~s (3-7)
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(3-10)
Dz,:(I,72m~t= 2 [,(r- f 2(r + m(r))t(r)dr (3-11)
DE(II72,~t) -2[/()- 2 2(r + m(r))t(r)dr (3-12)
DzE,E(m; t) = [(jlV m + 2puxl)tr + pl(u~ + vx)tv + pr(uz, x )tz ]dr (3-14)
DE,E(m; t) = [(lV -m +2pny,)ty + p(ul + vr)tz*l 1Vz + y z]dr (3-15)
D_,E (m; t) = [(AV m + 2pez)tz + pr(uz + x)tr #r(' y y~t ]dr (3-16)
Now, let Cc (R) denote the set of all infinitely differentiable functions with
compact support on s Then the Euler-Lagrange 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 3-8, 3-11,
3-12, and 3-13 by easy calculation:
DE,(fi,, 2,m) = 2[f,(r)- f2(r + m(r))] (3-17)
D E, (f,,f2, m) = -2[ f, (r)- f2( 2 ~r) (r + m(r)) (3-18)
DE E(II,2, fm) = -2[J;(r)- f2( 2 ~r) (r + m(r)) (3-19)
DE E(II,2, fm) = -2[f, (r) f2( 2 ~r) (r + m(r)) (3-20)
To obtain D2Ef;(, f2, m), We need to express the integral in Equation 3-9 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
JWxr) [det(I 3 (3-21)
where m'(r) =~ vx vy (-2
So, if a is the image of R under the operator (I3+m), fTOm Equation 3-9 we obtain
D2E:(fi72, ; t = 2 [ (IS F) (F)t(F dr(3-23)
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 3-21, 3-23 we obtain
D2E, (ft,, 2,m) = -2[f, (r)- f2(r)]/|I det(3 + m' (r))|1 (3-24)
where r = (I; -1)(r
The Frechet derivatives of L are given by
where h(r) [ ha (r) So, from Equations 3-17, 3-24, 3-25 and the fact that E does not
depend onf,,, we obtain the Frechet derivative of E with respect to fi and f2
g~'Hf, (i)(ih r
D2E( ,f,, fm )
We now determine the Frechet derivatives of Es with respect to u, v, w. The
Gateaux differentials in Equations 3-14, 3-15, 3-16 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 3-14, 3-15, 3-16, 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
2[ f (r 2 ()]/detI3 m'(r) + ~ r)- ()H (i)(ih, (r)1
(AZ + 2pu)uxx + pu(u, + uz)+ (/ + Pu)(xy xzw) (3-32)
(AZ + pu)(uxy+ yz)+ xxv, z) + (AZ + 2pu)v,, (3-33)
(AZ+ pu)(uxz + vyz)+ Pu(w + yy)+ (AZ+ 2pu)waz (3-34)
Hence we obtain the Frechet derivatives of Es
Dz,Es (m) = -(AZ + 2pu)uxx pu(u, + uzz) (/ + Pu)(xy xzw) (3-35)
DV,Es(m) = -(AZ + pu)(uxy + yz) uvxx zz) (A + 2pu)v,, (3-36)
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 (3-3 8)
2[ f (r) f2( 2 ~)l ( ~)
DE(f,, f2,m -[i u(xy yz) uvxx zz) (il+ 2U)vw]
2[J; (r) f2( 2 ~)l ( ~)
DwE(f,, f2,Il -[i ~ u(xz vyz uxx yyW) (jl+ 2p)waz i
2[f, (r) f2( 2 ~)l ( ~)
Finally, the Euler-Lagrange equations for E, satisfied by (f, f2 ,m ), are
2[ zr)-f2(~m~r)]+ h~)- h(r~ *()/Hf(i)= 0(3-41)
-2[ z(r)- f2(r)]/ det(I3 + m'(r) + hj(r)- h(r)rj()i/H2i (3-42)
2[ (r) -f2 2ml~f
B(A+2pU)ux + p(ul~ +uzz) (~ )(V xy xz1?) = 0
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 3-43, 3-44, 3-45 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 two-step 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 3-1 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))
= 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)] (3-47)
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)
Figure 3-1. Scheme of the RM algorithm.
3.3.2 Computation of R Step
Equation 3-46 can be regarded as a penalized maximum-likelihood 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
Green's one-step-late (OSL) algorithm  was proposed for the penalized
maximum-likelihood estimation of a single frame using a general class of priors. It was
used  for reconstructing gated frames, and convergence was obtained for a modified
version . We used a modified OSL for our R-step optimization.
Since EI and L are convex infi, f2, the Kuhn-Tucker optimality conditions are
necessary and sufficient for the minimizer. That is
Dm E(fi, f2," (k) ) > 0 (3-48)
and fmDmE( fi,f2, m(k)) = 0 for m =1,2. (3-49)
From Equations 3-26, 3-27, we define
af ha(r) g (i)Hf; (i)
Az (f,, f27M m> =1 (3-50)
ah7(r) + 2[ (r) f2( ~
A2 17, 2, _> 1=1 (3-51)
ah(r) 2[ (r)- f2 (r)] | det(I3+m'(^)
Hence the necessary and sufficient conditions for Equation 3-46 are
Am1 2mk) :, m =1,2 (3-52)
fAm(11,2,m~) mr''>=, m =1,2 (3-53)
So, it seems natural to estimate f,(ki+1) by the fixed-point iterations given by
my I 1+ = fml,t~l $1(f; 2f,17 (k)), l= 1,2,... (3-54)
and f,( k+) p ~kl)
If the image matching term EI is not present in the obj ective function in Equation
3-46, then Equation 3-54 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 3-54 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-~-l-A' zr) ,0):rf e ,, j=; 1,2)< 1 (3-55)
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
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 motion-estimation problem
considered by Klein and Huesman . 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
Hence, approximately we need only to compute
m'k+1) arg min Qk,/,m) + PEs (m)] (3-58)
where Qk,l ( k)=5 i')(r _i 2(k (k _~ic, _mr (kx)cr, (kI)c mi(k)~
From the Frechet derivatives of Qk,l and Es, the Euler-Lagrange equations for
m = mI are
2 i(k) 2i(k) (k) _I--I (k) )7 2 (k)( (k 2I ~ nd(k)(k
+ P(t/+ 2pU)ux + iuO,.+> uzz (;3 xyU),' xz,) = 0
2 (k,() 2() ( k m ") _m (k) 2(k (k) (k):(k
+ P + 2pl)v,, + pl(v, zzv)+(i + pl)(uxy ~ y)] = 0
2 () 2() k) (k 2(k) () ] 2(k) (k))
+ 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 3-60, 3-61, and 3-62 with the initialization
m k) (k) Appendix gives details of the numerical analysis for motion estimation.
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
semi-ellipsoid 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 semi-ellipsoid 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 4-1.
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
high-resolution (120 X 120 X 120) obj ect with the steps illustrated in Figure 4-2.
Two simulated data files are generated for each frame: a noise-added 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 PS-MLEM algorithm, and the true
image will be used to evaluate the reconstructed images.
Figure 4-1. 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 end-systolic and mid-systolic phases of the heart cycle, respectively. The source
obj ects of the 2 frames had different lengths of semi-axes in a manner that achieved a
constant myocardium volume across the 2 frames and an anatomically realistic
myocardial-to-chamb 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 . The lengths in mm of the
ellipsoid semi-axes for the 2 frames are shown in Table 4-1. Two of the three ellipsoid
semi-axes were equal in length and represented the circular, short-axis slices.
Figure 4-2. Simulation steps of the geometric phantom
Table 4-1. Ellipsoid semi-axes 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 mid-ventricular points
along the long axis, as indicated by the arrows in Figure 4-1.
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 Tc-99m
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  and
tagged MR images of the human heart . 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 short-axis slices (#5,
#15 and #20) are shown in Figure 4-3.
4.1.5 Forward Projection Simulation and Detector Response
Collimator-detector response was simulated assuming a low-energy high-resolution
parallel-hole 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 4-3. Ideal motion fields of the geometric phantom in 3 short-axis slices.
A) Slice #5. B) Slice #15. C) Slice #20.
The high-resolution proj section data were then collapsed to a low-resolution
proj section (30 X 60 X 30) by a procedure described in Figure 4-4 (only a 2D slice is
shown). While the final projection dimensions (30 X 60 X 30) are small relative to those
used by large field-of-view 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 end-systole 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 Bi-Value 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
bi-value Lame constants model: a fairly incompressible model for myocardium and a
compressible model for the other regions.
High-resolution 120 X 60 Proj section
Figure 4-4. Collapse high-resolution proj section to low-resolution 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=~ (4-1)
(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
. The ii term in the constraint equation penalizes non-zero divergence and the pu
term penalizes sharp discontinuities in the flow field .
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 4-2 gives the Lame constants we used in our study for some
specific Poisson ratios.
Table 4-2. 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 bi-value 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 threshold-based procedure. Firstly, we reconstructed the
image from the noisy data using the PS-MLEM algorithm with Nyquist cut-off frequency.
Then we compared each voxel intensity of the resulting image to a pre-chosen 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
The segmentation procedure and results are illustrated in Figure 4-5. 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
Figure 4-5. Segmentation for bi-value 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 4-6. The total
obj ective function E and the negative log-likelihood 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 4-7 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
4-8. The obj ective functional value in M step achieved with the PS-MLEM images is
larger than that with true images. This reflects the inconsistencies between the estimated
frames by PS-MLEM 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
2.4 x 105
Using ideal emotion
-2.7- U~sing: estimated motion
-32 10 20 30 40 50r
Figure 4-7. Convergence of R step
S- Using ideal images
10000 ."- Using PS MlLEM images
0 10 20 30 40 50
Figure 4-8. Convergence of M step
where f, and f2 are the true images, and in is the estimated motion vector field.
Regional motion evaluation is performed at some selected sub-regions. 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 4-2) and
hyper-parameter 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 4-9. 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 , that is, for the fairly incompressible
deformation, the motion-estimation 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.
O 0.02 0.04 0.06 0.OB 0.1
Hyper-pararneter in M step: P
Figure 4-9. Global motion error vs. B for several Poisson ratios.
4.4.3 Global Motion Evaluation
To compare the motion-estimation performance of RM and the individual
motion-estimation algorithm, we implemented the following 3 methods:
Method 1: M step using true images with change ofa.
Method 2: Individual method. M step using PS-MLEM images with change of
cut-off frequency and B;
Method 3: RM algorithm with change of hyper-parameter a and B.
Figure 4-10 shows the ideal motion superimposed on the true Frame 1 images of
the geometric phantom, in both long- and short- axis views. Figure 4-11 shows the
motion estimated by Method 1, superimposed on the true images. Figure 4-12 shows the
images and motion individually estimated by Method 2. Figure 4-13 shows the images
and motion simultaneously estimated by RM.
The minimum global motion error each method achieved is given in Table 4-3. 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 20-image ensemble. Method 1 and 2 used 40 iterations for M step, and Method 3
used 40 outer-loop iterations, each containing 1 M-step iteration and 1 R-step iteration.
Figure 4-10. Ideal motion vector fields superimposed on true images.
Figure 4-11. Estimated motion by M step using true images. (B = 0.02).
Figure 4-12. Individually reconstructed images by PS-MLEM and estimated motion by M
step. (cut-off = 0.3 cycle/pixel; B = 0.01).
Figure 4-13. Estimated motion and reconstructed images by RM. (a=0.02; J
Table 4-3. Global motion evaluation results of the geometric phantom
Method Optimum Conditions Minimum global motion error
1 B = 0.02. 21.4
2 cut-off = 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 4-3 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 sub-regions 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) mid-ventricle and in normal region. Each sub-region was a
2 X 2 X 2 voxel volume. Their locations are indicated by the squares in a slice view in
Figure 4-14. To obtain the magnitude and angle errors of motion estimate, the mean
magnitude and angle of the estimated motion in the sub-regions were computed and
compared with the mean magnitude and angle of the ideal motion in those sub-regions
accordingly. These errors were computed for a 20-image 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 4-4. 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
motion-estimation accuracy for RM compared with the individual method.
Figure 4-14. Chosen region to calculate regional motion error for the geometric phantom.
Table 4-4. Regional motion evaluation results of the geometric phantom
Individual method RM
Sub-region 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
PS-MLEM images. The evaluation methods include global image evaluation, regional
image evaluation, and defect contrast evaluation. Since the cut-off frequency of the
post-filtering used in the PS-MLEM algorithm can greatly affect the quality of the
resulting images, a range of cut-off frequencies was investigated.
The figure-of-merit for the global image evaluation was the Root Mean Squared
(RMS) image error defined as
RMS = fr-fr43
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 2-frame system,
the mean RMS error was computed across the 2 frames.
The regional image evaluation also considered the 3 sub-regions used for the
regional motion evaluation. For each sub-region, the mean bias and variance of the
intensity of the estimated image were computed across a 20-image ensemble. For the
PS-MLEM algorithm, a set of bias/variance points was computed across a range of
cut-off frequencies from 0. 1 to 1.0 cycle/pixel. The bias/variance curves were plotted for
each sub-region to show the trade-off between the absolute value of the bias and the
variance, as the cut-off frequency changed.
Defect contrast is defined as
Contrast = normal defect (4 -4)
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 20-image ensemble, and then the mean contrast
was computed. Like the regional image evaluation, a range of cut-off frequencies was
considered for the PS-MLEM algorithm.
4.5.2 Global Image Evaluation
The results of the global image evaluation for RM and PS-MLEM are given in
Figure 4-15. 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 PS-MLEM algorithm shows a
minimum RMS error with a cut-off frequency of approximately 0.4 cycles/pixel. This
minimum, however, is still substantially larger than the stable RMS error of the RM
algorithm (0.045 compared with 0.03).
The hyper-parameter effects on the RM images are illustrated in Figure 4-16. 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 hyper-parameter effect is shown in Figure 4-17. 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 4-16 (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
Figre4-1. MS rrrsof heRM ndPS-LE iage. ) R. ) P-MEM
A B C
Figure 4-16. Hyper-parameter 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).
3 0 04-
Figure 4-17. 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 4-18. With the
PS-MLEM algorithm for both defects, the defect contrast increases as the cut-off
frequency increases and, at the highest cut-off frequency, is slightly less than that of the
RM image. Although at this cut-off frequency the 2 methods have similar defect contrast,
the quality of the images is substantially different, as shown in Figure 4-19.
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
PS-MLEM cut-off frequency (cyclelpixel) A PS-MLEM cut-off frecluency (cyclelpixel) B
Figure 4-18. Defect contrast of the PS-MLEM and RM images. A) Defect 1. B) Defect 2.
Figure 4-19 (A) and (B) show the images of Frame 1 for RM (40 iteration) and
PS-MLEM (1.4 cycles/pixel cut-off). It is evident that while the images have similar
contrast, the RM image has substantially lower noise level. Also included in Figure 4-19
(C) are the PS-MLEM images with Nyquist cut-off 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 sub-regions (indicated by
the squares in Figure 4-20) were calculated in both 2 frames. Figure 4-21 shows the plots
of the absolute bias versus variance for the PS-MLEM and RM images. In each plot, the
RM bias/variance point is downward and left to the PS-MLEM 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 PS-MLEM.
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 PS-MLEM reconstructed
images. Also we evaluated the image reconstruction by comparing the RM images to the
PS-MLEM 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 4-19. The trade-off between the contrast and noise level of the geometric phantom.
A) RM images. B) PS-MLEM images at similar defect contrast but higher
noise level. C) PS-MLEM images at similar noise level but lower contrast.
Figure 4-20. Sub-regions for regional image evaluation for the geometric phantom.
A) Frame 1. B) Frame 2.
O 0.05 0.1
0 0.05 0.1
- + RM
CO 0 045
O 0.05 0.1
Figure 4-21. The trade-off between regional intensity bias and variance of the geometric
phantom. A) Sub-region in basal normal region of Frame 1. B) Sub-region in
basal normal region of Frame 2. C) Sub-region in basal defect region of Frame
1. D) Sub-region in basal defect region of Frame 2. E) Sub-region in
mid-ventricle normal region of Frame 1. F) Sub-region in mid-ventricle
normal region of Frame 2.
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~URBS-Based 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 Non-Uniform Rational
B-Splines (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: end-diastole and
end-systole, 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 5-1 shows the long-axis and short-axis slice views of these 2 frames.
Figure 5-1. Long- and short- axis views of the NCAT cardiac phantom. A) Frame 1:
end-diastole. B) Frame 2: end-systole.
5.2 Motion Estimation Results
5.2.1 Global Motion Evaluation
As indicated in Klein' s study of the physical phantom , 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 PS-MLEM images with change of
cut-off frequency and B;
Method 3: RM algorithm with change of hyper-parameter a and B.
Figure 5-2 shows the ideal motion superimposed on the true frame 1 images of the
NCAT phantom, in both long- and short- axis views. Figure 5-3 shows estimated motion
by M step using true images, superimposed on the true images. Figure 5-4 shows the
images and motion individually estimated by Method 2. Figure 5-5 shows the images and
motion simultaneously estimated by RM.
Figure 5-2. Ideal motion vector fields superimposed on true images
Figure 5-3. Estimated motion by M step using true images. (B = 0.10).
Figure 5-4. Individually reconstructed images by PS-MLEM and estimated motion by M
step. (cut-off = 0.4 cycle/pixel; B = 0.05).
Figure 5-5. Estimated motion and reconstructed images by RM. (a=0.10; B = 0.05).
The global motion evaluation results are shown in Figure 5-6. The minimum global
motion errors of the 3 methods are shown in Table 5-1. The numbers in the parentheses
are the standard deviations of global motion error calculated from a 20-image ensemble.
Each method implemented 40 iterations.
Table 5-1. Global motion evaluation results of the NCAT phantom
Method Optimum conditions Minimum global motion error
1 B = 0.10. 267.7
2 cut-off = 0.4 cycle/pixel; 304.1 (11.3)
B = 0.05.
3 a=0.10; 288.6 (5.8)
B = 0.05.
O 02 O 04 0.06 O0 OBO1 0 12 0.14 O 16
r~~ : 300
a 03 200
0.12 1 1
P 0.04 0.4 OO0
0.02 0.2 cut-off
Figure 5-6. Global motion errors of three methods. A) global motion error vs. a for
Method 1. B) global motion error vs. cut-off 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 sub-region within the defect region in slice 7 (indicated by the square in
Figure 5-7) was chosen for the regional motion evaluation. The results are shown in
Table 5-2. 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 motion-estimation accuracy for the RM
algorithm compared with the individual method.
Figure 5-7. Chosen region to calculate regional motion error for the NCAT phantom
Table 5-2. 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 image-reconstruction 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
PS-MLEM images over a range of cut-off frequencies.
5.3.1 Global Image Evaluation
For PS-MLEM algorithm, the mean RMS error of 2 frames changed with the
cut-off frequency, as indicated in Figure 5-8 (A). The minimum RMS error (0.284) was
achieved at approximately 0.6 cycles/pixel cut-off. For RM algorithm, a surface
representing the change of the mean RMS error with the hyper-parameters (a, a) is shown
in Figure 5-8 (B). The minimum RMS error of RM (0.268) was achieved when a=0.13
Figure 5-9 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,
then stayed flat up to iteration 60. This demonstrates again the stability of the RM image
with the outer-loop iteration.
2 0.3 .1
0. e 01 .0.
O 0. .015 .
~~f0. 0. .2
0 ~ ~ 20 400 6050
Itertio nurnberc of RMM oute loo
Figure 5-9. RMS errorvs.o ithe rato anume ofRM outger. loop. rrr s ct-
532DfrqetCntras EvaPS-luation 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 5-10 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 5-10. Normal and defect regions in the NCAT phantom. A) Frame 1. B) Frame2.
As indicated in Figure 5-11i, PS-MLEM using a higher cut-off frequency produced
a higher contrast. The highest contrast of PS-MLEM images (0.53) was achieved with the
cut-off of 1.4 cycle/pixel. RM images had a lower contrast (0.5) compared with the
highest contrast achieved by PS-MLEM. However, RM images presented better trade-off
between the contrast and noise than PS-MLEM. That is, RM had lower noise level while
their contrasts were similar, and higher contrast while their noise levels were similar.
Table 5-3 illustrates this trade-off between the contrast and noise of RM and PS-MLEM.
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 PS-MLEM
image had comparable noise level at the cut-off of 0.4 cycle/pixel, but the contrast was
substantially reduced to 0.43. At a higher cut-off 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 5-12.
Figure 5-11i. Contrast vs. cut-off frequency of PS-MLEM
Figure 5-12. The trade-off between the contrast and noise level of the NCAT phantom.
A) RM images. B) PS-MLEM images with similar contrasts but higher noise
levels. C) PS-MLEM images with similar noise levels but lower contrasts.
Table 5-3. The contrast and noise level of the NCAT phantom
Method Contrast Std. Dev.
RM 0.50 0.15
PS-MLEM 0.50 0.25
PS-MLEM 0.43 0.15
5.3.3 Regional Image Evaluation
Figure 5-13 shows the 2 sub-regions of interest in Frame 1 and Frame 2: one in the
normal region and the other in the defect region. Figure 5-14 shows the trade-off between
intensity bias and variance at these sub-regions in Frame 1 and Frame 2. We investigated
5 sets of hyper-parameters (see Table 5-4) 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 hyper-parameters presented different bias/variance property, as
indicated in Figure 5-14. For the sub-region in the defect of Frame 1, all of the Hyve RM
images presented close or slightly worse trade-off between intensity bias and variance
than PS-MLEM. For the sub-region in the defect of Frame 2 or the normal of Frame 1
and 2, all of the Hyve RM images presented better trade-off between intensity bias and
variance than PS-MLEM. The one out of these 5 sets that achieved the best bias/variance
trade-off 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 5-4. Five sets of hyper-parameters 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
\+ RM: set 1
+ RM: set 2
RM: set 3
\x RM: set 4
RM: set 5
Figure 5-13. Sub-regions for regional image evaluation for the NCAT phantom.
A) Frame 1. B) Frame 2.
\+ RM: set 1
\+ RM: set 2
\RM: set 3
x RM: set 4
i~RM: set 5
0.05 0.1 0.15
0.05 0.1 0.16
O 0.1 0.2 0.3 0.4
0.05 0.1 0.15 0.2
Figure 5-14. The trade-off between regional intensity bias and variance of the NCAT
phantom. A) Sub-region in the defect region of Frame 1. B) Sub-region in the
defect region of Frame 2. C) Sub-region in the normal region of Frame 1. D)
Sub-region in the normal region of Frame 2.
OBSERVER EVALUATION OF IMAGE QUALITY
6.1.1 Single-Slice CHO Model
Our study of observer evaluation of image quality for the RM algorithm started
with a single-slice CHO (SCHO) model. A number of channel models have been
presented since Myers et al.  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 band-pass cut-off 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, non-overlapping channels
with cut-off frequencies of (0.03 125, 0.0625), (0.0625, 0. 125), (0. 125, 0.25), (0.25, 0.5)
cycles/pixel. Choice of cut-off frequencies was somewhat arbitrary. Figure 6-1 shows the
RSC-model channels in the frequency domain.
The channelization procedure was described in Equation 2-19 as v = Uf To
compute the vector v from the image ffor a location-known defect, we Fourier
transformed each frequency-domain channel to a spatial-domain template centered at the
defect location. The ith COmponent of the vector v was then calculated by integration of
the pixel-by-pixel product of the image and the ith spatial-domain template. This is to
avoid excessive Fourier transformation for the large number of images.
Figure 6-1. The RSC-model channels used for channelization in single-slice CHO model.
(Note: Images are shown in the frequency domain.)
Figure 6-2 shows the 2 defects in the geometric phantom and the defect in the NCAT
phantom. Figure 6-3 shows the spatial-domain templates of the 4 RSC-model channels
used in our study. Three rows indicate the different defect locations in the two phantoms.
Figure 6-2. Defect locations in the phantoms. A) Geometric phantom. B) NCAT
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 N-channel 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 4-channel 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 PS-MLEM, 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 (SCHO-d,) 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 SCHO-d,. 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 SCHO-d,.
6.1.2 Multi-Frame Multi-Slice CHO-HO 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.  presented a multi-slice multi-view CHO-HO model for ungated SPECT
myocardial perfusion imaging. We modified it to a multi-frame multi-slice CHO-HO
(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 multi-frame multi-slice 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 (MMCHO-d,).
Figure 6-4 shows the scheme for this 2-step procedure.
Figure 6-3. The RSC-model channels represented as spatial-domain 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
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 test-statistic vectors as the training set to determine the HO discriminant
vector, and then applied the discriminant to the remaining test-statistic vectors to
compute the MMCHO-d,. The mean and standard deviation of MMCHO-d, were
computed by the same procedure used in the SCHO model.
3 slices 3 slices
S~A- CHO HO O SP~- CHO -HO
Multi-frame multi-slice CHO-HO detectability index
Figure 6-4. Two-step procedure of the multi-frame multi-slice CHO-HO model.
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 PS-MLEM and RM are shown in Figure 6-5 and Figure 6-6, respectively. In each
figure, the image slices for each frame are grouped in the order of #14, #15, and #16. The
cut-off frequency for PS-MLEM was 0.3 cycle/pixel. The RM algorithm used the
hyper-parameters (a=0.02, 8=0.04), which produced minimum RMS image error. The
PS-MLEM images present more smoothness than the RM images, particularly at the
6.2.1 Single-Slice CHO Results
We performed SCHO model to the 2 defects in the 3 long-axis slices of the 2
frames. For PS-MLEM, we used a range of cut-off frequency from 0. 1 to 1.0 cycle/pixel
and found the minimum SCHO-d,. The results of PS-MLEM and RM are shown in Table
6-1. Standard deviations of the SCHO-d, are given in parentheses.
The RM images had significantly improved SCHO-d, 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 multi-frame multi-slice observer model that produces a single
detectability index for a series of gated myocardial ET images.
6.2.2 Multi-Frame Multi-Slice CHO-HO Results
We applied the MMCHO model to the 2 frames and 3 slices shown in Figure 6-5
and Figure 6-6. The results are shown in Table 6-2. The numbers in parentheses are the
Figure 6-5. SP images of the geometric phantom reconstructed by PS-MLEM.
A) Frame 1. B) Frame 2.
Figure 6-6. SP images of the geometric phantom reconstructed by RM.
A) Frame 1. B) Frame 2.
standard deviations of MMCHO-d,. The results show that the RM images presented
significantly higher MMCHO detectability than the PS-MLEM images. Figure 6-7 shows
how the MMCHO-d, of the PS-MLEM images changed with the cut-off frequency. For
both defects, the maximum detectability indices of PS-MLEM were achieved with the
cut-off of 0.4 cycle/pixel. The MMCHO-,' s of the RM images are also plotted.
Table 6-1. SCHO detectability index results of the geometric phantom
SCHO-d, 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)
SCHO-d, 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)
0 0. 2 0. 4 0. 6 0. 8 1 1. 2
cut-off frequency (cycle,'pixel)
0 0.2 0.4 0.6 0.8 1 1.2
cut-off frequency (cycle,'pixel)
Figure 6-7. 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
Table 6-2. MMCHO detectability index results of the geometric phantom
MMCHO-d, for Defect 1 MMCHO-d, for Defect 2
PS-MLEM 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 6-8 shows the true image slices of two frames. Figure
6-9 shows the images reconstructed by the PS-MLEM algorithm using the cut-off
frequency of 0.3 cycle/pixel. Figure 6-10 shows the images reconstructed by the RM
algorithm using the hyper-parameters 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 Single-Slice CHO Results
The results of single-slice CHO performance for the PS-MLEM and RM images
are shown in Table 6-3. Again, the minimum SCHO-da for PS-MLEM was found by
using a range of cut-off frequencies from 0. 1 to 1.0 cycle/pixel.
Figure 6-8. SP true images of the NCAT phantom.
A) Frame 1. B) Frame 2.
Figure 6-9. SP images of the NCAT phantom reconstructed by PS-MLEM.
A) Frame 1. B) Frame 2
Figure 6-10. SP images of the NCAT phantom reconstructed by RM.
A) Frame 1. B) Frame 2
Table 6-3. SCHO detectability index results of the NCAT phantom
Slice Algorithm SCHO-d, for Frame 1 SCHO-d, for Frame 2
7 PS-MLEM 11.85(0.23) 12.81(0.64)
RM 20.35(0.80) 19.80(1.11)
8 PS-MLEM 10.71(0.30) 13.01(0.81)
RM 18.94(0.69) 18.45(0.45)
9 PS-MLEM 8.01(0.26) 9.05(0.46)
RM 8.32(0.35) 16.66(0.49)
6.3.2 Multi-Frame Multi-Slice CHO-HO 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 6-4. It is evident that for both 50% and 25%
defect phantoms, the RM images presented significantly higher MMCHO detectability
than the PS-MLEM images. Figure 6-11 shows how the MMCHO-d, of the PS-MLEM
images changed with the cut-off frequency.
Table 6-4. MMCHO detectability index results of the NCAT phantom
MMCHO-d, for 50% defect MMCHO-d, for 25% defect
PS-MLEM 16.98(0.77) 8.14 (0.36)
RM 33.41(0.66) 14.59 (0.52)
6.4 Multi-Frame Multi-Slice Multi-View CHO-HO 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
multi-frame multi-slice CHO-HO model to a multi-frame multi-slice multi-view
CHO-HO (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 test-statistic 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 test-statistic
vectors to compute a Einal MMMCHO detectability index (MMMCHO- da).
35 ~- -PS-MLEM __-A-PS-MLEM
0 0.2 0.4 0.6 0.8 1 20 0. 2 0. 4 0. 6 0. 8 1 1. 2
cut-off frequency (cycle/pixel) Acut-off frequency (cycle/pixel) B
Figure 6-1 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 6-12, 6-13 and
6-14, presenting the PS-MLEM algorithm (cut-off = 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 test-statistic 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 test-statistic 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
test-statistic vectors (50 SP and 50 SA) to estimate the MMMCHO detectability index
The MMMCHO-d, versus cut-off frequency of PS-MLEM is shown in Figure 6-15.
The highest MMMCHO-d, of PS-MLEM was achieved with the cut-off of 0.4 cycle/pixel.
The MMMCHO-,' s and their standard deviations of the RM images and the PS-MLEM
images are shown in Table 6-5.
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 MMMCHO-d, as a function of a and B.
Figure 6-12. Multi-frame multi-slice multi-view SP images reconstructed by PS-MLEM.
Figure 6-13. Multi-frame multi-slice multi-view SP images reconstructed by RM
Figure 6-14. Multi-frame multi-slice multi-view SP images reconstructed by RM
.~16 -r -4-RM 2
0 0. 2 0. 4 0. 6 0. 8 1 1. 2
cut-of ffrequency (cycle/pixel)
Figure 6-15. MMMCHO detectability index results of the NCAT phantom.
Table 6-5. MMMCHO detectability index results of the NCAT phantom.
PS-MLEM 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 single-slice CHO model and two multi-image
CHO-HO 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 PS-MLEM algorithm.
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
non-rigid 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 motion-estimation 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
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 PS-MLEM.
Our results showed that the RM images presented lower global image error and regional
image error, relative to PS-MLEM. 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 contrast-noise trade-off over the PS-MLEM images. The RM images
had a substantially higher defect contrast compared to the PS-MLEM images except for
the highest cut-off frequency PS-MLEM images. However, at this high cut-off frequency,
the PS-MLEM 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 PS-MLEM 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 single-slice CHO model, we proposed a multi-image CHO-HO model for
gated myocardial ET. Results with both the single-slice CHO model and the multi-image
CHO-HO model demonstrated that the RM images had significantly improved
detectability than the PS-MLEM images.
Our future work includes
1. Multi-frame 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, multi-frame processing will expand the current two-frame 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 motion-estimation accuracy for gated myocardial ET, relative to
the independent motion-estimation and image-reconstruction algorithms. More studies
are required to investigate the performance of the RM algorithm in clinical application.
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 first-order derivative is a difference quotient:
du .u(x + x) u(x)
= hm (A-1)
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)
du u(x) u(x &x)
Backward difference: 4 (A-3)
Central difference: 22
and the central difference approximation for the second-order derivative:
d2u u(x + &) 2(x)+ u(x h)
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 motion-estimation problem, we rewrite the
Equations 3-60, 3-61, 3-62 as
p[(A + 2~)u_ + p~(uY + uzz) + (I C P)(xy xzW ) 2m Vf2(k)( k) 2(,k) k
= 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 ''
= -2 fi*(k) 2(k) + ik)) Ilk)) V2(k) ( lik)) k)+II
$(A+ 2~w,+ pw +w )+ ( + )(u yz m -Vf2k) (k)) 2(k) +(k))
= f(k 2() (k) k) 2(k (k)1 2(k) (k)I
where k denotes the kth outer loop, I denotes the lth M-step inner loop, and m k) is
current motion estimation. The solution m of Equations A-6, A-7, A-8 is the motion
estimation for the (Il+)th M-step inner loop. Equations A-6, A-7, A-8 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
u u+1,],k -1Ik, +1,k 1, ,],k+1 ,,- A 9
x2 2 a2
uxx = (uz+,l,,k ],k ,, 2-1], l+,,,k 2ur], 2-1] (A-10)
uW = (u",J+1,k U,,,k (U,,,k 2,J-1,k Ui,J+1,k 2uI,,,k + 2,J-1,k (A-11)
uzz = (ul,],k+1 U, ,k (U,,,k 2,,,k-1 Ul,],k+1 2uI,,,k + 2,,,k-1 (A-12)
u (u+1,+1, 2-,J+,k +,J-,k 1,J1,k(A-13)
u ( rJ1k1 2J1k+ J1k1 ,-,- (A-14)
For each voxel (i, j, k), Equations A-6, A-7, A-8 generate 3 rows of entries in the
matrix A. The coefficients of (uz, u27"'>, N) for voxel (i, j, k) from Equation A-6 can be
filled into a chart as follows, where i-1, i, i+1, j-1,j +, k-1, 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.
i-1 0 0 0 0 (a +2pu)P 0 0 0
i pf 0 upf df2 (k) 2 pf 0 up
i+1 0 O2)
l+pu i+pu 0
i-1 0 0 0 P 0 p 0 0
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
.- 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
j-1 j j+1
j-1 J j+1
Similarly, the coefficients of (v,, v27"'>, N) in Equation A-6 for (i, j, k) can be expressed
using the following chart:
j-1 j j+1
j-1 j j+1
The coefficients of (w,,w27"' N,) in Equation A-6 for voxel
(i, j, k)
j-1 j j+1
j-1 j j+1
The coefficients of (uz, u27"'>, N) for Equation A-7 for voxel
(i, j, k) are
(3/ (k) (f(k) 0 0
i 0 0 0 0 -2( 2 O
i+ o -p o~ 0 0
k-1 k k+1
The coefficients of (v, ,v27"'>v N) for Equation A-7 for voxel (i, j, k) are
j-1 j j+1 j-1 j j+1 j-1 j j+1
i-1 0 0 0 0 ##p 0 0 0 0
pup (a + 2 p) f(k (a + 2 p) f
i 0 0 of 2o 0 #
k-1 k k+1
The coefficients of (w ,w27"' N,) for Equation A-7 for voxel (i, j, k) are
j-1 j j+1 j-1 j j+1 j-1 j j+1
i-1 0 0 0 0 0 0 0 0
il+u p A+ p af2 (k) d2(k) lu
c az2( )(/ 4
j-1 j j+1
j-1 j j+1
The coefficients of (u, ,u27"'>u N) for Equation A-8 for voxel (i, j, k) are
i- Op o o o -p0
8f~:k) ~S(k) 0
i 0 0 0 0 ~- 2(2 2) 0 0 0
i+1 l+pu i+pu 0
k-1 k k+1
The coefficients of (v,, v27"'>, N) for Equation A-8 for voxel (i, j, k) are
j-1 j +1 j-1 j +1 j-1 j +1
i-1 0 0 0 0 0 0 0 0 0
i +u p f, (k) df(k) lu
4 y 8 4 0 4
k-1 k k+1
The coefficients of (w,,w27"' N,) for Equation A-8 for voxel (i, j, k) are
j-1 j +1 j-1 j +1 j-1 j +1
i-1 0 0 0 0 ##p 0 0 0 0
(2il + 8pu)P
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
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
(rrm)= (r )m)+& (r +m1)+ Sv (r +m1)+ &v (r +m1)
(r m =(r + m1) + At (r + m1) + Sv (r + m1) + Av (r + m1)
The differential of f (r + m1) can be approximated using central difference.
d, f,(i + u+1, j +v, k +w)- f,(i + u-1, j +v, k +w)
(r +m) =(A-20)