Optimal Observation Selection in Magnetic Resonance Spectroscopic
Imaging
Yun Gao
Certificate of Approval:
Thomas S. Denney, Jr.
Associate Professor
Electrical and Computer Engineering
Stanley J. Reeves, Chair
Associate Professor
Electrical and Computer Engineering
Jitendra K. Tugnait
Professor
Electrical and Computer Engineering
Amnon J. Meir
Associate Professor
Mathematics Department
John F. Pritchett
Dean
Graduate School
Optimal Observation Selection in Magnetic Resonance Spectroscopic
Imaging
Yun Gao
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
August 14, 2000
Vita
Yun Gao, daughter of Xin Gao and Shuhua Song, was born in Henan province, P.
R. China on December 2, 1965. She received her Bachelor degree and Master degree
in Powering Engineering from Wuhan University of Hydraulic and Electric Engineering
in July, 1985, and July, 1988, respectively. She worked in China Institute of Water
Resources and Hydropower Research from July 1988 to August 1996. She entered the
Ph.D. program in the Department of Electrical and Computer Engineering in September
1996. Her research interests are in the areas of digital signal/image processing, medical
imaging, optimal acquisition/image compression, image reconstruction, and multiband
digital signal reception.
iii
Dissertation Abstract
Optimal Observation Selection in Magnetic Resonance Spectroscopic
Imaging
Yun Gao
Doctor of Philosophy, August 14, 2000
(M.S., Wuhan University of Hydrolic and Electric Engineering, P. R. China, July 1988)
(B.S., Wuhan University of Hydrolic and Electric Engineering, P. R. China, July 1985)
133 Typed Pages
Directed by Stanley J. Reeves
Magnetic resonance spectroscopic imaging (MRSI) is a completely non-invasive
method for obtaining quantitative information regarding biochemical parameters [1].
It shows great promise for use in basic physiological research and for clinical imaging of
metabolic function [2,3]. The information obtained with MRSI can be used to assess
regional metabolic abnormalities in various pathologies. However, MRSI requires a great
deal of time to gather the data necessary to achieve satisfactory resolution. When prior
information about the image is available, it may be possible to reconstruct the image
from a subset of k-space samples. Therefore, we desire to choose the best possible combi-
nation of a small number of k-space samples to guarantee the quality of the reconstructed
image by using the available prior knowledge. In this thesis, we assume prior knowledge
only of the region of support (ROS) of the spatial-domain image.
Sequential forward selection (SFS) is appealing as an optimization method because
the previously selected sample can be observed while the next sample is selected. We
develop an e?cient computational strategy for this algorithm that allows SFS to be
iv
applied to this problem when the image region has a known region of support (ROS).
The combined algorithm e?ciently selects a reduced set of k-space samples from which
the ROS can be reconstructed with minimal noise amplification. Furthermore, if there
is no noise, the minimum density can be reached with this algorithm.
Hexagonal sampling gives a 13.4% sampling density reduction compared to rectan-
gular sampling for images with a circular ROS. However, nonuniform sampling patterns
are more e?cient than hexagonal sampling for the same ROS. To reduce selection time
and achieve higher resolution, we develop a sequential backward selection (SBS) algo-
rithm from samples on a hexagonal grid. Simulation results show that more e?ciency
and reduced selection time can be achieved with the proposed method in comparison
withSBSonarectangulargrid.
We develop two e?cient algorithms for optimizing the dithering pattern so that an
image can be reconstructed as reliably as possible from a periodic nonuniform set of
samples, which can be obtained from a dithered rectangular-grid array.
One algorithm is SBS of sample arrays. Taking into account the ROS of the im-
age, we sequentially eliminate the least informative array recursively until the minimal
number of arrays remain. In this scheme, we provide an e?cient update formula for the
criterion based on convolution kernels rather than large, non-sparse matrices, an e?cient
update formula for the convolution kernels based on the deleted array, and an e?cient re-
construction method based on convolutions. The proposed method dramatically reduces
storage and computational complexity.
The other algorithm is SFS of sample arrays. Based on the ROS image, we sequen-
tially select the array that minimizes the noise amplification recursively until the desired
v
number of arrays are selected. To avoid the singularity of the criterion when extending
the selection procedure to more samples than unknowns, we propose a modified criterion
for the case when the number of unknowns is more than the number of selected samples
and the complementary case. We also propose an e?cient method to update the crite-
rion based only on the deleted array in the previousstep to greatly reduce computational
time and avoid the inversion of a huge matrix. This method has great practical potential
because it can finish the selection process within half a minute for practical sizes.
The proposed schemes in this dissertation e?ciently optimize the MRSI observa-
tion in di?erent ways. In general, they will reduce observation time and overcome the
problems in various available optimization methods.
vi
Acknowledgments
I would like to especially thank my research advisor, Dr. Stanley J. Reeves, for his
generous guidance during my Ph.D. study. I would not have finished this dissertation
without his valuable ideas on research topics and generous help in real life. I would also
like to especially thank Dr. Kai-Hsiung Chang for his generous help during my graduate
studies. I would like to thank Dr. Thomas S. Denney for being my committee member
and giving me opportunities to use his machine. I would like to thank Dr. Jitendra K.
Tugnait and Dr. Amnon J. Meir for being members of my committee. I would like to
thank Dr. Donald B. Twieg and Dr. Brad Newcomer at the University of Alabama at
Birmingham for providing me real MR data. I would like to thank Charles W. Hayes for
sketching a preliminary version of the proof in appendix A. I would also like to thank
Dr. Zhi Ding at the University of Iowa for his suggestions.
I would like to thank the Department of Electrical and Computer Engineering,
my research advisor, and the Whitaker Foundation for providing financial support. I
would like to thank Dr. Nikolaos V. Tsekos and the cardiovascular imaging center in the
Washington University School of Medicine for their patience with me.
I would like to especially thank Debbie Reeves, Hsiu-zhu Chang, and Cynthia Rouse
for their generous help in my di?culties. I would also like to thank Betty J. Kelly, Jo
Ann Loden, Shelia M. Collis, and Jo Ann Smith for their assistance.
Finally, I would like to especially thank my beloved parents, Xin Gao and Shuhua
Song, my beloved husband, Heping Jiao, my cherished sons, Gaoliang Jiao and David
vii
A. Jiao, as well as my brothers and sisters whose support and love have made this work
possible. This dissertation is dedicated to them.
viii
Style manual or journal used Discrete Mathematics (together with the style known
as ?auphd?). Bibliograpy follows van Leunen?s A Handbook for Scholars.
Computer software used The document preparation package T
E
X (specifically
L
A
T
E
X2
?
) together with the departmental style-file auphd.sty.
Table of Contents
List of Figures xi
List of Tables xiv
1 Introduction 1
1.1 MRI,MRS,andMRSI............................. 3
1.1.1 MRI................................... 3
1.1.2 MRS................................... 8
1.1.3 MRSI .................................. 10
1.2 MRSIAcquisitionMethodsandExistingProblems............. 12
1.2.1 RectangularSampling......................... 12
1.2.2 MultiSpinEchoImaging ....................... 13
1.2.3 Hexagonal Sampling . . . . ...................... 16
1.2.4 OtherExistingObservationOptimizationStrategies ........ 16
1.3 ThesisOverviewandContributions...................... 19
1.3.1 Overview ................................ 19
1.3.2 Contributions.............................. 20
2 Efficient Sequential Forward Selection 23
2.1 Introduction................................... 23
2.2 Criteria ..................................... 24
2.2.1 SignalModel .............................. 24
2.2.2 StandardCriterion........................... 24
2.2.3 ModifiedCriterion........................... 25
2.3 SFSAlgorithm ................................. 27
2.3.1 StandardSFSalgorithm........................ 27
2.3.2 ModifiedSFSalgorithm........................ 28
2.4 E?cientComputationForSFSAlgorithm.................. 30
2.5 ExtendingtheSFSProcess .......................... 33
2.6 Experiments................................... 36
2.7 Discussion.................................... 51
3 Efficient Backward Selection of Hexagonal Samples in k-space 55
3.1 Introduction................................... 55
3.2 Selectionalgorithm............................... 56
3.3 Regular Hexagonal Sampling . . . ...................... 58
3.4 GeneralizedMultidimensionalDiscreteFourierTransform ......... 62
ix
3.5 Experiments................................... 65
3.6 Conclusion ................................... 68
4 Sequential Backward Array Selection 70
4.1 Introduction................................... 70
4.2 TheMathematicsofSampling......................... 70
4.3 BasicOptimizationAlgorithm......................... 74
4.4 Imagereconstructionfromnonuniformsamples............... 79
4.5 Experiments................................... 80
4.6 Conclusion ................................... 83
5 Sequential Forward Array Selection 85
5.1 Introduction................................... 85
5.2 Notation..................................... 85
5.3 Criterion..................................... 87
5.4 E?cientrecursivealgorithm.......................... 88
5.5 E?cientcomputationalmethod ....................... 92
5.6 Experiments................................... 94
5.7 Conclusion ................................... 99
6 Conclusion and Discussion 102
6.1 E?cientsequentialforwardselection.....................102
6.2 E?cient backward selection of hexagonal samples . . ............103
6.3 E?cientselectionofsampleblock.......................103
6.4 E?cientsequentialforwardselectionofarray ................103
6.5 Futurework...................................104
Bibliography 110
x
List of Figures
1.1 (a) The spin magnetization M
0
without the application of the rf pulse,
(b) The spin magnetization M
0
with the rf field H
1
(t)applied. ...... 4
1.2 ImageweightingasafunctionofTRandTE ................ 6
1.3 The replication of elliptical ROS for rectangular sampling of k-space at
Nyquist density(a and b aresemimajorandminorrespectively). ..... 14
1.4 The replication of elliptical ROS for hexagonal sampling of k-space at
Nyquistdensity................................. 17
2.1 (a) The real image, (b) The ROS image (white indicates the support region) 37
2.2 (a) The criterion curve versus the number of selected samples with the
128x128 ROS image, (b) The flop curve versus selection step with the
128x128 ROS . . . . . ............................. 40
2.3 (a) 5297 optimized k-space sample distribution (the 5168 black dots repre-
sent the optimized minimum density, larger gray dots indicates 129 more
optimized samplesthan theminimumdensity), (b) Absolutedi?erence be-
tween the optimized image reconstruction and the fully sampled 128?128
image. (Dynamic range is scaled by a factor of 5 for greater visibility). . . 42
2.4 (a) The image from the minimum density k-space data(5168 points). (b)
The image from 5297 k-space samples (the ROS has been imposed as a
constraintinthereconstruction)........................ 43
2.5 (a) The point spread function from optimized k-space data (b) The point
spreadfunctionfrominverseFFT. ...................... 44
2.6 (a) Lowpass image from FFT, (b) Lowpass recursive solution at iteration
3, (c) Lowpass recursive solution at iteration 20, (d) Lowpass recursive
solution at iteration 100. . . . . . . ...................... 45
2.7 (a) Water image from the full 64 ? 64 real MRSI data, (b) Fat image
from the full 64?64 real MRSI data, (c) The 64?64 ROS image (white
indicates the support region), (d) 48?48 ROS image extracted from the
64?64 ROS (white indicates the support region), (e) 48?48 water image
extracted fromthe 64?64 real water image, (f) 48?48 fat image extracted
from the 64?64realfatimage......................... 47
xi
2.8 (a) 1513 optimized k-space data distribution with 64?64 FOV (the min-
imum optimized 1439 points indicated by black dots, larger gray dots
indicate 74 more samples than the minimum density), (b) 1513 optimized
k-space samples with a reduced 48 ? 48 FOV (the minimum optimized
1439 points indicated by black dots, larger gray dots indicate 74 more
samplesthantheminimumdensity)...................... 48
2.9 (a) Noise criterion curve versus the number of selected samples with the
48 ?48 ROS, (b) Flops per selection step versus the number of selected
samples with the 48?48ROS......................... 50
2.10 (a) OptimizedPSFinsidetheROSforareduced48?48 FOV, (b)Lowpass
PSF for 48?48 FOV with 1513 lowpass k-space samples, (c) Lowpass PSF
for 64?64 FOV with 1513 lowpass k-space samples. ............ 51
2.11 (a) Water image reconstructed from the 1513 optimized k-space samples
with 64 ? 64 FOV with ROS constraint imposed, (b) Fat image recon-
structed from the 1513 optimized k-space samples with 64?64 FOV with
ROSconstraintimposed. ........................... 52
2.12 (a) Absolute di?erence between the optimized water image reconstruction
and the fully sampled 64?64 image (b) Absolute di?erence between the
optimized fat image reconstruction and the fully sampled 64?64 image.
(Dynamic range is scaled by a factor of 5 for greater visibility). ..... 52
2.13 48?48 FOV reconstruction: (a) Water image reconstructed from the 1513
optimized k-space samples with ROS constraint imposed, (b) Fat image
reconstructed from the 1513 optimized k-space samples with ROS con-
straint imposed, (c) Water image from the 1513 lowpass k-space samples,
(d) Fat image from the 1513 lowpass k-space samples. . . . . . . ..... 53
2.14 (a) Absolute di?erence between the optimized water image reconstruc-
tion and the 48 ? 48 extracted image, (b) Absolute di?erence between
the optimized fat image reconstruction and the 48?48 extracted image.
(Dynamic range is scaled by a factor of 5 for greater visibility). ..... 54
3.1 (a) Hexagonal sampling pattern, (b) Hexagon . . . ............ 61
3.2 Smithformdecompositionflowchart ..................... 64
3.3 (a) 96 ?96 image, (b) 96?96 ROS, (c) Reconstructed image from 5550
rectangular k-space samples, (d) Reconstructed image from 5550 hexago-
nalk-spacesamples. .............................. 66
xii
3.4 Criterioncurvesversussamplesremaining................... 67
3.5 Timeversussamplesremaining......................... 67
4.1 typicalsamplingpattern. ........................... 71
4.2 Hexagonal sampling pattern. . . . ...................... 72
4.3 Periodicnonuniformsamplingexample. ................... 73
4.4 (a) The 64 ? 64 region of support (3205 unknowns) (b) 3584 selected
samples in the case 1 (c) 3328 selected samples in the case 2 (d) 3328
selected samples in the case 3 (black markers indicate selected samples) . 82
4.5 (a) The 128?128 region of support(white points indicate 3205 unknowns)
(b) 3328 selected samples in Case 4 ...................... 83
5.1 (a) The 64 ? 64 region of support (3205 unknowns) (b) 3328 selected
samples in Case 1 (c) 3328 selected samples in Case 2 (d) 3328 selected
samplesinCase3(blackmarkersindicateselectedsamples)........ 96
5.2 (a) The 128?128 region of support(white points indicate 3205 unknowns)
(b) 3328 selected samples in Case 4 (c) 4096 selected samples in Case 4. . 97
5.3 (a) The semilog criterion curve in Case 2 (b) The semilog criterion curve
inCase4 ....................................100
xiii
List of Tables
2.1 Comparisonofthecriterionandtimeversustol............... 39
2.2 Comparison of the selection criterion and time versus FOV with MRSI
ROSimage ................................... 46
4.1 Comparisonofthecriterioninthreecases.................. 81
5.1 Comparisonofthecriterioninfourcases................... 98
5.2 Comparisonofthecriterioninfourcases................... 98
xiv
Chapter 1
Introduction
Magnetic resonance spectroscopic imaging (MRSI) is an extension of localized mag-
netic resonance spectroscopy (MRS), which allows collection of in-vivo spectroscopic
information from multiple regions simultaneously. MRSI is a completely non-invasive
methodforobtainingquantitativebiochemicalinformation[1], whichshowsgreatpromise
for use in basic physiological research and for clinical imaging of metabolic function [2,3].
The information acquired with MRSI is shown as either spectral line plots at a particu-
lar location or images of metabolite distributions, which can be used to assess regional
metabolic abnormalities in various pathologies. In recent years, spectroscopic imaging
studies have been performed on a large number of pathologies [3]. 1H MRSI studies have
indicated both focal and global metabolic changes in a variety of diseases including brain
tumors [4?9], sub-acute and acute cerebral infarction [10,11], multiple sclerosis [12,13],
AIDS dementia [14,15], Alzheimer?s disease [16?18], and epilepsy [19?21]. Most of these
diseases present challenges to neuronal viability, which relate to a reduction in the NAA
concentration. For instance, a reduction in the NAA concentration was shown in severely
a?ected areas in stroke. Lactate is often increased in the same regions, indicating is-
chemia [22]. The observed regional reduction of NAA could mark onset of Alzheimer?s
disease in an early stage. Observed reductions in brain myo-inositol could indicate alter-
ations in liver biochemistry [18]. In addition, an increase of choline is observed in most
malignant brain tumors. For example, the choline increase has been found to correlate
with tumor grade in astrocytomas [11]. On the other hand, applications of phosphorus
1
2
MRSI have focused on energy metabolism in the human brain, in skeletal muscle and
cardiac muscle [23?25]. For example, metabolic changes have been observed in the brain
due to tumors [26] and epilepsy [21]. ATP depletion was found to indicate infarcted
regions in the cardiac muscle [25].
MRSI requires dedicated acquisition techniques because it combines elements from
MRS and magnetic resonance imaging (MRI) [1,3,27,28]. Despite its potential, MRSI
has remained largely in a research setting because a relatively long time is required to
acquire the images. In order forMRSIto reach its full potential in a clinical environment,
the acquisition time must be reduced. When prior knowledge of the image is available,
it may be possible to reconstruct the image from a subset of k-space samples. Examples
of such prior knowledge are region of support, locations of anatomical boundaries, and
relative smoothness of various image areas. Given this potential for reducing the number
of k-space samples, we desire to choose the best possible combination of these samples
to guarantee the quality of the reconstructed image. Towards this end, we are pursuing
strategies where the imaging process is optimized so that only the most significant por-
tions of the data set are acquired. The strategies we are developing have demonstrated
significant reductions in imaging time with minimal noise amplification and no reduction
in resolution.
3
1.1 MRI, MRS, and MRSI
1.1.1 MRI
MRI is a completely non-invasive method for imaging soft tissues and other struc-
tures in the body. MRI makes use of the nuclear magnetic resonance (NMR) phe-
nomenon, which was discovered in 1946 [29], to obtain a proton density image of a given
object. The history of MRI based on the principle of NMR goes back to the early 1970s.
At that time both Lauterbur [30] and Damadian [31] proposed that NMR spectroscopic
techniques could be applied to human imaging. Andrew et al. [32,33] demonstrated
that NMR imaging could show the details of the human anatomy. Subsequent images
obtained by Moore et al. [34] and Holland et al. [35] proved NMR tomography is indeed
capable of performing diagnostic imaging. Since MRI became clinically viable about 15
years ago, it has proven itself as an important diagnostic tool for a wide variety of health
problems. MRI is unique in that it exploits a magnetic moment called spin, which is a
physical property possessed by certain atomic nuclei such as those of hydrogen (protons).
It is completely noninvasive and appears to be safe for children and probably even for
pregnant women.
The Physics of MRI [1,29,36?41]
All materials contain nuclei, which are protons, neutrons, or a combination of both
[39]. Nuclei that have an odd number of protons, neutrons, or both in combination,
possess anuclear spin or a magnetic moment. Under normal circumstances, these atomic
nuclei in a sample are aligned randomly, and therefore their magnetic moments cancel
each other and no signal can be measured. However, if a sample is placed into a powerful
static magnetic field H
0
, spins align themselves with the external field either in a parallel
4
x x
y
y
Ho Ho
Mo
Mo
?
z z
H1
?=??
o
?
p
H1(t )dt 
( a ) ( b )
Figure 1.1: (a) The spin magnetization M
0
without the application of the rf pulse, (b)
The spin magnetization M
0
with the rf field H
1
(t) applied.
direction (in low-energy state) or in an antiparallel direction (in high-energy state) and
startto precessattheLarmorprecession frequency, whichisproportionalto themagnetic
field strength. At thermal equilibrium, the net spin magnetization vector M
0
is along
the external magnetic field H
0
. When an rf signal is applied to the system in the form
of an rf magnetic field H
1
, those protons in the lower-energy state tend to be excited to
the higher-energy state and as a result, M
0
can be rotated or flipped by any angle ?.
After cessation of the rf magnetic field H
1
, the excited protons tend to return to their
low-energy state by emitting a well-defined rf frequency (i.e. the same frequency as the
applied rf) by spin relaxation mechanisms such as spin-lattice and spin-spin relaxations,
and M
0
lines up with H
0
again. This emission of rf signals is then picked up by an rf
coil placed close to the excited object. This signal is called a free induction decay (FID),
which is the central part of NMR imaging. The spin magnetization in the rotating frame
with and withoutan rfpulseis shownin Figure 1.1. Theflipangle ? of the magnetization
M
0
is given by ? = ?
integraltext
t
p
o
H
1
(t)dt,whereH
1
(t) is the time-varying rf field intensity and
?
p
is the application length of the rf pulse. The angle ? is typically set to 90
o
or 180
o
.
5
The two relaxation mechanisms are the transverse or spin-spin relaxation known as
T
2
relaxation andthelongitudinalorspin-lattice relaxation knownasT
1
relaxation. Both
of these relaxations and their respective relaxation times (T
1
and T
2
)arequitesensitive
to the molecular structures and the environments surrounding the nuclei. For example,
the mean T
1
and T
2
values of normal tissues and those of many malignant tissues di?er
significantly, therebyallowing usto di?erentiate malignant tissue from normal tissues [1].
The imaging capabilities of these two important parameters, T
1
and T
2
, together with
the spin densities and the flow-dependent phase information of the objects, make NMR
imaging a unique, versatile, and powerful technique in diagnostic imaging.
Three widely used spin-echo techniques play an essential role in data acquisition
for NMR and NMR imaging: the HAHN spin-echo technique [42], the CPMG technique
[43,44], and the stimulated echo techniques [42,45]. Both Hanh and CPMG techniques
have wide applications in all phases of NMR and NMR imaging to overcome several
adverse e?ects that often arise in experimental situations, such as the e?ects of the
gradient pulse rise time and the field inhomogeneity. The same goal can be achieved
with stimulated echo but without su?ering additional T
2
decay.
In a typical spin echo sequence, the RF pulses are repeated many times; the interval
between the repeated pulses is called the repetition time (TR). The interval in the pulse
sequence between the end of the RF pulse and the detection of the returning MR signal
is called the echo time (TE). By adjusting propercombinations of TR and TE, the image
will show contrast related mainly to di?erences in proton density, T1 relaxation time, or
T2 relaxation time of the tissue ? in other words, the image is proton density weighted,
T1 weighted or T2 weighted. If long TR and long TE are chosen so that T1 di?erences
6
TE
TR
    Short
Long
Short
Long
T1-weighted
Spin density-
    weighted
T2-weighted
Figure 1.2: Image weighting as a function of TR and TE
are responsible for the observed contrast, the image is said to be T1-weighted. If short
TR and short TE are chosen so that T2 di?erences are responsible for the observed
contrast, the image is said to be T2-weighted. If long TR and short TE are selected so
that both the T1 and T2 e?ects are minimized, then the image is said to be spin density-
or proton density-weighted. Figure 1.2 shows these weighting rules.
In NMR imaging or MRI, a field gradient or a set of gradients is deliberately added
to resolve the spatial distribution of spins by Fourier encodings. The basic form of the
signal obtained can be expressed as [46?48]
s(t)=M
0
integraldisplayintegraldisplayintegraldisplay
?
??
?(x,y,z)exp{?i?
integraldisplay
t
0
[xG
x
(t
prime
)+yG
y
(t
prime
)+zG
z
(t
prime
)]dt
prime
}d
x
d
y
d
z
(1.1.1)
The resulted FID or echo is in e?ect a Fourier transform-domain representation of
the spatially distributed spin density. Many imaging algorithms can be derived from this
basic form of the 3-D equation.
From the basic imaging equation in (1.1.1), a 3-D discrete imaging equation can be
obtained bydiscrete Fourier encoding. Fourier-domain scanning usesboth frequencyand
phase encodings. Discrete Fourier encoding can be achieved by changing the gradient
7
amplitude (or its length):
S(g
x
,t
y
,g
z
)=M
0
integraldisplayintegraldisplayintegraldisplay
?(x,y,z)exp[?i?(xg
x
T
x
+yG
y
T
y
+zg
z
t
z
)]d
x
d
y
d
z
(1.1.2)
where g
x
= ntriangleG
x
and g
z
= ntriangleG
y
. triangleG
x
and triangleG
z
are the increments of the x-
gradientamplitudeandz-gradientamplitude,andndenotesthestepnumberofthephase
encodinggradient. Thepulselength can also bevaried providedthat thepulseamplitude
is kept constant. However, the former phase encoding method is more widely used due
to the fact that the constant time or pulse width T
x
employed in the former eliminates
NMR parameter-dependent signal variations such as T
2
decay [1]. The y-gradient G
y
is used as the readout gradient, which is also the frequency-encoding gradient. This is
distinguished from the phase-encoding gradients which in this case are the x-gradient
and z-gradient. The number of phase encoding steps determines the x-directional and
z-directional resolution, while the y-directional resolution remains constant regardless of
the number of encoding steps employed. Theoretically the resolution of the y-directional
resolution can be infinitely high. However, in most cases the resolution in the y-direction
is limited by the bandwidth of the NMR data, the sampling interval, the frequency
response of the electronics, and the SNR. In practical NMR imaging, the resolution
is therefore largely determined by the number of phase encoding steps employed (e.g.,
the x-gradient and the z-gradient in this case) provided that the signal has a su?cient
SNR.Other important physical parameters involved in NMR imaging arethe spin-lattice
and spin-spin relaxation times T
1
and T
2
. The increased phase encoding steps will
improve resolution but cause long imaging time. If we let w
x
= ?g
x
T
x
, w
y
= ?G
y
t
y
,and
8
w
z
= ?g
z
T
z
, (1.1.2) can be rewritten as
S(w
x
,w
y
,w
z
)=M
0
integraldisplayintegraldisplayintegraldisplay
?(x,y,z)exp[?i(w
x
x+w
y
y +w
z
z)]d
x
d
y
d
z
(1.1.3)
The 3-D density function can be obtained by taking a 3-D inverse Fourier transform of
S(w
x
,w
y
,w
z
):
?(x,y,z)=c
integraldisplayintegraldisplayintegraldisplay
S(w
x
,w
y
,w
z
)exp[?i(w
x
x+w
y
y +w
z
z)]d
w
x
d
w
y
d
w
z
(1.1.4)
where c is a constant. 2-D Fourier imaging can be derived directly by setting z = z
0
as
S(w
x
,w
y
)=M
0
integraldisplayintegraldisplay
?(x,y,z = z
0
)exp[?i(w
x
x+w
y
y)]d
x
d
y
(1.1.5)
and
?(x,y,z = z
0
)=c
integraldisplayintegraldisplay
S(w
x
,w
y
)exp[?i(w
x
x+w
y
y)]d
w
x
d
w
y
(1.1.6)
1.1.2 MRS
Because the magnetic field in the vicinity of a nucleus is influenced by the local
conditions, the Larmor frequencies for two di?erent nuclei of the same species may be
shifted apart. This phenomenon is also known as chemical shift. The Larmor frequencies
for two di?erent nuclei of the same species are given by
v
k
= ?(1??
k
)B
0
(1.1.7)
v
ref
= ?(1??
ref
)B
0
(1.1.8)
9
where ?
k
and ?
ref
are the screening (or shielding) constants of the sample and the
reference materials, respectively, due to the diamagnetic, paramagnetic, and interatomic
current e?ects. The chemical shift is defined as the di?erence of these two equations:
w
k
? v
k
?v
ref
= ?(?
ref
??
k
)B
0
(1.1.9)
To eliminate the field dependency of w
k
, a relative chemical shift ?
k
is defined as
?
k
=
?
ref
??
k
1??
ref
?10
6
ppm (1.1.10)
In proton NMR spectroscopy, tetramethylsilane (TMS) is a typical reference mate-
rial. A positive value of ?
k
indicates that the nuclei of the sample resonate at a frequency
higher than the reference material. In a high magnetic field the spectrum of a partic-
ular substance may be resolvable into several peaks due to the chemical shift, and the
frequency di?erence between the peaks becomes larger with increasing magnetic field
strength B
0
. Therefore, chemical shift plays a key role in NMR spectroscopy.
Magnetic resonance spectroscopy is a unique tool in that it can noninvasively obtain
quantitative biochemical information. Proton magnetic resonance spectroscopy plays an
important role for the study of in vivo metabolites in animals and humans [49?55]. Phos-
phorus spectroscopy plays a key role in NMR spectroscopy especially in in vivo human
studies [56,57]. However, most of the spectroscopy techniques critically depend upon
the quality of the spectral resolution that can be obtained. Good resolution can be
achieved by optimizing the magnetic field homogeneity over the sample volume of inter-
est by localization techniques. A number of volume localization techniques have been
10
developed such as rotating-frame zeugmatography [58,59], DRESS [60], SPARS [61,62],
STEAM [50,63] and other techniques using high-order gradient fields and multidimen-
sional localization [49,51?56,64?77].
1.1.3 MRSI
MRSI is a combination of MRS and MRI. MRSI di?ers from conventional spec-
troscopy in that the spectroscopic information is presented in the form of an image
instead of in the form of a graph. Consequently, MRSI gives results which are readily
interpretable in a qualitative sense and allows collection of in vivo spectroscopic data
from multiple regions simultaneously. It shows great promise for use in basic physiolog-
ical research and for clinical imaging of metabolic function [2]. Spectroscopic imaging
studies have been performed on a large number of pathologies [3]. 1H MRSI studies have
indicated bothfocaland globalmetabolic changes inanumberofdiseasesincludingbrain
tumors [4?9], sub-acute and acute cerebral infarction [10,11], multiple sclerosis [12,13],
AIDS dementia [14,15], Alzheimer?s disease [16?18], and epilepsy [20]. Most of these
diseases cause a reduction in the NAA concentration and challenges to neuronal vi-
ability [3]. Applications of phosphorus MRSI have focused on energy metabolism in
the human brain, skeletal muscle and cardiac muscle [23?25]. For instance, Metabolic
changes have been observed in the brain due to brain tumors [26] and epilepsy [21]. ATP
depletion was found in the cardiac muscle to indicate infarcted regions [25].
MRSI can be expressed in the form of four-dimensional Fourier NMR imaging [27],
which o?ered new possibilities of doing spatially resolved high-resolution spectroscopy.
11
This 4-D NMR technique is a variation of the three-dimensional Fourier imaging tech-
nique with an additional spectral dimension. The new spectral dimension indicates the
chemical shift information. The acquired data in the frequency domain is then trans-
formed into the spatial domain by an inverse four-dimensional Fourier transform to gen-
erate the NMR spectra as well as the 3-D image ? in other words, the NMR spectrum
at each voxel in the spatially resolved 3-D volume image. The observed NMR signal is
given by
S(w
x
,w
y
,w
z
,t
x
)=
integraldisplay
w
k
integraldisplay
z
integraldisplay
y
integraldisplay
x
?(x,y,z,w
k
)exp[?i(w
x
x+w
y
y+w
z
z+w
k
t
k
)]d
x
d
y
d
z
d
w
k
(1.1.11)
where w
x
= ?g
x
T
x
, w
y
= ?g
y
T
y
,andw
z
= ?g
z
T
z
, g
x
, g
y
,andg
z
are the amplitude
variables of the x, y,andz gradients, and T
x
, T
y
,andT
z
are constant time intervals. A
three-dimensionally resolved volume spectral image can be reduced to a two-dimensional
slice spectral image by converting any one of the three coding gradients into the selection
gradient.
The standard MRSI protocol is to acquire a set of k-space samples on a rectangular
grid and perform an inverse Fourier transform on the acquired data to obtain the spatial-
domain image and the NMR spectrum. Full four-dimensional Fourier imaging requires
an imaging time of N
x
N
y
N
z
T
R
seconds for spectrally resolved 3-D volume imaging. N
x
,
N
y
, N
z
and T
R
are the number of encoding steps of the x, y, z gradients and repetition
time, respectively. Then N
x
, N
y
,andN
z
will represent the number of image pixels in
the spatial domains x, y,andz. Quite a large number of coding steps will be required
even for the case of 2-D slice-selected spectroscopic imaging, and they are often too long
for human in vivo imaging. This limitation results in truncated k-space data, which
12
causes ringing and limits the resolution of the reconstructed image. However, increased
resolution is particularly important for making MRSI a clinically viable technique. For
example, theuseofMRSItolocalizeandassessbraintumors[7]andmultiplesclerosis[78]
as well as to determine the seizure focus of temporal lobe epilepsy [19] will all require
increased resolution. A number of image modeling techniques have been proposed to
incorporate a priori information from a proton-density MR imaging study to improve
the resolution of the MR spectroscopic image of the same object [79?81]. When these
image models are imposed on the image reconstruction process, the combination of k-
space samples acquired can have a considerable e?ect on the quality of the reconstructed
image. In this thesis, we consider images with a limited region of support (ROS). When
the image has a limited region of support (ROS), it is possible to reconstruct the image
from a subset of k-space samples. Therefore, we desire to develop strategies to choose the
best possible combination of a reduced set of k-space samples to guarantee the quality
of the reconstructed image while reducing the imaging time required. This is the goal of
the thesis.
1.2 MRSI Acquisition Methods and Existing Problems
1.2.1 Rectangular Sampling
Rectangular sampling acquires a set of k-space samples on a rectangular grid. The
spatial-domain image is reconstructed by performing an inverse Fourier transform on the
acquired data. This is the standard MRSI protocol. If the spatial extent of the image
is limited, rectangular sampling at the Nyquist rate will permit accurate reconstruction
with no spatial aliasing (overlapping of the periodic images). Rectangular sampling
13
may be formally defined as sampling on a regular lattice whose basis vectors form an
orthogonal set. The basis vectors need not be of equal magnitude. However, the time
required to gather the data while maintaining an adequate signal-to-noise ratio (S/N)
limitsthenumberofk-spacesamplesthatcanbeacquired, whichwehavediscussedinthe
previoussection. Thislimitation istypically addressedbylimitingthek-spaceacquisition
to low spatial-frequency samples, which causes ringing and limits the resolution of the
reconstructed image. Ontheother hand,intwo ormoredimensionsrectangular sampling
is ine?cient if the image is isotropically spatially limited. Consider rectangular sampling
ofk-space foranimage witha2-D circularorelliptical ROSat Nyquistdensity. Asshown
in Figure 1.3, gaps exist among the replications. In the 2-D case, the mean sampling
density is
4
?
as dense as the minimum sampling density required to avoid aliasing. This
translates into 27% more time to obtain an image than is necessary, which can be critical
in MRSI. Therefore, various other sampling schemes have been proposed to optimize the
sampling strategy.
1.2.2 Multi Spin Echo Imaging
Multi-spin echo imaging is a class of rapid imaging techniques that includes echo
planar imaging (EPI) and spiral scan echo planar imaging (SEPI).
MansfieldwasthefirsttoproposeEPI.Thismethod, inprinciple, allows ustoobtain
2-D image data by a single excitation [82] and can be extended to a 3-D form. First, a
slice of thickness trianglez is selected by applying G
z
together with a 90
o
selective rf pulse.
Then an oscillating G
x
and a small constant gradient G
y
are applied. The generated
14
y
x
a
b
Figure 1.3: The replication of elliptical ROS for rectangular sampling of k-space at
Nyquist density(a and b are semi major and minor respectively).
FID and echoes can be written as
f(t)={
?
summationdisplay
n=0
g
n
(t?2n?)exp(?
t
T
2
) (1.2.1)
where g
n
(t?2n?) is the echo signal:
g
n
(t?2n?)=s(t) (1.2.2)
The entire frequency or k-space can be covered with a single 90
o
rf pulse excitation
by adjusting G
x
and G
y
. However, there exist some problems in EPI. One problem is
that T
2
decay requires a rapidly alternating gradient pulse of large amplitude, which is
often di?cult to obtain in practice. Another problem in EPI is the poor resolution in the
15
phase encoding direction due to the limitation in the number of phase encoding steps or
the number of echoes that can be recalled due to the signal decay by T
2
.
Some of the di?culties in the original EPI such as the resolution limits in the phase
encoding direction can be partially resolved by the application of a spiral scan in k-space
known as SEPI [75]. The spiral scan represents a close approximation of the concentric
circles with which the entire frequency range of interest can be covered uniformly [83],
thereby achieving approximately circularly symmetric resolution. The required k-space
radial and angular sampling rates can be defined as
trianglek
r
=
2?
2N
r
(1.2.3)
trianglek
?
=
2?
N
?
(1.2.4)
where N
r
and triangler denote the number of rotations in the k-domain and the image resolu-
tion in the spatial domain, respectively. The total required image data acquisition time
T
D
is given by multiplication of angular incremental time triangleT, the number of rotations,
and the number of samples in one complete revolution.
T
D
= N
r
N
?
triangleT (1.2.5)
The main advantage of SEPI over standard EPI is the conjugate symmetry property
of the spiral scan. This property allows the sampling rate to be increased by a factor
of 2. One of the drawbacks of SEPI is that the complex 2-D interpolation of the polar
coordinatedatatoCartesian dataisnecessaryforusinganFFTtoreconstructthespatial
image to display. Another drawback is that the samples on the circles near the origin are
16
too close while the samples on the outermost circles are too wide apart, causing either
wasting of time or aliasing and increasing the complexity of the interpolation.
Multi-spin echo imaging techniques are inappropriate for standard MRSI, since the
standard MRSI protocol uses the readout direction for the spectroscopic information and
each point in k-space must be acquired by a separate excitation sequence. But they can
be adapted to MRSI with some modifications [84?88].
1.2.3 Hexagonal Sampling
Hexagonal sampling [89] allows an image with an elliptical ROS to be reconstructed
with no aliasing from fewer data points than rectangular sampling. Hexagonal sampling
at the Nyquist rate has a 13.4% e?ciency gain over rectangular sampling at the Nyquist
rate and yields the most densely packed spatial replications for space-limited signals
[90,91]. However, Nyquist density sampling cannot reach the minimum density for an
elliptical ROS image. Consider hexagonal sampling of k-space for an image with a 2-D
circular or elliptical ROS at Nyquist density. As shown in Figure 1.4, gaps exist among
the replications. Minimum density sampling yields samples linearly independent of each
other, which cannot be reached through this uniform sampling strategy.
1.2.4 Other Existing Observation Optimization Strategies
The problem of optimizing the choice of k-space samples has also been addressed
in various other ways. Von Kienlin et al. [92] introduced a criterion for optimizing the
k-space samples using the SLIM [79] model for the image. This model assumes that
the image is made up of homogeneous regions. The proposed criterion is based on the
17
b
a
  y
x
Figure1.4: Thereplication ofelliptical ROSforhexagonal samplingofk-space atNyquist
density
resulting point-spread function. However, their final criterion must be numerically opti-
mizedbyamultidimensionalnonlinearoptimization algorithm, andthecriterioncontains
multiple local minima. Furthermore, the criterion has only an indirect relation to image
quality. Cao and Levin [93] used an image database to constrain the reconstructed image
to a limited set of basis images and optimized the choice of k-space samples to minimize
mean square error (MSE) in the reconstruction. The problem with this approach is that
it may overly constrain the reconstruction when one is attempting to detect anomalies in
the image. Furthermore, the authors did not propose an e?cient optimization method.
Plevritis andMackovski [94] proposeda heuristick-spaceselection method basedon find-
ing the highest-energy k-space locations of a binarized scout image. While this method
may improve reconstruction results as compared to lowpass sampling, it has no clear
connection to any measure of optimality in the reconstructed image.
18
Recognizing the combinatoric problem inherent in the previous approaches, two
independent groups have proposed an alternative sampling scheme called coset sampling
[95?97]. In this scheme, k-space is divided into rectangular blocks and the same irregular
pattern of samples is selected within each block so that the sampling scheme becomes
block-periodic. This approach is quite e?ective and reduces the optimization problem
considerably under certain conditions. However, to achieve a significant reduction in
the sampling density, the ROS must be well approximated by a collection of rectangles.
Reeves et al. proposedsequential backwardselection (SBS)[98,99] andsequentialforward
selection (SFS) [100] based on minimizing MSE in the reconstructed image. However,
SBS cannot be implemented in real time because SBS eliminates the worst candidate
rather than selects the best one. Therefore, observation cannot begin until the whole
selection process finishes. The SFS algorithm based on MSE does not apply when there
are more unknowns than selected samples, because the criterion is undefined in this case.
Thus, the SFS algorithm cannot begin with no samples selected.
In addition to optimizing MRSI data acquisition, the techniques proposed here can
be used for 3-D MRI. In the case of 3-D MRI the readout direction is used to sample
alongalineinthethirdk-spacedimension,sotheselineswillappearasindividualk-space
points in the first two k-space dimensions. Therefore, if the 3-D image has a limited ROS
in every slice along a particular dimension, the proposed algorithm can be used without
modification to select the k-space distribution. The most e?cient acquisition would
use echo-volume imaging to traverse several selected k-space lines for every excitation
period. The only restriction is that in the optimization algorithm we must use an ROS
that completely contains the ROS at each slice. However, a relative shift of the ROS
19
from slice to slice will not a?ect the optimal set of k-space samples for that slice, since
the criteria in (2) and (8) are una?ected by shifts of the ROS. We have proved this
assertion in the appendix.
1.3 Thesis Overview and Contributions
1.3.1 Overview
Nyquist density sampling yields the most densely packed spatial replications for
space-limited signals, while minimum density sampling yields samples linearly indepen-
dentofeach other butnotnecessarily uniformlydistributed. Itcan beshownthat for1-D
space-limited functions, the two densities (or rates) are identical. For higher-dimensional
space-limited functions, Nyquist densities can besignificantly higher than minimum den-
sities [101]. This is the rationale for optimized k-space sampling with a limited ROS. If
the ROS is smaller than the image, Nyquist sampling represents a redundant sampling
of k-space. The multidimensional extension of Papoulis? generalized sampling expansion
states that the minimum frequency-sampling density of an M-D space-limited function
is equal to the area of the ROS of the function [101]. That is, the minimum number of
samples in k-space equals the area of the ROS. However, this theorem does not tell us
how to achieve this minimum density sampling without losing any signal information.
In thisthesis, we develop four MRSIobservation optimization strategies to overcome
the existing problems in di?erent ways: 1) e?cient SFS [102?104], 2) e?cient SBS ap-
plied to hexagonal samples [105], 3) e?cient SBS for block-periodic sampling [106?108],
and 4) e?cient SFS for block-periodic sampling. In each case, we develop an e?cient op-
timization criterion based on minimizing the sum of squared errors in the reconstructed
20
image samples, an e?cient computational strategy for the optimization procedure, and
in one case an e?cient image reconstruction method. These strategies tend to minimize
the error criterion, but they do not guarantee the optimal observation selection.
1.3.2 Contributions
The contributions we propose are summarized as follows:
? A modified SFS algorithm [102?104] is developed to overcome the problems with
standard SFS [100]. This entails the following contributions:
? We present a modified form of the criterion that overcomes the existing prob-
lem for the standard SFS algorithm and develop an SFS algorithm for the
new criterion.
? We develop an e?cient computational strategy for the new algorithm as well
as the standard SFS algorithm, which drastically reduces storage space and
computational time.
? The combined algorithm can e?ciently select a reduced set of k-space samples
from which we can reach the minimum density without losing any signal
information. Furthermore, if the data is noisy as is the case in MRSI, we can
reduce noise amplification considerably by selecting a few more samples than
the number of unknowns in the ROS. Simulations with real MR data show
that a good-quality image can be reconstructed with only a few more samples
than the minimum density with this method.
21
? It is more applicable in practice than SBS in that the best sample can be
observed as soon as it is selected, making possible real-time selection and
acquisition.
? We present an e?cient SBS algorithm from hexagonal samples when the ROS is
approximated by an ellipse. In this method, we make the following contributions:
? A general periodic matrix FFT is implemented based on a rectangular FFT
with the Smith form decomposition. A structured Smith form decomposition
is implemented.
? By choosing an appropriate periodic matrix N, the hexagonal samples in the
frequency domain can be transformed to rectangular samples in the spatial
domain without interpolation directly with the generalized FFT. Therefore,
artifacts and extra computational time are avoided.
? An e?cient SBS algorithm is applied to hexagonal samples. The combined
algorithm e?ciently optimizes a reduced set of k-space samples and saves
a great deal of selection time in comparison with SBS applied to common
rectangular samples. The savings results from the fact that the algorithm
begins with fewer samples than with a rectangular grid.
? We develop an SBS algorithm for periodic sample blocks. In this scheme, we make
the following contributions:
? It provides an e?cient update formula for the criterion based on convolution
kernels rather than large, non-sparse matrices;
22
? It provides an e?cient update formula for the convolution kernels based on
the deleted array.
? It provides an e?cient reconstruction method based on convolutions rather
than matrix inverses. The proposed method dramatically reduces storage and
computational complexity.
? We develop an e?cient SFS algorithm for periodic sample blocks. In this scheme,
we make the following contributions:
? It provides a modified criterion to avoid the singularity problem when extend-
ing the selection process.
? It provides an e?cient update formula for the criterion based only on the
previous selected array.
? It decomposes the block sample algorithm into a combination of single sample
algorithms to e?ciently solve the complicated array selection problem.
? It seamlessly integrates the modified SFS algorithm for fewer samples than
unknowns with the standard algorithm for more samples than unknowns.
? It has great practical potential in that it can finish the selection in half a
minute for a practical size image.
Chapter 2
Efficient Sequential Forward Selection
2.1 Introduction
Sequential forward selection (SFS) is appealing as an optimization method because
the previously selected sample can be observed while the next sample is selected. SFS
tends to minimize the error criterion, but it does not guarantee the optimal observation
selection. When the number of selected k-space samples is less than the number of un-
knowns at the beginning of the selection process, the optimality criterion is undefined
and the resulting SFS algorithm cannot be used. In this chapter, we present a modified
form of the criterion that overcomes this problem and develop an SFS algorithm for the
new criterion. Then we develop an e?cient computational strategy for this algorithm
as well as for the standard SFS algorithm. The combined algorithm e?ciently selects a
reduced set of k-space samples from which the ROS can be reconstructed with minimal
noise amplification. With the proposed algorithm, we can reach the minimum density
without losing any signal information. If the samples are noisy as are the case in MRSI,
we can reduce noise amplification considerably with a few more samples than the number
of unknowns in the ROS. Furthermore, the proposed algorithm brings real-time imple-
mentation closer to reality because the previously selected sample can be observed while
the next sample is selected.
23
24
2.2 Criteria
2.2.1 Signal Model
We model the observed k-space samples y as a vector resulting from a linear trans-
formation of the original image x in the presence of additive noise w;thatis
y = Ax+w (2.2.1)
where w is zero-mean, independently identical distributed (i.i.d.) noise and A ? C
m?n
,
m and n are integers. Given the observed signal y, the goal is to reconstruct a good
estimate of x. We want to select the limited set of observations that will yield the
best possible reconstruction of x by using the known mapping from the original to the
observations. This is equivalent to choosing rows for A that correspond to the best set
of observations y. For MRI and MRSI, A represents a Fourier transform matrix with
columns removed corresponding to the voxels outside the ROS.
2.2.2 Standard Criterion
If the noise is zero-mean, i.i.d. and the reconstruction of x is performed via least
squares, we can define a criterion E(A) that is proportional to the mean square error in
the reconstruction [109] as
E(A)=tr(A
H
A)
?1
(2.2.2)
Obviously, the criterion is not defined when A has fewer rows than columns. A
critical limitation of the SFS algorithm based on this criterion is that it only applies
25
when the initial matrix has a column rank of n. A necessary condition, therefore, is that
A have at least n rows to begin the update process [100]. Consequently, this criterion
can not be used directly for an SFS algorithm that starts with no selected observations.
Therefore, in the next section we present a modified form of the selection criterion
that uses a pseudo-inverse in place of the inverse to make the initial selection process
possible. We optimize the selection process based on the modified criterion until A is
square. Afterwards, we switch the selection process to the standard algorithm.
2.2.3 Modified Criterion
When A has fewer rows than columns, we propose to use the following modified
criterion:
E(A)=tr(A
H
A)
+
(2.2.3)
where (A
H
A)
+
is the pseudo-inverse of A
H
A.LetU
H
AV = ? be the singular value
decomposition (SVD) of A ? C
m?n
with r = rank(A). Then (A
H
A)
+
= V(?
H
?)
+
V
H
,
where
(?
H
?)
+
= diag(
1
?
2
1
,...,
1
?
2
r
,0,...,0), (2.2.4)
and ?
i
is the ith diagonal entry of ?. Now, we consider the validity of the modified
criterion and its relation to the standard criterion. The minimum-norm least-squares
solution is given by
x
LS
=(A
H
A)
+
A
H
y
=
r
summationdisplay
i=1
?
i
?
H
i
x+(A
H
A)
+
A
H
w (2.2.5)
26
and
x = VV
H
x =
n
summationdisplay
i=1
?
i
?
H
i
x (2.2.6)
where x is the ideal solution with A full rank and no noise and ?
i
is the ith column
vector of V. Then,
E||x?x
LS
||
2
= ||
n
summationdisplay
i=r+1
?
i
?
H
i
x||
2
+?
2
w
tr(A(A
H
A)
+
(A
H
A)
+
A
H
)
=
n
summationdisplay
i=r+1
||x
T
?
i
||
2
+?
2
w
tr(A
H
A)
+
(2.2.7)
When A is underdetermined, there are two parts in the error. One part is due to
the loss of signal components; the other is due to noise. The error due to the loss of
these signal components is not reflected in the modified criterion. However, the criterion
does reflect the noise error. Thus, the criterion still provides a useful indicator of the
degree to which the available components suppress or amplify noise. When a new row
is selected, the number of zeroed signal components (n ? r) will decrease by one, and
the corresponding noise will be reflected in the criterion and increase the criterion value.
Therefore, we can choose the best new row based on the modified criterion such that the
criterion increases by the minimum value. When at least n rows are ultimately to be
selected, the error in the reconstructed signal due to the zeroed signal components will
eventually disappear, and the pseudo-inverse criterion will be equivalent to the original
inverse criterion. This justifies our neglect of the signal error component in using the
modified error criterion.
27
Having modified the criterion, we can no longer use the standard SFS algorithm to
optimize it. Note that AA
H
= U(??
H
)U
H
. Therefore,
E(A)=tr(A
H
A)
+
=
n
summationdisplay
i=1
1
?
2
i
=tr(AA
H
)
?1
(2.2.8)
Now, the form of the criterion given by (2.2.8) no longer needs a pseudo-inverse and
consequently is much more tractable as the basis of a selection algorithm.
2.3 SFS Algorithm
2.3.1 Standard SFS algorithm
Standard SFS begins with a matrix A with the number of rows at least equal to the
number of columns, or equivalently, with the number of selected observations at least
equal to the number of unknowns. It then sequentially selects the row that reduces the
MSE criterion the most at each step until the desired number of rows are selected. Let
A
+
represents the new matrix after the ith row a
i
is added to matrix A.Usingthe
Sherman-Morrison matrix inversion formula [110] gives
(A
H
+
A
+
)
?1
=(A
H
A)
?1
?
(A
H
A)
?1
a
H
i
a
i
(A
H
A)
?1
1+a
i
(A
H
A)
?1
a
H
i
(2.3.1)
Taking the trace of both sides gives
tr{(A
H
+
A
+
)
?1
} =tr{(A
H
A)
?1
}?
tr{(A
H
A)
?1
a
H
i
a
i
(A
H
A)
?1
}
1+a
i
(A
H
A)
?1
a
H
i
(2.3.2)
28
and using the property that tr{CDE} =tr{DEC},weobtain
tr{(A
H
+
A
+
)
?1
} =tr(A
H
A)
?1
?
a
i
(A
H
A)
?2
a
H
i
1+a
i
(A
H
A)
?1
a
H
i
(2.3.3)
Therefore, the criterion is minimized at each step by selecting the row (observation)
that maximizes the second term of (2.3.3).
2.3.2 Modified SFS algorithm
In this section, we derive an e?cient SFS algorithm for the modified criterion.
Lemma 2.3.1 [111] Let R denote an invertible partitioned matrix,
R =
?
?
?
?
?
?
?
?
P
.
.
. Q
??? ???
S
.
.
. T
?
?
?
?
?
?
?
?
(2.3.4)
where P and T are invertible. Then
R
?1
=
?
?
?
?
?
?
?
?
E
?1
.
.
. FH
?1
??? ???
H
?1
G
.
.
. H
?1
?
?
?
?
?
?
?
?
(2.3.5)
where
E = P ?QT
?1
S, (2.3.6)
F = ?P
?1
Q (2.3.7)
G = ?SP
?1
(2.3.8)
29
and
H = T ?SP
?1
Q (2.3.9)
The matrix E is called the Schur complement of P,andthematrixH is called the Schur
complement of T.
Lemma 2.3.2 [111] Let E denote the Schur complement of P,
E = P ?QT
?1
S. (2.3.10)
Then the inverse of E is
E
?1
= P
?1
+FH
?1
G (2.3.11)
Now denote A
+
as the new matrix after the ith row a
i
isaddedtomatrixA.Then
we have a partitioned expression of A
+
A
H
+
:
A
+
A
H
+
=
?
?
?
?
?
?
?
?
AA
H
.
.
. Aa
H
i
??? ???
a
i
A
H
.
.
. a
i
a
H
i
?
?
?
?
?
?
?
?
(2.3.12)
From (2.3.5), (2.3.9), and (2.3.11)
tr{(A
+
A
H
+
)
?1
} =tr(E
?1
)+tr(H
?1
)
=tr(AA
H
)
?1
+
1+a
i
A
H
(AA
H
)
?2
Aa
H
i
a
i
a
H
i
?a
i
A
H
(AA
H
)
?1
Aa
H
i
(2.3.13)
30
Therefore, the criterion is minimized at each step by selecting the row that minimizes
the second term in (2.3.13) [103].
2.4 E?cient Computation For SFS Algorithm
From(2.3.13), weknowthatweneedtocalculatea
i
A
H
(AA
H
)
?2
Aa
H
i
anda
i
A
H
(AA
H
)
?1
Aa
H
i
for every candidate row a
i
to determine the best one among all the candidate rows. If
we have many candidate rows, then the computational cost can be enormous. How-
ever, by making use of the previous calculations for the next choice, we can reduce the
computational burden.
Let v
i
=(AA
H
)
?1
Aa
H
i
and ?
i
= a
i
a
H
i
?a
i
A
H
(AA
H
)
?1
Aa
H
i
.Then
(A
+
A
H
+
)
?1
=
?
?
?
?
?
?
?
?
(AA
H
)
?1
+
1
?
i
v
i
v
H
i
.
.
. ?
1
?
i
v
i
??? ???
?
1
?
i
v
H
i
.
.
.
1
?
i
?
?
?
?
?
?
?
?
(2.4.1)
A
H
+
(A
+
A
H
+
)
?1
A
+
= A
H
(AA
H
)
?1
A
+
1
?
i
(A
H
v
i
v
H
i
A?a
H
i
v
H
i
A?A
H
v
i
a
i
+a
H
i
a
i
) (2.4.2)
Let u
i
= A
H
v
i
.Then
a
j
A
H
+
(AA
H
)
?1
A
+
a
H
j
= a
j
A
H
(AA
H
)
?1
Aa
H
j
+
1
?
i
|a
j
u
i
?a
j
a
H
i
|
2
(2.4.3)
(A
+
A
H
+
)
?2
=(A
+
A
H
+
)
?1
(A
+
A
H
+
)
?1
31
=
?
?
?
?
PQ
ST
?
?
?
?
(2.4.4)
P =(AA
H
)
?2
+
1
?
i
v
i
v
H
i
(AA
H
)
?1
+
1
?
i
(AA
H
)
?1
v
i
v
H
i
+
1
?
2
i
(v
i
v
H
i
)
2
+
1
?
2
i
v
i
v
H
i
(2.4.5)
Q = ?
1
?
i
(AA
H
)
?1
v
i
?
1
?
2
i
v
i
v
H
i
v
i
?
1
?
2
i
v
i
(2.4.6)
S = ?
1
?
i
v
H
i
(AA
H
)
?1
?
1
?
2
i
v
H
i
v
i
v
H
i
?
1
?
2
i
v
H
i
(2.4.7)
T =
1
?
2
i
v
H
i
v
i
+
1
?
2
i
(2.4.8)
Let z
i
=(AA
H
)
?1
v
i
and w
i
= A
H
z
i
.Then
A
H
+
(A
+
A
H
+
)
?2
A
+
= A
H
(AA
H
)
?2
A+
1
?
i
(u
i
w
H
i
+w
i
u
H
i
)
?
1
?
i
(a
H
i
w
H
i
+w
i
a
i
)
+
1+bardblv
i
bardbl
2
?
2
i
(u
i
u
H
i
?a
H
i
u
H
i
?u
i
a
i
+a
H
i
a
i
) (2.4.9)
a
j
A
H
+
(A
+
A
H
+
)
?2
A
+
a
H
j
= a
j
A
H
(AA
H
)
?2
Aa
H
j
+
2
?
i
[Re((a
j
w
i
)(a
j
u
i
)
H
)?Re((a
j
w
i
)(a
j
a
H
i
)
H
)]
+
1+bardblv
i
bardbl
2
?
2
i
|a
j
u
i
?a
j
a
H
i
|
2
(2.4.10)
32
Now, when we choose the bestrow from the remaining candidate rows to add to A
+
,
we can recursively update a
j
A
H
+
(A
+
A
H
+
)
?1
A
+
a
H
j
and a
j
A
H
+
(A
+
A
H
+
)
?2
A
+
a
H
j
for every
candidate row a
j
using (2.4.3) and (2.4.10). The initial values of these quantities can be
shown to be zero. Once u
i
and w
i
have been computed, the update formulas in (2.4.3)
and (2.4.10) only require vector inner products. The computational cost will be reduced
greatly in this way.
To calculate u
i
and w
i
, we need to calculate v
i
and z
i
first. To calculate v
i
=
(AA
H
)
?1
Aa
H
i
and z
i
=(AA
H
)
?1
v
i
, we need to deal with (AA
H
)
?1
. However, in this
problem, AA
H
is huge and cannot be stored or inverted. Therefore, instead of finding
the inverse matrix explicitly, we compute the solution by finding the minimizer of a
quadratic function of the form:
?(f)=||g ?Df||
2
(2.4.11)
where g is a known vector and D is a known matrix. We use the conjugate gradients
method [112] which is defined as follows:
p
0
= r
0
= D
H
g (2.4.12)
f
0
= 0 (2.4.13)
r
k
= ?
1
2
??(f
k
)=D
H
(g ?Df
k
) (2.4.14)
p
k
= r
k
+(||r
k
||
2
/||r
k?1
||
2
)p
k?1
(2.4.15)
f
k+1
= f
k
+?
k
p
k
(2.4.16)
33
?
k
=
(p
k
,r
k
)
||Dp
k
||
2
(2.4.17)
where p
k
is called the direction vector, and ?
k
is chosen to minimize ?(f)inthep
k
direction. Note that the conjugate gradients iteration does not require the storage of D
as long as the operation Dxcan be performed on an arbitrary vector x in some way other
than by a matrix-vector multiply. In this setting, Dx can be computed using FFT?s and
other simple operations. Since f has finite dimensions, the minimum can be reached
in finite steps. The vectors v
i
and u
i
will be the solutions corresponding to di?erent
definitions of g and D. Therefore, we can greatly reduce the computational time and
storage space by using the conjugate gradients algorithm to solve the inversion problem
instead of calculating the inverse directly.
2.5 Extending the SFS Process
When more samples than unknowns are desired, the selection can be continued
with the standard SFS algorithm in Section 2.3.1 based on (2.2.2). For computational
e?ciency, we derive a recursive method for updating a
i
(A
H
A)
?2
a
H
i
and a
i
(A
H
A)
?1
a
H
i
,
which are needed to evaluate the second term in (2.3.3). For completeness, we reproduce
the derivation in [100], along with our own extension that allows usto compute the initial
values of the recursion. Using the Sherman-Morrison matrix inversion formula, we have
the update formula for (A
H
+
A
+
)
?1
as follows:
(A
H
+
A
+
)
?1
=(A
H
A)
?1
?
(A
H
A)
?1
a
H
i
a
i
(A
H
A)
?1
1+a
i
(A
H
A)
?1
a
H
i
=(A
H
A)
?1
?
1
?
i
g
i
g
H
i
(2.5.1)
34
where g
i
=(A
H
A)
?1
a
H
i
.Leth
i
=(A
H
A)
?1
g
i
.Then
(A
H
+
A
+
)
?2
=(A
H
A)
?2
?
1
?
i
g
i
h
H
i
?
1
?
i
h
i
g
H
i
+
1
?
2
i
(g
i
g
H
i
)
2
(2.5.2)
a
j
(A
H
+
A
+
)
?1
a
H
j
= a
j
(A
H
A)
?1
a
H
j
?
1
?
i
|a
j
g
i
|
2
(2.5.3)
a
j
(A
H
+
A
+
)
?2
a
H
j
= a
j
(A
H
A)
?2
a
H
j
?
2
?
i
Re(a
j
g
i
(a
j
h
i
)
H
)+
bardblg
i
bardbl
2
?
2
i
|a
j
g
i
|
2
(2.5.4)
Now we need to find the initial value of a
i
(A
H
A)
?2
a
H
i
and a
i
(A
H
A)
?1
a
H
i
.WhenA is
square (the initial A for the extended SFS algorithm), we have the following relations
after some matrix manipulation:
a
j
(A
H
A)
?1
a
H
j
= a
j
A
H
(AA
H
)
?2
Aa
H
j
(2.5.5)
a
j
(A
H
A)
?2
a
j
= a
j
A
H
(AA
H
)
?3
Aa
H
j
(2.5.6)
The term a
j
A
H
(AA
H
)
?2
Aa
H
j
is already available from the initial SFS process from
(2.4.10). We also need an e?cient recursive formula to update a
j
A
H
(AA
H
)
?3
Aa
H
j
.
Let c
i
=(AA
H
)
?1
z
i
and d
i
= A
H
c
i
. Using (2.4.1) and (2.4.4), after some matrix
manipulation we have
(A
+
A
H
+
)
?3
=(A
+
A
H
+
)
?2
(A
+
A
H
+
)
?1
=
?
?
?
?
EF
GH
?
?
?
?
(2.5.7)
35
where
E =(AA
H
)
?3
+
1
?
i
(v
i
c
H
i
+c
i
v
H
i
)
+
1
?
i
z
i
z
H
i
+
1+bardblv
i
bardbl
2
?
2
i
(v
i
z
H
i
+z
i
v
H
i
)
+
v
H
i
z
i
?
2
i
v
i
v
H
i
+
(1+bardblv
i
bardbl
2
)
2
?
3
i
v
i
v
H
i
(2.5.8)
F = ?
1
?
i
c
i
?
v
H
i
z
i
?
2
i
v
i
?
1+bardblv
i
bardbl
2
?
2
i
z
i
?
(1+bardblv
i
bardbl
2
)
2
?
3
i
v
i
(2.5.9)
G = ?
1
?
i
c
H
i
?
v
H
i
z
i
?
2
i
v
H
i
?
1+bardblv
i
bardbl
2
?
2
i
z
H
i
?
(1+bardblv
i
bardbl
2
)
2
?
3
i
v
H
i
(2.5.10)
H =
1
?
2
i
v
H
i
z
i
+
(1+bardblv
i
bardbl
2
)
2
?
3
i
(2.5.11)
A
H
+
(A
+
A
H
+
)
?3
A
+
= A
H
(AA
H
)
?3
A+
1
?
i
w
i
w
H
i
+
1
?
i
(u
i
d
H
i
+d
i
u
H
i
?a
H
i
d
H
i
?d
i
a
i
)
+
1+bardblv
i
bardbl
2
?
2
i
(u
i
w
H
i
+w
i
u
H
i
?a
H
i
w
H
i
?w
i
a
i
)
+(
v
H
i
z
i
?
2
i
+
(1+bardblv
i
bardbl
2
)
2
?
3
i
)(u
i
u
H
i
?a
H
i
u
H
i
?u
i
a
i
+a
H
i
a
i
) (2.5.12)
36
a
j
A
H
+
(A
+
A
H
+
)
?3
A
+
a
H
j
= a
j
A
H
(AA
H
)
?3
Aa
H
j
+
1
?
i
|a
j
w
i
|
2
+
2
?
i
(Re((a
j
d
i
)(a
j
u
i
)
H
)?Re((a
j
d
i
)(a
j
a
H
i
)
H
))
+
2(1+bardblv
i
bardbl
2
)
?
2
i
(Re((a
j
w
i
)(a
j
u
i
)
H
)?Re((a
j
w
i
)(a
j
a
H
i
)
H
))
+(
v
H
i
z
i
?
2
i
+
(1+bardblv
i
bardbl
2
)
2
?
3
i
)|a
j
u
i
?a
j
a
H
i
|
2
(2.5.13)
Now, we can use this recursive formula to update a
j
A
H
+
(A
+
A
H
+
)
?3
A
+
a
H
j
until we
have the same number of samples as unknowns. Then we use this value as the initial
value of a
j
(A
H
A)
?2
a
j
to begin the recursion for the extended SFS algorithm.
2.6 Experiments
To demonstrate the value of our method, we applied it to an actual MR data set.
Figure 2.1(a) shows a 128 ? 128 MR scout image (courtesy of the Center for Nuclear
Imaging Research, University of Alabama at Birmingham) reconstructed from the full
128?128 k-space sampleset. Thek-space samplesofthescoutimage areconsidered to be
the full image data set for selection purposes. This data set was a rapidly acquired scout
image, which contains some significant artifacts outside the ROS. The presence of these
artifacts tests the ability of the algorithm to provide a good selection and reconstruction
in spite of violating the assumptions of the ROS.
Figure 2.1(b) shows the ROS identified by hand. It contains 5168 possibly nonzero
voxels. Using this ROS, the k-space data were selected optimally by the proposed algo-
rithm. Since the initial criterion with a single selected sample in the frequency domain
37
(a) (b)
Figure 2.1: (a) The real image, (b) The ROS image (white indicates the support region)
is the same regardless of the sample, we arbitrarily choose the DC sample as the initial
point for the recursive method. We have shown in the appendix that any circular shift
of the selection pattern in k-space yields the same criterion value. Therefore, a shift
of the initial point only causes a circular shift of the resulting optimal selection, which
is also optimal. To illustrate how the criterion changes as more samples are selected,
we continued the SFS process until all 16,384 k-space samples were selected. First, we
chose 5168 k-space samples (matching the number of points in the ROS) from 16,384
candidate points with the modified algorithm in Section 2.4. Then we switched to the
algorithm in Section 2.5 to continue the SFS process.
The criterion curve for the whole selection process is shown in Figure 2.2(a). When
matrix A is underdetermined, the criterion increases with each selected sample and
reaches a maximum value when A is square. This is because the pseudo-inverse solution
38
reduces the degrees of freedom for the noise, and the corresponding error resulting from
zeroing out some signal components is not included in the criterion. Equation (2.2.7)
shows that the MSE criterion for the minimum-norm least-squares solution includes two
terms when A is underdetermined ? one term represents the error due to the missing
signal componentsandtheother representingtheerror dueto observation noise. Asmore
samples are selected, the error from missing signal components decreases. When the
number of observations equals the number of unknowns, this term disappears entirely.
Since we know we will select the minimum density so that no signal components are
missed, we ignore this term in our modified criterion and only keep the noise term that
reflects the noise amplification level during the selection process. For each row of A
that is added, the number of degrees of freedom for noise increases by one, causing the
criterion to increase until A becomes square. Since the criterion only represents the noise
portion of the error and the missing signal portion is unknown, we must select rows for A
untiltheunknownsignalerrordisappears?thatis, whenAbecomessquare. Thereafter,
the criterion decreases with each selected sample because the average noise will decrease
when more samples are observed while the number of the unknowns is fixed. The flop
curve in Figure 2.2 (b) shows the relation of the criterion to the condition number
indirectly, which has the same shape as the criterion curve. The average number of flops
for the whole selection process is 2.1759?10
8
. The average number of flops is 3.086?10
8
before the minimum density and 1.7566?10
8
after it.
An interesting result can be observed from the criterion curve. The curve drops
sharply right after the maximum value and then becomes smooth. Therefore, the recon-
structed image will improve greatly if we select only a few more samples than the number
39
Table 2.1: Comparison of the criterion and time versus tol
tol=1e-5 tol=1e-3
#samples 5168 5297 16384 5168 5297 16384
MSE 10.8543 2.3362 0.3153 10.8542 2.3421 -0.3816
time(hr) 18 22 41 12 15 30
of unknowns. After this, selecting more samples yields diminishing returns in terms of
image quality improvement as a function of additional acquisition time. In Table 2.1, we
see that the criterion value is 10.8542 when the selected number equals the number of the
unknowns in the ROS (5168). Then with only 129 more samples selected, the criterion
drops drastically from 10.8542 to 2.3362; i.e, the net drop is 8.52. However, with 11216
more samples selected, the criterion decreases only by 2. Image quality improvement is
purely a function of the noise level decrement in the reconstructed image. There is no
direct relation of the quality criterion to the spatial resolution since the spatial resolu-
tion in the reconstruction is fixed; However, the samples optimized with the proposed
criterion covers the whole k-space nonuniformly. With this approach, we achieve the
same spatial resolution as that obtained with the full uniform k-space data.
The optimal sampling pattern with 5168 samples is shown in Figure 2.3(a) by the
black dots . The image reconstructed from these 5168 k-space samples is shown in Figure
2.4(a). To demonstrate the improvement obtained by selecting a few more samples than
the ROS, we reconstructed the image by selecting 129 more noisy samples than the ROS
(5297 k-space points), where the criterion drops from 10.85 to 2.34. The corresponding
sample pattern is shown in Figure 2.3(a) by 129 gray dots. The reconstructed image
is shown in Figure 2.4(b). The di?erence image between the optimized image from
40
0 2000 4000 6000 8000 10000 12000 14000 16000 18000
0
2
4
6
8
10
12
MSE criterion
Modified
(a)
0 2000 4000 6000 8000 10000 12000 14000 16000 18000
0
1
2
3
4
5
6
7
8
9
10
x 10
9
(b)
Figure 2.2: (a) The criterion curve versus the number of selected samples with the
128x128 ROS image, (b) The flop curve versus selection step with the 128x128 ROS .
41
5297 k-space samples and the original image is shown in Figure 2.3 (b). The artifacts
outside the ROS are the systematic di?erence between them. When k-space samples are
acquired independently as in MRSI, the noise variance at each k-space location is the
same. However, our data set was acquired by a rapid scan through k-space in which
the samples are not acquired independently. In addition, the artifacts outside the ROS
are low-frequency noise, which is not considered in the model. Thus, the noise is larger
at the lowest, most powerful frequencies in the acquisition process. Therefore, once we
foundtheoptimal sampledistribution, wecircularly shiftedit sothatDC and afew other
adjacent frequencies were unselected. Since according to the theorem in the appendix
the criterion is invariant to periodic shifts of the selection, this step will have no negative
e?ects on reconstruction quality.
The discrete impulse response for the optimized set is an impulse for every location
in the ROS. The impulse response from the optimized k-space data is shown in Figure
2.5(a). That explains why we can reconstruct the image exactly at the minimum den-
sity if there is no noise outside the ROS. For comparison purposes, we also used 5297
low-frequency noisy samples in a circle around DC to reconstruct the image from a ze-
ropadded FFT and a least-squares method. The impulse response from the inverse FFT
of the lowpass filter is shown in Figure 2.5(b). Only about 1/3 of the energy is concen-
trated on the center point with 2/3 of the energy spread out to form a lowpass impulse
response. The lowpass image from a zero padded FFT has about half the resolution
of the optimized image. There are obvious rings in the reconstructed image from a ze-
ropadded FFT as expected. Therefore, the zeropadded FFT reconstruction of the image
is not included here. Instead, iterative reconstructions of the lowpass image at 3 di?erent
42
20 40 60 80 100 120
20
40
60
80
100
120
(a) (b)
Figure 2.3: (a) 5297 optimized k-space sample distribution (the 5168 black dots represent
the optimized minimum density, larger gray dots indicates 129 more optimized samples
than the minimum density), (b) Absolute di?erence between the optimized image recon-
struction and the fully sampled 128?128 image. (Dynamic range is scaled by a factor
of 5 for greater visibility).
43
(a) (b)
Figure 2.4: (a) The image from the minimum density k-space data(5168 points). (b)
The image from 5297 k-space samples (the ROS has been imposed as a constraint in the
reconstruction).
iterations are shown in Figure 2.6. Figure 2.6(a) shows the reconstructed image from
an inverse FFT of 5297 lowpass samples. Figure 2.6(b) shows the reconstructed images
from a least-squares solution at iteration 5. Figure 2.6(c) shows the reconstructed im-
ages from a least-squares solution at iteration 20. Figure 2.6(d) shows the reconstructed
images from a least-squares solution at iteration 100.The resulting system is extremely
ill-conditioned and gets progressively worse as the iterations continue.
The tolerance for ending the conjugate gradients algorithm has some e?ect on the
accuracy of the computational criterion and the computational time which is shown in
Table 2.1. Here, we define the tolerance as the gradient vector norm of a quadratic
function that needs to be minimized. The function reaches the minimum when the
gradient vector norm equals zero. In theory, the MSE with all 16,384 points observed is
44
5
10
15
5
10
15
0
0.5
1
5
10
15
5
10
15
0
0.5
1
(a) (b)
Figure 2.5: (a) The point spread function from optimized k-space data (b) The point
spread function from inverse FFT.
0.3154. When tol =10
?5
, the final computational criterion is 0.3153, almost the same as
the theoretical value. When tol =10
?3
, there is an obvious accumulating computational
error to the criterion when a huge number (16384) of samples is selected as shown in
Table 2.1, and the sample selection pattern is a little di?erent from that for tol =10
?5
.
When the selection number is not quite as large, (5297 in Table 2.1 and 4096 or 2304
in Table 2.2), the accumulating error is acceptable. To finish the whole process takes 30
hours when tol =10
?3
and41hourswhentol =10
?5
(running on a 250MHz Sun Ultra
Sparc processor). The smaller the FOV is, the faster the selection process will be. Here,
we chose to use tol =10
?5
to assure reliability of the results.
To evaluate the method for MRSI, we applied the method to an actual 1H MRSI
data set. Figure 2.7(a) and (b) shows the water and fat image from a full 64 ?64 real
k-space MRSI data set (courtesy of the Center for Nuclear Imaging Research, University
of Alabama at Birmingham) and the ROS identified by hand. It contains 1439 possibly
45
(a) (b)
(c) (d)
Figure 2.6: (a) Lowpass image from FFT, (b) Lowpass recursive solution at iteration 3,
(c) Lowpass recursive solution at iteration 20, (d) Lowpass recursive solution at iteration
100.
46
Table 2.2: Comparison of the selection criterion and time versus FOV with MRSI ROS
image
64?64 FOV 48?48 FOV
#samples 1439 1513 4096 1439 1513 2304
tol= MSE 6.4186 1.8906 0.3513 8.4318 2.2886 0.6245
1e-5 time(hr) 1.13 1.39 2.47 0.77 0.92 1.15
tol= MSE 6.1709 1.9027 0.3184 8.4318 2.2887 0.6440
1e-3 time(hr) 0.72 0.9 1.64 0.497 0.598 0.757
nonzerovoxels. UsingthisROS,thek-spacedatawereselected optimallybytheproposed
algorithm. Figure 2.7(c) shows the ROS identified by hand. Using this ROS, the k-space
data were selected optimally by the proposed algorithm. To illustrate how the criterion
changes asmore samples are selected, we continued the SFS processuntil all 4096 k-space
samples were selected. First, we chose 1439 k-space samples (matching the number of
points in the ROS) from 4096 candidate points with the modified algorithm in Section
2.4. Then we switched to the algorithm in Section 2.5 to continue the SFS process.
To compare the results with a reduced FOV, we applied the proposed method and the
traditional lowpass method to a reduced 48?48 FOV. Figure 2.7 (d) shows the 48?48
ROS image extracted from the 64 ? 64 ROS image. Figure 2.7 (e) shows the 48 ? 48
water image extracted from the 64?64 real image. Figure 2.7 (f) shows the 48?48 fat
image extracted from the 64?64 real image. The results are summarized in Table 2.2.
The set of 1513 optimized k-space samples with the 64?64 FOV is plotted in Figure
2.8(a). Figure 2.8(b) shows the 1513 optimized k-space data set distribution with the
48?48 FOV. The 1439 black dots represent the optimized minimum density sampling.
The larger gray dots represent an additional 74 optimized samples to reduce the noise.
Both distributions cover all of k-space with the larger interval in the 48?48 case.
47
(a) (b)
(c) (d)
(e) (f)
Figure 2.7: (a) Water image from the full 64 ?64 real MRSI data, (b) Fat image from
the full 64?64 real MRSI data, (c) The 64?64 ROS image (white indicates the support
region), (d) 48 ? 48 ROS image extracted from the 64 ? 64 ROS (white indicates the
support region), (e) 48?48 water image extracted from the 64 ?64 real water image,
(f) 48?48 fat image extracted from the 64?64 real fat image.
48
10 20 30 40 50 60
10
20
30
40
50
60
5 10 15 20 25 30 35 40 45
5
10
15
20
25
30
35
40
45
(a) (b)
Figure 2.8: (a) 1513 optimized k-space data distribution with 64?64 FOV (the minimum
optimized 1439 points indicated by black dots, larger gray dots indicate 74 more samples
than the minimum density), (b) 1513 optimized k-space samples with a reduced 48?48
FOV (the minimum optimized 1439 points indicated by black dots, larger gray dots
indicate 74 more samples than the minimum density).
49
The noise criterion curve versus the number of selected samples with the reduced
48?48 ROS image is shownin Figure 2.9(a). Thecriterion value dropssharplyrightafter
the maximum value at the minimum density. From Table 2.2, we see that the selection
time reduced greatly with smaller FOV, making it possible for real-time implementation.
The flop curve is shown in Figure 2.9(b). The average number of flops per recursive step
for the whole selection process is 8.0730?10
7
. The average number of flops per recursive
step is 8.6022 ? 10
7
before the minimum density and 7.1926 ? 10
7
after the minimum
density.
The optimized PSF with a reduced 48?48 FOV is shown in Figure 2.10(a), which is
an impulse. The optimized PSF with the 64?64 FOV is also an impulse (not include).
The PSF for 1513 lowpass samples with the reduced 48 ? 48 FOV is shown in Figure
2.10(b). The PSF for 1513 lowpass samples with the 64 ? 64 FOV is shown in Figure
2.10(c). From Figure 2.10(b) and Figure 2.10(c), we see that the percent of the energy on
the point is proportional to the percent of the lowpass area in the FOV for the lowpass
PSF, while the width of the lowpass PSF is inversely proportional to the area of the
lowpass region. However, with the circular or elliptical ROS, the lowpass image cannot
reconstruct the image exactly.
The water and fat reconstruction images from 1513 optimized k-space samples with
a64?64 FOV are shown in Figure 2.11 (a) and (b). Thecorrespondingdi?erence images
between the optimal image and the original image are shown in Figure 2.12. The water
and fat reconstruction images from 1513 optimized k-space samples with the reduced
48?48 FOV are shown in Figure 2.13 (a) and (b). The corresponding di?erence images
50
0 500 1000 1500 2000 2500
0
1
2
3
4
5
6
7
8
9
MSE criterion
Modified
(a)
0 500 1000 1500 2000 2500
0
2
4
6
8
10
12
14
x 10
8
(b)
Figure 2.9: (a) Noise criterion curve versus the number of selected samples with the
48 ? 48 ROS, (b) Flops per selection step versus the number of selected samples with
the 48?48 ROS.
51
0
5
10
15
0
5
10
15
0
0.2
0.4
0.6
0.8
1
0
5
10
15
0
5
10
15
0
0.2
0.4
0.6
0.8
1
0
5
10
15
0
5
10
15
0
0.2
0.4
0.6
0.8
1
(a) (b) (c)
Figure 2.10: (a) Optimized PSF inside the ROS for a reduced 48?48 FOV, (b) Lowpass
PSF for 48?48 FOV with 1513 lowpass k-space samples, (c) Lowpass PSF for 64?64
FOV with 1513 lowpass k-space samples.
between the optimal image and the original image are shown in Figure 2.14. The lowpass
water and fat image from 1513 k-space samples are shown in Figure 2.13 (c) and (d).
2.7 Discussion
The work presented here is important for optimizing MRSI data acquisition and
can also be used for 3-D MRI. The new SFS algorithm presented in this paper has the
following advantages: 1) it overcomes the limitation of the SFS algorithm based on the
standard MSE criterion in that it can start with an empty initial matrix; 2) it drastically
reducesstoragespaceandcomputationaltimewiththeproposedSFSrecursivealgorithm
and the conjugate gradients algorithm; 3) it is more applicable in practice than SBS in
that the best sample can be observed as soon as it is selected, making possible real-time
selection and acquisition. Simulation with real MR data shows a good-quality image
can be reconstructed with only a few more samples than the minimum density by our
method. Further work is clearly necessary to accelerate the algorithm if it is to be useful
52
(a) (b)
Figure 2.11: (a) Water image reconstructed from the 1513 optimized k-space samples
with 64?64 FOV with ROS constraint imposed, (b) Fat image reconstructed from the
1513 optimized k-space samples with 64?64 FOV with ROS constraint imposed.
(a) (b)
Figure 2.12: (a) Absolute di?erence between the optimized water image reconstruction
and the fully sampled 64?64 image (b) Absolute di?erence between the optimized fat
image reconstruction and the fully sampled 64?64 image. (Dynamic range is scaled by
a factor of 5 for greater visibility).
53
(a) (b)
(c) (d)
Figure 2.13: 48?48 FOV reconstruction: (a) Water image reconstructed from the 1513
optimized k-space samples with ROS constraint imposed, (b) Fat image reconstructed
from the 1513 optimized k-space samples with ROS constraint imposed, (c) Water image
from the 1513 lowpass k-space samples, (d) Fat image from the 1513 lowpass k-space
samples.
54
(a) (b)
Figure 2.14: (a) Absolute di?erence between the optimized water image reconstruction
and the 48?48 extracted image, (b) Absolute di?erence between the optimized fat image
reconstruction and the 48?48 extracted image. (Dynamic range is scaled by a factor of
5 for greater visibility).
in a real-time context on an ordinary computer. However, the algorithm described here
brings real-time implementation closer to reality.
Chapter 3
Efficient Backward Selection of Hexagonal Samples in k-space
3.1 Introduction
Hexagonal sampling at the Nyquist rate has a 13.4% e?ciency gain over rectangular
sampling at the Nyquist rate for an elliptical ROS and yields the most densely packed
spatial replications for space-limited signals [90,91]. However, Nyquist density sampling
cannot reach the minimum density for an elliptical ROS image. Minimum density sam-
pling yields samples linearly independent of each other, which can be reached through a
nonuniform sampling strategy [101].
Reeves et al. proposed sequential backward selection (SBS) on a rectangular grid
[98,99,103] based on minimizing mean squared error (MSE) in the reconstructed image.
SBS begins with a full set of candidate samples and sequentially eliminates one sample
at a time until the desired number remain. SBS represents a suboptimal approach to
determine an optimized set of sample locations. However, SBS can guarantee a given
level of performancebefore actually performingthe optimization [100,113]. Furthermore,
in extensive simulations we have found that SBS provides nearly optimal results in every
case. By comparison, random or heuristic selection methods may yield wildly erratic
results. To reduce selection time, we present sequential backward selection (SBS) on a
hexagonal grid.
The motivation for using a hexagonal grid is that we can begin the elimination
process with fewer candidate samples than in the rectangular case. Thus, we have fewer
samples to eliminate to reach the desired number of remaining samples. It would be
55
56
possible simply to eliminate randomly some samples from a rectangular grid candidate
set so that there are fewer samples to eliminate. However, there is no guarantee that this
initial candidate set would be a good one from which to begin the backward selection
process. Reeves andZhaohave shownthattheupperboundonMSEof thefinalselection
is proportional to the MSE resulting from the initial candidate set [100]. The initial
candidate set defined by a hexagonal grid yields the lowest possible MSE for the same
number of samples. Thus, the upper bound on the MSE of the final selection resulting
from a hexagonal candidate grid is as low as possible when beginning with the same
number of candidate samples. As such, it is an ideal starting point for the backward
selection process.
Standard fast Fourier transform (FFT) algorithms cannot be used directly for non-
rectangularly sampled data. To make use of the standard fast Fourier transform, we use
a Smith normal form decomposition to transform the hexagonal DFT into a rectangular
DFT. Since standard display devices require images to be displayed on a rectangular
grid, the proposed algorithm yields an image sampled on a rectangular grid directly
without interpolation by properly setting the periodic matrix in the frequency domain.
3.2 Selection algorithm
Sequential backward selection begins with a candidate matrix and sequentially elim-
inates the least important row at each step until the desired number of rows remain. Let
a
i
represent row i of A. If we eliminate a
i
from A, the inverse matrix (
?
A
H
?
A)
?1
is given
57
by the Sherman-Morrison matrix inversion formula [110] as
(
?
A
H
?
A)
?1
=(A
H
A)
?1
?
(A
H
A)
?1
a
H
i
a
i
(A
H
A)
?1
1+a
i
(A
H
A)
?1
a
H
i
(3.2.1)
Taking the trace of both sides gives
tr{(
?
A
H
?
A)
?1
} =tr{(A
H
A)
?1
}+
tr{(A
H
A)
?1
a
H
i
a
i
(A
H
A)
?1
}
1?a
i
(A
H
A)
?1
a
H
i
(3.2.2)
and using the property that tr{CDE} =tr{DEC},weobtain
tr{(
?
A
H
?
A)
?1
} =tr(A
H
A)
?1
+
a
i
(A
H
A)
?2
a
H
i
1+a
i
(A
H
A)
?1
a
H
i
(3.2.3)
Therefore, the criterion is minimized at each step by selecting the row (observation) that
minimizes the second term of (3.2.3). With this recursive formula, we need only evaluate
a
i
(A
H
A)
?1
a
H
i
and a
i
(A
H
A)
?2
a
H
i
for all i to determine the row to be eliminated rather
than computing the matrix inverse (A
H
A)
?1
with each row a
i
in turn being eliminated.
We can follow a procedure similar to the one in [103] and described in Chapter
2 to recursively compute the quantities needed to evaluate the second term in (3.2.3).
Let g
i
=(A
H
A)
?1
a
H
i
, h
i
=(A
H
A)
?1
g
i
,and?
i
=1? a
i
(A
H
A)
?1
a
H
i
. To avoid storing
a huge matrix, we compute g
i
and h
i
iteratively by minimizing ?(v
i
)=||e
i
? Av
i
||
2
and ?(w
i
)=||v
i
? A
H
Aw
i
||
2
with the conjugate gradients method [114]. With these
quantities in hand, we can e?ciently compute
a
j
(
?
A
H ?
A)
?1
a
H
j
= a
j
(A
H
A)
?1
a
H
j
+
1
?
i
|a
j
g
i
|
2
(3.2.4)
58
a
j
(
?
A
H ?
A)
?2
a
H
j
= a
j
(A
H
A)
?2
a
H
j
+
2
?
i
Re(a
j
g
i
(a
j
h
i
)
H
)+
|g
i
|
2
?
2
i
|a
j
g
i
|
2
(3.2.5)
Note that v
i
and w
i
only need to be computed for the row being eliminated in each elim-
ination step, since every remaining row can be evaluated from the recursively computed
terms in (3.2.4) and (3.2.5).
To use the recursion formulas, we must have initial conditions. Assume that A
represents a Fourier transform matrix with columns removed corresponding to voxels
outsidetheROS.Lettheimage sizebemandtheROSsizeben. A
H
A = mI becausethe
columns of A are orthogonal to each other. Therefore, the initial value of a
j
(A
H
A)
?1
a
H
j
is n/m for all j. The initial value of a
j
(A
H
A)
?2
a
H
j
is n/m
2
for all j. Since the criterion
increment is the same for all j, we may arbitrarily choose any sample as the first to be
eliminated.
3.3 Regular Hexagonal Sampling
All MR imaging methods use discrete data that has been sampled from an idealized
continuous data function. The data sampling pattern determines the image replication
pattern. If X
c
(d) represent the continuous data function, the sampled data X at an
integer sequence of points n may be expressed as
X(n)=X
c
(Vn) (3.3.1)
where V is the sampling matrix whose columns are the basis vectors of the sampling
pattern [115]. The hexagonal sampling pattern is shown in Figure 3.1(b). The axis
59
direction is drawn in an orientation that conforms to an image matrix. For rectangular
sampling, V is a diagonal matrix. For hexagonal sampling,
V =
?
?
?
?
2d
1
?d
1
0 d
2
?
?
?
?
(3.3.2)
where d
1
= d
2
/
?
3. The sampling matrix is not unique for a given sampling pattern
because di?erent basis vectors can be chosen. The sampling matrix in (3.3.2) has the
advantage of creating the rectangular image pixels directly from a discrete Fourier trans-
form without interpolation. To keep pixels on a square array, d
2
/d
1
must be a rational
fraction rather than the optimal value of
?
3 [116]. Use of d
2
/d
1
=7/4 yields a near-
optimal e?ciency gain of 12.5% compared to 13.4% for d
2
/d
1
=
?
3. When sampling
in the frequency domain, the reconstructed image from a discrete Fourier transform is
replicated periodically in the spatial domain. The periodic matrix in the spatial domain
is
N =2?(V
T
)
?1
=
?
?
?
?
2?
2d
1
0
2?
2d
2
2?
d
2
?
?
?
?
=
?
?
?
?
7
8
M 0
M
2
M
?
?
?
?
(3.3.3)
where M ? 2?/d
2
is an integer divisible by 8 that represents the length of one side of the
square image. On the other hand, discrete sampling of the image in the spatial domain
causes periodic replication of the sample pattern in the frequency domain. The periodic
60
matrix in the frequency domain is
N
T
=
?
?
?
?
7
8
M
M
2
0 M
?
?
?
?
(3.3.4)
It should be noted that the periodic matrix in (3.3.4) should be an integer matrix. It
defines a rectangular replication pattern in the frequency domain, although its sampling
pattern is hexagonal. The periodic matrix in (3.3.3) defines a hexagonal replication
pattern in the spatial domain. Figure 3.1(a) shows as an example of the rectangular
replication of the 4 ? 6 representative hexagonal samples within the region defined by
the diamond symbols. It is known that a rectangular replication pattern in one domain
corresponds to a rectangular sampling pattern in the transformed domain, and a hexago-
nal replication pattern in one domain correspondsto a hexagonal sampling pattern in the
transformed domain regardless of the sampling pattern in the original space. Therefore,
the periodic matrix in (3.3.4) defines image pixels on a rectangular grid in the spatial
domain.
For the image to be exactly recoverable from the hexagonal sequence defined by
(3.3.1), it must be space-limited with a hexagonal ROS as shown in Figure 3.1(b) with
d
1
<
?
b
,d
2
<
4?
2a+c
(3.3.5)
When a = c =
2?
?
3
b, the ROS is a regular hexagon. However, the horizontal and
vertical sampling interval are not the same even for a regular hexagon. Consider a space-
limited waveform with a circular ROS. If b denotes the radius of this region, then for this
waveform sampled rectangularly the maximum sampling periods are d
1
= ?/b, d
2
= ?/b.
61
c
b
a
( a )
( b )
v1
)(
0
4
N
T
= [
4
0
3
6
]
0
3
N = [
4
6
]
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
2d2
2d1
( )
3
6
Figure 3.1: (a) Hexagonal sampling pattern, (b) Hexagon
62
For sampling this same waveform hexagonally we can choose d
1
= ?/b, d
2
=2?/
?
3b.
Themean sampling densityis lessfor hexagonal sampling. In general, the mean sampling
density is proportional to the area of the assumed ROS. For a circularly space-limited
waveform, hexagonal sampling involves 13.4% fewer samples than rectangular sampling.
If we use SBS beginning with a full hexagonal grid rather than a full rectangular grid, the
selection time will be reduced because fewer samples must be eliminated. The resolution
will not be reduced.
3.4 Generalized Multidimensional Discrete Fourier Transform
An FFT from a hexagonal grid is necessary both for the selection algorithm as
well as for reconstruction once the samples are acquired. Many FFT algorithms have
been developed for the evaluation of rectangular multidimensional DFTs (MDFTs). We
present a procedure to make use of these algorithms for general MDFTs. A generalized
MDFT pair is
x(n)=
1
|detN|
summationdisplay
k
X(k)e
jk
T
(2?N
?1
)n
(3.4.1)
X(k)=
summationdisplay
n
x(n)e
?jk
T
(2?N
?1
)n
(3.4.2)
whereN is an integer periodic matrix in (3.3.3). By using the Smith form decomposition
[117,118], N can be decomposed into
N = UDV (3.4.3)
63
where U and V are unimodular integer matrices (determinant is ?1) and D is diago-
nal with |det(D)| = |det(N)|. To get desired Smith forms, two steps are performed.
First, using elementary (regular unimodular) matrices, multiply the given integer ma-
trix recursively on the left and right until the matrix is reduced to a diagonal form with
the diagonal entry values in increasing order. Elementary row and column operations
correspond to multiplication of the matrix on the left and right respectively by suitable
unimodular matrices. Initially, U = I, D = N,andV = I. The flowchart of the Smith
form decomposition is shown in Figure 3.2. The diagonal matrix D obtained in the first
step may not be what we desire. However, we can impose structure on the available
Smith forms by pivoting factors between diagonal elements of D [119]. Specifically, ma-
trix product PDS can move an integer factor i from one diagonal entry to the other
diagonal entry with det(P)=?1anddet(S)=?1.
With Smith forms in (3.4.3), the MDFT pair in (3.4.1) and (3.4.2) can be expressed
as
X(k)=
summationdisplay
n
x(n)e
?jk
T
(2?V
?1
D
?1
U
?1
)n
=
summationdisplay
n
x(n)e
?j(V
?T
k)
T
D
?1
(U
?1
n)
=
summationdisplay
hatwiden
x(Uhatwiden)e
?j
hatwide
k
T
D
?1
hatwiden
=
hatwide
X(
hatwide
k)
=
hatwide
X(V
?T
k) (3.4.4)
x(n)=
summationdisplay
k
X(k)e
jk
T
(2?V
?1
D
?1
U
?1
)n
64
set the minimum at  x(1,1)
entries on k row and
and k column except
x(k,k)?
zero out all
k = N - 1 ?
   end
add column j
to column k
if x(i,j) not
divisible
get quotient
of each entriy
with x(k,k)
the minimum
at x(k,k)
set the minimum
at x(k,k)
x(k,k)?
entries divisible 
by
k = 1
k = k+1
No
Yes
Yes
No
Yes
Yes
No
No
 all remaining
Figure 3.2: Smith form decomposition flowchart
65
=
summationdisplay
k
X(k)e
j(V
?T
k)
T
D
?1
(U
?1
n)
=
summationdisplay
hatwide
k
X(V
T
hatwide
k)e
j
hatwide
k
T
D
?1
hatwiden
= hatwidex(hatwiden)
= hatwidex(U
?1
n) (3.4.5)
With (3.4.4) and (3.4.5), a matrix-N DFT becomes a matrix-D DFT. The general
MDFT problem converts to a rectangular MDFT problem. The matrix-N MDFT al-
gorithm involves three steps. First, scramble the input sequence x(n)accordingtothe
relation n =(Uhatwiden)modN. Second, compute a rectangular matrix-D DFT of the re-
sulting sequence. Third, unscramble the output sequence X(
hatwide
k) according to the relation
hatwide
k =(V
?T
k)modD.Thematrix-N inverse MDFT involves three similar steps.
3.5 Experiments
To demonstrate the value of our method, we applied it to an actual MR data set.
Figure 3.3(a) shows a 96?96 MR image from a 128?128 rapidly acquired scout image
(courtesy of the Center for Nuclear Imaging Research, University of Alabama at Birm-
ingham) reconstructed from the full 128?128 k-space sample set. Figure 3.3(b) shows
a96?96 ROS identified from Figure 3.3(a). It contains 5168 possibly nonzero voxels.
For the ROS in 3.3(b), the k-space data were selected by SBS from a 96 ? 96 rectan-
gular k-space grid, which correspond to samples of the rectangular DFT coe?cients of
the image in Figure 3.3(a). Likewise, k-space samples were selected by SBS from an
84 ?96 hexagonal k-space grid, which correspond to samples of the hexagonal DFT of
the equivalent 84?96 image.
66
(a) (b)
(c) (d)
Figure 3.3: (a) 96?96 image, (b) 96?96 ROS, (c) Reconstructed image from 5550 rect-
angular k-space samples, (d) Reconstructed image from 5550 hexagonal k-space samples.
67
5000550060006500700075008000850090009500
0
5
10
15
# of remain samples
MSE
rectangular
hexagonal
Figure 3.4: Criterion curves versus samples remaining.
5000550060006500700075008000850090009500
0
1
2
3
4
5
6
# of remain samples
hours
rectangular
hexagonal
Figure 3.5: Time versus samples remaining.
68
Figure 3.4 shows the criterion comparison between SBS on a rectangular grid and
SBSonahexagonal grid. Figure3.5 showstheselection timecomparison between SBSon
a rectangular grid and a hexagonal grid. We observed that when the selection approaches
5168 samples, the criterion curves and selection time curves increase sharply but increase
very little when stopping with only a few more samples than 5168. The two criterion
curves are very close throughout the process. However, the selection time for SBS from
a hexagonal grid is always significantly less than that for SBS from a rectangular grid.
When selecting more than 6500 samples, the selection time in the hexagonal case takes
less than half of the selection time of the rectangular case. These observations tell us we
can sample at a rate lower than the Nyquist rate on a hexagonal grid with minimal noise
amplification as long as we acquire a few more samples than the number of nonzero pixels
in the ROS. As a further comparison, we reconstructed the image from 5500 selected
samples for both the rectangular and hexagonal cases. The reconstructions were done
iteratively using a conjugate gradients algorithm. The results shown in Figure 3.3(c)
and (d) demonstrate the good results that can be obtained with either selection method.
Although there is no di?erence in resolution for the two cases and little di?erence in
noise amplification, the selection time is much smaller for SBS on a hexagonal grid as
showninFigure3.5.
3.6 Conclusion
We proposed a new method for reducing selection time for MRSI samples by apply-
ing SBS on a hexagonal grid rather than a rectangular grid. To make use of a standard
69
rectangular FFT, we computed a generalized MDFT based on a Smith form decompo-
sition. To avoid interpolation artifacts, we chose an appropriate periodic matrix in the
frequency domain such that we can compute rectangular-grid image pixels directly from
hexagonal-grid frequency domain samples without interpolation. A lower than Nyquist
rate on a hexagonal grid can be achieved with the proposed method. Furthermore, SBS
on a hexagonal grid takes significantly less selection time than SBS on a rectangular
grid.
Chapter 4
Sequential Backward Array Selection
4.1 Introduction
Inthischapter, wedevelop ane?cient algorithmforoptimizingtheditheringpattern
so that the image can be reconstructed as reliably as possiblefrom a periodicnonuniform
setof samples, which can beobtained fromaditheredrectangular-grid array. Takinginto
account the frequency support of the image, we sequentially eliminate the least informa-
tive array recursively until the minimal number of arrays remain. The advantage of the
proposed algorithm is that convolution and downsampling take the place of operations
involving a huge non-sparse matrix. Therefore, the proposed algorithm shows high e?-
ciency in computational time and memory space. The resulting algorithm can be useful
not only in magnetic resonance imaging but is also extremely important in some sensor
array optimization problems, in digital reception of multiband signal transmissions, and
in PMMW imaging [106,107].
4.2 The Mathematics of Sampling
Typical MRSI observations are laid out in a rectangular grid pattern as shown
in Figure 4.1. The heavy dots represent the locations of the sensor elements in the
unshifted sensor array. The light dots represent other locations to which the sensor
array can be shifted. Assuming a circular ROS, if the radius of the ROS is R m, the
Nyquist sample spacing in each dimension considered individually is
?
R
. This is the
maximum uniform rectangular sample spacing that avoids aliasing. However, it is well
70
71
Figure 4.1: typical sampling pattern.
known that a circularly bandlimited image can be sampled uniformly on a hexagonal
grid with a 13.4% lower sampling density without aliasing [90,91], where ?x =2?/
?
3R
and ?y = ?/R. (See Figure 4.2.)
In fact, if we allow nonuniform sampling we can reduce the average sampling density
even further [101]. In particular, we can sample an image using a periodic replication of a
nonperiodicsampling pattern to capture the information in a bandlimited image without
aliasing. An example of this kind of sampling pattern is shown in Figure 4.3. In this
example, each 3?3 block sampling pattern is periodically replicated over the frequency
image plane. This pattern results from shifting the array in Figure 4.1 according to each
o?set sample in one of the 3?3blocks.
If we represent the discretized frequency image as a vector y, the unknown dis-
cretized spatial support by a vector x, and the mapping from the spatial-domain samples
to the frequency-domain image byF, the fullysampledfrequency image can beexpressed
72
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
 
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?       
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?    ?    ?    ?    ?    ?    ?    ?    ?    ?    ?
?
?
?
?
?
?
?x
?y
ky 
kx , n1
n2
Figure 4.2: Hexagonal sampling pattern.
as
y = Fx+u (4.2.1)
where u is zero-mean, i.i.d. noise. If we sample the image y using a single position of the
sensor array with o?set indexed by i, we can represent this by
y
i
= Q
i
y (4.2.2)
where Q
i
downsamples the fully sampled image and orders the result into a vector. Then
we can rewrite y in a rearranged form y
r
as
y
r
=[y
H
1
y
H
2
... y
H
n
]
H
(4.2.3)
=[Q
H
1
Q
H
2
... Q
H
n
]
H
(Fx+u)
= Q
r
Fx+u
r
73
Figure 4.3: Periodic nonuniform sampling example.
where u
r
is the similarly rearranged version of u and Q
r
=[Q
H
1
Q
H
2
... Q
H
n
]
H
.Ifwe
choose a subset of k of the n shifted arrays, we obtain
?y
r
=[y
H
i
1
y
H
i
2
... y
H
i
k
]
H
(4.2.4)
=[Q
H
i
1
Q
H
i
2
... Q
H
i
k
]
H
(Fx+u)
=
?
Q
r
Fx+?u
r
where ?u
r
is the corresponding subset of u
r
and
?
Q
r
is the corresponding subset of Q
r
.
As long as the {Q
i
} are properly chosen and enough subsets are selected, the un-
known spatial samples x can be reconstructed from ?y
r
in a least-squares sense by
hatwidex =(F
H
?
Q
H
r
?
Q
r
F)
?1
F
H
?
Q
H
r
?y
r
(4.2.5)
Thepreciseconditions underwhichaleast-squares reconstruction is definedarediscussed
in [101] and [95]. A necessary condition is that the product of the number of samples
74
in the array and the number of selected array positions equals or exceeds the number
of unknown frequency samples. Assuming that u has unit variance, the sum of squared
errors (SSE) in the reconstructed spatial samples is given by
E(
?
Q
r
)=tr(F
H
?
Q
H
r
?
Q
r
F)
?1
(4.2.6)
4.3 Basic Optimization Algorithm
We must optimize (4.2.6) with respect to the entries of
?
Q
r
to provide the most
useful information for reconstruction of images in the minimal necessary acquisition
time. Unfortunately, optimization of (4.2.6) is a large-scale combinatoric optimization
problem, and an optimal solution is prohibitive to find. We use sequential backward
selection (SBS) to optimize (4.2.6). The key challenge in this optimization task is that
the available SBS algorithms do not apply directly, since they are designed to eliminate
one sample location at a time [98]. Instead, we must eliminate a whole periodic array
of samples at one time here. This severely complicates the relevant equations and the
resulting computational complexity of an already complex algorithm.
Without loss of generality, we assume that F is scaled so that it has orthonormal
columns. Therefore, we have the following from [98]:
E(Q
r
)=tr(F
H
?
Q
H
r
?
Q
r
F)
?1
(4.3.1)
=tr(F
H
?
Q
H
r
?
Q
r
F)
?1
F
H
F (4.3.2)
=trF(F
H
?
Q
H
r
?
Q
r
F)
?1
F
H
(4.3.3)
75
To manipulate this matrix, we use the following theorem.
Theorem 4.3.1 The matrix B = F(F
H
?
Q
H
r
?
Q
r
F)
?1
F
H
can be represented in the follow-
ing form:
B =[H
1
H
2
...H
n
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
(4.3.4)
where H
i
,i=1,...,n,are convolution matrices.
The proof is in Appendix B. Let Z = F
H
?
Q
H
r
?
Q
r
F and Y be the matrix that results after
Q
j
is deleted. Then Y = Z ? F
H
Q
H
j
Q
j
F. By using the matrix inversion formula, we
have
Y
?1
= Z
?1
+Z
?1
F
H
Q
H
j
(I ?Q
j
FZ
?1
F
H
Q
H
j
)
?1
Q
j
FZ
?1
(4.3.5)
FY
?1
F
H
= FZ
?1
F
H
+FZ
?1
F
H
Q
H
j
(I ?Q
j
FZ
?1
F
H
Q
H
j
)
?1
Q
j
FZ
?1
F
H
(4.3.6)
Suppose that we delete the array corresponding to Q
j
from Q
r
so that B becomes
?
B and Q
r
becomes
?
Q
r
. Then from (4.3.6) and the definition of B
?
B = B +BQ
H
j
(I ?Q
j
BQ
H
j
)
?1
Q
j
B
76
= B +BQ
H
j
G
j
Q
H
j
B (4.3.7)
where G
j
=(I?Q
j
BQ
H
j
)
?1
. This matrix inverse can be evaluated e?ciently as a result
of the following theorem:
Theorem 4.3.2 Assume that the {H
i
} are circulant matrices. Then G
j
is a circulant
matrix.
The proof is in Appendix B.
Because G
j
is the inverse of a circulant matrix, it can be computed e?ciently using
FFT?s.
Using the fact that B is symmetric, Q
j
Q
H
j
= I and Q
i
Q
H
j
=0withi negationslash= j,the
updating term of Eq. (4.3.7) becomes
?
j
= BQ
H
j
G
j
Q
j
B
H
= H
j
Q
H
j
Q
j
Q
H
j
G
j
Q
j
Q
H
j
Q
j
H
H
j
= H
j
Q
H
j
G
j
Q
j
H
H
j
(4.3.8)
From the property tr(ABC)=tr(BCA), we have
tr(?
j
)=tr(G
j
Q
j
H
H
j
H
j
Q
H
j
) (4.3.9)
In the form of (4.3.9), the criterion increment can be e?ciently evaluated for all
remaining j, since the matrix in the right-hand trace can be evaluated much more e?-
ciently than can ?
j
. The array j that gives the minimum value of (4.3.9) will cause the
77
minimum SSE increase and therefore should be deleted. Only after we determine the
optimal j are we required to compute ?
j
.
Once the optimal j is chosen, we must update the convolution kernels h
i
,i =
1,2,...,n, that construct the convolution matrices H
i
,i=1,2,...,n, so that they are
available for the next elimination step. First, we need initial values for the {H
i
} to begin
the SBS process. If we begin with all arrays selected, Q
H
r
Q
r
= I, and we have
B = F(F
H
F)
?1
F
H
(4.3.10)
= FF
H
(4.3.11)
Thus, the initial B represents an operation in the frequency domain that is equivalent
to zeroing the spatial-domain signal outside the ROS. Consequently, the initial {H
i
}
matrices are chosen to be the equivalent convolution operation with the ROS as the
spatial response. These matrices are stored and updated as impulse responses rather
than in their full matrix representation.
Because both B and
?
B in (4.3.7) are ofthe form(4.3.4), theupdateterm?
j
=
?
B?B
is also of the same form. Thus, the convolution matrices that define
?
B can be updated
by updating H
i
, i =1,...,n. The form of the update is given by
?
j
=[?H
1,j
?H
2,j
... ?H
n,j
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
(4.3.12)
78
More specifically, these update terms are given by
?
j
= BQ
H
j
G
j
Q
j
B
H
= H
j
Q
H
j
G
j
Q
j
[H
1
H
2
... H
n
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
(4.3.13)
Thus,
?
B =[H
1
+?H
1,j
H
2
+?H
2,j
... H
n
+?H
n,j
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
=[
?
H
1
?
H
2
...
?
H
n
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
(4.3.14)
Since the matrices {?H
i,j
} are circulant, we only need to calculate one column of each.
Let h
i
be the first column of H
i
. Then we use the following procedure to calculate the
update:
For i =1,...,n:
1. r = Q
j
h
i
. This is accomplished by downsampling h
i
with o?set j ?1.
79
2. s = G
j
r. This is a filtering operation using FFT?s.
3. t = Q
H
j
s. This is an upsampling operation with o?set j ?1.
4. ?h
i
= H
j
t. This is a filtering operation using FFT?s.
Thus, we can implement the algorithm with convolution kernels rather than huge
non-sparse matrices. We can e?ciently evaluate the criterion increase with (4.3.9) and
update the convolution kernels with (4.3.12) and (4.3.14) using only convolution, down-
sampling, and upsampling operations. The above updating process can be continued
recursively until the minimal number of arrays remains. The proposed method dramat-
ically reduces storage and computational complexity.
4.4 Image reconstruction from nonuniform samples
Once the best array set ?y
r
is selected, the unknown spatial samples x (and thus
the entire spatial-domain image) can be reconstructed from ?y
r
in a least-squares sense.
However, the least-square solution in (4.2.5) involves inversion of a huge matrix. There-
fore, directly using (4.2.5) to reconstruct the image is probably prohibitive. However,
if we store the resulting impulse responses from the SBS process that correspond to
the selected array o?sets, we can e?ciently compute the least-squares reconstruction as
follows:
hatwidey = F
H
(F
?
Q
H
r
?
Q
r
F
H
)
?1
F
?
Q
H
r
?y
r
=
?
B
?
Q
H
r
?y
r
80
=[
?
H
1
?
H
2
...
?
H
n
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
?y
r
(4.4.1)
where
?
H
1
?
H
2
...
?
H
n
has been computed with the e?cient algorithm in Section 4.3.
Then the spatial samples can be reconstructed as follows:
hatwidex = Fhatwidey (4.4.2)
With (4.4.1) and (4.4.2), we can e?ciently reconstruct the whole spatial image with
convolutions.
4.5 Experiments
Figure 4.4(a) shows a circular region of support (ROS) containing 3205 unknowns in
a64?64 grid. We experimented with three cases for this ROS, in which the number of
array elements decreased from case to case but the number of possible dithered positions
increased, so that the number of potential sample locations was fixed. In each case we
began with all possible array o?sets and sequentially eliminated one array o?set at a
time.
Table 4.1 lists the simulation results for these three cases. In the first case, the
criterion matrix was undefinedafter three arrays were eliminated. (The inverse in (4.3.3)
does not exist.) To keep the criterion defined, we needed to keep 14 of the 16?16 arrays,
corresponding to 3584 samples. In the second and third cases, 52 of the 8?8 arrays and
81
Table 4.1: Comparison of the criterion in three cases
Case 1 # of arrays 16 15 14 13
16x16 criterion 3205 4244 5612 undefined
Case 2 # of arrays 64 60 56 52
8x8 criterion 3205 4178 5357 8542
Case 3 # of arrays 256 240 224 208
4x4 criterion 3205 4171 5358 7658
Case 4 # of arrays 64 60 56 52
8x8 criterion 3779 4202 4990 8445
#ofsamples 4096 3840 3584 3328
208 of the 4?4 arrays were necessary to keep the criterion defined, corresponding to 3328
samples. The latter sample pattern yielded a lower error criterion than the former with
the same total number of samples. Figure 4.5(a) shows the same ROS in a 128?128 grid.
This allowed us to consider array o?sets spaced more closely together. We experimented
with one case (Case 4) for this ROS. In comparison with Case 2, the array size is the
same but there are 256 possible array locations in Case 4. The same number of arrays
and samples were necessary to keep the criterion defined, but the criterion achieved was
lower. The comparison of Case 2 and Case 4 suggests that with the same array size, the
more possible o?sets the better the potential performance. However, this seems to be a
function of where the SBS procedure is stopped, which is evident from a comparison of
the various intermediate stopping points.
Figure 4.4(b), 4.4(c), and 4.4(d) show the selected sample patterns for Case 1, Case
2, and Case 3. Figure 4.5(b) shows the selected sample pattern for Case 4.
82
10 20 30 40 50 60
10
20
30
40
50
60
10 20 30 40 50 60
10
20
30
40
50
60
(a) (b)
10 20 30 40 50 60
10
20
30
40
50
60
10 20 30 40 50 60
10
20
30
40
50
60
(c) (d)
Figure 4.4: (a) The 64?64 region of support (3205 unknowns) (b) 3584 selected samples
in the case 1 (c) 3328 selected samples in the case 2 (d) 3328 selected samples in the case
3 (black markers indicate selected samples)
83
20 40 60 80 100 120
20
40
60
80
100
120
20 40 60 80 100 120
20
40
60
80
100
120
(a) (b)
Figure 4.5: (a) The 128?128 region of support (white points indicate 3205 unknowns)
(b) 3328 selected samples in Case 4
4.6 Conclusion
We have developed an e?cient sampling pattern optimization technique that pro-
vides the most useful information for reconstructing images with a limited ROS in the
minimal necessary acquisition time. The proposed algorithm has the following advan-
tages: 1) it provides an e?cient update formula for the criterion based on convolution
kernels rather than large, non-sparse matrices; 2) it provides an e?cient update formula
for the convolution kernels based on the deleted array; 3) it provides an e?cient re-
construction method based on convolution. The proposed method dramatically reduces
storage and computational complexity. We found that the lowest SSE is obtained when
the number of possible dithered positions is largest, given the same total number of
samples.
84
We have not accounted for some practical concerns. That is, we have assumed
throughout that all filtering operations are circular convolutions. Furthermore, one must
take into account the fact that the array and the image have a finite size, and therefore
boundary e?ects become an issue.
Chapter 5
Sequential Forward Array Selection
5.1 Introduction
In this chapter, we develop an e?cient algorithm di?erent from that in Chapter
4 for optimizing the dithering pattern so that the image can be reconstructed from a
periodic nonuniform set of samples as reliably as possible. In this algorithm, we begin
the development with the modified criterion from Chapter 2 and sequentially select the
least noisy sample recursively based on a given ROS until the desired number of arrays
is selected. To simplify the complicated block selection algorithm, we decompose the
block selection algorithm into a sum of single sample selection algorithms. To avoid
the singularity at the switching point, we further modify the criterion for extending the
selection processfromthenumberofunknownsmorethanthenumberofselected samples
to less than the number of selected samples. The advantage of this method compared
to the SFS of single samples is time saved from selecting one block each time. The
disadvantage is that it may be impossible to reach the minimum density. The advantage
compared to SBS of sample blocks is that it brings real-time implementation into the
realm of possibility because the previously selected sample block can be observed while
the next sample block is selected.
5.2 Notation
If we represent the discretized frequency image with size m?n as a column-ordered
vector y, the unknown discretized spatial support by a vector x, and the mapping from
85
86
the spatial-domain samples to the frequency-domain image by F, the fully sampled
frequency image can be expressed as
y = Fx+w (5.2.1)
where w is zero-mean, i.i.d. noise. Suppose the array size is m
1
? n
1
and the possible
dithered positions of the array are m
2
?n
2
.Wehavem = m
1
?m
2
and n = n
1
?n
2
.
If we sample the frequency image y using a single position of the sensor array with
o?set indexed by 1 ? i ? n, we can represent this by
y
i
= Q
i
y (5.2.2)
where Q
i
downsamples the fully sampled image and orders the result into a vector. Then
we can rewrite y in a rearranged form y
r
as
y
r
=[y
H
1
y
H
2
... y
H
n
]
H
(5.2.3)
=[Q
H
1
Q
H
2
... Q
H
L
1
]
H
(Fx+w)
= Q
r
Fx+w
r
where w
r
is the similarly rearranged version of w and Q
r
=[Q
H
1
Q
H
2
... Q
H
n
]
H
.Ifwe
choose a subset of k of the N shifted arrays, we obtain
?y
r
=[y
H
k
1
y
H
k
2
... y
H
k
k
]
H
(5.2.4)
=[Q
H
k
1
Q
H
k
2
... Q
H
k
k
]
H
(Fx+w)
87
=
?
Q
r
Fx+?w
r
where ?w
r
is the corresponding subset of w
r
and
?
Q
r
is the corresponding subset of Q
r
.
5.3 Criterion
Let A =
?
Q
r
F. In chapter 2, we have shown that the criteria is the trace of (AA
H
)
?1
when the number of selected samples is less than the number of unknowns in in the ROS
and the trace of (A
H
A)
?1
when the number of selected samples is more than the number
of unknowns in the ROS. To unite these two di?erent criteria into one for selection
purposes, we propose a modified criterion as follows
E(
?
Q
r
)=tr(AA
H
+?I)
?1
=tr(
?
Q
r
FF
H
?
Q
r
H
+?I)
?1
(5.3.1)
where ? is a very small number.
Let n
r
be the number of unknowns in the ROS. When m
1
n
1
k ? n
r
and A
k
A
H
k
is
nonsingular,
lim
??0
tr(A
k
A
H
k
+?I
k
)
?1
=tr(A
k
A
H
k
)
?1
(5.3.2)
Let the singular value decomposition of A
k
= U
k
?
k
V
H
k
.Then
tr(A
k
A
H
k
+?I)
?1
=tr(U
k
?
k
?
H
k
U
H
k
+?I
k
)
?1
=trU
H
k
(U
k
?
k
?
H
k
U
H
k
+?I
k
)
?1
U
k
88
=tr(?
k
?
H
k
+?U
H
k
U
k
)
?1
=tr(?
k
?
H
k
+?I
k
)
?1
(5.3.3)
When m
1
n
1
k ? n
r
,?
k
has only n
r
nonzero elements.
tr(?
k
?
H
k
+?I
k
)
?1
=
n
r
summationdisplay
i=1
1
|?
i
|
2
+?
+(m
1
n
1
k?n
r
)
1
?
(5.3.4)
tr(A
H
k
A
k
)
?1
= lim
??0
n
r
summationdisplay
i=1
1
|?
i
|
2
+?
= lim
??0
tr(A
k
A
H
k
+?I)
?1
?(m
1
n
1
k ?n
r
)
1
?
(5.3.5)
Since the second term in 5.3.5 is a constant at each step, it has no e?ect on the
selection procedure. We only need this term to evaluate the true criterion after the
selection. Therefore, 5.3.3 and 5.3.5 show the reasonability of the selection criterion in
5.3.1.
5.4 E?cient recursive algorithm
To simplify the computational method, we first transform the block matrix in 5.3.1
into a diagonal block matrix after some modulation.
AA
H
=
?
QFF
H
?
Q
H
=[Q
H
n
1
FF
H
Q
H
n
2
???Q
H
n
k
]
H
FF
H
[Q
n
1
Q
n
2
???Q
n
k
]
89
=
?
?
?
?
?
?
?
?
?
?
?
?
Q
n
1
FF
H
Q
H
n
1
Q
n
1
FF
H
Q
H
n
2
??? Q
n
1
FF
H
Q
H
n
k
Q
n
2
FF
H
Q
H
n
1
Q
n
2
FF
H
Q
H
n
2
??? Q
n
2
FF
H
Q
H
n
k
.
.
.
.
.
. ???
.
.
.
Q
n
k
FF
H
Q
H
n
1
Q
n
k
FF
H
Q
H
n
2
??? Q
n
k
FF
H
Q
H
n
k
?
?
?
?
?
?
?
?
?
?
?
?
(5.4.1)
Q
i
FF
H
Q
H
j
,i,j=1,2,???,n, are easily seen to be circulant matrices, since FF
H
is circu-
lant and Q
i
FF
H
Q
i
are the same for i =1,2,???,nbecause FF
H
is a circulant matrix.
Therefore, AA
H
is conjugate symmetric matrix with circulant blocks [120] (circulant-
block-circulant blocks in the 2-D case).
Let F
d
represent the downsampled Fourier matrix corresponding to the dithered ar-
ray. ThenF
H
d
F
d
= I. F
d
Q
i
FF
H
Q
j
F
H
d
diagonalizesthecirculantmatrixQ
i
FF
H
Q
j
,i,j=
1,2,???,n[120]. Let diag(d
ij1
,d
ij2
,???,d
ijg
)=F
d
Q
i
FF
H
Q
j
F
H
d
. Finally, let F
D
be a
block-diagonal matrix formed by replicating F
d
along the diagonal k times,
(AA
H
)
?1
= F
D
F
H
D
(AA
H
)
?1
F
D
F
H
D
= F
D
(F
D
AA
H
F
H
D
)
?1
F
H
D
= F
D
?
?
?
?
?
?
?
?
?
?
?
?
F
D
Q
n
1
FF
H
Q
H
n
1
F
D
d
H
F
D
Q
n
1
FF
H
Q
H
n
2
F
H
D
??? F
D
Q
n
1
FF
H
Q
H
n
k
F
H
D
F
D
Q
n
2
FF
H
Q
H
n
1
F
H
D
F
D
Q
n
2
FF
H
Q
H
n
2
F
H
D
??? F
D
Q
n
2
FF
H
Q
H
n
n
FD
H
.
.
.
.
.
. ???
.
.
.
F
D
Q
n
k
FF
H
Q
H
n
1
F
H
D
F
D
Q
n
k
FF
H
Q
H
n
2
FD
H
??? F
D
Q
n
k
FF
H
Q
H
n
k
F
H
D
?
?
?
?
?
?
?
?
?
?
?
?
?1
F
H
D
90
= F
D
?
?
?
?
?
?
?
?
?
?
?
?
diag(d
111
???d
11m
) diag(d
121
???d
12m
) ??? diag(d
1n1
d
1n2
???d
1nm
)
diag(d
211
???d
21m
) diag(d
221
???d
22m
) ??? diag(d
2n1
???d
2nm
)
.
.
.
.
.
. ???
.
.
.
diag(d
k11
???d
k1m
) diag(d
k21
???d
k2m
) ??? diag(d
nn1
???d
nnm
)
?
?
?
?
?
?
?
?
?
?
?
?
?1
F
H
D
= F
D
?
?
?
?
?
?
?
?
?
?
?
?
diag(
hatwide
d
111
???
hatwide
d
11m
) diag(
hatwide
d
121
???
hatwide
d
12m
) ??? diag(
hatwide
d
1n1
hatwide
d
1n2
???
hatwide
d
1nm
)
diag(
hatwide
d
211
???
hatwide
d
21m
) diag(
hatwide
d
221
???
hatwide
d
22m
) ??? diag(
hatwide
d
2n1
???
hatwide
d
2nm
)
.
.
.
.
.
. ???
.
.
.
diag(hatd
k11
???
hatwide
d
k1m
) diag(
hatwide
d
k21
???
hatwide
d
k2m
) ??? diag(
hatwide
d
nn1
???
hatwide
d
nnm
)
?
?
?
?
?
?
?
?
?
?
?
?
F
H
D
(5.4.2)
where
D
?1
g
=
?
?
?
?
?
?
?
?
?
?
?
?
d
11g
d
12g
??? d
1kg
d
21g
d
22g
??? d
2kg
.
.
.
.
.
. ???
.
.
.
d
k1g
d
k2g
??? d
kkg
?
?
?
?
?
?
?
?
?
?
?
?
?1
=
?
?
?
?
?
?
?
?
?
?
?
?
hatwide
d
11g
hatwide
d
12g
???
hatwide
d
1kg
hatwide
d
21g
hatwide
d
22g
???
hatwide
d
2kg
.
.
.
.
.
. ???
.
.
.
hatwide
d
k1g
hatwide
d
k2g
???
hatwide
d
kkg
?
?
?
?
?
?
?
?
?
?
?
?
(5.4.3)
where g =1,2,???,m. To calculate the vector d
ij
, we do the following steps, which
involves only FFT?s instead of matrix operations [120,121].
91
1. Take the FFT of the ROS image to get the impulse response h,whichisthefirst
column of FF
H
. FF
H
is a circulant matrix. Only the first column needs to be
stored.
2. Shift h by j, then downsample the shifted h by Q
i
to get the first column h
d
of the
circulant matrix Q
i
FF
H
Q
j
,i,j=1,2,???,n.
3. TaketheFFTofh
d
togetthediagonalvector ofthediagonalmatrixF
d
Q
i
FF
H
Q
j
F
H
d
.
By using the property tr(ABC)=tr(CAB), we have
tr(AA
H
)
?1
=tr
?
?
?
?
?
?
?
?
?
?
?
?
F
H
F
d
?
?
?
?
?
?
?
?
?
?
?
?
diag(
hatwide
d
111
hatwide
d
112
???
hatwide
d
11m
) ??? diag(
hatwide
d
1k1
hatwide
d
1k2
???
hatwide
d
1km
)
diag(
hatwide
d
211
hatwide
d
H
212
???
hatwide
d
H
21m
) ??? diag(
hatwide
d
2k1
hatwide
d
2k2
???
hatwide
d
2km
)
.
.
. ???
.
.
.
diag(
hatwide
d
H
k11
hatwide
d
H
k12
???
hatwide
d
H
k1m
) ??? diag(
hatwide
d
kk1
hatwide
d
kk2
???
hatwide
d
kkm
)
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
=tr
?
?
?
?
?
?
?
?
?
?
?
?
diag(
hatwide
d
111
hatwide
d
112
???
hatwide
d
11m
) ??? diag(
hatwide
d
1k1
hatwide
d
1k2
???
hatwide
d
1km
)
diag(
hatwide
d
211
hatwide
d
212
???
hatwide
d
21m
) ??? (
hatwide
d
2k1
hatwide
d
2k2
???
hatwide
d
2km
)
.
.
.
.
.
. ???
.
.
.
diag(
hatwide
d
k11
hatwide
d
k12
???
hatwide
d
k1m
) ??? diag(
hatwide
d
kk1
hatwide
d
kk2
???
hatwide
d
kkm
)
?
?
?
?
?
?
?
?
?
?
?
?
(5.4.4)
Let
?
A represent the new A after the ith array is added to A.
?
A
?
A
H
=[Q
H
Q
H
i
]
H
FF
H
[QQ
i
]
=
?
?
?
?
QFF
H
Q
H
QFF
H
Q
H
i
Q
H
i
FF
H
Q
H
Q
i
FF
H
Q
H
i
?
?
?
?
(5.4.5)
92
Let
?
D
g
represent the corresponding new D
g
after adding the ith array.
?
D
g
=
?
?
?
?
D
g
S
H
i
g
S
i
g
T
i
g
?
?
?
?
(5.4.6)
With Lemma 2.3.5 and 2.3.11 in Chapter 2, we have
?
D
?1
g
=
?
?
?
?
P
?1
+P
?1
S
H
i
?
i
S
i
p
?1
?P
?1
S
H
i
??
i
S
i
P
?1
?
?1
i
?
?
?
?
(5.4.7)
where ?
i
= T
i
?S
i
P
?1
Q
i
is a scalar, T
i
are the same for i =1,2,???,mand P = D
g
.
tr(
?
D
?1
g
)=tr(P
?1
)+
1+S
i
P
?2
S
H
i
T
i
?S
i
P
?1
R
i
(5.4.8)
Let ?
g
represent the second term in (5.4.8) and ? represent the criterion increment
for each step. Then
?
g
=
1+S
i
P
?2
S
H
i
T
i
?S
i
P
?1
R
i
(5.4.9)
? =
m
summationdisplay
g=1
?
g
(5.4.10)
where g =1,2,???,m.
We select the array to minimize ? each step.
5.5 E?cient computational method
To find the minimizer of ? in 5.4.10, which we denote as the ith array, we need to
calculate ?
g
,g=1,2,???,nfor every possible j =1,2,???,n. That is, we must calculate
93
S
j
P
?1
S
H
j
and S
j
P
?2
S
H
j
for every possible j =1,2,???,n. The computational burden
will be huge. To reduce the computational cost, we derived recursive formulas for the
above two terms such that the updated parts are calculated based only on the previous
selected array.
Let
?
S
j
represent the new S
j
after the ith array added. Note that the elements of
?
S
j
come from S
j
, but
?
S
j
only includes the elements corresponding to the o?sets already
selected, s
ji
.Let
?
P =
?
D
g
, ?
i
= H
i
= T
i
?S
i
P
?1
Q
i
, R
i
= P
?1
S
H
i
.Then
?
S
j
?
P
?1
?
S
H
j
=[S
j
s
ji
]
?
?
?
?
P
?1
+
1
?
i
P
?1
S
H
i
S
i
P
?1
?P
?1
S
H
i
?
1
?
i
S
i
P
?1 1
?
i
?
?
?
?
[S
j
s
ji
]
H
= S
j
P
?1
S
H
j
+S
j
P
?1
S
H
i
H
?1
i
S
i
P
?1
S
H
j
?S
j
P
?1
S
H
i
H
?1
i
s
H
ji
+s
ji
H
?1
i
s
ji
= S
j
P
?1
S
H
j
+
1
?
i
(S
j
R
i
)(S
j
R
i
)
H
?
1
?
i
s
ji
(S
j
R
i
)
H
?
1
?
i
(S
j
R
i
)s
H
ji
+
1
?
i
s
ji
s
H
ji
= S
j
P
?1
S
H
j
+
1
?
i
(||S
j
R
i
||?2Re((S
j
R
i
)s
H
ji
)+|s
ji
|) (5.5.1)
Let F
i
= P
?1
R
i
.
?
S
j
?
P
?2
?
S
H
j
=[S
j
s
ji
]
?
?
?
?
P
?1
+
1
?
i
P
?1
S
H
i
S
i
P
?1
?P
?1
S
H
i
?
1
?
i
S
i
P
?1 1
?
i
?
?
?
?
2
[S
j
s
ji
]
H
= S
j
P
?2
S
H
j
+S
j
P
?2
S
H
i
H
?1
i
S
i
P
?1
S
H
j
?S
j
P
?2
S
H
i
H
?1
i
s
H
ji
+S
j
P
?1
S
H
i
H
?1
i
S
i
P
?2
S
H
j
+S
j
P
?1
S
H
i
H
?1
i
S
i
P
?2
S
H
i
H
?1
i
P
?1
S
H
j
?S
j
P
?1
S
H
i
H
?1
i
S
i
P
?2
S
H
i
H
?1
i
s
j
i
H
?s
ji
H
?1
i
S
i
P
?2
S
H
j
94
?s
ji
H
?1
i
S
i
P
?2
S
H
i
H
?1
i
S
i
P
?1
S
H
j
+s
ji
H
?1
i
S
i
P
?2
S
H
i
H
?1
i
s
H
ji
+S
j
P
?1
S
H
i
H
?2
i
S
i
P
?1
S
H
j
?S
j
P
?1
S
H
i
H
?2
i
s
H
ji
?s
ji
H
?2
i
S
i
P
?1
S
H
j
+s
ji
H
?2
i
s
H
ji
= S
j
P
?2
S
H
j
+
2
?
i
Re(S
j
R
i
(S
j
F
i
)
H
)?(
2
?
i
((S
j
F
i
)s
H
ji
)+S
j
F
i
s
H
ij
)
+
1+||R
i
||
2
?
2
i
||S
j
R
i
||
2
+
1+||R
i
||
2
?
2
i
|s
ji
|
2
?
1+||R
i
||
2
?
2
i
2Re((S
j
R
i
)s
H
ji
)
= S
j
P
?2
S
H
j
+
2
?
i
(Re(S
j
R
i
(S
j
F
i
)
H
)?(S
j
F
i
)s
H
ji
)
+
1+||R
i
||
2
?
2
i
(||S
j
R
i
||
2
?2Re((S
j
R
i
)s
H
ji
)+|s
ji
|
2
) (5.5.2)
With (5.5.1) and (5.5.2), the inverses only operate on the previous selection array
rather than all possible arrays. We use Gaussian elimination to solve for R
i
and F
i
.The
computational cost is drastically reduced.
5.6 Experiments
Figure 5.1(a) shows a circular region of support (ROS) containing 3205 unknowns in
a64?64 grid. We experimented with three cases for this ROS, in which the number of
array elements decreased from case to case but the number of possible dithered positions
increased, so that the number of potential sample locations was fixed. In each case we
began with all possible array o?sets and sequentially eliminated one array o?set at a
time.
In Case 1, the array size is 4?4, and the number of possible dithered positions is
16 ? 16. In Case 2, the number of dithered positions is 8 ? 8,thearraysizeis8? 8.
In Case 3, the number of possible dithered positions is 16?16, the array size is 4 ?4.
Figure 5.1(b), (c), (d) shows the sampling pattern with 13 arrays selected in Case 1,
95
52 arrays selected in Case 2, and 208 arrays selected in Case 3 (corresponding to 3328
points).
To make a comparison with the same array size but di?erent number of potential
dithered positions we considered Case 4. In Case 4, the array size is 8?8, the same as
that in Case 2. However, the number of possible dithered positions is 16?16, di?erent
than that in Case 2. In Case 4, we set ? =10
?2
to reduce accumulated error due to
ill-conditioning. Figure 5.2(a) shows a 128 ? 128 ROS with 3205 unknowns. Figure
5.2(b) shows the selection pattern with 52 arrays (3328 points) selected.
Table 5.1 shows the criterion values at various steps for the four cases. In the
experiments, we used ? =10
?4
to avoid the singularity in the selection process. When
all 4096 samples are selected, the criterion values for Case 1, Case 2, and Case 3 are
approximately equal to 891/?. This is consistent with the analysis. Let A = FS,Sis
diagonal matrix with 3205 diagonal elements corresponding to nonzero unknowns as 1
and 891 diagonal elements corresponding to zero locations as 0.
tr(AA
H
+?I)
?1
=(FSF
H
+?FF
H
)
?1
= F(S +?I)
?1
F
H
(5.6.1)
tr(AA
H
+?I)
?1
=tr(F
H
F(S +?I)
?1
)
=tr(S +?I)
?1
=
891
?
+
3205
1+?
? 8910000 +3205
96
10 20 30 40 50 60
10
20
30
40
50
60
10 20 30 40 50 60
10
20
30
40
50
60
(a) (b)
10 20 30 40 50 60
10
20
30
40
50
60
10 20 30 40 50 60
10
20
30
40
50
60
(c) (d)
Figure 5.1: (a) The 64?64 region of support (3205 unknowns) (b) 3328 selected samples
in Case 1 (c) 3328 selected samples in Case 2 (d) 3328 selected samples in Case 3 (black
markers indicate selected samples)
97
20 40 60 80 100 120
20
40
60
80
100
120
20 40 60 80 100 120
20
40
60
80
100
120
(a) (b)
20 40 60 80 100 120
20
40
60
80
100
120
(c)
Figure 5.2: (a) The 128?128 region of support (white points indicate 3205 unknowns)
(b) 3328 selected samples in Case 4 (c) 4096 selected samples in Case 4.
98
Table 5.1: Comparison of the criterion in four cases
Case 1 # of arrays 9 13 14 15 16
16x16 criterion 0.950431 1839970 3789951 6349934 8909999
Case 2 # of arrays 36 52 56 60 64
8x8 criterion 0.892832 1229985 3789957 6349921 8909938
Case 3 # of arrays 144 208 224 240 256
4x4 criterion 0.884639 1230000 3789985 6349982 8910071
Case 4 # of arrays 36 52 56 60 64
8x8 criterion 0.840235 12302.2 37901.3 63501.1 89100.9
#ofsamples 2340 3328 3584 3840 4096
Table 5.2: Comparison of the criterion in four cases
Case 1 #ofarrays 9 13 14 15 16
16x16 SFS time (s) 4.25 6.85 7.59 8.34 9.11
SBS time (s) N/A N/A 6.71 3.48 0
Case 2 #ofarrays 36 52 56 60 64
8x8 SFS time (s) 9.19 20.92 25.15 30.0 35.69
SBS time (s) N/A 1111.9 766 398 0
Case 3 #ofarrays 144 208 224 240 256
4x4 SFS time (s) 87.34 291.91 367.37 477.89 600.17
Case 4 #ofarrays 36 52 56 60 64
8x8 SFS time (s) 21.94 42.60 49.46 57.07 65.65
#ofsamples 2340 3328 3584 3840 4096
? 891/? (5.6.2)
Comparing the criterion values at 3328 samples ( 13 arrays in Case 1, 52 arrays
in Case 2 and Case 4 and 208 arrays in Case 3), the values are close for Case 2, 3,
and 4 (adjusting Case 4 by a factor of 10
2
to compensate for the di?erent choice of ?).
However, the criterion value for Case 1 is much higher. Taking into account the time
cost, the experimental result in Case 2 is desired. It takes 21 seconds to select 52 arrays
99
and half a minute to finish the whole selection process in Case 2. For the larger field
of view (FOV) in Case 4, it takes 42.60 seconds to select 52 arrays, 65.65 seconds to
select 64 arrays (corresponding to 4096 samples), and 1 hour 9 minutes 16 seconds to
finish the whole selection. Table 5.2 shows the time cost for the four cases and sequential
backward selection for Case 1 and Case 2. In comparison with the sequential backward
array selection, sequential forward array selection is of more practical value because it
can get the desired selection in a minute.
Figure 5.3(a) and (b) show the logarithmitic scale of criterion curve in Case 2 and
Case 4. We observed that the criterion has a sharp increase because of the singularity
of AA
H
. With a few more samples than the unknowns selected, the criterion curve
eventually becomes smooth. The criterion curve eventually increases uniformly because
adding the value ? makes the selection increment nearly uniform for the last phase of
the selection process. The increase due to the 1/? term dominates the increase.
5.7 Conclusion
We developed an e?cient sequential forward array selection technique that provides
the most useful information for restoring and superresolving bandlimited images in a
minute. We developed a modified criterion to unite the whole selection process. We
derived an e?cient recursive selection algorithm for the update criterion. We derived
a recursive algorithm for e?ciently computing inverse related terms such that they can
be updated based only on the previous selected array. We use Gaussian elimination to
find a set of small inverse solutions at each step. The proposed method has important
potential in practice because it can select a reduced set of the most useful samples in a
100
0 10 20 30 40 50 60 70
10
?2
10
?1
10
0
10
1
10
2
10
3
10
4
10
5
10
6
10
7
# of selected arrays
(a)
0 50 100 150 200 250 300
10
?2
10
?1
10
0
10
1
10
2
10
3
10
4
10
5
10
6
10
7
# of selected arrays
(b)
Figure 5.3: (a) The semilog criterion curve in Case 2 (b) The semilog criterion curve in
Case 4
101
minute and can reconstruct the image with fewer samples than regular sampling. This
may reduce the image acquisition time greatly and is very important in applications to
MRSI, PMMW imaging, and multiband digital signal reception.
Chapter 6
Conclusion and Discussion
Four e?cient optimal magnetic resonance observation selection schemes have been
proposed to overcome the problems in those existing methods, reduce observation time,
and improve the image resolution.
6.1 E?cient sequential forward selection
In Chapter 2, we proposed an e?cient SFS optimization scheme, which has the
following advantages: 1) it overcomes the limitation of the SFS algorithm based on the
standard MSE criterion in that it can start with an empty initial matrix; 2) it drastically
reducesstoragespaceandcomputationaltimewiththeproposedSFSrecursivealgorithm
and the conjugate gradients algorithm; 3) it is more applicable in practice than SBS in
that the best sample can be observed as soon as it is selected, making possible real-time
selection and acquisition. Simulation with real MR data shows a good-quality image
can be reconstructed with only a few more samples than the minimum density by our
method.
Further work is clearly necessary to accelerate the algorithm if it is to be useful in
a real-time context on an ordinary computer. However, the algorithm described here
makes the possibility of real-time implementation into the realm of feasibility appear
promising.
102
103
6.2 E?cient backward selection of hexagonal samples
In Chapter 3, we proposed a new method for e?ciently choosing MRSI samples
by applying SBS on hexagonal samples. We computed a generalized MDFT with Smith
formdecomposition andoutputrectangular image pixelsdirectlyfromhexagonal samples
without interpolation. A rate lower than the Nyquist rate on a hexagonal grid can be
achieved with the proposed method. Furthermore, SBS on a hexagonal sampling grid
takes less selection time than SBS on a rectangular sampling grid,
6.3 E?cient selection of sample block
In Chapter 4, we developed a method for selecting a periodic array of samples. It
used an e?cient update formula for the criterion based on convolution kernels rather
than large, non-sparse matrices, an e?cient update formula for the convolution kernels
based on the deleted array, and an e?cient reconstruction method based on convolution
and FFT?s. The proposed method dramatically reduces storage and computational com-
plexity. We found that the lowest SSE is obtained when the number of possible dithered
positions is largest, given the same total number of samples.
6.4 E?cient sequential forward selection of array
In Chapter 5, we developed an e?cient sequential forward array selection technique
thatrapidlyselectsareducedsetofthemostimportantsamplesforreconstructingimages
with a limited ROS. We developed a modified criterion applicable to the whole selection
process such that the selection smoothly switches from the case with the number of
samples less than the number of unknowns to the case with the number of samples more
104
than the number of unknowns. We derived an e?cient recursive selection algorithm for
the update criterion. We derived an e?cient recursive algorithm for computing inverse-
related terms such that they can be updated based only on the previous selected array.
We use Gaussian elimination to solve for a set of small matrix inversions at each step.
With this e?cient strategy, the proposed method can finish the selection process within
half a minute for a practical size of image and can reconstruct the image with fewer
samples than regular sampling. This advantage may reduce the image acquisition time
greatly and make it very promising in applications such as MRSI, PMMW imaging, and
multiband digital signal reception.
6.5 Future work
1. A remaining problem is to derive an upper bound on SBS for array selection. The
upper bound for SBS for single samples was derived by Reeves and Zhao [100].
Deriving the upper bound for SBS for array selection appears to be a much more
challenging problem.
2. Another problem remaining is to develop an algorithm for beginning with a stan-
dard ROS and adjusting the selection for a modified ROS. To decrease the selection
time for the SFS algorithm, it is possible to assume several standard ROS?s, choose
the ROS closest to the actual head as the standard ROS, and start selection of sam-
ples from the standard ROS selection.
Appendix A
We prove that circular shifts of the selected samples and the ROS yield the same
value of the MSE criteria. For ease of notation, we consider the 1-D case. The proof for
the 2-D case is exactly the same, although the notation is more involved.
Theorem 6.5.1 Let {s
i
} and {hatwides
i
} be two sets of size m such that hatwides
i
=((s
i
+l))
M
for
arbitrary l,where((n))
M
= n mod M. Also, let {r
i
} and {hatwider
i
} be two sets of size n
such that hatwider
i
=((r
i
+ k))
M
for arbitrary k.LetA and A
l,k
have elements defined as
a(i,j)=e
?j2?
s
i
r
j
M
and hatwidea(i,j)=e
?j2?
hatwides
i
hatwider
j
M
.Then,
tr(A
H
A)
?1
= tr(A
H
l,k
A
l,k
)
?1
if m ? n and (A
H
A)
?1
exists,
and
tr(AA
H
)
?1
= tr(A
l,k
A
H
l,k
)
?1
if m ? n and (AA
H
)
?1
exists.
Proof: The individual terms in A
l,k
are of the form,
e
?j2?
hatwides
i
hatwider
j
M
= e
?j2?
s
i
r
j
M
e
?j2?
lr
j
M
e
?j2?
ks
i
M
e
?j2?
kl
M
(6.5.1)
105
106
Therefore,
A
l,k
= e
?j2?
kl
M
?
?
?
?
?
?
?
?
?
?
?
?
e
?j2?
ks
1
M
0 ??? 0
0 e
?j2?
ks
2
M
??? 0
.
.
.
.
.
.
.
.
.
.
.
.
00??? e
?j2?
ks
m
M
?
?
?
?
?
?
?
?
?
?
?
?
A
?
?
?
?
?
?
?
?
?
?
?
?
e
?j2?
lr
1
M
0 ??? 0
0 e
?j2?
lr
2
M
??? 0
.
.
.
.
.
.
.
.
.
.
.
.
00??? e
?j2?
lr
n
M
?
?
?
?
?
?
?
?
?
?
?
?
= ?H
s
AH
r
(6.5.2)
For convenience, we have denoted the left-hand matrix above as H
s
and the right-
hand matrix as H
r
,and? = e
?j2?
kl
M
. It is apparent that both matrices are unitary,
which will be useful below.
If m ? n,
tr(A
H
l,k
A
l,k
)
?1
=tr(?H
s
AH
r
)
H
?H
s
AH
r
)
?1
= ??
?
tr(H
?1
r
(A
H
A)
?1
H
r
)
=tr(H
r
H
?1
r
(A
H
A)
?1
)
=tr(A
H
A)
?1
(6.5.3)
If m ? n,
tr(A
l,k
A
H
l,k
)
?1
=tr(?H
s
AH
r
(?H
s
AH
r
)
H
)
?1
= ??
?
tr(H
s
(AA
H
)
?1
H
?1
s
)
=tr(H
?1
s
H
s
(AA
H
)
?1
)
=tr(A
H
A)
?1
(6.5.4)
Appendix B
Proof of Theorem 4.3.1:
For simplicity, we present a proof for the 1-D case; the 2-D case is conceptually identical
but requires more complex notation. Define
?
S
j
as a matrix implementing a k-space shift
and S
j
as the corresponding spatial-domain matrix operator made up of linear-phase
weights; that is,
?
S
j
F = FS
j
(6.5.5)
We observe that
?
S
H
j
?
S
j
= I and S
H
j
S
j
= I.Leti = j mod N,andk = Nfloorleftj/Nfloorright,whereN
is the spacing between samples in the sensor array; thus, j = i+k. Define e
j
as the unit
sample vector shifted by j?1 samples. Then the jth column of F
H
(F
?
Q
H
r
?
Q
r
F
H
)
?1
F is
given by
F(F
H
?
Q
H
r
?
Q
r
F)
?1
F
H
e
j
= F(F
H
?
Q
H
r
?
Q
r
F)
?1
F
H
?
S
k
e
i
(6.5.6)
= FS
k
(S
H
k
F
H
?
Q
H
r
?
Q
r
S
k
FS
k
)
?1
S
H
k
F
H
?
S
k
e
i
(6.5.7)
=
?
S
k
F(F
H
?
S
H
k
?
Q
H
r
?
Q
r
?
S
k
F)
?1
Fe
i
(6.5.8)
Because
?
S
k
shifts the sensor array by an integer number of periods of the sensor array
pattern,
?
Q
r
S
k
=
?
Q
r
. In this case we have from (6.5.8)
F
H
(F
?
Q
H
r
?
Q
r
F
H
)
?1
Fe
i+k
=
?
S
k
F
H
(F
?
Q
H
r
?
Q
r
F
H
)
?1
Fe
i
(6.5.9)
107
108
Thus, the response to an impulse at location i+k is a copy shifted by k of the response
to an impulse at location i. That is, there is an impulse response associated with each
shifted position of the sensor array.
Using the matrix form in (4.3.4) and post-multiplying by e
j
,wehave
B =[H
1
H
2
...H
n
]
?
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
Q
1
Q
H
2
Q
2
.
.
.
Q
H
n
Q
n
?
?
?
?
?
?
?
?
?
?
?
?
e
j
(6.5.10)
= H
i+1
e
j
(6.5.11)
=
?
S
k
H
i+1
e
i
(6.5.12)
Thus, by defining a Toeplitz matrix H
i+1
such that F
H
(F
?
Q
H
r
?
Q
r
F
H
)
?1
Fe
i
= H
i+1
e
i
,we
obtain the structure in (4.3.4).
Proof of Theorem 4.3.2:
Observe that Q
j
Q
H
j
= I and Q
i
Q
H
j
=0fori negationslash= j.Then
Q
j
BQ
H
j
= Q
j
[H
1
H
2
...H
n
]
?
?
?
?
?
?
?
?
?
?
?
?
Q
H
1
(Q
1
Q
H
j
)
Q
H
2
(Q
2
Q
H
j
)
.
.
.
Q
H
n
(Q
n
Q
H
j
)
?
?
?
?
?
?
?
?
?
?
?
?
= Q
j
H
j
Q
H
j
(6.5.13)
Thus,
G
j
=(I ?Q
j
H
j
Q
H
j
)
?1
(6.5.14)
109
Using the fact that a circulant matrix with downsampled rows and columns is a circulant
matrix, Q
j
H
j
Q
H
j
corresponds to a circulant matrix. Furthermore, the sum of circulant
matrices and the inverse of a circulant matrix are still circulant matrices. Therefore, G
j
is a circulant matrix.
Bibliography
[1] Z. H. Cho, J. P. Jones, and M. Singh, Foundations of Medical Imaging. Wiley,
1993.
[2] T. F. Kirn, ?Magnetic resonance spectroscopy may hold promise in studying
metabolites, tissues,? Journal of the American Medical Association, vol. 261,
p. 1103?., 1989.
[3] J. H. Duyn, ?Spectroscopic imaging techniques.? ISMRM 1998 spectroscopy teach-
ing course, 1998.
[4] S. J. Nelson, D. B. Vigneron, and W. P. Dillon, ?Serial evaluation of patients with
brain tumors using volume MRI and 3D 1H MRSI,? NMR in Biomedicine, vol. 12,
pp. 123?138, May 1999.
[5]C.M.Segebarth,D.F.Baleriaux,P.R.Luyten,andJ.A.denHollander,?Detec-
tion of metabolic heterogeneity of human intracranial tumors in vivo by 1H NMR
spectroscopic imaging,? Magn Reson Med, vol. 13, pp. 62?76, Jan 1990.
[6] M. J. Fulham, A. Bizzi, M. J. Dietz, H. H. Shih, R. Raman, G. S. Sobering,
J. A. Frank, A. J. Dwyer, J. R. Alger, and G. D. Chiro, ?Mapping of brain tumor
metabolites with MR spectroscopic imaging: clinical relevance,? Magn Reson Med,
vol. 13, pp. 62?76, Jan. 1990.
[7]P.R.Luytenet al., ?Metabolic imaging of patients with intracranial tumors: H-1
MR spectroscopic imaging and PET,? Radiology, vol. 176, pp. 791?799, 1990.
[8] K. Herholz, W. Heindel, P. R. Luyten, J. A. D. Hollander, U. Pietrzyk, J. Voges,
H. Kugel, G. Friedman, and W. D. Heiss, ?In-vivo imaging of glucose consumption
and lactate concentration in human gliomas,? Ann. Neurol, vol. 31, p. 319, 1992.
[9] M. J. Fulham, A. Bizzi, M. J. Dietz, H. L. Shih, R. Raman, G. S. Sobering, J. A.
Frank, A. J. Dwyer, J. R. Alger, and G. D. Chiro, ?Mapping of brain tumor
metabolites with proton MR spectroscopic imaging: clinical relevance,? Radiology,
vol. 185, p. 675, 1992.
[10] J. H. Duijn, G. B. Matson, A. A. Maudsley, J. W. Hugg, and M. W. Weiner,
?Human brain infarction: proton MR spectroscopy,? Radiology, vol. 183, pp. 711?
8, June 1992.
110
111
[11] P. B. Barker, J. H. G. JH, P. C. M. van Zijl, B. J. Soher, D. F. Hanley, A. M.
Agildere, S. M. Oppenheimer, and R. N. Bryan, ?Acute stroke: evaluation with
serial proton MR spectroscopic imaging,? Radiology, vol. 183, pp. 711?718, 1992.
[12] R. I. Kuzniecky, J. W. Hugg, H. P. Hetherington, E. Burrerworth, E. Biler,
E. Faught, and F. Gilliam, ?Relative utility of 1H spectroscopic imaging and
hippocampal volumetry in the lateralization of mesial temporal lobe epilepsy,?
Neurology, vol. 51, pp. 66?71, July 1998.
[13] D. L. Arnold et al., ?Proton magnetic resonance spectroscopic imaging for
metabolic characterization of demyelinating plaques,? Annals of Neurology, vol. 31,
no. 3, pp. 235?241, 1992.
[14] D. J. Meyerho?, S. MacKay, L. Bachman, N. P. W. P. Dillon, M. W. Weiner,
and G. Fein, ?Reduced brain n-acetylaspartate suggests neuronal loss in cogni-
tively impaired human immunodeficiency virus- seropositive individuals: in vivo
1H magnetic resonance spectroscopic imaging,? Neurology, vol. 43, pp. 509?515,
1993.
[15] P. B. Barker, R. R. Lee, and J. C. McArthur, ?AIDS dementia complex: evaluation
with proton MR spectroscopic imaging,? Radiology, vol. 195, pp. 58?64, 1995.
[16] D. J. Meyerho?, S. MacKay, J. M. Constant, D. Norman, C. V. Dyke, G. F.
C, and M. W. Weiner, ?Axonal injury and membrane alterations in Alzheimer?s
disease suggested by in vivo proton magnetic resonance spectroscopic imaging,?
Ann. Neurol, vol. 36, pp. 40?47, 1994.
[17] G. Tedeschi, A. Bertolino, N. Lundbom, S. Bonavita, N. J. Patronas, J. H. Duyn,
L. V. Metman, T. N. Chase, and G. D. Chiro, ?Cortical and subcortical chem-
ical pathology in Alzheimer?s disease as assessed by multislice proton magnetic
resonance spectroscopic imaging,? Neurology, vol. 47, pp. 696?704, 1994.
[18] B. L. Miller, R. Moats, T. Shonk, T. Ernst, S. Woolley, and B. D. Ross,
?Alzheimer disease: depiction of increased cerebral myo-inositol with proton MR
spectroscopy,? Radiology, vol. 187, pp. 433?439, May 1993.
[19] W. J. Chu, H. P. Hetherington, R. I. Kuzniecky, T. Simor, G. F. Mason, and
G. A. Elgavish, ?Lateralization of human temporal lobe epilepsy by 31P NMR
spectroscopic imaging at 4.1 T,? Neurology, vol. 51, pp. 472?479, August 1998.
[20] F. Cendes, Z. Caramanos, F. Andermann, F. Dubeau, and D. L. Arnold, ?Proton
magnetic resonance spectroscopic imaging and magnetic resonance imaging vol-
umetry in the lateralization of temporal lobe epilepsy: a series of 100 patients,?
Ann. Neurol, vol. 42, pp. 737?746, 1997.
112
[21] J. W. Hugg, K. D. Laxer, G. B. Matson, A. A. Maudsley, C. A. Husted, and
M. W. Weiner, ?Lateralization of human focal epilepsy by 31P magnetic resonance
spectroscopic imaging,? Neurology, vol. 42, pp. 2011?2018, 1992.
[22] O. A. Petro?, J. W. Prichard, T. Ogino, and R. G. Shulman, ?Proton magnetic res-
onance spectroscopic studies of axonal carbohydrate metabolism in rabbit brain,?
Neurology, vol. 38, pp. 1569?1574, 1988.
[23] P. A. Bottomley, C. J. Hardy, and P. B. Roemer, ?Phosphate metabolite imaging
and concentration measurements in human heart by nuclear magnetic resonance,?
Magn. Reson. Med., vol. 14, pp. 425?434, 1990.
[24] H. P. Hetherington, D. J. Luney, J. T. Vaughan, J. W. Pan, S. L. Ponder,
O. Tschendel, D. B. Twieg, and G. M. Pohost, ?31P spectroscopic imaging of
the human heart at 4.1 T,? Magn. Reson. Med., vol. 33, pp. 427?431, 1995.
[25] A. J. Farrall, R. T. Thompson, G. Wisenberg, C. M. Campbell, and D. J. Drost,
?Myocardial infarction in acanine modelmonitored bytwo-dimensional 31P chem-
ical shift spectroscopic imaging,? Magn. Reson. Med., vol. 38, pp. 577?584, 1997.
[26] J. S. Taylor, D. B. Vigneron, J. Murphy-Bosch, S. J. Nelson, H. B. Kessler, l. Coia,
W. Curran, and T. R. Brown, ?Free magnesium levels in normal brain and brain
tumors: 31P chemical-shift imaging measurements at 1.5 T,? Proc. Natl. Acad.
Sci. USA, vol. 88, pp. 6810?6814, 1991.
[27] A. A. Maudsley, S. K. Hilal, W. H. Perman, and H. E. Simon, ?Spatially resolved
high-resolution spectroscopy by four-dimensional NMR,? J. Magn. Reson., vol. 51,
pp. 147?152, 1983.
[28] T.R.Brown, B.M.Kincaid, andK.Ugurbil, ?NMRchemical shiftimaging in three
dimensions,? Proceedings of the National Academy of Sciences, vol. 79, pp. 3523?
3526, 1982.
[29] F. Bloch, ?Nuclear induction,? Phys. Rev., vol. 70, pp. 460?474, 1946.
[30] P. C. Lauterber, ?Image formation by induced local interactions: Examples em-
ploying nuclear magnetic resonance,? Nature, vol. 242, pp. 190?191, 1973.
[31] R. Damadian, ?Tumor detection by nuclear magnetic resonance,? Science, vol. 171,
pp. 1151?3, 1971.
[32] E. Andrew, P. Bottomley, W. Hinshaw, G. Holland, and W. Moore, ?NMR images
by the multiple sensitive point method: application to larger biological systems,?
Physics in Medicine and Biology, vol. 22, pp. 971?974, Sept. 1977.
113
[33]W.S.Hinshaw,E.R.Andrew,P.A.Bottomley,G.N.Holland,andW.S.Moore,
?Display of cross sectional anatomy by nuclear magnetic resonance imaging,?
British journal of Radiology, vol. 51, pp. 273?80, Apr. 1978.
[34] W. Moore, G. Holland, and L. Kreel CT, vol. 4, p. 1, 1980.
[35] G. Holland, R. Hawkes, and W. Moore, ?Nuclear magnetic resonance NMR to-
mography of the brain: coronal and sagittal sections,? J. Comput. Assist. Tomogr,
vol. 4, pp. 429?33, 1980.
[36] D. D. Stark and W. G. Bradley, ?Magnetic resonance imaging,? C. V. Mosby
Company, St. Louis, 1988.
[37] B. R. Friedman, J. P. Jones, G. C. M. noz, A. P. Salmon, and C. R. B. Merritt,
?Principles of MRI,? McGraw-Hill Information Services Company, 1989.
[38] E. M. Haacke, R. W. Brown, M. R. Thompson, and R. Venkatesan, ?Physical
principles and sequence design,? John Wiley and Sons, Inc. , New York, 1989.
[39] D. Shaw, ?Fourier transform NMR spectroscopy.,? New Yor: Elsevier Scientific,
1971.
[40] C. P. Slichter, ?Principles of magnetic resonance,? Berlin: Springer-verlag,vol.1,
1978.
[41] A. Abraham, ?The principles of nuclear magnetism,? Oxford: Oxford University
Press, 1961.
[42] E. L. Hahn, ?Spin echoes,? Phys. Rev., vol. 80, pp. 580?594, 1950.
[43] H. Y. Carr and E. M. Purcell, ?E?ects of di?usion on free precession in nuclear
magnetic resonance experiments,? Phys. Rev., vol. 94, pp. 630?638, 1954.
[44] S. Meiboom and D. Gill, ?Modified spin-echo method for measuring nuclear relax-
ation times,? Rev. Sci. Instrum., vol. 29, pp. 668?691, 1958.
[45] D. E. Woessner, ?E?ects of di?usion in nuclear magnetic resonance spin-echo ex-
periments,? J. Chem. Phys., vol. 34, pp. 2057?2061, 1961.
[46] Z. H. Cho, H. S. Kim, H. B. Song, and J. Cumming, ?Fourier transform nuclear
magnetic resonance tomographic imaging,? IEEE Proceedings, no. 70, pp. 1152?
1173, 1982.
[47] W. S. Hinshaw and A. Lent., ?An introduction to NMR imaging: from the Bloch
equation to the imaging equation,? IEEE proceedings, vol. 71, pp. 338?350, 1983.
[48] A. Kumar, D. Welti, and R. Ernst., ?NMR Fourier zeugmatography,? J. Magn.
Reson, vol. 18, p. 69, 1975.
114
[49] P. R. Luyten, A. J. H. Marien, B. Sijtsma, and J. A. D. Hollander., ?Solvent-
suppressed spatially resolved spectroscopy. an approach to high-resolution NMR
on a whole body system,? J. Magn. Reson., vol. 67, p. 148, 1986.
[50] J. Frahm, H. Bruhn, M. L. Gyngell, K. D. Merboldt, W. H?anicke, and R. Sauter.,
?Localized proton NMR spectroscopy in di?erent regions of the human brain in
vivo. Relaxation times and concentrations of cerebral metabolites,? Magn. Reson.
Med., vol. 11, pp. 47?67, July 1989.
[51] W. P. Aue, S. Mueller, T. A. Cross, and J. Seelig., ?Volume-selective excitation.
A novel approach to topical NMR,? J. Magn. Reson., vol. 56, no. 2, pp. 350?354,
1984.
[52] S. Mueller, W. P. Aue, and J. Seelig., ?NMR imaging and volume-selective spec-
troscopy with a single surface coil,? J. Magn. Reson., vol. 63, no. 3, pp. 530?543,
1985.
[53] P. R. Luyten and J. A. D. Hollander. Proc. SMRM, vol. 2, p. 1021, 1985.
[54] P.J.Hore., ?SolventsuppressioninFouriertransformnuclearmagneticresonance,?
J. Magn. Reson, vol. 55, no. 2, pp. 283?300, 1983.
[55] J. Frahm, K. D. Merboldt, and W. H?anicke., ?Localized proton spectroscopy using
stimulated echoes,? J. Magn. Reson., vol. 72, no. 3, pp. 502?507, 1987.
[56] C. Y. Rim, J. B. Ra, and Z. H. Cho., ?Radial scanning technique for volume
selective 31P spectroscopy,? Magn. Reson. Med., vol. 24, pp. 100?8, Mar. 1992.
[57] K. D. Merboldt, D. Chien, W. H?anicke, M. L. Gyngell, and H. Bruhn., ?Localized
31p NMR spectroscopy of the adult human brain in vivo using stimulated-echo
(STEAM) sequences,? J. Magn. Reson., vol. 89, no. 2, pp. 343?361, 1990.
[58] P. Styles, C. A. Scott, and G. K. Radda, ?A method for localizing high-resolution
NMR spectra from human subjects,? Magn. Reson. Med., vol. 2, p. 402, 1985.
[59] M. Garwood, T. Schliech, G. B. Matson, and G. Acosta, ?Spatial localization of
tissue metabolites by phosphorus-31 NMR rotating frame zeugmatography,? J.
Magn. Reson, vol. 60, no. 2, pp. 268?279, 1984.
[60] P. A. Bottomley, T. H. Foster, and R. D. Darrow., ?Depth-resolved surface-coil
spectroscopy (DRESS) for in vivo 1H, 31P, and 13C NMR,? J. Magn. Reson.,
vol. 59, no. 2, pp. 338?342, 1984.
[61] P. R. Luyten and J. A. D. Hollander, ?1H MR spatially resolved spectroscopy of
human tissues in situ,? Magnetic. Resonance Imaging, vol. 4, no. 3, pp. 237?239,
1986.
115
[62] P. R. Luyten, C. M. Anderson, and J. A. D. Hollander, ?1H NMR relaxation mea-
surements of human tissues in situ by spatially resolved spectroscopy,? Magnetic
Resonance in Medicine, vol. 4, pp. 431?40, May 1987.
[63] J. Frahm, H. Bruhn, M. L. Gyngell, K. D. Merboldt, W. H?anicke, and R. Sauter,
?Localized high-resolution proton NMR spectroscopy using stimulated echoes: ini-
tial applications to human brain in vivo,? Magn. Reson. Med., vol. 9(1), pp. 79?93,
1989.
[64] S. Lee and Z.Cho, ?Localized volume selection technique using an additional radial
gradient coil,? Magnetic Resonance in Medicine, vol. 12, pp. 56?63, October 1989.
[65] C. J. Hardy and H. E. Cline, ?Spatial localization in two dimensions using NMR
designer pulses,? J. Magn. Reson, vol. 82, p. 647, 1989.
[66] P. A. Bottomley and C. J. Hardy, ?Two dimensional spatially selective spin inver-
sion and spin-echo refocusing with a single nuclear magnetic resonance pulse,? J.
Appl. Phys., vol. 62, pp. 4284?90, 1987.
[67] P. A. Bottomley and C. J. Hardy, ?Progress in e?cient three-dimensional spatially
localized in vivo 31p NMR spectroscopy using multidimensional spatially selective
(?) pulses,? J. Magn. Reson, vol. 74, no. 3, pp. 550?556, 1987.
[68] C. J. Hardy, P. A. Bottomley, and P. B. Roemer, ?O?-axis spatial localization with
frequency modulated nuclear magnetic resonance rotating (?) pulses,? J. Appl.
Phys., vol. 63, p. 4741, 1988.
[69] C. J. Hardy, P. A. Bottomley, , M. O?Donnel, and P. B. Roemer, ?Optimization of
two-dimensional spatially selective NMR pulses by stimulated annealing by stim-
ulated,? J. Magn. Reson, vol. 77, no. 2, pp. 233?250, 1988.
[70] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, ?Optimization by stimulated an-
nealing,? Science, vol. 220, pp. 671?680, 1983.
[71] N. Metropolis, A. W. Rosenbluth, A. H. Teller, and E. Teller, ?Equation of state
calculations by fast computing machines,? J. Chem. Phys., vol. 21, p. 1087, 1953.
[72] J. Pauly, D. Nishimura, and A. Macovski, ?A k-space analysis of small-tip-angle
excitation,? J. Magn. Reson., vol. 81, pp. 43?56, 1989.
[73] J. Pauly, D. Nishimura, and A. Macovski Proc. SMRM, vol. 7, p. 654, 1988.
[74] J. B. Ra, C. Y. Rim, and Z. H. Cho, ?Application of single-shot spiral scanning
for volume localization,? Magn. Reson. Med., vol. 17, pp. 423?433, Feb. 1991.
[75] C. B. Ahn, J. H. Kim, and Z. H. Cho, ?High speed spiral-scan echo planar NMR
imaging,? IEEE Tran. Med. Imag., vol. 5, no. 1, pp. 2?7, 1986.
116
[76] C. Y. Rim, J. B. Ra, and Z. H. Cho Proc. SMRM, vol. 8, p. 30, 1989.
[77] C. J. Hardy and H. E. Cline Proc. SMRM, vol. 8, p. 26, 1989.
[78] D. L. Arnold, P. M. Matthews, G. S. Francis, J. O?Connor, and J. P. Antel, ?Pro-
ton magnetic resonance spectroscopic imaging for metabolic characterization of
demyelinating plaques,? Ann. Neurol, vol. 31, pp. 235?241, 1992.
[79] X. Hu, D. N. Levin, P. C. Lauterbur, and T. A. Spraggins, ?SLIM: Spectral local-
ization by imaging,? Magnetic Resonance in Medicine, vol. 8, pp. 314?322, 1988.
[80] Z. Liang and P. C. Lauterbur, ?A generalized series approach to MR spectroscopic
imaging,? IEEE Transactions on Medical Imaging, vol. 10, pp.132?137, June1991.
[81] E. M. Stokely and D. B. Twieg, ?Functional image reconstruction enhancement
(FIRE)forMRspectroscopicandnuclear medicineimages,? in SPIE Medical Imag-
ing VI, 1992.
[82] P. Mansfield, ?Multi-planar image formation using NMR spin echoes,? J. Phys. C:
Solid State Phys., vol. 10, pp. 155?l58, 1977.
[83] S. Ljunggren, ?A simple graphical representation of Fourier-based imaging meth-
ods,? J. Magn. Reson, vol. 54, no. 2, pp. 338?343, 1983.
[84] J. H. Duyn and C. T. Moonen, ?Fast proton spectroscopic imaging of human brain
using multiple spin-echoes,? Magn Reson Med, vol. 30, pp. 409?422, Oct. 1993.
[85] J. H. Duyn, J. A. Frank, and C. T. Moonen, ?Incorporation of lactate measure-
ment in multi-spin-echo proton spectroscopic imaging,? Magn Reson Med, vol. 33,
pp. 101?107, Jan. 1995.
[86] E. Adalsteinsson, P. Irarrazabal, D. M. Spielman, and A. Macovski, ?Three-
dimensional spectroscopic imaging with time-varying gradients,? Magn. Reson.
Med, vol. 33, pp. 461?466, 1995.
[87] D. G. Norris and W. Dreher, ?Fast proton spectroscopic imaging using the sliced
k-space method,? Magn. Reson. Med. 30, 641-645, vol. 30, pp. 641?645, 1993.
[88] S. Posse, G. Tedeschi, R. Risinger, and D. L. R. Ogg, ?High speed 1H spectroscopic
imaging in human brain by echo planar spatial-spectral encoding,? Magn. Reson.
Med, vol. 33, pp. 34?40, 1995.
[89] J. C. Ehrhardt, ?MR data acquisition and reconstruction using e?cient sampling
schemes,? IEEE Transactions on Medical Imaging, vol. 9, pp. 305?309, September
1990.
117
[90] R. M. Mersereau, ?The processing of hexagonally sampled two-dimensional sig-
nals,? Proceeding of IEEE, vol. 67, pp. 930?949, 1979.
[91] D. P. Petersen and D. Middleton, ?Sampling and reconstruction of wave-number
limited functions in n-dimensional euclidean spaces,? Information and Control,
vol. 5, pp. 279?323, 1962.
[92] M. von Kienlin, T. Ceckler, R. Mejia, and R. S. Balaban, ?Numerical optimization
of the SLIM experiment for localized spectroscopy,? in Proceedings of the Society
of Magnetic Resonance, (New York, NY), 1990.
[93] Y. Caoand D. N.Levin, ?Usingan image database to constrain theacquisition and
reconstruction of MR images of the human head,? IEEE Transactions on Medical
Imaging, vol. 14, pp. 350?361, June 1995.
[94] S. K. Plevritis and A. Macovski, ?Alternative k-space sampling distributions for
MR spectroscopic imaging,? in Proceedings of the 1994 IEEE International Con-
ference on Image Processing, vol. III, pp. 11?14.
[95] R. Venkataramani and Y. Bresler, ?Further results on spectrum blind sampling of
2D signals,? in Proceedings of the 1998 IEEE International Conference on Image
Processing, vol. II, pp. 752?756, 1998.
[96] S.K.Nagle andD.N.Levin, ?AnewclassofsamplingtheoremsforFourierimaging
of multiple regions,? in Proceedings of the 1998 IEEE International Conference on
Image Processing, vol. II, 1998.
[97] S. K. Nagle and D. N. Levin, ?Multiple region MRI,? Magnetic Resonance in
Medicine, vol. 41, pp. 774?786, 1999.
[98] S. J. Reeves and L. P. Heck, ?Selection of observations in signal reconstruction,?
IEEE Transactions on Signal Processing, vol. 43, pp. 788?791, March 1995.
[99] S. J. Reeves, ?Selection of k-space samples in localized spectroscopy of arbitrary
volumes of interest,? Journal of Magnetic Resonance Imaging, vol. 5, pp. 245?247,
March/April 1995.
[100] S. J. Reeves and Z. Zhao, ?Sequential algorithms for observation selection,? IEEE
Transactions on Signal Processing, vol. 47, pp. 123?132, January 1999.
[101] K. F. Cheung, Advanced Topics in Shannon Samping and Interpolation Theory,
ch. A Multidimensional Extension of Papoulis? Generalized Sampling Expansion
with the Application in Minimum Density Sampling, pp. 85?119. New York:
Springer-Verlag, 1993.
118
[102] Y. Gao and S. J. Reeves, ?An e?cient computation for sequential forward obser-
vation selection in image reconstruction,? in Proceedings of the 1998 IEEE Inter-
national Conference on Image Processing, vol. 3, pp. 380?384, 1998.
[103] Y. Gao and S. J. Reeves, ?Optimal k-space sampling in MRSI for images with a
limited region of support,? submitted to IEEE Transactions on Medical Imaging.
[104] Y. Gao and S. J. Reeves, ?E?cient forward selection of k-space samples for images
with a limited region of support,? in BMES-EMBS Conference, 1999.
[105] Y. Gao and S. J. Reeves, ?E?cient backward selection of k-space samples in MRI
on a hexagonal grid,? Circuits, Systems, and Signal Processing, to appear Aug
2000.
[106] Y. Gao and S. J. Reeves, ?Optimal dithering of focal plane arrays in passive
millimeter-wave imaging,? Optical Engineering, vol. 40, pp. 2179?2187, October
2001.
[107] Y. Gao and S. J. Reeves, ?E?cient algorithm for optimizing dithering pattern in
passive millimeter-wave imaging,? the SPIE?s 14th Annual International Sympo-
sium on Aerospace/Defense, Apr. 2000.
[108] Y. Gao and S. J. Reeves, ?Optimal sampling in array-based image formation,?
accepted by 2000 IEEE International Conference on Image Processing.
[109] S. J. Reeves and L. P. Heck, ?Selection of observations in signal reconstruction,?
in Proceedings of the 1993 IEEE International Conference on Acoustics, Speech,
and Signal Processing, vol. III, pp. 444?447.
[110] G. H. Golub and C. Van Loan, Matrix Computations. Baltimore, MD: Johns
Hopkins University Press, 2nd ed., 1989.
[111] Standard Mathematical Tables. CRC Press, Inc., 30th ed., 1995.
[112] R.L.Lagendijk, R. M.Mersereau, andJ.Biemond, ?On increasingthe convergence
rate of regularized iterative image restoration algorithms,? in Proceedings of the
1987 IEEE International Conference on Acoustics, Speech, and Signal Processing,
pp. 1183?1186.
[113] S. J. Reeves and Z. Zhao, ?New results on observation selection in signal recon-
struction,? in Proceedings of the 1996 IEEE International Conference on Acoustics,
Speech, and Signal Processing, vol. III, pp. 1676?1679, 1996.
[114] W. H. Press et al., Numerical Recipes in C. Cambridge, 1988.
[115] D. E. Dudgeon and R. M. Mersereau, Multidimensional Digital Signal Processing.
New Jersey: Prentice-Hall, 1984.
119
[116] J. C. Ehrhardt, ?MR data acquisition and reconstruction using e?cient sampling
schemes,? IEEE Transactions on Medical Imaging, vol. 9, pp. 305?309, September
1990.
[117] A. Kaufmann and A. Henry-Labord`ere, Integer and Mixed Programming: Theory
and Applications. New York: Academic, 1977.
[118] A. Guessoum, Fast Algorithms for the Multimensional Discrete Fourier Transform.
PhD thesis, Georgia Institute of Technology, Mar. 1984.
[119] B. L. Evans, T. R. Gardos, and J. H. McClellan, ?Imposing Structure on Smith-
Form Decompositions of Rational Resampling Matrices,? IEEE Transactions on
Signal Processing, vol. 42, pp. 970?973, April 1994.
[120] A. K. Jain, Fundamentals of Digital Image Processing. Englewood Cli?s, NJ:
Prentice-Hall, 1989.
[121] A.V.OppenheimandR.W.Schafer, Discrete-Time Signal Processing. NewJersey:
Prentice-Hall, 1989.