Image reconstruction in dynamic contrast enhanced magnetic resonance imaging

Dynamic contrast enhanced magnetic resonance imaging (DCE-MRI) is a powerful tool to detect cardiac diseases and tumors, and both spatial resolution and temporal resolution are important for disease detection. Sampling less in each time frame and applying sophisticated reconstruction methods to overcome image degradations is a common strategy in the literature. In this thesis, temporal TV constrained reconstruction that was successfully applied to DCE myocardial perfusion imaging by our group was extended to three-dimensional (3D) DCE breast and 3D myocardial perfusion imaging, and the extension includes different forms of constraint terms and various sampling patterns. We also explored some other popular reconstruction algorithms from a theoretical level and showed that they can be included in a unified framework. Current 3D Cartesian DCE breast tumor imaging is limited in spatiotemporal resolution as high temporal resolution is desired to track the contrast enhancement curves, and high spatial resolution is desired to discern tumor morphology. Here temporal TV constrained reconstruction was extended and different forms of temporal TV constraints were compared on 3D Cartesian DCE breast tumor data with simulated undersampling. Kinetic parameters analysis was used to validate the methods.

IMAGE RECONSTRUCTION IN DYNAMIC CONTRAST ENHANCED MAGNETIC RESONANCE IMAGING by Liyong Chen A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Bioengineering The University of Utah May 2012 Copyright ©Liyong Chen 2012 All Rights Reserved Th e Un i v e r s i t y o f Ut a h Gr a d u a t e S c h o o l STATEMENT OF DISSERTATION APPROVAL The dissertation of has been approved by the following supervisory committee members: , Chair Date Approved , Member Date Approved , Member Date Approved , Member Date Approved , Member Date Approved and by , Chair of the Department of and by Charles A. Wight, Dean of The Graduate School. Liyong Chen Edward V.R. DiBella 2/23/2012 Frederic Noo 2/23/2012 Sarang Joshi 2/28/2012 Edward W. Hsu 2/28/2012 Eugene Kholmovski 2/23/2012 Patrick A. Tresco Bioengineering ABSTRACT ABSTRACT Dynamic contrast enhanced magnetic resonance imaging (DCE-MRI) is a powerful tool to detect cardiac diseases and tumors, and both spatial resolution and temporal resolution are important for disease detection. Sampling less in each time frame and applying sophisticated reconstruction methods to overcome image degradations is a common strategy in the literature. In this thesis, temporal TV constrained reconstruction that was successfully applied to DCE myocardial perfusion imaging by our group was extended to three-dimensional (3D) DCE breast and 3D myocardial perfusion imaging, and the extension includes different forms of constraint terms and various sampling patterns. We also explored some other popular reconstruction algorithms from a theoretical level and showed that they can be included in a unified framework. Current 3D Cartesian DCE breast tumor imaging is limited in spatiotemporal resolution as high temporal resolution is desired to track the contrast enhancement curves, and high spatial resolution is desired to discern tumor morphology. Here temporal TV constrained reconstruction was extended and different forms of temporal TV constraints were compared on 3D Cartesian DCE breast tumor data with simulated undersampling. Kinetic parameters analysis was used to validate the methods. iv 2D imaging with serial acquisition of different slices is regularly used for myocardial perfusion imaging. 3D imaging has potential advantages including robustness to through plane motion, and accuracy of sizing ischemia. Here 3D stack-of-stars sampling with spatiotemporal TV constrained reconstruction is developed and is shown to be a promising alternative for myocardial perfusion imaging. Other groups proposed a number of reconstruction algorithms for undersampled MRI recently, including HYPR-LR, PR-FOCUSS, k-t BLAST/k-t SENSE, k-t FOCUSS and regularized iterative SENSE. The work here reveals the relationships among these methods by incorporating these algorithms into a generalized reference image framework. Reconstruction of simulated data, as well as undersampled myocardial cine datasets and perfusion datasets, showed that the superiority of x-t and x-f reference image is sensitive to the data characteristics and baseline images. All of the above efforts will lead to improvements in the diagnosis of diseases like myocardial ischemia and breast tumors, through improving image quality and better quantifying kinetic parameters. This is dedicated to my wife, Haiyan, my wonderful son, Isaac, and my parents. CONTENTS CONTENTS ABSTRACT ....................................................................................................................... iii LIST OF TABLES ............................................................................................................. ix LIST OF ACRONYMS ...................................................................................................... x ACKNOWLEDGEMENTS ............................................................................................. xiii 1. INTRODUCTION ....................................................................................................... 1 1.1 Organization of Thesis ......................................................................................... 4 2. MAGNETIC RESONANCE IMAGING BACKGROUND ....................................... 5 2.1 Nuclear Magnetic Resonance Effect .................................................................... 5 2.1.1 The Underlying Quantum Mechanics Explanation of Nucleus .................... 5 2.1.2 RF Excitation ................................................................................................ 6 2.1.3 Relaxation Mechanism .................................................................................. 7 2.2 Signal Localization ............................................................................................... 7 2.2.1 Slice Excitation ............................................................................................. 8 2.2.2 Spatial Encoding ........................................................................................... 8 2.2.3 Sampling Pattern ........................................................................................... 9 2.2.4 FOV and Resolution ................................................................................... 11 2.3 MRI Scanner Hardware Architecture ................................................................. 12 2.4 Contrast Mechanisms ......................................................................................... 13 2.5 DCE MR Imaging and Cine Imaging ................................................................. 14 2.5.1 DCE Breast Tumor Imaging ....................................................................... 14 2.5.2 DCE Myocardial Perfusion Imaging .......................................................... 15 2.5.3 Cardiac Cine Imaging ................................................................................. 16 3. RECONSTRUCTION BACKGROUND .................................................................. 18 3.1 Non-Cartesian Reconstruction ........................................................................... 18 3.2 Parallel Imaging ................................................................................................. 20 3.3 Constrained Reconstruction ............................................................................... 21 3.3.1 The Generalized Series (GS) Model ........................................................... 21 3.3.2 Compressed Sensing ................................................................................... 22 3.3.3 Total Variation ............................................................................................ 24 vii 3.3.4 The Extensions of Total Variation .............................................................. 26 3.4 Summary ............................................................................................................ 28 4. 3D DCE BREAST TUMOR IMAGING WITH TCR ............................................... 30 4.1 Introduction ........................................................................................................ 30 4.2 Image Reconstruction ......................................................................................... 32 4.3 Data Acquisition and Kinetic Parameter Fitting ................................................ 34 4.3.1 Data Acquisition and Simulation ................................................................ 34 4.3.2 Implementations .......................................................................................... 35 4.3.3 Kinetic Parameters from Breast Data ......................................................... 38 4.4 Experimental Result ........................................................................................... 38 4.4.1 Comparison of Different Methods .............................................................. 38 4.5 Discussion .......................................................................................................... 44 4.6 Conclusion .......................................................................................................... 50 5. 3D STACK-OF-STARS MYOCARDIAL PERFUSION IMAGING ...................... 51 5.1 Introduction ........................................................................................................ 51 5.2 k-Space Acquisition ........................................................................................... 54 5.2.1 3D Stack-of-stars Acquisition ..................................................................... 54 5.2.2 Numerical Simulation ................................................................................. 54 5.2.3 Phantom Study ............................................................................................ 57 5.2.4 Human Study .............................................................................................. 58 5.2.5 Comparison of 3D SOS and 2D Radial ...................................................... 58 5.3 Reconstruction and Analysis .............................................................................. 59 5.3.1 Image Reconstruction ................................................................................. 59 5.3.2 Image Analysis (Perfusion data: SNR/CNR) .............................................. 60 5.4 Experiment Results ............................................................................................ 60 5.4.1 Numerical Simulations ................................................................................ 60 5.4.2 Phantom Studies .......................................................................................... 61 5.4.3 Effect of the Approach to Steady State on the Image Quality .................... 63 5.4.4 3D Stack-of-Stars Images in Human Subjects ............................................ 63 5.4.5 Comparison of 3D SOS and 2D Radial ...................................................... 65 5.5 Discussion .......................................................................................................... 66 5.6 Conclusion .......................................................................................................... 69 6. GENERALIZEDREFERENCEIMAGEFRAMEWORK ......................................... 70 6.1 Introduction ........................................................................................................ 70 6.2 The Extended GS Model .................................................................................... 72 6.2.1 Generalized Reference Framework ............................................................. 72 6.2.2 Relation to PR-FOCUSS, k-t FOCUSS, and k-t BLAST/k-t SENSE ........ 75 6.2.3 Relation to HYPR-LR ................................................................................. 76 6.2.4 Sparsity in x-f and x-t Domains .................................................................. 77 viii 6.2.5 Inclusion of Coil Sensitivity Profiles .......................................................... 79 6.3 Experiment Results ............................................................................................ 80 6.3.1 HYPR-LR Simulation ................................................................................. 80 6.3.2 k-t FOCUSS and x-t FOCUSS Comparison ............................................... 81 6.4 Discussion .......................................................................................................... 89 6.5 Conclusion .......................................................................................................... 91 7. SUMMARYAND CONCLUSIONS ......................................................................... 92 REFERENCES ................................................................................................................. 96 LIST OF TABLES LIST OF TABLES 4.1 The linear relationship between pharmacokinetic parameters determined from the images reconstructed using different methods and those determined from the fully sampled data. The data in the bracket are the correlation coefficients. ............................ 46 6.1 Comparison of different GS model derived methods ................................................. 80 6.2 RMSE of six myocardial perfusion imaging dataset using k-t FOCUSS and x-t FOCUSS with the temporal average image as baseline image ......................................... 90 LIST OF ACRONYMS LIST OF ACRONYMS BOLD Blood oxygen level dependent. A mechanism in fMRI that MRI T2* weighted signal is dependent on blood oxygen level that is believed to be related to neural activity CAIPIRINHA Controlled aliasing in parallel imaging results in higher acceleration. A multiband technique that aliases multiple slices differently by slice dependent phase modulation CNR Contrast to noise ratio DCE Dynamic contrast enhanced DCE-MRI Dynamic contrast enhanced magnetic resonance imaging FLASH MRI Fast low angle shot magnetic resonance imaging GRAPPA Generalized autocalibrating partially parallel acquisitions. A parallel imaging technique that calculates the missing k-space points by linearly weighting neighbor points of multiple coils; the weighting factors are calibrated using extra acquired low frequency k-space HYPR Highly constrained back projection xi HYPR-LR Local HYPR reconstruction k-t BLAST k-t Broad use linear acquisition speedup technique k-t FOCUSS k-t Focal underdetermined system solver k-t SENSE k-t Sensitivity encoding MRI Magnetic resonance imaging Multiband A technique exciting multiple slices simultaneously NUFFT Nonuniform fast Fourier transform POCS Projection onto convex sets PR-FOCUSS Projection reconstruction focal underdetermined system solver RIGR Reduced encoding imaging by generalized-series reconstruction SENSE Sensitivity encoding SNR Signal to noise ratio SPGR Spoiled gradient recalled sequence SOS Stack-of-stars TCR Temporal constrained reconstruction xii TE Echo time TR Repetition time TV Total variation UNFOLD Unaliasing by Fourier encoding the overlaps using the temporal dimension ACKNOWLEDGEMENTS ACKNOWLEDG EMENTS I would like to thank many professors, friends, relatives, and supporters who have made this happen. Firstly, I thank my advisor Prof. Edward V.R. Di Bella. His patience and earnest scientific attitude is a good example for me to do research. Without his generous support for the past several years, I could not have done research continuously. He always patiently listens to others and rarely interrupts others, which is also a great example for me. He organizes the group members very well, and all of us are inspired to do fast imaging by applying sophisticated image reconstruction algorithms and other modeling techniques. I would like to thank Dr. Parker, director of the Utah Center for Advanced Imaging Research (UCAIR), Department of Radiology. UCAIR wins a lot of respect worldwide for its wide range of research in medical imaging. When I came to Utah, I did research on dynamic molecule simulations, which was boring to me. So I wrote an email to Dr. Parker to apply for a position at UCAIR. He forwarded my email to other professors at UCAIR. That is how I came to UCAIR. I would like to thank Dr. Noo who gave me the interview after Dr. Parker forwarded the email, and it lasted for about an hour, much longer than expected. He also taught me 3D image reconstructions, which is the most mathematically intensive course. It reminded me of high school mathematical contests. I had a lot of fun doing the exercises. xiv I would like to thank Dr. David Feinberg, adjunct professor of University of California, Berkeley, also, president of Advanced MRI Technologies. He is a great pulse sequence designer targeting real applications of new techniques. I did three months of internship and several months' part-time work in his group, which allowed to me to learn multiband excitation and simultaneous echo refocusing (SER) techniques. This experience taught me that pulse sequence design is an important avenue for fast MRI. He always told me to learn more physics. I would like to thank my committee members. I took Dr. Hsu's "Principles of MRI" course. He assigned difficult exercises covering broad areas, which made me struggle when doing homework; however, I learned a lot while pursuing the correct answers. Eugene's Matlab program that transfers raw data from dat format to Matlab format is the first step for our reconstruction. I also got a lot of other help from him, especially pulse sequence programming and ICE programming. Dr. Joshi works at Scientific Computing Institute (SCI) and I do not see him very often. Every time I see him, I have the strong impression that he is a very responsible person and takes everything earnestly. I would like to thank all UCAIR people. Professor Larry Zeng is a good helper in everything. He brings us a lot of fun with his special way of teaching and presenting. UCAIR has good seminars and discussions, which is a great learning opportunity. Lastly, but not least, I would like to thank my family. My wife gave me a lot of support, especially when I did the internship. My son's smile and funny behavior always comfort me when I feel stressed. My parents are far away from here and have not seen me for several years, but they never bother me and always let me focus on my studies. I CHAPTER 1 INTRODUCTION 1. INTRODUCTION Magnetic resonance imaging (MRI) is widely used in the field of healthcare due to its ability to detect cancer and accurately diagnose several other diseases noninvasively. Compared to Computed Tomography (CT), it does not have the risk of radioactive harm to people and provides much more flexibility in image contrast by using different imaging sequences and scanning protocols. MRI can be designed to show T1, T2 and proton density contrast in images and to measure other physical parameters, such as velocity, temperature, or diffusion coefficients. Several imaging techniques have been proposed for individual diseases by introducing contrast agents. Among them dynamic contrast enhanced MRI (DCE-MRI) is a well-known MRI technique that monitors the enhancement of a tissue or organ continuously by acquiring a series of MRI images after injecting contrast agent. This helps to show tissue perfusion to identify tumors or to make movies of the heart to obtain anatomical and functional information of cardiac disease so as to diagnose these diseases. Unlike optical imaging in which a whole image is acquired at the same time, MRI collects raw data in frequency domain in a pixel-by-pixel scheme which limits its acquisition speed. Since the invention of MRI over 30 years ago, image acquisition speed and quality have greatly improved as the result of endeavors of investigators worldwide. These advances in image speed and quality have been achieved through a) hardware 2 improvements, b) innovative pulse sequence designs and c) more efficient sampling techniques. Examples of these types of improvements include: a) multiple RF receiver coils, high performance gradient coil design and high static magnetic field. b) echo planar imaging (EPI) (1), fast spin-echo (2) and GRASE (3). c) partial k-Space, radial sampling, propeller sampling and spiral sampling. The efficient sampling methods often require specific reconstruction techniques. Typically fast techniques take advantage of several of the above methods to speed up acquisition and keep good image quality. One example is multiplexed EPI (4), which acquires multiple images in one EPI echo train by interleaving signal from several slices using simultaneous echo refocusing (SER) technique (5) and exciting several slices simultaneously using multiband technique (6). It requires high performance hardware, such as multiple RF receiver coil and high static magnetic field, and often applies partial Fourier methods. All these types of techniques have been applied to speed up DCE-MRI. Since DCE-MRI captures each image in a short acquisition window, it assumes that the image remains unchanged during the readout. This assumption may not hold well, especially for 3D acquisitions, which typically take longer, and in the presence of normal physiological body motions, such as respiratory and cardiac movement. For 3D DCE breast tumor imaging, the image contrast changes especially during the contrast agent uptake and washout. Better temporal resolution may reduce the violation of the above assumption of static image for each time frame. At the same time high spatial resolution is desirable for discerning the tumor morphology. Although there is still controversy as to what spatial and temporal resolution should be and which has a high priority, undersampled Cartesian sampling with sophisticated reconstruction is a good 3 way to better balance the tradeoff between spatial and temporal resolution so as to accurately track the tissue enhancement to tell the difference between malignant and benign tumors. One of the reconstruction methods that have been successfully applied to dynamic MRI is spatiotemporal total variation (TV) constrained reconstruction proposed by our group and applied to some 2D applications, such as myocardial perfusion imaging (7), and temperature imaging (8). Since MRI images are complex-valued, the TV constraint can be of different forms, such as complex form, separate real and imaginary form, and separate magnitude and phase term. One contribution of this thesis is to compare different forms of temporal TV constraints on3D DCE breast tumor datasets and verify the resulting images with pharmacokinetic parameter analysis. For DCE myocardial perfusion imaging, currently multislice 2D imaging can provide only 4 slices with about 3mm in-plane resolution even with parallel imaging technique in clinical setting (9). 2D radial sampling has been proposed to be superior to 2D Cartesian sampling for its robustness to motion and undersampling (7,10). 3D myocardial perfusion imaging has several potential advantages to 2D imaging, such as contiguous coverage of the left ventricle and high SNR (11,12), although it is limited by its longer acquisition window. Thus for 3D myocardial perfusion, undersampled k-space data is especially desirable to shorten the long acquisition window and reduce the effect of heart motion. A 3D form of radial sampling, ‘3D hybrid radial' (also known as ‘stack-of-stars') is applied to 3D myocardial perfusion imaging with spatiotemporal TV constrained reconstruction to verify its feasibility, and this is the second contribution. 4 A third contribution of this thesis is to generalize the generalized series (GS) model and create a framework to include several recent algorithms. Many sophisticated reconstruction algorithms with different names have been proposed by the MRI image reconstruction community. The GS model and compressed sensing are two well-known reconstruction algorithm families from which many algorithms can be derived. In this thesis, some of these algorithms are shown to be derived from the extension of the GS model. This helps to better understand these algorithms and the relationship between the GS model and compressed sensing. 1.1 Organization of Thesis This thesis is organized as follows: Chapter 2 gives a basic overview of the principles of MRI. Chapter 3 gives a background of MRI reconstruction that introduces three important topics in MRI image reconstruction field so as to provide a background to understand much of the research described in the subsequent chapters. Chapters 4-5 apply the spatiotemporal TV constrained reconstruction method to two 3D applications. Chapter 4 compares different forms of TV for 3D DCE breast tumor imaging with undersampled Cartesian SPGR sequence. Chapter 5 applies the complex form of spatiotemporal TV constrained reconstruction to 3D stack-of-stars myocardial perfusion imaging technique with an ECG gated saturated recovery turboFLASH sequence. Chapter 6 investigates several other groups' reconstruction methods mainly from theoretical level and presents a general framework that can include them, such as HYPR-LR, PR-FOCUSS, k-t FOCUSS, and regularized iterative SENSE. Chapter 7 summarizes the main achievements of the thesis. I CHAPTER 2 MAGNETIC RESONANCE IMAGING BACKGROUND 2. MAGNETIC RESONANCE IMAGING BACKGROUND This chapter introduces the principle of magnetic resonance imaging. The basic physical principle of magnetic resonance imaging lies in the nuclear magnetic resonance effect. 2.1 Nuclear Magnetic Resonance Effect In 1946, Purcell and Bloch discovered the nuclear magnetic resonance effect independently (13,14). Both were awarded the Nobel Prize for physics for this discovery. In this section, the NMR effect will be explained by going through quantum mechanics explanation of nucleus, RF excitation and relaxation mechanisms. 2.1.1 The Underlying Quantum Mechanics Explanation of Nucleus Atoms that have an odd number of protons and/or neutrons have angular spin momentum (called a spin), and they act as small magnetic dipoles. These dipoles are randomly aligned, so the net magnetization of an object is zero. When exposed to static magnetic field, the spin will align with the magnetic field and precess at frequency of 0 0  B where  is gyromagnetic ratio. Different elements have different gyromagentic ratio, for H1, r=42.58MHz/T. For simplicity, here we only focus on hydrogen, H1, which is commonly used for MRI due to its abundance in the body. When placed under external 6 magnetic field, some spins are aligned parallel to the magnetic field while others are antiparallel. The antiparallel state has a higher energy than the parallel state, and the spins at the lower energy state can transition to the higher energy state by absorbing energy 0 2 1   E  h where h is Planck's constant (h=6.62x10-34Js). Due to the preference for the lower energy state, based on the Boltzmann distribution, the ratio of the number of protons in the lower energy state to those in the higher energy state is given by e kT n n      , where  n is the number of spins in the lower energy state,  n is the number of spins in the higher energy state, 0 2 1     h Is the energy difference, k is the Boltzmann constant (k=1.38x10-23J/k), and T is the temperature. The total effect is that the object under static magnetic field will have nonzero net magnetizationM n n h  ( ) 4 1 0     . 2.1.2 RF Excitation The spins precess at the frequency of 0 0  B (for H1 at 3T, 0  =123.2MHz), which belongs to the spectrum of radio transmission. Resonance happens if an electromagnetic wave of the same frequency ( 0 0  B ) is applied, which is called "Radio Frequency excitation." In practice, the electromagnetic wave is generated by adjusting the electric current of the RF coils, which is analogous to applying time varying field B1+ in the plane orthogonal to static magnetic field. The governing equation for the RF excitation is given as: M B M   dt d where          z y x M M M M and          0 1 0 1 0 sin( ) cos( ) B B t B t   B . Due to the RF 7 excitation, the net magnetization M0 will tilt to transverse plane. The magnetization precessing in transverse plane sends out a signal that can be received with a receiver coil that has the same frequency. 2.1.3 Relaxation Mechanism The spins tend to recover to equilibrium state after RF excitation. Two different relaxation mechanisms, one is called spin-lattice or longitudinal relaxation, the other is called spin-spin or transverse relaxation, were found to affect the magnetization. To account for these two relaxations, Bloch extended the above equation by adding another relaxation term:                           1 0 2 2 T M M T M T M dt d z y x M B M  , where T1 is the longitudinal relaxation time, and T2 is the transverse relaxation time, as illustrated in Figure 2.1. Effective field inhomogeneity, which originates from static magnetic field inhomogeneity and susceptibility difference, causes the protons to dephase more quickly. In this case, T2* is used, which is the combination effect of T2 relaxation and effective field inhomogeneity. 2.2 Signal Localization The signal can be detected due to the NMR effect. However, it cannot indicate the resonating protons in specific regions, which is important for imaging. 8 Figure 2.1 Illustration of T1 and T2 relaxation. The governing equation for longitudinal relaxation is Mz(t)=M0(1-exp(-t/T1)); for transverse relaxation it is Mxy(t)=Mxy(0)exp(-t/T2). 2.2.1 Slice Excitation An RF pulse at the resonance frequency can be applied to excite the protons. In general, it is possible to excite the whole volume with 0 0  B without applying any gradient, which is named as "nonselective excitation." A certain bandwidth of radio frequency ( z G Z z    0 ( ) ) can be specified to excite a certain portion (slab or slice, for 3D imaging, called slab; for 2D imaging, called slice) in the slice direction by applying a gradient in the slice direction at the same time. Theoretically, a sinc function has a rectangle shaped spectrum which gives the desired rectangular slice profile. In practice, there are time limitations for the RF pulse which create an imperfect rectangle frequency box that results in an imperfect slice profile. 2.2.2 Spatial Encoding Without a spatial encoding gradient, the signal equation can be described as i t xy M ( , t) ( )e 0    r  r where  (r) is the excited object in the position of r and 0 0  B . By applying additional magnetic gradient fields, the magnetic fields varies spatially in x, 9 y, z directions. B t B t G t x G t y G t z x y z ( , ) ( , ) ( ) ( ) ( ) 0 ' 0 r  r    where            z y x r is the position with respect to the isocenter of the magnet, x G , y G and z G are the gradients of magnetic field in x, y, and z directions, respectively. The signal equation in the position of r becomes ( ( ) ) 0 0 ( , ) ( )     t i t d xy M t e     G τ r r r . The term 0  can be demodulated, so the signal detected can be simplified as     r r G τ r S t e d t i d 0 ( ) ( ) ( )    . This can be further described as    r r k rS t e d i2 (t ) ( ) ( )   where  t t d 2 0 1 ( )  ( )   k G τ . For a given time course of gradient fields applied after RF excitation, a series of sampled data, which is known as "k-space," can be generated based on the above equation. The excited object turns out to be an inverse Fourier transformation of k-space. Typical pulse sequences are shown in Figure 2.2. 2.2.3 Sampling Pattern In the conventional case, the k-space is acquired line-by-line, referred to as a Cartesian sampling scheme, which is the most popular pattern. It can be reconstructed simply by performing an inverse discrete Fourier transform. Non-Cartesian sampling patterns, like radial (15) or spiral patterned (16), have also been proposed. Radial acquisition is robust to motion and undersampling. Spiral sampling efficiently uses the gradient, and samples very fast. However, non-Cartesian sampling requires more complicated reconstruction algorithms, such as gridding, which will be discussed more in Chapter 3. Some 2D and 3D sampling patterns are shown in Figure 2.3. 10 Figure 2.2 Typical pulse sequence diagrams. (a) is the 2D Spoiled Gradient Recalled (SPGR) sequence. (b) is 3D SPGR sequence. Both sequences are composed of slice excitation, spatial slice/phase encoding using variable gradient amplitudes (hatched pulses on Gy axis) and readout, spoiling part. The slice excitation is to excite a portion of object in slice direction. The spatial slice/phase encoding is to encode object in slice/phase encoding direction so as to recovery the object. The difference between 2D SPGR and 3D SPGR is that 3D SPGR has slice phase encoding (hatched pulse on Gs axis) while 2D SPGR has only a slice refocusing gradient. The readout gradient is to encode the excited object in readout direction. The spoiling part is to get rid of in plane magnetization by dephasing, and this is specific for SPGR sequence. 11 Figure 2.3 Some common sampling patterns. Top row, left to right: Cartesian 2D, radial, spiral. Bottom row, left to right, Cartesian 3D, stack-of-star (or 3D hybrid radial), stack-of-spiral. 2.2.4 FOV and Resolution In k-space acquisition, the continuous Fourier transform of an object is sampled at discrete points. For simplification, only Cartesian sampling is considered here. The discrete sampling can be thought as multiplying k-space with a comb function with interval width k ; which means the convolution of excited object with the inverse Fourier transform of a comb function, which is another comb function with reciprocal interval width k 1 . This discrete sampling brings periodic duplicated object with adjacent distance of k 1 , which is described as field of view (FOV). The number of samples is denoted as base resolution n, the k-space ranges from k n   2 (denoted as max k ) to k n 1) 2 ( , and the spatial resolution is given by nk 1 . 12 2.3 MRI Scanner Hardware Architecture Figure 2.4 illustrates the system architecture of MRI scanner. The scanner tunnel contains built-in RF coil, gradient coil, and magnet, which are the basic components of a scanner. Custom RF coil and gradient coil can be used. A pulse sequence that runs on the host computer will control the operation of the switching of gradient coil and RF coil. The signal detected will be recorded and reconstructed into images. Figure 2.4 The system architecture of MRI scanner. The scanner is composed of magnet, gradient coil, and RF coil. The gradient coil and RF coil operation is controlled by pulse sequence. The detected RF signal can be reconstructed and shown on the host computer monitor. 13 2.4 Contrast Mechanisms Although many physical factors, such as velocity, diffusion coefficient, and temperature, play a role in image contrast and signal intensity, most MRI images can be categorized into three types, proton density weighted images, T1weighted images and T2(*) weighted images, for which the image contrast are dominated by three parameters: proton density (PD), T1 relaxation and T2 relaxation, respectively. Many diseases cause changes to at least one of these three parameters, which makes MRI very useful to diagnose disease. PD weighted images are acquired with long TR and short TE. The regions with more protons will have high magnitude, while regions with fewer protons will have low signal. In practice, a short TR with a very small flip angle can be used to acquire PD images more rapidly. T1 weighted images are acquired with short TR and short TE, and shorter T1 has larger signal. For dynamic contrast enhanced T1weighted MR imaging, the contrast agent is injected into a vein and the gadolinium contrast decreases the T1 value of its local environment. These T1 changes can be tracked by T1 weighted MRI signal intensity changes, so as to track the contrast agent concentration changes. T2(*) weighted images are acquired with long TR and long TE, and objects with longer T2(*) have larger signal. For functional MRI, the blood oxygen level dependent (BOLD) effect contributes to the signal changes, which can be used to track the neuron activity of functional area. 14 2.5 DCE MR Imaging and Cine Imaging Dynamic MRI, such as cine imaging and DCE applications in oncology, angiography, and perfusion of the heart and other organs, is an important and rapidly growing area in medical imaging. The dynamic MRI applications considered in the thesis include DCE breast tumor imaging, DCE myocardial perfusion imaging, and cine cardiac imaging. The strengths and limitations of these three applications will be discussed below. 2.5.1 DCE Breast Tumor Imaging X-ray mammography is the current standard method for the detection of breast tumors. It performs well in postmenopausal women and less well in perimenopausal women (17). It is not very sensitive for many cases and also exposure to X-ray is hazardous. DCE-MRI is an important routinely used MRI technique for detecting breast tumors. DCE breast tumor imaging is capable of acquiring contrast uptake patterns, which are used to distinguish malignant and benign tumors. DCE-MRI has been reported to have sensitivity (the fraction of patients with disease who test abnormal) approaching 100% and no radiation exposure is involved (18).The main limitation of DCE-MRI in the investigation of breast lesions lies in its low specificity(the fraction of patients without disease who test normal) (19). It was reported that multivariate models combining tumor morphology and contrast uptake dynamics have a superior diagnostic accuracy than that based on tumor architecture or contrast uptake pattern alone (20). This requires both high spatial resolution and temporal resolution; although there is still controversy on how much to prioritize spatial versus temporal resolution (19). The proposed image 15 reconstruction methods in this thesis with more extensive undersampling are expected to obtain higher spatial-temporal resolution, and to do this without SNR reduction. This should lead to better diagnostic accuracy of breast tumors. 2.5.2 DCE Myocardial Perfusion Imaging Myocardial perfusion imaging is important in evaluation of patients with coronary disease by providing functional and prognostic information. SPECT is a very widely used test to evaluate the myocardial perfusion. PET and stress echocardiography are also performed in clinical myocardial perfusion practices. However, all of the above methods have limitations. For SPECT, there are always tradeoffs between and specificity. PET offers better image quality than SPECT and it provides high sensitivity and specificity (21) but PET is still not widely available for cardiac perfusion imaging due to costly scanner and cyclotron operation and expensive radionuclide (22) and it also lacks the spatial resolution obtainable with MRI. Stress echocardiography with contrast agents can to some degree reflect myocardial perfusion but it requires adequate skill of the operator. MRI has the potential to become a widely used tool for myocardial perfusion measurement. Compared to SPECT and PET, it is more realizable for MRI to get high spatial-resolution, temporal resolution and volumetric coverage. The spatial resolution of MRI makes it possible to differentiate between subendocardial and subepicardial regions (23) which is not possible with clinical SPECT and PET. Subendocardial perfusion defects can be a more sensitive indicator of ischemia (23). The study of perfusion and MPR (myocardial perfusion reserve, the ratio of stress to rest perfusion) distribution which is a research focus requires high temporal resolution to get signal intensity-curves especially under stress condition (24). 16 Current DCE imaging methods with MRI cannot provide full spatial coverage of the heart while at the same time provide images with high spatial and temporal resolution and the necessary SNR. Only a few 2D slices can be acquired per heartbeat, especially in stress condition when the heart rate is high. It was reported that four slices/beat with ~3mm in-plane resolution was possible in a general clinical application with multicoil methods (9,25). Approximately 10 short axis slices (6mm thick) and 1-3 long axis slices are desired to give full spatial coverage of the left ventricle. The proposed image reconstruction methods in this thesis are expected to improve spatial coverage and spatial-temporal resolution without compromised SNR and significant artifact, for myocardial perfusion imaging. This development could lead to improved diagnostic accuracy of coronary artery disease. As well, accurate sizing of ischemic regions could improve predictions of how the patient will do in the future and enable optimal treatment selection. 2.5.3 Cardiac Cine Imaging Cardiac cine MRI imaging is a basic technique to assess the contractile cardiac function. FLASH and SSFP sequences are typically used for cardiac cine imaging, and SSFP is reported to be superior to FLASH in terms of SNR and CNR in both 1.5 and 3.T although it contains some artifact (26). For this technique, one or several slices are imaged at each stage or "phase" of the cardiac cycle, and the images acquired at different stages can be viewed as a movie, so termed as "cine." Due to the short acquisition window of each stage, typically a portion of k-space lines of each image are acquired in each heartbeat and the lines from multiple heartbeats can be combined as full k-space to recover the image at each stage. The number of the k-space lines acquired in each 17 heartbeat is termed as lines per segment. Given the spatial resolution, the temporal resolution is proportional to the lines per segments, and the acquisition time is inverse proportional to the lines per segments. Furthermore, the image quality is dependent on heart rate regularity and motion consistency. To gain better spatiotemporal resolution while keep good image quality, several algorithms have been applied to cardiac cine imaging with undersampling dataset gaining an acceleration factor of about 4 to 6 without much image degradation (27,28). Current cardiac cine imaging can provide one or two slice with both high spatial and temporal resolution in a reasonable breath-hold time. The main limitation is that multiple breath-holds are needed to acquire stacks of cardiac slices which result in long acquisition time and inaccurate cardiac volume due to inconsistent respiratory motion (29). The k-space undersampling combined with sophisticated reconstruction techniques make is possible to acquire more slices so as to mitigate or overcome this limitation. In this thesis, the SSFP cardiac cine imaging datasets are used to test the superiority of some algorithms that can be derived from the extension of the Generalized Series (GS) model in Chapter 6. I CHAPTER 3 RECONSTRUCTION BACKGROUND 3. RECONSTRUCTION BACKGROUND As stated in Chapter 2, the MRI data acquired by scanner give values in k-space. The frequency domain k-space data need to be transformed to get an image, and this process is termed as "reconstruction." In this chapter, an overview of reconstruction algorithms will be presented. Three topics will be covered: non-Cartesian reconstruction, parallel imaging, and constrained reconstruction. These three topics are only enough to cover the main aspects of reconstruction techniques, but will suffice to provide a background to understand much of the research described in the subsequent chapters. For many applications, other specific reconstruction procedures are required to get good images, such as off resonance correction and motion correction. These topics will not be covered here. 3.1 Non-Cartesian Reconstruction Radial sampling, spiral sampling, and other more arbitrary sampling patterns have been proposed in literature and have gained great popularity due to robustness to motion, undersampling and efficiency, although Cartesian sampling is the most widely used in clinical practice. In this thesis, the 3D form of radial sampling, 3D stack-of-stars sampling pattern is applied to myocardial perfusion imaging. It has several potential advantages to 2D imaging, such as contiguous coverage of left ventricle, through-plane 19 motion and high SNR (11,12). Noteworthy, the robustness to through-plane motion is due to the fact that the slab thickness of 3D imaging is much larger than the slice thickness of 2D imaging. Thus motion out of plane for 3D imaging will be more negligible than 2D imaging. For conventional Cartesian sampling, the reconstruction can be easily and efficiently implemented by simple inverse Fourier transform. However, for non-Cartesian sampling, inverse Fourier transform is no longer applicable. There are several options for non- Cartesian reconstruction, such as projection reconstruction (30), conjugate phase reconstruction (31,32) and resampling (33-35). Projection reconstruction, which does filtered back projection of 1D inverse Fourier transform of each line, can be applied only to radial sampling. Conjugate phase reconstruction, which calculates an integral for each pixel separately, is extremely computationally expensive. One feasible and efficient solution is to sample the non-Cartesian data to Cartesian data, then do inverse Fourier transform. One of the most commonly used resampling methods is called gridding (also regridding) (35). There are many variations of gridding. The mathematical description of one typical of gridding algorithm is given here: ˆ ( , ) [( ( , ) ( , ) ( , )) ( , )] ( , ) y y x x x y x y x y x y x y k k k k M k k M k k S k k w k k C k k      where ( , ) x y w k k is density compensation function, and ( , ) x y C k k is convolution kernel, ( , ) x y S k k is the measured data sampling, ( , ) y y x x k k k k    is the Cartesian grid sampling. One simple and efficient resampling algorithm used by our group is to sample data from non- Cartesian point to the nearest integer point using triangle based interpolation (7). In the undersampled case, the gridded data cannot be inverse Fourier transformed to get the 20 final image due to aliasing artifact. Instead it is incorporated into a constrained reconstruction framework as the fidelity term 2 2 ~ ~ d m WF ) m (    where m ~ is the image estimate, F is Fourier transform, W is a sampling mask (W is a diagonal matrix, and it is the identity matrix if there is no undersampling), which will be covered in Chapter 4. 3.2 Parallel Imaging Parallel imaging is a method that acquires the data from multiple receiver coils that have different spatial sensitivities in order to increase the speed of MRI acquisition. For the past few decades, many different parallel imaging reconstruction techniques have been proposed. They can be categorized into two types, image based reconstruction, such as PILS (36) or SENSE (37) and k-space based reconstruction, such as SMASH (38) or GRAPPA (39). Various algorithms have been extended from both types, such as image based TSENSE (40), kSPA(41), and PARS (42); k-space based TGRAPPA (43), iGRAPPA(44), and SPIRiT(45). In this chapter, only the most basic and widely used methods, SENSE and SMASH/GRAPPA will be explained. SENSE reconstruction represents the signal of each pixel of each coil image as    R l k l kl I S 1  , where k is the coil index, kl S is the sensitivity profile of the kth coil at location l, l  is the signal value at location l , l ranges from 1 to R and specifies the pixel location and its aliased pixel location, and R is the acceleration factor. This can be written as matrix formI  Sρ, and can be solved asρ S S S I H 1 H ( )  . SMASH is based on the assumption that the missing k-space steps y y im k e  can be modeled as linear combination of coil sensitivities, and it represents the k-space data of 21 composite image as     Nc k k y m y y k comp S k m k n S k 1 ( ) ( ) where k is the coil index, m is the skipped k-space lines, c N is the coil number. m k n is calculated by fitting coil sensitivity profile to y y im k e  . GRAPPA extends SMASH, by representing the k-space of each coil, rather than the composite, as     Nc k k y m l y y k S k m k n S k 1 ( ) ( ) . m k n is fitted by acquiring extra autocalibration lines. Image from multiple coils are reconstructed separately, and combined using sum of square method. 3.3 Constrained Reconstruction Constrained reconstruction was proposed several decades ago in MRI reconstruction. There are so many different kinds of constraints, including implicit and explicit, in the literature that a thorough discussion of constrained reconstruction is out of the scope of this work. Here we discuss several seminal and review papers which help to sketch the roadmap. 3.3.1 The Generalized Series (GS) Model Early constrained reconstruction work has been reviewed by Liang Z-P (46), and the constraint was defined as a priori information, bounds, or parametric models. Partial Fourier reconstructions that incorporate phase information, extrapolation algorithm based on the assumption of finite image support were reviewed there. In addition, several parametric models, including autoregressive moving average model, localized polynomial approximation, and the generalized series (GS) model, were reviewed there. 22 Here only the generalized series model is explained. The GS model is a general mathematical framework to handle prior constraints, and image was represented as  l gs l l  (x) a (θ, x)where l  is parameterized basis function and l a is the series coefficients for which the number is much smaller than image pixel number. Tsao et al. (47) extended generalized series model as     Nterms l static dynamic l l R R c 1  (x) (x) (x)  (x)where  (x) is the reconstructed image, (x) static R and (x) dynamic R are static and dynamic reference images, l c are the unknown basis coefficients, terms N is the number of basis coefficients, and (x) l  is the basis function. This model is reported to be able to incorporate at least 14 algorithms. 3.3.2 Compressed Sensing Compressed sensing is hot topic in the signal processing area and it is a technique that recover signal from underdetermined linear systems by minimizing L1 norm of the sparse signal and/or its transformation (48). MRI reconstruction is one of many applications that compressed sensing has gained much popularity in recent years. Compressed sensing is a great improvement over classic sampling requirements enforced by Shannon sampling theorem. Shannon sampling theorem states that the sufficient condition to recover a band limited function G(f) with band limit of B is to sample data at a rate higher than 2B 1 , which is illustrated in Figure 3.1. In the scenario of MRI, MRI images can be inverse Fourier transformed from k-Space data as noted in Chapter 2. To avoid overlapping (aliasing), the k-space interval k should be less than FOV 1 . Here k 23 Figure 3.1 Illustration of Shannon sampling theorem. (a) is the band limited signal or function G(f) with band limit of B. (b) is the signal recovered with sampling rate of 2B 1 . (c) is the signal recovered with sampling rate higher than 2B 1 . (d) is the signal recovered with sampling rate lower than 2B 1 with aliasing showing up. 24 is the sampling rate and FOV is the field-of-view of the image which corresponds to 2B in Figure 3.1. Using compressed sensing, the MRI images can be recovered from measurements that are drastically fewer than those required by Shannon sampling theory by constraining the L1 norm of images and/or transformed images. Several different forms of optimization schemes and many algorithms have been proposed to solve this problem, and a list of software can be found in (49). Here we introduce the constrained optimization scheme that can be formulated as min( ) 1 1 m s.t.   2 F m d u where is the image transformation, m is the image estimate, d is the k-space measurement, u F is the undersampled Fourier transform and is the noise level.  can be identity matrix if MRI images are sparse , and many different kinds of transformations have been proposed to enforce sparsity, such as wavelet transform (50), finite difference operation (the constraint term in this case is "total variation") (7,51) and curvelet transform (52). 3.3.3 Total Variation One of the most used constraint terms in this thesis is total variation (TV), which can be used as the sparsity for compressed sensing. Generally the total variation of a real valued or complex valued function f , defined on an interval a,b , is       1 0 1 ( ) sup ( ) ( ) np i i i P a b TV f f x f x , where the supremum runs over the set of all partitions { ,..., } 0 np p  x x , p is a partition of a,b (53). For a function of n dimensional real variables defined on  which is an open subset of n  , the TV norm of the function is ( ) sup{ : ( ; ), 1, } 1            TV f  f wdx w C w x n c where ( ) 1  c C is the set of 25 smooth functions on  that vanish on the boundary of  and f (x) can be real or complex(54). Correspondingly, w(x) can be real or complex functions. This definition is applicable to nondifferentiable f (x) . This definition helps to formulate the primal-dual algorithm for TV minimization(55). To make the total variation more intuitive and understandable, here the total variation of the differential function of one variable is considered specifically, and the definition gives   b a b a TV ( f ) f (x) dx ' .The discrete form of total variation is       1 0 ( ) ( 1) ( ) n n TV f f n f n with its complex form as          1 0 ( ) ( 1) ( ) ( ( 1) ( )) n n R R I I TV f f n f n i f n f n , where f (n) R and f (n) I are the real part and imaginary part of function, respectively. In this thesis, gradient descent (or time marching) method is used to minimize the total variation(56). To avoid singularities in the derivative of the functional as shown by Acar et al (57),  is a small positive constant is added     2 ( ) dx df TV f . The gradient descent method gives     2 2 2 ( ( )) dx df dx d f dx d TV f with its discrete implementation, denoted here as S , as              2 2 ( 1) ( ) ( 1) ( ) ( 1) ( ) ( 1) ( ) f n f n f n f n f n f n f n f n S . The gradient update term can be written as n n n n f  f  S 1 , where n  is the step size of nth iteration. Other 26 algorithms have been proposed to solved TV minimization problems, such as fixed points (58) and primal-dual method (55). 3.3.4 The Extensions of Total Variation There are several extensions of TV minimization problems. One extension is to extend L1 norm to homotopic L0 norm. In the scenario of dynamic MR imaging, the total variation is the sum of the absolute value (L1 norm) of signal intensity difference of adjacent image pixels in spatial or temporal dimension. The homotopic L0 norm of the constraints become ( ) ( ( 1) ( ) ) 1 0       x x H f h f x f x , here h is the homotopic L0 norm function, such as  | | ( ) 1 x h x e    ,   | | | | ( ) x x h x (59). Interestingly, anisotropic diffusion and robust statistics that have been reported to be closely related to each other (60) are related to compressed sensing (including L1 norm and homotopic L0 norm of transformed image with finite difference operation) (59). Investigating total variation from anisotropic diffusion perspective, the flux function (which is influence function in robust statistic) is    2 | | ( ) ' ( ) x x g x h x .The incorporation of into TV term is making it analogous to Huber function, L2 norm in the low difference value area, L1 norm in the high difference value area; which is desirable for denoising in noisy areas while keeping the sharp edges. In the scenario of DCE-MRI applications, for signal intensity change curve of each pixel, the use of the total variation will keep its sharp enhancement while smooth it when the signal intensity change is not sharp (61) which is illustrated in Figure 3.2. This is a valuable property so that the data do not get smoothed in sharp transitions as with L2 27 Figure 3.2 The illustration of signal intensity change curve of one typical pixel. The dashed line with square is the simulated tissue enhancement curve. The stair-wise area of the thick line indicates the area that needs to keep sharp edges. The smooth area of the thick line indicates the area that needs to be smooth. norm. The same principle applies to the total variation in spatial domain, and this has been widely used for images and natural scenes (51,56). There are other good models such as wavelets (50) but the extensive discussions are beyond the scope of this work. Another extension is to iteratively refine the TV constraint term, such as Bregman iteration algorithm (62,63) and reweighted L1 algorithm (64). Here we introduce only the Bregman iteration algorithm. The Bregman iteration algorithm is reported to better preserve fine structures than the TV regularization method by adaptively refining the regularization term, and the minimization problem becomes: argmin( ( , )) 1 1) 2 2      f Ef d D f f k k f k  where ( , ) k1) D f f is the Bregman distance between f and k1 f , defined 28 by k TV k TV k k TV D f f f f f f f 1) 1 1 1 ( , ) ,          , where , is inner product operation, k TV f 1  is the subgradient of the TV norm at point k1 f , and Ef  d is the encoding equation. A simplified implementation is as below: 0 0 argmin( ) 0 1 2 1 2                 k v v d Ef f d v Ef f k k k k TV f k . The TV minimization in each Bregman iteration can be solved with various algorithms as indicated before, such as gradient descent (or time marching)(56), fixed points (58) and primal-dual method (55). 3.4 Summary Non-Cartesian reconstruction, parallel imaging and constrained reconstruction are independent of each other. Many algorithms have been proposed to incorporate several of them to gain additive benefits. For example, k-t SENSE (65) incorporates both parallel imaging and constrained reconstruction; when k-t SENSE is applied to non-Cartesian sampling, it incorporates all of the above three topics. Radial GRAPPA (66,67) incorporates non-Cartesian sampling and parallel imaging. In the following chapters, different constrained reconstruction algorithms with Cartesian/non-Cartesian sampling will be investigated and applied to clinical data. In Chapter 4, different types of temporal TV constraints are included into the POCS framework, and implemented for 3D Cartesian DCE breast tumor application. In Chapter 5, a 3D stack-of-stars pattern was applied to myocardial perfusion imaging. Here only complex TV was tried based on the results in Chapter 4 that complex TV performs best. 29 In Chapter 6, the generalized series model constrained reconstruction methods were extended to include several recent algorithms from the literature, and were applied to cardiac cine imaging and DCE myocardial perfusion imaging. CHAPTER 4 3D DCEBREAST TUMOR IMAGING WITH TCR 4. 3D DCE BREAST TUMOR IMAGING WITH TCR In this chapter, it is demonstrated that application of TCR to DCE breast tumor imaging may help to achieve better image quality. Different forms of temporal constraints are presented and projection onto convex set (POCS) framework is introduced to include these constraints into reconstruction. DCE breast tumor data are tested using these algorithms, and the resulting images are analyzed with kinetic parameter models for verification. The results are published in Magnetic Resonance Imaging, Vol. 28, Page 637-645, 2010 (68), and reproduced here with permission. 4.1 Introduction Dynamic MRI plays an important role in a number of clinical MRI applications, such as in dynamic contrast enhanced MRI and functional MRI. For such applications, a common strategy used to balance tradeoffs between spatial resolution and temporal resolution is to reduce sampling of k-space data at each time frame. A variety of reduced k-space data acquisition and reconstruction techniques have been proposed to do this. Examples include sliding window, UNFOLD (69), keyhole (70), RIGR (71), k-t BLAST/k-t SENSE (27), k-t FOCUSS (28), compressed sensing (72,73), and HYPR (74). Most of these methods use constraints (also known as "prior information") to compensate 31 for the information loss from reduced sampling. Sophisticated methods are typically needed to reconstruct the images with the constraints. For dynamic MRI applications, the images of adjacent time frames are often assumed to be similar, especially when motion is minimal, in which case temporal TV is a reasonable regularization term (7,75). In this paper, two tools were applied to the implementation of constrained reconstruction. One powerful tool is the Projection onto Convex Sets (POCS) formalism, which can include prior information flexibly and has been extensively used in MRI reconstruction applications (46,76-78). Another tool is the gradient descent method, which is regularly used for the minimization of an objective function, and can be considered as a type of projection to be included in the POCS framework. For constrained MRI image reconstruction, the regularization term is typically used in its complex form (27,28,79). However, separate real and imaginary TV regularization, and separate magnitude and phase regularization terms, have also been investigated by several investigators. Fessler et al. (80) reported that the L2 norm of the spatial derivative in separate magnitude and phase form worked better than that in complex form on simulated phantom data. He et al. (81) reported that separate real and imaginary constraints produced results similar to the use of the complex form of regularization on phantom data. In this paper, different forms of temporal TV terms are compared for reconstruction of undersampled DCE-MRI data acquired in breast cancer patients. The paper is organized as follows: in the first section, we present the form of the fidelity term and of various temporal TV terms, including complex TV, real and imaginary TV, and magnitude TV. Next, specific implementation details for the serial 32 and parallel POCS methods are presented. Finally, we present and discuss results of undersampled breast DCE-MRI reconstructed using the different POCS methods. 4.2 Image Reconstruction Two types of POCS, serial POCS (also known as sequential POCS) and parallel POCS, as summarized in (78), are used in this paper. For serial POCS, different projections are sequentially applied to update the data term at each iteration, while for parallel POCS different projections are weighted to update the data term each iteration. Both the L2 norm of the fidelity term and the temporal TV term are convex functions. The gradient descent forms of both can be viewed as projections. The fidelity term is defined as the L2 norm:    Nf t (m) W t Fm t d t 1 2 2  ~ ( ) ~(r, ) ( ) , [4.1] where m~(r, t) is the image estimate including all time frames, t is from 1 to time frame number f N , F is the (2D) spatial Fourier transform applied on each time frame in the dynamic sequence, W(t) is a binary undersampling pattern of time frame t (W is diagonal matrix, and it is identity if no undersampling), that changes each time frame and matches the undersampling pattern of the acquired k-space, and d(t) is the undersampled data in k-space of one image slice of time frame t (see acquisition section). The projection corresponding to the fidelity constraint term can be denoted as n n n m  m   ~ ~ 1 with the updating term n  given by         Nf t H H Nf t H W t F W t Fm t d t F W t Fm t F W t d t dm d 1 1 ~ ( ( ) ) ( ( ) ~(r, ) ( )) ( ( ) ~(r, ) ( ) ( ))  [4.2] 33 where H is the Hermitian transpose operator and H F is the 2D inverse Fourier transform applied on each time frame. Besides the mandatory data fidelity convex set described above, any other convex sets can be used to regularize the solution. Three forms of the temporal TV constraint are considered here. The first and most widely used form is the complex temporal TV term: ) ε dt dm ( dt dm ψ(m)    2 1 ~ ~ ~ [4.3] where ε is a small positive constant to avoid singularities in the derivative of the functional as shown by Acar et al (57). The gradient descent method gives     2 2 2 ~ / ~ ~ dt dm dt d m dm d with its discrete implementation, denoted here as S, as              2 1 1 2 1 1 ~ ~ ~ ~ ~ ~ ~ ~ t t t t t t t t m m m m m m m m S . The projection corresponding to the complex temporal TV term can be written as n n n n m  m  S  ~ ~ 1 , where n  can be set to a constant step size. The second form of TV constraint is to use separate real and imaginary temporal TV terms, defined as: ) ε dt dR m ( dt dR m (m)    2 1 ~ ( ~) ( ~)  where ) ~ (m R is the real part ofm ~ [4.4] ) ε dt dI m ( dt dI m (m)    2 1 ~ ( ~) ( ~)  where ) ~ (m I is the imaginary part of m ~ [4.5] 34 The gradient descent projection gives ' ' 1 ~ ~ n n n n m  m  S  and '' '' 1 ~ ~ n n n n m  m  i S  where n n n dm d m S ~ ( ~ ) '   n n n dm d m S ~ ( ~ ) ''   and ' n  and '' n  can be Polyak's step size (82) 2 2 ' ' ( ~ ) n n n S  m   and 2 2 '' '' ( ~ ) n n n S  m   [4.6] that make it unnecessary to fit the step size or can be useful to provide a reference for step size selection. The third form of temporal TV constraint is to use separate magnitude and phase terms. It was found that use of temporal magnitude and phase TV terms gave only slightly better reconstructions than temporal magnitude temporal TV alone (see the Discussion section), and so magnitude alone was used in this work: ) ε dt dM m ( dt dM m (m)    2 1 ~ ( ~) ( ~)  [4.7] where ) ~ (m M is the magnitude part ofm ~ . The gradient descent projection gives ''' ''' 1 ( ~ ) ( ~ ) n n n n M m  M m  S  where n n n dm d m S ~ ( ~ ) '''   and ''' n  can be Polyak's step size [19] 2 2 ''' ''' ( ~ ) n n n S  m   (4.8). 4.3 Data Acquisition and Kinetic Parameter Fitting 4.3.1 Data Acquisition and Simulation Breast DCE-MRI data were acquired using a 3D spoiled gradient echo pulse sequence with the following imaging parameters: TR=2.35-3.16 msec, TE=0.99-1.24 msec, flip angle=10-15º using a seven channel dedicated breast coil. Temporal resolution per frame was 12-15 seconds with data acquired with 6/8 reduced Fourier space in the 35 phase and slice directions and elliptical acquisition in the kx-ky plane. The acquisition matrix for the breast data varied between 256 x (80-104) x 80 of 42-60 time frames. The acquisition was bilateral, with the read direction left to right. The fast inverse Fourier transform (IFT) was performed in the read (kx) direction, and the ky-kz datasets were extracted from each slice in the x dimension. Four datasets from three study participants with histopathologically confirmed breast cancer were obtained under an IRB approved protocol. One subject was imaged on two separate occasions. Undersampled k-space data were simulated by deleting a portion of the acquired phase encodes in the ky and kz directions. In the outer areas in the ky direction, one in every two points was sampled; while in the kz direction, one in every three points was sampled. An example of the sampling pattern over a series of time frames is shown on Figure 4.1(a). In the k-space center, a 6x6 window of ky-kz phase encodes were fully sampled for every time frame, as shown in Figure 4.1(b). The net acceleration was R=6 (17% of the k-space data were used). The elliptical partial Fourier acquisition of the original data further increased the undersampling to an acceleration factor of R=10, though note that the "true" data used to compare with the constrained reconstruction methods had only six times as many samples as the undersampled version, so our results are reported as using an acceleration of R=6. 4.3.2 Implementations Four reconstruction methods were implemented: parallel POCS with complex temporal TV term, serial POCS with complex temporal TV term, serial POCS with separate real/imaginary temporal TV term, and serial POCS with magnitude TV alone. These were denoted as parallel+complex, serial+complex, serial+real/imaginary, and 36 Figure 4.1 Undersampling pattern for breast DCE data. The gray circles were not sampled, and the black circles were sampled. (a) The outer k-space sampling of eight adjacent time frames. (b) A typical k-space sampling pattern of one time frame. serial+magnitude, respectively. Reconstruction using a simple "sliding window" method was presented for comparison and for algorithm initialization. Sliding window was implemented by inverse Fourier transformation of k-space data after filling missing measurements in k-space using the corresponding measurement from the most recent time frame in which it was acquired. This is not technically a "sliding window," but this method gave better results than interpolating based on all of the data within a window. 37 Parallel POCS requires weighting factors to be chosen, and serial POCS accomplishes similar weighting by the number of iterations each convex set is performed before going on to the next convex set. For parallel POCS, 150 iterations were used and the fidelity weighting was set to be 1, and temporal TV weighting was set by trying a range of different parameters from 0.01 to 1.2 on the data set of each coil, and the weighting factor of 0.1 was selected based on minimizing the root mean squared error (RMSE) in a test dataset. This weighting was then used for all of the coils and all of the datasets. The RMSE value was calculated by square root of the mean square difference between the reconstructed images and the true images that were reconstructed by inverse Fourier transform of the elliptical partial k-space data. For serial POCS, the weighting of different terms are affected by both the step size of the projection term and the iteration number ratios among different convex set projections. For simplicity, the iteration ratio of 1:1:1 or 1:1 was applied for all serial POCS methods. For the fidelity term, the step size was set to be 1, and could be viewed as replacing the measured k-space data in the corresponding k-space points of the current estimate. For the temporal TV term, Polyak's step size (Equation [4.6]) was used for the initial estimation of the constant step size. There were two reasons for not using Polyak's step size to adapt the step size at each iteration: one was that from our tests, the constant step size converges faster; the other is that computation of Polyak's step size takes some time during each iteration. The POCS methods were applied independently to sparse data obtained from each of the seven coils. The reconstructions from each coil were then combined using the square root of the sum of squares. 38 4.3.3 Kinetic Parameters from Breast Data After dynamic images were reconstructed and the baseline precontrast signal subtracted, the signal intensity difference curve of every pixel was fitted to the extended Tofts-Kety two compartment model for tissue contrast agent (CA) concentration: C (t) K C (t) e v C (t) p p k t p trans t ep     where Ktrans and kep are the transfer constant and rate constant respectively,  is the convolution operator, vp is the blood plasma volume fraction, and Cp(t) is the concentration of CA in the blood plasma (83). The linearized regression method described in (84) was used to perform curve fitting, and a population averaged arterial input function (AIF) was used for Cp(t) (85). To quantify the linear relationship between the kinetic parameters generated from constrained methods and that generated from the true data, L1 regression, where the sum of absolute difference is minimized, was used. This type of analysis was used due to its robustness to outliers, rather than least square regression that minimizes the sum of squares difference. 4.4 Experimental Result 4.4.1 Comparison of Different Methods The images reconstructed from undersampled data of one subject using the undersampling pattern described in Figure 4.1 are shown in Figure 4.2. Figure 4.2(a) shows a time frame in a typical DCE sequence obtained from full k-space data using IFT. Figure 4.2(b) shows the corresponding time frame reconstructed using IFT on the undersampled data. Figure 4.2(c-f) shows the corresponding time frame reconstructed 39 Figure 4.2 Comparison of reconstructions from full data and R=6 (using pattern shown in Figure 4.1) data using different methods with all coils. (a) The 22nd time frame reconstructed from full k-space data using IFT. The corresponding time frame reconstructed from undersampled data using IFT method is shown in panel (b), the parallel+complex in (c), serial+complex in (d), serial+magnitude in (e), and the sliding window method in (f). using the sliding window (SW) method, serial+complex, parallel+complex, and serial+magnitude, respectively. The RMSE plots of images reconstructed from this subject with the simulated undersampling using the different POCS methods are shown in Figure 4.3. The separate real and imaginary TV does not work as well as complex TV and magnitude TV in terms of RMSE. Figure 4.3 shows that for all time frames, serial+complex, parallel+complex, and serial+magnitude constrained reconstructions had reduced RMSE as compared to the SW method. The serial+magnitude method consistently had the lowest RMSE. 40 Figure 4.3 RMSE plot for each time frame computed for different methods with one data set of all coils. The black line with dots is parallel+complex, the blue line is sliding window method, the green dash line with plus is serial+complex, the red line with circles is serial+magnitude, the cyan line with circles is serial+real/imaginary, and the blue dash line with circle is serial+magnitude/phase. Figure 4.4 compares the mean signal intensity time curves from one breast lesion region using different methods. Figure 4.4(a) shows the region of interest (ROI) in the breast. Figure 4.4(b-c) compares the mean signal intensity curves for the region. The time curves obtained from POCS methods matched with the full data reconstructions closely. POCS methods were applied independently to the sparse R=6 data obtained from each of the seven coils. The reconstructions from each coil were then combined using the square root of the sum of squares and the results are shown in Figure 4.5. Figure 4.5(a) shows the images reconstructed from full k-space data using IFT. Figure 4.5(b-e) shows the images reconstructed from parallel+complex, serial+complex, serial+magnitude, SW, 41 Figure 4.4 Comparison of dynamics of reconstructions from undersampled data (R=6) in two different breast lesion regions using different methods. (a) Images showing the one ROI in the breast lesion, indicated by the small black rectangle. (b) Comparison of mean signal intensity time curves for the lesion region shown in (a), and (c) is the magnified images of (b). The magnified image shows the signal intensity curve of SW methods have a larger deviation from that of the true images. The red line is the full sampled reference, the blue line is parallel+complex, the cyan line is serial+complex, the black isserial+magnitude, and the dash green is sliding window method 42 Figure 4.4 continued. 43 Figure 4.5 Reconstructions results from all coils. a-e (left column): the 12th time frame of reconstructed images (from top to bottom, the left column is full sampled image, parallel+complex, serial+complex, serial+magnitude, sliding window). f-i (right column): the difference image between the corresponding left image and (a). Larger residual errors of the images reconstructed with the serial+magnitude and SW methods are evident in the bottom two rows of the right column. The left column images are scaled to [0,30]; and the right column images are scaled to [0,2]. 44 respectively. Figure 4.5(g-i) shows the difference images with fully sampled image (a) of corresponding time frame. The residual error of image reconstructed using the complex constraint was smaller than that of the other methods. The relationship of kinetic parameters from the reconstruction of the R=6 data and that from fully sampled data are shown in Figure 4.6 for parallel POCS with complex temporal TV, the reconstruction method that correlated with full data best in the tumor area. Pharmacokinetic parameters determined from the images reconstructed using this method showed strong linear correlation with those determined from the fully sampled data. The results are summarized in Table 4.1, with the exception of p v . The values of p v determined from all methods are close to zero ( p v =0.027±0.051 for fully sampled data, p v =0.014±0.033forparallel+complex). 4.5 Discussion Three types of temporal TV terms were used for reconstruction of undersampled breast data. A POCS method with gradient descent method was used for implementation. From signal intensity curves (Figure 4.4), difference images (Figure 4.5) and kinetic parameters (Figure 4.6), it can be seen that parallel+Complex, serial+Complex, and serial+magnitude are capable of accurately reproducing the measured signal intensity curves and pharmacokinetic parameters. In Fessler's work (80), it was demonstrated that the L2 norm of the spatial derivative of separate magnitude and phase performed better than that of the complex form in one set of simulated data. It is known that most MRI data have relatively smooth phase in spatial dimension, and that the spatial phase L2 is a strong constraint. In the breast DCE 45 Figure 4.6 The correlation plots between kinetic parameters (Ktrans, kep) generated from images using parallel+complex and that using IFT of fully sampled k-space data, with Ktrans plot shown in (a), kep plot shown in (b). The kinetic parameters data sets came from all four datasets' lesion areas. 46 Table 4.1 The linear relationship between pharmacokinetic parameters determined from the images reconstructed using different methods and those determined from the fully sampled data. The data in the bracket are the correlation coefficients. Ktrans kep parallel+complex Y=0.97X+0.00 (0.98) Y=0.95X+0.00 (0.85) serial+complex Y=0.97X+0.00 (0.98) Y=0.96X+0.00 (0.84) serial+magnitude Y=0.93X+0.00 (0.98) Y=0.90X+0.00 (0.92) sliding window Y=0.97X+0.00 (0.97) Y=0.94X+0.00 (0.80) MRI data used here, serial+magnitude/phase worked approximately the same as serial+magnitude without the temporal phase TV term. The phase TV term did not help significantly to get better images (see Figure 4.3), possibly due to the good initialization of the phase images. Keyhole techniques have been used for quantitative dynamic contrast enhanced breast MRI (86) and it was reported that the minimum keyhole size should be restricted by the approximate minimum size of the expected lesions. Parallel imaging and generalized series model have been combined to accelerate dynamic contrast enhanced breast cancer imaging with an acceleration factor of 3-4 for the 2D case (87). The techniques used here enabled an acceleration factor of 6 for 3D acquisitions while maintaining good correlation with the kinetic parameters in the tumor. Polyak's step size is useful to adjust the step size range, which is necessary for L1 norm minimization. However, it was found that the optimal step size, in terms of efficiently reaching images that gave minimum RMS errors compared to the full data reconstruction, was not close to Polyak's step size. However, all of the POCS methods performed robustly to perturbations (0.5α) in the optimal step size of α (75), which 47 implies relatively few trial step sizes are needed to get the optimal step size. As well, the same step size performed well for different datasets, indicating that it is likely not a parameter that needs to be found for each dataset. The parallel POCS methods used here can also be termed temporally constrained reconstruction (TCR) (75). TCR was performed by iteratively minimizing a cost function of a data fidelity term and constraint terms. The cost function was defined asC  argmin{ (m~)  (m~)}. The gradient descent method gives n n n n n m  m C m  m  Wd WFm S  ~ ~ ' ( ~) ~ ( ~ ) 1 . Interestingly, this gradient descent implementation of TCR is the same as parallel POCS with complex temporal TV constraints used here. A spatial TV term can also be added to the parallel (or serial) POCS methods (7). Compared to parallel POCS, serial POCS makes it easier to add other constraints. In (78), several convex sets and associated projection operators pertinent to MRI data reconstruction were defined, such as fixed phase and limited object support. Other types of convex functionals can also be included in the POCS framework, such as a prior image constraint (88). For the kinetic parameter analysis, the signal intensity curves were not converted to the contrast agent concentration, as is sometimes done (89). Most of the signal intensity difference curves were expected to remain linear with the contrast agent concentration, and since the truth was computed in the same manner, it was not essential to perform this extra step. The extra step of conversion to concentration would have made the comparison less direct. 48 The sliding window method also worked reasonably well. This can be explained in part through Figure 4.7: although the sliding window method is often biased more than the temporally constrained method when the intensity increases sharply, the fitted line is not so sensitive to this area because of sharp onset of the population AIF. Thus similar kinetic parameters similar to truth can be found even when the curves from the undersampled reconstructions have a slower onset. Figure 4.7 The delta signal intensity values and model fits for images reconstructed from parallel+complex, sliding window and fully sampled data (denoted as "TCR," "SW," "True," respectively). The delta signal intensity value is the intensity value with the mean value of the first 10 time frames subtracted off. The plot demonstrates that in particular for the sliding window curve, the fit is much closer to the fully sampled data when the tissue curves increases sharply. This is because the sharp onset of the AIF used does not permit an exact fit to the curve. 49 Another limitation of this study is that the signal intensity curves of two of the data sets show relatively slow and steady uptake in the lesions, and this is particularly amenable to undersampling. Also, vp can be difficult to analyze with voxel wise curves in general, and did not correlate well here (r=0.55 for parallel POCS with complex TV). The results here were based on simulated undersampling, in order to have a measure of truth. The simulated undersampling may not have as high temporal resolution as actual undersampled acquisitions. Actual undersampled data will likely be more robust since temporal and/or spatial resolution can be increased, and the effective rate of change of the contrast will be slower and easier to reconstruct. It is also possible that the current acquisitions were undersampled temporally and that the time curve will vary more rapidly when temporal resolution is improved. In this case, there may be greater differences between the reconstruction methods. Spatial and temporal resolution are crucial for MRI breast cancer detection and characterization (90). The proposed method can be used to increase temporal resolution without compromising spatial resolution and SNR loss. High spatial resolution is required for detection of small lesions and for assessment of lesion morphology. Thus, this approach may increase the detectability of small lesions. It is also possible that the high temporal resolution can make it possible to track the tissue enhancement curve more accurately and thus increase specificity for diagnosing malignancy (91). The computation time is demanding, especially when the dataset size is large. In a Matlab (The Mathworks, Natick, MA, USA) implementation on a desktop PC, it takes approximately 40s to reconstruct one slice. Considering that 20-40 slices of 5-12 coils will have to be reconstructed in a clinically acceptable time span of 30-60 seconds, the 50 computation time will have to be improved by a factor of between 60 and 640. An efficient C++ implementation on a more powerful computer will provide improvement in computation time. Recently published papers have shown that computationally intensive medical imaging tasks can be processed on a graphics processing unit to increase computation speed by a factor of 85-100 (92,93). Taking advantage of these techniques, clinical implementation would be feasible. 4.6 Conclusion We have demonstrated that temporal TV could be successfully employed for dynamic MRI breast perfusion applications. Complex TV or magnitude TV constraints could be used to give good results at an acceleration factor R of 6, which can translate into improved spatial and temporal resolution for DCE breast scans without a cost to image quality. In the tumor area, the best method, parallel POCS with complex temporal TV, gave kinetic parameter trans,R6 K =0.97 trans,R1 K +0.00 with correlation coefficient r=0.98, ep,R6 k =0.95 ep,R1 k +0.00 (r=0.85). These promising methods warrant further study to determine how increasing spatial or temporal resolution affects clinical assessment and management of breast cancer and other cancers. CHAPTER 5 3D STACK-OF-STARS MYOCARDIAL PERFUSION IMAGING 5. 3D STACK-OF- STARS MYOCARDIAL PERFUSION IMAGING In this chapter, the 3D stack-of-stars sequence with spatiotemporal TV constrained reconstruction is demonstrated to be feasible for 3D myocardial perfusion MRI. The stack-of-stars sampling pattern is presented and simulations are undertaken to select the optimal saturation recovery time (SRT) and flip angle value. Then, the acquired 3D myocardial perfusion data are reconstructed and compared with 2D myocardial perfusion data. The results have been composed in a paper entitled "Myocardial perfusion MRI with an undersampled 3D stack-of-stars sequence," which is being prepared for submission. 5.1 Introduction MR myocardial perfusion imaging is an effective method to evaluate perfusion defects and detect cardiac ischemia. Current methods typically provide three to four 2D slices per heartbeat at stress with parallel imaging (94-96). An echo planar readout can provide more than 10 slice spatial coverage with in-plane spatial resolution as high as 1.5mm (97). However, echo-planar is sensitive to chemical shift and susceptibility effects, which thus far have prevented its use in clinical practice. Large spatial coverage 52 of the heart with high spatial and temporal resolution and good SNR is important to improve the utility of cardiac MRI perfusion. Greater spatial coverage makes it less likely to miss ischemic areas and allows for better sizing of ischemia. High spatial resolution can reduce the dark rim artifact (98-100) which can mimic subendocardial defects (23). High temporal resolution can also be important to reduce dark rim effects and to accurately track signal intensity changes. Besides parallel imaging techniques, undersampling with sophisticated reconstructions have been proposed to obtain more spatial coverage and higher spatial and temporal resolution for 2D perfusion scans. k-t SENSE methods using Cartesian undersampling have been reported to give good results for three to four slices with a net acceleration factor of 3 to 4 by acquiring 23-33 phase encoding lines (101,102). Compressed sensing combined with parallel imaging was reported to gain an acceleration factor of 8 by acquiring 16-24 lines (103). Radial undersampling patterns have been explored due to their robustness to motion and undersampling. A constrained reconstruction method with temporal and spatial total variation constraints was reported to acquire 10 slices at rest using 24 rays per slice and five slices were acquired after each saturation pulse and gave image quality comparable to 68 phase encodes with Cartesian data (7) . SW-CG-HYPR was proposed using 16 rays per slice and six to eight slices were acquired per beat (10,104). Some of these accelerated methods are sensitive to motion, or focus on high spatial resolution and do not achieve high coverage. 3D perfusion MRI might be advantageous compared to 2D in terms of volume coverage and a consistent contrast for all slices. 3D also may be more robust to inter-frame motion and may permit greater undersampling, although the longer 3D readout 53 could be sensitive to cardiac motion. Undersampled 3D Cartesian myocardial perfusion imaging with SENSE reconstruction was reported to provide whole left ventricle (LV) coverage with 10 slices and relatively poor spatial resolution of 3x(4.3-4.5)x10 mm3with an acceleration factor of 6, acquiring 110-115 phase encoding lines (11). This 3D method was shown to perform better than 2D multislice imaging in terms of the accuracy of estimating the size of perfusion defects in a phantom (11). However, the limited spatial resolution may make it hard to detect subendocardial ischemia as well as making the acquisition more prone to dark rim artifact. Recently, an undersampled 3D acquisition was reconstructed with the k-t PCA method and was reported to obtain an acceleration factor of 7, acquiring 125 phase encoding lines and providing 10 slices with spatial resolution of 2.3x2.3x10 mm3 (12) and matrix size of 168x133x10 with partial Fourier and elliptical sampling in ky-kz plane. A similar approach but using k-t SENSE was reported to give an acceleration factor of 6.3 and shown to be useful for detection of ischemia in patients (105). In this paper, a 3D sampling pattern with radial sampling in the kx-ky plane and Cartesian encoding in the kz direction is used. This sampling pattern has been termed 3D hybrid radial sampling or 3D stack-of-stars (3D-SOS) sampling. Due to the relatively long acquisition, 3D imaging can have more signal variations for different readouts than 2D, which may result in image artifacts. In this paper, simulations were performed to show the dependence of the signal transients on flip angle and saturation recovery time. Phantom studies were used to analyze the effect of flip angle on image quality. Human studies were performed to further assess the approach. 54 5.2 k-Space Acquisition 5.2.1 3D Stack-of-stars Acquisition An ECG-triggered, 3D turboFLASH sequence with SOS k-space sampling and saturation recovery preparation as shown in Figure 5.1(a-c) was used. Figure 5.1(a) shows an example of the sampling pattern. The 3D-SOS pattern was chosen instead of 3D radial to obtain a cylindrical field of view that better matched the heart. For stack-of-stars myocardial perfusion imaging, inconsistent projections can cause severe streaking artifacts in-plane (7) and, as with 3D Cartesian imaging, there can be crosstalk artifacts in the slice direction. To reduce the effect of inconsistent projections in-plane, the k-space data were acquired by sampling all in-plane radial lines of one partition (one kz encode) with an interleaved pattern, then sampling other partitions as illustrated in Figure 5.1(b). Centric ordering was applied in the slice (kz) direction. The radial sampling was rotated in the temporal dimension and the slice encoding direction so that data sharing can make an evenly distributed fully-sampled 3D-SOS sampling. Such a change in sampling pattern over time is essential for the reconstruction method to be effective. An example of the sampling pattern over a series of time frames is shown in Figure 5.1(c). 5.2.2 Numerical Simulation In order to minimize the signal inconsistencies for 3D-SOS imaging, simulation studies were done to determine the optimal acquisition parameters. The signal of the n-th readout of saturation recovery turboFLASH is (106) a a M n M e a M e n T TR T n SRT xy           1 1 ( ) (1 ) (1 ) 1 1 1 1 [5.1] 55 Figure 5.1 Illustration of pulse sequence. (a) Illustration of stack-of-stars. (b) Schematic diagram for the 3D stack-of-stars acquisition with ECG gating. The centric reordering in slice direction is applied. (c) The sampling pattern of adjacent 3 time frames that interleaving both slice encoding and temporal direction with an interleave factor of 3, and the partial Fourier sampling is on in slice direction. The black circle is the sampled, and the gray circle is not sampled. 56 where 1 * sin( ) 2 , cos( ) 0 T TR T TE M M e a e       , and SRT is the saturation recovery time between the saturation pulse and the first readout radio-frequency (RF) pulse. From equation [5.1], it can be demonstrated that, when                  ) 1 cos ( 1 1 1 1 1 T SRT T SRT T TR T TR e e e  e , [5.2] Mxy(n) is independent of n and the transverse magnetization Mxy immediately reaches its steady-state value. This is an important insight reported in (107) - the readout of a saturation recovery prepared signal can be obtained immediately at steady-state, if SRT and TR and T1 are known and α is selected by equation [5.2]. This expression has been given for the case of 2D spiral-based sequences (107)and similar work has been done in another context - to use saturation pulses to bring spoiled gradient echo sequences to steady-state more rapidly(108). Since T1 is not known a priori, we evaluated the effect of varying T1 on the flip angle given by equation [5.2]. Simulations with physiologically relevant parameters were used to study the effect of non-steady-state readouts in more detail. TR was fixed to 2.5msec to keep the acquisition time short. For each set of T1(ranging from 100 ms to 2000 ms with interval steps of 100 ms), SRT(from 50 ms to 300 ms with interval steps of 1 ms), and flip angle (from 2º to 30º with interval steps of 0.1º), a signal intensity-readout index curve was determined by equation [5.1]. The coefficient of variation (CV), the standard deviation divided by the mean value, was then calculated to evaluate how much the signal varied over the readouts. CV is a measure of the consistency of the signal intensity relative to the readout index, so a perfectly steady-state set of readouts would have CV=0. 57 5.2.3 Phantom Study For comparison to the simulation results, a phantom was imaged on a 3T Trio scanner with a 3D saturation recovery turboFLASH sequence with slice encoding turned off. The DC term (the sum of signal intensity over the excited volume) of an 8cm slab covering the center portion of the cylinder phantom was recorded for 160 readouts after the saturation pulse with saturation recovery time (SRT)=150ms, TR=2.5ms,TE=1.39 ms, flip angle α=8º, 10º, 12º, 14º, and 25º, FOV=220x220mm2. The 160 readouts were composed of 8 sets of 20 readouts which are each composed of 4 sets of interleaved rays (flip angles [0 36º 72º 108º 144º], [18º 54º 90º 126º 162º], [9º 45º 81º 117º 153º], [27º 63º 99º 135º 171º]). To analyze the effect of the transient approach to steady-state on the image quality, the same phantom was also imaged with slice encoding turned on. Image acquisition parameters were: SRT=150ms, TR=2.5ms, TE=1.39 ms, flip angle α=10º and25º, FOV=220x220mm2, number of rays per slice=20 in an interleaved fashion with an interleave factor of 5, 8 slices with partial Fourier factor in slice direction=6/8, slice oversampling factor=25%, spatial resolution=1.7x1.7x10mm3, total readout time≈300ms for one time frame. Imaging was performed twice to evaluate the random spoiling effect that has been reported to show better spoiling in 2D radial imaging (109). The first time was with random RF spoiling. The second was with the standard built-in RF spoiling that uses a phase increment between RF pulses of 50⁰. Since there were no concentration changes between time frames in the phantom, a sliding window reconstruction method was used to reconstruct the images to compare the effect of flip angles and different spoiling patterns. 58 5.2.4 Human Study To determine the feasibility of stack-of-star sampling in-vivo, experiments were performed using a 3T Trio or Verio Siemens scanner under an IRB-approved protocol with an ECG-gated, SOS saturation recovery turbo-FLASH sequence and a 12-element coil array in three subjects. A dose of 0.015-0.05 mmol/kg of contrast agent (Gd-BOPTA or gadofovesettrisodium), was injected at a rate of 5 ml/s followed by a 25 ml saline flush at the same rate. Based on the simulation results, SRT was set to 140-160ms, and flip angle was specified to be 10-14º to compensate for the B1+ inhomogeneity to obtain an actual flip angle of ~8-12º (110). Other image acquisition parameters were as follows: TR=2.1-2.9ms, TE=1.1-1.4ms, FOV=(260-360)x(260-360)mm2, number of rays per slice=20-24 in an interleaved fashion, 8-10 slices with partial Fourier factor in slice direction=6/8, spatial resolution=(1.8-2.8)x(1.8-2.8)x(6-10)mm3, total readout time≈300ms for one time frame. 5.2.5 Comparison of 3D SOS and 2D Radial To compare the SNR, both 3D SOS and 2D radial imaging were performed on a cylindrical phantom. The following parameters were used: SRT=140ms, TR=2.6ms, TE=1.43ms, flip angle=14º, FOV 220x220mm2, the number of projections 20, interleave factor=5, slice thickness=10mm. For 3D, the slice number was 8, and 25% oversampling was performed. For 2D, one slice was acquired. A sliding window reconstruction method was used. A 2D multislice myocardial perfusion imaging dataset with radial sampling was also acquired in one subject with the same dose of 0.015mmol/kg of gadofovesettrisodium as with the 3D-SOS imaging for comparison. Image acquisition parameters for the 2D 59 sequence were: SRT=20ms, TR=2.3ms, TE=1.4ms, flip angle=14º, FOV=360x360mm2, matrix size=144x144, slice thickness=10mm, the number of projections=30, 10 slices were acquired in one heartbeat with five slices after each saturation pulse (the SRTs are 20, 89, 158, 227, 296 ms for each of the five slices in a set). The slices with SRT=158ms were used for the SNR/CNR comparison. 5.3 Reconstruction and Analysis 5.3.1 Image Reconstruction After acquiring the 3D data, the images were reconstructed using spatiotemporal total variation (TV) constrained reconstruction (7,68), with the cost function: ( ) ( ) ( ) ( ) ( ) ( ) 1 2 2 C m W t F Gm t d t TV m TV m temporal spatial N t z f       [5.3] where m(t) represents complex image estimate of time frame t, t ranges from 1 to the total time frame number Nf, G is a nonuniform FFT applied to all slices (34) that transforms images from the x-y-z domain to the kx-ky-z domain, Fz is a Fourier transform in the slice encoding direction that transforms data from the kx-ky-z domain to the kx-ky-kz domain, W(t) is the undersampled binary pattern of time frame t as shown in Figure 5.1(b),d(t) is the measured k-space data of time frame t and α, β are the weighting factor of the temporal and spatial TV constraint term. The gradient descent method was used to minimize the cost function. Different weighting factors for TV constraints were tried on one dataset, and α=0.7 and β=0~0.2 were empirically determined after setting the k-space center (the mean image value) to be ~102. These weights were used to reconstruct other datasets based on the assumption that the reconstruction method was robust to small changes of the weights (75). The image was initialized with an inverse nonuniform FFT 60 of the undersampled radial data, which is similar to doing filtered backprojection of the undersampled projections. The number of iterations was empirically chosen to be 50 because the reconstructed images changed little after 50 iterations based on visual assessment. The reconstruction was applied independently to the data obtained from each coil and the reconstructions from each coil were then combined using the square root of the sum of squares. 5.3.2 Image Analysis (Perfusion data: SNR/CNR) For the phantom experiments, SNRs were calculated by the ratio of the mean value of a 3x3 block from the center area of signal in the images with the standard deviation of the signal intensities from a background area. For the in-vivo experiments, the reconstructed images were evaluated using SNR and CNR. SNR was calculated by the ratio of the mean and standard deviation of the signal intensities from a uniform region in the myocardium of a postcontrast time frame. CNR was computed by (Myopost-Myopre)/σ, where Myopost is the mean of the signal intensities from a uniform region in the myocardium in a postcontrast time frame, Myopre and σ are the mean and standard deviation of the signal intensities from a similar region in the myocardium in a precontrast timeframe. 5.4 Experiment Results 5.4.1 Numerical Simulations Figure 5.2 shows the flip angle-SRT plot calculated using equation [5.2] for three different T1 values. This plot shows that for a given saturation recovery time, the flip angle depends only weakly on T1. 61 Figure 5.2 The flip angle-SRT plot calculated using equation [5.2] for three different T1 values. The flip angle that gives steady state readouts is relatively insensitive to T1 changes. Figure 5.3 shows the results of CV values obtained from simulations with different SRT and flip angle values using T1=700ms. The sets composed of SRT and flip angle, such as (110ms, 12º), (150ms, 10º), and (220ms, 8º), provide the smallest CV values - meaning those readouts were closest to steady-state. Similar results were found for T1=300ms and 1200ms (not shown here). 5.4.2 Phantom Studies Figure 5.4 shows the measured signal intensity (DC term) plotted against the readout number. The signal intensity-readout curves are obtained with SRT=150ms, TR=2.5ms using a 3D saturation recovery turboFLASH sequence with slice encoding turned off and the flip angles α specified as 8º, 12º, and 25º. The periodic fluctuations are consistent 62 Figure 5.3 The coefficients of variation with different SRT and different flip angle when T1=700ms and TR=2.5ms. The sets composed of SRT and flip angle, such as (110ms, 12º), (150ms, 10º), and (220ms, 8º), provide the smallest CV values. Figure 5.4 The signal intensity changes with readout index acquired with TR=2.5, SRT=150ms and the specified α of 8º, 12º, and 25º are shown by the solid lines, and the signal intensity is calculated by doing linear interpolation to obtain the k-space center for each ray. The dashed lines are calculated from equation [5.1] and manually fitted to the solid line. Lower flip angles: 6º, 9º, 18º were used to give better fits. These flip angles are closer to the actual flip angles (112). 63 with in-plane radial angle (period=5) described in the methods where sets of five rays over 180 degrees are repeated with different angular offsets. The periodic signal fluctuations that are consistent with the flip angle changes are due to the gradient delay effect (111). The effect of the signal fluctuations is negligible as described in the discussion section. 5.4.3 Effect of the Approach to Steady State on the Image Quality Figure 5.5 shows the comparison of phantom images acquired with SRT=150ms and flip angles of 8º, 14º and 25º with random (109)and the standard built-in RF spoiling with a phase increment between RF pulses of 50º. The five center slices are shown here. The images acquired with flip angle of 25º show more crosstalk and smearing artifact as indicated by the red and blue arrows, respectively. Images with random RF spoiling have less smearing artifact than that with the standard 50º increment RF spoiling. 5.4.4 3D Stack-of-Stars Images in Human Subjects Figure 5.6 shows three time frames of 3D-SOS images acquired from one typical subject, at precontrast, RV enhancement and LV enhancement phases after reconstruction with spatiotemporal TV constraints. The different slices show the similar contrast and the edge slices show some crosstalk artifact. 64 Figure 5.5 A set of phantom image acquired with SRT=150ms and different flip angles with different RF spoiling pattern. The reconstruction is without gradient delay correction. The top 3 rows are images with Siemens built-in RF spoiling with flip angle of 8º, 14º,25º (from top to bottom) and the bottom 3 rows are images with random RF spoiling with flip angle of 8º, 14º and 25º. The arrows indicate the artifact, including crosstalk artifact (red arrows) and smearing artifact (blue arrows). The center five slices are shown here. The nonuniform images are due to coil effect that the coil closer to the phantom gives more signals. 65 Figure 5.6 One set of 8 slices (left to right) and three time frames at precontrast, RV enhancement and LV enhancement phases, of the representative 3D myocardial perfusion images from another subject, each in a different row. A total of 8 slice encodings were acquired. Partial Fourier factor=6/8 in slice direction was used so 10 slices were acquired. The two edge slices with the most aliasing artifacts were not used. 5.4.5 Comparison of 3D SOS and 2D Radial The SNR of cylindrical phantom using 3D-SOS and 2D radial imaging of the same slice are 64.7±1.42 and 46.8±1.8, respectively. Figure 5.7 shows the comparison of myocardial perfusion images using 3D-SOS and 2D multislice imaging with contrast agent injection using spatiotemporal TV constrained reconstruction. In this case, 3D-SOS provides an SNR of 21.5±3.0 and a CNR of 7.7±1.0 compared with an SNR of 19.8±2.5 and a CNR of 7.0±0.8 for the slice with SRT=158ms of 2D multislice imaging. 66 Figure 5.7 Image comparison of the myocardial perfusion imaging using 3D stack-of- stars (left) and multislice 2D imaging (right) reconstructed with spatiotemporal TV constraints. Both of the images have high SNR. 5.5 Discussion This paper demonstrated the feasibility of 3D myocardial perfusion imaging using 3D-SOS sampling reconstructed with spatiotemporal TV constrained reconstruction to achieve large coverage with high spatial resolution. Simulation and phantom studies were performed to show that the magnetization transient is a function of flip angle and saturation recovery time, and incorrect selection of flip angle and poor spoiling may degrade images. The use of a small flip angle and random spoiling is helpful to reduce image artifacts. Compared to 2D multislice myocardial perfusion imaging, 3D myocardial perfusion imaging requires a longer temporal acquisition window. However, it provides volume excitation which is more robust to through-plane motion and offers contiguous volume coverage, which is reported to be advantageous for sizing perfusion defects (11). The 3D readout is also advantageous because a single, relatively long saturation recovery time can be used for high SNR. For 2D imaging it is not practical to have a long saturation recovery time unless multiple slices are acquired after a single saturation pulse, in which 67 case the saturation recovery time and image contrast is variable. This issue may be manageable (9) but having the same contrast for all of the slices may be an advantage of 3D imaging. SNR measured in-vivo is only a relative indicator in this work, because spatiotemporal TV reconstruction may change the noise characteristics of the images. The phantom studies use fully sampled data with a noniterative reconstruction method and thus reflect the standard SNR measurement. . Compared to a 3D Cartesian acquisition, 3D-SOS inherits the robustness to undersampling and motion of 2D radial acquisition. However, 3D-SOS is more restricted in terms of requiring in-plane isotropic resolution with evenly distributed undersampled projections. For myocardial perfusion imaging, in-plane spatial resolution (maximum kx and ky) is desired to be similar while the resolution in the slice direction (kz) is much coarser, which makes it reasonable to apply a 3D-SOS sampling pattern. The dependence of signal intensity on readout number is determined by flip angle, T1, SRT, and TR. For a saturation recovery spoiled gradient echo pulse sequence with any given set of these parameters, there exists a null point in flip angle where steady-state magnetization is reached immediately (at the first readout). Thus, the degradation of the point spread function (PSF) that arises from readouts that are not at steady state vanishes at the null point, providing the potential for substantially improved image quality. While this is a larger effect with radial imaging due to the repeated sampling of the k-space center, the different weighting of phase encodes in Cartesian readouts also degrades PSF (113,114). Spatial variation in T1 and flip angle makes it impossible to image at the exact null point for all voxels, but the sensitivity to T1 is weak near the null point (Figure 5.2), making it possible to obtain nearly optimal consistency across readouts by flip angle 68 optimization. The existence of this optimal flip angle was tested using measured data, assuming that the T1 and spatial flip angle variation can be ignored (Figure 5.4). Several artifacts arise in the phantom study. The images of 8º and 25º (Figure 5.5) show more crosstalk than that of 14º. The larger flip angle shows more smearing artifact. The greater crosstalk in the slice direction is due to the greater signal variation in the approach to steady state as simulated in Figure 5.3. The smearing artifact seen in Figure 5.5 may result from imperfect spoiling that has more effect on large flip angle images. This is supported by the experiment that random RF spoiling helps to attenuate the artifact(115) as is also shown in Figure 5.5. Gradient delays may bring streaking artifacts to radial sampling, which can be compensated through calibration (Figure 5.8), although this artifact is not obvious in the myocardial datasets. For 3D myocardial perfusion imaging, the slice encoding number is small due to the short acquisition window, resulting in crosstalk or Fourier leakage (116). Also, the slab-excitation profile (especially when a fast RF pulse is applied with a small time-bandwidth Figure 5.8 One slice of phantom image with and without k-space center offset, and the difference image of them. The offset of about .25 is measured using method described in (117), and corrected by adjusting this offset in k-space sampling of NUFFT reconstruction. 69 product) is not perfect which will also degrade edge slices. From our results when reconstructing 6 or 8 kz encodes that were offset (partial Fourier), into 8 or 10 slices, the outermost slice at each edge of the slab showed significant aliasing. Discarding two slices at each edge left approximately six central slices that appeared to be free of aliasing. It is also possible to shorten the acquisition time by reducing the number of readout lines in the higher slice encoding planes. This could enable more oversampling in the slice direction for no net cost in acquisition time. The reconstruction time is demanding, especially when the dataset size is large. In a Matlab (The Mathworks, Natick, MA) implementation on a desktop PC, it takes approximately 10 minutes to reconstruct one slab of 50-60 time frames for one coil. Recently published papers have shown that computationally intensive medical imaging tasks can be processed on a graphics processing unit (GPU) to increase computation speed by a factor of 85 to 100 (92). Taking the advantage of these techniques, clinical implementation could be feasible. 5.6 Conclusion The contributions of this paper include showing the dependence of the transients on flip angle and saturation recovery time, and analyzing the effect of the flip angle on image qualities for 3D SOS perfusion imaging. The initial evaluations show that 3D stack-of-stars myocardial perfusion imaging with spatiotemporal TV constrained reconstruction is a promising alternative to provide images with consistent contrast and contiguous volume coverage of the heart. CHAPTER6 GENERALIZED REFERENCE IMAGE FRAMEWORK 6. GENERALIZEDREFERENCEIMAGEFRAMEWORK In this chapter, it is demonstrated that a generalized series framework proposed decades ago can be extended to include several recent reconstruction algorithms, like HYPR-LR, PR-FOCUSS, k-t FOCUSS and regularized iterative SENSE. First, the mathematical derivations of the generalized series model are given. Then, the relationships of GS model with different algorithms are clarified. Finally, different algorithms are tested on cine imaging dataset and myocardial perfusion dataset. The results are composed in a paper entitled "A generalized framework for reference image reconstruction methods including HYPR-LR, PR-FOCUSS, and k-t FOCUSS," and are accepted and in press by the Journal of Magnetic Resonance Imaging, and reproduced here with permission. 6.1 Introduction A Generalized Series model (GS) was proposed to improve image quality by including prior or reference images for dynamic MRI applications (118). Various related methods have been proposed since then, and it has been shown that many of them, such as the keyhole method (70) and the Reduced encoding Imaging by Generalized-series 71 Reconstruction (RIGR) (119), can be considered as special cases of a formulation similar to the GS method (47). Recently, Highly Constrained Back Projection (HYPR) and related methods such as HYPR-LR and I-HYPR were proposed for dynamic reconstruction of undersampled radial k-space magnetic resonance imaging (74,120-122). Note that two types of HYPR have been proposed. The first type of HYPR is similar to the first iteration of the maximum likelihood method (122) and weights the composite image by the backprojection of the ratio of the 1D inverse Fourier transform of each ray of k-space data, and the projections of the reference image: ) P(r(x)) F (d) ρ r(x) B(   1 1 ~ 1    d p N where p N is the number of projections, r(x)  is the reference (or composite) image, P and B are projection and backprojection operators, 1 1  d F is the 1D inverse Fourier transform along the radial direction that transform the k-space measurements d into a sinogram. The second type of HYPR weights the composite image by the ratio of the backprojection of the 1D inverse Fourier transform of each ray of k-space data and the backprojection of the composite or reference image's projection B(P(r(x))) B(F (d)) ρ r(x)   1 ~ 1    d . In this paper, only the second type HYPR is investigated and this type of HYPR is a specific case of HYPR-LR (121) (refer to next section). Another type of approach for radial k-space reconstruction is the Projection Reconstruction FOCal Underdetermined System Solver (PR-FOCUSS) (123), which was designed to minimize the L1 image norm to constrain the image sparsity. This is a compressed sensing type of method. 72 HYPR-LR and PR-FOCUSS deal with each time frame separately (here denoted as "x-t methods"), and use a composite image (or images) reconstructed from multiple time frames as a reference image. Methods such as k-t BLAST/k-t SENSE (27) include spatiotemporal correlations in a different way, by working in the x-f domain (denoted as "x-f method"). k-t FOCUSS (28,124), which is the reweighted L2 norm implementation of compressed sensing, has been shown to be an extension of k-t BLAST/k-t SENSE for general sampling patterns. k-t FOCUSS constrains the image sparsity in the x-f domain. In this work, it is shown that HYPR-LR, PR-FOCUSS, k-t BLAST/k-t SENSE and k-t FOCUSS can be included in a unified multiplicative correction framework. As a consequence, they are all susceptible to errors in the reference image caused by signal zeroing. Previously, k-t SENSE and k-t FOCUSS have been presented in a unified framework (124), but HYPR-LR and PR-FOCUSS have not previously been shown to have specific relationships to each other. The theory is presented first. Demonstrations of the differences of these methods are then presented as well as similarities. Real data examples are also used to show the tradeoffs of different methods, in particular, if the reconstruction operates in the x-t or x-f domain. 6.2 The Extended GS Model 6.2.1 Generalized Reference Framework The MRI signal equation relates acquired k-space samples d(k ) to image domain object (x)   as follows: d k x e dx j k x            2 ( ) ( ) [6.1] In matrix form, this can be written as 73 d  Eρ [6.2] where E is the encoding matrix, and d, ρ are column vectors of k-space data and image pixel values, respectively. The elements of the encoding matrix are given by:   2 , i i j k k x e    x E The generalized series (GS) approach suggests that the image may be obtained by multiplicative correction of the reference or composite image r(x)  . The following model for image reconstruction was proposed (47,118,119,125):      terms i N i j k ic e 1 2 x ρ(x) r(x) m(x) r(x)        [6.3] where the reference image r(x)  is typically a magnitude only image (119), multiplication  is a pixel-by-pixel operation here, and ci is the ith GS basis coefficient. The basis functions can be general (118), but a Fourier basis is typically used (119). Nterms is the number of basis coefficients and typically matches the number of k-space samples (125). Note that this formulation (with an additional additive reference image term) can represent at least 14 different reconstruction methods such as feature-recognizing MRI (126) and RIGR, as shown in (125). Written in matrix form, equation [6.3] becomes: ρ Rm RE c H   [6.4] where R is the diagonal matrix with r on the main diagonal; and superscript H is the Hermitian transpose operator. In the original GS method and other methods listed in (125), r is limited to the image domain (x-t space). Here r is generalized to other domains, such as the x-f domain; correspondingly, the encoding matrix E is generalized to encode in the appropriate 74 domain. This new more general framework is referred to here as "generalized reference framework." For Cartesian acquisitions undersampled in the phase encode direction, E WF where W is a binary diagonal matrix choosing which k-space points that are sampled, andFis the forward Fourier transformation matrix. To adapt the formulation for non-Cartesian sampling, interpolation should be included, E=TF, H H H E  F T , where T is the interpolation and resampling matrix. To include coil sensitivity profiles, the encoding matrix may be extended to include coil sensitivity profiles as follows:              Nc 2 1 TFS TFS TFS E  , where i S are the diagonal matrix with values of ith coil sensitivity on the main diagonal ( c i 1,...,N ). For x-f domain methods, the encoding matrix is k t E  E F , where Ft is used to Fourier transform images from the x-f domain to the x-t domain, and Ek is used to transform images from the x-t domain to the k-t domain. Substituting equation [6.4] into [6.2], we have ERE c  d H . If a regularization term is added to provide more resistance to noise (127), (ERE I)c  d H , or equivalently c (ERE I) d    H [6.5] In the original GS paper, no regularization term was added, and c was found using a Toeplitz solver (118). Other methods have been used for calculating c, such as a singular value decomposition (SVD) method (128). Here we include a regularization term and solve for c using conjugate gradient methods. Substituting back into equation [6.4] gives 75 ρ RE (ERE I) d ~ 1   H H [6.6] This has a related mathematical form as noted in (129) ρ (E E R ) E d ~ H 1 1 H   [6.7] The next sections will show how these equations [6.5-6.7] relate to more recent reconstruction techniques. 6.2.2 Relation to PR-FOCUSS, k-t FOCUSS, and k-t BLAST/k-t SENSE PR-FOCUSS is a specific case of equation [6.6] when the reference image r is in x-t space and is updated at each iteration. PR-FOCUSS uses radial sampling and the inverse Fourier transform (IFT) and projection operations for encoding. Also very similar to equation [6.6] is k-t FOCUSS, although for this reconstruction method, r is in the x-f domain and is updated at each iteration. k-t FOCUSS has been shown to generalize k-t BLAST/k-t SENSE (124), and thus only k-t FOCUSS is referred to in this section. From equation [6.6], when a baseline term is added, we get ρ ρ RE (ERE I) (d Eρ ) 0 0       ~ H H [6.8] Equation [6.8] is the same as the k-t FOCUSS equation [17] of (124), when p in (124) is set to be ½ to provide minimization of the L1 norm of x-f space and the pseudo inverse is solved using the conjugate gradient method. 76 6.2.3 Relation to HYPR-LR In order to express HYPR-LR in the generalized reference framework, we multiply both sides of equation [6.2] by H E with a density compensation function D and substitute equation [6.3], which gives E Dd E DEρ E DE(r(x) m(x))      H H H ,here the multiplication operator is a pixel-by-pixel operation. E DE H can be denoted as a kernel h E Dd h (r(x) m(x))      H [6.9] When h is a good approximation to the Dirac delta function (upon proper choice of density compensation function D) or m is a constant: h(rm)  (hr)m If m does not change significantly over the effective spatial support of h, this relationship forms a good approximation, as has been exploited for homodyne detection (130) and fast RIGR(131). That is, h(r m)  (hr)m [6.10] From this approximation and equation [6.7], E Dd (h r(x)) m(x) (E DEr(x)) m(x)          H H Substituting into equation [6.3], E DEr E Dd ρ r(x) m(x) r(x) H H     ~    [6.11] where the multiplication and division operators are pixel-by-pixel operations. From equation [6.11], when r is set to be the composite image(s) which has high resolution, and d is the undersampled radial data, the equation is similar to the HYPR-LR method. The only difference is that equation [6.11] is written with gridding or inverse gridding with a density compensation function, rather than the notation in HYPR-LR where a projection or filtered backprojection image convolved with a low pass filter 77 (121) is used. Compared to the fast RIGR method (131), h is slightly different: HYPR-LR uses the filtered backprojection of projection convolved with a lowpass filter; while fast RIGR uses the inverse Fourier transform of the low resolution (truncated in k-space) data. These are essentially the same, but fast RIGR was originally created for keyhole-type undersamplings, while HYPR-LR was designed for an undersampling pattern spread throughout k-space. In order for results from the multiplicative correction method of equation [6.6] and HYPR-LR to be most similar, equation [6.10] should be a good approximation. In order to obtain a locally smooth m over the effective spatial support of h, the terms E Dd H and E DEr H should be locally smooth. This is consistent with the fact that HYPR-LR uses a filtered backprojection with a low frequency filter; in order to reduce the effective spatial support of h, a Gaussian filter was applied in (121). Note that the phase of m is usually locally smooth and can be included in m, and r is typically set to be the magnitude image without a phase term for HYPR-LR. When | r | 1 h  , HYPR-LR is the same as the second type of HYPR as mentioned above. 6.2.4 Sparsity in x-f and x-t Domains From compressed sensing theory, the image or its transformation should be sparse. It was reported that a sparse image is necessary for HYPR-LR and k-t FOCUSS to be effective (124,132). Strategies can be applied to make the images sparser. Current literature has applied a DC baseline image to enhance the sparsity, such as in (27,124). Other methods have also been proposed to get the baseline, such as RIGR (133), and a motion estimation/compensation scheme (28). 78 k-t FOCUSS essentially minimizes the L1 norm of the image in the x-f domain by minimizing the reweighted L2 norm iteratively, which assumes image sparsity in the x-f domain. The question arises as to how well this works when the image in the x-t domain is constrained instead. Here we compare the x-f method and the x-t method (the x-f method is denoted as "k-t FOCUSS," and the x-t method is denoted as "x-t FOCUSS," which is the same as PR-FOCUSS when the sampling pattern is radial) from the generalized reference framework perspective. For k-t FOCUSS, equation [6.3] can be written as ( ) x f x f t k k t ρ r F E c      H H [6.12] where H t F is the Fourier transform in time dimension, which will change x-t domain to x-f domain, and k E is the transformation operator that changes the x-t domain to the k-t domain. This equation sets up the minimization problem for the first iteration of k-t FOCUSS. Taking the inverse Fourier transform in the time dimension, we get an alternate expression for k-t FOCUSS: ρ r (E c ) xt xt k kt   H [6.13] On the other hand, for<