University of Southern California Dissertations and Theses

Dynamic cardiovascular magnetic resonance imaging for improved assessment of ischemic heart disease

(USC Thesis Other)

Dynamic cardiovascular magnetic resonance imaging for improved assessment of ischemic heart disease

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content DYNAMIC CARDIOVASCULAR MAGNETIC RESONANCE IMAGING FOR
IMPROVED ASSESSMENT OF ISCHEMIC HEART DISEASE
by
Taehoon Shin
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(ELECTRICAL ENGINEERING)
December 2009
Copyright 2009 Taehoon Shin
Acknowledgements
Many years have passed since I joined USC, and they now seem to have gone by in a
night’s dream. It is my great pleasure to recall my memories from those years, and to
express gratitude to the people to whom I am indebted.
I would like to thank my advisor, Professor Krishna Nayak, for his tremendous support
during my PhD studies. He introduced me to the field of medical imaging: a field that
I had not previously considered, but I now hope to continue to study in my professional
career. His keen intuition and fresh research experience have guided my research and
facilitated its progress. I have also been touched by his incredible optimism in the face of
any difficulties in research and in life.
I thank Professor Jerry Mendel and Dr. Sanza Kazadi for their guidance when I
was mentoring high school students at Jisan Research Institute. Professor Mendel has
shown attentive leadership and patience in working with the JRI students. Dr. Kazadi
demonstrated the importance of educating and inspiring potential young scientists.
I have been fortunate to have had Dr. Gerald Pohost as my clinical advisor. His
profound knowledge of physiology and medicine, and rich experience on clinical research
have been essential to my study of cardiac MRI. I also thank Dr. Ramdas Pai and
ii
Dr. Padmini Varadarajan for their generous support for clinical studies at Loma Linda
Hospital.
I have met many wonderful colleagues at the Magnetic Resonance Engineering Lab-
oratory. Previous post-docs, Dr. Chia-Ying Liu and Dr. Jon-Fredrik Nielsen have been
always generous with their advice. I have also been fortunate to have Dr. Houchun Hu
who has shown a genuine enthusiasm for research, which inspired me to greater efforts. I
thank earlier students, Dr. Kyunghyun Sung, Dr. Joao Carvalho, and Dr. Hsu-Lei Lee for
helping me get started, and for sharing their knowledge of MRI implementation. I thank
Zungho Zun for his time for productive and enjoyable discussion. I wish all younger stu-
dents the best with their research: Yoon-chul Kim, Mahender Makhijani, Samir Sharma,
and Travis Smith.
It is impossible to say enough about the unconditional love and support that I have
received from my family. My parents dedicated themselves to ensuring that I obtained
the finest education, and always put their children ahead of themselves. I also would like
to express my deep gratitude to my sister, Eunyoung for her constant help and advice.
iii
Table of Contents
Acknowledgements ii
List of Tables vii
List of Figures viii
Abstract xv
Chapter 1 : Introduction 1
1.1 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
Chapter 2 : MRI Background 6
2.1 Basic Principles of MRI . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.1 MR Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1.2 Selective Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.3 Spatial Encoding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Advanced Reconstruction of MR Images . . . . . . . . . . . . . . . . . . . 14
2.2.1 Non-Cartesian Reconstruction . . . . . . . . . . . . . . . . . . . . . 14
2.2.2 Parallel Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.3 Temporal Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.4 Compressed Sensing . . . . . . . . . . . . . . . . . . . . . . . . . . 21
Chapter 3 : Accelerated Spiral Imaging of Ventricular Function 24
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2 Reconstruction Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.1 Aliasing in Dynamic Spiral Imaging . . . . . . . . . . . . . . . . . 26
3.2.2 Algebraic Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.3 Rank Analysis and Zero Assumption . . . . . . . . . . . . . . . . . 32
3.2.4 Pseudoinverse and Regularization . . . . . . . . . . . . . . . . . . . 36
3.2.5 Iterative Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Numerical Phantom Simulation . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4 In-vivo Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4.1 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
iv
3.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Chapter 4 : Three-Dimensional Myocardial Perfusion Imaging: Feasibility 52
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2 Pulse Sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.3 Phantom experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.3.1 Phantom Description . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.3.2 Scan Protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.3.3 Image Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.4 In-vivo Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.4.1 Scan Protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.4.2 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.4.3 Image Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
Chapter 5 : Three-Dimensional Myocardial Perfusion Imaging: Systole vs Diastole 74
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.2.1 Pulse Sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.2.2 Imaging Protocols . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.2.3 Image Reconstruction and Analysis . . . . . . . . . . . . . . . . . . 78
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Chapter 6 : Three-Dimensional Myocardial Perfusion Imaging with Regularization 90
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.3.1 Pulse Sequence and Sampling Scheme . . . . . . . . . . . . . . . . 93
6.3.2 Scan Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.3.3 Image Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Chapter 7 : Summary and Future Work 105
7.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
7.1.1 Accelerated Spiral Cardiac Imaging . . . . . . . . . . . . . . . . . . 106
7.1.2 Three-Dimensional Myocardial Perfusion Imaging . . . . . . . . . . 107
7.1.3 Highly Accelerated 3D Myocardial Perfusion Imaging . . . . . . . . 108
v
References 109
Appendix A: Modulation of Spiral Trajectories 118
vi
List of Tables
4.1 Imaging parameters of 3D and 2D multi-slice scans. . . . . . . . . . . . . . 60
5.1 Imaging parameters of dual 3D MPI. . . . . . . . . . . . . . . . . . . . . . 77
5.2 Systolic and diastolic myocardial thicknesses, measured in pixel units. . . 79
6.1 Imaging parameters used for in-vivo experiments. . . . . . . . . . . . . . . 97
vii
List of Figures
1.1 Ischemiccascadeintheheart. Thenarrowingorblockageofcoronaryartery
causes reduction in blood flow in the myocardium, followed by metabolic
change, impairment of diastolic and systolic function, and ECG change.
This dissertation focuses on MRI methods for detecting impaired cardiac
function (ventricular function imaging), and blood flow reduction (myocar-
dial perfusion imaging) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Nuclear spins with and without a static field. (a): In the absence of static
field B
0
, spins align in a random way such that the net magnetization is
zero. (b): In the presence of B
0
, spins are lined up either in the direction
of B
0
or in the reverse direction. . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Precession of a nuclear spin around static field. In the presence of B
0
spins not only spin around their own magnetization axes, but this axis also
precesses around the direction of static field B
0
. . . . . . . . . . . . . . . 8
2.3 Excitation of a nuclear spin by a RF field. RF field B
1
(t) tuned to the
resonance frequency, tips down spins from the equilibrium state. This be-
havior is seen as (a): spiral motion in the absolute coordinate, and (b): as
simple rotation around the axis of B
1
(t) in the relative coordinate rotating
at the resonance frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4 Selective excitation using a RF field and gradient fields. Amplitude modifi-
cation of B
1
(t) field creates limited bandwidth such that only spins within
the band can be excited. With linear gradient turned on, resonance fre-
quencyofspinslinearlychangesacrossspatiallocation. Asaresult, spinsin
afinitespatialregionareselectivelyexcited. Theslicethicknessisincreased
by either decreasing G
z
or increasing the bandwidth of B
1
(t) . . . . . . . 12
2.5 Two examples of k-space trajectories. k-space is traversed by the time
integral of gradients. (a): 2DFT trajectory is the most commonly used,
where k-space is scanned in a raster fashion. (b): spiral trajectory, one of
the non-Cartesian trajectories, can cover k-space in more time-efficient way. 15
viii
2.6 Gridding of one-dimensional non-Cartesian data. Data from each non-
Cartesian sample denoted by “x” are convolved with a kernel, and re-
sampled onto neighboring Cartesian grids denoted by circles. In the image
domain, this procedure is equivalent to convolution with an impulse re-
sponse and apodization by the transform of the gridding kernel. . . . . . 16
2.7 Pictorial description of multi-coil acquisition in Cartesian two-fold under-
sampling. The true image is weighted by the sensitivity profile of individual
coil element, and is convolved by a point spread function. . . . . . . . . . 19
2.8 UNFOLDmethodintwo-foldCartesianundersampling. UNFOLDacquires
odd-order and even-orderk-space encoding lines alternately over time. (a):
The phase of the aliasing will alternate, and (b): the aliasing will be shifted
to the Nyquist frequency by a Fourier transform along temporal dimension.
The true signal can be separated from the aliasing by filtering. . . . . . . 20
2.9 Importance of incoherent aliasing in compressed sensing reconstruction.
(a): A Shepp Logan phantom (top) has only a few non-zeros total varia-
tion (TV) coefficients (bottom). (b): If k-space is uniformly undersampled
(top), a true TV coefficient is not distinguishable from its alias since its
contribution to the aliased signal is the same (bottom). (c): If k-space
is randomly undersampled (top), large TV coefficients are distinguishable
from their noise-like aliases (bottom). Smaller TV coefficients can be itera-
tively detected by subtracting the aliased signal by the contribution of the
detected coefficients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1 Proof of opposite sidelobe phases in the PSFs from interleaved spiral un-
dersampling. The top row shows spiral sampling function and the bottom
row shows the corresponding point spread functions (PSF). Sampling odd
ordered or even ordered spiral interleaves results in ring-shaped sidelobe as
well as mainlobe while the PSF from full sampling contains mainlobe only.
Since full sampling is the summation of the two undersampling functions
and Fourier transform is linear, the two aliased PSFs should add up to the
PSF from full sampling. Therefore, the sidelobe phases of the two aliased
PSF should be opposite to each other. . . . . . . . . . . . . . . . . . . . . 28
3.2 Illustration of dynamic undersampled spiral imaging in xyf space. From
two-fold undersampling, ring-shaped sidelobes are shifted to the Nyquist
frequency. Since the mainlobe and sidelobes are confined to DC and the
Nyquist frequency respectively, the solutions of only two temporal frequen-
cies are coupled to generate aliasing at one temporal frequency. . . . . . . 29
ix
3.3 Illustration of the system matrix in Cartesian two-fold undersampling with
the solution for two coupled temporal frequencies. The rank is half the ma-
trix dimension since the first and fourth quarter-columns (gray) are iden-
tical, and the second and third quarter-columns (white) are identical. To
condition rank deficiency, zeros can be assumed in certain regions of the
solution. P
k
and Q
k
(k =1;2) represent the region of zero assumption.
Accordingly, P
k
c
and Q
k
c
(k =1;2) represent support regions of the solu-
tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4 Diastolicandsystolicphasesofthedynamicnumericalphantom. Thephan-
tom was designed to resemble the heart in a short axis view, with the size
of the ventricle changing sinusoidally. . . . . . . . . . . . . . . . . . . . . . 40
3.5 Cross-sectional view of the zero assumption in the numerical phantom
study. The support region is cylindrical in xyf space. . . . . . . . . . . . . 41
3.6 Normalized RMS error and maximum pixel-by-pixel error from simulated
phantom experiments. The proposed method has smaller error than sliding
window reconstruction. Compared to filtering, slight improvements are ob-
served because of regularization effects. The large periodic increase in error
in sliding window reconstruction occurs when motion is the fastest (twice
per cycle). The periodic increase in errors for filtering and the proposed
method is related to the size of the dynamic region approaching the edge
of the support region at diastolic phases (once per cycle). . . . . . . . . . 43
3.7 Representative reconstructed images and residue images from three consec-
utive frames around early systolic phase. Residue images are the difference
between the reconstructed images and the known true images, and are dis-
played using harsh windowing. . . . . . . . . . . . . . . . . . . . . . . . . 44
3.8 Cross-sectionalviewofthezeroassumptioninthein-vivoexperiments. The
support region resembles a truncated cone in xyf space. . . . . . . . . . . 45
3.9 Evaluation of the depiction of fine cardiac structure in in-vivo images from
sliding window reconstruction, UNFOLD-like filtering and the proposed
method. (a): Signalprofilesacrossthemoderatorbandoftherightventricle
(black cross-hairs) were measured (20
th
time frame). (b): Intensity profiles
indicate sharp depiction of the moderator band in the order of the proposed
method, filtering and the sliding window. . . . . . . . . . . . . . . . . . . 47
x
3.10 Evaluation of motion artifacts in the three sets of in-vivo images. (a):
To visualize motion artifact reduction, short axis images from the three
methods were windowed harshly (59
th
time frame), and signal in a region
outside of the body (white box) was measured. (b): Mean signal within the
box is plotted as a function of time illustrating that the proposed method
reduces motion artifacts as well as filtering. . . . . . . . . . . . . . . . . . 48
4.1 Pulse sequence diagram and sampling pattern for 3D myocardial perfu-
sion imaging. The pulse sequence consists of a 90
±
BIR-4 adiabatic RF
pulse for myocardial saturation, followed by a saturation recovery time and
undersampled 3DFT GRE readout. Three-dimensional k-space was under-
sampled by factor of three and two in k
y
and in k
z
, respectively, resulting
in a net undersampling factor of six. The k-space sampling locations were
shifted in a cyclic pattern such that coil sensitivity map could be obtained
from the data in six neighboring R-R intervals. SENSE reconstruction of
undersampled data from the k
th
cardiac cycle utilized coil sensitivity maps
derived from the (k-2)
th
to (k+3)
th
cardiac cycles. . . . . . . . . . . . . . 55
4.2 Waveforms of saturation pulse and 3DFT imaging. BIR-4 RF pulse (8 ms)
was used for B1 insensitive 90
±
saturation. In 3DFT imaging, Gz rewinder
was combined with spoiler to reduce TR (2.3 ms). . . . . . . . . . . . . . . 56
4.3 Transverse component of magnetization over k-space encoding lines with
T
SR
=130 ms and ®=10
±
. The k-space response shows little variation over
encoding lines in the wide range of T
1
values. . . . . . . . . . . . . . . . . 57
4.4 Cardiac phantom and rubber insets. Left ventricular and myocardial cavi-
ties were filled with Gadolinium-doped saline of 2 mmol/L and 0.6 mmol/L.
The rubber insets ranging in size from 5-70 % of the myocardial volume,
mimic transmural perfusion defects. . . . . . . . . . . . . . . . . . . . . . . 58
4.5 Representative heart phantom images with a 20.3 % rubber inset. (a):
images from 3D scan, (b): images from 2D 4-slice scan, and (c): images
from 2D 3-slice scan. In 2D multi-slice case, images from the three sets of
cross-sectional scan planes were shown. In all phantom images, the dark
rimaroundtheLVwascreatedbythe»2mmthickseparationwallbetween
the LV and the myocardial cavities. . . . . . . . . . . . . . . . . . . . . . . 61
4.6 Distribution of errors in the estimation of phantom defect size. The hori-
zontal axis and vertical axis represent the true volume fraction of the de-
fect, and the estimation errors, respectively. In each plot, the solid and
dotted lines represent the mean and §1.96£standard deviation, respec-
tively. 3D imaging resulted in a significanty lower bias bias/ standard de-
viation (-0.44/1.49%) compared to 2D 4-slice (2.23/2.97%) and 2D 3-slice
(2.59/3.18%). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
xi
4.7 Representative 3DFT perfusion images selected every other heart beat (1
th
-
13
th
time frames). Blood pools and myocardium are well suppressed by
BIR-4 saturation pulse at the beginning. LV and RV get bright sequentially
by the arrival of contrast agent. . . . . . . . . . . . . . . . . . . . . . . . . 66
4.8 Representative 3DFT perfusion images selected every other heart beat
(15
th
-27
th
time frames). After LV enhancement, myocardium reaches peak
enhancement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.9 SNR at peak myocardial enhancement and CNR values for the first eight
slices. The SNR and CNR were 57.6§22.0 and 37.5§19.7 across all seg-
ments. Note that the two values are higher in slices at center than slices at
edges due to the slab excitation profile of parabolic shape. . . . . . . . . . 68
4.10 Time intensity curves in a healthy volunteer. Representative time intensity
curves (TIC) are shown from the first nine slices on a segment basis. Six
segments were used for the first six slices (basal to mid-short axis levels),
and four segments were used for the rest (apical level). All TICs show
homogenous myocardial enhancement, consistent with normal perfusion.
Flickering artifacts caused by errors in coil sensitivity map estimation were
observed (see an arrow in the 8
th
and 9
th
slices). The artifacts are most
severe at the start of LV enhancement, but have no significant effect on the
TIC upslope. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.11 Colormap visualization of time-intensity curve upslope in a healthy volun-
teer. Upslope value of time intensity curve was computed by linear fitting
of the curve during signal enhancement, and visualized by a color map.
The average upslope value of the whole myocardium is scaled to 100 %. . 71
5.1 Schematic of pulse sequence for dual 3D perfusion imaging. The pulse
sequence consists of two sets of saturation RF pulse followed by a delay
(T
SR
) and 3DFT GRE readout. The centers of two data acquisitions were
located at end-systole and mid-diastole. The largest T
SR
was used with a
maximum threshold of 120 ms, depending on a given heart rate. . . . . . . 77
5.2 Representativesystolic3Dperfusionimagesfromahealthyvolunteeratpre-
contrast, RVenhancement, LVenhancement, andmyocardialenhancement.
The passage of contrast agent is clearly captured with complete myocardial
coverage. The images suffer slightly from blurring and dark-rim artifacts
presumably due to lower in-plane spatial resolution. . . . . . . . . . . . . . 81
5.3 Representative diastolic 3D perfusion images from the same volunteer as in
Fig.5.2atpre-contrast,RVenhancement,LVenhancement,andmyocardial
enhancement. Dark rim is thinner than that in systolic images, but its ratio
to myocardial thickness is comparable to the ratio in systolic images . . . 82
xii
5.4 Representative time intensity curves based on (a) systolic and (b) diastolic
perfusion images from a healthysubject. Homogeneous signal enhancement
is observed in nearly all myocardial segments. Diastolic slice #3 suffered
from flickering artifacts caused by errors in coil sensitivity map estimation
(see gray arrow). Some endo-cardial TICs exhibit a signal dip (see black
arrows) that is temporally aligned with peak LV enhancement, and is likely
to be related to dark rim artifact. . . . . . . . . . . . . . . . . . . . . . . . 83
5.5 Correlation coefficients between systolic and diastolic TICs. Mean and
standard deviation of the coefficient were 0.9841 and 0.0166 across a to-
tal of 240 segments in five subjects (48 segments per subject), and mean
coefficient was larger than 0.97 in all five subjects, which validates strong
similarity between overall dynamics in the two sets of TICs. . . . . . . . . 84
5.6 Scatter plot of time intensity curve (a) upslope and (b) peak enhancement
for diastolic and systolic acquisitions. Linear fitting across all myocardial
segments results in y=0.9768x (P < .001) and y=0.9728x (P < .001) for
upslope and peak enhancement, respectively. . . . . . . . . . . . . . . . . . 86
5.7 Comparisonoftimecurveupslopeandpeakenhancementfromendo-cardial
andepi-cardiallayers: (a)systolicupslope, (b)diastolicupslope, (c)systolic
peak enhancement, and (d) diastolic peak enhancement. There was no
significant difference in either upslope or peak enhancement between the
endo- and epi-cardial layers. This was true for both systolic and diastolic
data across all subjects (at the significance level = 0.05). The ratio of
endo- to epi-cardial upslopes was 1.0539 and 1.0741 in the systolic and
diastolic TICs, respectively. The ratio of average endo- to epi-cardial peak
enhancement values was 0.9557 and 0.9580 in the systolic and diastolic
TICs, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.1 Schematicofpulsesequenceforhighlyaccelerated3Dperfusionimagingvia
TV regularized SENSE. The pulse sequence consists of saturation RF pulse
followed by a delay (T
SR
) and eight-fold undersampled 3DFT readout.
The center of k-spacewas always sampled, and the remaining region was
undersampled randomly, but in a disjoint way across nine cardiac cycles
for coil map estimation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.2 Flow chart of sampling function design. The center of k-space is initially
chosen, and is included in all nine sampling functions. For i
th
sampling
function, sampling points are randomly chosen from k-space positions that
are not included in 1
st
-(i¡1)
th
sampling functions. The final 9
th
sampling
function is automatically determined as all residual k-space points. N
cent
= number of samples in k-space center. N
outer
= number of samples in
outer k-space. S
cent
= k-space center that is always sampled. S
(i)
outer
=
outer k-space in i
th
sampling function. S
(i)
= i
th
sampling function. . . . 95
xiii
6.3 Point spread function of the nine chosen sampling function with (a) normal
and (b) harsh display window. The energy of sidelobe in all nine PSFs is
evenly diffused over the FOV, which improves the incoherence of aliasing
caused by k-space undersampling. . . . . . . . . . . . . . . . . . . . . . . . 96
6.4 Representative3Dperfusionimagesfromahealthyvolunteeratpre-contrast,
RV enhancement, LV enhancement, and myocardial enhancement. The
passage and arrival of contrast agent is clearly captured with high SNR.
Compared to the images in Chapter4, blurring is reduced due to the in-
creased spatial resolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.5 Effect of regularization coefficient on total variation and data fidelity. Ten
incremental values of the regularization coefficient was used, ranging from
0.00001 to 0.01. An increase in the coefficient reduces TV while sacrificing
data fidelity. The L curves suggest using 0.0005 with which TV is minimal
while the loss of data fidelity is negligible. The overall shape of the curve
is similar in the two data sets. . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.6 Effectofregularizationcoefficientonimagereconstruction. Eachrowshows
fourrepresentativepartitionslicesreconstructedusingdifferentvalueofreg-
ularizationcoefficients. Theimagesfromtoosmallregularizationcoefficient
suffer noisy aliasing artifacts while the images from too heavy regulariza-
tionsufferover-smoothing. ¸=0.0005produceswellbalancedimages, which
agrees with the findings from L-curve analysis. . . . . . . . . . . . . . . . 102
6.7 Representative myocardial time intensity curves (TIC) from a healthy vol-
unteer. TIC was generated from the images reconstructed using ¸=0.0005.
Uniform signal enhancement is seen across all myocardial segments. . . . 103
A.1 PSF profiles of the original concentric rings and the shifted rings by half of
radial spacing. Two profiles have the same magnitudes over whole re gion.
Phases of two profiles are alike at mainlobe and the second sidelobe but
diffrent from each other at the first sidelobe by ¼ rad. . . . . . . . . . . . 119
xiv
Abstract
Magnetic resonance imaging (MRI) is one of the most powerful non-invasive imaging
modalities whose strengths include flexible tissue contrast, sensitivity to flow, and the
lack of any harmful radiation. Ischemic heart disease (IHD) is a leading cause of death
in the western world, and as such, MR approaches for the assessment of IHD have gained
considerable attention. This dissertation will describe two improvements in cardiac MRI
methodology that are aimed at improving the assessment of IHD: (i) Real-time cardiac
imaging and (ii) First-pass myocardial perfusion imaging.
The first method, real-time cardiac MRI continuously captures the beating heart with-
out the need for cardiac gating and breathhold, and is a promising approach for the as-
sessment of ventricular function. Efficient k-space sampling trajectories (e.g., spiral or
echo-planar) are typically used to reduce acquisition time, but further improvement in
spatio-temporal resolution is desirable. A novel reconstruction method was developed for
accelerating real-time spiral imaging, based on an algebraic reconstruction framework that
can be applied to undersampled data. The method utilizes knowledge of the aliasing pat-
tern caused by spiral undersampling, and the fact that temporal bandwidth of dynamic
cardiac images is a function of spatial position.
xv
The second method, first-pass myocardial perfusion imaging (MPI) is a tool for early
detection of IHD by evaluating blood flow to the myocardium. The most widely used
protocol for MR MPI involves a 2D multi-slice acquisition, which has demonstrated high
sensitivity and specificity, but suffers from incomplete spatial coverage. We developed
a three-dimensional MR MPI approach that allows contiguous and complete coverage of
the left ventricle. We demonstrate its accuracy in a cardiac phantom, and its feasibility
in several healthy volunteers. We also compare the performance of diastolic acquisition,
which has a longer quiescent period and enables higher spatial resolution, versus systolic
acquisition, which is less sensitive to heart-rate variation or arrhythmia. Achieving both
complete coverage and high spatial resolution is challenging in 3D MPI. The last part of
this dissertation focuses the use of regularized reconstruction, and how it may be used to
improve spatial resolution.
xvi
Chapter 1 :
Introduction
Magnetic resonance imaging (MRI) has emerged as a powerful tool for the detection and
evaluation of cardiovascular disease because it is non-invasive and involves no ionizing
radiation, can provide many types of tissue contrast, is sensitive to flow, and can examine
arbitrary scan volumes. MRI has already been established as the gold standard modality
for examining left ventricular structure and function, and the presence of myocardial
scar [18]. Active research continues to advance MR in myocardial perfusion imaging,
valvular imaging and coronary artery imaging [51,73,78]. Many believe that MRI will
become a primary imaging modality for the comprehensive cardiac examination in the
near future [50].
Ischemic heart disease (IHD) is a leading cause of mortality in the world, account-
ing for over 7 million deaths in 2002 and is predicted to cause over 9 million deaths in
2030 [60]. IHD is also the number one cause of death in the United States, and was re-
sponsible for more than 4 million mortalities in 2005, and $165.4 billion direct and indirect
cost estimated in 2009 [63]. Figure 1.1 shows the series of events in IHD, beginning with
blockage (stenosis) of coronary arteries [66]. There are a variety of available treatments,
1
Blood flow reduction
Impaired diastolic
relaxation
Impaired systolic
function
ECG change
Angina pectoris
Metabolic change
Coronary Atherosclerosis
Time
Risk of mortality
Myocardial perfusion
Imaging (MPI)
Ventricular function imaging
Figure 1.1: Ischemic cascade in the heart. The narrowing or blockage of coronary artery
causes reduction in blood flow in the myocardium, followed by metabolic change, impair-
ment of diastolic and systolic function, and ECG change. This dissertation focuses on
MRI methods for detecting impaired cardiac function (ventricular function imaging), and
blood flow reduction (myocardial perfusion imaging)
.
involving lifestyle (diet, exercise, smoking, stress), pharmaceuticals, and surgical interven-
tion. However, early detection is critical to reduce mortality. This dissertation focuses on
MRI methods for assessing ventricular function and myocardial perfusion.
Ventricular function study is already an established MR examination, providing ex-
cellent blood-myocardium contrast using balanced steady-state free precession (SSFP)
sequence [11,30,68]. The current protocol involves a segmented data acquisition which is
synchronized by electrocardiogram (ECG) or plethysmographic gating [52,98]. This so-
called cine imaging requires a good breathhold and an assumption that cardiac motion is
purely periodic over many cardiac cycles (typically >10 cycles). However, subjects often
exhibit aperiodic patterns, and may also have difficulty in maintaining breathhold long
2
enough. Moreover, rapid switching of gradient pulses and radio frequency (RF) power
may corrupt the acquired ECG signal, which results in data inconsistency artifacts [21].
Real-time imaging continuously captures a dynamic object as it is without the need
for cardiac gating [43,65,106]. Thus, real-time imaging is free from the breathhold re-
quirement and the assumption on periodic cardiac motion. Typically, real-time imaging
requires time-efficient k-space trajectories to reduce acquisition time and thus avoid mo-
tion artifacts. Spiral and echo-planar readouts are popular examples of such trajectories,
but further reduction in data acquisition time is demanded for higher spatial and tempo-
ral resolution in many applications. The first part of this dissertation deals with image
reconstruction from undersampled real-time spiral imaging to improve spatio-temporal
resolution.
Myocardial perfusion imaging is a non-invasive way to assess blood flow into the my-
ocardium, which enables early detection of IHD [16,24,32,62,86,88]. MPI is a promising
alternative to invasive X-ray coronary angiography that is the current gold standard for
diagnosis of IHD. Presently, single photon emission computed tomography (SPECT) is
routinely used in clinical myocardial perfusion studies, but suffers from low spatial reso-
lution, attenuation artifacts, and exposure to ionizing radiation.
Magnetic resonance imaging is an ideal modality for myocardial perfusion imaging
(MPI), due to high spatial resolution, lack of radiation exposure, and relatively short scan
time. The most popular protocol of MR MPI involves 2D multi-slice acquisition which
enables the acquisition of 3-6 slices per one or two cardiac cycles [40]. While this 2D multi-
slice approach has demonstrated high sensitivity and specificity [4,75,104], its incomplete
spatial coverage leaves room for potential improvements in diagnosis and prediction of
3
outcomes via increased coverage. Three-dimensional MR MPI is a promising alternative
to2Dmulti-sliceMPIduetoitscontiguousspatialcoverageandhighSNR,yetsuffersfrom
long data acquisition time. The second part of this dissertation deals with the feasibility
and acceleration of whole heart 3D MR MPI.
1.1 Outline
The structure of the dissertation is as follow:
Chapter 2: MRI Background
This chapter contains an overview of MR imaging concepts, including MR physics, acqui-
sition and reconstruction principles. Three types of undersampling techniques, parallel
imaging, temporal filtering and compressed sensing, are also explained.
Chapter 3: Accelerated Spiral Imaging of Ventricular Function
A reconstruction method from undersampled dynamic spiral data is presented. Based
on the description of the aliasing from spiral undersampling, an iterative reconstruction
method is proposed, and validated using numerical simulation and in-vivo experiments.
Chapter 4: Three-Dimensional Myocardial Perfusion Imaging: Feasibility
The development of 3D first-pass MPI at 3 Tesla is presented, including pulse sequence,
image reconstruction and analysis. The feasibility of 3D MPI with complete LV coverage
is demonstrated in healthy volunteers.
4
Chapter 5: Three-Dimensional Myocardial Perfusion Imaging: Systole vs Di-
astole
Systolic 3D MPI is developed as an alternative to diastolic MPI to reduce sensitivity to R-
R variability. Semi-quantitative perfusion indices from systolic and diastolic acquisitions
are compared in healthy volunteers.
Chapter 6: Three-Dimensional Myocardial Perfusion Imaging with Regular-
ization
Highly accelerated 3D MPI is developed via total variation regularized sensitivity encod-
ing, and its feasibility is tested in healthy volunteers.
Chapter 7: Summary and Future Work
Contributions of this dissertation are summarized, and recommendations for future work
are presented per each topic.
5
Chapter 2 :
MRI Background
This chapter contains a brief description of the basic MRI principles. It will also provide
introduction to a few advanced MRI reconstruction methods, which are preliminaries for
the remaining chapters.
2.1 Basic Principles of MRI
2.1.1 MR Physics
Nucleus with odd number of proton and/or odd number of neutron possess angular mag-
netic moment, and create electromagnetic field, as discovered by Bloch [6,81]. Such nuclei
are referred to as "spins", and their magnetic moment are represented by a magnetization
vector for convenience (Fig. 2.1). The hydrogen nucleus has a single proton, and is the
most frequently imaged in MRI because of abundant H
2
O in human body. Another pop-
ular atoms include
31
P,
31
C, and
19
F which are useful for studies on cellular metabolism.
MR images are formed based on the interactions between “spins” and three external mag-
netic field components: strong static field B
0
, radio frequency (RF) field B
1
(t), and linear
6
B
0
a b
Figure 2.1: Nuclear spins with and without a static field. (a): In the absence of static
field B
0
, spins align in a random way such that the net magnetization is zero. (b): In the
presence of B
0
, spins are lined up either in the direction of B
0
or in the reverse direction.
gradient fieldsG(t). In absence of static field B
0
, the axes of spin magnetic moments are
arranged in a random way so that their net moment is zero. When B
0
is applied, two
phenomena are observed. First, all spins are lined up either in the direction of B
0
or in
the reverse direction. The difference in spins aligned in the two directions will create net
magnetic moment, which is called polarization (Fig. 2.1). This phenomenon is similar to
the alignment of compass needle by the earth’s magnetic field. Second, each magnetic spin
not only spins around its own magnetization axis, but also precesses around the axis of
B
0
. The angular velocity of the precession around B
0
is called Larmor frequency, written
as
!
0
=°B
0
: (2.1)
where ° is called the gyromagnetic ratio, and is a nucleus specific constant. For example,
° = 4.26 kHz/G for hydrogen.
The purpose of the RF field B
1
(t) is to tip the longitudinal magnetization at equilib-
rium to produce a transverse component, which is detectable by a receiver coil. If B
1
(t)
7
B
0
x
y
Figure 2.2: Precession of a nuclear spin around static field. In the presence of B
0
spins
not only spin around their own magnetization axes, but this axis also precesses around
the direction of static field B
0
.
tuned to the Larmor frequency, is applied in the transverse plane, the magnetization not
only precesses around the axis of B
0
field, but also rotates around the axis of B
1
(t) field.
The behavior of magnetization vector in the presence of B
0
and B
1
(t) looks like spiral
motion in the absolute coordinate system (Fig. 2.3(a)). If we look at the same motion in
a coordinate system rotating at the Larmor frequency, the motion will be simple rotation
around the axis of B
1
(t) field (Fig. 2.3(b)).
When only B
0
and B
1
(t) fields are applied, all magnetic spins in the range of the two
fields will experience global polarization and resonance. In practice, only a fraction of
whole object is of interest for MR imaging, and spatially selective excitation and encoding
are needed. Linear gradient field G(t) provides spatial selectivity, by generating mag-
netic fields linearly varying proportional to spatial location. Then spatial position can be
8
x
M
B
1
z
y
x’
M
B
1
z
y’
a b
Figure 2.3: Excitation of a nuclear spin by a RF field. RF field B
1
(t) tuned to the
resonance frequency, tips down spins from the equilibrium state. This behavior is seen as
(a): spiral motion in the absolute coordinate, and (b): as simple rotation around the axis
of B
1
(t) in the relative coordinate rotating at the resonance frequency.
mapped, proportionally, to the resonance frequency by the Larmor equation (2.1). Utiliz-
ing spatial variation in resonance frequency for selective excitation and spatial encoding
will be explained in the subsequent subsections.
Once spins are tipped down by B
1
(t) field, they start going back to the original
equilibrium state as in Fig. 2.1(b), which is called relaxation. Since the equilibrium state
contains purely longitudinal magnetization, any reduced longitudinal component of a spin
will grow up to the equilibrium level, and any increased transverse component will decay
to zero. Longitudinal relaxation is described by,
M
z
=M
0
+(M
z
(0)¡M
0
)e
¡t=T
1
(2.2)
9
where T
1
, M
0
, M
z
(0) are the time constant of longitudinal relaxation, longitudinal mag-
netization at the equilibrium, and longitudinal magnetization at the initial state, respec-
tively. Transverse relaxation is described by,
M
xy
=M
xy
(0)e
¡t=T
2
(2.3)
where T
2
and M
xy
(0) are a time constant of transverse relaxation and initial transverse
magnetization. ThetworelaxationtimeconstantsT
1
andT
2
aretissuespecificparameters,
and therefore, are critical parameters in determining tissue contrast in MRI. Since pulse
sequencedesigncanadjustthe degreeofdependenceofeachtimeconstantontheresulting
signal level, different tissue contrast can be obtained from the same subject.
The behavior of magnetization vector in the presence of an applied magnetic field B,
is described phenomenologically by the Bloch equation, written as,
dM
dt
=(M£°B)¡
M
x
^ x+M
y
^ y
T
2
¡
(M
z
¡M
0
)^ z
T
1
; (2.4)
where ^ x, ^ y, and ^ z are unit vectors in the x, y, and z directions respectively, andB includes
the various magnetic fields applied such as static field B
0
RF field B
1
(t) and gradient
field G(t).
2.1.2 Selective Excitation
As explained in the previous section, RF field B
1
(t) tuned to the resonance frequency
of nuclear spins will tip all spin magnetization. In this case, all spins have the same
resonancefrequency, andthereisnoselectivityduringexcitationprocess. Spatialencoding
10
of large 3D volume is time-consuming, and is impractical. Therefore, the excitation of
smaller region of interest (ROI) is desirable, which can be achieved by modifying two
magnetic field components. First, B
1
(t) whose baseline frequency is tuned to the Larmor
frequency, should have time-varying amplitude such that it contains limited bandwidth.
Then, only nuclear spins with resonance frequencies within the bandwidth, will be excited.
Second, linear gradient fieldG(t) is applied in the direction of desired spatial selectivity, z
direction for example. In the presence of G
z
, local resonance frequency !(z) is described
by °(B
0
+G
z
¢z), i.e., a function of spatial location z. As a result, only spins in a finite
region corresponding to the resonance bandwidth will be selectively excited.
2.1.3 Spatial Encoding
After the selective excitation explained in the previous section, all magnetic spins within
the slice will contribute to the net signal captured by receiver coils. The net signal can
not be decomposed based on the spatial position, and a question about how to encode
each position arises. MRI performs spatial encoding by applying Fourier kernel over the
spin population.
The transverse magnetization within a voxel can be modeled as a complex number
with some magnitude m(r) and with some phase Á(r;t). Here, r represents the spatial
position and the magnitude is assumed to be time-invariant. The signal from each voxel
at position r can be expressed as m(r)e
¡iÁ(r;t)
. Then the total received signal is the
summation of transverse magnetizations within the plane of interest.
s(t)=
Z
m(r)e
¡iÁ(r;t)
dr (2.5)
11
ω ω
0
z
B (t)
1
Fourier transf.
( ) z G B
z
⋅ + =
0
γ ω
Figure 2.4: Selective excitation using a RF field and gradient fields. Amplitude modifi-
cation of B
1
(t) field creates limited bandwidth such that only spins within the band can
be excited. With linear gradient turned on, resonance frequency of spins linearly changes
across spatial location. As a result, spins in a finite spatial region are selectively excited.
The slice thickness is increased by either decreasing G
z
or increasing the bandwidth of
B
1
(t)
From the Larmor equation ! =°B, the phase of magnetization can be written as
Á(r;t)=
Z
t
0
!(r;¿)d¿ =°
Z
t
0
B(r;¿)d¿ (2.6)
In the presence of both B
0
and gradient,
B(r;t)=B
0
+G(t)¢r (2.7)
12
Then the signal equation is written as
»(t)=
Z
m(r)e
¡i!
0
t
e
¡i°
R
t
0
G(¿)¢rd¿
dr (2.8)
After demodulation of the baseline frequency, we get
s(t) =
Z
m(r)e
¡i2¼[k(t)¢r]
dr
= M(k(t)) (2.9)
where k is defined as
k(t)=
°
2¼
Z
t
0
G(¿)d¿ (2.10)
and M(k) is the Fourier transform of m(r). Equation (2.9) states that the received signal
can be interpreted as Fourier transform of transverse spin magnetization at a specific
spatial frequency, which is determined by the time integral of gradient. In other words,
gradient waveforms can be designed in such a way that Fourier domain can be traversed
as we want. By convention, the Fourier domain M(k) is called k-space.
The most commonly used k-space trajectory is 2DFT readout which covers k-space in
a raster-scan fashion (Fig. 2.5(a)). The 2DFT trajectory conforms to the Cartesian grid,
which makes the reconstruction simple, and is insensitive to errors due to off-resonance
and gradient delay. The drawback of 2DFT readout is the large number of encoding
lines required, which increases acquisition time. Several non-Cartesian trajectories have
been designed, which cover k-space in more time-efficient way. One example called spiral
trajectory, is shown in Fig. 2.5(b). Larger portion of k-space can be traversed using
13
single spiral interleaf, which reduces the required number of encoding lines for fullk-space
coverage, and therefore decreases total imaging time.
2.2 Advanced Reconstruction of MR Images
2.2.1 Non-Cartesian Reconstruction
Non-Cartesiank-spacetrajectoriessuchasspiralandechoplanarenablerapidMRimaging
due to efficient k-space coverage. In contrast to conventional 2DFT readout, the fast
Fourier transform (FFT) cannot be directly applied to k-space data that do not fall on a
Cartesian grid. Non-uniform Fourier transform or its approximate versions are required,
but are computationally heavy [2,22].
Gridding is the most widely used method for MRI reconstruction from non-Cartesian
data [37,69,82,87]. The non-Cartesian samples are interpolated onto a Cartesian grid,
followed by the conventional FFT. Gridding with a Kaiser-Bessel kernel provides high
interpolation accuracy at reasonable computational cost [37,100]. Each data point on the
non-Cartesian grid is convolved with a gridding kernel, and the convolution result is re-
sampled on the Cartesian grid (Fig. 2.6). This procedure can be mathematically written
as,
^
M(k)=[fM(k)S(k)g¤C(k)]comb
µ
k
¢
k
¶
(2.11)
where M(k) is the true k-space data, S(k) =
P
i
±(k¡k
i
) is a non-Cartesian sampling
function, C(k) is a gridding kernel, and comb(x)´
P
i
±(x¡i) is the Cartesian sampling
14
−4
−2
0
2
4
G
x
[G/cm]
0 2 4
−4
−2
0
2
4
t [ms]
G
y
[G/cm]
−4
−2
0
2
4
G
x
[G/cm]
0 1.5
−4
−2
0
2
4
t [ms]
G
y
[G/cm]
a
b
a b c
a
b c
a
a
b
b
c
c
Figure 2.5: Two examples ofk-space trajectories. k-space is traversed by the time integral
of gradients. (a): 2DFT trajectory is the most commonly used, where k-space is scanned
in a raster fashion. (b): spiral trajectory, one of the non-Cartesian trajectories, can cover
k-space in more time-efficient way.
15
x
k
Figure 2.6: Gridding of one-dimensional non-Cartesian data. Data from each non-
Cartesian sample denoted by “x” are convolved with a kernel, and re-sampled onto neigh-
boring Cartesian grids denoted by circles. In the image domain, this procedure is equiv-
alent to convolution with an impulse response and apodization by the transform of the
gridding kernel.
function. The effect of the sampling function, convolution with the gridding kernel, and
re-sampling on the Cartesian grid can been seen in the reciprocal image domain as,
^ m(x)=[fm(x)¤s(x)gc(x)]¤comb(¢
k
x) (2.12)
The true image is convolved with an impulse response (inverse Fourier transform of
the non-Cartesian sampling function), and is weighted by apodization function (inverse
Fourier transform of the gridding kernel). The impulse response should ideally be an
impulse but is not due to uneven density of the non-Cartesian samples. To correct for
the impulse response distortion, the non-Cartesian data is pre-multiplied by the sample
density which can be estimated geometrically or numerically [34,74]. The weighting by
apodization function can be also corrected for by dividing the resulting image by the
transform of the apodization function.
16
2.2.2 Parallel Imaging
Parallel imaging requires the use of receiver coil arrays, and exploits sensitivity differences
between individual coil elements. Simultaneous acquisition of MR data modulated by
different coil sensitivities provides additional capacity for spatial encoding, which reduces
the required number of gradient encoding steps. The saving of gradient encoding can
reduce scan time and/or increase spatial resolution. Many parallel imaging methods have
been developed in the past decade, and can be divided into k-space based approaches and
image based approaches.
Generalized auto-calibrating partially parallel acquisitions (GRAPPA) is the most
widely used k-space based method, and extracts sensitivity information from obtained
k-space samples to estimate missing k-space data [25]. GRAPPA reconstruction is sub-
optimal in the sense that it uses only partialk-space to extract coil sensitivity information.
However, the self-calibrating nature of GRAPPA makes it robust in some applications
where coil sensitivity estimation is difficult.
Sensitivity encoding (SENSE) is based on the exact description of spatial encoding by
both coil sensitivities and gradients in image domain [80]. With accurate coil sensitivity
maps, SENSE provides the best possible reconstruction performance [5]. Main limita-
tion of SENSE is its reliance on accurate knowledge on coil sensitivity maps, which is
difficult particularly in abdominal and dynamic imaging. Adaptive sensitivity encoding
incorporating temporal filtering (TSENSE) is a variation of SENSE suited for dynamic
imaging [41]. Dynamick-space data are undersampled in time interleaved fashion, and are
temporally low-pass filtered to form full FOV images used to estimate coil sensitivities.
17
In the following, we explain the basic principle of SENSE that will be used in subsequent
chapters 4-6.
Fig. 2.7 shows how differences in coil sensitivities can be utilized in an example of two-
fold undersampling. The true image is weighted by the coil sensitivity, and is convolved
with the point spread function (PSF) from two-fold undersampling, which consists of main
lobe at the center, and sidelobe at the top and the bottom. For simplification, we focus on
two image pixels which are aliased together, denoted by x
1
and x
2
. Sensitivity weighting
and aliasing in each coil element can be described as,
r
k
=s
k1
x
1
+s
k2
x
2
(k =1;2;¢¢¢ ;4) (2.13)
where r
k
are aliased pixels, and s
k
l represent sensitivity of k
th
coil element. Applying the
similar relation to all aliased pixels, we get
Sx=r (2.14)
where S is the encoding matrix that represents sensitivity weighting and aliasing, and
x and r are true image, and aliased images in coil elements, respectively. The inversion
problem can be solved using conventional least-square criterion.
2.2.3 Temporal Filtering
Temporalcorrelationinamovingobjectisanothersourceofredundancytobeexploredfor
the reconstruction from undersampled k-space data. Early examples of such approaches
18
GJX GJY GJZ GJ[ G G U U U U U U U U U U U U U U X Y XX XY YX YY ZX ZY [X [Y X Y Z [ Q U U U GG Figure 2.7: Pictorial description of multi-coil acquisition in Cartesian two-fold undersam-
pling. The true image is weighted by the sensitivity profile of individual coil element, and
is convolved by a point spread function.
19
t=0 t=1 t=2
UUU x x x
f
Nyquist
DC
f
-f
Nyquist
a b
Figure 2.8: UNFOLD method in two-fold Cartesian undersampling. UNFOLD acquires
odd-order and even-order k-space encoding lines alternately over time. (a): The phase of
the aliasing will alternate, and (b): the aliasing will be shifted to the Nyquist frequency
by a Fourier transform along temporal dimension. The true signal can be separated from
the aliasing by filtering.
include keyhole methods [39,99], reduced-encoding MR imaging by generalized-series re-
construction (RIGR) [101], and reduced field-of-view (FOV) reconstruction [10,35,90].
Recently, the concept of utilizing temporal correlation was advanced by UNaliasing by
Fourier-encoding the Overlaps using the temporaL Dimension (UNFOLD) [58,96] and k-t
Broad-use Linear Acquisition Speed-up Technique (k-t BLAST) [97].
UNFOLD undersamplesk-space in a time-interleaved fashion (e.g., odd order and even
order encoding lines for two-fold undersampling case). Consequently, aliasing will occur
along the phase encoding direction in the reconstructed images, but the phase of aliasing
will alternate over time (Fig. 2.8(a)). If each spatial point is Fourier transformed along
the time axis, the aliasing will be shifted to the Nyquist frequency by the modulation
property (Fig. 2.8(b)). Since the spectral support of the true signal is finite, the overlap
of the two signals will be small enough to be separated by filtering.
In addition to the time-interleaved undersampling, k-t BLAST uses training data as
a priori information on the solution. Low resolution images are obtained from a separate
pre-scan or during an actual scan using variable density k-space sampling, and serve as
20
training data. The signal power in the training image provides an approximation of true
signal distribution, which facilitate resolving aliased voxels. Mathematically, k-t BLAST
reconstruction looks for the minimum norm solution weighted by the squared magnitude
of training image.
x=M
2
E
H
¡
EM
2
E
H
¢
¡1
½ (2.15)
where M is the magnitude of training data, E is the encoding matrix obtained from
the PSF of a given undersampling pattern, and ½ is the aliased signal. Note that in
k-t BLAST, all images are first Fourier-transformed along the temporal dimension, and
therefore operations in (2.15) are in the hybrid domain of image space and temporal
frequency.
2.2.4 Compressed Sensing
Transform sparsity of an image itself is another potential source of redundant information
to be explored for higher acceleration. MR images, like any other natural images, are
sparse in terms of the number of non-zero coefficients in an appropriate transform domain.
For instance, images from contrast enhanced angiography will be sparse in the image
domain since only branches of blood vessel are bright amidst a background of nulled
tissue signal. If a certain image is ideally piece-wise constant, its finite difference will have
only a few non-zero values.
The concept of utilizing image sparseness has inspired theoretical studies on the recon-
struction from insufficient data samples, which was named compressed sensing (CS) [17],
[12]. One remarkable result of CS theory states that [12,95]: For an unknown discrete
21
and sparse object f(t), t=t
1
;¢¢¢ ;t
d
, and Fourier transform known at a randomly chosen
incomplete set of frequencies F(!)j
!2
, f(t) can be recovered by solving the following
optimization
^
f(t)=argmin
g
kªgk
l
1
subject to G(!)=F(!)j
!2
(2.16)
whereF(!)andG(!)aretheFouriertransformsoff(t)andg(t)respectively. represents
the set of randomly undersampled frequencies. ª andk¢k
l
1
denotes an appropriate sparse
transform and l
1
norm, respectively.
Two requirements for the successful application of CS to MRI should be emphasized.
First, the image to be reconstructed should be transform sparse, and finding the transform
providing the largest sparseness is necessary. Finite difference and wavelet transform
are mostly widely used transforms, and the effectiveness of each transform depends on
specific applications. Second, k-space should be undersampled in a random way, so that
the aliasing caused by the undersampling is incoherent. The incoherence requirement
is important to easily distinguish the sparse true signal from the aliasing (Fig. 2.9). A
perturbed spiral trajectory, radial trajectory, and variable density Cartesian trajectory
have been shown to produce incoherent aliasing [7,56,57].
22
a b
c
Figure 2.9: Importance of incoherent aliasing in compressed sensing reconstruction. (a):
A Shepp Logan phantom (top) has only a few non-zeros total variation (TV) coefficients
(bottom). (b): Ifk-space is uniformly undersampled (top), a true TV coefficient is not dis-
tinguishable from its alias since its contribution to the aliased signal is the same (bottom).
(c): If k-space is randomly undersampled (top), large TV coefficients are distinguishable
from their noise-like aliases (bottom). Smaller TV coefficients can be iteratively detected
by subtracting the aliased signal by the contribution of the detected coefficients.
23
Chapter 3 :
Accelerated Spiral Imaging of Ventricular Function
3.1 Introduction
Dynamic magnetic resonance imaging (MRI), which captures temporal variations of ob-
jects of interest, is receiving increasing interest as an effective diagnostic tool. It has
been proven to be particularly useful for cardiac imaging since the human heart involves
a high degree of dynamic motion [73], [19]. Cardiac gated scanning during a breathhold
can achieve high spatial and temporal resolutions, but is affected by arrhythmias and the
breathhold requirements. In some cases, the aperiodic pattern itself may represent actual
diagnostic information. Recent advances in MR hardware and pulse sequences have made
it possible to perform dynamic cardiac imaging without gating or breath holding [65].
However, many applications still demand higher temporal resolution, which has remained
an important research topic.
One way to increase the temporal resolution of dynamic MRI is by sampling only
a fraction of k-space. The fundamental principle of these methods is to exploit infor-
mation redundancy, and they can be divided into two groups based on the source of
24
the redundancy. The first group includes partial k-space reconstruction [67] and parallel
imaging [80], [92], [25] where the redundancy source is inherent properties of the image or
external hardware. Dynamic properties of moving objects are irrelevant to these methods,
and they can be applied to both static and dynamic imaging. The second group utilizes
information redundancy from temporal correlations of dynamic images. Keyhole methods
and reduced-encoding MR imaging by generalized-series reconstruction (RIGR) assume
that temporal variation is primarily from the central portion of k-space [99], [39], [101].
Similarly, reduced field of view (FOV) reconstruction utilizes the fact that only a lim-
ited region of an image undergoes substantial variation over time [35], [90], [10]. That
is, the information from relatively stationary regions in object space is considered to be
redundant.
The second group, introduced above, was advanced by a novel reconstruction based on
temporal filtering called UNaliasing by Fourier-encoding the Overlaps using the temporaL
Dimension (UNFOLD) [58], [96]. By undersamplingk-space in an interleaved fashion over
time, UNFOLD displaces aliased signals along the frequency axis and recovers the original
signal using a temporal low-pass filter. Further acceleration is achievable in k-t Broad-use
Linear Acquisition Speed-up Technique (k-t BLAST), which uses training data as prior
information along with the interleaved sampling function [97]. Conceptually, UNFOLD
and k-t BLAST reduce the information redundancy by economically using k-t space (the
combination of k-space with a time axis [105]). Until now, these two methods have been
used mainly with Cartesian k-space trajectories, due to the simplicity of describing and
correcting for the aliasing produced by recti-linear sampling patterns (refer to [27] as a
25
non-Cartesian example). However, spiral readouts are appropriate for certain applica-
tions because they have good flow properties [61], and are time efficient, which enables
further increases in spatial and/or temporal resolution when used in conjunction with
undersampling techniques.
In this chapter, we describe the aliasing process for spiral sampling trajectories and
formulate unaliasing as an inverse problem in 3Dxyf space, the transform ofk
x
k
y
t space.
Ingeneral, anyaliasingcausedbytheviolationoftheNyquistconditioncannotberesolved
under arbitrary condition, which manifests, in our case, as numerically under-determined
system matrices. We show that the aliasing is formed such that only a few temporal
frequency components are coupled, and rank deficiency problem is isolated in each set of
temporal frequencies. From the fact that the support region is bounded in xyf space,
a zero assumption in the solutions is used to condition each linear system. We inspect
how the size and the location of assumed zeros affect the degree of rank deficiency by
borrowing ideas from the Cartesian undersampling case. Finally, the proposed method is
validated through numerical simulations and in-vivo spiral cardiac MRI experiments.
3.2 Reconstruction Model
3.2.1 Aliasing in Dynamic Spiral Imaging
During dynamic MRI, k-space data of moving objects are acquired according to a specific
sampling function at different time points. The k
x
k
y
t space should be considered as the
26
working space of data acquisition in which a sampled object R
kxkyt
is expressed as an
original object I
kxkyt
multiplied by a sampling function S
kxkyt
I
kxkyt
(k
x
;k
y
;t)¢S
kxkyt
(k
x
;k
y
;t)=R
kxkyt
(k
x
;k
y
;t) (3.1)
In the reciprocal domain called xyf space, the above relation becomes
I
xyf
(x;y;f)¤
3
S
xyf
(x;y;f)=R
xyf
(x;y;f) (3.2)
where S
xyf
(x;y;f) = F
¡1
3
¡
S
kxkyt
¢
, the 3D inverse Fourier transform of the sampling
function, is a point spread function (PSF), and ¤
3
denotes the 3D convolution operator.
In the rest of this paper, all operations are in xyf space and the subscript xyf will be
omitted.
When a fraction of k-space is scanned in a time-interleaved fashion, the PSF sidelobe
has different phase at successive time points, with the pattern depending on the acquisi-
tion order. Aliased signals are shifted along the temporal frequency axis, which makes it
easier to separate the original signal from the aliased ones. The modulation phenomenon
of the PSF can be proved in a straightforward way in Cartesian undersampling [58], but
the rigorous proof of it in spiral sampling is challenging due to the difficulty in mathe-
matically handling spiral trajectories. Instead, concentric ring trajectories can be used as
an approximation model and the relevant proof is described in Appendix. An intuitive
proof, particularly in two-fold undersampling, can be provided as follow (Fig. 3.1). The
point spread function (PSF) from sampling even ordered or even ordered spiral interleaves
27
+ =
k
x
k
y
x
y
+ =
Full G Odd ordered G Even ordered G Figure 3.1: Proof of opposite sidelobe phases in the PSFs from interleaved spiral under-
sampling. The top row shows spiral sampling function and the bottom row shows the
corresponding point spread functions (PSF). Sampling odd ordered or even ordered spiral
interleaves results in ring-shaped sidelobe as well as mainlobe while the PSF from full
sampling contains mainlobe only. Since full sampling is the summation of the two under-
sampling functions and Fourier transform is linear, the two aliased PSFs should add up to
the PSF from full sampling. Therefore, the sidelobe phases of the two aliased PSF should
be opposite to each other.
consists of ring-shaped sidelobe as well as mainlobe while PSF from full sampling contains
mainlobe only. Since full sampling is the summation of the two undersampling functions,
and Fourier transform is linear, the summation of the two aliased PSFs should be equal
to the PSF from full sampling, e.g. mainlobe only. Therefore, the sidelobe phases of the
two aliased PSF should be opposite to each other.
The top of Fig. 3.2 shows an example object, the PSF associated with spiral two-fold
undersampling, and the corresponding aliasing pattern in 3D xyf space. It is assumed
that only the central part of the object is moving. Once the PSF is described in this 3D
space, the resulting aliasing can be understood as a 3D convolution of the original object
with its mainlobe and sidelobes.
28
Figure 3.2: Illustration of dynamic undersampled spiral imaging inxyf space. From two-fold undersampling, ring-shaped sidelobes
are shifted to the Nyquist frequency. Since the mainlobe and sidelobes are confined to DC and the Nyquist frequency respectively,
the solutions of only two temporal frequencies are coupled to generate aliasing at one temporal frequency.
29
3.2.2 Algebraic Formulation
The 3D convolution (3.2) can be formulated as a linear matrix system. However, solving
the inversion problem as a whole is impractical on a conventional workstation due to
memory requirements. Also, the sparseness of the PSF remains unexplored in the 3D
convolution formulation. Because of the temporal modulation of the PSF in two-fold
undersampling, the mainlobe and sidelobe are strictly confined to DC and the Nyquist
frequency, respectively, as shown in Fig. 3.2. Then we can express the 3D convolution as:
R(x;y;f) =
L¡1
X
l=0
N¡1
X
p=0
N¡1
X
q=0
[I(p;q;l)S(x¡p;y¡q;f¡l)]
=
L¡1
X
l=0
2
4
N¡1
X
p=0
N¡1
X
q=0
I
(l)
(p;q)S
(f¡l)
(x¡p;y¡q)
3
5
=
L¡1
X
l=0
h
I
(l)
¤
2
S
(f¡l)
i
= I
(f)
¤
2
S
(0)
+I
(f¡
L
2
)
¤
2
S
(
L
2
)
(3.3)
where I
(l)
and S
(l)
denote the original signal and the PSF at temporal frequency f = l,
and ¤
2
denotes 2D linear convolution. In the third line of (3.3), the initial 3D convolu-
tion is expressed as the summation of 2D convolutions. Furthermore, only two temporal
frequency terms (S
(0)
and S
(
L
2
)
) are non-zero, and the equation therefore reduces to the
last expression in (3.3).
30
If the 2D convolution matrices associated with S
(0)
and S
(
L
2
)
are A
1
and A
2
, the
aliasing from two-fold undersampling can be written as
2
6
6
4
A
1
A
2
A
2
A
1
3
7
7
5
2
6
6
4
x
l
x
l+
L
2
3
7
7
5
=
2
6
6
4
b
l
b
l+
L
2
3
7
7
5
l =0;1;¢¢¢ ;L=2¡1 (3.4)
where x
l
=
¡
I
(l)
¢
N£N!N
2
£1
and b
l
=
¡
R
(l)
¢
N£N!N
2
£1
, i.e., N
2
£ 1 column vectors in
which the elements ofI
(l)
andR
(l)
are stacked in a raster-like fashion. Note that only two
temporal frequenciesI
(l)
andI
(l+
L
2
)
are coupled into the same linear system. For example,
the aliasing at the Nyquist frequency is generated from the DC and the Nyquist frequency
components which are convolved with the mainlobe and the sidelobe, respectively (see
bottom of Fig. 3.2).
The generalization of this algebraic formulation into higher order undersampling is
straightforward. The PSF from m-fold undersampling contains non-zero signals at only
m temporal frequencies f = 0;
L
m
;¢¢¢ ;
m¡1
m
L. If we represent the corresponding 2D con-
volution matrices as A
1
;A
2
;¢¢¢ ;A
m
, then the generalized version of (3.4) is written as
2
6
6
6
6
6
6
6
6
6
6
4
A
1
A
m
::: A
2
A
2
A
1
::: A
3
.
.
.
.
.
.
.
.
.
.
.
.
A
m
A
m¡1
::: A
1
3
7
7
7
7
7
7
7
7
7
7
5
2
6
6
6
6
6
6
6
6
6
6
4
x
l
x
l+
L
m
.
.
.
x
l+
m¡1
m
L
3
7
7
7
7
7
7
7
7
7
7
5
=
2
6
6
6
6
6
6
6
6
6
6
4
b
l
b
l+
L
m
.
.
.
b
l+
m¡1
m
L
3
7
7
7
7
7
7
7
7
7
7
5
Ax
(l)
= b
(l)
l =0;1;¢¢¢ ;
L
m
¡1 (3.5)
31
The solution of every
¡
L
m
¢
th
frequency is coupled into the same linear system, and a total
of
L
m
such systems can be solved independently.
In general, the aliasing described by (3.2) cannot be easily resolved since the cor-
responding linear equations are ill-posed. The separability of the solution in (3.4) and
(3.5) indicates the insufficient rank problem is isolated in each set of aliased temporal
frequencies. The next subsection discusses a method to condition each set of equations.
3.2.3 Rank Analysis and Zero Assumption
The existence of a unique solution to (3.5) is directly related to the rank of the system
matrixA. The solution can be uniquely determined only if the rank of the system matrix
is equal to the number of unknowns. The direct rank analysis of the system matrix
associated with spiral undersampling is quite complicated due to the complexity of its
PSF. To obtain intuition about rank deficiency, we will borrow ideas from Cartesian two-
fold undersampling and show that the degree of rank deficiency depends on the overlapped
area between the non-zero signals at two aliased temporal frequencies.
Fig. 3.3 shows the image view of the system matrix defined in (3.4) associated with the
Cartesian two-fold undersampling and the corresponding solutions I
(l)
and I
(l+
L
2
)
. A
1
is
diagonal andA
2
has two off-diagonal lines sinceS
(0)
has one impulse at the object domain
origin and S
(
L
2
)
has two impulses at the edge in the phase encoding direction. The rank
ofA is exactly half the matrix dimension. Assuming I
(l)
has the size of N£N and thus
the dimension of A is 2N
2
£2N
2
, the first quarter columns of dimension 2N
2
£
N
2
2
are
the same as the last quarter columns, and the second quarter columns are the same as the
third quarter. These repeated columns indicate that the overlap between the upper (or
32
Figure 3.3: Illustration of the system matrix in Cartesian two-fold undersampling with
the solution for two coupled temporal frequencies. The rank is half the matrix dimension
since the first and fourth quarter-columns (gray) are identical, and the second and third
quarter-columns (white) are identical. To condition rank deficiency, zeros can be assumed
in certain regions of the solution. P
k
and Q
k
(k =1;2) represent the region of zero
assumption. Accordingly,P
k
c
andQ
k
c
(k =1;2)representsupportregionsofthesolution.
lower) half of I
(l)
and the lower (or upper) half of I
(l+
L
2
)
occurs in both R
(l)
and R
(l+
L
2
)
.
Consequently, any aliased position forms a pair with another location to produce single
effective rank.
A good candidate for mitigating this rank deficiency is assuming zeros in a certain
region of the solution since the solution’s support region is bounded in xyf space. Then
the columns of the system matrix corresponding to the position of zero assumption will be
discarded. A key question about this assumption is the effect of the size and the location
of assumed zeros on the rank deficiency. For notational convenience, we define four sets
33
of zero locations P
k
and Q
k
(k =1;2) and four sets of the corresponding support regions
P
k
c
and Q
k
c
(k =1;2) (Fig. 3.3).
P
1
´
½
(x;y)jI
(l)
(x;y)=0;0·y·
N
2
¡1
¾
P
2
´
½
(x;y)jI
(l+
L
2
)
(x;y)=0;
N
2
·y·N¡1
¾
Q
1
´
½
(x;y)jI
(l)
(x;y)=0;
N
2
·y·N¡1
¾
Q
2
´
½
(x;y)jI
(l+
L
2
)
(x;y)=0;0·y·
N
2
¡1
¾
(3.6)
P
1
c
´
½
(x;y)jI
(l)
(x;y)6=0;0·y·
N
2
¡1
¾
P
2
c
´
½
(x;y)jI
(l+
L
2
)
(x;y)6=0;
N
2
·y·N¡1
¾
Q
1
c
´
½
(x;y)jI
(l)
(x;y)6=0;
N
2
·y·N¡1
¾
Q
2
c
´
½
(x;y)jI
(l+
L
2
)
(x;y)6=0;0·y·
N
2
¡1
¾
(3.7)
Then the matrix rank r and the number of unknowns y are expressed as
r = N
2
¡(jP
1
\P
2
j+jQ
1
\Q
2
j)
y = 2N
2
¡(jP
1
j+jP
2
j+jQ
1
j+jQ
2
j) (3.8)
34
The degree of rank deficiency can be written as
r¡y = jP
1
[P
2
j+jQ
1
[Q
2
j¡N
2
= j(P
1
c
\P
2
c
)
c
j+j(Q
1
c
\Q
2
c
)
c
j¡N
2
= ¡jP
1
c
\P
2
c
j¡jQ
1
c
\Q
2
c
j (3.9)
InCartesiantwo-foldundersampling, theupperhalfandthelowerhalfofI
(l)
aresuper-
imposed with the lower half and the upper half of I
(l+
L
2
)
, respectively. Since jP
1
c
\P
2
c
j
and jQ
1
c
\Q
2
c
j represent the area of the overlap between non-zero signals of I
(l)
and
I
(l+
L
2
)
, (3.9) tells us that the degree of rank deficiency depends on the size of this over-
lapped region.
In this study these properties are assumed to carry over into the spiral undersampling
withoutarigorousproof. Sincethesidelobeofthespiralundersamplingisring-shaped, the
radial extents of the support regions inI
(l)
andI
(l+
L
2
)
are critical to the overlapping area.
Hence, a support region of circular shape will be used for specifying the zero assumption
in each I
(l)
.
The zero assumption, in its function, is similar to the temporal filtering used in UN-
FOLD. It should be noticed that the solutions are obtained from the algebraic processes,
which enables incorporation of additional algebraic operation into the proposed frame-
work. The regularization mentioned in the next subsection can be considered as a simple
form of such incorporation.
35
3.2.4 Pseudoinverse and Regularization
Withthezeroassumption, thelinearsystembecomesoverdeterminedandtheconventional
inverse does not exist. The optimal solution in minimum norm least square sense can be
obtained by solving the normal equation A
H
Ax = A
H
b.
¡
A
H
A
¢
¡1
can be computed
only if all columns of A are independent. However, A
H
A will still be ill-posed when the
size of the zero assumption is insufficient.
We need regularization which changes the original problem such that the modified
system becomes stable while still being reasonably close to the original one. Tikhonov
filtering is the most widely used regularization method due to the existence of closed form
solution. The Tikhonov regularized linear system is written as
¡
A
H
A+½I
¢
x=A
H
b (3.10)
where the regularization coefficient ½ balances the tradeoff between solution stability and
accuracy. The addition of ½I is equivalent to replacing the original cost function with
kAx¡bk
2
+½kxk
2
which adds weighting on the norm of the solution vector.
3.2.5 Iterative Solution
With a reconstructed image of N£N pixels, the dimension of the system matrix with
two-fold undersampling will be
¡
2N
2
¢
£
¡
2N
2
¢
as seen from (3.4). Moreover, the number of
operations required for direct computation of its inverse is of order O
¡
N
6
¢
. One common
way to overcome this obstacle is using an iterative method to solve the inverse problem.
In this study, the conjugate gradient (CG) method, a classical iterative algorithm, is
36
used. It is well known that CG guarantees an exact solution only if the system matrix
is symmetric and positive definite, which can be easily proved for the regularized system
matrix A
H
A+½I.
The most computationally demanding part in each k
th
iteration of CG is the matrix-
vector multiplication
¡
A
H
A+½I
¢
x
(k)
which requires the storage of
¡
2N
2
£2N
2
¢
com-
plex matrices and O
¡
N
4
¢
operations. Since the product Ax originates from the com-
bination of several 2D convolutions as described in (3.3) and (3.4), it can be rewritten
as
Ax =
2
6
6
4
A
1
x
l
+A
2
x
l+
L
2
A
2
x
l
+A
1
x
l+
L
2
3
7
7
5
=
2
6
6
4
³
S
(0)
¤
2
I
(l)
+S
(
L
2
)
¤
2
I
(l+
L
2
)
´
N£N!N
2
£1
³
S
(
L
2
)
¤
2
I
(l)
+S
(0)
¤
2
I
(l+
L
2
)
´
N£N!N
2
£1
3
7
7
5
(3.11)
Each of the four 2D convolutions involved in (3.11) can be computed by
S¤
2
I =F
¡1
2
[F
2
(S)F
2
(I)] (3.12)
where F
2
and F
¡1
2
denote the forward and the backward 2D Fourier transforms and all
superscripts are discarded for brevity.
Additional product with A
H
can be written in a similar way
37
A
H
Ax =
2
6
6
4
A
H
1
y
1
+A
H
2
y
2
A
H
2
y
1
+A
H
1
y
2
3
7
7
5
=
2
6
6
4
³
^
S
(0)
¤
2
Y
1
+
^
S
(
L
2
)
¤
2
Y
2
´
N£N!N
2
£1
³
^
S
(
L
2
)
¤
2
Y
1
+
^
S
(0)
¤
2
Y
2
´
N£N!N
2
£1
3
7
7
5
(3.13)
where y
1
= A
1
x
l
+A
2
x
l+
L
2
, y
2
= A
2
x
l
+A
1
x
l+
L
2
and Y
j
is the reshaped version of
y
j
, i.e., y
j
= (Y
j
)
N£N!N
2
£1
. Note that the product with A
H
1
or A
H
2
is equivalent to the
convolution with the flipped and conjugated PSF. That is,
^
S
(l)
(i;j) ´ S
(l)
(N¡i;N¡j)
and thereforeF
2
³
^
S
(l)
´
=F
2
¡
S
(l)
¢
. Then each of four convolutions can be computed by
^
S¤
2
Y =F
¡1
2
h
F
2
(S)F
2
(Y)
i
(3.14)
Since 2D FFT has the complexity of O
¡
N
2
logN
¢
, the computational load at each itera-
tion can be reduced by a factor of O
³
N
2
logN
´
.
One should pay attention to the boundary effects of linear convolutions denoted by¤
2
above. The PSF S
(0)
and S
(
L
2
)
have infinitely many sidelobes along the radial direction.
If only the nominal FOV (of sizeN£N) is taken and convolved with the solutionsI
(l)
and
I
(l+
L
2
)
, the signal produced by the sidelobe outside the FOV will be ignored. Assuming
that all non-zero signals exist within the FOV, including the sidelobes within 2£FOV will
be enough for handling all aliased signals introduced in the FOV. For the implementation
of the 2D convolutions, the PSF of 2£FOV (2N £ 2N) is prepared and the solution
I
(l)
is zero-padded to generate a 2£FOV image. After computing the extended 2£FOV
38
convolutions using (3.12), desired convolutions free of boundary effects can be obtained
by cropping the 2£FOV convolution results.
3.3 Numerical Phantom Simulation
3.3.1 Methods
A dynamic numerical phantom is used for the validation of the proposed method applied
to spiral two-fold undersampling. The phantom, which resembles a short axis cardiac
image, includes a chest wall, right and left ventricles, a myocardium and papillary muscles
represented by two dots (Fig. 3.4). The ventricles and the myocardium keep contracting
and expanding periodically based on the cosine function with period 1 second. The outer
and inner diameters of the myocardium change such that the area of the myocardium
should be conserved. A total of 6 spiral interleavesare used with a matrix size of128£128.
Each set of odd and even order interleaves is assumed to be sampled simultaneously so
that motion artifacts within each set is being neglected. The time interval between each
set of interleaves is set to 48 ms which produces, for a heart rate of 60 beats per minute,
roughly 21 frames per cardiac cycle from two-fold acceleration. For error quantification,
we use the maximum error and the normalized root-mean-square (RMS) error expressed
as
E
(t)
RMS
=
v
u
u
u
u
t
P
N¡1
i=0
P
N¡1
j=0
¯
¯
¯r
(t)
i;j
¡o
(t)
i;j
¯
¯
¯
2
P
N¡1
i=0
P
N¡1
j=0
¯
¯
¯o
(t)
i;j
¯
¯
¯
2
(3.15)
whereE
(t)
RMS
denotes RMS error at timet. r
(t)
i;j
ando
(t)
i;j
represent the(i;j)
th
pixel intensity
of the reconstructed image and the true image at time t respectively.
39
Figure3.4: Diastolicandsystolicphasesofthedynamicnumericalphantom. Thephantom
was designed to resemble the heart in a short axis view, with the size of the ventricle
changing sinusoidally.
The cross section view of the zero assumption is shown in Fig. 3.5. Since we know that
most dynamic part of the phantom lies inside half the FOV, the remaining outer region
is assumed to be zero, which can be applied to all pairs of I
(l)
and I
(l+
L
2
)
, except DC and
the Nyquist frequency (l = 0). The whole region of I
(0)
must be non-zero and I
(
L
2
)
is
likely to have small non-zero values in its support region. In this case the area of assumed
zeros is insufficient and the corresponding system matrix will be ill-posed. In this study,
we simply assume that I
(
L
2
)
is zero and take the DC of aliased signal R
(0)
as the solution
of I
(0)
.
For comparison purposes, reference images are generated using sliding window recon-
struction [84] and UNFOLD-like filtering which does not involve algebraic inversion. To
make the reconstruction time points of the sliding window reconstruction coincide with
those of UNFOLD-like filtering and the proposed method, the image at each time point
is generated from the raw data at the current time point and the average of the raw data
40
Figure 3.5: Cross-sectional view of the zero assumption in the numerical phantom study.
The support region is cylindrical in xyf space.
from the previous and next time points. UNFOLD-like filtering is implemented using the
same support region in xyf space shown in Fig. 3.5.
3.3.2 Results
Normalized RMS errors, maximum errors, reconstructed images, and residue images from
the sliding window reconstruction, UNFOLD-like filtering and the proposed method are
shown in Fig. 3.6 and Fig. 3.7. The proposed method produces smaller errors than the
sliding window reconstruction in terms of both RMS and maximum values. The proposed
method produces comparable but smaller errors compared to UNFOLD-like filtering.
The largest improvement is in maximum errors because of regularization effects. The
large error fluctuation of the sliding window reconstruction arises from the fact that the
ventricle and myocardium signals change sinusoidally. The motion artifacts are small
around the time points of the maximum and the minimum of the cosine, and large at
phases between them. The image of frame number 29 that corresponds to the highest
41
slope of the cosine function shows the largest motion blurring, which is verified by the
corresponding residue image. Using filtering and the proposed method, the temporal blur-
ring of the myocardium and the ventricle is reduced due to increased temporal resolution.
Filtering and the proposed reconstruction produce small periodic error fluctuations
which have local maxima at each diastole. In diastole, when the simulated heart is at
its widest, the extent of the dynamic region exceeds the specified support region. In this
simulation, the periodic error increase is the result of an incorrect zero assumption which
designates some non-zeros pixels outside the support region.
3.4 In-vivo Experiments
3.4.1 Methods
Experiments were performed on a 1.5 T GE Signa EXCITE system, with 22 mT/m
maximum gradient amplitude and 77 T/m/s maximum slew rate. A phased array coil
with four channels was used for signal reception, and the raw data from two channels were
used for reconstruction. Real-time (free breathing and ungated) spiral balanced SSFP
cardiac imaging was conducted [64]. The imaging protocol used 36 spiral interleaves,
an in-plane resolution of 2.4 mm£2.4 mm, a slice thickness of 10 mm, a 60 degree flip
angle, and a repetition time (TR) of 6.4 ms. The gradient waveforms were designed such
that both the zeroth and first moments are zero between RF pulses to ensure steady-
state signal coherence and to avoid in-plane flow effects. A total of 100 image frames were
reconstructed from the proposed method, which used a set of odd or even order interleaves
alternately to form one image frame.
42
W XW YW ZW [W \W ]W W WUWY WUW[ WUW] WUW_ WUX WUXY frames
RMS error
sliding window
filtering
proposed method
W XW YW ZW [W \W ]W W WUY WU[ WU] WU_ X frames
max error
sliding window
filtering
proposed method
Figure 3.6: Normalized RMS error and maximum pixel-by-pixel error from simulated
phantom experiments. The proposed method has smaller error than sliding window re-
construction. Compared to filtering, slight improvements are observed because of regu-
larization effects. The large periodic increase in error in sliding window reconstruction
occurs when motion is the fastest (twice per cycle). The periodic increase in errors for
filtering and the proposed method is related to the size of the dynamic region approaching
the edge of the support region at diastolic phases (once per cycle).
43
sliding window filtering proposed method
Y^ Y_ Y` Figure 3.7: Representative reconstructed images and residue images from three consec-
utive frames around early systolic phase. Residue images are the difference between the
reconstructed images and the known true images, and are displayed using harsh window-
ing.
44
Figure 3.8: Cross-sectional view of the zero assumption in the in-vivo experiments. The
support region resembles a truncated cone in xyf space.
Since the numerical rank of the system matrix depends on the radial extents of the
support region in I
(l)
, as mentioned in the section 3.2.3, it is reasonable to use a stack of
circles for the zero assumption inxyf space. Then, the design parameter to be determined
is the profile of the diameters along the frequency axis. In this study, a ramp-shaped
profile was used, assuming the temporal variation of images to be dominant in the central
region (Fig. 3.8). Note that high frequency region, which amounts to 16 % of the whole
bandwidth, is filtered out since in this band there will be strong aliased signals arising
from the sidelobe. The linear systems associated with the zero assumption displayed in
Fig. 3.8 will be slightly ill-posed. We observe that this loose zero assumption produces
more robust results than a tight support region which enables better conditioning, but
fails to truly contain all non-zero signals. Reference images were generated using the
sliding window reconstruction and UNFOLD-like filtering which used the same support
region shown in Fig. 3.8.
45
3.4.2 Results
Fig.3.9(a)showssystolicshortaxisimages(20
th
frame)reconstructedfromslidingwindow
reconstruction, UNFOLD-like filtering, and the proposed method. The image difference
is subtle, but slight improvement can be found from the examination of fine cardiac
structures. Fig. 3.9(b) shows the signal intensity profiles through the moderator band in
the right ventricle. The paths of the signal profiles are the missing parts of the horizontal
and vertical lines shown in Fig. 3.9(a). The band is most sharply captured using the
proposed method, followed by filtering, and finally sliding window, which matches the
expected order of least motion-blur to most motion-blur. Another set of systolic short
axis images (59
th
frame) are shown in Fig. 3.10(a) using harsh windowing to visualize
motion artifacts which appear as background whirling. The average signal intensity in
the region of interest (ROI), specified by the dotted box, is plotted as a function of time
(Fig.3.10(b)). Themeanvaluesfluctuateattheheartrate, whichimpliesthattheartifacts
in the ROI relate to cardiac motion. The proposed method reduces motion artifacts as
wellasUNFOLD-likefilteringcomparedtoslidingwindowreconstruction. Reconstruction
was performed in MATLAB
R °
on a Dell 8400 desktop computer (3.6 GHz Intel processor,
1 GB RAM). Proposed reconstruction of 100 temporal frames, using two coils, required
approximately 5 minutes of CPU time.
3.5 Discussion
The most critical component in the implementation of the proposed method is determin-
ing the shape of the zero assumption in xyf space. If the zero assumption is insufficient
46
horizontal
intensity profiles
vertical
sliding window proposed method filtering
sliding window
filtering
proposed method
a
b
Figure3.9: Evaluationofthedepictionoffinecardiacstructureinin-vivoimagesfromslid-
ing window reconstruction, UNFOLD-like filtering and the proposed method. (a): Signal
profiles across the moderator band of the right ventricle (black cross-hairs) were measured
(20
th
time frame). (b): Intensity profiles indicate sharp depiction of the moderator band
in the order of the proposed method, filtering and the sliding window.
(i.e., the diameters of the circular windows in Fig. 3.8 are large), the numerical rank of the
corresponding linear systems will be smaller. Excessive zero assumption forces incorrect
zeros into certain regions of the solution. For example, since in-vivo data are always noisy,
applying a zero assumption even in the stationary region of the FOV will cause noise am-
plification in support regions [77]. In reality, it is difficult to achieve sufficient conditioning
and a completely accurate zero-assumption. From our experience, using an insufficient
47
0 10 20 30 40 50 60 70 80 90 100
mean of ROI
frames
a
b
sliding window proposed method filtering
sliding window
filtering
proposed method
Figure 3.10: Evaluation of motion artifacts in the three sets of in-vivo images. (a):
To visualize motion artifact reduction, short axis images from the three methods were
windowed harshly (59
th
time frame), and signal in a region outside of the body (white
box) was measured. (b): Mean signal within the box is plotted as a function of time
illustrating that the proposed method reduces motion artifacts as well as filtering.
48
zero assumption with insufficient conditioning produces more robust results than using an
excessive assumption with sufficient conditioning. During implementation, we manually
specified the shape of the zero assumption by referring to the signal distribution in xyf
space obtained from the sliding window reconstruction. Support regions were made larger
than the distribution from the sliding window reconstruction which underestimates high
temporal frequency components. Automatically determining proper shapes of support
regions is still an open question, and will be application dependent.
The computational complexity of the CG method employed in the proposed method
is of practical importance. In the in-vivo studies, we observe that CG converges at around
8 iterations for one set of linear equations, which should be repeated L=2 times for whole
sets of systems. One design parameter affecting the convergence rate is the extent of the
support region inI
(l)
. If it becomes large, the corresponding system will be more ill-posed,
which will increase the required number of iterations. The regularization coefficient ½ is
another parameter that influences the convergence rate. A large ½ increases the values
of the diagonal elements of the regularized system matrix (refer to (3.10)) decreasing the
matrix condition number, and vice versa. In the in-vivo studies we used the coefficient
value ½ = 0.01, which is 0.46% of the temporal DC values averaged over the FOV. We
observe that when ½ = 0.001 is used instead, the required number of iterations becomes
around 16 with little change in image quality. When a larger value ½ = 0.1 is used, the
iteration number is reduced to 4, but the restored images are smoothed more.
One assumption used in the analysis of the aliasing from spiral undersampling is that a
group of interleaves are collected at the same time, which is not true. The analytic model
used in this study will be less accurate when larger numbers of interleaves are combined
49
to produce a single time frame. In fact, the true distribution of aliased signals in xyf
space can be obtained by considering a single spiral acquisition as a single time frame.
Combining a set of spiral interleaves is equivalent to applying sinc weighting followed by
aliasing along the frequency dimension. The dilation of the sinc function is determined
according to the number of interleaves combined into a single time point. Of course, the
system equations derived without combining interleaves will be the most accurate, but
in this case the highest temporal frequency to be resolved becomes larger without the
increase in system rank. The tradeoff between the accuracy of system equations and the
increased dimensionality of the solution needs to be explored in future works.
The two-fold acceleration covered in this study might be a trivial example considering
the involved numerical complexity of the proposed approach. The same problem can
be solved via UNFOLD-like filtering which is free from an algebraic inversion, as shown
in the Results section. The numerical simulation showed that filtering produces subtle
increases in RMS error, and more apparent increases in maximum error, which can be
explained by the regularization effects of the proposed method. The in-vivo study showed
that the reduction of motion-blur and background whirling with the proposed method is
comparable to filtering. Still, two-fold acceleration is not the ultimate goal, and should
be considered as a simple example to test the proposed algebraic model. One advantage
of the proposed approach is its capability of being fused with any acceleration method
based on algebraic operation, for example non-Cartesian sensitivity encoding [79]. Any
arbitrary sampling scheme (arbitrary k-space trajectory and arbitrary view ordering) can
be handled within the proposed algebraic framework.
50
3.6 Conclusions
We have developed a reconstruction method for accelerating dynamic spiral MRI, which
is based on the algebraic formulation of aliasing in xyf space. Aliasing from m-fold
undersamplingisformedsuchthatonlymtemporalfrequenciesarecoupledtothesameset
of ill-posed linear equations. Each ill-posed system can be solved independently assuming
afinitesupportregioninthesolution. Numericalsimulationandin-vivoexperimentsusing
spiral two-fold undersampling demonstrate that the proposed method reduces motion blur
and motion artifacts due to increased temporal resolution.
Thereductionofmotionartifactsandmotionblurachievedusingtheproposedmethod
is comparable to but slightly larger than UNFOLD-like filtering which is computation-
ally more efficient. However, the proposed algebraic framework enables the potential to
incorporate additional acceleration techniques based on linear operations, such as non-
Cartesian sensitivity encoding. The proposed method also benefits from the flexibility
to handle arbitrary k-space trajectories, arbitrary acquisition orders, and arbitrary data
dimensionality.
51
Chapter 4 :
Three-Dimensional Myocardial Perfusion Imaging:
Feasibility
4.1 Introduction
The goal of myocardial perfusion imaging (MPI) is to evaluate regional blood flow to the
heart, which is essential for early detection of ischemic heart disease (IHD) (see Fig. 1.1).
While catheter based coronary angiography can detect coronary stenosis that is the ba-
sis for IHD, MPI directly assesses physiological significance of the reduced blood supply.
Nuclear medicine imaging approaches such as single photon emission computed tomogra-
phy (SPECT) and positron emission tomography (PET) have long been used in clinical
studies, but are limited in spatial resolution, and require the exposure to ionizing radia-
tion [36,71,89].
First-passmagneticresonance(MR)myocardialperfusionimaging(MPI)ispotentially
an ideal tool for the evaluation of cardiac perfusion, due to its high spatial resolution
52
and lack of the exposure to ionizing radiation. MR MPI uses contrast agent such as
Gadolinium-DTPA which shortens the T
1
of tissue by the following relation.
1
T
1
=
1
T
1;0
+¹[C] (4.1)
where T
1
, T
1;0
, ¹ and [C] are the resulting T
1
, original T
1
in the absence of contrast
agent, relaxivity and the concentration of contrast agent, respectively. Time resolved
and T1-weighted pulse sequences are used to track the passage of contrast agent. Two
dimensional (2D) multi-slice imaging is a typically used data acquisition mode in MR
MPI, and its diagnostic value has been demonstrated with high sensitivity (>85%) and
specificity (>75%) [4,48,75,104].
Besidesitsdiagnosticability, theprognosticutilityofMPIisimportantforappropriate
treatment of symptomatic patients. Patient follow-up studies based on nuclear medicine
and echocardiography have shown that the extent of myocardial perfusion defects is vital
to the prediction of cardiac outcomes [9,26,59]. However, its accurate estimation is
restricted in conventional 2D multi-slice MPI, due to the limited tradeoff between spatial
coverage and inter-slice spacing.
Three dimensional (3D) MPI is an attractive alternative since it supports contigu-
ous spatial coverage and therefore has the potential to more accurately estimate defect
size (volume). Moreover, 3D encoding inherently retains high signal-to-noise ratio (SNR)
efficiency and high capacity for parallel imaging acceleration [102], and has been demon-
strated at 1.5 T [42]. The purpose of this study was (i) to evaluate the performance
of 3D MPI and 2D multi-slice MPI for estimating defect size by a phantom study, and
53
(ii) to demonstrate the feasibility of in-vivo 3D MPI accelerated by sensitivity encoding
(SENSE) with whole left ventricular (LV) coverage on a 3 Tesla platform.
4.2 Pulse Sequence
The pulse sequence used in this study is a 3DFT gradient echo (GRE) acquisition that
is preceded by a 90
±
global saturation pulse and a saturation delay time T
SR
(Fig. 4.1).
We specifically adopted an adiabatic saturation pulse (BIR-4) due to its insensitivity to
transmit B1 non-uniformity [45,93]. For slab selective excitation, a sinc pulse with a
time-bandwidth-product (TBW) of four was used. ECG gating was used to synchronize
data acquisition at diastole to minimize the effects of cardiac motion. Figure 4.2 shows
actual waveforms of BIR-4 saturation RF pulse, and 3DFT GRE readout. BIR-4 pulse
enables robust saturation insensitive to B1 variation. In 3DFT readout, Gz rewinder was
combined with spoiler to reduce TR (2.3 ms).
Three-dimensional k-space for a given spatial resolution and FOV was undersampled
at a net acceleration factor of six. Since the readout encoding line (k
x
) is always fully
sampled, undersampling in the other two encoding lines (k
y
, k
z
) can be viewed in the
2D k
y
¡ k
z
plane (Fig. 4.1 ). Six-fold undersamped data were obtained by three-fold
undersampling in k
y
and two-fold undersampling in k
z
. The position of sampled k-space
encoding lines was shifted in a cyclic pattern in the k
y
¡k
z
plane such that acquired data
from six neighboring cardiac cycles could be combined to form a coil sensitivity map for
subsequent SENSE reconstruction.
54
B
I
R
4
Six-fold
undersampled
3DFT readout
T
SR
Figure 4.1: Pulse sequence diagram and sampling pattern for 3D myocardial perfusion imaging. The pulse sequence consists of
a 90
±
BIR-4 adiabatic RF pulse for myocardial saturation, followed by a saturation recovery time and undersampled 3DFT GRE
readout. Three-dimensional k-space was undersampled by factor of three and two in k
y
and in k
z
, respectively, resulting in a net
undersampling factor of six. The k-space sampling locations were shifted in a cyclic pattern such that coil sensitivity map could
be obtained from the data in six neighboring R-R intervals. SENSE reconstruction of undersampled data from the k
th
cardiac
cycle utilized coil sensitivity maps derived from the (k-2)
th
to (k+3)
th
cardiac cycles.
55
RF
Gz
Gx
Gy
imaging sequence (3DFT)
|RF|
Gz
saturation sequence (BIR4)
RF
Figure 4.2: Waveforms of saturation pulse and 3DFT imaging. BIR-4 RF pulse (8 ms)
was used for B1 insensitive 90
±
saturation. In 3DFT imaging, Gz rewinder was combined
with spoiler to reduce TR (2.3 ms).
Saturation recovery time (T
SR
) and flip angle (®) are critical parameters which deter-
mines k-space response and the contrast between different T
1
values. Underperfusion can
be differentiated from normal perfusion by the T
1
contrast, and the need for large contrast
is obvious. Due to the gated data acquisition at diastole only, the magnetization remains
in a transient state during the acquisition. Variation in the magnitude of transverse mag-
netization over different k-space encoding lines can cause artifacts in object domain by
distorting PSF. The two imaging parameters are determined by the following two steps.
^ ®(T
SR
)=argmin
®
E
T
1
[V
k
[M
xy
(®;T
SR
;T
1
;k)]] (4.2)
^
T
SR
=argmax
T
SR
M
xy
(^ ®(T
SR
);T
SR
;T
1s
;N
t
=2)¡M
xy
(^ ®(T
SR
);T
SR
;T
1l
;N
t
=2) (4.3)
where M
xy
represents transverse component of magnetization as a function of flip angle,
T
SR
, T
1
, and an index to phase encoding linek, E
T
1
(¢) and V
k
(¢) denote expectation over
56
20 40 60 80 100 120
0
0.05
0.1
0.15
0.2
0.25
0.3
view index
M
xy
/M
0

T
1
=200 ms
T
1
=400 ms
T
1
=600 ms
T
1
=800 ms
T
1
=1000 ms
Figure 4.3: Transverse component of magnetization over k-space encoding lines with
T
SR
=130 ms and ®=10
±
. The k-space response shows little variation over encoding lines
in the wide range of T
1
values.
T
1
and variance over encoding line index k, respectively. T
1s
and T
1l
are representative
small and large T
1
values, set to 300 ms and 1100 ms. N
t
is the total number of phase
encoding lines, and thus N
t
=2 is the index to the center of k-space assuming a linear view
order. Fig. 4.3 shows normalized transverse magnetization over k-space encoding lines
using T
SR
=130 ms and ®=10
±
selected from (4.2) and (4.3). It can be seen that k-space
response is reasonably flat in a range of T
1
=[ 200, 1000] ms.
All subsequent experiments were performed on a General Electric 3 Tesla scanner with
40 mT/m gradient amplitude and 150 T/m/s gradient slew rate, using a body-coil for RF
transmission and an eight-channel cardiac coil array for signal reception.
57
Figure 4.4: Cardiac phantom and rubber insets. Left ventricular and myocardial cavities
were filled with Gadolinium-doped saline of 2 mmol/L and 0.6 mmol/L. The rubber insets
ranging in size from 5-70 % of the myocardial volume, mimic transmural perfusion defects.
4.3 Phantom experiments
4.3.1 Phantom Description
Aheartphantom(ModelRH-2, CapintecInc.) wasscannedtocomparetheaccuracyof3D
and 2D multi-slice MPI in estimating defect size (Fig. 4.4) [49]. The phantom contained
left ventricular cavity of 132 mL volume and myocardial cavity of 193 mL volume, 10.5
mm thickness and 11 cm long-axis length. The two cavities were separated by a» 2mm
thick wall, and were filled with Gadolinium-doped saline of 2 mmol/L and 0.6 mmol/L, to
mimic peak enhancement in LV and myocardium, respectively [44]. Seven rubber insets
ranging in size from approximately 5-60 % of the myocardial volume were used to mimic
transmural perfusion defects.
4.3.2 Scan Protocol
3D scanning was performed using the pulse sequence illustrated in Fig. 4.1. Full k-space
data was acquired in six segments with NEX=10, requiring 6£10= 60 data acquisitions.
58
2D multi-slice scanning was performed using a BIR-4 saturation pulse followed by 2DFT
GRE acquisition. Full k-space data was acquired in two segments per slice with NEX=10,
requiring a total of 2£10 = 20 data acquisitions per slice.
Specific imaging parameters used for 3D and 2D 4-slice/3-slice scans are represented
in Table 4.1. The long-axis spatial coverage was 10.0 cm, 8.2 cm and 7.0 cm in 3D, 2D
4-slice, and 2D 3-slice scans, respectively. The gap between neighboring slices was 14 mm
and 20 mm in 2D 4-slice and 2D 3-slice scans, respectively. Three sets of cross-sectional
scan planes each of which was shifted by 10 mm along the partition encoding direction,
were used for 2D multi-slice imaging to examine the sensitivity of scan plane locations on
the estimation of defect size.
4.3.3 Image Analysis
The borders of defects in reconstructed images were manually outlined, and the volume of
a defect was expressed as the percent of the whole myocardium. To accurately estimate
partial inclusion of the insets in a voxel, the phantom image with a defect was normalized
by the image without a defect, and the fraction of defect was linearly interpolated in each
normalized voxel.
4.3.4 Results
Figure 4.5 shows representative images of the phantom with a 20.3% rubber inset obtained
from 3D, 2D 4-slice and 2D 3-slice scans. Figure 4.6 contains scatter plots of errors in
the estimation of defect size using the three methods. The horizontal axis and vertical
axis represent the true volume fraction of the defect, and estimation errors, respectively.
59
The mean and §1.96£standard deviation are denoted by solid and dotted lines, respec-
tively. In 2D cases, three measurements were made for each phantom defect from three
shifted scan orientations. 3D imaging resulted in a significantly smaller bias/ standard
deviation (-0.44/1.49 %) compared to 2D 4-slice imaging (2.23/2.97 %) and 2D 3-slice
imaging (2.59/3.18 %). 2D multi-slice methods show a slight bias towards overestimation,
presumably due to the fact that the majority of defects in this study were positioned at
mid-short axis level and were partially visible in all slices. 2D multi-slice methods also
show a larger measurement standard deviation due to incomplete spatial coverage.
4.4 In-vivo Experiments
4.4.1 Scan Protocol
In-vivo experiments using 3D MPI were performed in three volunteers. Written informed
consent was obtained from all participants. Imaging parameters were the same as those
in Table 4.1 except T
SR
=130 ms and matrix size= 100£66-69£10. Data acquisition was
located at the center of diastole, and the acquisition time was 304 ms. A proton density
Table 4.1: Imaging parameters of 3D and 2D multi-slice scans.
3D 2D (4-slice/3-slice)
T
SR
100 ms 100 ms
TR 2.3 ms 2.9 ms
TE 0.9 ms 1.3 ms
flip angle 12
±
12
±
matrix size 100£ 66£ 10 120£ 80
FOV 30£ 30£ 10 cm
3
30£ 30 cm
2
slice thickness 10 mm 10 mm
slice gap 0 14/ 20 mm
long-axis coverage 100 mm 82/ 70 mm
NEX 10 10
60
a
b c
+10mm shift in z
-10mm shift in z
reference scan plane
Figure 4.5: Representative heart phantom images with a 20.3 % rubber inset. (a): images from 3D scan, (b): images from 2D
4-slice scan, and (c): images from 2D 3-slice scan. In 2D multi-slice case, images from the three sets of cross-sectional scan planes
were shown. In all phantom images, the dark rim around the LV was created by the »2 mm thick separation wall between the
LV and the myocardial cavities.
61
0 10 20 30 40 50 60
0
2
4
6
8
10
True fractional volume (%)
Error (%)
0 10 20 30 40 50 60
0
2
4
6
8
10
True fractional volume (%)
Error (%)
0 10 20 30 40 50 60
0
2
4
6
8
10
True fractional volume (%)
Error (%)
3D 2D 4-slice 2D 3-slice
Figure 4.6: Distribution of errors in the estimation of phantom defect size. The horizontal axis and vertical axis represent the
true volume fraction of the defect, and the estimation errors, respectively. In each plot, the solid and dotted lines represent
the mean and§1.96£standard deviation, respectively. 3D imaging resulted in a significanty lower bias bias/ standard deviation
(-0.44/1.49%) compared to 2D 4-slice (2.23/2.97%) and 2D 3-slice (2.59/3.18%).
62
weighted data set was obtained using 4
±
flip angle with the saturation pulse turned off
during the first six cardiac cycles.
Contrast media (0.1 mmol/kg Gd-DTPA, Magnevist) was injected at a rate of 5 ml/s
followed by 20 ml saline flush at the same rate. Subjects were instructed to hold their
breath as long as possible.
4.4.2 Reconstruction
All image reconstruction was performed off-line using MATLAB
R °
(Mathworks, Natick,
MA). Proton density weighted images were obtained by simply combining the k-space
data from the first six cardiac cycles. 2D SENSE was used for the image reconstruction
from undersampled 3DFT perfusion data [102], [41]. The data from six neighboring heart
beats are combined to generate an aliasing free image and to form coil sensitivity maps.
In contrast to 2DFT undersampling, aliasing in 3DFT undersampling occurs in both
phase encoding and partition encoding directions. Considering three-fold and two-fold
undersamplings in phase encoding and partition encoding directions, aliased points can
be resolved by the following equation.
Sx=r (4.4)
63
S =
2
6
6
6
4
S
(1)
(y;z) S
(1)
(y;z+
M
2
) S
(1)
(y+
N
3
;z) S
(1)
(y+
N
3
;z+
M
2
) S
(1)
(y+
2N
3
;z) S
(1)
(y+
2N
3
;z+
M
2
)
S
(2)
(y;z) S
(2)
(y;z+
M
2
) S
(2)
(y+
N
3
;z) S
(2)
(y+
N
3
;z+
M
2
) S
(2)
(y+
2N
3
;z) S
(2)
(y+
2N
3
;z+M=2)
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
S
(C)
(y;z) S
(C)
(y;z+
M
2
) S
(C)
(y+
N
3
;z) S
(C)
(y+
N
3
;z+
M
2
) S
(C)
(y+
2N
3
;z) S
(C)
(y+
2N
3
;z+
M
2
)
3
7
7
7
5
x
T
=
·
I(y;z) I(y;z+
M
2
) I(y+
N
3
;z) I(y+
N
3
;z+
M
2
) I(y+
2N
3
;z) I(y+
2N
3
;z+
M
2
)
¸
r
T
=
·
I
(1)
alias
(y;z) I
(2)
alias
(y;z) ::: I
(C)
alias
(y;z)
¸
(4.5)
where S
(k)
(y;z) represents coil sensitivity of k
th
coil element, and I(y;z) and I
(k)
(y;z)
denote true image and aliased image observed byk
th
coil element. C,N andM are a total
number of coil elements, the dimension of phase encoding and the dimension of partition
encoding respectively. The index of readout dimension (x) is skipped in all arguments for
brevity. Note that the number of coil elements greater than six will render the system
matrix over-determined.
4.4.3 Image Analysis
Raw perfusion images were normalized by proton density images to remove variations in
the receiver coil sensitivities. Corrected images were then segmented into six (basal and
mid-short axis levels) or four (apical level) myocardial sectors, according to the American
Heart Association 17-segment model [13]. Within each myocardial sector, a time intensity
curve (TIC) was generated and the corresponding upslope value was computed by linear
fitting of the data during signal enhancement.
Noise standard deviation is typically calculated from background region-of-interest
(ROI)forthecomputationofSNR,assumingspatiallyuniformnoisedistribution[31],[15],
but this assumption is no longer valid in parallel imaging reconstruction due to spatially
64
varying reconstruction noise. Hence in this work, SNR was calculated by the “difference
method” [83], where noise standard deviation was estimated from the difference in two
consecutive images as follow.
SNRj
ROI
=
¹ x
ROI
¾
ROI
=
mean(I
1
+I
2
)j
ROI
=2
jstd(I
1
¡I
2
)j
ROI
=
p
2
(4.6)
4.4.4 Results
Fig. 4.7 and Fig. 4.8 show representative perfusion images from one subject, at every other
time frame. Overall image quality is excellent due to high SNR and effective SENSE
reconstruction, clearly showing the arrival and the passage of contrast agent. A small
amount of signal fluctuation is seen during the course of RV and LV enhancement. This
flickering artifact is due to errors in sensitivity map estimation where combining data
from neighboring six cardiac cycles may be inaccurate during rapid signal changes. Dark
rim artifacts around LV are seen on most of the volunteer images, presumably due to
low spatial resolution [20]. The extent of these rims may be reduced by improving the
spatial resolution, which requires acceleration of higher rate. Potential techniques for
higher acceleration are described in the Discussion section.
Figure 4.9 shows SNR at peak myocardial enhancement and CNR values in the first
eight slices. A total number of 12-18 segments (number of volunteers (3)£number of
segments per slice (4 to 6)) were used to compute the average and standard deviation of
SNR and CNR values in each slice. Collectively, the SNR and CNR were 57.6§22.0 and
37.5§19.7 across all segments. Note that the two values are higher in slices at center than
slices at edges due to the slab excitation profile of parabolic shape.
65
WUWG XU_[G ZU]_G \U\YG ^UZ[G `UYWG XXUX[G Figure 4.7: Representative 3DFT perfusion images selected every other heart beat (1
th
-13
th
time frames). Blood pools and
myocardium are well suppressed by BIR-4 saturation pulse at the beginning. LV and RV get bright sequentially by the arrival of
contrast agent.
66
XZUW^G X[U`\G X]U__G X_U_]G YWU_`G YYU`_G Y\UX[G Figure 4.8: Representative 3DFT perfusion images selected every other heart beat (15
th
-27
th
time frames). After LV enhancement,
myocardium reaches peak enhancement.
67
Peak SNR CNR
0
20
40
60
80
100
0
20
40
60
80
1 2 3 4 5 6 7 8 All
Slice index
Figure 4.9: SNR at peak myocardial enhancement and CNR values for the first eight
slices. The SNR and CNR were 57.6§22.0 and 37.5§19.7 across all segments. Note that
the two values are higher in slices at center than slices at edges due to the slab excitation
profile of parabolic shape.
68
Figure4.10showsrepresentativesegment-basedTICsobtainedfromthefirstnineslices
in a volunteer. All TICs show homogeneous myocardial enhancement, consistent with
normal perfusion. Flickering artifacts caused by the errors in sensitivity map estimation,
manifest as fluctuations in the TIC, as denoted by an arrow. The artifacts are most
severe at the start of LV enhancement, but have no significant effect on the TIC upslope.
Figure 4.11 depicts the regional TIC upslope from the first nine slices in this study. The
average upslope value of the whole myocardium is scaled to 100 % in the map. Only small
variation of upslope values was seen over the whole myocardium.
4.5 Discussion
The heart phantom study has demonstrated the superiority of 3D MPI over 2D multi-
slice MPI in sizing perfusion defects. The limitation of this phantom study was that only
transmural defects were tested, due to the lack of a rubber model for a subendocardial
defect. The capability of sizing various types of defects will be investigated by perfusion
scans in patients with known CAD.
Whenviewedinavideoformat, thein-vivoperfusionimagesshowedflickeringartifacts
due to errors in sensitivity map used for SENSE reconstruction. This artifact appeared
only during the start of RV and LV enhancement, which barely affected the upslope of
myocardial TIC. The flickering can be reduced by employing variable density k-space
undersampling. The artifacts will be reduced and diffused over the FOV by taking central
region of k-space only from current cardiac cycle, and combining outer part of k-space
from neighboring cardiac cycles. More actively, strong quadratic regularization can be
69
W XW YW ZW W WU X WU Y WU Z WU [ GOP G GOUUP W XW YW ZW W WU X WU Y WU Z WU [ GOP W XW YW ZW W WU X WU Y WU Z WU [ GOP W XW YW ZW W WU X WU Y WU Z WU [ GOP W XW YW ZW W WU X WU Y WU Z WU [ GOP s} h u{ h u{ Gz l w z l w p u m w vz sh { W XW YW ZW W WU X WU Y WU Z WU [ GOP G GOUUP W XW YW ZW W WU X WU Y WU Z WU [ GOP W XW YW ZW W WU X WU Y WU Z WU [ GOP W XW YW ZW W WU X WU Y WU Z WU [ GOP s} h u { z lw p u m sh { ]G [G Figure 4.10: Time intensity curves in a healthy volunteer. Representative time intensity curves (TIC) are shown from the first nine
slices on a segment basis. Six segments were used for the first six slices (basal to mid-short axis levels), and four segments were
used for the rest (apical level). All TICs show homogenous myocardial enhancement, consistent with normal perfusion. Flickering
artifacts caused by errors in coil sensitivity map estimation were observed (see an arrow in the 8
th
and 9
th
slices). The artifacts
are most severe at the start of LV enhancement, but have no significant effect on the TIC upslope.
70
W YW [W ]W _W XWW XYW Figure4.11: Colormapvisualizationoftime-intensitycurveupslopeinahealthyvolunteer.
Upslope value of time intensity curve was computed by linear fitting of the curve during
signal enhancement, and visualized by a color map. The average upslope value of the
whole myocardium is scaled to 100 %.
used for the separate reconstruction of each coil map from smaller number of neighboring
k-space data sets, which will involve less signal change in blood pools.
Six-fold acceleration was achieved by applying SENSE in the in-vivo studies, but
further acceleration will be explored for the following reasons. First, spatial resolution of
3.0£4.3¡4.5£10 mm
3
used in this study may be insufficient for visualizing subendocardial
defects. Higher spatial resolution may also help to avoid the dark rim artifacts which
were seen on most of volunteer images. Gibbs ringing is one of critical reasons for the
artifacts, and the use of high spatial resolution has been shown to restrict its transmural
extent [20,76]. Second, the data acquisition time (304 ms) need to be shortened to avoid
potential motion artifacts in subjects with high heart rates. An eight-channel receiver
coil was used in this study, and using larger array of coil elements would be one way to
increase acceleration rate. Recently introduced compressed sensing concept is another
promising source for further acceleration [56]. The addition of a regularization of l1-norm
71
such as total variation to SENSE reconstruction will enable image reconstruction from
more highly undersampled k-space data.
SNR in a 3D encoded image varies along the partition encoding direction in proportion
to the shape of excitation profile (refer to Fig. 6). With TBW of slab excitation= 4 used
in this study, the lowest peak SNR in the basal slice (slice #1) was 64 % of the highest
SNR in the mid short axis slice (slice #4). Evenly high SNR over partition slices should
be desired for consistent diagnostic capacity. Variable-rate selective excitation (VERSE)
is a technique that modifies the original RF and gradient waveform such that either peak
RF value or total RF duration can be reduced with the slice profile unchanged [14,29].
We will explore the use of VERSE technique to increase TBW of slab excitation RF for
a sharper profile without increase in RF duration.
In addition to high SNR and contiguous coverage, the large convenience of image
registration is a potential merit of 3D encoding. The registration of time-resolved perfu-
sion images is important for subsequent semi-quantitative and/or absolute quantitative
perfusion analysis such as TIC upslope. Typically, the registration of 2D images should
be performed by nonrigid methods that allow image deformation to accurately model
through-plane motion. Since through-plane motion can be registered as accurately as in-
plane motion in 3D images, respiratory motion can be corrected for by relatively simple
3D rigid-body model.
72
4.6 Conclusion
We have demonstrated that 3D MR MPI is superior to 2D multi-slice MR MPI in sizing
transmural perfusion defects. Mean/standard deviation of errors in estimating the volume
fraction across the seven defects were -0.44/1.49 %, 2.23/2.97 %, and 2.59/3.18 % in 3D,
2D 4-slice, and 2D 3-slice imaging, respectively. We have also demonstrated the feasibility
of in-vivo 3D MPI accelerated by sensitivity encoding at 3T. 3D MR MPI at 3T can
provide complete LV coverage with excellent image quality and high peak SNR and CNR
(57.6§22.0 and 37.5§19.7, respectively). Evaluation of this technique in patients with
known CAD, and the use of regularized parallel imaging reconstruction for higher spatial
resolution and speed, are under investigation.
73
Chapter 5 :
Three-Dimensional Myocardial Perfusion Imaging: Systole
vs Diastole
5.1 Introduction
The previous chapters demonstrated that 3D MR MPI is a promising alternative to 2D
multi-slice approach, due to its contiguous spatial coverage. Three-dimensional MPI is
more accurate than 2D multi-slice technique in estimating the size of defects, which is
crucial for prognosis in symptomatic patients, as validated by patient follow-up studies
[9,26,59,108]. Additional advantages of 3D MPI include high signal-to-noise ratio (SNR),
imaging of all myocardium at a single cardiac phase, a potentially long saturation recovery
time, and the capacity for high rate parallel imaging. Typically, the acquisition window is
positioned at mid-diastole where the heart is the most stationary within an R-R interval.
In healthy subjects, previous 3D MPI studies have used »300 ms acquisition, although
the optimal duration and timing of the acquisition tends to vary over subjects [38].
End-systole is another cardiac phase well-suited for data acquisition in MR MPI. It is a
shorter quiescent period compared to mid-diastole, but is less sensitive to R-R variability
74
and arrhythmia. In normal subjects, the duration and timing of diastasis is known to
vary more significantly than the duration and timing of end-systole, as a function of heart
rate [103]. For example, an increase in heart rate from 80 to 90 beats per minute (BPM)
results 5.6% and 16.4% decreases in the duration of systole and diastole, respectively.
Systolic data acquisition is likely to be more robust in patients with atrial fibrillation
(AF),themostcommoncardiacarrhythmiawithaprevalenceof9%inelderlypeople(>65
years) with cardiovascular disease [23]. In patients with chronic AF, mean variations in
systolic and diastolic time intervals were reported to be 10.6% and 37.9%, respectively [3].
A change in the duration and timing of a stable period induces either motion artifacts
when data is acquired as scheduled, or loss of temporal resolution when data is rejected.
It should be noted that the myocardium is thicker at the systole. Thus, lower in-plane
spatial resolution may be tolerated, which, to some extent, compensates for the reduction
in available acquisition time compared to diastole.
In this chapter, we investigate the feasibility of 3D MPI with data acquisition at end-
systole, and compare it to 3D MPI with data acquisition at mid-diastole. We implemented
a MPI pulse sequence that acquires 3D perfusion data set at both end-systole and mid-
diastole, and performed a comparison of perfusion images and region-based time-intensity
curves from the two cardiac phases. In order to achieve complete LV coverage, we used a
relatively low spatial resolution, which is more likely to suffer from dark rim artifacts and
maybeincapableofidentifyingsmallperfusiondefects. Inthispilotstudy, westudiedonly
healthy subjects with no known perfusion defect. We hypothesize that systolic 3D MPI
anddiastolic3DMPIwillprovidecomparablemyocardialresolutionandsemi-quantitative
75
perfusion indices, and thus systolic imaging may be favorable in subjects with severe R-R
variability.
5.2 Methods
All subsequent experiments were performed on a General Electric 3 Tesla scanner with 40
mT/m gradient amplitude and 150 T/m/s gradient slew rate, using a body-coil for RF
transmission and an eight-channel cardiac coil array for signal reception.
5.2.1 Pulse Sequence
The pulse sequence consisted of two saturation recovery 3DFT acquisition modules, one
at end-systole and the other at mid-diastole per each R-R interval (Fig. 5.1). For the
preparatory saturation, we used tailored hard-pulse train that is robust in the presence
of B
0
and B
1
inhomogeneities and involves relatively low RF power [94]. A sinc-shaped
RF pulse with a time-bandwidth-product of four was used for slab selective excitation.
Acquisitions were triggered using vector-ECG. A cine balanced SSFP scan was used to
identify the start and end of end-systole and the start and end of mid-diastole relative
to the R-R interval. Trigger delays were prospectively adjusted in real-time based on
the heart rate. Three-dimensional k-space for given FOV and spatial resolution was
undersampled by factor of three and two along k
y
and k
z
directions, e.g., by a net factor
of six, in both systolic and diastolic imaging. The position of sampled k-space encoding
lines was shifted in a cyclic pattern such that the data acquired from six neighboring
cardiac cycles can be combined to generate a coil sensitivity map for sensitivity encoding
(SENSE) reconstruction.
76
S
A
T
3DFT GRE

254~305 ms
T
SR,dias
T
SR,sys

3DFT GRE

116~145 ms
S
A
T
Figure5.1: Schematicofpulsesequencefordual3Dperfusionimaging. Thepulsesequence
consists of two sets of saturation RF pulse followed by a delay (T
SR
) and 3DFT GRE
readout. Thecentersoftwodataacquisitionswerelocatedatend-systoleandmid-diastole.
The largest T
SR
was used with a maximum threshold of 120 ms, depending on a given
heart rate.
5.2.2 Imaging Protocols
The imaging parameters are summarized in Table 5.1. We used a shorter data acquisition
window for systolic imaging due to its shorter quiescent period. We used a lower in-plane
resolution and a smaller number of partition slices. This choice was motivated by the
fact that, in general, systolic myocardium is transmurally thicker and is shorter along the
long axis. The number of partition slices also depended on the subject’s heart rate with
a threshold of 85-90 bpm. Saturation recovery time T
SR
, defined as the time interval
between the end of saturation RF pulse and the start of data acquisition, was limited by
the R-R interval. The pulse sequence used the longest possible T
SR
within a maximum
Table 5.1: Imaging parameters of dual 3D MPI.
Systolic imaging Diastolic imaging
Low HR / High HR Low HR / High HR
Matrix size 62£42£8 / 62£42£6 100£66£10 / 100£66£8
FOV (cm
2
) 28£28£8 / 28£28£6 28£28£10 / 28£28£8
Acquisition time (ms) 144.8 / 115.8 304.7 / 253.9
TR (ms) 2.1 2.3
TE (ms) 0.9 1.0
flip angle 12
±
12
±
77
value of 120 ms. The acquisition of systolic and diastolic data in the same R-R, limits the
saturation recovery time for diastolic imaging, resulting in slightly reduced CNR. Based
on Bloch simulation, the CNR reduction is less than 20% for heart rates between 50 and
120 bpm. The systolic scan volume is shifted toward the apex by 8-10 mm relative to the
diastolic scan volume.
3D MPI experiments were performed in five healthy subjects (Sex: 3F/2M, Age: 21-
55). Prior to perfusion scans, SSFP cine images were obtained for the determination of
the centers of the two acquisition windows. Contrast agent (0.05 mmol/kg Gd-DTPA,
Magnevist) was injected at a rate of 4 ml/s 0-5 second before the start of scan. Subjects
were instructed to hold their breath as long as possible.
5.2.3 Image Reconstruction and Analysis
All image reconstruction and subsequent analysis were implemented offline in MATLAB
(Mathworks, Natick, MA). Two-dimensional TSENSE was used for image reconstruction
from undersampled 3DFT perfusion data from the two cardiac phases [41,102]. The data
from six neighboring heart beats were combined to form coil sensitivity maps.
All perfusion images were normalized by pre-contrast images obtained by combining
data from the first six cardiac cycles, to correct for the effects of different T
SR
and re-
ceiver coil sensitivities. From the normalized images, endocardial and epicardial contours
were manually segmented. Myocardium was then automatically subdivided into endo-
and epi-cardial layers and six azimuthal sectors to generate time intensity curves (TICs)
from a total of twelve segments per slice. Upslope and peak enhancement values were
78
calculated for each myocardial segment, and were normalized by an upslope from the LV
to compensate for inter-subject variation.
Systolic and diastolic data sets were manually registered, and the four partition slices
closest to the mid-short axis level were used for TIC similarity analyses on a segment-
by-segment basis. First, correlation coefficient between systolic and diastolic TICs was
calculated to measure the similarity of overall dynamics in TICs. Second, TIC upslope
andpeakenhancementwerecomparedtomeasurethesimilaritybetweensemi-quantitative
perfusion indices. A total of 240 segments were analyzed in five subjects (48 segments per
subject).
5.3 Results
Table 5.2 contains the systolic and diastolic myocardial thicknesses measured in pixel
units. The transmural spatial resolution is higher in the anterior and inferior walls com-
pared to the septal and lateral walls, because the acquisition has asymmetric spatial
resolution, and the anterior-inferior direction is more closely aligned with the (higher
resolution) frequency encoding direction. The transmural resolution in the systolic my-
ocardium was comparable to that in the diastolic myocardium (3.35§0.52 vs 3.30§0.49 in
Table 5.2: Systolic and diastolic myocardial thicknesses, measured in pixel units.
Systolic myocardium Diastolic myocardium
Anterior/Inferior Septal/Lateral Anterior/Inferior Septal/Lateral
Subject1 3.10/3.41 1.89/2.10 3.00/4.00 1.65/2.31
Subject2 3.72/4.34 2.94/2.52 3.50/3.50 1.98/2.64
Subject3 3.72/3.41 2.94/2.31 2.50/4.00 2.97/2.97
Subject4 3.41/3.10 1.89/2.52 3.00/3.50 1.98/2.64
Subject5 2.64/2.64 1.89/1.68 2.75/3.25 1.82/1.98
All 3.35§0.52 3.35§0.52 3.35§0.52 3.35§0.52
79
anteriorandinferiorwalls, 2.27§0.45vs2.29§0.48inseptalandlateralwalls), presumably
due to the increased myocardial thickness at end-systole.
Figure 5.2 and Fig. 5.3 show representative systolic and diastolic 3D perfusion images
from one subject, at pre-contrast, RV enhancement, LV enhancement, and myocardial
enhancement. Perfusion images from both cardiac phases clearly show the arrival and
the passage of contrast agent with complete myocardial coverage and high image quality.
Systolic images suffer slightly from blurring and dark-rim artifacts presumably due to
lower in-plane spatial resolution. Dark rim in diastolic images is thinner than that in
systolic images, but its ratio to myocardial thickness is comparable to the ratio in systolic
images.
Figure 5.4 shows representative systolic and diastolic TICs from the four partition
slices closest to mid-short axis level. All TICs show uniform enhancement, except ones
in a diastolic basal slice which suffers slightly from fluctuation indicated by the gray
arrow. This fluctuation temporally corresponds to rapid signal elevation in LV, which
causesaliasingartifactsincoilsensitivitymapestimatedbyslidingwindowreconstruction.
Many endo-cardial TICs show a signal dip (see black arrows in Fig. 5.4) that is temporally
aligned with peak LV enhancement, and is likely to be related to dark rim artifact.
Figure 5.5 shows the average and standard deviation of the correlation coefficients
between systolic and diastolic TICs. The mean and standard deviation of the correlation
coefficient were 0.9841 and 0.0166 across all 240 segments, and mean values were larger
than 0.97 in all five subjects, indicating excellent correlation between temporal dynamics
of the two TICs.
80
Figure 5.2: Representative systolic 3D perfusion images from a healthy volunteer at pre-contrast, RV enhancement, LV enhance-
ment, and myocardial enhancement. The passage of contrast agent is clearly captured with complete myocardial coverage. The
images suffer slightly from blurring and dark-rim artifacts presumably due to lower in-plane spatial resolution.
81
Figure 5.3: Representative diastolic 3D perfusion images from the same volunteer as in Fig. 5.2 at pre-contrast, RV enhancement,
LV enhancement, and myocardial enhancement. Dark rim is thinner than that in systolic images, but its ratio to myocardial
thickness is comparable to the ratio in systolic images
82
Cardiac cycles
Signal intensity [a.u.] Signal intensity [a.u.]
slice #2 slice #3 slice #4 slice #5
Cardiac cycles Cardiac cycles Cardiac cycles
a
b
0 10 20 30 40
0
2
4
6
8
LV
ANT ENDO
SEP ENDO
POS ENDO
ANT EPI
SEP EPI
POS EPI
0 10 20 30 40
0
2
4
6
8
0 10 20 30 40
0
2
4
6
8
0 10 20 30 40
0
2
4
6
8
0 10 20 30 40
0
2
4
6
8
0 10 20 30 40
0
2
4
6
8
0 10 20 30 40
0
2
4
6
8
0 10 20 30 40
0
2
4
6
8
slice #3 slice #4 slice #5 slice #6
Figure 5.4: Representative time intensity curves based on (a) systolic and (b) diastolic perfusion images from a healthy subject.
Homogeneous signal enhancement is observed in nearly all myocardial segments. Diastolic slice #3 suffered from flickering artifacts
caused by errors in coil sensitivity map estimation (see gray arrow). Some endo-cardial TICs exhibit a signal dip (see black arrows)
that is temporally aligned with peak LV enhancement, and is likely to be related to dark rim artifact.
83
subj. 1 subj. 2 subj. 3 subj. 4 overall
Correlation coefficient
0.8
0.85
0.9
0.95
1
1.05
subj. 5
Figure 5.5: Correlation coefficients between systolic and diastolic TICs. Mean and stan-
dard deviation of the coefficient were 0.9841 and 0.0166 across a total of 240 segments in
five subjects (48 segments per subject), and mean coefficient was larger than 0.97 in all
five subjects, which validates strong similarity between overall dynamics in the two sets
of TICs.
Figure 5.6 contains scatter plots of the TIC upslope and peak enhancement values of
thesystolicanddiastolicdata, andresultsoflinearfitting. Thelinearcorrelationsbetween
the two sets of upslope and peak enhancement values were significant (P <0.001 in both).
Figure 5.7 compares the TIC upslope and peak enhancement values for endo- and epi-
cardial layers. There was no significant difference in either upslope or peak enhancement
between the endo- and epi-cardial layers. This was true for both systolic and diastolic
data across all subjects (at the significance level= 0.05). The average upslope was slightly
larger in the endo-cardial layer, likely due to dark rim artifact that is reduced as LV signal
decreases (see Fig. 5.4). The ratio of average endo- to epi-cardial upslopes was 1.0539, and
1.0741 in the systolic and diastolic TICs, respectively. The average peak enhancement
was slightly smaller in the endo-cardial layer, likely due to the residual effects of elevated
84
LV blood pool signal. The ratio of average endo- to epi-cardial peak enhancement values
was 0.9557 and 0.9580 in the systolic and diastolic TICs, respectively.
5.4 Discussion
We have demonstrated the feasibility of systolic 3D MPI, and the equivalence of TICs
derived from systolic versus diastolic acquisitions, in healthy subjects. Region-based time
intensity curves from systolic images were shown to be highly correlated with those from
diastolic images in terms of overall temporal dynamics, upslope, and peak enhancement.
Wefoundnostatisticallysignificantdifferencebetweentheupslopeandpeakenhancement
values obtained from endo- and epi- cardial segments. On average, upslope was slightly
larger and peak enhancement was slightly lower in endo-cardial layer than in epi-cardial
layer, likely due to dark rim artifact. This inter-layer difference was similar for both
cardiac phases, presumably due to similar myocardial thickness in pixel units.
One limitation of this study is the use of relativelylowspatial resolution, which is more
likely to suffer from dark rim artifacts and may be incapable of identifying small perfusion
defects. This is an important consideration, and needs to be investigated in patients
with known perfusion defects. At this time, dark rim artifacts are the most significant
limitation to the clinical acceptance of MR MPI, which underscores the need to achieve
higher spatial resolution. This may be achieved through the use of higher acceleration,
including use of high-density coil arrays and high-rate parallel imaging, and/or regularized
reconstruction [7,55].
85
Upslope
diast
Upslope
sys
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
0
0.05
0.1
0.15
0. 2
0.25
0. 3
0.35
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
subject 1 (Endo-cardial)
subject 2
subject 3
subject 4
subject 5
subject 1 (Epi-cardial)
subject 2
subject 3
subject 4
subject 5
Peak enhancement
diast
Peak enhancement
sys
a b
Figure 5.6: Scatter plot of time intensity curve (a) upslope and (b) peak enhancement for diastolic and systolic acquisitions.
Linear fitting across all myocardial segments results in y=0.9768x (P < .001) and y=0.9728x (P < .001) for upslope and peak
enhancement, respectively.
86
Upslope (a.u.)
a
c
0
0.05
0.1
0.15
0.2
0.25
0
0.05
0.1
0.15
0.2
0.25
Endocardial
Epicardial
Systolic TIC Diastolic TIC
0
0.5
1
1.5
2
2.5
3
0
0.5
1
1.5
2
2.5
3
Peak enhancement (a.u.)
subj. 1 subj. 2 subj. 3 subj. 4 overall subj. 5
subj. 1 subj. 2 subj. 3 subj. 4 overall subj. 5
subj. 1 subj. 2 subj. 3 subj. 4 overall subj. 5
subj. 1 subj. 2 subj. 3 subj. 4 overall subj. 5
Endocardial
Epicardial
b
d Systolic TIC Diastolic TIC
Upslope (a.u.) Peak enhancement (a.u.)
Figure 5.7: Comparison of time curve upslope and peak enhancement from endo-cardial and epi-cardial layers: (a) systolic upslope,
(b) diastolic upslope, (c) systolic peak enhancement, and (d) diastolic peak enhancement. There was no significant difference in
either upslope or peak enhancement between the endo- and epi-cardial layers. This was true for both systolic and diastolic data
across all subjects (at the significance level = 0.05). The ratio of endo- to epi-cardial upslopes was 1.0539 and 1.0741 in the systolic
and diastolic TICs, respectively. The ratio of average endo- to epi-cardial peak enhancement values was 0.9557 and 0.9580 in the
systolic and diastolic TICs, respectively.
87
The rationale for using lower in-plane resolution for systolic imaging was based on my-
ocardial wall thickening and shorter quiescent duration. Across the five healthy subjects,
the average systolic and diastolic wall thicknesses were 15.1 mm and 9.5 mm, respectively.
However, patients with severe irreversible defects may suffer from reduced systolic wall
thickening [72], which will decrease transmural spatial resolution in affected segments.
The incidence and degree of wall thickening impairment should be investigated in pa-
tients with ischemic heart disease. The duration of acquisition time for the two cardiac
phases was empirically chosen (see Table 5.1). Despite the use of the two most stable
cardiac phases, these acquisitions are longer than individual image acquisition times used
with the traditional 2D multi-slice approach, and may experience more significant motion
artifacts. Further investigation is needed to characterize and minimize motion artifacts in
3D MPI.
Flickering artifacts were seen in both systolic and diastolic perfusion images, especially
when viewed in video format. This artifact is caused by errors in coil sensitivity maps
that are derived from six neighboring cardiac cycles, and thus suffers from aliasing during
rapid signal change. Across temporal frames, the flickering is the most severe during the
arrival of contrast agent into the RV and LV. Across partition slices, it was the most
severe in apical slices whose coil maps suffer from aliasing of the large blood pool in
basal slices, as seen by fluctuation in TICs (Fig. 5.4). This had a negligible affect on
the overall myocardial TIC upslopes, due to the high frequency nature of the flickering.
Variable density k-space trajectories and key-hole reconstruction for coil map generation
may reduce and diffuse the artifacts over the FOV [39,99].
88
The TIC upslope has shown variation across myocardial segments and subjects in both
systolic (0.1699§0.0223) and diastolic data (0.1717§0.0212). The corresponding percent-
age standard deviation (13.1% and 12.3% for systolic and diastolic upslope, respectively)
was lower than standard deviation reported in previous studies (18.1%-25.8%, [70,91]).
Potential reasons for this variation include spatial variation in image intensity due to
receiver coil sensitivity, transmit RF inhomogeneity, and the slab excitation profile. Nor-
malizationbypre-contrastimagesshouldhavecorrectedformostofthevariation, yetwere
not perfect. The diagnostic value of the TIC upslope as a semi-quantitative perfusion in-
dex will be valid only if the upslope distribution of normal and defected myocardium is
well separated. This separability should be further confirmed in patients with perfusion
defects.
5.5 Conclusion
The feasibility of systolic 3D MPI has been demonstrated in healthy subjects. Myocardial
signal enhancement captured by systolic 3D imaging has been shown to be comparable
to diastolic 3D imaging by region-based TIC analysis. The mean and standard deviation
of the correlation coefficients between systolic and diastolic TICs were 0.9841 and 0.0166,
and the linear correlation between the two sets of TIC upslope and peak enhancement
values was statistically significant (P <0.001). Systolic 3D MPI may be advantageous in
subjects with severe R-R variability.
89
Chapter 6 :
Three-Dimensional Myocardial Perfusion Imaging with
Regularization
6.1 Introduction
Three-dimensional MR MPI is a promising alternative to 2D multi-slice approach, due
to its potential to accurately estimate the volume of perfusion defects. One challenging
requirement in 3D MR MPI is the need for high rate acceleration technique to limit data
acquisition time to 250-300 ms at end-diastole or 150-200 ms at end-systole. In previ-
ous chapters, six-fold acceleration was achieved using 2D sensitivity encoding (SENSE)
reconstruction, yet may not be sufficient for the following reasons.
First, in-plane spatial resolution of 3.0£4.5 mm
2
in diastolic images and 4.5£6.7 in
systolic images may be insufficient for visualizing subendocardial defects. Higher spatial
resolution may also help to avoid the dark rim artifacts which were seen on most of
volunteer images. Gibbs’ ringing is one of critical reasons for the artifacts, and the use of
high spatial resolution has been shown to restrict its transmural extent [20,76]. Second,
the data acquisition time need to be shortened to avoid potential motion artifacts in
90
subjects with high heart rates. The motion artifacts will be problematic particularly in
stress perfusion tests where the injection of Adenosine typically increases the heart rate.
Sparsity of a natural image itself is an important source of redundancy to be explored
for higher acceleration. This prior information can be utilized by regularization method
which improves the conditioning of an inverse problem by penalizing solution with rough
variation in pixel values. L
2
norm of difference between neighboring pixel values has been
the most frequently used penalty due to the existence of closed form solution and its
fast implementation [47,53,54], but typically involves image blurring. L
1
norm of finite
difference, so-called total variation (TV) has been recently introduced in MR community
[57,85], and has been shown to be effective in reconstruction from under-sampled k-space
while preserving edge information. Three-dimensional MPI is more advantageous than
2D multi-slice MPI in using the regularization method since sparsity increases in higher
dimensional images.
Sensitivity encoding (SENSE) is a well established acceleration method that has been
successfully applied to 3D MPI in previous chapters. SENSE reconstruction is formulated
as an inverse problem to which regularization penalty can be easily incorporated. In this
chapter, TV regularization SENSE reconstruction is developed for high rate acceleration
of 3D MPI to increase spatial resolution with whole heart coverage.
6.2 Theory
Total variation penalized SENSE reconstruction from under-sampled time-resolved 3D
k-space data can be formulated as the following optimization problem.
91
^ x=argmin
x
kF
u
Sx¡dk
2
+¸kxk
TV
(6.1)
where F
u
, S, x, d, and ¸ are under-sampled Fourier encoding, coil sensitivity weight-
ing, true images, under-sampled k-space measurement, and a regularization coefficient,
respectively. TV is l
1
norm of finite difference, defined as
kxk
TV
=
X
m;n
jD
y
(m;n)j+jD
z
(m;n)j (6.2)
where D
y
and D
z
denote finite difference in y and z directions, respectively.
D
y
(m;n)=x(m;n)¡x(m¡1;n)
D
z
(m;n)=x(m;n)¡x(m;n¡1) (6.3)
Since k-spaceis fully sampled along readout direction, all equations above need to be
solved independently per readout pixel.
Two requirements should be met for successful reconstruction. First, TV values of an
object to be solved should be truly sparse. Second, k-spaceshould be undersampled in a
random way, so that the aliasing artifacts due to the undersampling are incoherent.
6.3 Methods
All subsequent experiments were performed on a General Electric 3 Tesla scanner with 40
mT/m gradient amplitude and 150 T/m/s gradient slew rate, using a body-coil for RF
transmission and an eight-channel cardiac coil array for signal reception.
92
6.3.1 Pulse Sequence and Sampling Scheme
The pulse sequence consisted of saturation recovery 3DFT acquisition modules centered
at mid-diastole(Fig. 6.1). For the preparatory saturation, we used tailored hard-pulse
train that is robust in the presence of B0 and B1 inhomogeneities and involves relatively
low RF power [94]. A sinc-shaped RF pulse with a time-bandwidth-product of four was
used for slab selective excitation.
Three-dimensional k-spacewas undersampled in a pseudo-random fashion such that
aliasing is diffused over FOV while the data from neighboring cardiac cycles fully fill the
k-space. That is, the center of k-space(» 2 % of full k-space) was always sampled, and
the remaining k-spacewas randomly undersampled in a disjoint way across nine cardiac
cycles. The procedure of sampling function design is described in Fig. 6.2. The fully
sampled k-spacedata from neighboring nine cardiac cycles was used for coil sensitivity
map estimation.
Sampling functions should be determined such that resulting aliasing is incoherent as
much as possible. The design of the nine sampling functions was randomly repeated 1,000
times, and the best set of sampling functions was chosen such that the maximum value
of sidelobes in point spread function (PSF) is minimal [46]. Figure 6.3 shows PSFs of
the nine chosen sampling functions using a normal (top) and a harsher (bottom) display
windows. Signal intensities of sidelobes are evenly small throughout the FOV in all nine
PSFs. The maximum intensity of sidelobe was 17 % of the mainlobe.
Right before the perfusion scan, proton density weighted images were obtained using
the same imaging parameters except the saturation pulse turned off and low readout flip
93
S
A
T
Psuedo-randomly
undersampled
3DFT readout
T
SR
Figure6.1: Schematicofpulsesequenceforhighlyaccelerated3DperfusionimagingviaTVregularizedSENSE.Thepulsesequence
consists of saturation RF pulse followed by a delay (T
SR
) and eight-fold undersampled 3DFT readout. The center of k-spacewas
always sampled, and the remaining region was undersampled randomly, but in a disjoint way across nine cardiac cycles for coil
map estimation.
94
Y
cent cent cent
N S t s S = . . Choose
cent
S S =
) 0 (
Set
outer
i
outer
k
i
k
i
outer
N S t s
S U S
=
− ⊂
−
=
) (
) (
1
0
) (
. .
choose Randomly
U
9 if = i
) 9 , , 1 ( , 0
) (
L = = = n S i
n
φ
) (
1
0
) ( k
i
k
i
outer
S U S
U
−
=
− =
N
9 if < i
Y
N
1 + =i i
END
) ( ) ( i
outer cent
i
S S S U =
Figure 6.2: Flow chart of sampling function design. The center of k-space is initially
chosen, and is included in all nine sampling functions. Fori
th
sampling function, sampling
points are randomly chosen from k-space positions that are not included in 1
st
-(i¡1)
th
sampling functions. The final 9
th
sampling function is automatically determined as all
residual k-space points. N
cent
= number of samples in k-space center. N
outer
= number
of samples in outer k-space. S
cent
= k-space center that is always sampled. S
(i)
outer
= outer
k-space in i
th
sampling function. S
(i)
= i
th
sampling function.
95
a
b
¡ Figure6.3: Pointspreadfunctionoftheninechosensamplingfunctionwith(a)normaland
(b) harsh display window. The energy of sidelobe in all nine PSFs is evenly diffused over
the FOV, which improves the incoherence of aliasing caused by k-space undersampling.
angle (4
±
). The proton density images were reconstructed using the same key-hole type
method as in coil sensitivity estimation. Both proton weighted imaging and T1 weighted
perfusion imaging were performed during breath-hold at exhalation.
6.3.2 Scan Parameters
The imaging parameters used for in-vivo experiments are summarized in Table 6.1. The
number of acquired samples in the ky-kz plane was 116, corresponding to undersampling
factor of 7.9. The data acquisition time was 116*TR = 289 ms.
96
Perfusion experiments were performed in two healthy volunteers. Gadolinium-based
Contrast agent (0.05 mmol/kg) was injected at a rate of 5 ml/s 0¡5 second before the
start of scan. Volunteers were instructed to hold their breath as long as possible.
6.3.3 Image Reconstruction
All image reconstruction was performed off-line using MATLAB
R °
(Mathworks, Natick,
MA) on a Desktop computer equipped with 3.6 GHz Intel processor and 2GB RAM. Raw
k-spacedata was normalized by the maximum intensity of temporally averaged images.
This ensures consistent signal level throughout different scans so that the same regular-
ization coefficient can bring the same effects on the reconstruction. The 3D k-spacewas
1D Fourier transformed alongk
x
, and 2D reconstruction was performed ink
y
¡k
z
domain
independently per readout sample to reduce computational and memory load.
Full sampledk-space data for coil sensitivity estimation was obtained by key-hole type
reconstruction. For k
th
cardiac cycle, k-space data acquired at the current cardiac cycles
(S
(k)
), and disjoint samples in outer k-space from neighboring nine cycles ([
i=4
i=¡4
S
(k+i)
outer
),
were combined. Data inconsistency due to CA enhancement and motion can cause alias-
ing, yet will be limited to high frequency component since thek-space center was acquired
Table 6.1: Imaging parameters used for in-vivo experiments.
T
SR
120 ms
TR 2.5 ms
TE 1.2 ms
flip angle 12
±
matrix size 116£ 76£ 10
FOV 28£ 28£ 10 cm
3
97
at current cycle only. Resulting high frequency error in sensitivity estimation can be mit-
igated by the use of TV regularization.
The optimization problem written as Eq. 6.1 was iteratively solved using non-linear
conjugate gradient (CG) method [33]. Computing gradient of the cost function is a key
step in CG, but is not doable since the absolute value function in the cost function is not
differentiable. Thus the absolute function was approximated byjxj¼
p
¹ xx+¹ where ¹ x is
complex conjugate ofx, and¹ is a small positive parameter. Back-tracking line search was
used to expedite convergence of CG. Zero images were used for initial guess. Maximum
iteration number was set to 150.
After normalization by proton density weighted images, the 3D perfusion images were
manually segmented into six or four myocardial sectors. Time intensity curves (TIC) were
generated per segment.
6.4 Results
Figure 6.4 shows representative perfusion images reconstructed using regularization co-
efficient ¸ = 0.0005. The four rows correspond to pre-contrast, RV enhancement, LV
enhancement, and myocardial enhancement, from top to bottom, respectively. Aliasing
from high rate undersampling was effectively suppressed due to the combination of sen-
sitivity encoding and TV regularization. Blocky artifacts which would appear with too
large regularization coefficient is not seen.
Coilsensitivityestimatedusingthedatafromneighboringninecardiaccyclesisinaccu-
rate during rapid signal changes in blood pools, as mentioned in Chapter 4. Compared to
98
regular undersampling of lattice structure, variable density random undersampling results
in such reconstruction errors more diffused over the FOV.
Figure 6.5 shows L-curve analysis on the effect of regularization coefficient in two
volunteers. Ten incremental values of the coefficient was used, ranging from 0.00001
to 0.01. An increase in the regularization coefficient reduces TV values (or improves
smoothness) while sacrificing data fidelity. It is desired to choose a coefficient which
produces low values of both TV and the norm of data fit term. Graphically, this optimal
coefficient corresponds to the elbow of the L-shaped curve, and 0.0005 was chosen in this
study. It should be noted that the optimal coefficient is similar in the two volunteer data
sets although absolute values of TV and the norm of data fidelity are different.
Figure 6.6 shows representative perfusion images at peak myocardial enhancement,
reconstructed using different regularization coefficients (0.00005, 0.0002, 0.0005, 0.001,
0.01). The image from too small coefficient (¸=0.00005) suffer noisy artifacts, due to high
frequency aliasing from variable density undersampling, and errors in sensitivity estima-
tion. On the other hand, the image from too big coefficient (¸=0.01) over-smoothed, ex-
hibiting “blocky” artifacts. The images from¸=0.0005 are well balanced between smooth-
ness and noisy aliasing. This qualitative decision on the regularization coefficient agrees
with the findings based on L-curve analysis (Fig. 6.5).
Figure 6.7 shows representative myocardial time intensity curve (TIC) from 2
nd
,
4
th
,6
th
, and 8
th
partition slices. TICs were generated from the images reconstructed
using ¸=0.0005. Uniform signal enhancement is seen across all myocardial segments.
Using MATLAB on a Desktop computer with 3.6 GHz and 2 GB RAM, the non-linear
CG required approximately 3.5 hours for the reconstruction of 30 time frames.
99
Figure 6.4: Representative 3D perfusion images from a healthy volunteer at pre-contrast, RV enhancement, LV enhancement, and
myocardial enhancement. The passage and arrival of contrast agent is clearly captured with high SNR. Compared to the images
in Chapter4, blurring is reduced due to the increased spatial resolution.
100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
400
600
800
1000
1200
data fidelity
Total variation
subj #1
subj #2
λ = 0.00001
λ = 0.0002
λ = 0.0005
λ = 0.001
λ = 0.01
Figure 6.5: Effect of regularization coefficient on total variation and data fidelity. Ten
incremental values of the regularization coefficient was used, ranging from 0.00001 to
0.01. An increase in the coefficient reduces TV while sacrificing data fidelity. The L
curves suggest using 0.0005 with which TV is minimal while the loss of data fidelity is
negligible. The overall shape of the curve is similar in the two data sets.
6.5 Discussion
Three-dimensional first-pass MPI has been eight-fold accelerated using total variation reg-
ularized SENSE reconstruction. Aliasing from highly undersampledk-space data was well
suppressed by the synergetic combination of sensitivity encoding and TV regularization.
The regularization coefficient retrospectively chosen by L-curve analysis has reconstructed
images well balanced between data fidelity and smoothness.
Undersampled k-spacedata were acquired by fully sampling a pre-defined central re-
gion of k-spaceand randomly sampling the remaining outer part in a disjoint way across
101
λ = 0.00005
λ = 0.0002
λ = 0.0005
λ = 0.001
λ = 0.01
slice #3 slice #4 slice #5 slice #6 slice #7
Figure 6.6: Effect of regularization coefficient on image reconstruction. Each row shows
four representative partition slices reconstructed using different value of regularization
coefficients. The images from too small regularization coefficient suffer noisy aliasing
artifacts while the images from too heavy regularization suffer over-smoothing. ¸=0.0005
produces well balanced images, which agrees with the findings from L-curve analysis.
102
time (sec)
signal intensity (a.u.)
time (sec)
time (sec) time (sec)
6 segments
4 segments
ANT
ANT SEP
SEP
INF
POS
LAT
ANT
SEP
INF
LAT
slice #2 slice #4
slice #6 slice #8
0 10 20 30
0
0.2
0.4
0.6
0.8
1
0 10 20 30
0
0.2
0.4
0.6
0.8
1
0 10 20 30
0
0.2
0.4
0.6
0.8
1
0 10 20 30
0
0.2
0.4
0.6
0.8
1
signal intensity (a.u.)
Figure 6.7: Representative myocardial time intensity curves (TIC) from a healthy volun-
teer. TIC was generated from the images reconstructed using ¸=0.0005. Uniform signal
enhancement is seen across all myocardial segments.
cardiac cycles. In such way, full sampled data for coil map estimation could be obtained
by combining data from neighboring cardiac cycles while maintaining “variable density”
sampling pattern. The central region of k-space to be fully sampled was determined only
empirically. The optimal choice of it, if any, will depend on target acceleration rate and
matrix size, and is an open research question.
Finding the optimal regularization coefficient is challenging. We chose the coefficient
value retrospectively from L-curve analysis, and obtained the same value from the per-
fusion data of two volunteers [28]. The goodness of the chosen coefficient value was
investigated only qualitatively in healthy volunteers. The same investigation should be
103
performedinpatientswithperfusiondefects, includingacomparisonwithreferenceimages
which could be obtained by an additional 2D scan without acceleration.
The combination of data from nine cardiac cycles for coil sensitivity estimation is
another limitation of the proposed reconstruction method. In the presence of severe
breathing motion, the sensitivity estimation will be inaccurate, resulting in degradation
of image quality. To reduce the number of cardiac cycles required, strong quadratic
regularization can be used for the reconstruction of each coil image, due to the smooth
nature of coil sensitivities. Another possible approach would be joint estimation of coil
map and perfusion images [107], at the expense of much longer reconstruction time.
6.6 Conclusion
Highly accelerated 3D first-pass MPI was developed via TV regularized SENSE recon-
struction, and tested in healthy volunteers. Ten partition slices with in-plane spatial
resolution of 2.4 £ 3.6 mm
2
was achieved from 8-fold acceleration. The regularization
coefficient chosen by L-curve analysis produced perfusion images well balanced between
TV penalty and data fidelity. The effect of and optimal determination of the coefficient
should be further investigated in patients with perfusion defects.
104
Chapter 7 :
Summary and Future Work
Dynamic cardiac MRI is a promising non-invasive imaging tool for the detection and
assessment of ischemic heart disease. This dissertation presented MR imaging and recon-
struction methods for improved assessment of cardiac function and perfusion. The key
contributions of the dissertation are summarized as follow.
² The aliasing process caused by dynamic undersampled spiral imaging was described
both algebraically and graphically. Based on such description of the dynamic alias-
ing, time efficient iterative reconstruction method was developed to increase tem-
poral resolution. The proposed method was validated by numerical phantom study
and in-vivo scan, and was proven to be superior to conventional sliding window
reconstruction.
² Three-dimensional myocardial perfusion imaging (MPI) with complete left ventric-
ular coverage has been developed. Cardiac phantom study has shown that 3D MPI
outperforms conventional 2D multi-slice MPI in estimating the size of defects. The
feasibility of 3D MPI has been demonstrated in healthy volunteers, by showing
excellent diagnostic image quality with high CNR, and complete spatial coverage.
105
² The feasibility of systolic 3D MPI has been investigated with a hypothesis that
systolic acquisition may be more robust in the presence of severe R-R variability.
We have conducted the comparison of systolic and diastolic data acquisitions for
3D MPI, and have shown extremely high cross correlation between time intensity
curves (TIC) from the two acquisitions, and significant linear correlation between
the two TIC upslopes.
² Highlyaccelerated3DMPIviatotalvariationregularizedSENSEreconstructionhas
been developed and tested on healthy volunteers. Ten partition slices with in-plane
resolution of 2.4£ 3.6 mm
2
was achieved from 8-fold acceleration.
7.1 Future Work
The MR imaging and reconstruction methods developed in the previous chapters are
ultimately targeted towards clinical utility. Since all of the proposed methods have been
tested on healthy volunteers only, potential problems associated with patient scanning
need to be explored and resolved. This section summarizes potential improvements of the
proposed techniques and important clinical validation, organized by previous chapters.
7.1.1 Accelerated Spiral Cardiac Imaging
² Replace the assumption on finite support region with sensitivity encoding: The sup-
portregionofsolutionwasempiricallychosen, whichmaydegradethereconstruction
106
accuracy depending on subjects and imaging parameters. Substitution of coil sen-
sitivity for the finite support assumption will achieve the same goal of conditioning
while improving the stability of the proposed reconstruction.
² Apply the method to Dobutamine stress testing in patients with ischemic heart
disease: Dobutamine typically increases patient’s heart rate up to »150 beats per
minute, which will require imaging with high temporal resolution to avoid motion
artifacts.
7.1.2 Three-Dimensional Myocardial Perfusion Imaging
² Compare with conventional 2D multi-slice MPI in patients with perfusion defects:
The most important advantage of the proposed 3D MPI is its capability of accurate
estimation of defect size. The estimated extent of perfusion defects by 3D and 2D
multi-slice MPI should be compared in patients.
² Compare systolic and diastolic 3D MPI in patients with severe R-R variability or
arrhythmia: The superiority of systolic acquisition in 3D MPI should be validated
in patients with severe R-R variability. Potential degradation of diagnostic perfor-
mance due to reduced spatial resolution should be also investigated.
² Perform a joint analysis of whole heart MPI and late Gadolinium enhancement
(LGE) imaging: Complete left ventricular coverage of 3D MPI enables ideal com-
bination with LGE imaging that is a gold standard method for scar imaging. The
combined analysis can completely determine the distribution of normal, reversibly
ischemic, and infarcted myocardium.
107
7.1.3 Highly Accelerated 3D Myocardial Perfusion Imaging
² Investigate the effect of total variation regularization on the depiction of perfusion
defects: A range of regularization coefficient values should be used for the proposed
reconstruction in patients with perfusion defects. Potential diagnostic degradation
caused by the regularization should be carefully examined.
108
References
[1] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions with For-
mulas, Graphs, and Mathematical Tables. Dover, Newyork, 1972.
[2] S. Bagchi and S. Mitra. The nonuniform discrete Forier transform and its applica-
tions in signal processing. Kluwer, 1999.
[3] H. Benjelloun, R. Itti, L. Philippe, J. M. Lorgeron, and M. Brochier. Beat-to-
beat assessment of left ventricular ejection in atrial fibrillation. Eur. J. Nucl. Med.,
8:206–210, 1983.
[4] P. Bernhardt, T. Engels, B. Levenson, K. Haase, A. Albrecht, and O. Strohm.
Prediction of necessity for coronary artery revascularization by adenosine contrst-
enhanced magnetic resonance imaging. Int. J. Card., 112(2):184–190, 2006.
[5] M. Blaimer, F. Breuer, M. Mueller, R. M. Heidemann, M. A. Griswold, and R. M.
Jakob. SMASH, SENSE, PILS, GRAPPA: How to choose the optimal method. Top
Magn Reson Imaging, 15:223–236, 2004.
[6] F.Bloch, W.W.Hansen, andM.E.Packard. Nuclearinduction. Phys. Rev., 69:127,
1946.
[7] K. T. Block, M. Uecker, and J. Frahm. Undersampled radial MRI with multiple
coils.Iterativeimagereconstructionusingatotalvariationconstraint. Magn. Reson.
Med., 57:1086–1098, 2007.
[8] R. N. Bracewell. The Fourier Transform and Its Applications. McGraw-Hill,
Boston, 1999.
[9] K. A. Brown, C. A. Baucher, R. D. Okada, T. E. Guiney, J. B. Newell, H. W. Straus,
and G. M. Pohost. Prognostic value of exercise thallium201 imaging in patinets
presenting for evaluation of chest pain. J. Amer. Coll. Cardiol., 1:994–1001, 1983.
[10] M. E. Brummer, D. Moratal-Perez, C.-Y. Hong, R. I. Pettigrew, J. Millet-Roig, and
W. T. Dixon. Noquist: Reduced field-of-view imaging by direct fourier inversion.
Magn. Reson. Med., 51:331–342, 2004.
[11] J. Bundy, O. Simonetti, G. Laub, and J. P. Finn. Segmented true FISP cine imaging
of the heart. In Proc., ISMRM, 7th Annual Meeting, page 1282, Philadelphia, 1999.
109
[12] E. Candes, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal
reconstruction from highly incomplete frequency information. IEEE Trans. Inf.
Theory, 52:489–509, 2006.
[13] M. D. Cerqueira, N. J. Weissman, V. Dilsizian, A. K. Jacobs, S. Kaul, W. K. Laskey,
D.J.Pennell,J.A.Rumberger,T.Ryan,andM.S.Verani. Standardizedmyocardial
segmentation and nomenclature for tomographic imaging of the heart: a statement
for healthcare professionals from the cardiac imaging committee of the council on
clinical cardiology of the american heart association. Circulation, 105(4):539–542,
2002.
[14] S. Conolly, D. G. Nishimura, A. Macovski, and G. Glover. Variable-rate selective
excitation. J. Magn. Reson., 78(3):440–458, 1988.
[15] C. D. Constantinides, E. Atalar, and E. R. McVeigh. Signal-to-noise measurements
in magnitude images from NMR phased arrays. Magn. Reson. Med., 38:852–857,
1997.
[16] A. Cuocolo, W. Acampa, M. Imbriaco, N. De Luca, G. L. Iovino, and M. Salvatore.
The many ways to myocardial perfusion imaging. Q. J. Nucl. Med. Mol. Imaging,
49:4–18, 2005.
[17] D. L. Donoho. Compressed sensing. IEEE Trans. Inf. Theory, 52:1289–1306, 2006.
[18] J. P. Earls, V. B. Ho, T. K. Foo, E. Castillo, and S. D. Flamm. Cardiac MRI: recent
progress and continued challenges. J. Magn. Reson. Imaging, 16:111–127, 2002.
[19] J. P. Earls, V. B. Ho, T. K. Foo, E. Castillo, and S. D. Flamm. Cardiac MRI: recent
progress and continued challenges. J. Magn. Reson. Imaging, 16:111–127, 2002.
[20] D.L. Parker E.V.R. Di Bella and A.J. Sinusas. On the dark rim artifact in dynamic
contrast-enhanced MRI myocardial perfusion studies. Magn. Reson. Med., 54:1295–
1299, 2005.
[21] J. Felblinger, C. Lehmann, and C. Boesch. Electrocardiogram acquisition during mr
examinations for patient monitoring and sequence triggering. Magn. Reson. Med.,
32:523–529, 1994.
[22] J. A. Fessler and B. P. Sutton. Nonuniform fast Forier transform using min-max
interpolation. IEEE Trans. Signal Proc, 51:560–574, 2003.
[23] C. D. Furberg, B. M. Psaty, T. A. Manolio, J. M. Gardin, V. E. Smith, and P. M.
Rautaharju. Prevalence of atrial fibrillation in elderly subjects (the cardiovascular
health study). Am. J. Cardiol., 74:236–241, 1994.
[24] B. L. Gerber, S. V. Raman, K. S. Nayak, F. H. Epstein, P. Ferreira, L. Axel,
and D. L. Kraitchman. Myocardial first-pass perfusion cardiovascular magnetic
resonance: history, theory, and current state of art. J. Cardiovasc. Magn. Reson.,
10:18, 2008.
110
[25] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang,
B. Kiefer, and A. Haase. Generalized autocalibrating partially parallel acquisitions
(GRAPPA). Magn. Reson. Med., 47:1202–1210, 2002.
[26] R. Hachamovitch, D. S. Berman, L. J. Shaw, H. Kiat, I. Cohen, J. A. Cabico,
J. Friedman, and G. A. Diamond. Incremental prognostic value of myocardial per-
fusion single photon emission computed tomography for the prediction of cardiac
death. Circulation, 97:535–543, 1998.
[27] M. S. Hansen, , C. Baltes, J. Tsao, S. Kozerke, K. P. Pruessmann, and H. Eggers.
k-t BLAST reconstruction from non-Cartesian k-t space sampling. Magn. Reson.
Med., 55:85–91, 2006.
[28] P. C. Hansen. The L-curve and its use in the numerical treatment of inverse prob-
lems. In: P. Johnston, editor. Computational Inverse problems in Electrocardiology.
WIT Press, Southampton, UK, 2001.
[29] B.A.Hargreaves, C.H.Cunningham, D.G.Nishimura, andS.M.Conolly. Variable-
rate selective excitation for rapid MRI sequences. Magn. Reson. Med., 52:590–597,
2004.
[30] O. Heid. True FISP cardiac fluoroscopy. In Proc., ISMRM, 5th Annual Meeting,
page 320, Vancouver, 1997.
[31] R. M. Henkelman. Measurement of signal intensities in the presence of noise in MR
images. Med. Phys., 12:232–233, 1985.
[32] B. Hesse, K. Tagil, A. Cuocolo, C. Anagnostopoulos, M. Bardies, J. Bax, F. Bengel,
E. Busemann Sokole, G. Davies, M. Dondi, L. Edenbrandt, P. Franken, A. Kjaer,
J. Knuuti, M. Lassmann, M. Ljungberg, C. Marcassa, P. Y. Marie, F. McKiddie,
M. O’Connor, E. Prvulovich, R. Underwood, and B. van Eck-Smit. EANM/ESC
procedural tuideline for myocardial perfusion imaging in nuclear cardiology. Nucl.
Med. Mol. Imaging, 32:855–897, 2005.
[33] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear
systems. 49:409–436, 1952.
[34] R. D. Hoge, R. K. S. Kwan, and G. B. Pike. Density compensation functions for
spiral MRI. Magn. Reson. Med., 38:117–128, 1997.
[35] X. Hu and T. Parrish. Reduction of field of view for dynamic imaging. Magn.
Reson. Med., 31:691–694, 1994.
[36] P. Hunold, T. Schlosser, and J. Barkhausen. Magnetic resonance cardiac perfusion
imaging-a clinical perspective. 16:1779–1788, 2006.
[37] J. Jackson, C. H. Mayer, D. G. Nishimura, and A. Macovski. Selection of a con-
volution function for Fourier inversion using gridding. IEEE Trans. Med. Imaging,
10:473–478, 1991.
111
[38] K. R. Johnson, A. Whigham S. J. Patel, A. Hakim, R. I. Pettigrew, and J. N.
Oshinski. Three-dimensional, time-resolved motion of the coronary arteryies. J.
Magn. Reson. Imaging, 6:663–673, 2004.
[39] R. A. Jones, O. Haraldseth, T. B. Muller, P. A. Rinck, and A. N. Oksendal. K-space
substitution: a novel dynamic imaging technique. Magn. Reson. Med., 29:830–834,
1993.
[40] P. Kellman and A. E. Arai. Imaging sequences for first pass perfusion-a review. J.
Cardiovasc. Magn. Reson., 9:525–537, 2007.
[41] P. Kellman, F. H. Epstein, and E. R. McVeigh. Adaptive sensitivity encoding
incorporating temporal filtering TSENSE. Magn. Reson. Med., 45:846–852, 2001.
[42] P. Kellman, Q. Zhang, A. C. Larson, O. P. Simonetti, E. R. McVeigh, and A. E.
Arai. Cardiac first-pass perfusion MRI using 3D trueFISP parallel imaging using
TSENSE. In Proc., ISMRM, 12th Annual Meeting, page 310, Kyoto, 2004.
[43] A. B. Kerr, J. M. Pauly, B. Hu, K. C. Li, C. J. Hardy, C. H. Meyer, A. Macovski,
and D. G. Nishimura. Real-time interactive MRI on a conventional scanner. Magn.
Reson. Med., 38:355–367, 1997.
[44] D. Kim and L. Axel. Multislice, dual-imaging sequence for increasing the dynamic
range of the contrast-enhanced blood signal and CNR of myocardial enhancement
at 3T. J. Magn. Reson. Imaging, 23:81–86, 2006.
[45] D. Kim, O. Gonen, N. Oesingmann, and L. Axel. Comparison of the effectiveness
of saturation pulses in the heart at 3T. Magn. Reson. Med., 59:209–215, 2008.
[46] Y. C. Kim and K. S. Nayak. Optimization of undersampled variable density spiral
trajectories based on incoherence of spatial aliasing. page 422, 2008.
[47] K. F. King. SENSE image quality improvement using matrix regularization. page
1771, 2001.
[48] I. Klem, J. F. Heitner, D. J. Shah, M. H. Sketch, V. Behar, J. Weinsaft, P. Cawley,
M. Parker, M. Elliott, R. M. Judd, and R. J. Kim. Improved detection of coronary
artery disease by stress perfusion cardiovascular magnetic resonance with the use of
delayed enhancement infarction imaging.
[49] K. Knesaurek and J. Machac. Comparison of
18
f spect with pet in myocardial
imaging: A realistic thorax-cardiac phantom study. BMC Nucl. Med., 6:5, 2006.
[50] C. M. Kramer. Cardiovascular MRI: The future is now. Electromedica, 69(2):1–6,
2001.
[51] R. Y. Kwong. Cardiovascular Magnetic Resonance Imaging. Humana Press, 2007.
[52] G. W. Lenz, E. M. Haacke, and R. D. White. Retrospective cardiac gating: a review
of technical aspects and future directions. Magn. Reson. Med., 7:445–455, 1989.
112
[53] Z.-P. Liang, R. Bammer, J. Ji, N. Pelc, and G. Glover. Making better SENSE:
wavelet de-nosing, Tikhonov regularization, and total-least squares. page 2388,
2002.
[54] F.-H. Lin, K. K. Kwong, J. W. Belliveau, and J. J. Wald. Parallel imaging recon-
struction using automatic regularization. Magn. Reson. Med., 51:559–567, 2004.
[55] B. Liu, K. King, M. Steckner, J. Xie, J. Sheng, and L. Ying. Regularized sensitiv-
ity encoding (sense) reconstruction using bregman iterations. Magn. Reson. Med.,
61:145–152, 2009.
[56] M. Lustig, D. L. Donoho, and J. M. Pauly. Sparse MRI: The application of com-
pressed sensing for rapid MR imaging. Magn. Reson. Med., 58:1182–1195, 2007.
[57] M. Lustig, J. H. Lee, and D. L. Donoho. Faster imaging with randomly perturbed,
undersampled spiral and l
1
reconstruction. In Proc., ISMRM, 13th Annual Meeting,
page 685, Miami, 2005.
[58] B. Madore, G. H. Glover, and N. J. Pelc. Unaliasing by Fourier-encoding the
overlaps using the temporal dimension (UNFOLD), applied to cardiac imaging and
fMRI. Magn. Reson. Med., 42:813–828, 1999.
[59] T. H. Marwick, C. Case, S. Sawada, C. Rimmerman, P. Brenneman, R. Kovacs,
L.Short,andM.Lauer. Predictionofmortalityusingdobutamineechocardiography.
J. Amer. Coll. Cardiol., 37:754–760, 2001.
[60] C. D. Mathers and D. Loncar. Projections of global mortality and burden of disease
from 2002 to 2030. 3:e442, 2006.
[61] C. H. Meyer, B. S. Hu, D. G. Nishimura, and A. Macovski. Fast spiral coronary
artery imaging. Magn. Reson. Med., 28:202–213, 1992.
[62] E. Nagel, C. Klein, I. Paetsch, S. Hettwer, B. Schnackenburg, K. Wegscheider, and
E. Fleck. Magnetic resonance perfusion measurements for the noninvasive detection
of coronary artery disease. Circulation, 108:432–437, 2003.
[63] E. Nagel, C. Klein, I. Paetsch, S. Hettwer, B. Schnackenburg, K. Wegscheider,
and E. Fleck. Heart disease and stroke statistics 2009 update. a report from the
american heart association statistics committee and stroke statistics subcommittee.
Circulation, 119:e1–e161, 2009.
[64] K. S. Nayak, B. A. Hargreaves, B. S. Hu, D. G. Nishimura, J. M. Pauly, and C. H.
Meyer. Spiral balanced steady-state free precession cardiac imaging. Magn. Reson.
Med., 53:1468–1473, 2005.
[65] K. S. Nayak and B. S. Hu. The future of real-time cardiac magnetic resonance
imaging. 7:45–51, 2005.
[66] R. W. Nesto and G. J. Kowalchuk. The ischemic cascade: Temporal sequence of
hemodynamic, electrocardiographic and symptomatic expressions of ischemia. Am.
J. Cardiol., 59:C23–C30, 1987.
113
[67] D. C. Noll, D. G. Nishimura, and A. Macovski. Homodyne detection in magnetic
resonance imaging. IEEE Trans. Med. Imaging, 10(2):154–163, June 1991.
[68] A. Oppelt, R. Graumann, H. Barfuss, H. Fischer, W. Hartl, and W. Shajor. FISP-a
new fast MRI sequence. Electromedica, 54:15–18, 1986.
[69] J. D. O’Sullivan. A fast sinc function gridding algorithm for Fourier inversion in
compter tomography. IEEE Trans. Med. Imaging, 4:200–207, 1985.
[70] J. R. Panting, P. D. Gatehouse, G.-Z. Yang, F. Grothues, D. N. Firmin, P. Collins,
and D. J. Pennell. Abnormal subendocardial perfusion in cardiac stndrome x de-
tected by cardiovascular magnetic resonance imaging. N. Engl. J. Med., 346:1948–
1953, 2002.
[71] J. R. Panting, P. D. Gatehouse, G. Z. Yang, F. Grothues, D. N. Firmin, P. Collins,
and D. J. Pennell. Abnormal subendocardial perfusion in cardiac syndrome X de-
tected by cardiovascular magnetic resonance imaging. N. Engl. J. Med., 346:1948–
1953, 2002.
[72] P. Perrone-Filardi, L. Bacharach, V. Dilsizian, S. Maurea, J. A. Frank, and R. O.
Bonow. Regional left ventricular wall thickening. relation to regional uptake of
18Fluorodeoxyglucose and 201tl in patients with chronic coronary artery disease
and left ventricular dysfunction. Circulation, 86:1125–1137, 1992.
[73] R. I. Pettigrew, J. N. Oshinski, G. Chatzimavroudis, and W. T. Dixon. MRI tech-
niques for cardiovascular imaging. J. Magn. Reson. Imaging, 10:590–601, 1999.
[74] J. G. Pipe and P. Menon. Sampling density compensation in MRI: rationale and
an iterative numerical solution. Magn. Reson. Med., 41:179–186, 1999.
[75] S. Plein, A. Radjenovic, J. P. Ridgway, D. Barmby, J. P. Greenwood, S. G. Ball,
and M. U. Sivananthan. Coronary artery disease: myocardial perfusion mr imaging
with sensitivity encoding versus conventional angiography. Radiology, 235(2):423–
430, 2005.
[76] S. Plein, S. Ryf, J. Schwitter, A. Radjenovic, P. Boesiger, and S. Kozerke. Dynamic
contrast-enhanced myocardial perfusion MRI accelerated with k-tSENSE. Magn.
Reson. Med., 58:777–785, 2007.
[77] S.K.PlevritisandA.Macovski. Spectralextrapolationofspatiallyboundedimages.
IEEE Trans. Med. Imaging, 14(3):487–497, Sep 1995.
[78] M. Poustchi-Amin, F. R. Gutierrez, J. J. Brown, S. A. Mirowitz, V. R. Narra,
N. Takahashi, and P. K. Woodard. Performing cardiac MR imaging: an overview.
Magn. Reson. Imaging Clin. North Am., 11:1–18, 2003.
[79] K. P. Pruessmann, M. Weiger, P. Bornert, and P. Boesiger. Advances in sensitivity
encoding with arbitrary k-space trajectories. Magn. Reson. Med., 46(4):638–651,
2001.
114
[80] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger. SENSE: sensi-
tivity encoding for fast MRI. Magn. Reson. Med., 42:952–962, 1999.
[81] E. M. Purcell, H. C. Torrey, and R. V. Pound. Resonance absorption by nuclear
magnetic moments in a solid. Phys. Rev., 69:37, 1946.
[82] V. Rasche, R. Proska, R. Sinkus, P. Boernert, and H. Eggers. Resampling of data
betweenarbitrarygridsusingconvolutioninterpolation. IEEE Trans. Med. Imaging,
18:385–392, 1999.
[83] S. B. Reeder, B. J. Wintersperger, O. Dietrich, T. Lanz, A. Greiser, M. F. Reiser,
G. M. Glazer, and S. O. Schoenberg. Practical approaches to the evaluation of
signal-to-noise ratio performance with parallel imaging: application with cardiac
imaging and 32-channel cardiac coil. Magn. Reson. Med., 54:748–754, 2005.
[84] S.J.Riederer, T.Tasciyan, F.Farzaneh, J.N.Lee, R.C.Wright, andR.J.Herfkens.
MR fluoroscopy: technical feasibility. Magn. Reson. Med., 8(1):1–15, 1988.
[85] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal
algorithms. Physica D, 60:259–268, 1992.
[86] N. K. Sabharwal and A. Lahiri. Role of myocardial perfusion imaging for risk
stratification in suspected or known coronary artery disease. Heart, 89:1291–1297,
2003.
[87] H. Schomberg and J. Timmer. The gridding method for image reconstruction by
Fourier transformation. IEEE Trans. Med. Imaging, 14:596–607, 1995.
[88] J. Schwitter. Myocardial perfusion. J. Magn. Reson. Imaging, 24:953–963, 2006.
[89] J. Schwitter, D. Nanz, S. Kneifel, K. Bertschinger, M. Buchi, P. R. Knusel, B. Mar-
incek, T. F. Luscher, and G. K. von Schulthess. Assessment of myocardial perfusion
in coronary artery disease by magnetic resonance: a comparison with positron emis-
sion tomography and coronary angiography. Circulation, 103:2230–2235, 2001.
[90] H. Sedarat, A. B. Kerr, J. M. Pauly, and D. G. Nishimura. Partial-FOV reconstruc-
tion in dynamic spiral imaging. Magn. Reson. Med., 43:429–439, 2000.
[91] P. Sipola, K. Lauerma, M. Husso-Saastamoinen, J. Kuikka, E. Vanninen, T. Laiti-
nen, H. Manninen, P. Niemi, K. Peuhkurinen, P. Jaaskelainen, M. Laakso, J. Kuu-
sisto, and H. Aronen. First-pass mr imaging in the assessment of perfusion impair-
ment in patients with hypertrophic cardiomyopathy and the asp175asn mutation of
the alpha-tropomyosin gene. Radiology, 226:129–137, 2003.
[92] D. K. Sodickson and W. J. Manning. Simultaneous acquisition of spatial harmon-
ics (SMASH): Fast imaging with radiofrequency coil arrays. Magn. Reson. Med.,
38:591–603, 1997.
[93] R. S. Staewen, A. J. Johnson, B. D. Ross, T. Parrish, H. Merkle, and M. Garwood.
3-D FLASH imaging using a single surface coil and a new adiabatic pulse, BIR-4.
Invest. Radiol., 25(5):559–567, May 1990.
115
[94] K. Sung and K. S. Nayak. The design and use of tailored hard-pulse trains for
uniform saturation of myocardium at 3 tesla. Magn. Reson. Med., 60:997–1002,
2008.
[95] J. Trzasko and A. Manduca. Highly undersampled magnetic resonance image recon-
struction via homotopic l
0
-minimization. IEEE Trans. Med. Imaging, 28:106–121,
2009.
[96] J. Tsao. On the UNFOLD method. Magn. Reson. Med., 47:202–207, 2002.
[97] J. Tsao, P. Boesiger, and K. P. Pruessmann. k-t BLAST and k-t SENSE: dynamic
MRI with high frame rate exploiting spatiotemporal correlations. Magn. Reson.
Med., 50:1031–1042, 2003.
[98] J. A. Utz, R. J. Herfkens, and J. A. Heinsimer et al. Cine MR determination of left
ventricular ejection fraction. AJR, 148:839–843, 1987.
[99] J. J. van Vaals, M. E. Brummer, W. T. Dixon, H. H. Tuithof, H. Engels, R. C.
Nelson, B. M. Gerety, J. L. Chezmar, and J. A. den Boer. Keyhole method for
accelerating imaging of contrast agent uptake. J. Magn. Reson. Imaging, 3:671–
675, 1993.
[100] F. T. A. W. Wajer, R. Lethmate, J. A. C. van Osch, D. Graveron-Demilly, M. Fud-
erer, and D. van Ormondt. Interpolation from arbitrary to cartesian sample posi-
tions: Gridding. In ProRISC/IEEE Workshop, pages 571–577, 2000.
[101] A.G.Webb,Z.-P.Liang, R.L.Magin, andP.C.Lauterbur. Applicationsofreduced-
encoding MR imaging with generalized-series reconstruction (RIGR). J. Magn.
Reson. Imaging, 3:925–928, 1993.
[102] W. M Weiger, K. P. Pruessmann, and P. Boesiger. 2D SENSE for faster 3D MRI.
Magma, 14:10–19, 2002.
[103] A. M. Weissler, W. S. Harris, and C. D. Schoenfeld. Systolic time intervals in heart
failure in man. Circulation, 37:149–159, 1968.
[104] S. D. Wolff, J. Schwitter, and R. Coulden et al. Myocardial first-pass perfusion mag-
netic resonance imaging a multicenter dose-ranging study. Circulation, 110(6):732–
737, 2004.
[105] Q. S. Xiang and R. M. Henkelman. K-space description for MR imaging of dynamic
objects. Magn. Reson. Med., 29:422–428, 1993.
[106] P. C. Yang, A. B. Kerr, A. C. Liu, D. H. Liang, C. J. Hardy, C. H. Meyer, A. Ma-
covski, J. M. Pauly, and B. S. Hu. New real-time interactive magnetic resonance
imaging complements echocardiography. J. Amer. Coll. Cardiol., 32:2049–2056,
1998.
[107] L. Ying and J. Sheng. Joint image reconstruction and sensitivity estimation in
SENSE (JSENSE). Magn. Reson. Med., 42:952–962, 1999.
116
[108] H.Yokota, S.Heidary, C.K.Katikireddy, P.Nguyen, P.M.Pauly, M.V.McConnell,
and P. C. Yang. Quantitative characterization of myocardial infarction by cardio-
vascular magnetic resonance predicts future cardiovascular events in patients with
ischemic cardiomyopathy. J. Cardiovasc. Magn. Reson., 10:17, 2008.
117
Appendix A :
Modulation of Spiral Trajectories
In this appendix the modulation property of spiral trajectories is proven using concentric
rings as an approximation model. The sampling function for the concentric rings with a
radial spacing ¢
r
in k-space can be written as
S(k
r
)=
M
R
X
n=1
±(k
r
¡n¢
r
) (A.1)
where M
R
is the number of rings. Since the sampling function is circularly symmetric
its Fourier transform, i.e., PSF is also circularly symmetric and can be written as a
functionofsinglevariableinradialdirectionr. Fromthebasictransformpair±(k
r
¡a)$
2¼aJ
0
(2¼ar) [8], the PSF P(r) can be obtained as,
P (r)=F
¡1
fS(k
r
)g=2¼¢
r
M
R
X
m=1
mJ
0
(2¼rm¢
r
) (A.2)
118
0 5 10 15 20 25
−0.5
0
0.5
1
r
Amplitude
PSF of original rings
PSF of shifted rings
Figure A.1: PSF profiles of the original concentric rings and the shifted rings by half of
radial spacing. Two profiles have the same magnitudes over whole re gion. Phases of two
profiles are alike at mainlobe and the second sidelobe but diffrent from each other at the
first sidelobe by ¼ rad.
whereJ
0
(r) is Bessel function of the1
st
kind of zero order. When the set of rings is shifted
in the radial direction by half the radial spacing, its PSF will be modified into,
P
shift
(r)=2¼¢
r
M
R
X
m=1
(m¡0:5)J
0
(2¼r(m¡0:5)¢
r
) (A.3)
The PSF profile of this shifted ring samples together with the original one is shown in
Fig. A.1 where the number of rings M
R
= 10 and the radial spacing ¢
r
= 0:1 are used.
The phase difference is exactly ¼ rad around the first sidelobe while two profiles are in
phase around the mainlobe. Therefore when the shifted samples are used together with
the original ones alternatively the aliasing stemming from the sidelobe will be modulated
to Nyquist frequency. The sidelobes keep occurring in the radial direction with smaller
119
magnitudes and in general two profiles are in and out of phase at the even and odd order
sidelobes respectively but in practice the second and higher order sidelobes are likely to
be out of FOV and deserve little attention.
The desired property regarding the phase of PSF can be also shown in an analytic way.
Generally the asymptotic approximation of the Bessel function of k
th
order is represented
by [1],
J
k
(t)¼
8
>
>
>
<
>
>
>
:
1
¡(k+1)
µ
t
2
¶
k
for t¿1 (A.4a)
r
2
¼t
cos
µ
t¡
k¼
2
¡
¼
4
¶
for tÀ1 (A.4b)
where ¡(k) is Gamma function. The behavior at mainlobe, i.e., the same magnitude
in both P(r) and P
shift
(r) can be easily seen using (A.4a), since for zeroth order Bessel
function (k =0),
P(r) = P
shift
(r)
¼ 2¼¢
r
M
R
1
¡(1)
near the origin (r¿1) (A.5)
Using (A.4b), P(r) and P
shift
(r) far from the origin can be expressed as follow,
P(r) =
M
R
X
m=1
m
s
2
¼(2¼rm¢
r
)
cos
³
2¼rm¢
r
¡
¼
4
´
=
1
¼
M
R
X
m=1
r
m
x
cos
³
2¼mx¡
¼
4
´
(A.6)
120
P
shift
(r) =
M
R
X
m=1
m
s
2
¼(2¼r¢
r
)(m¡0:5)
cos
³
2¼r(m¡0:5)¢
r
¡
¼
4
´
=
1
¼
M
R
X
m=1
s
m
2
(m¡0:5)x
cos
³
2¼mx¡¼x¡
¼
4
´
=
8
>
>
<
>
>
:
1
¼
P
M
R
m=1
q
m
2
(m¡0:5)x
cos
¡
2¼mx¡
¼
4
¢
for even x
¡
1
¼
P
M
R
m=1
q
m
2
(m¡0:5)x
cos
¡
2¼mx¡
¼
4
¢
for odd x
(A.7)
where x =r¢
r
is a dimensionless variable such that the PSF shows peaks at its interger
values. Comparing (A.6) and (A.7) we can see that for a large value of M
R
, P(r) and
P
shift
(r) have almost same amplitude with the same phase at even order sidelobes and
opposite phase at odd order sidelobes.
121

Abstract (if available)

Abstract Magnetic resonance imaging (MRI) is one of the most powerful non-invasive imaging modalities whose strengths include flexible tissue contrast, sensitivity to flow, and the lack of any harmful radiation. Ischemic heart disease (IHD) is a leading cause of death in the western world, and as such, MR approaches for the assessment of IHD have gained considerable attention. This dissertation will describe two improvements in cardiac MRI methodology that are aimed at improving the assessment of IHD: (i) Real-time cardiac imaging and (ii) First-pass myocardial perfusion imaging.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Assessment of myocardial blood flow in humans using arterial spin labeled MRI

PDF

Radio-frequency non-uniformity in cardiac magnetic resonance imaging

PDF

Improving the sensitivity and spatial coverage of cardiac arterial spin labeling for assessment of coronary artery disease

PDF

Improved myocardial arterial spin labeled perfusion imaging

PDF

Improving sensitivity and spatial coverage of myocardial arterial spin labeling

PDF

Fast upper airway MRI of speech

PDF

Radiofrequency pulse performance for myocardial ASL

PDF

Velocity-encoded magnetic resonance imaging: acquisition, reconstruction and applications

PDF

Fast flexible dynamic three-dimensional magnetic resonance imaging

PDF

Dynamic functional magnetic resonance imaging to localize activated neurons

PDF

Estimating liver iron non-invasively with high-field MRI

PDF

Visualizing and modeling vocal production dynamics

PDF

Characterization of lenticulostriate arteries using high-resolution black blood MRI as an early imaging biomarker for vascular cognitive impairment and dementia

PDF

High-dimensional magnetic resonance imaging of microstructure

PDF

Functional real-time MRI of the upper airway

PDF

Modeling of cardiovascular autonomic control in sickle cell disease

PDF

Quantitative MRI for oxygenation imaging in cerebrovascular diseases

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Theory and simulation of diffusion magnetic resonance imaging on brain's white matter

PDF

Shift-invariant autoregressive reconstruction for MRI

Asset Metadata

Creator Shin, Taehoon (author)

Core Title Dynamic cardiovascular magnetic resonance imaging for improved assessment of ischemic heart disease

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Publication Date 08/27/2009

Defense Date 07/16/2009

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag cardiac MRI,ischemic heart disease,magnetic resonance imaging,myocardial perfusion imaging,OAI-PMH Harvest,real time imaging

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Nayak, Krishna S. (committee chair), Mendel, Jerry (committee member), Pohost, Gerald M. (committee member), Valero-Cuevas, Francisco J. (committee member)

Creator Email shinage@gmail.com,taehoons@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m2578

Unique identifier UC1194419

Identifier etd-Shin-3045 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-254626 (legacy record id),usctheses-m2578 (legacy record id)

Legacy Identifier etd-Shin-3045.pdf

Dmrecord 254626

Document Type Dissertation

Rights Shin, Taehoon

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu