E cient Acquisition and Reconstruction for Magnetic Resonance Spectroscopic Imaging by Wenting Deng A dissertation submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama May 5th 2013 Keywords: Magnetic Resonance Spectroscopic Imaging, Sequential Data Selection, Nonlinear Reconstruction Copyright 2013 by Wenting Deng Approved by Stanley Reeves, Chair, Professor of Electrical and Computer Engineering Thomas Denney, Professor of Electrical and Computer Engineering Bogdan Wilamowski, Professor of Electrical and Computer Engineering Abstract Magnetic resonance spectroscopic imaging (MRSI), combining both magnetic resonance spectroscopy (MRS) and magnetic resonance imaging (MRI) techniques, has proven to be a powerful approach to reveal information about metabolite distributions and discriminate multiple resonant frequencies in the spectrum. Due to its nondestructive nature and sensi- tivity to the molecular environment of individual atoms [1], MRSI is widely applied in the clinical community. The typical applications include mapping abnormal tissues in the brain [2, 3, 4], in the prostate [3, 5], in the breast [6, 7], as well as pathologic analysis after resection operations [8]. All of these clinical studies require satisfactory resolution in both spatial and spectral dimensions, which in conventional MRSI demands a great deal of acquisition time. Unfortunately, patient discomfort, motion artifacts and cost will signi cantly increase when the acquisition time lengthens. To overcome these problems, we propose to implement an imaging protocol that only acquires an optimal subset of data to accelerate the collection process without sacri cing spatial and spectral accuracy in reconstruction. In most MRSI applications, the spectral domain has great sparsity, which raises the pos- sibility of reconstructing spectral information with limited time series data [9]. Therefore, an e cient sequential backward selection (SBS) [10] technique is proposed to select a limited set of but the most informative echo-time values, which are then applied to echo-planar imaging (EPI) acquisition [11]. By exploiting multi-echo EPI, multiple k-space frames can be ac- quired within one excitation to further reduce the acquisition time. To achieve this purpose, we modify the selection method to a more e cient approach. Instead of selecting echo-time value one by one, the modi ed algorithm selects multiple echo-time values simultaneously, which would then be used in one excitation acquisition. ii For the EPI technique, a k-space frame cannot be collected instantly. In other words, every k-space sample will have a di erent time delay even in the same k-space frame. Con- sequently, selecting one echo-time value for a whole k-space frame might not be accurate enough. We then extend the data selection method to both k-space and time domains. In addition, an advanced EPI strategy is introduced. On the other hand, if SBS algorithm is the only restriction for k-t data selection, the acquisition e ciency might be reduced, which leads to longer observation procedure. Considering this, sequential k-t selection with constraint will be studied for a better selection result in a more e cient way. Due to the selection method and the EPI acquisition technique, the collected data are normally nonuniform and time varying. Therefore, fast Fourier Transforms (FFTs) are not capable reconstructing the spatial and spectral information directly. Besides, an FFT cannot separate spatial information of di erent metabolite resonances. On the other side, conventional optimization methods [12], such as conjugate gradients (CG), require very high computational e ort and large memory storage to nd the matching parameters of the images. A fast reconstruction method combining polynomial approximation with FFTs is investigated, which can greatly accelerate the reconstruction process without sacri cing reconstruction quality. During the reconstruction procedure, the nite data set raises a practical problem: frequency local minima, which means convergence to the global minimum is not guaranteed. In order to overcome this problem, we study the origin of the local minima, and propose an easily implemented method: varying estimated decays during the optimization. However, this method has its own limitations, especially in reconstruction e ciency and accuracy. Therefore, another advanced technique, applying weighted scalars to the cost function model, is presented. The results show that the second method can e ciently attenuate frequency local minima e ect, while avoiding the reconstruction speed and accuracy problems. iii Acknowledgments It has been a great pleasure to be a Ph.D student working closely with Dr. Stanley J. Reeves. I would like to give my heartfelt thanks for his guidance and encouragement over the years. He led me into the medical imaging world and supported me throughout my research process with his patience and knowledge. This work would never have existed without him. I would like to thank my advisory committee members, Dr. Thomas Denney and Dr. Bogdan Wilamowski for their generous comments and suggestions on this dissertation. I also thank Dr. Tin-Yau Tam for serving as the outsider reader. I would like to thank Dr. Donald B. Twieg from the University of Alabama at Birm- ingham, Dr. Nouha Salibi from Seimens for their kindness help of subjects preparation and data acquisition. In addition, my appreciation goes to all the previous and present group members and collaborators, especially Weidong Tang, Tao Ma, Kevin Perry, Chun G. Shiros, Wei Zha, Bharath Ambale Venkatesh and Ningzhi Li, for their help and friendship throughout the course of this research. Finally, I want to thank my beloved family for their unconditional support to me in every possible way. Their trust and con dence in me is my forward power. My wonderful husband, Dayu Yang, makes me become a better person. Without his endless love, care and patience, I would not be who I am now. My beloved parents, Yucheng Deng and Lizhong Yang, give their hearts to me. They do everything good for me and everything they can to support me in every step in my life. I would like to dedicate this dissertation to my dearest family. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Organization of the Proposal . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Physical Bases of Nuclear Magnetic Resonance . . . . . . . . . . . . . . . . 2 1.2.1 Chemical Shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 NMR Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 MRSI Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Imaging Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.2 Echo-Planar Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2 E cient Data Acquisition with Echo-Time Selection . . . . . . . . . . . . . . . 13 2.1 Echo-Time Selection Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.1 Observation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.2 Selection Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Phantom Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.1 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3 Alternative Optimal Data Acquisitions . . . . . . . . . . . . . . . . . . . . . . . 36 3.1 Optimal k-t Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.1.1 k-t Selection Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 v 3.1.2 k-t Selection with Constraint . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Echo-Time Selection with Overlapped EPI . . . . . . . . . . . . . . . . . . . 43 3.3 Experiments and Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3.1 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4 Fast Reconstruction for Nonlinear Model . . . . . . . . . . . . . . . . . . . . . . 56 4.1 Introduction to Iterative Reconstruction Method . . . . . . . . . . . . . . . . 56 4.1.1 Parameter Gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1.2 Line Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.2 Fast Reconstruction Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2.1 Reconstruction Model Directly Utilizing FFTs . . . . . . . . . . . . . 65 4.2.2 Reconstruction model based on polynomial approximations . . . . . . 67 4.3 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.1 Simulation experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.2 Phantom Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.3.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5 Optimization Techniques to Attenuate Frequency Local Minima E ect . . . . . 80 5.1 Origins of Local Minima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.2 Optimization Protocol with Varying Decay . . . . . . . . . . . . . . . . . . . 85 5.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2.2 Experiments and Results . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.3 Optimization Protocol with Weighted Scalars . . . . . . . . . . . . . . . . . 92 5.3.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.2 Experiments and Results . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.4 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 vi 6.1 Summary of the Contributions of This Dissertation . . . . . . . . . . . . . . 100 6.2 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 vii List of Figures 1.1 Stationary states of nuclear spins in a static magnetic eld . . . . . . . . . . . . 3 1.2 Molecular structure of (a) ethanol (b) acetone . . . . . . . . . . . . . . . . . . . 6 1.3 Molecular structure of butanone . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 (a) EPI pulse sequence (b) EPI trajectory . . . . . . . . . . . . . . . . . . . . . 10 1.5 EPI images with di erent echo-time values . . . . . . . . . . . . . . . . . . . . . 11 2.1 (a) k-space sample selection, (b) echo-time sample selection . . . . . . . . . . . 14 2.2 (a) Global frequency shift (b) local variations . . . . . . . . . . . . . . . . . . . 19 2.3 Basic multi-echo EPI sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 (a) Elimination pattern (b) periodic nonuniform elimination . . . . . . . . . . . 22 2.5 Phantom structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.6 1H spectrum of methanol (upper plot) 1H spectrum of ethanol (lower plot) . . . 24 2.7 Number of remaining echo-time values vs. tr (AHA) 1 (which is proportional to MSE) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.8 Spatial images reconstructed from 500 echo-time values ( rst row), 48 optimized echo-time values (second row), 48 randomly selected echo-time values (third row), 48 equally spaced echo-time values (fourth row). Each column from left to right represents spatial distributions of: hydroxyl, ethyl, methyl of methanol, methyl of ethanol. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 viii 2.9 Spectra in the tube of water and ethanol mixture (upper plot), in the tube of water and methanol mixture (middle plot), in the big glass cylinder of pure water (lower plot). Each plot from lower left to upper right from: 500 echo-time values, 48 optimized echo-time values, 48 randomly selected echo-time values, and 48 equally spaced echo-time values. Plots are o set for easier viewing. . . . . . . . 29 2.10 Experiment ow: Parameter o sets a ect the selection performance . . . . . . . 31 2.11 MSE increase as a function of (a) magnitude o set (b) decay o set (c) frequency o set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.12 Spectrum of equally-spaced selection (a) from a small echo-time interval (b) from a large echo-time interval. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.13 MSE as a function of the number of remaining echo-time values (a) with equally spaced selection (b) with random selection . . . . . . . . . . . . . . . . . . . . . 35 3.1 (a) Standard EPI with uniform k-space acquisition (b) Modi ed EPI with nonuni- form k-space acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 ky-t data selection strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Trajectory of (a) multi-echo EPI selection (b) k-t selection without constraint . 42 3.4 Flow of the updated method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.5 (a) Standard EPI with uniform k-space acquisition (b,c) Modi ed EPI with coarse ky resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.6 (a) Traditional multi-echo EPI with uniform k-space acquisition for two excitation (b) Modi ed multi-echo EPI with coarse ky resolutions for two excitation . . . 48 3.7 Data distribution of conventional phase-encoding imaging . . . . . . . . . . . . 48 ix 3.8 Phantom structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.9 Spatial images on the left picture are reconstructed from full k-space data ( rst row), 3072 optimized (ky;t) samples out of 65,536 (second row), 48 optimized echo-time values using regular echo-time selection (third row), 48 optimized echo- time values using overlapped echo-time selection (fourth row). The rst column represents spatial distributions of hydroxyl; the second column represents spatial distributions of methyl. Spatial images on the right picture are reconstructed from 48 equally spaced echo-time values. The upper image represents spatial distribution of hydroxyl and the lower image represents spatial distribution of methyl. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.10 Spectra in the big glass cylinder and lower left tube full of pure water (left plot), in the tubes of water and methanol mixture (right plot). Each plot from lower left to upper right from: full data set, 3072 optimized (ky;t) samples out of 65,536, 48 optimized echo-time values using regular echo-time selection, 48 optimized echo-time values using overlapped echo-time selection, and 48 equally spaced echo-time values. Plots are o set for easier viewing. . . . . . . . . . . . . . . . . 54 4.1 Original vs. quadratic approximation of line search . . . . . . . . . . . . . . . . 62 4.2 Quadratic approximation of line search based on di erent intersections . . . . . 63 4.3 Original exponential function vs. polynomial approximation . . . . . . . . . . . 69 4.4 Absolute error varies with polynomial order . . . . . . . . . . . . . . . . . . . . 70 x 4.5 Original simulated spatial images ( rst row), spatial images reconstructed from 48 optimized echo-time values with regular conjugate gradient method (second row), spatial images reconstructed from 48 optimized echo-time values with polynomial approximation method (third row), spatial images reconstructed from 48 opti- mized echo-time values with modi ed conjugate gradient method (fourth row). Each column from left to right represents spatial distributions of: hydroxyl, ethyl, methyl of methanol, methyl of ethanol. . . . . . . . . . . . . . . . . . . . . . . . 76 4.6 Spatial images reconstructed from full k-space data with regular conjugate gra- dient method ( rst row), spatial images reconstructed from 48 optimized echo- time values with regular conjugate gradient method (second row), spatial images reconstructed from 48 optimized echo-time values with polynomial approxima- tion method (third row), spatial images reconstructed from 48 optimized echo- time values with modi ed conjugate gradient method (fourth row). Each column from left to right represents spatial distributions of: hydroxyl, ethyl, methyl of methanol, methyl of ethanol. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.1 Cost function varying with frequency and decay o sets (upper plot) frequency and magnitude o sets (middle plot) magnitude and decay o sets (lower plot) . . 83 5.2 Cost function varying with frequency at (a) di erent magnitude o sets (b) dif- ferent decay o sets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.3 Overview of the proposed method . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.4 Cost function with (a) decay=-2 (b) decay=-20 (c) decay=-50 . . . . . . . . . . 88 5.5 Cost function with (a) decay=-2 (b) decay=-20 (c) decay=-50 . . . . . . . . . . 90 5.6 Cost function with estimated (a) decay=-2 (b) decay=-20 (c)decay=-50 . . . . . 92 xi 5.7 Typical cost function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.8 Cost function using estimated decay =-2 (a) without weighted scalars (b) with weighted scalars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.9 Cost function using estimated decay =-2 (a) without weighted function (b) with weighted function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.10 Cost function using estimated decay =-2 (a) without weighted function (b) with weighted function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 xii List of Tables 2.1 Selection Time Comparison (s) . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1 Selection Time(s) and MSE Comparison . . . . . . . . . . . . . . . . . . . . . . 49 4.1 MSEPCT of reconstructed magnitude (%) . . . . . . . . . . . . . . . . . . . . . . 75 4.2 MSEPCT of reconstructed decay (%) . . . . . . . . . . . . . . . . . . . . . . . . 75 4.3 MSEPCT of reconstructed frequnecy (%) . . . . . . . . . . . . . . . . . . . . . . 75 4.4 Comparison of reconstruction time for 100 iterations (s) . . . . . . . . . . . . . 77 5.1 Comparison of o set range (Hz)and number of iterations . . . . . . . . . . . . . 89 5.2 Comparison of o set range (Hz) and number of iterations . . . . . . . . . . . . 91 5.3 Comparison of o set range (Hz) and number of iterations . . . . . . . . . . . . 91 5.4 Comparison of o set range (Hz) and number of iterations . . . . . . . . . . . . 95 5.5 Comparison of o set range (Hz) and number of iterations . . . . . . . . . . . . 98 5.6 Comparison of o set range (Hz) and number of iterations . . . . . . . . . . . . 98 xiii Chapter 1 Introduction Magnetic resonance spectroscopic imaging, which is a combination of magnetic reso- nance spectroscopy (MRS) and magnetic resonance imaging (MRI), is an important appli- cation of the phenomenon of nuclear magnetic resonance (NMR). It can produce metabolite distributions and spatially localized spectra inside subjects. Because of its relative safety, MRSI is widely used to diagnose health conditions of human and animals. 1H, 13C, 31P MRSI studies have many applications in clinical practice. They can indicate information about cellular activities in a variety of diseases including brain tumors [2], breast tumors [6], prostate tumors [5], epilepsy [13], as well as abnormalities in various pathologies [14, 15]. Despite its attractive potential, MRSI has remained largely in the realm of health care research because long acquisition time is required in a standard acquisition scheme. Factors like patient comfort, motion artifacts and expense limit the time available to obtain good image resolutions in MRSI. To accelerate the data acquisition and keep the quality of the reconstructed image, we are pursuing strategies where the imaging procedure is optimized so that only an optimal subset of data is collected and used for reconstruction. On the other hand, these acquired data are normally laid out nonuniformly. Fast reconstruction methods, like the Fourier transform, are not applicable anymore. Therefore, we propose to develop a new optimization method to nd the best matching parameters of the images. 1.1 Organization of the Proposal In this chapter, we brie y review the physics of nuclear magnetic resonance and basic principles of magnetic resonance spectroscopic imaging. 1 In Chapter 2, a new e cient acquisition technique is presented. We propose a sequential selection method to choose limited but optimal echo-time values which are then used in echo-planar imaging acquisition. 1H phantom experiments demonstrate that our approach achieves similar results to standard MRSI while using much less acquired data. In Chapter 3, we address limitations of echo-time selection in data acquisition. We propose two alternative data selection strategies: one is k-t selection, which extends the sequential selection technique in both k-space and time domains; another one is overlapped echo-time selection, which makes the acquisition time for each k-space frame shorter and may lead to a fast reconstruction process. In Chapter 4, we summarize the standard reconstruction method used for MRSI. Then a new iterative reconstruction method based on a polynomial model is detailed discussed. This technique decomposes the exponential time function and utilizes fast Fourier transform to speed up the optimization procedure. In Chapter 5, a practical problem of parameter tting in reconstruction is described: local minima. Two methods, varying estimated decays / applying weighted scalars to the cost function model during the optimization, are prsented to solve the local minima problem. Experiments and comparisons indicate their e ciency and simplicity of implementation. Finally, conclusions are given in Chapter 6. The innovative ideas are summarized, and possible future work is discussed. 1.2 Physical Bases of Nuclear Magnetic Resonance Nuclear magnetic resonance (NMR) is a phenomenon when the nuclei [16] of certain atoms are immersed in a magnetic eld, absorbing and re-emitting electromagnetic (EM) energy. This energy owns a speci c resonant frequency that depends on the strength of the external magnetic eld and the chemical environment of the nuclei. Many scienti c techniques exploit NMR phenomena to probe molecule structures as well as detailed internal forms of the observed subjects. 2 Figure 1.1: Stationary states of nuclear spins in a static magnetic eld In the absence of an external magnetic eld, the nuclei have random orientations. The vector sum of these orientations will be zero. When an external magnetic eld is applied, randomly oriented nuclei experience an external force that aligns the nuclei either in a parallel or an anti-parallel direction in reference to the applied magnetic eld. The lower-energy state with the nuclei aligned with the eld is called the -spin state; the higher-energy state with the nuclei aligned against the eld is called the -spin state. This situation is depicted in Figure 1.1. In fact, the energy di erence between two states is proportional to the strength of the external magnetic eld. E = rh2 B0 (1.1) where E is the energy di erence between and states; h is Planck?s constant; B0 is the strength of the external magnetic eld andr is the gyromagnetic ratio (26:753sec 1Gauss 1for a proton). Due to Boltzmann?s law, the lower-energy state always has a larger population of 3 spins than the higher-energy state. However, when a nucleus is irradiated with an extra EM radiation with energy E = rhB0=2 , the nucleus can ip from the lower-energy state to the higher-energy state. This energy is usually supplied by application of an rf pulse signal to the system. After the rf signal, the excited nuclei tend to return to their low-energy state by emitting a well-de ned resonant frequency (the same frequency as the applied rf). This emission of rf signals is then detected by rf coils placed close to the excited object. This is the origin of NMR and the emitted signal is called a free induction decay (FID). Besides (1.1), nuclear energy can also be represented as E = hv (1.2) where h is Planck?s constant and v is the frequency of the EM wave. When (1.1) combines with (1.2) E = rh2 B0 = hv (1.3) The resonant frequency v is proportional to the applied eld B0 v = r2 B0 = 26;753sec 1Gauss 1 2 B0 = 4257:8sec 1Gauss 1 B0 (1.4) For the elds of currently available magnets, the resonant frequencies occur in the radio- frequency range of the spectrum. Because hydrogen is a major component of organic com- pounds and historically NMR was rst used to study protons, we limit our discussion of NMR to proton magnetic resonance. 1.2.1 Chemical Shift In nature, protons do not exist independently. Instead, they are surrounded by electrons. In the presence of an external magnetic eld, the electrons will partially shield the proton 4 from the external eld. Their negative charges generate an opposite magnetic eld to the externally applied eld. As a result, the actual eld strength is weaker than the external eld strength. This is known as chemical shielding. Bactual = Bexternal Bshielding (1.5) vactual = 4257:8sec 1Gauss 1 Bactual (1.6) In di erent molecular environments, protons are shielded by di erent amounts of electrons. Therefore, the resonant frequencies of individual protons vary with the chemical surround- ings. For example, in methanol, the hydroxyl proton is not shielded as much as the methyl proton, so the actual eld strength applied to the hydroxyl proton is stronger. Consequently the hydroxyl proton has a higher resonant frequency than the methyl proton. Generally speaking, chemical shielding is very small compared to the external eld. Therefore, it is not accurate enough to utilize absolute di erences of resonant frequencies to distinguish individual protons. A more reasonable way to express resonant frequency di er- ences is in parts per million (ppm), which de nes frequency di erences between a reference frequency and the observed resonant frequency relative to a reference frequency. = (vsignal vref) 10 6 vref (1.7) The relative frequency di erences expressed in (1.7) are known as chemical shifts. For a given proton, the chemical shift in ppm is the same regardless of the external eld, which makes the chemical shift values easier to interpret. 1.2.2 NMR Signals As we mentioned above, di erent chemical environments lead to di erent chemical shielding of protons. Thus, the number of NMR signals of a molecule corresponds to the number of di erent types of protons present in the molecule. For instance, ethanol, as shown 5 Figure 1.2: Molecular structure of (a) ethanol (b) acetone in Figure 1.2a, has one methyl group (1), one ethyl group (2) and one hydroxyl group (3). The protons from di erent groups have their own chemical shifts. Consequently, three NMR signals can be observed from ethanol. Acetone, as shown in Figure 1.2b, for another instance, has two identical methyl groups (4), where all protons from these two groups have the same chemical shift and are said to be chemically equivalent. In practice, however, there may be fewer signals in the NMR spectrum than the number of types of protons in a molecule, especially a large molecule. The reason is that large molecules may have similar functional groups, which lead to similar resonant frequencies. Association of chemical shifts with di erent types of protons must be done carefully. Unlike the number of resonant peaks, the strength of resonant peaks represented by the area under each peak is only determined by the number of protons contributing to that peak. For example, butanone, as shown in Figure 1.3, has two di erent methyl groups (1)(3). These two groups create two NMR peaks with di erent resonant frequencies; but these two signals have the same area under the peaks. One thing we need to mention here: the number of protons is only proportional to the area under the peak, not the peak height. Therefore, when we reconstruct the magnitudes of resonances, the number of protons cannot be used for accurate proportional estimation of magnitudes, may only be used for the start estimation. 6 Figure 1.3: Molecular structure of butanone 1.3 MRSI Basics MR imaging and MR spectroscopy are two important applications of NMR, which make use of the property of NMR to get spatial and spectral information inside the subjects. The combination of these two techniques is called magnetic resonance spectroscopic imaging (MRSI). For MR spectroscopy, spectral information is only presented at one spatial location, which makes MRS unable to visualize the detailed internal form of subjects. On the other hand, in conventional MR imaging, spatial information is non-time-varying, which makes MRI unable to separate di erent resonant frequencies. Thus, MRI can only map the proton distribution of water. MRSI, instead, gives results that are readily interpretable. At a given resonant frequency, there is a detailed spatial map; while at a given spatial location, a spectral plot is available (examples are shown in experiment section of Chapter 2). MRSI allows collection of spectroscopic data from multiple regions simultaneously. 1.3.1 Imaging Model The basic 2-D form of MRI can be expressed as: s(kx;ky) = Z Z (x;y) exp[ j2 (kxx+kyy)]dxdy (1.8) where 7 kx = 2 Z T 0 Gx( )d ky = 2 Z T 0 Gy( )d Gx( ), Gy( ) are the time-dependent eld gradients along x and y axes, which are applied to the system to distinguish the spatial distribution of spins by Fourier encodings [1]. (x;y) represents the spatial distribution. Notice that the observed s(kx;ky) is in e ect a Fourier transform domain representation of the spatially distributed spin density of (x;y) at the spatial-frequency coordinate (kx;ky). The di erence between MRI and MR spectroscopic imaging is the inclusion of the time domain (or equivalently the spectral domain) to the spatial MR data in MRSI. Therefore, in order to extend the 2-D MRI model to 3-D MRSI model, the original spatial distribution (x;y) needs to include an extra time coordinate becoming a time-varing spatial distribution ~ (x;y;t) (t is later used for spectral reconstruction). In consequence, the observations of MRSI can be expressed as s(kx;ky;t) = Z Z ~ (x;y;t) exp[ j2 (kxx+kyy)]dxdy (1.9) As we mentioned in our previous NMR physical section, after the extra rf signal is turned o , the spin system returns to the low-energy state with the emission of NMR signals. Two relaxation mechanisms are associated with this process: longitudinal(T1) relaxation, which realigns the spins along the original external B0 eld direction; another called transverse or T2 relaxtion, which is a decay process of the FID. In practice, because of magnetic eld in- homogeneity, the observed FID decays even more rapidly than the inherent T2. The e ective relaxation time T 2 is given as 1 T 2 = 1 T2 + B 2 (1.10) 8 where B represents the eld inhomogeneity [17]. Because the readout time is very short in conventional MR imaging, it is not necessary to include these relaxation e ects in imag- ing model. However, MRSI requires relatively long readout time to get the spectroscopic information. In this case, relaxation times can not be ignored. Beacuse the FID signal is the major interest of this research work, in our study, we only consider the e ect of trans- verse relaxation. Now the time-varing spatial distribution ~ (x;y;t) in MRSI model can be represented as ~ (x;y;t) = (x;y)e 1 T 2 (x;y)t (1.11) Normally, in MRSI, we want to observe more than one type of proton. Since the relax- ations and their respective relaxation times are quite sensitive to the chemical enviroments surrounding the nuclei, each type of proton has its own T 2 . In addition to the decay of the FID, each type of proton has its own resonant frequency (in this research work, we treat water proton having 0Hz resonant frequency). Hence, for multiple resonance MRSI, (1.11) is updated as ~ (x;y;t) = X i i(x;y)e 1 T 2i(x;y)te j!i(x;y)t (1.12) where i represents di erent resonant peaks and !i is the frequency o set from water peak. Now, the observed MRSI signal can be rewritten as s(kx;ky;t) = Z Z X i i(x;y)e 1 T 2i(x;y)te j!i(x;y)te j2 [kxx+kyy]dxdy (1.13) 1.3.2 Echo-Planar Imaging The standard MRSI protocol is to acquire a set of k-t samples on uniform grids and perform Fourier transforms on the acquired data to obtain the spatial distributions and the resonant spectra. A large number of coding steps are required for spectroscopic imaging, 9 Figure 1.4: (a) EPI pulse sequence (b) EPI trajectory and they are often too time-consuming for in-vivo imaging. Therefore, \high-speed imaging" has been introduced to the MRSI eld [17]. Due to developments in gradient hardware and processing methods, an imaging technique called echo-planar imaging (EPI) has become popular in MRSI [18, 19]. In this technique, one rf excitation is followed by a series of gradient reversals, thereby collecting a complete data set of one k-space image. Figure 1.4 shows a standard EPI trajectory and pulse sequence. Echo-planar imaging can be further divided into single-echo EPI and multi-echo EPI. In single-echo EPI, only one complete data set of a k-space image is observed in one excitation, while in multi-echo EPI, several k-space frames can be observed in one excitation, which will accelerate the imaging procedure further. As shown in Figure 1.4a, there is an echo-time (TE) parameter de ned as the time between the start of the rf pulse and the center of k-space. By changing the echo-time value, we will get a set of k-space images with di erent time o sets and therefore the 3-D (kx;ky;t) MRSI data. Figure 1.5 shows 21 k-space frames with di erent echo-time shifts. For instance, 10         Figure 1.5: EPI images with di erent echo-time values TE10 is the echo-time value describing the time between the start of the rf pulse and the center of 10th k-space image and the observed data at TE10 is represented as s(TE10) = Z Z ~ (x;y;TE10) expf j2 [kxx+kyy]gdxdy (1.14) In practice, a k-space image cannot be collected instantly in EPI. In other words, k- space images acquired in EPI are evolving images, not static images as in conventional MRI methodology. Each k-space sample is time dependent, so imaging model (1.13) can be to be updated for EPI with only one varible: s(t) = Z Z ~ (x;y;t) expf j2 [kx(t)x+ky(t)y]gdxdy (1.15) 11 Besides of its high speed, EPI is able to capture speci ed echo-time values and k-space trajectory during the acquisition process. Thus, in our study, we choose to work with variations of echo-planar imaging in MRSI acquisition. 12 Chapter 2 E cient Data Acquisition with Echo-Time Selection Long acquisition time is one of the major challenges of magnetic resonance spectro- scopic imaging. Standard MRSI collects data of both k-space and time domains to obtain spatial and spectral information of subjects. In order to maintain satisfactory spatial and spectral accuracy while reducing the acquisition time, several selection techniques have been developed to choose the best possible data subset to guarantee the reconstruction quality. Gao and Reeves [20, 21] and Plevritis and Macovski [22] investigated the optimal subset of k-space samples. The authors proposed to collect a limited number of k-space samples rather than all k-space data to reconstruct the image with a limited region of support, while keeping full data information in the time series. Conventional phase encoding imaging [17] is used in these methods, which requires one excitation to record an individual sample in k-space through the time domain. In this way, if accurate spatial details are needed, the acquisition time is still relatively long. Alternatively, within a comfortable acquisition time, only a limited number of k-space samples can be gathered for the subject, and the reconstructed images will have relatively poor spatial resolution. In some clinical applications, especially in locating abnormal biochemistry in speci c tissues, spatial details may be more useful than spectal information to indicate disease states. Therefore, other researchers emphasize the spatial distribution of key metabolites, acquiring high-resolution k-space coverage with limited chemical-shift encoding echo-time values. S. B. Reeder et al. [23] discussed optimization of echo-time spacing for maximum noise per- formance. However, their proposed method was only applied to a relatively small image assuming the spectral peaks are ideal at known locations; no T2=T2 [17] e ects are present; and each chemical metabolite only has one resonant frequency. These assumptions are not 13         Figure 2.1: (a) k-space sample selection, (b) echo-time sample selection always realistic for in vivo imaging. Figure 2.1 shows two types of selection methods de- scribed above: in the upper plot, data on each blue line is collected in one excitation, so the total acquisition time is proportional to the number of k-space sample locations; in the lower plot, data on each k-space is collected in one excitation, so the total acquisition time is proportional to the number of echo-time values. In MRSI, the spectral domain normally has greater sparsity than the spatial domain. Thus, we can reconstruct spectra with limited time series data more easily than reconstruct- ing the spatial domain with limited k-space data. In addition, as we mentioned in Chapter 1, EPI can easily capture speci ed echo-time values during the acquisition procedure. There- fore, in this chapter we focus on an e cient selection technique to choose the best echo-time values, which are then applied in EPI acquisition. 14 2.1 Echo-Time Selection Theory 2.1.1 Observation Model To understand the method for selecting echo-time values, consider a linear observation model in matrix notation: = Aq +e (2.1) where is a set of observations, q is the unknown vector of interest, A is a matrix that describes the transformation between the unknown q and the observations , and e represents zero-mean white Gaussian noise. The least-squares solution of q is: ^q = (AHA) 1AH (2.2) where the superscript H represents Hermitian transpose. Since e is independently and identically distributed, the mean squared error (MSE) in the reconstruction is proportional to Efk^q qkg= tr (AHA) 1 (2.3) where trfg means the trace of the inside matrix. This error criterion only depends on the transformation matrix A. Once we determine the A matrix, the MSE is xed and consequently the reconstruction quality of q is determined. For most types of MRSI, a desired time-varying image (x;y;t) can be described as (x;y;t) = MX i=1 ai(x;y) expf[ di(x;y) +j!i(x;y)]tg (2.4) where M is the maximum number of metabolite resonances contained in a voxel. The com- plex weight ai is proportional to the contribution from di erent resonances. The parameter 15 !i is the resonant frequency and di is the rate of damping of each resonance. After sampling, we can rewrite the basic image model at a given (x;y) location as (tn) = MX i=1 ai expf[ di +j!i]tng+e(tn) (2.5) where e(tn) is zero-mean, independent, identically distributed Gaussian noise, and tn is the echo time used to encode the nth EPI frame. The primary di erence of (2.5) from (2.1) is that the equations are nonlinear in the parameters of interest. To utilize the error criteria shown in (2.3), we propose to linearize (2.5) and use the linearized model to develop the error criterion. Note that the linearized model is only used in the selection procedure; during reconstruction, parameters of interest are still estimated using the more accurate nonlinear model. The linearized model is given by: (tn) = MX i=1 a0i exp [( d0i +j!0i)tn] + (ai a0i) exp [( d0i +j!0i)tn] + [( di +j!i) ( d0i +j!0i)]tna0i exp [( d0i +j!0i)tn] +e(tn) = MX i=1 ai exp [( d0i +j!0i)tn] + ( di +j!i)tna0i exp [( d0i +j!0i)tn] ( d0i +j!0i)tna0i exp [( d0i +j!0i)tn] +e(tn) = MX i=1 ai exp [( d0i +j!0i)tn] + ( di +j!i)tna0i exp [( d0i +j!0i)tn] +e (tn) (2.6) Experiments show that the the linearized model is adequate enough for the selection process (more details in Section 2.3 ). Here ai, di, !i are unknown signal parameters, while a0i, d0i, !0i are known from prior information. Equation (2.6) can be rewritten in matrix format like (2.1) with bi = di + j!i and b0i = d0i + j!0i, as the unknown parameters q = 16 [a1;a2;:::;aM;b1;b2;:::;bM]T. The transformation matrix is then given by A = 0 BB BB @ eb01t1 eb0Mt1 a01t1eb01t1 a0Mt1eb0Mt1 ... ... ... ... ... ... eb01tn eb0Mtn a01tneb01tn a0Mtneb0Mtn 1 CC CC A (2.7) At this point, A has only one set of variables, the choice of echo-time values tn (a0i;d0i;!0i are xed values from prior knowledge). Each row of the A matrix corresponds to a speci c tn. Hence, when we optimize the combination of rows in the A matrix, we are selecting the optimal combination of echo-time values that determines the reconstructed image quality. The MSE value is not merely a function of the number of rows removed from A. The speci c choice of echo-time values is more important than the number of echo times that are selected. In some cases fewer acquired echo times may result in better reconstruction. 2.1.2 Selection Algorithm Careful selection of the echo-time values to be applied in EPI acquisition, represented by rows of A, is required to minimize the MSE and hence optimize the reconstruction perfor- mance. Unfortunately, exhaustive search for the optimal row combination is prohibitive for a large matrix. A suboptimal sequential backward selection algorithm (SBS) is therefore used to optimize echo-time selection. We begin with all possible echo-time values and evaluate all possible choices of eliminating a single row from the current A, which corresponds to a speci c echo-time value. Let ri represent row i of A. If ri is eliminated from A, using the Sherman-Morrison matrix inversion formula [10], the modi ed (AHA) 1 is given by: ( ~AH ~A) 1 = (AHA) 1 + (A HA) 1rH i ri(A HA) 1 1 ri(AHA) 1rHi (2.8) Taking the trace of both sides and using the property that trfABg = trfBAg, we can formulate an updated error criterion as 17 trf( ~AH ~A) 1g= trf(AHA) 1g+ ri(A HA) 2rH i 1 ri(AHA) 1rHi (2.9) The rst term in (2.9) is the same for every choice of an eliminated row. The error criterion will yield the smallest increase by eliminating the row that minimizes the second term of (2.9). Then the least important row will be removed and A updated to ~A. This process is repeated until the desired number of echo-time values remains. This approach can achieve a satisfactory selection result in an e cient manner. Moreover, SBS has other attractive advantages for echo-time selection: The actual values of the parameters of interest do not enter into the linear criterion. Therefore, the optimization process can be accomplished before MRSI acquisition. The a0i values are scalars in the criterion and have only a small e ect on the optimiza- tion process. Therefore, even if the prior amplitude information is not accurate, the inaccuracy will only have a small e ect on the optimality of the echo-time selection. If the di erences among resonant frequencies are xed, a global frequency shift due to eld inhomogeneity will not a ect the optimality of the selected samples. The global frequency shift we mention here is not spatially global. It means that at a given (x;y) spatial location, due to eld inhomogeneity, the same frequency shift occurs for all resonant peaks. Compared to the global shifts, in general, the local variation for each resonance is small enough to ignore. Figure 2.2 shows the global shifts ! and local variations !1; !2; !3 at one spatial location. At di erent locations, the values of global shifts can vary. 18 w D w w 1w D 2w D 3w D w D w D Figure 2.2: (a) Global frequency shift (b) local variations Recall (2.7), if a global shift ! exists, the A matrix will update as A = 0 BB BB @ eb01t1ej !t1 eb0Mt1ej !t1 a01t1eb01t1ej !t1 a0Mt1eb0Mt1ej !t1 ... ... ... ... ... ... eb01tnej !tn eb0Mtnej !tn a01tneb01tnej !tn a0Mtneb0Mtnej !tn 1 CC CC A (2.10) let D = 0 BB BB BB B@ ej !t1 0 0 0 0 0 ej !t2 0 0 0 ... ... ... ... ... ... 0 0 0 0 ej !tn 1 CC CC CC CA then A = A D (2.11) 19 Notice that the updated A can be represented as the original A multiplying by a diagonal matrix with terms ej !ti. When we apply the updated A in (2.9), tr ( AH A) 1 = tr (DHAHAD) 1 = tr (ADDHAH) 1 = tr (AHA) 1 (2.12) These two diagonal matrices multiplied together are identity. Therefore, the mean squared error will not change with a global shift ! and so does selection result. Thus in MRSI, the echo-time selection is minimally a ected by global frequency shifts due to eld inhomogeneity. Sometimes, the number of remaining echo-time values is required to be relatively large. For EPI acquisition, if we only scan one k-space frame, corresponding to one echo-time value, within one excitation, the total acquisition time will still be long. By exploiting multi-echo EPI, one can acquire multiple k-space images within one excitation to further reduce the acquisition time. The basic timing scheme of Multi-Echo EPI is shown in Figure 2.3. The phase gradient Gy(t) is rewound to the original starting position after each k-space image to ensure identical k-space trajectories for all echo-images. Furthermore, the adjacent echo-time values used in one excitation should be just long enough to accommodate one k-space frame acquisition so that scanning time is not wasted between k-space frames. The time di erence between adjacent echoes cannot be smaller than one k-space frame acquisition time, or the k-space frames cannot be completely collected. To achieve this purpose, we modify the selection algorithm. Instead of removing one echo-time value at each elimination step, the modi ed algorithm removes multiple echo-time values simultaneously. We consider candidate echo-time values for individual k-space frames to be laid out uniformly in time as shown in Figure 2.4a. The square red dots represent the elimination pattern, while the round black dots represent other time choices to which the elimination array can be shifted. Figure 2.4b shows an example of elimination using a 20 Figure 2.3: Basic multi-echo EPI sequence periodic pattern to remove the least useful twelve echo-time points, four at a time. That is to say: the red, green and blue frames are not collected in the acquisition process; only black frames are actually acquired. In general, the selection algorithm allows the elimination pattern to be arbitrarily spaced. However, if the elimination pattern is unequally spaced, sometimes the remaining echo-time points cannot be e ciently used in multi-echo EPI. Therefore, we restrict the elimination pattern to be equally spaced and the echo interval to be a little longer than one k-space frame acquisition time. The MSE criterion for multi-echo selection can be generalized as trf( ~AH ~A) 1g= trf(AHA) 1g+ trf Rj(A HA) 2RH j 1 Rj(AHA) 1RHj g (2.13) 21               Figure 2.4: (a) Elimination pattern (b) periodic nonuniform elimination where Rj represents the jth combination of several rows in A, which corresponds to echo-time values that might be removed simultaneously in one elimination step. The combination that minimizes the term on the right is selected for elimination. 2.2 Phantom Experiment 2.2.1 Materials and Methods Data from a 1H phantom intended to verify the selection algorithm was acquired on a 4.7T 60cm-vertical-bore Varian primate MRI system (Varian Inc., Palo Alto, CA) at the University of Alabama at Birmingham, courtesy of Prof. Donald B. Twieg. The phantom, shown as Figure 2.5, was constructed from four cylindrical test tubes. One of them (right upper corner) held a methanol and water mixture with a 1:1 volume ratio. The other three tubes held ethanol and water mixtures with a 1:1 volume ratio. The four tubes were arranged in a rectangular con guration and placed inside a larger glass cylinder lled with water. 22 Figure 2.5: Phantom structure Figure 2.6 shows spectra of the phantom obtained with a conventional phase-encoding technique. Chemical shifts are relative to hydroxyl at 4.7T. In the methanol-water spectrum, the methyl (CH3) peak is 298 Hz down eld from hydroxyl (OH), while the signal contribu- tion of CH3 : OH 1 : 1:8. In the ethanol-water spectrum, the ethyl (CH2) peak is 241 Hz and the methyl peak is 728 Hz down eld from hydroxyl, while the signal contribution of CH2 : CH3 : OH 2 : 3 : 7:5. These parameters are treated as prior information which is used in the selection algorithm. 3-D MRSI data were observed by echo-planar imaging with a traditional bottom-up uniform trajectory as shown in Figure 1.4b. 500 echoes separated by 0.3788ms were used in the EPI acquisition (from 17.81ms to 206.83ms). These represent the distribution of echo times available for optimal echo-time selection. Note that in practice only the optimized echo-time EPI frames would be acquired, whereas here all EPI frames were acquired so that di erent sampling patterns could be studied. In this work, the time required for one EPI 23 Figure 2.6: 1H spectrum of methanol (upper plot) 1H spectrum of ethanol (lower plot) frame acquisition is 32.77ms. Therefore, if we want to utilize multi-echo EPI, the adjacent k- space images acquired within one excitation must have at least 32.77 ms echo-time di erence. In addition, the phase gradient needs some time to rewind to the original starting position after each k-space image. Thus, we chose the elimination pattern with a 34 ms interval. The spectral bandwidth was 2640Hz, and spatial resolution was 2mm 2mm with a 64 64 matrix. Since the EPI technique is used in this work, each k-space frame cannot be collected instantly. In other words, every k-space sample will have its own time delay even in the same frame. Therefore, we cannot simply use fast Fourier Transforms to reconstruct the images. Furthermore, an FFT cannot separate spatial information from di erent metabolite resonances. Therefore, we use conjugate gradients (CG) to nd the best matching parameters of the images. (details are discussed in Chapter 4). 24 0 50 100 150 200 250 050100150 Number of remaining echo?time values tr[(A H A) ?1 ] (proportional to MSE) multi?echo selection single?echo selection number =48 Figure 2.7: Number of remaining echo-time values vs. tr (AHA) 1 (which is proportional to MSE) Here all programming used MATLAB running on a single core of a workstation equipped with 2.4 GHz AMD Opteron 880 dual-core processors (@Xi Computer Corp). 2.2.2 Results In order to compare the speed and accuracy of single/multi-echo algorithms, we use both methods to eliminate echo-time values (from 500 available echoes). The multi-echo algorithm removes four values at a time. Figure 2.7 shows how the error criteria vary with the number of remaining echo-time values. We observe that the single-echo selection has better MSE performance. However, the MSE performance of the multi-echo selection is almost as good. In both cases, the MSE 25 increases signi cantly when the number of remaining echo-time values is less than 48 but increases slowly when the remaining number is greater than 80. Nyquist density is not neces- sary for our selection technique. Because the remaining echoes are nonuniformly spaced and we have prior information about the spectral peaks, the Nyquist criterion does not apply. The selection time (average of 100 repetitions) comparison of the two selection methods is shown in Table 2.1. For the same number of remaining echoes, the multi-echo selection is much faster. In both selection processes, eliminating all but 32 echo times can be accom- plished in less than 5 seconds. Even if real-time selection is necessary, the time required is feasible for clinical applications. However, su cient prior information may be available to do the selection o ine before the imaging session begins. Table 2.1: Selection Time Comparison (s) Number of remaining echo-time values 256 128 96 64 48 32 Single-echo selection time 3.0352 3.8885 4.0268 4.0325 4.0455 4.0657 Multi-echo selection time 0.4031 0.4990 0.5152 0.5248 0.5298 0.5354 Since MSE from the multi-echo technique is comparable to the single-echo techique and it has advantages in the acquisition process, we used the multi-echo algorithm (eliminat- ing echo-time values four at a time) to demonstrate the value of the selection technique. Forty-eight out of ve hundred echo-time values were selected for EPI acquisition, and these collected EPI data were used for image reconstruction. For comparison purposes, images reconstructed from 500 echoes were treated as ground truth. We also reconstructed images from 48 out of 500 randomly selected echoes and 48 out of 500 equally spaced echoes. Figure 2.8 shows reconstructed spatial details at di erent resonant frequencies. When water mixes with ethanol or methanol, their hydroxyls exchange very quickly. Therefore, we can treat all hydroxyls in the phantom as having the same resonant frequency. Hence, the hydroxyl resonance should exist at every spatial location of the cylindrical container and four tubes (Figure 2.8, rst column). Only ethanol has an ethyl part; when the spatial distribution 26 Figure 2.8: Spatial images reconstructed from 500 echo-time values ( rst row), 48 optimized echo-time values (second row), 48 randomly selected echo-time values (third row), 48 equally spaced echo-time values (fourth row). Each column from left to right represents spatial distributions of: hydroxyl, ethyl, methyl of methanol, methyl of ethanol. 27 of ethyl is shown, the resonant signals should only appear in the three ethanol tubes (Figure 2.8, second column). Both ethanol and methanol have a methyl part; however, the molecular environments around the methyl are very di erent from each other. Consequently, there are two resonances of methyl, one from ethanol and another one from methanol. In other words, when the spatial distribution of methyl from methanol is shown, there should be no signals in the ethanol tube (Figure 2.8, third column); when the spatial distribution of methyl from ethanol is shown, there should be no signals in the methanol tube (Figure 2.8, fourth column). Figure 2.9 shows reconstructed spectra at di erent spatial locations. In the water- ethanol tube, there should be three resonant peaks from hydroxyl, ethyl and methyl (Figure 2.9, upper plot); in the water-methanol tube, there should be two resonant peaks from hydroxyl and methyl (Figure 2.9, middle plot); in the big glass cylinder lled with water, there should be only one resonant peak from hydroxyl (Figure 2.9, lower plot). As shown in the spatial and spectal reconstruction results, both multi-echo SBS and full selection distinguished these four di erent resonances well. The spatial distribution and spectra display good agreement between the full data set and the optimized echo-time data set. In contrast, random selection and equally spaced selection yielded much greater reconstruction errors and artifacts, and they did not separate di erent resonant frequencies clearly. 2.3 Discussion In this chapter, we presented a novel technique to accelerate the data acquisition in MRSI. An e cient SBS method was developed to optimize the echo-time selection in the context of multi-echo EPI. We analyzed the computation time and reconstruction MSE of this selection algorithm and demonstrated its feasibility in MRSI acquisition. The phantom experiments demonstrate that spatial maps of di erent resonances of interest can be sepa- rated with a limited number of echo-time values if those echo times are chosen properly. In 28 Figure 2.9: Spectra in the tube of water and ethanol mixture (upper plot), in the tube of water and methanol mixture (middle plot), in the big glass cylinder of pure water (lower plot). Each plot from lower left to upper right from: 500 echo-time values, 48 optimized echo-time values, 48 randomly selected echo-time values, and 48 equally spaced echo-time values. Plots are o set for easier viewing. 29 addition, spatial maps of di erent spectra can also be reconstructed properly with a limited number of echo-time values. In our experiments, the volume of acquired data for the opti- mized selection was only 9.6% that of the full echo-time set. Thus, the proposed method has great potential to reduce overall acquisition time for obtaining spatial maps of spectral peaks of interest. As mentioned in the Materials and Methods subsection, linearization is an important step of the selection algorithm. If the prior values a0i;d0i;!0i have the same values as un- known signal parameters ai;di;!i, (2.6) is equal to (2.5). Therefore, accuracy of the prior information is a factor in the selection performance. In order to understand the impact of each kind of parameter on the selection performance, we did separate tests with param- eter o sets of magnitude, decay and frequency. Each study adheres to the following the experiment ow, as shown in Figure 2.10. Using original parameters instead of estimated parameters in the selection, the speci c combination of selected echo-time values can be totally di erent. However, that does not means the echo-time selection relying on the esti- mated parameters is not good. A more reasonable way to evaluate the selection performance is to examine how much the MSE increases if the parameters used in the selection process vary from the true parameters. Thus, a selection is made with parameters that are per- turbed from the true parameters, and the reconstruction MSE for that selection is evaluated assuming the true parameters. Meanwhile, the reference MSE for the selection from the true parameters is approximated using the linearized matrix equation around the true pa- rameters, which is asymptotically accurate as noise variance decreases. Here we treat the parameters reconstructed from 500 echo-time values as original values. For each kind of parameter, we tested 10 o set ranges. For each o set range, we repeated the experiment ow 1000 times with random o sets in that range to nd the minimum, maximum and average MSE increases. Figure 2.11 shows how the MSE increase was a ected by di erent parameter o sets. In Figure 2.11a, when estimated magnitudes are randomly perturbed by a maximum o set of 10% of the original magnitudes, the minimum MSE 30 ! " " # $ " # % ! # # ! " % # & ? ? $ $ ( ! " % ! " ) ! " ! * $ $ & ? + ! * $ $ & # # $ Figure 2.10: Experiment ow: Parameter o sets a ect the selection performance increase is 0%, the maximum MSE increase is 0.2% and the average MSE increase within 10% o set range is 0.01%. Similarly, when estimated magnitudes are randomly perturbed by a maximum of 15% to 50% from the original magnitudes, the maximum increase in MSE goes from 0.25% to 4.3%. The average increase in MSE goes from 0.04% to 0.37%. In Figure 2.11b, when estimated decays are randomly perturbed by a maximum of 10% to 50% from the original decays, the maximum increase in MSE goes from 0.7% to 9.3%, and the average increase in MSE goes from 0.23% to 1.6%. In Figure 2.11c, when estimated frequencies are randomly perturbed by a maximum of 1 to 10 Hz o set from the original frequencies, the maximum increase in MSE goes from 0.2% to 23%, and the average increase in MSE goes from 0.05% to 3.5%. From these observations, errors in prior knowledge of magnitudes have the smallest e ect on echo-time selection performance. If the magnitude error is limited to 50%, the MSE increase is less than 5%. Errors in prior knowledge of decay have a larger e ect on 31 Figure 2.11: MSE increase as a function of (a) magnitude o set (b) decay o set (c) frequency o set 32 echo-time selection. However, if we can keep the maximum decay error within 50%, the MSE increase can be limited to 10%. Compared to the resonant frequencies themselves, the frequency o set is relatively small. But the maximum MSE ampli es rapidly with the increase of frequency o set. Therefore, frequency accuracy is the most critical factor in echo-time selection. In practice, the chemical shifts of common metabolites used in MRSI are well known. In addition, normally the local variations of the resonant frequencies change very little, while frequency global shifts have no e ect on the optimality of the echo-time selection as shown in Section 2.1.2. For the proposed method, the MSE decreases when the number of selected echoes in- creases (Figure 2.7), which means the reconstruction quality improves with more data acqui- sition. For equally spaced selection and random selection, however, this property is not true. Figure 2.13 shows the MSE variation based on the number of remaining echo-time values. All the parameters come from prior information. For these sampling methods, MSE does not decrease monotonically when the number of echoes increases. In other words, equally spaced selection and random selection cannot guarantee better reconstruction results with more acquisition data. For the equal sampling pattern, we know that the frequency band- with is the reciprocal of time resolution [24]. When the bandwidth is larger than any of the frequency di erences, which means the time resolution (the echo-time interval) is ne enough, there will be no aliasing (Figure 2.12a ! < bw1). However, once the di erence between two resonant frequencies is an integral multiple of the spectral bandwidth (Figure 2.12b ! = 2bw2), reconstruction from equally spaced selection cannot separate resonant peaks due to aliasing. In our phantom experiment, the time di erence between adjacent echoes for equally spaced selection is 4.1668 ms and the spectral bandwidth is 239.9 Hz. The frequency di er- ence between hydroxyl and ethyl is 241 Hz, and the frequency di erence between hydroxyl and methyl from ethanol is 487 Hz. Both di erences are close to an integral multiple of 239.9 Hz, hence the reconstruction from equally spaced selection cannot separate these three 33 w D Figure 2.12: Spectrum of equally-spaced selection (a) from a small echo-time interval (b) from a large echo-time interval. resonances properly. For random selection, the echo-time values are not optimized, and the performance is-not surprisingly-random. Thus, random selection is not reliable for data acquisition either. In the selection procedure, a global frequency shift of resonannces has no e ect on the optimality of the selection. However, in the reconstruction procedure, large global shifts must be compensated by a eld map (This will be discussed in more detail in Chapter 5). Otherwise, spectral peaks may be initialized too far from the correct values, resulting in a failure to nd these values due to local minima. Consequently, magnitudes and decays of metabolite resonances will also be incorrect. In conclusion, we have developed an e cient echo-time selection technique, which can be easily implemented in MRSI acquisition. A lower-than-Nyquist density can be achieved with the proposed method. It greatly reduces the data requirements (and thus acquisition time) without sacri cing spatial and spectral accuracy in reconstruction. 34 Figure 2.13: MSE as a function of the number of remaining echo-time values (a) with equally spaced selection (b) with random selection 35 Chapter 3 Alternative Optimal Data Acquisitions In the previous chapter, we focused on selection techniques which choose the best echo- time values for echo-planar imaging acquisitions, while keeping high-resolution of k-space coverage. In fact, using EPI, a high-resolution k-space image cannot be acquired instantly at a single echo time; instead EPI requires a relatively long readout time. Every k-space sample has its own time delay even in the same k-space frame. Therefore, the echo-time selection might be optimal for a subset of k-space data but may not be the best choice for every k-space datum. Thus, in EPI, acquiring full coverage of k-space at each selected echo- time is somewhat imprecise. In addition, we cannot simply ignore the time required to scan each echo and use FFTs to reconstruct the images, because of this relative long acquisition time for high-resolution k-space frames. With these issues in mind, we introduce two alternative selection strategies: one is k-t selection, which use the sequential backward selection in both k-space and time domains; another one is echo-time selection that explicitly utilizes overlapped EPI, which can reduce the acquisition time for each k-space image and allows the possibility of using FFTs to shorten the reconstruction procedure. 3.1 Optimal k-t Selection Unlike the echo-time selection method, the k-t selection algorithm considers the k-space and time domain selection simultaneously. We try to choose the most useful k-t samples, which are then acquired for image reconstruction. As shown in Figure 3.1a, for standard echo-planar imaging, the k-space data are dis- tributed uniformly in the ky direction. This is because the pulse blips in the Gy direction 36 Figure 3.1: (a) Standard EPI with uniform k-space acquisition (b) Modi ed EPI with nonuni- form k-space acquisition 37 have the same amplitudes and durations, leading to equal spacing. Similarly, the k-space data are also distributed uniformly in the kx direction. The reason is the pulses in the Gx direction have the same amplitudes and durations. Here we can see that, for the same ky value, the kx value can vary from N=2 to N=2, which means the data acquisition in the kx direction is much faster than in the ky direction. Therefore, we still maintain full coverage in the kx direction and limit our k-t data selection only in ky-t domain. If we change the pulse amplitude or duration of Gy, the sample density in the ky direction changes. Normally, there is no data collection during the ky jump, so we keep the duration of each pulse in the Gy direction short. When the pulse amplitude increases, the trajectory will jump by a bigger step in the ky direction. When the pulse amplitude decreases, the trajectory will shrink to a smaller jump in the ky direction. In consequence, the nonuniform k-space frame can be acquired, as shown in Figure 3.1(b). Once the locations in the ky domain are determined, the pulse sequence for Gy can be implemented in echo-planar imaging. 3.1.1 k-t Selection Theory Figure 3.2 shows the main idea of k-t selection theory: all potential samples are supposed to lie on a uniform grid (black dots) in both ky and t domains. And the red dots represent the selected ky-t locations. As the gure shows, at some time points several ky samples are selected, while at some time points no ky sample is chosen. Similarly, at di erent ky locations, the number of selected time samples varies. On the other hand, at each selected ky-t location, we still fully sample the kx dimension. Similar to the echo-time selection, k-t selection also utilizes the sequential backward selection algorithm. Consider a linear observation model in matrix notation: = Aq +e (3.1) 38 . uni004Buni0079 uni004Buni0078 uni0054uni0069uni006Duni0065 Figure 3.2: ky-t data selection strategy where is a set of observations, q is the unknown vector of interest, A is the transformation matrix and e represents zero-mean white Gaussian noise. Then the mean squared error in the reconstruction is proportional to E k^q qk2 = tr (AHA) 1 (3.2) This error criterion only depends on the transformation matrix A. Once we determine the A matrix, the MSE is xed and consequently the reconstruction quality of q is determined. Since we fully sample the kx direction, an inverse FFT can be directly applied to the kx direction to get x spatial information. Here we only consider ky in k-space and y in the spatial domain. For most types of MRSI, a ky-t acquired image model can be described as: (ky;t) = NX y=1 MX i ai(y) expf[ di(y) +j!i(y)]tgexp [ j2 (kyy)] (3.3) 39 where M is the number of metabolite resonances contained in a voxel, N is the spatial coverage in the y direction, and ai, !i and di are the complex weight, resonant frequency and decay rate of each resonance. In order to utilize the error criterion shown in (3.2), we need to modify (3.3) and make it linear in the parameters of interest. After a linearized approximation, (3.3) becomes: (ky;t) = NX y=1 MX i=1 ai(y) exp [b0i(y)t] exp [ j2 (kyy)] +bi(y)ta0i(y) exp [b0i(y)t] exp [ j2 (kyy)] b0i(y)ta0i(y) exp [b0i(y)t] exp [ j2 (kyy)] +e(t) (3.4) where bi = di + !i. ai and bi are unknown signal parameters, while a0i and b0i are known from prior information. The noise term e(t) is zero-mean, white Gaussian noise. Equation (3.4) can be rewritten in matrix format, as the unknown parametersq = [a1(y1);:::;aM(y1);a1(y2);:::; aM(y2);:::;a1(yN); :::;aM(yN);b1(y1);:::;bM(y1);b1(y2);:::;bM(y2);:::;b1(yN);:::;bM(yN)]T. The transforma- tion matrix is then given by A = 0 BB BB BB BB BB BB BB BB B@ P1:M;y1(ky(1);t(1)) P1:M;yN(ky(1);t(1)) Q1:M;y1(ky(1);t(1)) Q1:M;yN(ky(1);t(1)) P1:M;y1(ky(2);t(1)) P1:M;yN(ky(2);t(1)) Q1:M;y1(ky(2);t(1)) Q1:M;yN(ky(2);t(1)) ... ... ... ... ... ... P1:M;y1(ky(N);t(1)) P1:M;yN(ky(N);t(1)) Q1:M;y1(ky(N);t(1)) Q1:M;yN(ky(N);t(1)) ... ... ... ... ... ... ... ... ... ... ... ... P1:M;y1(ky(N);t(n)) P1:M;yN(ky(N);t(n)) Q1:M;y1(ky(N);t(n)) Q1:M;NyN(ky(N);t(n)) 1 CC CC CC CC CC CC CC CC CA (3.5) 40 where P1:M;yl(ky(j);t(h)) = [eb01t(h)e j2 ky(j) l; ;eb0Mt(h)e j2 ky(j) l] Q1:M;yl(ky(j);t(h)) = [t(h)a01eb01t(h)e j2 ky(j) l; ;t(h)a0Meb0Mt(h)e j2 ky(j) l] At this point, A has two sets of variables, the choice of ky and t. Each row of A corresponds to a speci c (ky;t) sample. Hence, optimizing the combination of rows in A means selecting the optimal combination of (ky;t) samples which determines the quality of the reconstructed images. The MSE value is not merely a function of the number of rows remaining in A. The speci c choice of (ky;t) samples is more important than the number of samples that are selected. Compared to individual echo-time selection, the transformation matrix is much bigger in (ky;t) data selection. Thus, even using the sequential technique, the k-t selection requires longer computation time. The good thing is, in many clinical applications, su cient prior information may be available to do the selection o ine before the imaging session begins. In the selection procedure, we begin with all possible (ky;t) samples and evaluate all possible choices of eliminating a single row from the currentA, which corresponds to a speci c (ky;t) sample. Then the least important row, which causes the smallest MSE increase, will be removed from A. This process is repeated until the desired number of (ky;t) samples remain. Experiments shows that if there is no constraint in k-t selection procedure, the remaining (ky;t) samples scatter in a seemingly random pattern. 3.1.2 k-t Selection with Constraint As we mentioned at the beginning of this chapter, using EPI, the entire k-space cannot be acquired instantly at a single echo time; every k-space sample has its own time delay even in the same k-space frame. Therefore, we can modify our echo-time selection into the ky-t domain in a more accurate way, as shown in Figure 3.3a. If time points 1, 9, 17 and 25 41     Figure 3.3: Trajectory of (a) multi-echo EPI selection (b) k-t selection without constraint are selected in echo-time selection, this implies that all samples along the four red lines are chosen in k-t selection. In Figure 3.3, we have 8 36 potential (ky;t) samples in total (black dots), and our goal is to select the optimal 96 of them. Utilizing muti-echo time selection (acquiring 4 k-space images in one shot), only three excitations are needed to acquire the 96 (ky;t) samples. If non-constrained k-t selection is used, the acquired ky-t samples have no speci c pattern as shown in Figure 3.3b. At time point 1, ve ky samples are selected (red dot, blue dot, green dot, yellow dot and purple dot); at time point 2, three ky samples are selected; at time point 4, four ky samples are selected; at time point 6, two ky samples are selected; at time point 13, no ky sample is selected; and at time point 15, only one ky sample is selected. For one excitation, we can only pick one ky sample at one time point. Therefore, we need at least ve excitations (shown as ve color trajectories) to acquire 96 (ky;t) data. None of these trajectories can collect a selected point at every available sample time, so the overall scheme 42 is ine cient. For example, the red trjactory misses ky samples at time point 13 and point 25. Other color trajectories miss more ky samples during one excitation. Consequently, even though this k-t selection can achieve smaller MSE with the same amount of data (details in Experiment section) than echo-time selection, it still somewhat wastes acquisition time. In order to avoid wasting acquisition time, we update the k-t selection technique with a constraint. Assume the the potential range of the ky-t domain is m n and the total number of desired (ky;t) samples is N. In order to get fully e cient acquisition in one excitation, at each time point one ky sample must be selected. Then q = N=n excitations are needed to acquire the total N samples. Figure 3.4 shows an overview of the updated method: (1) create an accumulator for each time point and set the initial value to zero; (2) use the SBS technique to eleminate the least useful (ky;t) sample, which make the smallest increment of MSE; (3) determine the time point (i) of the selected sample and increment the ith accumulator by one; (4) compare the value of accumulator to (m q), if it is smaller than (m q); then repeat Step 2 and 3; if it is larger than (m q), shrink the potential k-t selection range to exclude all samples at time point i; (5) repeat Steps 2, 3 and 4 until the desired number of (ky;t) samples remain. Of course, the k-t selection with constraint has less freedom in selection, but it greatly enhance the acquisition e ciency compared to the original k-t selection method. More comparison details are shown in the experiment section. 3.2 Echo-Time Selection with Overlapped EPI For the k-t selection with or without constraint, the pulse sequence design would be a challenge, because it requires di erent k-space trajectories for each excitation. Hence, we introduce another selection method: echo-time selection with overlapped EPI, which is still a form of time domain selection and much easier in the pulse sequence design. In EPI, every k-space sample has its own time delay within the same k-space frame. A high-resolution k-space image (including many k-space samples) requires long scanning time. Therefore, the original echo-time selection technique is somewhat imprecise because it 43 ! " # $ % & ? ( ) ) " ) # ) % & $ * ) ( ( $ + & ) + & Figure 3.4: Flow of the updated method 44 assumes full coverage of k-space at each selected echo-time. However, if we can acquire the entire k-space in a negligible period, then the echo-time values we choose would be optimal for more samples of the selected k-space frame. For traditional echo-planar imaging, the observed data are distributed uniformly in the whole k-space domain as shown in Figure 3.5a. TR1 is the total readout time for a whole k-space image. If we want to decrease the readout time, we can simply reduce the resolution either in the kx domain, ky domain or both. Since for the same ky value, the kx value varies (data acquisition in the kx direction is much faster than in the ky direction), we here only consider reducing the resolution in the ky direction. When we double the pulse amplitude in the Gy direction shown in Figure 3.5b , the jump size in the ky direction also doubles and the sample density in the ky direction becomes half of the original one. Scanning the same k-space range, the total number of k-space samples in this frame will be halved and the readout time, TR2, is consequently half of TR1. On the other hand, in some clinical applications, especially in locating abnormal biochemistry in speci c tissues, high-resolution in k-space is very important. We want to shrink the readout time but without sacri cing the spatial resolution. Hence, in the next excitation, another half-resolution k-space image is acquired at the same echo-time as shown in 3.5c. The only thing di erent between these two low-resolution k-space frames is the initial setup of Gy: GY1 is larger than GY2, so there is an o set in the ky direction between the blue trajectory and red trajectory. In consequence, we can combine these two images into one with full k-space coverage, and the image readout time shrinks to half. For one full-coverage k-space image, the orignal EPI requires one excitation and there is one echo-time value to represent the whole k-space image with TR1 readout time; while the overlapped EPI requires two excitations and also only one echo-time value to represent this k-space image with smaller (TR2) readout time. That is, the original EPI is more e cient in data acquisition, while the overlapped EPI is more accurate in echo selection. 45 Figure 3.5: (a) Standard EPI with uniform k-space acquisition (b,c) Modi ed EPI with coarse ky resolution 46 In order to increase the e ciency of overlapped EPI acquisition, we try to observe more echoes during one excitation than in original EPI. For regular multi-echo EPI, as shown in 3.6a, four full-coverage k-space frames (four di erent echoes) are acquired in one excitation. These four images are distributed equally in the time domain (= TR1) for high e ciency of data acquisition. So we choose these four echoes as a set in the echo-time selection process. During the next excitation, another set (four full-coverage) of k-space frames are acquired with an optimal time shift, which is calculated in SBS echo-time selection. For overlapped multi-echo EPI, as shown in 3.6b, eight half-coverage k-space frames (eight di erent echoes) are acquired in one excitation. The time interval of these echoes is TR2 for data acquisition e ciency. Here we choose these eight echoes as a set in the echo-time selection process. During the next excitation, another eight half-coverage k-space frames are acquired using the same echo time values. For one excitation process, both selection methods acquire the same amount of k-space data. For the same acqusition e ciency, regular EPI has more exibility in choosing echo- time values, while echo time values selected for the overlapped EPI are more accurate in representing one k-space image. Because the k-space images acquired in overlapped EPI have a shorter readout time, FFTs may be su cient for reconstructing spatial information. Similarly, if we want to reduce the readout time further for one k-space image, we can reduce the resolution of ky more but with even less exibility in echo-time selection. 3.3 Experiments and Comparisons 3.3.1 Materials and Methods To compare the performance of di erent selection methods, 3-D MRSI data were ob- served by conventional phase-encoding imaging. During one excitation, data is acquired uniformly along the time domain at one k-space location at the red-dot line shown in Figure 3.7. For the next excitation, the pulse sequence changes the k-space location and rewinds to the original time point to start acquiring data at the blue-dot line shown in Figure 3.7. 47   a0a1 a2   a0a1 a2   a0a1 a3   a0a1 a3 Figure 3.6: (a) Traditional multi-echo EPI with uniform k-space acquisition for two excitation (b) Modi ed multi-echo EPI with coarse ky resolutions for two excitation     Figure 3.7: Data distribution of conventional phase-encoding imaging 48 Consequently, after a number of excitations, at each time point, we can get a full k-space image, in which the k-space date is distributed uniformly. In our experiment, 1024 time points separated by 0.5ms and 64 uniformly located ky samples acquired. These represent the distribution of (ky;t) samples available for optimal selection. Note that the phase-encoding imaging utilized, which has uniform distribution in ky and time domain, was used here for acquisition so that di erent sampling patterns could be studied under the same experiment conditions. In practice echo planar imaging would be used in the acquisition for all of these selection methods. The 1H phantom was acquired in a 4.7T 60cm-vertical-bore Varian primate MRI sys- tem (Varian Inc., Palo Alto, CA) at the University of Alabama at Birmingham. The phan- tom, shown in Figure 3.8, was constructed from four cylindrical tubes. Three of them held methanol and water mixture with a 1:4 volume ratio, while the other tube (left lower corner) and the large glass cylinder were lled with pure water. Chemical shifts are relative to hy- droxyl at 4.7T. In the methanol-water spectrum, the methyl (CH3) peak is 298Hz down eld from hydroxyl (OH), while the signal contribution of CH3 : OH 1 : 6:2. These parame- ters are treated as prior information, which is used in the selection algorithms. The spectral bandwidth was 2000Hz, and spatial resolution was 2mm 2mm with a 64 64 matrix. 3.3.2 Results Table 3.1: Selection Time(s) and MSE Comparison Selection time MSE k-t selection with constraint 7,565.4 8.8634 (3072 remaining (ky;t) samples out of 65,536) k-t selection without constraint 13,565.4 8.2845 (3072 remaining (ky;t) samples out of 65,536) Regular multi-echo time selection 4.413 11.3561 (48 remaining t samples out of 1024, eliminate 4 at a time) Overlapped multi-echo time selection (eight frames in one excitation) 1.096 12.6538 (48 remaining t samples out of 1024, eliminate 8 at a time) 49 Figure 3.8: Phantom structure Table 3.1 shows the comparisons of the selection time and mean-square errors of four selection methods mentioned in this chapter: k-t selection with constraint, k-t selection without constraint, regular muti-echo time selection (eliminate four least informative echoes at a time and at each of the remaining echoes, one full coverage k-space image is acquired) and overlapped multi-echo time selection (eliminate eight least informative echoes at a time and at each of the remaining echoes, two half coverage k-space image are acquired). The two k-t selection methods retained 3072 (ky;t) samples out of 65,536, while the two echo-time selections retained 48 time samples out of 1024. As discussed in Section 3.1.2, to compare these selection methods in a more resonable way, we should transfer the selection results from both echo-time selections to the k-t domain as shown in Figure 3.3a. For instance, if 48 echo-time values remain in echo-time selection, there should be 48 N (N is the resolution in ky domain, in this experiment N=64) samples left after k-t selection. Therefore, 48 echo-time values and 3072 (ky;t) samples represent the same amount of remaining k-space data. Two echo-time selection methods have a much faster selection process than the two k-t selection methods. This is because when we use the SBS algorithm to eliminate samples in these methods, the transformation matrix A (see (2.7) and (3.5)) has very di erent sizes. 50 Therefore, when we eliminate samples by calculating the MSE using tr (AHA) , the pro- cessing time will be much di erent. The k-t selection without constraint has an even slower selection process than the k-t selection with constraint. The reason is, for the k-t selection with constraint, once all the desired number of samples have been eliminated at one time point ti, we shrink the potential k-t selection range. This eliminates all samples at time point i, which reduces the size of transformation matrix A faster. The two k-t selection techinques have better MSE performance than the two echo- time selections. The MSE of the two k-t selection methods are comparable and so is the MSE of the two echo-time selection techniques. Nyquist density is not necessary for our selection techniques. Because the remaining samples from these selection methods are all nonuniformly spaced and we have prior information about the spectral peaks, the Nyquist criterion does not apply. Since MSE from k-t selection with constraint is comparable to k-t selection without constraint and the constrained version has higher scanning e ciency in the acquisition pro- cess, we used the k-t selection with constraint to demonstrate the value of the k-t selection technique. 3072 out of 65,536 (ky;t) samples were selected for modi ed EPI acquisition, and these collected EPI data were used for image reconstruction later. For regular muti-echo selection, 48 out of 1024 echo-time values were selected for regular EPI acquisition (elimi- nating echo-time values four at a time); for overlapped muti-echo selection, 48 out of 1024 echo-time values were selected for overlapped EPI acquisition (eliminating echo-time values eight at a time). Of course, when the remaining echo-time values are determined, we also need to transfer these echo-time values to (ky;t) samples to make the comparison fair. For comparison purposes, images reconstructed from full k-space data were treated as ground truth. We also reconstructed images from 48 out of 1024 equally spaced echoes. Figure 3.9 shows reconstructed spatial details at di erent resonant frequencies using di erent selection techniques. When water mixes with methanol, their hydroxyls exchange very quickly. Therefore, we can treat all hydroxyls in the phantom as having the same 51 Figure 3.9: Spatial images on the left picture are reconstructed from full k-space data ( rst row), 3072 optimized (ky;t) samples out of 65,536 (second row), 48 optimized echo-time values using regular echo-time selection (third row), 48 optimized echo-time values using overlapped echo-time selection (fourth row). The rst column represents spatial distributions of hydroxyl; the second column represents spatial distributions of methyl. Spatial images on the right picture are reconstructed from 48 equally spaced echo-time values. The upper image represents spatial distribution of hydroxyl and the lower image represents spatial distribution of methyl. 52 resonant frequency. Hence, a hydroxyl resonance should exist at every spatial location of the cylindrical container and four tubes (Figure 3.9, rst column of the left picture and the upper image of the right picture). Besides the hydroxyl part, methanol also has a methyl part. When the spatial distribution of methyl is shown, the resonant signals should only appear in the three methanol tubes: upper two tubes and right lower one (Figure 3.9, second column of the left picture and the lower image of the right picture). We observed that the images reconstructed from both full selection data and k-t selection data can distinguish these two resonances very well. Both echo-time selection methods have good reconstruction performance in hydroxyl images, although there are some small artifacts in the reconstructed methyl images. In contrast, equally spaced selection yielded much greater reconstruction errors and artifacts and it did not separate di erent resonant frequencies clearly. Figure 3.10 shows reconstructed spectra at di erent spatial locations. In the big glass cylinder and left lower tube lled with water, there should be only one resonant peak from hydroxyl (Figure 3.10, left plot). However, the reconstructed spectra from the two echo- time selection methods also have tiny spikes at methyl resonant frequency, which should not appear. Clearly the spectra reconstructed from k-t selection displays better agreement with the full data set than other selection methods. In the other three water-methanol tubes, there should be two resonant peaks from hydroxyl and methyl (Figure 3.10, right plot). All selection techniques except equally spaced echoes show good reconstruction results. 3.3.3 Conclusion This chapter introduces two new data selection methods. Compared to the original echo- time selection technique, each of them has advantages and limitations. k-t selection expands the selection range from only the time domain to both k-space and time domains. It makes the selection process more accurate, and the reconstructed images from selected samples have better quality. However, the k-t selection technique has more demanding requirements in pulse sequence design and requires longer computation time to select the most informative 53 Figure 3.10: Spectra in the big glass cylinder and lower left tube full of pure water (left plot), in the tubes of water and methanol mixture (right plot). Each plot from lower left to upper right from: full data set, 3072 optimized (ky;t) samples out of 65,536, 48 optimized echo-time values using regular echo-time selection, 48 optimized echo-time values using over- lapped echo-time selection, and 48 equally spaced echo-time values. Plots are o set for easier viewing. 54 (k;t) samples. On the other hand, the echo-time selection using overlapped EPI can be easily implemented in MSRI acquisition. The selection process can be done in seconds, so even if real-time selection is necessary, the selection time is feasible for many applications. Another advantage of overlapped echo-time selection is that the acquired k-space images has relative short readout time. In this case, FFTs might be good for spatial reconstruction, which leads to a much faster reconstruction procedure. From our experimental data, we can also observe that the reconstruction quality is comparable to the original echo-time selection but is not as good as the k-t selection method. In summary, choosing the selection methods depends on the application demands. If the application requires real-time selection, either the original echo-time selection or echo-time selection with overlapped EPI would be good. In addition, if a fast reconstruction process is desired, echo-time selection with overlapped EPI might be the best choice. In another aspect, if su cient prior information is available and selection o ine before the imaging session is acceptable, the k-t selection method would be preferred because it yields smaller reconstruction errors and artifacts. 55 Chapter 4 Fast Reconstruction for Nonlinear Model In our research work, echo-planar imaging is used for image acquisition, which allows us to obtain single/multiple k-space image(s) with one excitation. EPI models each datum as a sample in the (k;t) domains rather than a non-time-varying k-space image. In other words, every k-space datum has its own time delay (leading to local decay and phase variations) even in the same frame. Therefore, fast Fourier Transforms (FFTs) directly used in k-space can cause artifacts in the reconstruction. Additionally, for our proposed acquisition methods, the observed k-space data are selected based on sequential backward selection algorithm and normally are sampled nonuniformly in the time domain. Therefore, an FFT is not appropriate to separate spectrum information from di erent metabolite resonances either. In this chapter, we present a fast iterative reconstruction method by applying polynomial approximations to the exponential time functionand and then utilizing FFTs to nd the best matching parameters of the images. Simulation and phantom experiments demonstrate that our approach achieves satisfactory reconstruction quality while largely reducing computation time. 4.1 Introduction to Iterative Reconstruction Method As mentioned in Chapter 2, the observations of MRSI from echo-planar imaging can be represented as s(t) = Z Z (x;y;t) expf j2 [kx(t)x+ky(t)y]gdxdy (4.1) 56 and the discrete version of (4.1) is given by s(tn) = X x X y (x;y;tn) expf j2 [kx(tn)x+ky(tn)y]g (4.2) in which the time-dependent spatial distribution (x;y;tn) is (x;y;tn) = MX i=1 ai(x;y) expf[ di(x;y) +j!i(x;y)]tng (4.3) where M is the number of metabolite resonances contained in a voxel, ai(x;y), di(x;y) and !i(x;y) are signal parameters we want to estimate from the observed signal s(tn) by solving (4.2) and tn is the time point where we observe k-space data. From (4.2) and (4.3), one can understand EPI observation in another way: Fourier Transforms are taken at each time point tn for (x;y;tn), and only one sample at [kx(tn);kx(tn)] in each k-space frame is selected instead of the whole k-space image. Therefore, if the image size is N N, for traditional MRI, we only need to perform one FFT to get one k-space image, while in EPI we need N2 FFTs to get one k-space image. For conventional phase-encoding MRSI, the spatial information (x;y;tn) at each spe- ci c tn can be reconstructed by taking the inverse Fourier Transform of k-space at tn, because the trajectory within one k-space frame observes a single static image at the same time. The trajectory of EPI, however, samples the k-space of di erent spatial images that vary by an exponential time function exp ( d+j!)tn. Because the observation model is nonlinear in the parameters of interest and does not have the direct structure of an FFT, an iterative method is needed to reconstruct the spatial information. Besides the spatial details, the spectral information is also important in spectroscopic imaging. For conventional phase-encoding MRSI, after taking the inverse Fourier transform of k-space at each tn, one can do another direct FFT at every (x,y) location along the time domain to get the spectral information at each voxel. But for EPI acquisition with our proposed acquisition methods, the observed data normally are sampled nonuniformly in the 57 time domain and the average density is much lower that the Nyquist density. Therefore, one also needs to use an iterative method to reconstruct the spectral information. The most common way to solve this problem is to use an iterative technique to recon- struct the image parameters (the spatial distribution a, the decay rate d and the resonant frequency !) instead of reconstructing spatial images and spectra. We propose to use an iterative conjugate-gradients (GC) algorithm to minimize the cost function: C(p) = X tn jstn ^stn(p)j2 (4.4) with respect to p, where p represents the image parameters: a, d and !. Here stn is the observed data and ^stn is the estimated data from (4.2). The regular conjugate gradient algorithm used to reconstruct p is processed as follows [12]: 1. Initialize p = p0 2. Calculate the gradients of parameters g0 = rpC(p0). If g0 = 0, stop, else set p0 = g0 3. Line search k = arg min 0C(pk + k pk). Find the optimal k to minimize the cost function with updated parameters (pk + k pk) 4. Let pk+1 = pk + k pk 5. Calculate the gradients of parameters gk+1 = rpC(pk+1). If gk+1 = 0, stop 6. k = (gk+1)Tgk+1(g k)Tgk 7. Set pk+1 = gk+1 + k pk 8. Let k = k + 1; repeat from Step 3 to Step 8 58 In this procedure, Step 3 to Step 8 need to repeat as many times to achieve the optimal results. The line search calculation in Step 3 and gradients calculation in Step 5 are the two main computational parts of the operation. 4.1.1 Parameter Gradient The computation of parameter gradientsrpC includes three parts: @C=@a, @C=@d and @C=@!. C(p) = X tn jstn ^stn(p)j2 = X tn [^stn(p) stn][^stn(p) stn] @C @a = X tn @^stn @a [^stn(p) stn] + [^st n(p) stn] @(^stn) @a (4.5) Let ftn = [^stn(p) stn] . Combine with (4.2) plus (4.3), (4.5) can be modi ed as @C @a = X tn [@^stn@a ftn +f tn@(^stn) @a ] = 2Ref X tn @^stn @a ftng = 2Re (X tn e( d+j!)tne j2 [kx(tn)x+ky(tn)y]ftn ) (4.6) 59 Similarly, the other two types of gradients can be represented as @C @d = 2Ref X tn @^stn @d ftng = 2Re (X tn ae( d+j!)tn ( tn) e j2 [kx(tn)x+ky(tn)y]ftn ) (4.7) @C @! = 2Ref X tn @^stn @! ftng = 2Re (X tn ae( d+j!)tn (jtn) e j2 [kx(tn)x+ky(tn)y]ftn ) (4.8) For ith resonance at (x,y) spatial location, these parameter gradients are described as @C @ai(x;y) = 2Re (X tn e[ di(x;y)+j!i(x;y)]tne j2 [kx(tn)x+ky(tn)y]ftn ) (4.9) @C @di(x;y) = 2Re (X tn ai(x;y)e[ di(x;y)+j!i(x;y)]tn ( tn) e j2 [kx(tn)x+ky(tn)y]ftn ) (4.10) @C @!i(x;y) = 2Re (X tn ai(x;y)e[ di(x;y)+j!i(x;y)]tn (jtn) e j2 [kx(tn)x+ky(tn)y]ftn ) (4.11) 4.1.2 Line Search Line search is another critical step in the conjugate gradient method. After we calculate the gradients of parameters g and determine the directions of parameter optimization p, we need to gure out how large the adjustable step size should be to minimize the error cost (4.4) at the current iteration. Exhaustive search for the optimal is ine cient because we may need to try a large number of di erent values and for each value we try, we need to update the parameters and reevaluate (4.4), which is computationally intensive. Therefore, we desire to restrict the number of values that we evaluate to minimize the cost function C. We numerically analyzed the function of C( ) and found that the cost function along is a quadratic-like shape and the quadratic approximation can be used to reduce the computational complexity of the line search. 60 An overview of the quadratic approximation is described as follows (here \C0" represents the value of the cost function when = 0): 1. Initialize = 1 and update parameters as p1 = p + 1 p. 2. Evaluate the cost function C1 as shown in (4.4) with the updated parameters p1. 3. Update parameters as p2 = p + 2 1 p, and evaluate the cost function C2 based on updated parameters p2. 4. If C0 C1 and C2 >C1 , go to Step 7; else continue to the next step. 5. If C0 C1, set = 0:5 1 and C2 = C1; update p1 = p + 0:5 1 p; evaluate the cost function C1 with the updated parameters p1; and then go back to Step 4. Otherwise continue to the next step. 6. If C0 > C1 and C1 > C2, set = 2 1 and C1 = C2; update p2 = p + 2 2 1 p; evaluate the cost function C2 with the updated parameters p2; and then go back to Step 4. Otherwise stop with an error sign. 7. Approximate the cost function C as a quadratic function as C( ) = h 2 +k +l. 8. Utilize C(0) = C0, C( ) = C1 and C(2 ) = C2 to nd the values of h, k and l. 9. Calculate min = k=2h. Figure 4.1 shows the cost function varying with di erent values. We can see the original line search is very close to a quadratic approximation over the range of interest. At A ( = 0), B( = 4 10 7) and C( = 8 10 7) points, the two cost lines have intersections and these three points are used to create the quadratic shape and determine the min used in the current iteration. In theory, any three points on the C( ) function can determine a quadratic waveform shape. However, di erent points chosen from C( ) function will lead to a di erent quadratic approximation. Based on our numerical analysis, 61 Figure 4.1: Original vs. quadratic approximation of line search 62 Figure 4.2: Quadratic approximation of line search based on di erent intersections 63 we notice that if C point is lower than B point and B point is lower than A point (Figure 4.2, upper plot), or if C point is higher than B point and B point is higher than A point (Figure 4.2, lower plot), the quadratic approximation is not that close to the original line search. Therefore, in our line search process, we change the values to make A and C points higher than B point (try to make the B point close to the bottom of the original line search function) and then use the quadratic approximation to nd the min. Furthermore, in order to make this procedure more reliable, the value of C(p+ min p) will be compared to C(p). If C(p + min p) >C(p), it means the cost function along is not quadratic-like in the current searching range, and we need a full line search to be performed to nd the minimizer. In Chapter 2 and Chapter 3, we used the conjugate gradients algorithm as described above to reconstruct parameters from EPI acquisitions, and we found some drawbacks with this technique: Directly evaluating C(p) in the line search is computationally intensive. Directly evaluating rpC(p) is computationally intensive. The speed of convergence to optimal parameters is relatively slow. All of these problems lead to long reconstruction time. Therefore, we propose to develop an e cient parameter- tting algorithm that approximates the exponential time function as a polynomial model and utilizes FFTs to accelerate the reconstruction procedure. 4.2 Fast Reconstruction Algorithm Examining (4.2), (4.4), (4.6), (4.7) and (4.8), we notice that both data evaluations in the line search and gradient calculations include a exp[ 2 (kxx+kyy)] term, which is very similar to the form of a Fourier Transform. The reason we can not directly use FFTs in these equations is the exponential time function exp ( d +j!)tn. Unlike conventional MRI, which samples the k-space of a static image, EPI samples the k-space of di erent images that are 64 related by this exponential time function. Therefore, if we can utilize some approximations in the reconstruction model and replace the direct calculation of exp[ 2 (kxx + kyy)] by FFTs, the iteration process should accelerate. 4.2.1 Reconstruction Model Directly Utilizing FFTs In the previous chapters, we presented several acquisition methods for MRSI. All of them utilized EPI to shorten the acquisition duration. Echo-time selection methods use standard EPI, and the observed data is uniformly sampled in k-space while nonuniformly in the time domain. For k-t selection, a modi ed EPI technique is used and the observed data is only uniformly sampled in kx direction while nonuniformly sampled in the (ky;t) domain. In addition, no matter which EPI techniques we use, the data acquisition in the kx direction is very fast. For instance, in the phantom experiment of Chapter 2, the time needed to acquire a single kx line is only 0.5 ms. Therefore, in order to speed up the reconstruction process, we can ignore the time di erence among the data in the same kx line and use FFTs directly in the kx direction. The standard discrete observation model of MRSI from echo-planar imaging can be represented as s(tn) = X x X y MX i=1 ai(x;y) expf[ di(x;y) +j!i(x;y)]tngexpf j2 [kx(tn)x+ky(tn)y]g (4.12) After FFTs in the kx direction, the data used for reconstruction is given by s(xk;tn) = X y MX i=1 ai(xk;y) expf[ di(xk;y) +j!i(xk;y)]tngexp[ j2 ky(tn)y] (4.13) In consequence, when we use the same iterative conjugate-gradient idea mentioned in the rst section, the parameter gradients (4.5), (4.5) and (4.8) can be updated as 65 @C @a = 2Re (X tn e( d+j!)tne j2 ky(tn)yftn ) (4.14) @C @d = 2Re (X tn ae( d+j!)tn ( tn) e j2 ky(tn)yftn ) (4.15) @C @! = 2Re (X tn ae( d+j!)tn (jtn) e j2 ky(tn)yftn ) (4.16) Notice here, the reconstruction problem becomes parameter tting at one x location (xk), which means we need to repeat this iterative process at every x location. The parameters a, b and ! have no information in the x direction; they only vary with y location and di erent resonance i. After the whole iterative process, we combine the parameter set for eachxlocation to get the full parameter package. Even though we need to repeat the iterative process at each (x) location, the calculation time still decreases signi cantly. This is because, during each iteration, the sizes of calculation matrices in (4.13), (4.14), (4.15) and (4.16) shrink considerably, which is the main factor contributing to computational complexity. Recall that echo-time selection used overlapped EPI, the data acquisition for one k-space is relative fast and the k-space data is observed uniformly. So we may further accelerate the reconstruction process by using FFTs in both kx and ky domains. Consequently, the reconstruction problem becomes parameter tting at each (x;y) location. The parameter gradients (4.5), (4.5) and (4.8) can be updated as @C @a = 2Re (X tn e( d+j!)tnftn ) (4.17) @C @d = 2Re (X tn ae( d+j!)tn ( tn) ftn ) (4.18) @C @! = 2Re (X tn ae( d+j!)tn (jtn) ftn ) (4.19) 66 and s(tn) = PMi=1 ai expf[ di +j!i]tng used in line search calculations. The parameters a, b and ! have no information in the (x;y) direction; they only vary with di erent resonance i. To get the full parameter package, we need to repeat the iterative process at every (x;y) location. Sometimes, because the values of resonant frequency and decay rate are very large, even if the time di erence are very small, the local decay and phase di erence are still big. Therefore, ignoring the time di erences in k-space may not be accurate enough, and we need to nd a more accurate way to approximate the reconstruction model. 4.2.2 Reconstruction model based on polynomial approximations As we discussed in the previous sections, because of the exponential time function exp ( d +j!)tn, directly using FFTs in the gradient calculation and data evaluation in line search is sometimes not accurate enough. If we can separate the time variable tn from decay rate d and resonant frequency !, then a combination of FFTs can be used to speed up the reconstruction process. Tang and Reeves [25] proposed a polynomial approximation to separate desired image parameters and time in the exponential function. However, this method is limited to func- tional MRI with the SS-PARSE technique (Single-shot parameter assessment by retrieval from signal encoding). Unlike MRSI, functional MRI only includes one k-space image and only cares about one resonant frequency from water. Therefore, based on their methods, we develop a more general approximation algorithm for MRSI, which can accommodate more k-space images and more resonances. The common polynomial approximation of the exponential function is a Taylor series as shown below: ept = 1X n=0 (pt)n n! (4.20) 67 However, in order to get an acceptable approximation result, the number of the polynomial terms must be very large. Although we can separate the p and t, we still need to calculate FFTs many times, which means the Taylor approximation does not solve the problem of computational complexity. In the EPI technique, the acquisition time needed for one k-space image is relatively short (dozens of microseconds level). It is reasonable to assume that the exp[( d + j!)t] term is bounded within a limited small range. Here t represents the time points within one k-space frame. Based on this assumption, one can approximate the exponential function with limited polynomials and nd the coe cients in a least squares sense. Equation (4.20) can be updated as ept = cLpLtL +cL 1pL 1tL 1 +:::+c1pt+c0 (4.21) where ci is the polynomial coe cient for ascending powers and these coe cients are xed through the whole reconstruction process, not varying with time. The coe cients we look for need to provide a good approximation for a range of pt values, not just a single pt value. In other words, within a bounded range, no matter how p or t changes, each exp(pt) must be expressed with the same set of polynomial coe cients. The coe cient estimation process can be written in matrix form, if the length of vector pt is N argminc ept PTc 2 (4.22) where PT = 0 BB BB BB B@ (pt1)0 (pt1)1 (pt1)2 (pt1)L (pt2)0 (pt2)1 (pt2)2 (pt2)L ... ... ... ... ... (ptN)0 (ptN)1 (ptN)2 (ptN)L 1 CC CC CC CA (4.23) 68 0 1 2 3 4 5 6 7 8 9 ?1.5?1?0.5 00.51 Absolute Value of (p x t) Exponential Time Function Polynomial Approximation (real part) Original Time function (real part) Polynomial Approximation (imag Part) Original Time function (imag part) Figure 4.3: Original exponential function vs. polynomial approximation Figure 4.3 shows the precision of the polynomial approximation (order 17) for the ex- ponential time function exp(pt), where Re(p) varies within [ 10;0], Im(p) varies within [ 80 ;0] and t varies within [0;0:033]. Figure 4.4 shows the mean squared error (MSE) between the original exponential func- tion (ept, where p = 10+j80 and t2[0:0:033]) and the polynomial approximation varying with the polynomial orders. MSEL = X t e pt LX l=0 clpltl 2 69 10 11 12 13 14 15 16 17 18 19 20 10?2510?2010?1510?1010?5 Order Number (L) Mean Squared Error Figure 4.4: Absolute error varies with polynomial order 70 Apply the polynomial approximation to our exponential term exp[( d(x;y)+j!(x;y))tn] e[ d(x;y)+j!(x;y)]tn = e[ d(x;y)+j!0+j!1(x;y)]tn = ej!0tne[ d(x;y)+j!1(x;y)]tn = ej!0tne[ d(x;y)+j!1(x;y)]TEme[ d(x;y)+j!1(x;y)](tn TEm) = ej!0tne[ d(x;y)+j!1(x;y)]TEm X l cl[ d(x;y) +j!1(x;y)]l(tn TEm)l (4.24) where TEm is the echo-time (the time center) of the mth EPI frame. In this approximation, we make !(x;y) = !0 + !1(x;y) and !0 does not vary with spatial location (x;y). In our research, we assume that we know what kinds of resonances exist in our subject, so we can start our evaluation with the ideal resonant frequency of that resonance. For example, if we know the resonance is water, !0 here is equal to 0 Hz. Even if some spatial have no water resonance, we still let !0 = 0 everywhere while the magnitude of water resonance at these locations is 0. The term !1(x;y), which is caused by eld inhomogeneity and noise, can vary at di erent spatial locations. If there is more than one resonance in one voxel, there is a di erent ideal resonant frequency !0i for each resonance. And let pi(x;y) = di(x;y) +j!1i(x;y) for each resonance, the observed MRSI data (4.12) can be represented as s(tm) = X x X y X i ai(x;y)ejw0itmepi(x;y)tme j2 [kx(tm)x+ky(tm)y] = X x X y X i ai(x;y)ejw0itmepi(x;y)TEmepi(x;y)(tm TEm)e j2 [kx(tm)x+ky(tm)y] = X x X y X i ai(x;y)ejw0itmepi(x;y)TEm X l clpli(x;y)(tm TEm)le j2 [kx(tm)x+ky(tm)y] (4.25) = X i ejw0itm X l cl(tm TEm)l X x X y ai(x;y)epi(x;y)TEmpli(x;y)e j2 [kx(tm)x+ky(tm)y] (4.26) 71 where tm presents the time point in the mth k-space frame. If one wants to approximate the observed data in the (m+1)th k-space image, echo-time value TE needs to be changed to TEm+1. On the other hand, for observed data at di erent k-space frames, we still only need to do a polynomial approximation once to get a set of cl. This is because for di erent m values, the time variance (tm TEm) is bounded within the same small range (as in Chapter 2 phantom experiment, this time range is [-16ms 16ms]). In consequence, exp[pi(x;y)(tm TEm] can be approximated by the same set of cl. Similarly, for observed data of di erent resonance i, the same set of cl are used. The reason is the value of (tm TEm) is very small (dozens of microseconds level), even though the pi(x;y) varies with i (pi(x;y) = di(x;y) + j!1i(x;y) and nomally di(x;y) is no larger than 50 and j!1i(x;y)j is no larger than 150 ), the multiplication of [pi(x;y) (tm TEm)] is still bounded within a limited small range. Let Bil m(x;y) = ai(x;y)epi(x;y)TEmpli(x;y). Then (4.26) can be simpli ed as s(tm) = X i ejw0itm X l cl(tm TEm)l X x X y Bil m(x;y)e j2 [kx(tm)x+ky(tm)y] (4.27) We can calculate the inside term PxPyBil m(x;y)ej[kx(tm)x+ky(tm)y] of (4.27) with a 2D FFT, and the observed data s(tm) can be approximated with a linear combination of a relatively small number of FFTs. s(tm) = X i ejw0itm X l cl(tm TEm)l ~Bil m[k(tm)] (4.28) where ~Bil m[k(tm)] is the Fourier transforms of Bil m(x;y). Similarly, for the gradient calculation (4.9), (4.10) and (4.11) can be approximated as @C @ai(x;y) = 2Re (X t e[ di(x;y)+j!i(x;y)]te j2 [kx(t)x+ky(t)y]ft ) = 2Re (X m X tm ej!0itmepi(x;y)tme j2 [kx(tm)x+ky(tm)y]ftm ) 72 Here Ptm represents adding all observed time points in the mth k-space frame, while Pm represents adding all observed k-space frames. @C @ai(x;y) = 2Re (X m X tm ej!0itmepi(x;y)TEmepi(x;y)(tm TEm)e j2 [kx(tm)x+ky(tm)y]ftm ) = 2Re (X m X tm ej!0itmepi(x;y)TEm X l clpi(x;y)l(tm TEm)le j2 [kx(tm)x+ky(tm)y]ftm ) = 2Re (X l clpi(x;y)l X m epi(x;y)TEm X tm ej!0itm(tm TEm)le j2 [kx(tm)x+ky(tm)y]ftm ) (4.29) @C @di(x;y) = 2Re ( ai(x;y) X l clpi(x;y)l X m epi(x;y)TEm X tm ej!0itm(tm TEm)le j2 [kx(tm)x+ky(tm)y]( tm)ftm ) (4.30) @C @!i(x;y) = 2Re ( ai(x;y) X l clpi(x;y)l X m epi(x;y)TEm X tm ej!0itm(tm TEm)le j2 [kx(tm)x+ky(tm)y](jtm)ftm ) (4.31) The inside summations Ptm of (4.29), (4.30) and (4.31) can also be evaluated by FFTs. Therefore, the gradient calculations can be approximated with a linear combination of a relatively small number of FFTs. One thing we need to pay attention to: when FFTs are utilized, uniformly sampling in k-space is required. Therefore, k-t selection samples may not compatible with this polynomial approximation. 73 In this reconstruction protocol, even though an extra polynomial approximation step is added to the reconstruction process, the utilization of the FFTs and the reduction of the ma- trix size in the calculations give this algorithm high potential to accelerate the reconstruction speed. 4.3 Experiment Simulation and phantom experiments are presented in this section to validate the poly- nomial approximations used in the reconstruction process. Both simulation data and phan- tom data were acquired by the EPI technique. 500 echoes separated by 0.3788 ms were used in the EPI acquisition (from 17.81ms to 206.83ms). Utilizing SBS selection method, 48 out of 500 echo-time values were selected (48 k-space images) for the reconstruction procedure. For comparison, we used regular conjugate gradients, modi ed conjugate gradients (direct FFTs in kx direction) and polynomial approximation (order number=17) to reconstruct the spatial information. There are four di erent resonances: hydroxyl, ethyl, methyl from methanol and methyl from ethanol. The spectral bandwidth was 2640 Hz, and spatial resolution was 2 mm 2 mm with a 64 64 matrix. 4.3.1 Simulation experiment We evaluated the performance of di erent reconstruction techniques with a set of noisy simulation data. The rst column of Figure 4.5 shows the true spatial distributions at di erent resonant frequencies (0Hz-hydroxyl, 241Hz-ethyl, 298Hz-methyl from methanol and 728Hz-methyl from ethanol). This information is used to synthesize the decaying k-space images. We added 20 dB SNR white Gaussian noise to the synthesized k-space data and reconstructed the spatial information from these noisy data. The estimated parameters (a;d) we started with had 30% o set from the original ones and the estimated resonant frequencies started with 5 Hz o set. 74 Figure 4.5 shows reconstructed spatial details at di erent resonant frequencies. The reconstruction results from both regular conjugate gradient methods and polynomial ap- proximation display good agreement with the original spatial distributions. The MSEPCTs (mean squared errors in %) of the reconstructed parameters are listed in Table 4.1, 4.2 and 4.3. MSEPCT = P x;yjp(x;y) ^p(x;y)j 2 P x;yjp(x;y)j 2 100% Table 4.1: MSEPCT of reconstructed magnitude (%) hydroxyl ethyl methyl of methanol methyl of ethanol ordinary CG 3.72 13.73 15.90 7.56 polynomial approx 4.09 14.35 17.08 8.23 modi ed CG 5.36 18.89 22.49 10.25 Table 4.2: MSEPCT of reconstructed decay (%) hydroxyl ethyl methyl of methanol methyl of ethanol ordinary CG 0.86 14.90 11.63 3.48 modi ed CG 0.93 15.44 12.34 3.75 polynomial approx 0.98 21.57 15.45 5.63 Table 4.3: MSEPCT of reconstructed frequnecy (%) hydroxyl ethyl methyl of methanol methyl of ethanol ordinary CG 5.3957e-02 7.8051e-03 5.5608e-03 6.6698e-04 modi ed CG 5.5441e-02 7.8223e-03 5.7638e-03 6.6895e-04 polynomial approx 5.7344e-02 8.9078e-03 6.4547e-03 8.4903e-04 All of the reconstructions were implemented in MATLAB running on a single core of a workstation equipped with 2.4 GHz AMD Opteron 880 dual-core processors (@Xi Computer Corp). We compared their computation time for 100 iterations, shown in Table 4.4. One can see that using polynomial approximation in the reconstruction procedure, we can greatly shortens the calculation time without sacri cing the reconstruction quality. For the modi ed 75 Figure 4.5: Original simulated spatial images ( rst row), spatial images reconstructed from 48 optimized echo-time values with regular conjugate gradient method (second row), spatial images reconstructed from 48 optimized echo-time values with polynomial approximation method (third row), spatial images reconstructed from 48 optimized echo-time values with modi ed conjugate gradient method (fourth row). Each column from left to right represents spatial distributions of: hydroxyl, ethyl, methyl of methanol, methyl of ethanol. 76 conjugate gradient method (direct FFTs in kx direction), the reconstruction can also be accelerated and the reconstruction results are also acceptable, though not as good as the regular conjugate gradient method or polynomial approximation technique. Table 4.4: Comparison of reconstruction time for 100 iterations (s) ordinary CG modi ed CG polynomial approx 189.62k 16.23k 3.87k 4.3.2 Phantom Experiment The phantom data we used to demonstrate the value of the reconstruction techniques is the same set used in Chapter 2 phantom experiment. For comparison purposes, images reconstructed from 500 echoes with the regular conjugate gradient method were treated as ground truth. The reconstructed spatial results of the phantom experiment are compared in Figure 4.6. We observed that the images reconstructed from 48 optimized echo-time values with both the regular conjugate gradient method and polynomial approximation can distinguish di erent resonances very well. The reconstructed spatial distributions from these two methods yield good agreement with the ground truth. In contrast, the modi ed conjugate gradient method (direct FFTs in kx direction) also held acceptable reconstruction performance but there were some small artifacts, especially in the x direction. 4.3.3 Conclusion In this chapter, we focus on how to reduce the computational complexity to accelerate the reconstruction procedure of MRSI using the EPI technique. Sometimes, FFTs can be directly used in the reconstruction. In our case, we applied FFTs directly in the kx domain to obtain x direction information because the acquisition time in the kx direction is short enough to ignore. However, in many practical situations, this kind of approach may not be 77 Figure 4.6: Spatial images reconstructed from full k-space data with regular conjugate gra- dient method ( rst row), spatial images reconstructed from 48 optimized echo-time values with regular conjugate gradient method (second row), spatial images reconstructed from 48 optimized echo-time values with polynomial approximation method (third row), spatial images reconstructed from 48 optimized echo-time values with modi ed conjugate gradient method (fourth row). Each column from left to right represents spatial distributions of: hydroxyl, ethyl, methyl of methanol, methyl of ethanol. 78 accurate enough and the reconstruction time is still relatively long. We present a polynomial approximation used for the time exponential function in the observed MRSI model. In this case a linear combination of FFTs can be applied to reduce the computational complexity. Simulation and phantom experiments show that this method can greatly reduce the recon- struction time, and there is no signi cant di erence in reconstructed results compared to the regular conjugate gradient reconstruction method. However, this polynomial approxi- mation requires the observed data to be distributed uniformly in k-space; therefore, if the data is acquired by the k t selection method, we cannot use the polynomial approach to accelerate the reconstruction process. For future work, we intend to develop a more exible reconstruction algorithm that can accommodate more acquisition strategies and handle more parameter o sets. 79 Chapter 5 Optimization Techniques to Attenuate Frequency Local Minima E ect In MRSI, many reconstruction problems can be posed as image parameter estimation from a sequence of discrete, nite measurements. Unfortunately, frequency local minima are a common practical problem when signal model- tting procedures are implemented. Con- ventional, gradient-based optimization methods such as steepest descent, Newton?s method and conjugate gradients are subject to frequency local minima. As a result, convergence to the global minimum is not guaranteed. What is worse, if the estimated frequency is not correct, the estimated signal amplitude and decay rate will be inaccurate. 5.1 Origins of Local Minima Many signals used in MRSI are decaying sinusoids (cosines). Thus we only consider an exponential signal model de ned by amplitude, decay rate and frequency in our work: x(t) = ae( d+j!)t + (t) (5.1) where x(t) is the observed signal; a, d and ! are unknown signal parameters which need to be estimated; and (t) represents zero-mean, independent, identically distributed Gaussian noise. To understand the cause of local minima in the parameter tting procedure, consider a simple cost function with one resonant peak in 1-D: Z T 0 ae( d+j!)t a 0e( d0+j!0)t 2dt (5.2) 80 Here a0, d0 and !0 are original signal parameters; and T is the observed signal length. Let ~! = ! !0 and utilize ej = cos( ) +jsin( ), cos( ) = (ej +e j )2 and ej = 1. Then (5.2) can be updated as Z T 0 ej! 0t ae dtej~!t a 0e d0t 2dt = Z T 0 ae dt(cos ~!t+jsin ~!t) a 0e d0t 2dt = Z T 0 (ae dt cos ~!t a 0e d0t) +jae dt sin ~!t 2dt = Z T 0 (a2e 2dt +a20e 2d0t)dt Z T 0 2aa0e (d+d0)tcos~!tdt (5.3) = Z T 0 (a2e 2dt +a20e 2d0t)dt aa0 Z T 0 e (d+d0)te j~!tdt aa0 Z T 0 e (d+d0)tej~!tdt (5.4) We notice that the second and third terms of (5.4) are very similar to the form of a Fourier transform. However, in practice, T must be a nite value (none of the observed signals could be in nite length). In order to take advantage of the Fourier transform relationship, we extend the range of integration from (0;T) to ( 1;1) while adding a window function to the original expression. Therefore, (5.4) is modi ed as a02 2d0 (e 2d0T 1) + a 2 2d(e 2dT 1) aa0 Z 1 1 e (d+d0)tu(t)[u(t) u(t T)]e j~!tdt aa0 Z 1 1 e (d+d0)tu(t)[u(t) u(t T)]ej~!tdt (5.5) Now, the third term of (5.5) is the Fourier transform of e (d+d0)tu(t)[u(t) u(t T)], and the fourth term of (5.5) can be treated as a conjugate component of the third term. Now we use the fact that if two functions multiply in the time domain, their Fourier transform functions will convolve in the frequency domain. The extra u(t) here is used for easy Fourier transform 81 of the exponential term e (d+d0)t, which does not change the value of this expression. Let G(~!) =Ff[e (d+d0)tu(t)]g= 1(d 0 +d) +j~! H(~!) =Ff[u(t) u(t T)]g= Tsinc( ~!T2 ) e j~!T2 Then (5.5) will become a02 2d0 (e 2d0T 1) + a 2 2d(e 2dT 1) aa0 fG(~!) H(~!) + [G(~!) H(~!)] g = a0 2 2d0 (e 2d0T 1) + a 2 2d(e 2dT 1) 2aa0 <[G(~!) H(~!)] = a0 2 2d0 (e 2d0T 1) + a 2 2d(e 2dT 1) 2aa0 < [Tsinc( ~!T2 ) e j~!T2 ] 1(d 0 +d) +j~! (5.6) Examining (5.6), the rst term is a xed value based on the original signal parameters; the second term varies with estimated magnitude and decay does not oscillate; and the third term including sinc( ~!T2 ) has obvious local minima. Since the sinc( ~!T2 ) shape originates from the window function, the reason for local minima is the nite signal observations. Figure 5.1 shows the cost function as it varies with di erent parameter o sets. If there are only magnitude o set and decay o set, the cost function is non-oscillatory (Figure 5.1, lower plot). If frequency o set presents, local minima appear in the cost function (Figure 5.1, upper and middle plots). In practice, many signals in MRSI have more than one resonant peaks. In these cases, the general cost function is a combination of several oscillating functions and local minima still exist. For instance, if there are two resonant peaks, the cost function becomes Z T 0 X i=1;2 aie( di+j!i)t X i=1;2 a0ie( d0i+j!0i)t 2 dt (5.7) 82 Figure 5.1: Cost function varying with frequency and decay o sets (upper plot) frequency and magnitude o sets (middle plot) magnitude and decay o sets (lower plot) 83 After expansion, (5.7) updates as Z T 0 a 1e d1t cos!1t+a2e d2t cos!2t a01e d01t cos!01t a02e d02t cos!02t 2dt + Z T 0 a 1e d1t sin!1t+a2e d2t sin!2t a01e d01t sin!01t a02e d02t sin!02t 2dt = Z T 0 (a12e 2d1t +a22e 2d2t +a012e 2d01t +a022e 2d02t)dt + 2a1a2 Z T 0 e( d1 d2)t(cos!1tcos!2t+ sin!1tsin!2t)dt + 2a01a02 Z T 0 e( d01 d02)t(cos!01tcos!02t+ sin!01tsin!02t)dt 2a1a01 Z T 0 e( d1 d01)t(cos!1tcos!01t+ sin!1tsin!01t)dt 2a1a02 Z T 0 e( d1 d02)t(cos!1tcos!02t+ sin!1tsin!02t)dt 2a2a01 Z T 0 e( d2 d01)t(cos!2tcos!01t+ sin!2tsin!01t)dt 2a2a02 Z T 0 e( d2 d02)t(cos!2tcos!02t+ sin!2tsin!02t)dt (5.8) Utilizing the sum and di erence formulas of sin=cos cos!1tcos!2t+ sin!1tsin!2t = cos(!1 !2)t (5.8) becomes as Z T 0 (a12e 2d1t +a22e 2d2t +a012e 2d01t +a022e 2d02t)dt + Z T 0 [2a1a2e( d1 d2)t cos(!1 !2)t+ 2a01a02e( d01 d02)t cos(!01 !02)t]dt Z T 0 [2a1a01e( d1 d01)t cos(!1 !01)t+ 2a1a02e( d1 d02)t cos(!1 !02)t]dt Z T 0 [2a2a01e( d2 d01)t cos(!2 !01)t+ 2a2a02e( d2 d02)t cos(!2 !02)t]dt (5.9) 84 Normally, multi-resonant peaks in one MRSI signal have almost the same frequency o set, which means !1 !01 = !2 !02 = ~!, and !01 !02 = ! is known before observation. Therefore, (5.9) can be simpli ed as Z T 0 (a12e 2d1t +a22e 2d2t +a012e 2d01t +a022e 2d02t)dt + Z T 0 [2a1a2e( d1 d2)t + 2a01a02e( d01 d02)t] cos(!)tdt Z T 0 [2a1a01e( d1 d01)t + 2a2a02e( d2 d02)t] cos(~!)tdt Z T 0 2A1A02e( d2 d01)t cos(~! +!)tdt Z T 0 2A2A01e( d2 d01)t cos(~! !)tdt (5.10) Comparing (5.3) and (5.10), one can see that only the last three items including frequency o set(cos ~!t, cos(~! + !)t and cos(~! !)t) are the reasons for the oscillatory shape of the cost function. And the oscillating frequencies are all related to ~!, ~! +! and ~! !. 5.2 Optimization Protocol with Varying Decay In the previous section, we know that local minima come from the nite signal observa- tions (here, we focus on the case of one resonant peak), which is represented as the third term of (5.6). The estimated magnitude a is just a scalar in the third term of (5.6); it will not increase or decrease the number of local minima. Therefore, if we can decrease the oscillation from n [Tsinc( ~!T2 ) e j~!T2 ] 1(d0+d)+j~! o , the local minima e ect can be attenuated. 5.2.1 Theory Considering the various sectional views of Figure 5.1a and b, we can see that increasing the estimated a value does not change the shape of the cost function (Figure 5.2a), while increasing the estimated d value can spread out the main lobe of the cost function, and the local minima fade away. In the middle plot of Figure 5.2b, even though there are some local minima remaining, the main lobe that includes the global minimum is expanded. That 85 Figure 5.2: Cost function varying with frequency at (a) di erent magnitude o sets (b) di erent decay o sets means the potential of nding the global minimum is enhanced. On the other hand, when the decay rate decreases, the gradient of the main lobe becomes sharper, yielding faster convergence and more accurate estimates [12]. As a result, in order to converge to the global optimum accurately and e ciently, the following procedure is used. The initial estimated decay d is set large to widen the main lobe to ensure that the frequency o set value is inside the main lobe. Then conjugate gradient (GC) is used to modify the frequency parameters. Next, we decrease the decay value and repeat CG. The procedure becomes a narrower search, so that the global optimum is estimated more accurately. In addition, decreasing the decay value through iterations will also accelerate the convergence speed. In consequence, an accurate global optimum can be detected within a limited number of iterations. Here we need to be alert that in each 86 Figure 5.3: Overview of the proposed method iteration step, the frequency o set is always required to be inside the main lobe. Figure 5.3 shows an overview of the proposed method. 5.2.2 Experiments and Results To investigate the validity of the proposed method, three simulation experiments are presented. To focus on frequency local minima, we assume the signal magnitudes in the following experiments are one. Experiment 1: single resonant peak with uniform samples 87 Figure 5.4: Cost function with (a) decay=-2 (b) decay=-20 (c) decay=-50 x(t) = e( 1:3+j535 )t t2(0 : 0:18) (5.11) The cost function is C = Z 0:18 0 e( d+j!)t x(t) 2 (5.12) The observed x(t) has 64 uniform samples, d is the estimated decay and ! is the estimated frequency. We explore the frequency o sets (! 5352 ) from -50 Hz to 50 Hz. Figure 5.4 compares di erent estimated decays applied to the cost function: d = 2 (upper plot), d = 20 (middle plot) and d = 50 (bottom plot). If we initialize the estimated frequency 30 Hz away from the original: for the cost func- tion of Figure 5.4a, the reconstructed frequency is 28.5 Hz away from the original (converging to the nearest local minimum); while for the cost function of Figure 5.4b and 5.4c, the recon- structed frequency will be almost the same as the original (achieving the global minimum). 88 Table 5.1 shows the acceptable frequency o set range for cost functions with di erent decays and compares the number of iterations required for minimum convergence. Here acceptable frequency o set range indicates the region within which there are no frequency local minima, and the number of iterations indicates the convergence speed. Table 5.1: Comparison of o set range (Hz)and number of iterations Acceptable frequency Number of iterations for o set range converging from 30 Hz o set to a minimum Cost function as Figure 5.4a -4 $ 4 7 (Local minimum) Cost function as Figure 5.4b -45 $ 45 113 (Global minimum) Cost function as Figure 5.4c -45 $ 45 252 (Global minimum) Varying cost function as proposed method -45 $ 45 78 (Global minimum) Experiment 2: two resonant peaks with uniform samples x(t) = e( 1:3+j535 )t +e( 3+j580 )t t2(0 : 0:18) (5.13) The observed x(t) has 64 uniform samples, di is the estimated deacy and !i is the estimated frequency. We explore the frequency o sets from -50 Hz to 50 Hz. Figure 5.5 compares di erent estimated decays applied to the cost function: d = 2 (upper plot), d = 20 (middle plot) and d = 50 (bottom plot). Table 5.2 shows the acceptable frequency o set range for cost functions with di erent decays and compares the number of iterations required for minimum convergence. Experiment 3: one resonant peak with nonuniform samples Uniformly sampled signals were investigated in the previous two experiments, now we execute another experiment using the following signal, and the measured signal is selected 89 Figure 5.5: Cost function with (a) decay=-2 (b) decay=-20 (c) decay=-50 nonuniformly (64 samples) within (0;0:18): x(t) = e( 1:3+j535 )t (5.14) We explore the frequency o sets from -50Hz to 50Hz. Figure 5.6 compares di erent estimated decays applied to the cost function: di = 2 (upper plot), di = 20 (middle plot), di = 50 (bottom plot). Table 5.3 shows the acceptable frequency o set range for cost functions with di erent decays and compares the number of iterations required for minimum convergence. Comparing Experiment 1 and Experiment 3 (with the same signal parameters), one can see that the acceptable frequency o set range is di erent. Revising (5.4), Z T 0 (ae 2dt +a0e 2d0t)dt aa0 Z T 0 e (d+d0)te j~!tdt aa0 Z T 0 e (d+d0)tej~!tdt 90 Table 5.2: Comparison of o set range (Hz) and number of iterations Acceptable frequency Number of iterations for o set range converging from 25 Hz o set to a minimum Cost function as Figure 5.5a -4 $ 4 7 (Local minimum) Cost function as Figure 5.5b -26.5 $ 26.5 132 (Global minimum) Cost function as Figure 5.5c -27 $ 27 308 (Global minimum) Varying cost function as proposed method -27 $ 27 83 (Global minimum) Table 5.3: Comparison of o set range (Hz) and number of iterations Acceptable frequency Number of iterations for o set range converging from 15 Hz o set to a minimum Cost function as Figure 5.6a -3.5 $ 3.5 5 (Local minimum) Cost function as Figure 5.6b -15.5 $ 15.5 104 (Global minimum) Cost function as Figure 5.6c -16 $ 16 216 (Global minimum) Varying cost function as proposed method -16 $ 16 67 (Global minimum) if t is not sampled uniformly, we can not extended the range of integration from (0;T) to ( 1;1) by just adding a window function to the original expression. Besides the window function, there must be another sampling function to represent the sampling pattern. Now (5.4) can be modi ed as Z T 0 (ae 2dt +a0e 2d0t)dt 2aa0< Z 1 1 e (d+d0)tu(t)p(t)[u(t) u(t T)]e j~!tdt (5.15) where p(t) is the nonuniform sampling pattern. The second term of (5.15) becomes the Fourier transform ofe (d+d0)tu(t)p(t)[u(t) u(t T)], which is the convolution ofF[e (d+d0)tu(t)], F[u(t) u(t T)] andF[p(t)]. Normally, increasing the decay values to avoid local minima has better performance in the uniformly sampled case than in the nonuniformly sampled case. 91 Figure 5.6: Cost function with estimated (a) decay=-2 (b) decay=-20 (c)decay=-50 From these experiments, one can see that enlarging the estimated decay value is an easy way to extend the acceptable frequency o set range within which the global minimum can be detected by gradient methods. However, when we keep increasing the decay values (in these experiments, from -20 to -50), the acceptable frequency o set ranges rarely change. In other words, if the estimated decay exceeds some values, enlarging decay will not further help in solving local minimum problems. In addition, a cost function with large estimated decays converges slowly (needs more iterations). Therefore, we need to select the estimated decays very carefully so that we can fully utilize the acceptable frequency o set range while keeping a relatively sharp gradient of the cost function. 5.3 Optimization Protocol with Weighted Scalars Examing Figure 5.4, Figure 5.5 and Figure 5.6, we notice that when the estimated decays increase, the dynamic range of the cost function shrinks signi cantly. Especially 92 at the edge of the acceptable frequency o set range, the gradient of cost function is very small. In gradient-based method, if the gradient is very small, it will a ect the optimization performance and speed [12, 26]. Moreover, when we use varying estimated decays to extend the acceptable frequency o set range, there is no easy way to decide what is the estimated decay value to start with: if the estimated decay is not big enough, we may not nd the largest acceptable frequency o set range; but if the estimated decay is too big, it will sacri ce optimization speed and accurancy. Therefore, in this section, we introduce a new method: applying a sequence of weighted scalars to the signal, which could get the cost function having a better acceptable frequency o set range than the varying decay method, while avoiding these issues mentioned above. 5.3.1 Theory Figure 5.7 shows a typical cost function. Normally, it is reasonable to assume that the cost function C(~!) is symmetrical from ~! = 0, so here we focus on the positive side of the cost function. The solid waveform is the original shape of the cost function, the rst gradient of C(~!) in (0; ~!p) is always C0(~!) 0. Our goal is to change the shape of the cost function and let C0(~!) 0 in (0; ~!p1) as dashed line shown. In our case, the cost function is continuous, so is the gradient of the cost function. As shown in Figure 5.7, when the frequency o set = ~!p, then C0(~!p) = 0. In order to extend the acceptable frequency o set range to ~!p1, we need to make C0(~!p) > 0. In consequence, we utilize the following criterion: min ~!pX 0 [C0(~!)] (5.16) where [C0(~!)] = 8> >< >>: 1 C0(~!)+ ; if C 0(~!) 0: 1 C 0(~!); if C0(~!) < 0: 93 Figure 5.7: Typical cost function and is a very small positive number, which makes C0(~!) + not equal to zero. During each iteration step to minimize P~!p0 [C0(~!)], we need to update the cost function C(~!) and the rst gradient of the cost function C0(~!) with the varying weighted scalars. C(~!) = Z T 0 W(t) ae( d+j!)t a0e( d0+j!0)t 2dt = Z T 0 W(t)(ae 2dt +a0e 2d0t)dt Z T 0 W(t)2aa0e (d+d0)tcos~!tdt (5.17) C0(~!) = Z T 0 W(t)2aa0e (d+d0)ttsin~!tdt (5.18) Here W(t) is a function de ned by scalar weights that control the value of the function where the signal samples are observed. The following procedure is used to implement the proposed method: 1. Find rst positive ~! = ~!p, which makes C0(~!p) = 0. Hence, we can de ne the summa- tion range of criterion (5.16). 94 2. Update the cost function and the rst gradient of the cost function using (5.17) and (5.18) and the initial W(t) = 1. 3. Vary the values of W(t) to minimizeP [C0(~!)] in the range (0; ~!p), where ~!p is de ned in step 1. 4. Repeat step 1, 2 and 3 until the new ~!p changes only slightly from the previous one. 5.3.2 Experiments and Results In order to compare the performance between the weighted scalar method and varying decay method, three simulation experiments using the same orginal signal, are presented. Experiment 1: single resonant peak with uniform samples x(t) = e( 1:3+j535 )t t2(0 : 0:18) (5.19) The observed x(t) has 64 uniform samples. Exploring the frequency o sets from -50Hz to 50Hz, Figure 5.8 compares cost functions: without weighted scalars (upper plot), and with weighted scalars (bottom plot). Table 5.4 shows comparisons among the original cost function, the cost function applying varying decays and the cost function applying weighted scalars. Table 5.4: Comparison of o set range (Hz) and number of iterations Acceptable frequency Number of iterations for o set range converging from 30 Hz o set to a minimum Original cost function -4 $ 4 7 (Local minimum) Cost function applied varing decays -45 $ 45 78 (Global minimum) Cost function applied weighted scalars -50 $ 50 42 (Global minimum) 95 Figure 5.8: Cost function using estimated decay =-2 (a) without weighted scalars (b) with weighted scalars Experiment 2: two resonant peaks with uniform samples x(t) = e( 1:3+j535 )t +e( 3+j580 )t t2(0 : 0:18) (5.20) The observed x(t) has 64 uniform samples. We explore the frequency o sets from -50Hz to 50Hz. Figure 5.9 compares the cost functions with or without weighted scalars. Table 5.5 shows the acceptable frequency o set ranges and the number of iterations required for di erent cost functions. Experiment 3: one resonant peak with nonuniform samples x(t) = e( 1:3+j535 )t (5.21) The observed signal is selected nonuniformly (64 samples) within (0;0:18). Figure 5.10 compares the cost functions with or without weighted scalars for the frequency o sets from 96 Figure 5.9: Cost function using estimated decay =-2 (a) without weighted function (b) with weighted function -50Hz to 50Hz. Table 5.6 shows the acceptable frequency o set ranges and the number of iterations required for di erent cost functions. 5.4 Conclusion and Future Work In this chapter, we introduced two methods to alleviate local minima e ects. One is varying estimated decay during the reconstruction process; the other is optimizing a weighted function applied to the observed and estimated signal. Both methods intend to change the shape of the cost function and make the global minimum located in a wider main lobe. Based on the above experiments, one can see that enlarging the estimated decay value is an easy way to extend the acceptable frequency o set range. However, since it will signi cantly decrease the gradient of the cost function, the reconstruction might be not accurate and is relatively slow. On the other hand, applying the weighted scalars to the 97 Table 5.5: Comparison of o set range (Hz) and number of iterations Acceptable frequency Number of iterations for o set range converging from 25 Hz o set to a minimum Original cost function -4 $ 4 7 (Local minimum) Cost function applied varing decays -27 $ 27 83 (Global minimum) Cost function applied weighted scalars -33 $ 33 37 (Global minimum) Table 5.6: Comparison of o set range (Hz) and number of iterations Acceptable frequency Number of iterations for o set range converging from 15 Hz o set to a minimum Original cost function -3.5 $ 3.5 5(Local minimum) Cost function applied varing decays -16 $ 16 67 (Global minimum) Cost function applied weighted scalars -25 $ 25 26 (Global minimum) signal can also extend the acceptable frequency o set range, while avoiding problems with reconstruction speed and accuracy. In the reconstruction procedure, large global shifts due to eld inhomogeneity must be compensated by a eld map. Otherwise, spectral peaks may be initialized too far from the correct values, resulting in a failure to nd these values due to local minima. Future research will focus on making these methods are more robust to nonuniform measurements and accelerate the process of convergence. 98 Figure 5.10: Cost function using estimated decay =-2 (a) without weighted function (b) with weighted function 99 Chapter 6 Conclusion 6.1 Summary of the Contributions of This Dissertation In this dissertation, we studied several problems related to data acquisition and recon- struction in magnetic resonance spectroscopic imaging. To accelerate the acquisition time, we rst proposed an e cient acquisition scheme with echo-time selection. A criterion to estimate the mean square error in reconstruction was derived. Both single-echo SBS and multi-echo SBS were used to select a limited set of the most important echo-time values which were then applied to EPI acquisition. A lower- than-Nyquist density can be achieved with the proposed method. 1H phantom experiments demonstrate that our approach achieves similar results to standard MRSI while only using 9.6% as much acquired data. Due to the limitations of echo-time selection in data acquisition, two alternative data selection strategies were studied in Chapter 3. One is k-t selection, which extends the data selection in both k-space and time domains. Another one uses echo-time selection but with overlapped EPI, which can reduce the acquisition time for each k-space image. From the phantom experiment, we observe that k-t selection makes the selection process more accurate while requiring longer computation time to select the most informative data. In contrast, the selection process of echo-time selection using overlapped EPI can be done in seconds, but the reconstruction quality is not as good as k-t selection. Therefore, choosing which selection methods to use for acquisition depends on the application demands. In Chapter 4, fast reconstruction methods for the nonlinear imaging model were inves- tigated. We proposed to apply a polynomial approximation to the time exponential function 100 in the observed MRSI model and utilize a linear combination of FFTs to reduce the compu- tational complexity in reconstruction. Simulation and phantom experiments show that this technique can largely accelerate the reconstruction process without sacri cing reconstruction accuracy. The last part of this dissertation described two methods to alleviate frequency local minima e ects. One is varying estimated decay during the reconstruction process, which can be easily implemented but does not achieve very good performance in reconstruction e ciency and accuracy. Another advanced technique we developed is optimizing a weighted function applied to the observed and estimated signal. The experimental results show that the weighted function method can e ciently attenuate frequency local minima e ects, while avoiding the reconstruction speed and accuracy problems. 6.2 Future Works For the acquisition part, the future research will focus on how to simplify the selection process in the k-t domain while still keeping the selection accuracy. Since our proposed methods are validated by simulation and phantom experiments, if it is possible, we also plan to apply our optimized selection strategies to in vivo experiments. For the reconstruction part, we intend to develop a fast and more exible reconstruction algorithm, which can accommodate more acquisition strategies and handle more parameter o sets. 101 Bibliography [1] E. M. Haacke, R. W. Brown, and R. Venkatesan, Magnetic Resonance Imaging: Physical Principles and Sequence Design. Wiley-Liss, 1999. [2] 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, 1999. [3] J. Kurhanewicz, D. Vigneron, and S. Nelson, \Three-dimensional magnetic resonance spectroscopic imaging of brain and prostate cancer," Neoplasia, vol. 2, pp. 166{189, 2000. [4] N. Schu , F. Ezekiel, and etc., \Region and tissue di erences of metabolites in normally aged brain using multislice 1H magnetic resonance spectroscopic imaging," Magnetic Resonance in Medicine, vol. 45, pp. 899{907, 2001. [5] X. Wang, B. Wang, Z. Gao, J. Liu, and etc., \1H-MRSI of prostate cancer: The relation- ship between metabolite ratio and tumor proliferation," European Journal of Radiology, vol. 73, no. 2, pp. 345 { 351, 2010. [6] M. Fulham, A. Bizzi, M. Dietz, H. Shih, and etc., \Mapping of brain tumor metabolites with proton MR spectroscopic imaging: clinical relevance," Radiology, vol. 185, pp. 675{ 686, 1992. [7] M. Jacobs, P. B. DPhil, and etc., \Proton magnetic resonance spectroscopic imaging of human breast cancer: A preliminary study," Journal of Magnetic Resonance Imaging, vol. 19, pp. 68{75, 2003. [8] K. Zakian, K. Sircar, H. Hricak, and etc., \Correlation of proton MR spectroscopic imaging with gleason score based on step-section pathologic analysis after radical prosta- tectomy," Radiology, vol. 234, no. 3, pp. 804{814, 2005. [9] S. Hu, M. Lustigc, A. Chen, and etc., \Compressed sensing for resolution enhancement of hyperpolarized 13C yback 3D-MRSI," Journal of Magnetic Resonance, vol. 192, pp. 258{264, 2008. [10] S. J. Reeves and Z. Zhe, \Sequential algorithms for observation selection," IEEE Trans- actions Signal Processing, vol. 47, no. 1, pp. 123{132, 1999. [11] R. Pohmann, M. von Kienlin, and A. Haase, \Theoretical evaluation and comparison of fast chemical shift imaging methods," Journal of Magnetic Resonance, vol. 129, pp. 145{ 160, 1997. 102 [12] E. K. P. Chong and S. H. Zak, An Introduction to Optimization (Second Edition). Wiley-Interscience, 2001. [13] J. Hugg, K. Laxer, G. Matson, A. Maudsley, and etc., \Lateralization of human fo- cal epilepsy by 31P magnetic resonance spectroscopic imaging," Neurology, vol. 42, pp. 2011{2018, 1992. [14] J. Hugg, G. Matson, D. Twieg, and etc., \Phosphorus-31 MR spectroscopic imag- ing (MRSI) of normal and pathological human brains," Magnetic Resonance Imaging, vol. 10, pp. 227{243, 1992. [15] M. Swanson, D. Vigneron, and etc., \Proton HR-MAS spectroscopy and quantitative pathologic analysis of MRI/3D-MRSI-targeted postsurgical prostate tissues," Magnetic Resonance in Medicine, vol. 50, no. 5, pp. 944{954, 2003. [16] D. Shaw, Fourier Transform NMR Spectroscopy. New York: Elsevier Scienti c, 1971. [17] Z. H. Cho, J. P. Jones, and M. Singh, Foundations of Medical Imaging. Wiley, 1993. [18] F.-H. Lin, S.-Y. Tsai, R. Otazo, and etc., \Sensitivity-encoded (sense) proton echo- planar spectroscopic imaging (pepsi) in the human brain," Magnetic Resonance in Medicine, vol. 57, no. 2, pp. 249{257, 2007. [19] A. Chen, C. Cunningham, and etc., \High-speed 3T MR spectroscopic imaging of prostate with yback echo-planar encoding," Magnetic Resonance Imaging, vol. 25, pp. 1288{1292, 2007. [20] Y. Gao and S. J. Reeves, \Optimal k-space sampling in MRSI for images with a limited region of support," IEEE Transactions Medical Imaging, vol. 19, no. 12, pp. 1168{1178, 2000. [21] Y. Gao and S. J. Reeves, \Fast k-space sample selection in MRSI with a limited region of support," IEEE Transactions Medical Imaging, vol. 20, no. 9, pp. 868{876, 2001. [22] S. K. Plevritis and A. Macovski, \MRS imaging using anatomically based k-space sam- pling and extrapolation," Magnetic Resonance in Medicine, vol. 34, pp. 686{693, 1995. [23] S. Reeder, J. Brittain, T. Grist, and Y.-F. Yen, \Least-squares chemical shift separation for 13C metabolic imaging," Journal of Magnetic Resonance Imaging, vol. 26, no. 4, pp. 1145{1152, 2007. [24] S. J. Orfanidis, Introduction to Signal Processing. Pearson Education, Inc, 2009. [25] W. Tang, S. Reeves, and D. Twieg, \Fast joint estimation of local magnitude, decay, and frequency from single-shot mri," Proceedings of SPIE, vol. 6498, p. 649818, 2007. [26] C.T.Kelley, Iterative Methods for Optimization. SIAM, 1999. 103