Study of Phase Information in MRI with Applications to Fast Imaging and Inversion Recovery Imaging by ZHENG CHANG B.Sc., Shandong University, 2000 M . S c , The University of British Columbia, 2002 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF DOCTOR OF PHILOSOPHY OF SCIENCE In T H E F A C U L T Y OF G R A D U A T E STUDIES (PHYSICS) T H E UNIVERSITY OF BRITISH C O L U M B I A November 2005 © Zheng Chang, 2005 Abstract In Magnetic Resonance Imaging ( M R I ) , signal is complex valued and can be represented by a magnitude image and a phase map. more frequently Although magnitude images are used much in clinical diagnosis than phase maps, the latter should not be undervalued because a lot of valuable information is encoded only in phase, which has great potential in medical applications. M y doctoral thesis focuses on phase information in M R I and its applications in the areas of fast imaging and Inversion Recovery (IR) imaging. IR imaging is one of the most useful techniques in M R I contrast manipulation, but its use is limited because of the presence of phase errors. The thesis proposes a new method to correct phase errors that occur in IR imaging. The method models phase variation with a polynomial, whose coefficients are statistically determined by calculating relative vector rotations of complex fields after being shifted by n pixels. A s discovered in this thesis, increasing the pixel shift effectively enhances phase signal without amplifying the corresponding noise, and thereby improves phase correction. The method has been successfully demonstrated with 2 D in vivo IR imaging data. Phase information has great potential also in fast imaging in M R I . For example, a recently published fast imaging method named "Skipped Phase Encoding and Edge Deghosting ( S P E E D ) " conducts strategic spatial phase encoding and thereby accelerates ii M R I . S P E E D is promising but is demonstrated so far only in 2 D and only with a single coil. This thesis presents new developments based on the principle of S P E E D : First, S P E E D is extended from 2 D to 3 D to reduce scan time with more flexibility and efficiency; Second, S P E E D is combined with Half-Fourier imaging to accelerate M R I further by a factor of nearly 2; Third, S P E E D is simplified based on the sparseness of M R angiography data; Fourth, a new parallel imaging strategy named " S P E E D with Array C o i l Enhancement ( S P E E D - A C E ) " is proposed, which extends S P E E D from a single coil to multiple coils to further accelerate M R I . Finally, signal-to-noise ratio ( S N R ) in S P E E D is studied, and a novel approach for S N R improvement is proposed and demonstrated with computer simulations and in vivo data. iii Table of Contents Abstract ii Table of Contents iv List of Tables viii List of Figures ix List of Abbreviations xix Acknowledgement xxi Chapter 1 Imaging Goals of the Thesis and Fundamentals of Magnetic Resonance 1 1.1 The Outline and the Goals of the Thesis 2 1.2 A Brief history of N M R and M R I 4 1.3 Basic Physics of N M R 6 1.3.1 Nuclear Magnetization 6 1.3.2 Motion of the Magnetization: Quantum to Classical 11 1.3.3 R F Excitation and N M R Signal Detection 13 1.3.4 Relaxation Times 16 1.4 1.5 Basic Principles of Magnetic Resonance Imaging 19 1.4.1 Magnetic Field Gradient 20 1.4.2 Slice Selection 20 1.4.3 Frequency Encoding 22 1.4.4 Phase Encoding 23 1.4.5 The Concept of £-Space 24 1.4.6 Pulse Sequences 26 Summary 31 iv Chapter 2 Review of Applications of MRI Phase Information 32 2.1 M R I Phase Information and Phase Wrapping 32 2.2 Review of Applications of M R I Phase Information 34 2.2.1 F l o w Imaging 34 2.2.2 Chemical Shift Imaging 37 2.2.3 Fast Imaging in M R I 45 2.3 Objectives of the Thesis 49 Chapter 3 IR Imaging with Non-linear Phase Correction based on an Extended Statistical Algorithm 50 3.1 Introduction 50 3.2 Method 52 3.2.1 Review of the A C Method 52 3.2.2 N-Pixel-Shift Rotational Differential Field 56 3.2.3 Extended A C Method 59 3.2.4 Improved Weighted Least Squares Method 67 3.2.5 IR Image Reconstruction with Phase Correction 71 3.2.6 Experiments 72 3.3 Results 72 3.4 Discussion 79 Chapter 4 Simplified Skipped Phase Encoding and Edge Deghosting (SPEED) with Applications in MRA 82 4.1 Introduction 84 4.2 Method 85 4.2.1 Simplified S P E E D for Fast M R A 85 4.2.2 Extension of Simplified S P E E D to T w o P E Directions 93 4.2.3 Combination o f Simplified S P E E D and Half-Fourier Reconstruction. 95 v 4.2.4 Generalization of Simplified S P E E D 4.2.5 Experiments 96 101 4.3 Results 103 4.4 Discussion 117 Chapter 5 Accelerating MRI by Skipped Phase Encoding and Edge Deghosting (SPEED) and its Extensions 120 5.1 Introduction 120 5.2 Methods 120 5.2.1 Review o f the Original S P E E D : 120 5.2.2 Extension of S P E E D to T w o P E directions 129 5.2.3 Combination of S P E E D and Half-Fourier Imaging 132 5.2.4 Generalization of S P E E D 135 5.3 Results 136 5.4 Summary 142 Chapter 6 Highly Accelerated MRI by Skipped Phase Encoding and Edge Deghosting with Array Coil Enhancement (SPEED-ACE) 143 6.1 Introduction 143 6.2 Methods 144 6.2.1 Review of S E N S E 144 6.2.2 Review of S P E E D 146 6.2.3 SPEED-ACE 147 6.2.4 S P E E D - A C E with Multiple Skipped Samplings 156 6.2.5 Generalized S P E E D - A C E with Multi-layer Models 157 6.2.6 Experiments 158 6.3 Results 159 6.4 Discussion 171 vi Chapter 7 SNR Improvement for SPEED and Its Extensions 174 7.1 Introduction 175 7.2 Methods 176 7.2.1 Shepp-Logan phantom 176 7.2.2 Sampling with Skipped Phase Encoding 177 7.2.3 S N R Improved by Data Leakage of D F T in S P E E D 179 7.2.4 Proposed Method to Improve S N R for S P E E D and its Extensions 185 7.2.5 Experiments 188 7.3 Results 188 7.4 Summary 191 Chapter 8 Conclusions and Recommendations for Future Work 192 8.1 Conclusions o f the Thesis 192 8.2 Recommendations for Future W o r k 194 8.2.1 Other Applications of Phase Correction with N-Pixel-Shift R D F . . 194 8.2.2 Highly Accelerated M R A with Simplified S P E E D - A C E 195 8.2.3 Applications of S P E E D in Multiple Acquisitions 196 8.2.4 S N R Improvement by Using S P E E D 197 Bibliography 198 vii List of Tables Table 3.1: W R P V values for a series of R D F s are obtained by repeated application of E q . (3.10) with the total shifts of 125-pixel and 134-pixel, respectively, in x and y directions for the squared complex IR image from a transverse head scan 75 Table 3.2: Refocusing ratios for extended A C method applied to the squared complex IR image shown in Figure 3.4 (c) with the corresponding polynomials 75 Table 3.3: Refocusing ratios for various phase correction methods applied to the squared complex IR image from a transverse head scan shown in Figure 3.4 (c) 76 Table 3.4: Refocusing ratios for various phase correction methods applied to the squared complex IR image from a coronal head scan 78 Table 6.1: T R E values for the images reconstructed by generalized S P E E D - A C E shown in the bottom four rows of Figure 6.9 ; 170 viii List of Figures Figure 1.1: Zeeman splitting increases linearly with the magnetic field Bo. Spin up protons have lower potential energy than spin down protons; the former are aligned parallel with the magnetic field, while the latter are anti-parallel. The energy difference between the two states is AE = hco [11] 8 0 Figure 1.2: The motion of M is a precession about the static magnetic field Bo, just like the precession of a spinning top in a gravitational field 13 Figure 1.3: The magnetization vector M precesses about the effective magnetic field B\ i n the rotating frame [11,12] 15 Figure 1.4: Time evolution of the magnetization M: (a) is the T 2 decay of transverse component of M according to E q . (1.23); (b) is the recovery of longitudinal magnetization according to E q . (1.24) [11,13] 19 Figure 1.5: A typical pulse sequence of slice selection using a 90° R F pulse and a z directional linear magnetic field gradient G ....21 z Figure 1.6: A pulse sequence of frequency encoding using a F E gradient G after a typical slice selection 23 x Figure 1.7: A pulse sequence containing a slice selection, P E with a linear magnetic field gradient G , and F E with G 24 y x Figure 1.8: A typical spin-echo pulse sequence of M R I 27 Figure 1.9: Inversion recovery pulse sequence 28 Figure 1.10: A typical 3 D gradient-echo pulse sequence 29 Figure 1.11: The original E P I pulse sequence (a) and its k-space trajectory (b) [11]. 30 Figure 2.1: A n example obtained from a transverse gradient-echo head scan, (a) is the magnitude image representing the amplitude of the transverse magnetization, and (b) is the corresponding phase map that shows the direction of the transverse magnetization.. 33 Figure 2.2: A bipolar gradient used for velocity phase sensitization [11] ix 35 Figure 2.3: Chemical shift information of water and fat can be encoded with a spin-echo sequence (a) or a gradient-echo sequence (b) by manipulating the timing parameters: the refocusing time to, the echo time TE, and the time shift Af [30] 41 Figure 2.4: Half-Fourier is demonstrated with in vivo data from a spin-echo transverse head scan, (a) is the magnitude image reconstructed from full k-space data, and is therefore treated as the "gold standard", (b) is the recovered and desirable image obtained from 62.5% data by using Half Fourier imaging 48 Figure 3.1: (a) is real part of a simulated complex M R image denoted as S(x,y) that is obtained by imposing a noise-free phase error 1.2 + 0.03y rad onto a magnitude image, (b) is complex pixel distribution of S(x,y). (c) and (d) are complex pixel distributions of A y , and A y 3 0 o f S(x,y), respectively, (e) is noisy M R image denoted as C(x,y) that is the addition of S(x,y) and a complex Gaussian noise field [69]. (f-h) are repeated calculations of (b-d) for C(x,y). The angular estimation in (g) is less reliable because of relative large angular uncertainty, while the angular estimation in (h) is much more accurate because the relative angular uncertainty is greatly reduced by enhancing the angle by a factor of 30. This is the key idea of the proposed method 58 Figure 3.2: Flowchart of the extended A C method for non-linear phase error correction based on the n-pixel-shift R D F 66 Figure 3.3: Flowchart of the improved W L S method based on the n-pixel-shift R D F . . . . 70 Figure 3.4: (a) and (b) are the real part and phase of the complex experimental IR image of a transverse head scan, (c) and (d) are the real part and phase of the squared complex IR image, in which polarity factor is removed but the phase is doubled 73 Figure 3.5: In vivo real IR image from a transverse head scan and its corresponding phase corrections with various methods 76 Figure 3.6: In vivo real IR image and its corresponding phase corrected images with various methods from a coronal head scan 78 Figure 4.1: K-space sampling scheme and data flow of Simplified S P E E D with one P E direction. 2 D full k-space S (k) is sparsely sampled into two interleaved sparse k-space data sets: S,(k) and S (k), each with the same skip size of N and a different relative shift in P E . After inverse 2 D F T , 2 ghosted images —I,(r) and I (r)— are obtained and subsequently combined into a deghosted image In(r) 86 0 2 2 x Figure 4.2: Data flow of deghosting in Simplified S P E E D with one P E direction at N = 3. A L S E solution G„(r) and its associated ghost order map n(r) are firstly solved with two ghosted images: I,(r) and I (r). Then, G is sorted out according to n(r) and thereby yields N = 3 separate deghosted M R A images: G ( r ) , G , ( r ) , and G ( r ) . Finally, they are spatially registered to yield G ^ r ) , G , ( r ) , and G ( r ) and averaged into a deghosted image I (r) 91 2 n 0 r 2 r2 0 Figure 4.3: Flowchart of Simplified S P E E D with one P E direction used to achieve an acceleration factor of N/2 with a single coil 92 Figure 4.4: K-space sampling scheme and data flow of Simplified S P E E D with two P E directions. 2D cross-section of 3D k-space S (k) is sparsely sampled into two interleaved sparse k-space data sets: S (k) and S ( k ) , each with a skip size of M in horizontal direction and a skip size of N in vertical direction. After inverse 2D FT, two ghosted images I,(r) and I (r) are obtained and subsequently combined into a deghosted image Io(r) 94 0 { 2 2 Figure 4.5: The photo of the vasculature phantom (left) and the coverage of the 25 slices. 102 Figure 4.6: 2D reconstructed images from the phantom scan, (a) is one of the 25 slices reconstructed from full k-space data and is considered as the gold standard, (b) is a noise field where each pixel value is randomly taken from the background of the real and imaginary parts of the gold standard shown in (a), (c) is the summation of (a) and (b). (c) has a T R E of 0.10. (d) is one of the ghosted images, each obtained by sampling kspace with a skip size of N = 7 along vertical direction, (e) is the deghosted image yielded by Simplified S P E E D based on a single-layer model with an acceleration factor of 7/2 (f) is the residual map of (e). (g) is the deghosted image yielded by Simplified S P E E D based on a double-layer model with an acceleration factor of 7/3 (h) is the residual map of (g). Compared with the gold standard (a), both (e) and (g) have the errors comparable to the noise level of the gold standard quantified by T R E . (e) has a T R E of 0.17 while (g) has a T R E of 0.09. To allow more information to be observed, the contrast in (a), (c), (e) and (g) are reduced with a lower window and level setting (-100, 100) as shown in (i), (j), (k) and (1), respectively 105 Figure 4.7: Results of phantom scan reconstructed from three interleaved data sets with a skip size of N = 11 along vertical direction, (a) is the gold standard, which is the same as Figure 4.6 (a), (b) is ghosted image reconstructed by direct 2 D F T from one of the three skipped k-space data sets, (c) is the deghosted image yielded by Simplified S P E E D based on a double-layer model with an acceleration factor of 11/3 ~ 3.7. (d) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a double-layer xi model with an acceleration factor close to 6. (e-h) are corresponding M I P images of the 25 slices reconstructed in the same way as (a-d). Both (c) and (d) have the errors comparable to the noise level of the gold standard (a), and their corresponding M I P images (g) and (h) have quite good image qualities, (c) has a T R E of 0.18 while (d) has a T R E of 0.12. W i t h a lower window and level setting (-100, 100), the contrasts in (a-d) are reduced to allow more information to be observed as shown in (i-1), respectively.. 107 Figure 4.8: Results of phantom scan reconstructed from three interleaved data sets with a skip size of N = 5 along vertical direction and a skip size of M = 2 along horizontal direction. Figure 4.7 (a) and (e) are respectively the gold standard and the corresponding M I P image, (a) is one of the three ghosted images, (b) is the deghosted image reconstructed by Simplified S P E E D based on a double-layer model with an acceleration factor of 10/3. (c) is the deghosted image yielded by combination of Simplified S P E E D and Half-Fourier imaging based on a double-layer model with an acceleration factor of 5.5. (d-f) are corresponding M I P images of the 25 slices reconstructed in the same way as (a-c). The reconstructed images (b-c) and their corresponding M I P images (e-f) are all have reasonably good image qualities. Compared with the gold standard shown in Figure 4.7 (a), (b) and (c) are respectively measured to have T R E values of 0.13 and 0.20. For easy visualization, the contrasts in (a), (b) and (c) are reduced with a lower window and level setting (-100, 100) as shown in (g), (h) and (i), respectively 109 Figure 4.9: (a) is the gold standard image with full k-space sampling, (b) is what (a) becomes after being superimposed by noise points randomly taken from its background, resulting in a T R E of 0.67. T o reveal more information, the contrast in (a) and (b) are reduced with a lower window and level setting (-25, 25) as shown in (c) and (d), respectively 110 Figure 4.10: Results of the in vivo P C head scan reconstructed from various interleaved data sets, each sampled with a skip size of N = 5 and a different relative P E shift along vertical direction within 60% of k-space. (a) is the gold standard, which is identical to Figure 4.9 (a), (b) is the ghosted image reconstructed by direct 2 D F T from one of the interleaved k-space data sets, (c) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a single-layer model with an acceleration factor of 4. (d) is the deghosted image based on a double-layer model with an acceleration factor of 2.8. (c) and (d) are respectively measured to have T R E values of 0.62 and 0.51. (e-h) are corresponding axial M I P images of the 25 slices reconstructed in the same way as (a-d). (i-1) are the M I P images that are similar to (e-h) but at a different view angle. Compared with (a), (c) and (d) are both reasonably good while scan time has been highly reduced, (m-p) are the same as (a-d) but with a lower scale setting (-25, 25) 112 Figure 4.11: Results of the in vivo P C head scan reconstructed from various interleaved data sets, each sampled with a skip size of N = 7 and a different relative P E shift along xii vertical direction within 60% of k-space. Figure 4.10 (a) shows the gold standard; Figure 4.10 (e) and (i) are the corresponding M I P images of the gold standard at two different view angles, (a) is one of the ghosted images, (b) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a single-layer model with an acceleration factor close to 6. (c) is the deghosted image yielded based on a double-layer model with an acceleration factor close to 4. (b) has a T R E of 0.68 while (c) has a T R E of 0.59. (d-f) are corresponding axial M I P images of the 25 slices reconstructed in the same way as (ac). (g-i) are the M I P images that are similar to (d-f) but at a different view angle. Both (b) and (c) have the errors comparable to the noise level of the gold standard but with high compression factors. For easy visualization, the contrasts i n (a-c) are reduced with a lower window and level setting (-25, 25) as shown in (j-1), respectively 114 Figure 4.12: Results of the in vivo P C head scan reconstructed from various interleaved and relatively shifted data sets, each sampled respectively with a skip size of N = 3 along vertical direction and a skip size of M = 2 along horizontal direction within 60% of kspace i n vertical direction. Figure 4.10 shows the gold standard and the corresponding M I P images, (a) is one of the ghosted images, (b) is the deghosted image yielded by HalfFourier Simplified S P E E D based on a single-layer model with an acceleration factor of 5. (c) is the deghosted image yielded based on a double-layer model with an acceleration factor of 3.3. (b) has a T R E of 0.67 while (c) has a T R E of 0.56. (d-f) are corresponding axial M I P images of the 25 slices reconstructed i n the same way as (a-c). (g-i) are the M I P images that are similar to (d-f) but at a different view angle. A l l the aliasing ghosts in (a) have been successfully unfolded and the results shown i n (b) and (c) are both visually good, (j-1) are the same as (a-c) but with a lower scale setting (-25, 25) 116 Figure 5.1: K-space sampling scheme and data flow of 2D S P E E D [3]. 2 D full k-space S (k) is sampled into three interleaved sparse k-space data sets: S,(k), S (k) and S (k), each with a skip size of N and a relative shift i n P E . After inverse 2 D F T followed by a differential filtering, S ^ k ) , S (k), and S (k) are firstly turned into three ghosted images: I^r), I (r) and I (r), and then three ghosted edge maps: E,(r), E (r) and E (r). Finally, they are deghosted into a ghost-free edge map Eo(r), which is subsequently inversefiltered to yield a deghosted image Io(r) 122 0 2 2 2 3 3 3 2 3 Figure 5.2: Flowchart of the S P E E D with one P E direction used to achieve an undersampling factor of 1.7 with a single coil [3] 128 Figure 5.3: K-space sampling scheme and data flow of S P E E D with two P E directions. F u l l k-space S (k) is sampled into three interleaved sparse k-space data sets: S,(k), S (k) and S (k), each with a skip size of N along one P E direction and a different skip size of M along the other P E direction. After inverse 2 D F T followed by a differential filtering, S,(k), S (k) and S (k) are first turned into three ghosted images: I,(r), I (r) and I (r), and 0 2 3 2 3 2 xiii 3 then three ghosted edge maps: E,(r), E (r) and E (r). Finally, they are deghosted into a ghost-free edge map En(r), which is subsequently inverse-filtered into a deghosted image Io(r) 130 2 3 Figure 5.4: K-space sampling scheme and data flow of Half-Fourier-SPEED. Slightly more than half k-space S (k) is sampled into three interleaved k-space data sets: S,(k), S (k), and S (k), each with a skip size of N and a relative shift in P E . After inverse 2 D F T followed by a differential filtering, S^k), S (k) and S (k) are first turned into three ghosted images: Ij(r), I (r), and I (r), and then three ghosted edge maps: Ej(r), E (r), and E (r). They are deghosted into a ghost-free edge map En(r), which is subsequently inverse-filtered to yield a deghosted image Io(r). Finally, In(r) is processed by HalfFourier reconstruction to yield a final image If, ai(i") 134 0 2 3 2 2 3 3 2 3 n Figure 5.5: Results from S P E E D with one P E direction and a standard reconstruction for a spin-echo transverse head scan, (a) is the gold standard with full k-space sampling, (b) is one of the three ghosted images with P E skip size N = 5, having up to five layers of ghost overlapping, (c) is the corresponding edge map after high-pass filtering, with ghost overlapping layers being reduced to be only 2 in most pixels, (d) is the edge map after deghosting. (e) is the final deghosted image by inverse filtering (d). (e) is measured to have a T R E of 0.07. W i t h an undersampling factor of 5/3 ~ 1.7, the final image (e) looks comparable to the gold standard (a), (f) is a residual map for quality control 137 Figure 5.6: Results from S P E E D with two P E directions and a standard reconstruction for a spin-echo transverse head scan, (a) is the gold standard with full k-space sampling, (b) is one of the three ghosted images with a skip size N = 3 along vertical direction and another skip size M = 2 along horizontal direction, (c) is the corresponding edge map. (d) is the final deghosted image, which is comparable to the gold standard (a), (d) is measured to have a T R E of 0.09. In this case, the undersampling factor is 6/3 = 2...... 139 Figure 5.7: Results from Half-Fourier-SPEED and conventional Half-Fourier imaging for a spin-echo sagittal knee scan and an IR coronal head scan. Although scan time has been further reduced by a factor of 3/5, the images (b) and (d) reconstructed by Half-FourierS P E E D have no visible difference from the gold standard (a) and (c), each obtained by conventional Half-Fourier reconstruction, (b) and (d) are respectively measured to have T R E values of 0.03 and 0.05, which are both very small 141 Figure 6.1: K-space sampling scheme and data flow of 2 D S E N S E with four receiver coils. 2 D full k-space is selectively sampled with a skip size of N by using four receiver coils i n parallel. The sampled sparse data sets are reconstructed separately into four ghosted images: Ii(r), L:(r), h(r), and L;(r). Unfolding signal superposition yields a final deghosted image Io(r) 145 xiv Figure 6.2: K-space sampling scheme and data flow of 2 D S P E E D - A C E . B y using four receiver coils in parallel, k-space is selectively sampled with a skip size of N into four sparse data sets: S^k), S (k), S (k) and S (k). They are then reconstructed separately by direct inverse 2 D F T into four ghosted images: Ii(r), l2(r), l3(r), and U(r), with ghost overlapping up to N layers. After high-pass filtering, four ghosted edge maps — Ei(r), E (r), E (r), and E4(r) — are obtained and subsequently combined into a deghosted edge map Eo(r). Finally, Eo(r) is inverse-filtered into a deghosted image In(r) 148 2 2 3 4 3 Figure 6.3: (a) is an M R image from a typical spin-echo scan on a clinical 1.5 T scanner equipped with multiple receiver coils and is denoted as C(x,y). (b) is the corresponding coil sensitivity map for one receiver coil and is denoted as S(x,y), where signal region is surrounded by noise background. S(x,y) was slightly smoothed to improve signal-to-noise ratio [49]. If C(x,y) is sensitivity weighted by S(x,y), the sensitivityweighted image can be written as S(x, y)C(x, y). If the vertical direction is defined as the y-direction, the v-directional spatial differential of S(x,y)C(x,y) consists of two terms: S(x,y+\)C(x,y+\)S(x,y)C(x,y) = C(x,y+l)[S(x,y+l)-S(x,y)] + S(x,y)[C(x,y+\)-C(x,y)l W i t h the same scale, (c) shows the first term, and (d) shows the second term. Compared with (d), C(x,y+l)[S(x,y+l)-S(x,y)] shown in (c) can be assumed to be negligible. This is because coil sensitivity profile S(x,y) is varying slowly. For easy visualization, all images are shown here in magnitude although they are actually complex 150 Figure 6.4: Flowchart of S P E E D - A C E with single skipped sampling used to achieve an undersampling factor of 5 with four receiver coils 155 Figure 6.5: Results from the second slice of the pelvis scan. In the top row are individual low-resolution coil images obtained from 32 central k-space lines. In the middle row are individual sensitivity-weighted and ghosted images, each reconstructed from the data set sampled with a skip size of 5. Aliasing ghosts of various orders are seen along vertical P E direction. In the bottom row are corresponding ghosted edge maps of middle row and shows much less ghost overlapping 160 Figure 6.6: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with double skipped samplings, each accomplished by using a skip size of N = 5 steps and four receiver coils in parallel, (a) is the desired body image considered as the gold standard, (b) is what (a) becomes after being superimposed by noise points randomly taken from the background of the real and imaginary parts of the gold standard shown in (a), resulting in a T R E of 0.12. (c) is one of the eight ghosted images, (d) is the R M S of the four individual low-resolution coil images shown in the top row of Figure 6.5. (e) is the R M S of the coil images obtained by a direct inverse F T of all sampled data, (f) is the reconstructed images by S P E E D - A C E based on a bilayer model with an undersampling factor of 2.5. (f) has a T R E 0.082, and appears to be comparable to the gold standard (a) 162 xv Figure 6.7: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with single skipped sampling. S P E E D - A C E samples k-space with a skip size of N = 3 in P E by using four receiver coils in parallel and reconstructs deghosted images with multiple-layer models. Figure 6.6 (a) is the gold standard reconstructed from full kspace data, (a) is one of the four individual sensitivity-weighted and ghosted images, (b) is corresponding ghosted edge map of (a), (c) is the R M S of coil images obtained by a direct inverse F T of all sampled data, (d) and (e) are the reconstructed images by S P E E D A C E based respectively on a single-layer model and a double-layer one. (d) has a T R E 0.09 while (e) has a T R E of 0.11. (f) is the simple average of (d) and (e). W i t h an undersampling factor of 3, all the images shown in (d-f) are quite good compared with the gold standard shown in Figure 6.6 (a) 164 Figure 6.8: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with single skipped sampling. S P E E D - A C E samples k-space with a skip size of N = 5 in P E by using four receiver coils and reconstructs deghosted images with multi-layer models. Figure 6.6 (a) is the gold standard reconstructed from full k-space data, (a) is one of the four individual sensitivity-weighted and ghosted images, (b) is the R M S of the coil images obtained by direct inverse F T of all sampled data, (c-f) are the reconstructed images by S P E E D - A C E based respectively on zero-layer, single-layer, double-layer, and triple-layer models, (g) is the final image reconstructed by generalized S P E E D - A C E and has a T R E of 0.12. W i t h an undersampling factor of 5, the final image shown in (g) has the errors comparable to the noise level of the gold standard shown in Figure 6.6 (a), (h-k) are automatically generated model maps for (g), which indicate the number of overlapping layers in the model used at each pixel. M o r e specifically, (h) is the zerolayer model map, in which a white pixel has a value of unity representing a zero-layer model applied at that pixel. Likewise, (i-k) are respectively the single-layer, doublelayer, and triple-layer model maps, (i) is the residual map of (g) 168 Figure 6.9: Results from the 20-slice pelvis scan: In the top four rows are the gold standard from the 1 slice to the 2 0 slice, each reconstructed from full k-space data. In the bottom four rows are reconstructed images by generalized S P E E D - A C E with an undersampling factor of 5 169 st th Figure 7.1: Simulated M R images with a digital Shepp-Logan phantom: (a) is the ideal noise-free complex M R image; (b) is the noisy M R image imposed with simulated "white"- Gaussian noise 177 Figure 7.2: Ghosted images obtained with a skip size of N = 4 and their corresponding edge maps, (a) is noise-free ghosted image obtained by sampling k-space with a skip size of N = 4. (b) is the edge map of (a), (c) is simulated noisy ghosted image obtained with N = 4. (d) is the edge map of (c). It can be shown that every 1/4 of F O V is just the same as one another 178 xvi Figure 7.3: Ghosted images obtained with a skip size of N = 5 and corresponding edge maps, (a) is the ideal noise-free complex M R image, which is identical to Figure 4 (a), (b) is one of three ghosted images, (c) and (d) are respectively the edge maps of (a) and (b). (e) is profile of column N o . 128 of (d). It can be easily seen from (e) that every a 1/5 F O V is slightly different from one another 180 Figure 7.4: Noise-free deghosted images and edge maps yielded by S P E E D with a skip size of 5. (a) is the ideal noise-free complex M R image, which is the same as Figure 7.1 (a), (b) is the edge map of (a), (c) is one of the three ghosted edge maps, each obtained by sampling k-space with a skip size of 5. (d) is one of the five deghosted edge maps, (e) is the average of the five deghosted edge maps, (f) is the final deghosted image, which is obtained by inverse-filtering (e). (f) has a T R E of 0.007 183 Figure 7.5: Noisy deghosted images and edge maps yielded by S P E E D with a skip size of N = 5. (a) is the simulated noisy complex M R image, which is the same as Figure 7.1 (b) with an S N R of 11.5. (b) is the edge map of (a) with an S N R of 1.5. (c) is one of the five deghosted edge maps, and has an S N R of 1.3. (d) is the average of the five deghosted edge maps, and has an S N R of 1.5, which is the same as S N R of (b). (e) is obtained by inverse-filtering (c) and has an S N R of 9.0. (f) is the final deghosted image, which is obtained by inverse-filtering (d) and has an S N R of 11.5, which is the same as S N R of the original image in (a), (f) is measured to have a T R E of 0.08 184 Figure 7.6: Noisy deghosted images and edge maps yielded by S N R improved S P E E D with a skip size of N = 5. (a) is the simulated noisy complex M R image, which is the same as Figure 7.1 (b) with an S N R of 11.5. (b) is the edge map of (a) with an S N R of 1.5. (c) is one of the five deghosted edge maps, and has an S N R of 1.3. (d) is the average of the five deghosted edge maps, and has an S N R of 1.7, which is slightly greater than S N R of (b). (e) is obtained by inverse-filtering (c) and has an S N R of 9.0. (f) is the final deghosted image, which is obtained by inverse-filtering (d) and has an S N R of 12.0, which is slightly greater than the S N R of Figure 7.5 (f) obtained without randomizing (nl, n2). In addition, the T R E of (f) is measured to be 0.075, which is slightly less than that of Figure 7.5 (f) 187 Figure 7.7: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with single skipped sampling, which is accomplished by using a skip size of N = 4 steps and four receiver coils in parallel, (a) is the gold standard reconstructed from full kspace data., (b) is the image reconstructed by S P E E D - A C E based on a double-layermodel without S N R improvement, with an undersampling factor of 4. (c) is the image reconstructed by S P E E D - A C E based on a double-layer model with S N R improvement, with an undersampling factor of 4. (b) has a T R E of 0.18, and (c) has a T R E of 0.15. 189 xvii Figure 7.8: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with double skipped samplings, each accomplished by using a skip size of N = 8 steps and four receiver coils in parallel. Figure 7.7 (a) is the gold standard with full kspace sampling, (a) is one of the eight individual sensitivity-weighted and ghosted images, (b) is the corresponding ghosted edge map of (a), (c) is the reconstructed images by S P E E D - A C E based on a double-layer model with an undersampling factor of 4. (d) is reconstructed similar to (d) but with an S N R improvement, (c) has a T R E 0.14, while (d) has a T R E of 0.11, which is comparable to the noise level quantified by T R E in Figure 6.6 (b) 190 xviii List of Abbreviations This is a list of abbreviations that are used throughout the thesis. Contrast enhanced CE Cerebrospinal fluid CSF Chemical shift imaging CSI Chemical shift imaging with spectrum modeling CSISM Computed tomography CT Direct phase encoding DPE Echo planar imaging EPI Free-induction-decay FID Fourier transform FT Fast Fourier transform FFT Field-of-view FOV Ghost layer unfolding expansion GLUE Inversion recovery IR Least-square-error LSE Maximum-intensity-projection MIP M R angiography MRA Magnetic resonance imaging MRI Nuclear magnetic resonance NMR Phase-contrast PC Parallel imaging with localized sensitivities PILS Partially-opposed-phase POP Quadrature-sampling QS Phase encoding PE Root-mean-square RMS Projection reconstruction PR Radio-frequency RF xix Rotational differential field RDF Root-mean-square RMS Refocusing ratio RR Signal-to-noise ratio SNR Skipped Phase Encoding and Edge Deghosting SPEED S P E E D with array coil enhancement SPEED-ACE Inversion time TI Repetition time TR Echo time TE Time-resolved imaging of contrast kinetics TRICKS Total relative error TRE Weighted lease squares WLS Weighted relative phase variation WRPV xx Acknowledgement I would like to thank Dr. Qing-San X i a n g , my advisor, for encouraging me to explore the amazing world of M R I , inspiring me with his brilliant insights, and guiding me through my doctoral study. I would also like to thank Dr. A n n a Celler for sparking my interest in medical physics, helping me develop research skills, and offering constant encouragement. I am grateful to Drs. D a v i d F . Measday, A l e x L . M a c K a y , Vesna Sossi, Robert Rohling, and Brenda Clark for giving me invaluable guidance in medical physics. I want to give special thanks to Dr. Janis A . M c K e n n a for continuously helping me during my study at U B C . I thank all my colleagues who helped me in this work; in particular, I thank Wayne Patola and Ange Lalonde for their assistance in collecting in vivo data at St. Paul's hospital. I am also grateful to B C Children's Hospital for providing financial support. Finally, I wish to express my love and gratitude to my family. I thank my parents and my parents-in-law for their constant support and encouragement. Special thanks go to my wife, Peiyi Duan, for helping me with her wisdom and love. xxi Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging This chapter outlines the goals of the thesis, and then provides an overview of the phenomenon of nuclear magnetic resonance ( N M R ) and the fundamentals of magnetic resonance imaging ( M R I ) . The overview is intended to help readers understand the work presented in subsequent chapters. This chapter begins with the outline and the goals of the thesis. Then, a brief historical summary of N M R and M R I is provided. This is followed by quantum mechanical and classical descriptions of magnetization, which outline the origin of the N M R signal and the behavior of magnetization. Nest, the Bloch equations are introduced, R F excitation is described and N M R relaxation times are briefly reviewed. introduces an important concept: magnetic field gradient. Moreover, this chapter Based on the concept, M R I spatial localization techniques of slice selection, frequency encoding, and phase encoding are presented. Furthermore, spatially encoded M R signals are described in terms of "kspace"; that is, the Fourier domain used as a powerful analytical tool for understanding and developing M R I pulse sequences. Finally, several typical pulse sequences are briefly reviewed, which include gradient-echo imaging and echo planar imaging (EPI). pulse sequence of E P I is intuitively interpreted by using k-space. 1 The Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging 1.1 The Outline and the Goals of the Thesis M R I has become one of the most popular medical imaging modalities because it can provide an unparalleled view inside a human body without making any invasion into it. In M R I , signal consists of real and imaginary components, which can be used to derive a magnitude image and a phase map. The unique feature of M R I is phase, which often encodes important information such as field inhomogeneity, motion, and chemical shift. This implies great potential in medical applications such as flow imaging, water-fat imaging and fast imaging in M R I . M y doctoral thesis focuses on the study of phase information in M R I and explores its applications to fast imaging and Inversion Recovery (IR) imaging. IR is a useful technique for N M R T l - s i g n a l enhancement; however, IR imaging is limited by phase errors in M R I . Therefore, it is very desirable to conduct phase correction for IR imaging. The linear phase correction method proposed A h n and Cho ( A C method) [1] is efficient and robust; however, it is limited because phase often appears non-linear in M R I . In this work, I propose an extended A C method for nonlinear phase correction, which has been successfully demonstrated with in vivo IR data [2]. Phase information has great potential also i n fast imaging in M R I , which has been very desirable since M R I came into being. The recently proposed fast imaging method named "Skipped Phase Encoding and Edge Deghosting ( S P E E D ) " [3] is promising; however, it has been demonstrated so far only i n 2 D and only with a single coil. In this work, I present new development based on S P E E D [3]: First, S P E E D is extended from 2 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging 2 D to 3D and combined with Half-Fourier imaging for improved performance; Then, S P E E D is simplified based on the sparseness of M R angiography ( M R A ) data [4]; finally, a new parallel imaging strategy named " S P E E D with Array C o i l Enhancement ( S P E E D - A C E ) " is proposed, which extends S P E E D from a single coil to multiple coils to further accelerate M R I [5]. S P E E D - A C E can achieve an undersampling factor even greater than the number of receiver coils, which has so far not been for other existing parallel imaging methods [5]. To help readers understand my work, I outline the structures of the thesis as follows: Chapter 1 provides an overview of the phenomenon of N M R and the fundamentals of M R I . Chapter 2 introduces phase information in M R I and provides a brief review of its applications, revealing the motivation of the thesis. Chapter 3 presents a novel statistical non-linear phase correction method and demonstrates it with in vivo IR imaging data. Chapters 4, 5 and 7 present several fast imaging methods: In Chapter 4, Simplified S P E E D is proposed to accelerate M R A based on the sparseness of M R A data. In Chapter 5, the original S P E E D is reviewed, and its extension to 3 D and its combination with HalfFourier imaging have been explored. In Chapter 6, a new parallel imaging method named S P E E D - A C E is proposed by taking advantages of multiple coils available on modern M R scanners. Chapters 4, 5 and 6 focus on new fast imaging in M R I . As a continuation of this work, Chapter 7 studies noise behaviors and proposes a novel S N R improvement approach for the proposed fast M R I methods. 3 Finally, Chapter 8 briefly Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging reviews this thesis and proposes future work exploring further the potentials of phase information in M R I . 1.2 A Brief history of N M R and MRI In 1946, the phenomenon of N M R was first achieved in condensed matter by two independent groups of physicists using two different methods [6-8]; Bloch, Hansen, and Packard from Stanford University used the induction method [6,7], while Purcell, Torrey, and Pound from Harvard University used the absorption method [8]. The 1952 Nobel Prize in Physics was awarded to B l o c h and Purcell for their pioneering contributions to NMR. In the original experiment conducted by Bloch's group, a sample — water — was placed with two orthogonal coils in a strong static magnetic filed. After alternating electric current was introduced in one of the coils, an induction signal was detected in the other when the frequency of alternating current reached a specific value. Similarly, in the original experiment conducted by Purcell's group, a sample — paraffin — was wrapped by a coil and was placed in a strong static magnetic field. After alternating current was introduced i n the coil, significant power loss in the coil was detected when the frequency of current reached a specific value. Both the groups found a specific frequency denoted as OJo and currently known as the resonant frequency. The value of OJO was found to be 4 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging proportional to the strength of the static magnetic field used in the experiment, and can be described by the Larmor equation (o =yB 0 (1.1) 0 where yis the gyro-magnetic ratio depending only on the nucleus, and BQ is the strength of the static magnetic field. The gyro-magnetic ratio of a proton equals 26752 rad/sec/gauss, corresponding to a resonant frequency coo = 4 . 0 1 x l 0 rad/sec or 63.9 M H z 8 at a magnetic field BQ = 1.5 Tesla. N M R can analyze various proprieties of matter in a non-invasive way because its signal contains important information regarding the characteristic of the matter. However, it had been applied only to small and homogeneous samples until the first 2 D images were obtained with N M R by P . C . Lauterbur from the State University of N e w Y o r k at Stony Brook in 1973 [9]. Lauterbur applied a magnetic field gradient to an object — two tubes of water — and thereby made the Larmor frequencies of nuclei vary with spatial locations, encoding spatial information into N M R signal. The N M R signal acquired i n the presence of a magnetic field gradient represents a I D projection of the nuclei density distribution along the gradient direction. Applying the gradient in multiple directions, Lauterbur obtained a series of projections; he then reconstructed an image of the nuclei density distribution with the same reconstruction algorithm as used i n X - r a y computed tomography (CT). Almost at the same time, P. Mansfield from the University of 5 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging Nottingham demonstrated N M R ' s potential in mapping the 3 D structure of crystalline solids [10]. Lauterbur and Mansfield shared the 2003 Nobel Prize in Physiology or Medicine for their contributions to the development of M R I . A l o n g with the fast development of computer technology, N M R imaging has been commercialized and been made available in clinical practice since the early 1980's. Becoming increasingly popular, N M R imaging is now referred to as M R I without the somewhat scary word "nuclear". Today M R I is one of the most commonly used medical imaging modalities because of its non-invasive nature and powerful features [11]. 1.3 Basic Physics of N M R 1.3.1 Nuclear Magnetization The phenomenon of N M R is based on the fact that the nuclei are charged particles with a property known as "spin". Spin is purely a quantum mechanical construct loosely analogous to classical angular momentum. The nuclei of Hydrogen — also known as protons — are of particular importance in clinical M R I because Hydrogen (mostly in the form of water) is the most abundant atom in human body. Therefore, my discussion in this work focuses on protons. 6 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging According to quantum mechanics, a proton is a fermion with a spin of 1/2, and thus has 2x1/2+1 = 2 states, often referred to as "spin up" and "spin down". When an external static magnetic field BQ is applied, "spin up" protons preferentially line up parallel to the field, while "spin down" protons line up anti-parallel to the field. The interaction between the magnetic moment ju of each proton and the magnetic field Bo results in a potential energy E. E = -/lB (1.2) 0 In this thesis, all vector quantities are indicated with bold-faced characters. B y choosing the z direction to be along BQ and considering jU = y j , z z and j z = hs, the potential energy is given by E = -ju B =-yj B =-yhsB z 0 z 0 (1.3) 0 where y is the gyro-magnetic ratio, h is the Planck constant divided by 27t ( 1 . 0 5 x l 0 - 3 4 Joule-sec), and s is the quantum number for the angular momentum j and can only be z +1/2 or -1/2 for a proton [11]. Therefore, there are two possible energy levels for a proton in a magnetic field, corresponding respectively to the two possible quantum states: spin up and spin down. The potential energy of a spin up proton is (1.4) 7 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging Likewise, the potential energy of a spin down proton is E - ^ B 2 ° (1.5) This energy level splitting is the well known as nuclear Zeeman effect and is illustrated in Figure 1.1. Spin Down/ Anti-Parallel Spin U p / Parallel Figure 1.1: Zeeman splitting increases linearly with the magnetic field Bo- Spin up protons have lower potential energy than spin down protons; the former are aligned parallel with the magnetic field, while the latter are anti-parallel. The energy difference between the two states is AE = hco [11]. 0 A s discussed previously, the behavior of one single proton in a magnetic field is accurately described by using Eqs.(1.2-1.5). However, when a system with a large number of protons is considered, statistical mechanics should be used to describe the magnetization of the system. According to the Boltzmann distribution, the numbers of the protons in the parallel and anti-parallel states are given respectively by [11] 8 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging JV =Aexp(-^) kT (1.6) AL=Aexp(-|p) (1.7) + and where A is a constant of the system, T is absolute temperature, k is the Boltzmann constant (1.38x10 Joule/K), and N and is+are respectively the numbers of protons ± and their associated energies in the two states [11]. The population ratio between the parallel and anti-parallel protons can be easily estimated by using E q . (1.6) and E q . (1.7). N A C X P kP ( + N ~ Aexp -^) A£ k T yhB n k T ( Assume each unit volume contains N = N-+ N+ protons, each proton with a magnetic moment of ju. The magnitude of the total magnetic moment in the unit volume can therefore be calculated as [11] 9 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging M =H(N -N_) 0 + (N -N_) + = ju (N +N_) + (N +N_) + = uN-^—L (1.9) R+l exp(—) -1 exp(-) = LiN tanhf l + AF. 2kT ) In most systems, the factor of kT is usually much greater than the value of AE. For example, at a temperature of T = 293 K and a magnetic field of 1.5 Tesla, the value of AE/kT w i l l be as small as 1 . 0 4 x l 0 - 5 [11]. Therefore, E q . (1.9) can be Taylor expanded and simplified into M Q =vNtanh(—)~p:N— 2kT ^ V = Nju ^ 2 2kT = Ny h -^ 2 kT 2 4kT (1.10) A s seen from E q . (1.10), the magnitude of the magnetization Mo is proportional to the strength of the magnetic field Bo. The magnetization M o is called "the magnetization at thermal equilibrium value". E q . (1.10) also illustrates that M o i s inherently weak because it is formed by only a very small fraction of the total protons [11]. Although E q . (1.10) quantifies the amplitude of the magnetization, its dynamic behavior is yet to be 10 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance investigated. Imaging Therefore, in the following section, the motion of the magnetization is discussed in both quantum mechanical and classical sense. 1.3.2 M o t i o n o f the M a g n e t i z a t i o n : Q u a n t u m to C l a s s i c a l According to quantum mechanics, the expectation value of an observable F has the following equation of motion — <F>=-<H,F>=-<HF-FH> dt h h (1.11) where H is the energy operator known as the Hamiltonian. Considering F = fl and H = - / / • Z? , we have [11] 0 d i dt h — <//>=-<-// B ,/i> 0 = | <(-//• B )H - //(-// -B )> h 0 = h 0 (1.12) j<B x(fixju)> Q where fj. = yj is a product of gyro-magnetic ratio yaad angular momentum j . According to commutation properties of angular momentum operators i n quantum mechanics, j x j = ih j , the motion of a proton can therefore be described as [11] 11 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging d 7 — </i>=-y <B xj>=y<fi>xB dt (1.13) z 0 0 Consider a system with a magnetization M. The equation of motion for the magnetization vector M is therefore —M dt - yMxB (1.14) 0 Interestingly, E q . (1.14) can also be derived by using only classical physics. In a classical sense, the torque T exerted on M in the magnetic field BQ is T = MxB (1.15) 0 M and the time derivative of the corresponding angular momentum / = — is 7 — =T dt (1.16) Combining Eqs. (1.15) and (1.16), also considering the definition of % we have —M dt = yMxB (1.17) 0 which is exactly the same as E q . (1.14) and is known as the Bloch equation, which describes the behavior of the magnetization M in an external static magnetic field Bo [11] 12 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging as illustrated in Figure 1.2. The effects of relaxation are yet to be considered and w i l l be discussed in the section of Relaxation Times. Figure 1.2: The motion of M is a precession about the static magnetic field BQ, just like the precession of a spinning top in a gravitational field. 1.3.3 1.3.3.1 R F E x c i t a t i o n and N M R S i g n a l D e t e c t i o n Rotating Frame of Reference Because of Larmor precession and relaxation, it is difficult to describe the behavior of magnetization in a typical fixed reference frame; however, the difficulty can be easily overcome with a rotating frame of reference. In a frame of reference rotating at the Larmor frequency, magnetization is simply fixed in the frame i f the effects of relaxation 13 Chapter 1 Goals of the Thesis and Fundamentals are ignored. of Magnetic Resonance Imaging This is analogous to why we always choose a frame of reference rotating with the Earth rather than a frame of reference fixed on the Sun to simplify the motion on the Earth [11,12]. 1.3.3.2 R F Excitation A s discussed above, the magnetization vector M is constant in a frame of reference rotating at the Larmor frequency without taking relaxation into account. Therefore, the magnetization vector M can be tipped away from its initial position by using an orthogonal oscillating magnetic field tuned at the Larmor frequency. This process is referred to as " R F excitation" because the magnetic field applied usually rotates in the radio-frequency (RF) range [11,12]. A s illustrated in Figure 1.3, i f the magnetization M is initially at its thermal equilibrium Mo, a flip angle or of the magnetization M can be achieved by applying a rotating field B\ for a period of ras follows [11,12]: a=yB T (1.18) x 14 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging M(t) \^>y Figure 1.3: The magnetization vector M precesses about the effective magnetic field B\ in the rotating frame [11,12]. 1.2.3.3 Induction Signal Detection A s shown in Figure 1.3, there emerges a transverse component of the magnetization M after an R F excitation is applied. Because the precession of the transverse component causes magnetic flux change, a receiver coil placed nearby can detect an induction signal. Even after the R F field is turned off, induction signals can still be obtained because the transverse component of the magnetization M remains precessing until it decays out because of relaxation effects. The signal detected in this case is named "free-inductiondecay (FID)" [11-13]. 15 Chapter 1 Goals of the Thesis and Fundamentals 1.3.4 of Magnetic Resonance Imaging Relaxation Times According to the Bloch equation and Figure 1.3, once the magnetization M is tipped away by an R F field B\ from the main magnetic field Bo, M seems to be able to precess about BQ forever; however, this w i l l not happen because of relaxation effects. In reality, M w i l l eventually return to its initial thermal equilibrium position Mo along the direction of the magnetic field BQ. There are two major relaxation effects significantly influencing the behavior of the magnetization M. One is spin-lattice relaxation characterized by T l ; the other is spin- spin relaxation characterized by T2 [11-13]. The two relaxation effects are briefly described in the following section. A s shown in Figure 1.3, the magnetization M, having been tipped away from the main magnetic field Bo, can be treated as a combination of a longitudinal component along Bo and a transverse component orthogonal to Bo. After an R F excitation, the longitudinal component of M w i l l return to its thermal equilibrium value M o . The transition to the equilibrium state is caused by a transfer of energy from the nuclei to the nuclei's environment (lattice). Therefore, the relaxation along the longitudinal direction is named "spin-lattice relaxation", which is characterized with T l [13]. O n the other hand, the transverse component of M, after an R F excitation, w i l l dephase through interactions between different nuclei. A s a result, the net transverse magnetization appears to decay. Therefore, the transverse relaxation is referred to as "spin-spin relaxation", which is 16 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging characterized by T 2 [13]. Typically, different biological tissues have different T I and T 2 values. A s discussed previously, the Bloch equation depicted by E q . (1.14) and E q . (1.17) does not take the effects of relaxation into account. Including T I and T 2 into the Bloch equation gives a more comprehensive motion equation for the magnetization M, which is usually known as "the phenomenological B l o c h equation" [11,12]. dM dt where x, = MxB -—^ r (1.19) 0 2 1 y , and z are unit vectors along x, y and z directions, M , M , x y and M are the z components of the magnetization vector M along x, y and z directions, and BQ is the main magnetic field applied along z direction. 17 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging Solving E q . (1.19) along x, y and z directions respectively, we have [11] M (t) = [M (0) cos OJ t + M x x M (0 = [-M M (t) z 0 x (0) sin a^r] e x P(~ ^ ) ( (0) sin o^t + M (0) cos a^t] exp(——) T2 = M +[M 0 y - ) 2 0 (1.21) (0) - M ] e x p ( - ^ - ) z 1 (1.22) 0 where M (0), M (0), and M (0) are the initial components of the magnetization vector M x y z along x, y and z directions respectively, and 0)Q - }B is the Larmor frequency. 0 If we define transverse magnetization M (t) T magnetization M (t) L = M (t), z = ^M (t) 2 x +M (t) 2 y we have [11] M (0 = M (0)exp(-^) r M L (1.23) r (0 = M 0 and longitudinal + [M (0) - M L 0 ] exp(-^-) (1.24) where MT<0) and M (0) are respectively transverse and longitudinal components of M at L t = 0 . Eqs. (1.23) and (1.24) can be intuitively interpreted in Figure 1.4. 18 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance 'Mj(t) Imaging M (t) l L Mo Mj{0) M (0)\ L (a) (b) Figure 1.4: Time evolution of the magnetization M : (a) is the T 2 decay of transverse component o f M according to E q . (1.23); (b) is the recovery of longitudinal magnetization according to E q . (1.24) [11,13]. 1.4 Basic Principles of Magnetic Resonance Imaging A s discussed previously, a key step i n producing an M R image is to encode spatial information into the N M R signal. This section briefly reviews various techniques for spatial information encoding. The review firstly introduces the magnetic field gradient and describes the techniques of slice selection, frequency encoding, and phase encoding. Then several typical M R I pulse sequences are reviewed. Finally, the concept o f k-space is briefly discussed and is demonstrated with a fast imaging technique named E P I . 19 Chapter 1 Goals of the Thesis and Fundamentals 1.4.1 of Magnetic Resonance Imaging Magnetic Field Gradient According to E q . (1.1) known as the Larmor equation, the Larmor frequency coo varies with Bo- If a linear magnetic field gradient is applied along the z direction over the main magnetic field BQ, protons at different z locations w i l l have different Larmor frequencies. Thus, the z directional spatial information is successfully encoded frequencies. into Larmor Furthermore, this principle can be applied along an arbitrary direction so that the Larmor equation can be modified into [11,13] aXr) = y(B +Gr) (1.25) 0 where aXr) is the Larmor frequency of a proton at the location of r, G is the linear magnetic gradient field imposed. This equation indicates that the Larmor frequency OJ(r) is a function of the main magnetic field Bo, the magnetic field gradient G, and the location r [11,13]. 1.4.2 Slice Selection If a magnetic field Bo is superposed with a linear magnetic field gradient G, an R F excitation with a frequency of a; can selectively excite a slice that satisfies y(B +Gr) 0 = co (1.26) 20 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging The excited slice must be located in a plane perpendicular to the direction of G [11-13] because of its linearity. For example, selectively exciting a slice perpendicular to the z direction can be achieved by using a combination of an R F excitation and a linear magnetic field gradient along z direction. Figure 1.5 shows a typical sequence for the slice selection with a 90° R F pulse. j A 90° : RF A ! 1. If - i Slice selection I i \ ___x___j gradient G z \ J \~~~7 Figure 1.5: A typical pulse sequence of slice selection using a 90° R F pulse and a z directional linear magnetic field gradient G . z The negative gradient lobe applied after the R F pulse is used to refocus a linear phase dispersion caused by different Larmor frequencies across the slice thickness given by [11,13] (1.27) where d is the slice thickness and Aa) is the frequency bandwidth of the R F pulse. 21 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging According to the principle of slice selection, M R I can be used to image a slice in an arbitrary orientation, and the location of the slice can be shifted simply by offsetting the R F frequency without manipulating hardware, which is usually required by other imaging modalities such as X-ray C T . 1.4.3 Frequency Encoding The strategy of frequency encoding in M R I was first recognized and realized by P . C . Lauterbur from the State University of N e w Y o r k at Stony Brook in 1973 [9]. Lauterbur used a linear magnetic field gradient to distinguish the protons at different locations. His idea is based on the fact that protons at different locations precess with different Larmor frequencies in a linear magnetic gradient. The process of encoding spatial information into frequencies with a linear magnetic gradient is often referred to as "frequency encoding ( F E ) " [11-13]. A s shown in Figure 1.6, the frequency encoded N M R signal is collected in the presence of a gradient, which is often called the "readout gradient" or "frequency encoding gradient". 22 Chapter 1 Goals of the Thesis and Fundamentals A Slice selection j gradient G ' Frequency encoding gradient G \ j z \ / " J f " ! \~ : x ! I • - - | ; ' i j T ' n i Signal Imaging i 90° ii! A j K p RF of Magnetic Resonance I! I ! jl i 1 |( : — A • Figure 1.6: A pulse sequence of frequency encoding using a F E gradient G after a typical slice selection. x 1.4.4 Phase E n c o d i n g In addition to slice selection and frequency encoding, another strategy called phase encoding can be used for spatial information encoding. A s shown i n Figure 1.6, protons with the same x locations but different y locations have the same Larmor frequency, while protons with the same y locations but different x locations have different frequencies because of the application of G . If a linear magnetic gradient G is applied x y along the y direction for a short period between slice selection and signal acquisition, the protons with the same x locations but different y locations w i l l have the same Larmor frequency but different phase, depending on their y locations [11-13]. In this way, spatial information of protons can be successfully encoded by using a phase encoding gradient. 23 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging Such an encoding process is therefore named "phase encoding (PE)" [11-13]. Figure 1.7 shows a pulse sequence that uses slice selection, P E and F E to spatially localize protons in three dimensions. /h RF 90' _ Slice selection gradient G z Phase encoding gradient G v Frequency encoding gradient G x Signal V Figure 1.7: A pulse sequence containing a slice selection, P E with a linear magnetic field gradient G , and F E with G . y 1.4.5 x The Concept of £-Space A s discussed previously, after a gradient field G is applied, the Larmor frequencies of protons in a fixed reference frame can be described by E q . (1.25). In the rotating frame of reference, E q . (1.25) can be simplified into u)(r) = yG (1.28) r 24 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging If an arbitrary field gradient G(f) is applied to a sample represented by the magnetization M ( r ) , the phase of M(r) at a time t after an R F excitation can be expressed as <p(r,t) = \co(r,t)dt = (y\G(t)dt)-r o o = k(t)-r (1.29) t where we have defined k(t) - y\G(t)dt, o which is the position vector in k-space [11-14]. The time derivative of k(i) can therefore be calculated as follows: dk(t) dt 7G(t) (1.30) Clearly, G(i) is nothing but an instant velocity with a scaling factor of y in k-space. Since the induction signal is the total transverse magnetization in the rotating frame of reference, the induction signal can be expressed as 5 ( 0 = jM(r)exp[i0(r,t)]dr (1.31) B y combining Eqs. (1.29) and (1.31), we have S[k(t)] = SM(r)exp[ik(t)r]dr (1.32) 25 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging A s shown in the above equation, k-space is the Fourier counterpart of image space and as such is also named "spatial frequency space". It follows that reconstructing an M R I image from k-space raw data is just an operation of the Fourier transform. 1.4.6 P u l s e Sequences The previous section reviews magnetic field gradient, slice selection, frequency encoding, phase encoding, and the concept of k-space. The following section w i l l review briefly several typical pulse sequences and analyze the pulse sequence of E P I in k-space. 1.3.5.1 Spin Echo The pulse sequence of spin-echo is a basic M R I sequence, which was first proposed by Hahn in 1950 [15]. In his original experiment, a spin-echo was formed by using two 90° R F pulses. For simplicity, this work describes a standard spin-echo pulse sequence, which employs a 90° R F excitation pulse and a single 180° refocusing pulse to form a spin-echo. A s shown in Figure 1.8, the 90° R F pulse moves the magnetization from the longitudinal direction to the transverse plane. This is followed by an 180° refocusing pulse that reverses the dephasing of the transverse magnetization. 26 Chapter I Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging TR 180" 90" RF A TE Slice selection gradient G z Phase encoding gradient G v Frequency encoding gradient G s 4fa— Signal Figure 1.8: A typical spin-echo pulse sequence of M R I . In Figure 1.8, T R stands for "repetition time (TR)", which is the separation of two neighboring 90° R F pulses. T E is referred to as "echo time (TE)", which is the time separation between the 90° R F pulse and the peak of the echo signal [11-13]. Since T R and T E are both adjustable parameters of M R I pulse sequences, they can be adjusted to obtain various types of spin-echo images useful for clinical diagnosis, which include proton density images, T I weighted images, and T 2 weighted images [11-13]. 27 Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging 1.3.5.2 Inversion Recovery Inversion Recovery (IR) imaging is one of the most useful techniques i n M R I contrast manipulation [11-13,16,17]. A s shown in Figure 1.9, IR can be accomplished by applying an 180° R F pulse at a time of TI prior to the 90° R F pulse [11-13]. The 180° R F pulse is called "inversion pulse", and TI is referred to as "inversion time (TI)" [11-13]. Because longitudinal magnetization (M ) ranges from -Mo to Mo, various T l contrast can z be produced by IR imaging, leading to its two major applications: T l selective signal suppression by setting TI = T l l n 2 = 0.69T1, and T l contrast enhancement with doubled dynamic range (-Mo, Mo) of M , which is twice as much as that of the ordinary saturation z recovery experiment (0, Mo). For example, the IR sequence can be used to suppress strong signals from fatty tissue with a known T l value for the purpose of obtaining an improved image contrast [17]. I /!' 180° JV ; A RF \ / i \J ,< x j • Signal — —j— ; ; ' v +j i. TI : : I'll A.... A : 0 — - — : — [, j I' ^. i Figure 1.9: Inversion recovery pulse sequence. 28 — Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging 1.3.5.3 Gradient Echo A gradient-echo, as its name implies, is produced by switching gradients to refocus the dephasing transverse magnetization, and thus is also called "gradient-recalled-echo" [1113]. Unlike the spin-echo pulse sequence, the gradient-echo pulse sequence does not use a 180° R F refocusing pulse, but applies a series of gradients to turn the dephasing transverse magnetization into an in-phase state. Figure 1.10 shows a typical 3 D gradientecho sequence that uses two P E in both y and z directions. i i A RF —r\ 90° I i!' f Secondary phase encoding gradient G 2 Primary phase encoding gradient G — Y Frequency encoding gradient G — x I \ - J Signal • i : -M Figure 1.10: A typical 3 D gradient-echo pulse sequence. 29 IAT Chapter 1 Goals of the Thesis and Fundamentals of Magnetic Resonance Imaging 1.3.5.4 Echo Planar Imaging Echo planar imaging (EPI) is a fast M R I strategy, which obtains an entire 2 D k-space data after a single R F excitation, typically in 50 ms. E P I was originally proposed by P. Mansfield in 1977 [18]. N o w it is still one of most difficult pulse sequences to implement because of its great burden on hardware as well as biological considerations. The original version of E P I pulse sequence is shown in Figure 1.10 (a), associated with its k-space trajectory in (b) [11]. RF G, _ .r:.:::..—---^- G, Signal ••• 11 (a) (b) Figure 1.11: The original E P I pulse sequence (a) and its k-space trajectory (b) [11]. 30 c Chapter 1 Goals of the Thesis and Fundamentals 1.5 of Magnetic Resonance Imaging Summary In this chapter, the phenomenon of N M R is briefly described and the basics of M R I are selectively reviewed. A s the origin of N M R signal, magnetization and its behavior are described not only in a quantum mechanical but also in a classical sense. To quantitatively describe the motion of magnetization, the Bloch equation is introduced and its comprehensive version involving relaxation effects is also given in this chapter. To illustrate the principles of M R I , the concept of magnetic field gradient is introduced and its three applications in spatial localization are reviewed: slice selection, frequency encoding, and phase encoding. reviewed. After this, several typical pulse sequences are briefly Finally, the concept of k-space is introduced and is demonstrated with the pulse sequence of E P I . 31 Chapter 2 Review of Applications of MRI Phase Information Chapter 2 Review of Applications of MRI Phase Information The purpose of this chapter is to introduce phase information in M R I and to provide a brief review of its applications, revealing the motivation of the thesis. briefly describes M R I phase and the problem of phase wrapping. This chapter Then it selectively reviews the applications of phase information in the areas of flow imaging, chemical shift imaging, and fast imaging in M R I . Inspired by these applications, I explore further the potentials of M R I phase information i n inversion recovery imaging and fast imaging. 2.1 MRI Phase Information and Phase Wrapping A s discussed in the previous chapter, M R I maps the spatial distribution of the transverse magnetization in the rotating frame of reference. The M R I signal is complex valued [11- 13] and contains information about the amplitude and direction of the transverse magnetization. The directional information characterizes the rotation of the transverse magnetization and is hereafter referred to as " M R I phase information". M R I phase not only can be used to encode spatial information through P E , but also often conveys very important information such as field inhomogeneity, velocity of blood flow and chemical shift [19]. Usually, M R I phase information is represented by a phase map, while the amplitude of the transverse magnetization is described by a magnitude image. 32 A s an Chapter 2 Review of Applications of MRI Phase Information example, Figures 2.1 (a) and (b) show respectively the magnitude and phase images from a transverse gradient-echo head scan. Figure 2.1: A n example obtained from a transverse gradient-echo head scan, (a) is the magnitude image representing the amplitude of the transverse magnetization, and (b) is the corresponding phase map that shows the direction of the transverse magnetization. Although a lot of important information can be encoded into M R I phase, this is often counterbalanced by the fact that the value of a phase map can only be derived as modulo 2jt, which gives a principal value within the range of [-71, 7t) [20]. This implies that the M R I phase map may not necessarily represent the true phase information or the true rotational information of the transverse magnetization. For example, rotations of 7t/4 and -luJA are not distinguishable in a M R I phase map. This is because -77t/4 is beyond the principal value range [-7t, 7t), and is hence wrapped around to yield a principal value of 7t/4. Such phase wrapping is caused by adding or subtracting true phase by a multiple of 27C to yield a principal value in [-7E, 7t) [20]. 33 Chapter 2 Review of Applications Because of MRI Phase of phase wrapping, the Information information encoded in M R I phase may be misinterpreted, which limits its applications in clinical diagnosis. However, M R I phase information should not be undervalued, because neglecting M R I phase means losing the important information encoded only in M R I phase [19,20]. Therefore, the applications of phase information in M R I have been explored by many researchers. 2.2 Review of Applications of MRI Phase Information This section briefly reviews various techniques regarding M R I phase information encoding. This review begins with flow imaging, which is followed by chemical shift imaging. Then, fast imaging in M R I is explored. 2.2.1 Flow Imaging There are various types of physiological flows in a human body including the flows of blood, lymph and cerebrospinal fluid. Such physiological flows result in the motion of nuclei that can alter both magnitude and phase of M R signal [21-23]. Therefore, both magnitude and phase can be made sensitive to flows and thereby be used for flow imaging [11-13]. Since M R I phase information encoding is the focus of this thesis, magnitude change caused by flow w i l l not be discussed further. 34 Chapter 2 Review of Applications 2.2 A A of MRI Phase Information Velocity Sensitization in a Bipolar Gradient Field Consider a velocity sensitized bipolar gradient shown in Figure 2.2. First, the gradient with a negative value of -G is turned on at t = to for a time duration of T. After a time period of A , the gradient with a positive value of G is turned on at t = to +T+A for the same time period of T [11,12]. i Gradient G(t) 0 Jo Time t •G T Figure 2.2: A bipolar gradient used for velocity phase sensitization [11]. This bipolar gradient as a whole does not affect stationary nuclei, because the effect on the stationary nuclei caused by the first negative lobe is completely compensated by the subsequent positive one. In other words, the net effect on the stationary nuclei by the bipolar gradient is zero. However, it is interesting to see that the bipolar gradient causes non-zero net effect on the nuclei moving along the direction of the gradient. This is because the Larmor frequencies of the moving nuclei vary with time, which is a result of 35 Chapter 2 Review of Applications of MRI Phase Information their temporally varying locations in the gradient field. The net effect on the moving nuclei is usually interpreted as a phase shift, or a direction change of the transverse magnetization in the rotating frame of reference [11]. Therefore, the velocity of flow can be sensitized with a bipolar gradient field, and thereby is encoded into M R I phase. This is the basic principle of the phase encoding for velocity sensitization, or the so-called "phase-contrast ( P C ) " technique for flow imaging [11-13,22,24], which is demonstrated below with a specific example. Consider a tissue element moving with a constant velocity v along the direction of the bipolar gradient. After the bipolar gradient shown in Figure 2.2 is applied, the phase shift can be expressed as 0= rl' r(OG(t')dt'= to rl' (r + vt')G(t')df t0 0 = r^ \r +vO(-G)dt'+y!%^ (2.1) + 0 = yvGT(T + A) Being proportional to the velocity v and gradient strength G, the phase shift <f> does not depend on the initial location ro of the tissue element and the starting time to of the bipolar gradient [11]. 36 Chapter 2 Review of Applications 2.2.1.2 of MRI Phase Information Phase-Contrast M R Angiography A s described above, a bipolar gradient causes a phase shift i f nuclei are moving along the direction of the gradient. Based on this fact, the motion of nuclei or flow can be mapped with two interleaved complex images. One of them is a velocity-sensitized image acquired with a bipolar gradient, while the other is a reference image acquired with an identical sequence but without the bipolar gradient. These two images have different phases at the pixels where the nuclei are moving along the direction of the gradient. Therefore, flow velocities can be mapped simply by either subtracting the two complex images or calculating their phase differences. Such a method is used for M R Angiography ( M R A ) and is known as phase-contrast (PC) M R A [11-13]. 2.2.2 Chemical Shift Imaging A second typical technique using M R I phase information is chemical shift imaging (CSI). C S I is the production of spatially localized spectra and thus is considered a type of M R spectroscopy [25]. In C S I , spectroscopic data are obtained simultaneously and then are decoded into a set of spatially localized spectra. They can be presented either as a set of spectra representing the chemical composition of each pixel, or as a set of images mapping the spatial distribution of each chemical component. 37 Chapter 2 Review of Applications 2.2.2A of MRI Phase Information Chemical Shift A s the name of C S I suggests, it is based on the concept of chemical shift. Chemical shift characterizes the variation in the Larmor frequency of a given nucleus in different chemical environments [25]. The amount of the shift is proportional to magnetic field strength and is usually specified in parts per million (ppm) of the resonance frequency relative to a standard. It is the chemical shift that makes it possible to distinguish different chemicals in M R I . If we treat the chemical shift as an additional dimension, chemical shift images can be obtained by using Fourier transform (FT) after chemical shift information is frequency encoded [26-28] or phase encoded [29]. For example, signal can be acquired in the absence of gradient once spatial information is encoded with P E . Thus, chemical shift images can be reconstructed from the sampled data simply by using the F T . However, a lot of time is required for data acquisition in these C S I methods because of additional chemical shift dimension [30]. Amongst the most commonly used C S I techniques is water-fat imaging, which maps spatial distribution of water or fat i n a human body. It is well known that chemical shift between water and fat is 3.5 ppm [12]. This corresponds to 220 H z at a magnetic field of 1.5 T, making water and fat protons switch in and out of phase every 2.25 ms after an R F pulse is applied. 38 Chapter 2 Review of Applications 2.2.2.2 of MRI Phase Information Brief Review of Water-Fat Imaging W i t h the known chemical shift between water and fat, the proton spectra can be modeled as two peaks with a separation of 3.5 ppm [31]. Based on this double-peak model, many efficient approaches have been developed to achieve water-fat imaging with resolution, signal-to-noise ratio ( S N R ) and data acquisition time, which are comparable to conventional M R I [12, 30]. These approaches can be divided into two categories. Those in the first category usually apply a presaturation R F pulse to suppress fat (or water), allowing the other chemical component to be imaged [12]. The success of these approaches relies on a very high (< 3.5 ppm) homogeneity of the magnetic field, which is extremely difficult to obtain in reality [12,30]. The approaches in the second category are based on the idea that the chemical shift information of water and fat can be encoded into M R I phase information. This category includes the simple proton spectroscopic imaging known as the " D i x o n method" [31], the method with Quadrature-Sampling (QS) [32], and the extensions of D i x o n method [30, 33]. In this work, our discussion focuses on the second category. The following sections firstly discuss water-fat encoding with M R I phase by using a spin-echo sequence and a gradient-echo sequence. The D i x o n method [31] and the Q S method are then reviewed in a general perspective. This is followed by a brief description of a three-point direct phase encoding ( D P E ) method [30] for water-fat imaging, an extension of the D i x o n method. 39 Finally, chemical shift imaging with Chapter 2 Review of Applications of MRI Phase Information spectrum modeling ( C S I S M ) [32] is further explored as a generalization of the threepoint D P E method. 1) Water and Fat encoded by using M R I Phase Information Based on the water-fat double-peak model, each pixel of a typical M R image can be considered as a summation of two magnetization vectors representing respectively protons associated with water and fat. rotate at different angular Because of the chemical shift, the two vectors frequencies, and thereby interfere constructively or destructively, depending on the timing parameters of the pulse sequence. If the timing parameters of the pulse sequence are deliberately manipulated to control the interference of the two vectors, the resulting M R image can be expressed as [12,30,31,34] I = [W + exp(/a)F]exp(/0)exp(/0) (2.2) where W and F are respectively the magnitudes of the magnetization vectors of water and fat, a is the phase angle of the two vectors caused by chemical shift, 6 is the phase error caused by magnetic field inhomogeneity, and 0 is a static phase error caused by R F field inhomogeneity and data sampling window off-centering, etc. 40 Chapter 2 Review of Applications of MRI Phase Information 180" A 90" RF Slice selection gradient G [ z Phase encoding gradient G y Frequency encoding gradient G x Signal _ i /2 to/2 - • ( • » At TE (a) A .90" in RF Slice selection gradient G : TIT" Phase encoding gradient G v Frequency encoding gradient G ^47 v Signal Ar « re (b) Figure 2.3: Chemical shift information of water and fat can be encoded with a spin-echo sequence (a) or a gradient-echo sequence (b) by manipulating the timing parameters: the refocusing time to, the echo time TE, and the time shift Af [30]. 41 Chapter 2 Review of Applications of MRI Phase Information The complex image / can be acquired by using either a spin-echo sequence shown in Figure 2.3 (a) or a gradient-echo sequence as shown in Figure 2.3 (b) [30]. In E q . (2.2), relaxation effects are ignored or grouped into W and F, and a and 6 are both proportional to the time shift At in the pulse sequences shown in Figure 2.3. The mathematical expression of oris a = AcoAt (2.3) where A^yis the difference in Larmor frequencies between water and fat ( e.g. 220 H z at 1.5 T), and similar for 6 6 = yAB At (2.4) Q where ABo represents the local offset of the main magnetic field BQ because of field inhomogeneity. Thus, the phase angle a between the two magnetization vectors of water and fat vectors is successfully encoded by using the time shift At, which is measured from the refocusing time to to the echo time TE as shown in Figure 2.3. 2) The D i x o n Method In the D i x o n method, the time shift Af is adjusted to put the magnetization vectors of water and fat in-phase or opposed-phase, so that a pair of images can be obtained respectively with a= 0 and n [31]. According to E q . (2.3) and Figure 2.3, at 1.5 T, the 42 Chapter 2 Review of Applications of MRI Phase Information phase angle of a = 0 can be accomplished by using Af = 0 in the pulse sequence, while a = 71 is obtained by using Af = 2.25 ms. In this way, the two complex images can be expressed as [12,30] / , =oy + F)exp(ty) I =(W- (2.5) F) exp(/0) exp(i0) 2 (2.6) In the original D i x o n method, simply adding and subtracting I\ and h yields W and F i f phase errors 0 a n d 0are ignored. Because phase errors are clearly not negligible in most cases, phase error 0 should be determined in order to solve for W and F. Given the fact that W and F are both real and non-negative, 0can be estimated as [12,30] 20 = a r g [ ( / / * ) ] (2.7) 2 2 1 where arg [] returns the phase of the complex number i n the square bracket, and I* is the complex conjugate of I\. According to E q . (2.7), 0 can be directly estimated i f 29 is within the principal value range [-TC, TC). This, together with the fact that the 3.5 ppm water-fat chemical shift corresponds to a phase angle of TC, implies that the field inhomogeneity must be less than 1.75 ppm over the F O V . M a k i n g 26 within [-TC, TC) is practically difficult because even the subject induced field variation can easily exceed 1.75 ppm caused by magnetic 43 Chapter 2 Review of Applications of MRI Phase susceptibility variation of tissues [12,30]. Information Therefore, solving 6 has to depend on the success of 2 D phase unwrapping, which remains a problem yet to be completely solved [20]. 3) The Quadrature-Sampling (QS) Method Shortly after D i x o n method was published, a related single-point QS method was proposed, which acquires a single image when the vectors of water and fat are placed orthogonal to each other [32]. In the absence of phase errors, the real part of the complex image is the water image, and the imaginary part is the fat image; however, the method of quadrature sampling is limited because of the presence of phase errors in reality. 4) Three-Point D P E method In 1997, X i a n g and A n proposed a three-point method using direct phase encoding to solve water and fat on a pixel-by-pixel basis [30]. A s an extension of D i x o n method, this method acquires three complex images with different time shifts: Af = in, Xo+T, and Xo+2x, corresponding to an asymmetric sampling scheme (a&, (Xo+a, (Xo+2a). The threepoint D P E method encodes chemical shift information into M R I phase and thereby can solve for separate images of water and fat in an analytical approach [30]. 44 Chapter 2 Review of Applications of MRI Phase Information 5) Chemical Shift Imaging with Spectrum Modeling ( C S I S M ) Although the 3-point D P E method is proposed for water-fat imaging, its principle can be generalized to handle three or more chemical components. A s an extension of the 3-point D P E method, C S I S M proposed by A n and X i a n g models the chemical shift spectrum "as several delta-function shaped peaks with known relative Larmor frequencies but unknown amplitude"[33]. In this way, C S I S M can analytically solve separate images of multiple chemical components with only a few D P E spin-echo images [33]. 2.2.3 Fast I m a g i n g i n M R I Fast imaging is another application of M R I phase information, and is very desirable in M R I for improving its performance as well as increasing patient comfort level. The clinical significance of reducing scan time has driven the development of many fast imaging strategies since M R I came into being. Some strategies physically manipulate spin dynamics to make better use of available magnetization [18, 35-38], while others sparsely sample k-space and reconstruct a complete image through a non-standard reconstruction [3,39-52]. Half-Fourier imaging [39] is a typical example using M R I phase information to reduce scan time. 45 Chapter 2 Review of Applications 2.23 A of MRI Phase Information Half-Fourier Imaging Half-Fourier imaging is a fast acquisition and reconstruction technique, which accelerates M R I by recovering a desired image from only slightly more than half of k-space data. In a mathematical sense, the recovery is equivalent to deconvolving the truncation artifacts from image space [39]. The theory of Half-Fourier imaging can be illustrated as follows: 1) Basic Principle of Half-Fourier Imaging If M R signals are sampled to cover half of k-space, the acquired data can be Fourier transformed (FT) into a complex image associated with truncation artifacts. It is well known that sampling a half k-space data set is equivalent to multiplying full k-space data with a unit step function in k-space. According to the convolution theorem [12], i f kspace data are partially sampled along x direction, the reconstructed image l(x,y) from the partial data can be written as I(x,y) = where S(x) = — \S(x)- C(x,y)®S(x) (2.8) is the Fourier transform of the unit step function [12], and C(x,y) is the 2 D complex image reconstructed from full k-space data. If the phase map of C(x,y) is known to be <fi(x,y), multiplying exp[-iflx,y)] gives a real-valued image R(x,y) as follows, 46 and C(x,y) Chapter 2 Review of Applications of MRI Phase Information R(x,y) = C(x,y)exp[-i0(x,y)] (2.9) The real-valued image R(x,y) is the desired image. In order to obtain it, the reconstructed image I(x,y) is multiplied by exp[-i<p(x,y)] and as such yields [39] /(*, y) exp[-i0(x, y)] = R(x, y) ®S(x) = 2 - i[R(x, y) ® - L ] 2^c (2.10) It can be shown that the real part given by E q . (2.10) can be simply doubled to yield the desired image R(x,y). This is the basic theory of Half-Fourier imaging in an ideal situation. 2) Practical Implementation of Half-Fourier Imaging In practice, the implementation of Half-Fourier imaging can be divided into the following four steps, demonstrated below with in vivo data from a spin-echo transverse head scan. (1) First, slightly more than half of k-space data are sampled and are then Fourier transformed into a complex image, which is associated with truncation artifacts caused by undersampling. For example, Figure 2.4 (a) shows the magnitude image reconstructed from full 256x256 k-space data, which are acquired with a spin-echo transverse head scan. T o reduce scan time, the 256x256 k-space matrix is sampled as follows: 64 k-space lines are consecutively sampled near k-space center, and then half of k-space excluding the 64 central lines is fully sampled. 47 Chapter 2 Review of Applications of MRI Phase Information (2) Second, a phase estimate is obtained with the consecutively sampled data near the kspace center. For the above example, the phase estimate is determined by using the central 64 k-space data lines. (3) Third, the undersampled k-space data set is multiplied by a Harming filter [39] to reduce truncation artifacts. (4) Finally, the filtered k-space data are Fourier transformed into image space, and phase error is then removed by multiplying the complex conjugate of the phase estimate. The real part of phase corrected complex image is obtained as the final image. Compared with the "gold standard" shown in Figure 2.4 (a), the final image shown i n Figure 2.4 (b) is reasonably good but has a S N R reduction of ~V2 [39]. (a) (b) Figure 2.4: Half-Fourier is demonstrated with in vivo data from a spin-echo transverse head scan, (a) is the magnitude image reconstructed from full k-space data, and is therefore treated as the "gold standard", (b) is the recovered and desirable image obtained from 62.5% data by using Half Fourier imaging. 48 Chapter 2 Review of Applications 2.3 of MRI Phase Information Objectives of the Thesis A s discussed in Chapter 1 and Chapter 2, M R I signal is complex valued and can be used to derive a magnitude image and a phase map. Although magnitude images are used much more frequently than phase maps for clinical diagnosis, phase maps should not be undervalued because they often encode inhomogeneity, motion, and chemical shift. important information such as field Because a lot of valuable information is encoded only in M R I phase, the phase information in M R I has great potential in medical applications. M y doctoral thesis focuses on the study of phase information in M R I and its applications to fast imaging and inversion recovery imaging. Since phase correction is desirable in inversion recovery imaging, I also propose an efficient method for phase correction, which is described in the following chapter. 49 Chapter 3 IR Imaging with Non-linear Phase Correction Chapter 3 IR Imaging with Non-linear Phase Correction based on an Extended Statistical Algorithm 3.1 Introduction A s discussed in Chapter 1, inversion recovery (IR) imaging is used to enhance image contrast among tissues by making use of their differences in T l relaxation times [1113,16,17]. IR imaging has two major applications [11-13]: signal suppression of tissues with known T l values, and T l contrast enhancement with doubled dynamic range (-Mo, Mo) of M . The former has been successfully used in selective tissue signal suppression z such as fat suppression [17]; however, the use of latter has been limited because it requires real image reconstruction instead of the conventional magnitude image reconstruction [12]. In a perfect world where phase errors caused by magnetic field inhomogeneity and other factors are negligible, the image reconstruction for IR imaging is a simple display of the real part of the complex image. In reality, however, phase errors are always present in M R images and therefore must be corrected to reconstruct an IR image. In the past few years, many phase unwrapping and correction methods have been developed, which generally fall into two categories: path-following methods and minimum-norm methods 50 Chapter 3 IR Imaging with Non-linear Phase Correction [53]. Cut-line methods [20,54-55], disk growing method [56], region growing methods (including maximum [57] and minimum spanning tree methods [58]) and Flynn's minimum discontinuity approach [59] are all path-following methods, which are based on the idea of guiding the integration of the gradient map to avoid any error propagation [20]. Minimum-norm methods are used to minimize cost function to achieve phase unwrapping [53]. T o minimize cost function, the well-known least squares algorithm [60] or a generalized Z/'-norm algorithm [61] with data-dependent weights is generally used, leading to a direct minimum-norm solution to phase. Alternatively, phase can be modeled with a function by estimating its coefficients in a weighted least squares ( W L S ) sense [53]. Amongst the model-based methods, polynomial fitting is the most commonly used in M R I because phase variation in M R I can often be modeled with a polynomial [1,62-66]. Although polynomials may not necessarily represent all phase maps, especially when poles [20] are present in phase, polynomial fitting is of great use in the case of disconnected tissues in the field-of-view ( F O V ) [67], a typical example of which is an image showing two separated legs. A successful example of polynomial fitting is a linear phase correction method using autocorrelation, which was proposed by A h n and Cho ( A C method) [1]. W h i l e the A C method is efficient, its straightforward extension to handle higher order non-linear phase terms fails to yield good results because the signal in the estimated phase terms is substantially contaminated by noise as a result of 1-pixelshift phase differentials. 51 Chapter 3 IR Imaging with Non-linear Phase Correction In this work, we introduce the n-pixel-shift rotational differential field (RDF), which represents local vector rotations of a complex field relative to itself after being shifted by n pixels [2]. The amount of pixel shift is crucial, because we found that an increase in the amount significantly enhances the signal without changing the noise, resulting in a more accurate phase correction. Based on the n-pixel-shift R D F , we have successfully extended the A C method to handle higher order non-linear phase error terms, which are often important for polynomial expansion of phase variation in M R I [2]. The n-pixelshift R D F can also improve the weighted least squares phase unwrapping method proposed by Liang ( W L S method) [1,66]. The improvements brought by the n-pixel-shift R D F have been demonstrated with 2 D in vivo IR M R I data. 3.2 Method 3.2.1 R e v i e w o f the A C Method The A C method [1] statistically fits a linear polynomial or a plane to phase variation in M R I . Given a noise-free M R image with a zeroth-order phase error, represented by a and a first-order linear phase error bx along x direction, we can write the complex image as S(x,y) = S (x,y)exp[i(a 0 + bx)] where So(x,y) is a real-valued ideal image without phase error. model parameters " a " and "b", and to achieve phase correction. 52 (3.1) The goal is to estimate Chapter 3 IR Imaging with Non-linear Phase Correction Since So(x,y) and S(x,y) are both idealized images that do not account for other sources of phase error ( e.g. chemical shift, etc.), parameter b can be extracted by using a complex field defined as A%(x,y) = S(x,y)S\x-\,y) (3.2) where S* denotes the complex conjugate of S. The complex field defined by E q . (3.2) represents local vector rotations of the complex image S(x,y) relative to itself after being shifted by 1 pixel. In this work, the complex field is called "first-order x-directional 1-pixel-shift rotational differential field ( R D F ) " of S(x,y), and is denoted as A ^ , , given that superscript denoting the order of R D F is generally ignored in this work when the order is 1. Substituting E q . (3.1) for S in E q . (3.2) gives A *,i y) o (*. y)So* ( ~ = S x L y) exp(&) (3.3) Since So(x,y) is real-valued, the quantity of the statistical expectation of So(x,y)So*(x-l,y) is always real and positive for most images [68]. Therefore, b can be determined as b = rg{E[A (x,y)]} a (3.4) xl where E represents the statistical expectation operation over the F O V . 53 Chapter 3 IR Imaging with Non-linear Phase Correction If noise is present, the noisy M R image can be expressed as C(x,y) = S(x,y) + N(x,y) (3.5) where the noisy signal C(x,y) is modeled as a superposition of a true signal field S(x,y) defined by E q . (3.1) and a complex Gaussian noise field N(x,y) [69]. Similarly, the parameter b can be accurately estimated as the phase of the statistical expectation of C(x,y)C*(x—l,y), which automatically emphasizes signal over noise through a squared magnitude weighting. Furthermore, it has recently been reported that the A C method has its Fourier counterpart with complementary operations offering high computing efficiency for certain applications [70]. The Fourier counterpart of the A C method can be derived by using the Fourier power theorem [71], which can be mathematically expressed as 2Zf(x,y)g*(x,y) where F(k ,k ) x y and G(k ,k ) x y = 2ZF(k ,k )G*(k ,k ) x y x y are respectively Fourier transforms (3.6) of image domain functions f(x,y) and g(x,y). If we define j\x,y) = C(x,y) and g(x,y) = C ( * - l , y ) and consider the Fourier shift theorem [71], we have 54 Chapter 3 IR Imaging with Non-linear Phase Correction I C(x, y)C*(x-l,y) = £ | F(k , x k ) | exp(i2/tk 2 y x IN) (3.7) Therefore, the linear coefficient b can be estimated as the angle of the centroid of an angular distribution of |F(^,^ )| exp(/2^: /A0 in k-space as follows: 2 y A b = arg{£[|F(k ,k )\ x 2 y exp(i2nk /N)]} x (3.8) where F(k ,k ) is the Fourier transform of the complex image C(x,y). x y After the linear coefficient b is determined, the estimated first-order linear phase error bx is removed from C(x,y) by multiplying its complex conjugate exp(-ibx). In the original A C method [1], the zeroth-order phase error a is estimated by using histogram after the estimated linear phase error bx is removed. In this work, we found a can generally be determined by using E q . (3.9) because So(x,y) is positive real in most cases. a = a r g { £ [ C ( x , y)exp(-ifa)]} (3.9) Although E q . (3.9) cannot be used directly for certain cases such as IR imaging, a combination of E q . (3.9) and additional operations may be used to handle these cases. For example, phase doubling can be used to extend E q . (3.9) to handle IR imaging cases, and w i l l be addressed later in this chapter. 55 Chapter 3 IR Imaging with Non-linear Phase Correction Although the A C method and its Fourier counterpart have been demonstrated to be robust and efficient in linear phase correction [1,70], their applications are limited because phase error appears to be non-linear in many M R images. To correct non-linear phase error, it is desirable to extend the A C method to higher order. straightforward extension of the A C method However, we found a for higher order non-linear terms unsuccessful, especially where the complex images are abundant in noise. This is due to poor signal performance as a result of phase differential operations in calculating the 1pixel-shift R D F by the A C method. Although its performance can be improved by suppressing noise with low-pass or smoothing filters, we have found that the application of such filters is still insufficient to overcome this drawback. 3.2.2 N - P i x e l - S h i f t R o t a t i o n a l Differential F i e l d To overcome the drawback of 1-pixel-shift R D F , we developed n-pixel-shift R D F , which represents local vector rotations of a complex field relative to itself after the field is shifted by n pixels [2]. The shift is conducted by simply shifting out the complex field by n-pixels along the desired direction without wrapping it around. The key difference between the n-pixel-shift R D F and the original 1-pixel-shift R D F is a more general, and usually larger pixel shift. W e have found that the amount of pixel shift is very crucial because a larger shift significantly enhances the signal without amplifying the noise in 56 Chapter 3 IR Imaging with Non-linear Phase Correction the estimated phase terms, and as such allows a much more accurate phase correction [2]. The reasons are illustrated in both image space and k-space below. In image space, the noisy complex image can be expressed as E q . (3.5). Its first-order xdirectional n-pixel-shift R D F is denoted as A xn and is found to be x,n (*» y) = ( > y)C*(x - n, y) A c x = S(x, y)S*(x-n, y) + S(x, y)N*(x-n, + N(x, y)S*(x -n,y) + N(x, y)N*(x y) (3.10) -n,y) Because of the random nature of the noise field N(x,y), when a total summation is performed over the entire F O V , the contribution from the last three terms should not depend on the choice of pixel shift " n " , unlike that from the first term, which has a larger phase angle when n is bigger. This means that the phase signal contained in the first term has been enhanced without any change in the noise in other terms, leading to an improved phase correction. This can be illustrated intuitively in the example shown in Figure 3.1. A s shown in Figures 3.1(c) and (d), when there is no noise, the linear polynomial coefficient can be accurately estimated from the angularly concentrated complex pixel distribution i f either 1-pixel-shift R D F or n-pixel-shift R D F is used. However, when noise is present, the sharp angular distributions shown in Figures 3.1(c) and (d) are blurred by the last three terms of E q . (3.10). Therefore, the estimation of the phase angle of the complex pixel distribution in Figure 3.1(g) becomes less reliable because of angular uncertainty. In contrast, the angle estimation in Figure 3.1(h) is much more 57 Chapter 3 IR Imaging with Non-linear Phase Correction accurate because the angular uncertainty remains the same, but the angle itself is 30 times greater than the angle shown in Figure 3.1(g), resulting in a much reduced relative error. Im Re (d) "" •"...- • (e) (f) (g) GO Figure 3.1: (a) is real part of a simulated complex M R image denoted as S(x,y) that is obtained by imposing a noise-free phase error 1.2 + 0.03 y rad onto a magnitude image, (b) is complex pixel distribution of S(x,y). (c) and (d) are complex pixel distributions of A y jand A y 3 0 o f S(x,y), respectively, (e) is noisy M R image denoted as C(x,y) that is the addition of S(x,y) and a complex Gaussian noise field [69]. (f-h) are repeated calculations of (b-d) for C(x,y). The angular estimation in (g) is less reliable because of relative large angular uncertainty, while the angular estimation in (h) is much more accurate because the relative angular uncertainty is greatly reduced by enhancing the angle by a factor of 30. This is the key idea of the proposed method. In k-space, the Fourier counterpart of the A C method depicted in E q . (3.8) shows that the linear coefficient b corresponds to the angle of the centroid of an angular distribution of a k-space factor |F(fcj,fc )| exp(/2;*yW). 2 v Indeed, in this case, the k-space positions k and x 58 Chapter 3 IR Imaging with Non-linear Phase Correction k +N/n are assigned the same phase value rather than unique phase values. However, this x phase wrapping generally appears in high-frequency k-space and as such does not affect the k-space center, which has the most energy. Therefore, the statistical coefficient estimation is mainly determined by the heavily weighted k-space center. Therefore, accuracy w i l l not be compromised because of the phase wrapping. In other words, the shift n does not change the noise distribution due to phase wrapping of noise, but changes the phase of the k-space factor from 2nkx/N to 2rnikx/N, increasing the angle of the centroid by a factor of n as long as it is not wrapped. Although only first-order n-pixel-shift R D F is demonstrated, it is straightforward to extend n-pixel-shift R D F to any higher order. For example, m^-order x-directional n- pixel-shift R D F can be calculated as follows: A ^(x,y) { where A^™~ \x, l = A^ \x,y)A ^\x-n,y) X { n x (3.11) y) is ( m - l ^ - o r d e r x-directional n-pixel-shift R D F . 3.2.3 E x t e n d e d A C M e t h o d Based on the n-pixel-shift R D F , we propose a new statistical phase correction method, which extends the A C method to include non-linear phase correction (extended A C 59 Chapter 3 IR Imaging with Non-linear Phase Correction method) [2]. The extended A C method works in a fashion similar to the original A C method except that the former uses repeated n-pixel-shift R D F calculation to extract higher order phase errors. Because of signal enhancement resulting from the n-pixel-shift R D F , the extended A C method substantially overcomes the drawbacks of the original A C method, and as such is able to handle higher order non-linear terms in a polynomial expansion of phase variation in M R I . The principle of the extended A C method can be illustrated by an example of secondorder phase error correction. Given a noise-free 256x256 M R image So with a zerothorder phase error a, a first-order linear phase error bx, and a second-order phase error cx 2 along x direction, we can write the idealized image as S(x,y) - S (x,y)cxp[i(a 0 + bx + cx )] 2 (3.12) To correct the phase error, the extended A C method is divided into four steps as follows. (1) The first step is to choose a proper total pixel shift. Increasing the shift has been found to enhance signal substantially; however, the enhancement is counterbalanced by a greater loss of signal points due to the greater pixel shift and thus increased probability of failed statistics. This is because points at the boundaries between tissue and background noise in R D F may not fulfill the criterion that any useful point in R D F must be obtained from two signal points. A l s o , phase wrapping in R D F may occur as a result of using a larger shift. Therefore, an appropriate total shift size should be a compromise between the 60 Chapter 3 IR Imaging with Non-linear Phase Correction enhancement of signal strength and the loss of signal points. In our practice, an image is firstly processed by thresholding to yield a mask image. Then the mask image is used to calculate for x-directional n-pixel-shift R D F . The pixel-shift increases from 1-pixel to full F O V size until a cutoff shift is found, at which the product of the pixel-shift and the square-root of the number of signal points in R D F is maximized. The cutoff shift is used as the total shift size for x direction. Likewise, the total shift size for y direction is determined in a similar fashion. This is considered as a reasonable compromise for most spin-echo scans. Other shifts might be appropriate for other cases. (2) W i t h the total shift size being determined, a proper polynomial is chosen to model phase error by determining the highest order of the polynomial. For the example depicted by E q . (3.12), a series of R D F s are obtained by application of E q . (3.11) in x direction to estimate the order. B y measuring relative variations of the phase maps of these R D F s , the highest order of the polynomial needed w i l l be found to be 2 because the second-order R D F has a minimum relative phase variation, corresponding to a most uniform phase map. A n y other R D F has a greater relative phase variation because of either existing higher order polynomial phase error or dominating noise. For example, first-order R D F has a greater relative phase variation since it contains a linear phase ramp caused by cx , and third-order R D F has a greater relative phase variation because noise 2 becomes dominating. 61 Chapter 3 IR Imaging with Non-linear Phase Correction In this work, the relative phase variation is quantitatively measured by a parameter named "weighted relative phase variation ( W R P V ) " expressed by E q . (3.13). The proper cutoff order of the polynomial used is found when the W R P V is minimized. 2ZW(x,y)[AP(x,y)} 2 WRPV = ^ — (3.13) ZW(x,y)P x,y where P denotes the phase of complex mean value of normalized R D F A(x,y) and is expressed as P = arg(A), where A = X x, \ y ^ A(x, y) > AP(x, y) is a phase difference defined I by arg[A • A* (x, y)]; and W (x, y) = \A(x,y)\ represents a weighting coefficient defined as the magnitude of the R D F . (3) Once the second-order polynomial is chosen to model phase error, phase correction begins from the highest order. The second-order phase correction can be broken into three steps as follows: First, the first-order x-directional n-pixel-shift R D F denoted as A ^ „ and the second-order x-directional n-pixel-shift R D F denoted as A®\ are respectively calculated as A (x,y) xn * = SQ(x,y)S (x-n,y)exp[i(bn + 2ncx-n 0 62 2 )] (3.14) Chapter 3 IR Imaging with Non-linear Phase Correction * A S (*» 30 = x,n (*. y)A A xn = SQ(X, y ) S Second, compared 0 (x-n,y)S 0 with the corresponding angle of A^ (x,y) n (x - n, y) (x-n, angle y)5 (x-2n,y)exp(/2n c) 2 0 of statistical expectation of A^\(x, y ) , the shown in E q . (3.15) is enhanced by a factor of n , so 2 the coefficient c can be accurately determined by 2n c = a r g { £ [ A ^ , y ) ] } 2 (3.16) ( If the phase variation in S(x,y) is greater than ± jc/n, the use of n-pixel-shift w i l l surely cause phase wrapping in the first-order R D F A A (x,y)will xn x n ( x , y ) ; however, the phase wrapping in not affect the estimation of the second-order polynomial coefficient c because c is usually small enough to avoid phase wrapping i n the second-order R D F ^x!n ( > y) • M ° x r e specifically, c usually satisfies -n < 2n c < +7i. Therefore, i n practice 2 only the estimation of linear phase error involves a concern of phase wrapping because non-linear polynomial coefficients are often small. Third, the estimated second-order phase error is removed from the original image S(x,y) by multiplying it by the complex conjugate of the estimated phase error exp(-/cx ). 63 Chapter 3 IR Imaging with Non-linear Phase Correction After the second-order phase error is removed, bx and a can be estimated and removed to yield a final phase corrected image denoted as S'(x, y) in a fashion similar to the original A C method. (4) The quality of phase correction can be evaluated with a parameter named "refocusing ratio ( R R ) " defined by E q . (3.17). I 2ZS'(x,y)\ RR = - ^ f T.\S'(x,y)\ (3.17) x,y where S'(x, y) denotes the final phase corrected complex M R image. The refocusing ratio is the normalized total magnetization used to characterize the refocusing of M R signals. When phase error exists, M R signals are dephased. When M R signals are completely dephased, R R reaches a minimum possible value of 0. When phase correction is applied to remove the phase error, the M R signals begin to refocus with an increased R R value. A n ideal R R value of 1 (unity) implies complete phase error removal. Therefore, R R has a possible value in the range of [0,1]. Generally, the greater value R R has, the better the phase correction can be achieved. A s described previously, a polynomial is used for phase correction. The choice of the polynomial order is found to be crucial since it greatly affects the performance of phase correction. The proposed W R P V , as a statistical model, usually can successfully 64 Chapter 3 IR Imaging with Non-linear Phase determine the cutoff order. Correction In addition to the statistical model, we can use a practical strategy that tries several polynomials around the cutoff order to ensure the optimal performance of phase correction, by making use of the computing efficiency of the extended A C method. In this work, the practical strategy uses five polynomials in 2D and uses the R R to evaluate the various phase corrections achieved by the different polynomials. Finally, the best result with the maximum R R is output. In summary, W R P V is first used to search the cutoff order in x and the order in y separately. Then R R is used to search the best order for both x and y jointly around the orders found in the first step. Although only second-order x-directional phase correction is demonstrated, it is straightforward to apply the extended A C method to correct for higher order phase errors 3 2 such as x , other directional phase errors such as y and cross terms such as xy. For example, to correct for cross-term phase error dxy, the following y-directional R D F can be calculated after a first-order x-directional R D F A*„ is obtained. %n (x, y) = A „ (x, y)A\ A X; n (x, y-n) (3.18) After this, the coefficient d can be directly estimated from the second-order R D F A ^ „ i f there is no other higher order cross terms. The extended A C method is summarized by a flowchart shown in Figure 3.2. 65 Chapter 3 IR Imaging with Non-linear Phase Correction ^ I n p u t original complex image Determine a total shift size I Pre-select a small range o f 5 polynomials by using W R P V I Start phase correction from the highest order to the 0 order for each polynomial 1 7 lh Evaluate phase correction results by using R R ^ Output the optimal result Figure 3.2: Flowchart of the extended A C method for non-linear phase error correction based on the n-pixel-shift R D F . 66 Chapter 3 IR Imaging with Non-linear Phase Correction 3.2.4 Improved Weighted Least Squares Method In 1999, Liang proposed a polynomial-model-based method for phase unwrapping [66]. The proposed method consists of three major steps: (1) calculating the phase derivatives along both x and y directions using a fast Fourier transform (FFT) algorithm; (2) calculating the polynomial coefficients using a weighted least squares (WLS) method; (3) calculating a residual function [66]. The n-pixel-shift R D F can also be applied to improve the W L S method proposed by Liang [66], which fits polynomials to phase derivative fields. As stated by Liang, the phase derivatives dP(x,y)/dx and dP(x, y)/dy were calculated by using a FFT algorithm and the calculations may not necessarily yield the expected real numbers. From the resulting complex numbers, real parts were used as the approximate phase derivatives [66]. Obtaining phase derivatives in this way is found to be inaccurate, and therefore suboptimal. In this work, we developed an improved WLS method for M R I phase correction based on n-pixel-shift R D F (improved W L S method) [2]. This method fits polynomials to the RDFs of a complex M R image instead of the approximate phase derivatives used by Liang. In addition to the WLS method, Liang proposed a residual function to correct some minor phase error that cannot be represented by a smooth polynomial [66]. 67 However, the Chapter 3 IR Imaging with Non-linear Phase Correction residual function is assumed to be in the principal value range (-TC, TC]. This assumption is only acceptable when the polynomial used is appropriate for the phase variation in M R I . If the phase variation contains higher order terms than those in the polynomial used, the application of the residual function cannot account for the remaining error in the original W L S method. This is because the higher order polynomial terms are generally not in the principal value range, and as such cannot be represented by a residual function. Since this paper focuses on the phase correction using polynomial fitting, the residual function w i l l not be considered. Similar to the extended A C method, the improved W L S method can be divided into the following four steps illustrated with the example with x-directional phase error described by E q . (3.12). (1) A „ of S(x,y) is calculated by using E q . (3.10). x (2) In the same way as the extended A C method, a second-order polynomial is determined by using W R P V defined by E q . (3.13). (3) The order of the second-order polynomial is shifted down by one order in x direction into a linear polynomial that is actually fitted to A xn by applying magnitude W L S algorithm. The W L S algorithm minimizes the difference between the linear polynomial 68 Chapter and A xn 3 IR Imaging with Non-linear Phase Correction in a magnitude weighted least squares sense. In a fashion similar to the extended A C method, phase error is removed after polynomial coefficients are estimated. (4) The phase correction can be evaluated by using refocusing ratio defined by E q . (3.17) in the same way as the extended A C method. Because of the computing efficiency of the extended W L S method, we may consider the application of a practical strategy that tries a small range of polynomials around the cutoff order. The various phase corrections yielded by the different polynomials are evaluated with R R to yield a best result with the maximum refocusing ratio. The two-dimensional improved W L S method is summarized by a flowchart shown in Figure 3.3. 69 Chapter 3 IR Imaging with Non-linear Phase Correction Input original complex image F i n d two R D F s o f the complex image i n x and y direction respectively I Pre-select a small range o f 5 polynomials by using W R P V I Shift down every polynomial by one order in x or y direction and fit the shifted polynomial to the corresponding R D F by W L S I Evaluate phase correction results by using R R I Output the optimal result Figure 3.3: Flowchart of the improved W L S method based on the n-pixel-shift R D F . 70 Chapter 3 IR Imaging with Non-linear Phase Correction 3.2.5 IR Image Reconstruction with Phase Correction T o demonstrate the feasibilities of the extended A C method and the improved W L S method, they are applied to Inversion Recovery (IR) image reconstruction, which requires a phase correction. Generally, an original IR image is complex-valued and can be written as [72] I(x, y) = M (x, y)0(x, y) exp[/ 0(x, y ) ] (3.19) where M(x,y) is a magnitude factor that is always non-negative, 0(x,y) is a polarity factor that can only be 1 or - 1 , and exp[i0(x,y)] is a phase error factor. The purpose of IR image reconstruction is to find M(x,y)0(x,y). T o accomplish this, the phase error factor must be removed from the original IR image. The proposed method to remove the phase error factor can be broken down into three steps. First, the original complex IR image is squared. factor. This complex square operation removes the polarity A s a result, the phase error also gets doubled but can still be handled by the proposed phase fitting algorithm. Second, non-linear phase calculation is applied to determine the doubled phase error factor, which is subsequently divided by 2 to yield an estimate of the actual phase error factor. Finally, the estimated phase error is removed from the original IR image by multiplying it by the complex conjugate of the estimated phase error. 71 Chapter 3 IR Imaging with Non-linear Phase Correction 3.2.6 E x p e r i m e n t s The performance of various phase correction methods was evaluated by using the existing patient data from a 1.5 T clinical scanner ( E D G E , Picker International, Cleveland, O H , U S A ) acquired with a multi-slice IR spin-echo sequence ( T R / T I / T E : 2000/600/30 ms). A l l of the phase correction methods were implemented in C programming language on a personal computer with 256 M B R A M and a clock rate of 9 0 0 M H z . For each slice with 2 5 6 x 2 5 6 pixels, phase correction by our proposed methods is accomplished within seconds. In this work, I define horizontal as the x direction and vertical as the y direction in all images. 3.3 Results The proposed phase correction methods based on the n-pixel-shift R D F are demonstrated with the following two in vivo examples. One example is from a transverse head scan of a patient with massive hydrocephalus. The contrast between brain tissues and cerebrospinal fluid ( C S F ) should be enhanced substantially by IR reconstruction. Figure 3.4(a) and (b) are respectively the real part and the phase map of the complex experimental IR image. Figures 3.4(c) and (d) are the real part and the phase map of the squared complex IR image with a doubled phase error 72 Chapter 3 IR Imaging with Non-linear Phase Correction factor but no longer with the polarity factor. In this work, various methods are applied to estimate the actual phase error. (a) (b) (c) (d) Figure 3.4: (a) and (b) are the real part and phase of the complex experimental IR image of a transverse head scan, (c) and (d) are the real part and phase of the squared complex IR image, in which polarity factor is removed but the phase is doubled. In applying the extended AC method to this case, the total shift sizes for x and y directions are automatically calculated to be 125-pixel and 134-pixel, respectively. Then, a series of R D F s are obtained by repeated application of E q . (3.10) in both x and y directions. For each R D F , the W R P V is evaluated as shown in Table 3.1. The cutoff order of the polynomial is found to be 2 in both x and y directions, corresponding to minimum W R P V values. Next, the second-order polynomial and its four neighboring polynomials illustrated in Table 3.2 are used in phase calculation and correction to yield 73 Chapter 3 IR Imaging with Non-linear Phase Correction five results. Table 3.2 shows the corresponding refocusing ratios for them. Finally, the best result with the maximum refocusing ratio of 0.830, which is based on the squared complex IR image, is output. Results from different methods are shown in Figure 3.5 for comparison. Figure 3.5(a) is the real part of the original experimental IR image that is the same as Figure 3.4(a). Figure 3.5(b) is the real image after phase correction by the linear A C method based on the 1-pixel-shift R D F , which shows a slow varying remaining phase error across. Figure 3.5(c) is the real image after phase correction by the second-order A C method based on the 1-pixel-shift R D F , which gives biased polynomial coefficients because of noise and yields a worse result with extremely rapid phase change across. Figure 3.5(d) is the real image after phase correction by the second-order W L S method based on the 1-pixel-shift R D F , which shows similar contrast as in Figure 3.5(b). Figure 3.5(e) shows the real image after phase correction using the extended A C method that applies second-order phase correction in x direction and the third-order phase correction in y direction. The total shift size is 125-pixel for x direction while the total shift size is 134-pixel for y direction. The pixel shift sizes for each order respectively are 63-pixel and 45-pixel in x and y directions. Figure 3.5 (f) shows the real image after phase correction using the second-order improved W L S method based on 10-pixel shift R D F , which could avoid possible phase wrapping in R D F and give a considerable phase signal enhancement. These results show that the proposed methods offer much more accurate phase correction 74 Chapter 3 IR Imaging with Non-linear Phase Correction and are able to achieve the desired contrast of an IR image. The improvements shown i n the figures can be quantified consistently by the refocusing ratio as shown i n Table 3.3. Order of RDF WRPV in x-directional WRPV in y-dircctional RDF RDF 0 182.042 I 0.960' 8.651 2 0.481 0.59.1 3 14.150 9.730 182.042 Table 3.1: W R P V values for a series of R D F s are obtained by repeated application of E q . (3.10) with the total shifts of 125-pixel and 134-pixel, respectively, in x and y directions for the squared complex IR image from a transverse head scan. Along x ^direction Along > V direction Polynomial order: 1 Polynomial order: 2 Polynomial order: .1 Polynomial order: 2 Polynomial order: 3 RR: 0.746 RR: 0.685 RR: 0.818 Polynomial order: 3 RR: 0.824 RR: 0.830, Table 3.2: Refocusing ratios for extended A C method applied to the squared complex IR image shown in Figure 3.4 (c) with the corresponding polynomials. 75 Chapter 3 IR Imaging with Non-linear Phase ........... Correction ... (a) (b) (c) \ / (e) (d) J / (f) Figure 3.5: In vivo real IR image from a transverse head scan and its corresponding phase corrections with various methods. Phase calculation and correction method Refocusing ratio Squared complex IR image 0.013 Original AC method 0.223 2 order extended AC method based on 1-pixei-shift RDF 0.002 2 order improved WLS method based on 1-pixel-shift RDF 0.136 2 order x-directional and 3"' order ydirectional extended AC method based on 125-pixel-shift and 134-pixel-shift. 0.830 2 order improved WLS method based on 10-pixel-shift RDFs 0.813 nd nd Table 3.3: Refocusing ratios for various phase correction methods applied to the squared complex IR image from a transverse head scan shown in Figure 3.4 (c). 76 Chapter 3 IR Imaging with Non-linear Phase Correction Another example from a coronal head scan of a patient is shown in Figure 3.6. 3.6(a) is the real part of the original experimental IR image. Figure Figure 3.6 (b) is the real image after phase correction by the linear A C method based on the 1-pixel-shift R D F , which shows unsatisfactory phase correction in the head region of the image. Figure 3.6 (c) is the real image after phase correction by the second-order A C method based on the 1-pixel-shift R D F . Again, the contrast in the image is entirely distorted by an inaccurate phase correction. Figure 3.6 (d) is the real image after phase correction by the secondorder W L S method based on the 1-pixel-shift R D F , which shows large remaining phase error in the head. Figure 3.6(e) shows the real image after phase correction using the second-order extended A C method. The total pixel-shifts are 127-pixel and 184-pixel respectively for x and y directions. Figure 3.6 (f) shows the real image after phase correction using the second-order improved W L S method based on the 20-pixel shift RDF. Both Figures 3.6(e) and 3.6 (f) show greatly improved contrast across the image with minimum phase error and maximum refocusing ratios shown in Table 3.4. 77 Chapter 3 IR Imaging with Non-linear Phase (d) Correction (0 (e) Figure 3.6: In vivo real IR image and its corresponding phase corrected images with various methods from a coronal head scan. Phase calculation and correction method Squared c o m p l e x IR image Refocusing ratio 0.101 O r i g i n a l A C method 0.353 2™ order extended A C method based 1 o n 1-pixel-shift R D F T order i m p r o v e d W L S method based A o n 1-pixel-shift R D F 2** order extended A C method with total 0.033 0.375 1 shifts o f 127-pixel and 184-pixel 2 n<i order i m p r o v e d W L S method based o n 20-pixel-shift R D F 0.654 0.553 Table 3.4: Refocusing ratios for various phase correction methods applied to the squared complex IR image from a coronal head scan. 78 Chapter 3 IR Imaging with Non-linear Phase Correction 3.4 Discussion After introducing the concept of a general n-pixel-shift rotational differential field ( R D F ) , this work demonstrated that a larger pixel shift n can effectively enhance signal without amplifying the corresponding noise, allowing for an improved phase correction. This is the key idea of this work. Based on the n-pixel-shift R D F , we have extended the A C method to handle non-linear phase error, and have improved the W L S method. Each of the proposed methods can efficiently accomplish a 2 D phase correction within seconds. Although n-pixel-shift R D F loses some signal points at the boundaries between tissue and background, the signal loss can generally be negligible because usually an M R image has a very large number of signal points. Another potential complication is that a possible phase wrapping in n-pixel-shift R D F may be caused by the use of a large pixel shift. In our practice, all such undesirable phase wrapping occurs in the first-order n-pixel-shift R D F because the coefficients of the first-order phase error terms are usually much greater than those of higher order phase error terms. T o handle this complication, 1-pixel-shift can be used to estimate the first-order phase error coefficients, which is found very effective. The result of phase correction is known to be sensitive to the choice of the polynomial cutoff order. In this work, we have demonstrated how to use W R P V to determine the cutoff order. T o optimize phase correction, we used a small range of polynomials obtained by shifting the cutoff polynomial up and down by one order in both x and y 79 Chapter 3 IR Imaging with Non-linear Phase Correction directions instead of only one single polynomial. The final result is selected when its corresponding refocusing ratio is maximized, which ensures the best achievable result based on the n-pixel-shift R D F with a given total shift size. The W R P V and the optimization criterion were first developed for the extended A C method and have been found to be effective for the improved W L S method as well. This can be explained by the idea that the polynomial cutoff order is predominantly determined by the M R image itself when similar methods are used. Additionally, the use of a small range of polynomials can ensure the quality of phase correction against the minor errors caused by the W R P V order estimation. Finally, this work introduced the refocusing ratio as an ultimate measure for the quality of phase correction. Although the extended A C method and the improved W L S method are presented as two different methods in this paper, they are both categorized as polynomial fitting and share the n-pixel-shift R D F in common. The key difference between them lies in the order of the n-pixel-shift R D F applied. For example, to fit an m^-order polynomial to an M R I phase map, the extended A C method applies a series of n-pixel-shift R D F from first-order R D F up to m -order R D F and fits a zeroth-order polynomial to each R D F , while the th improved W L S method uses two first-order n-pixel-shift R D F s and fits an ( m - l ) - o r d e r th polynomial to each R D F . The above analysis implies that the extended A C method and the improved W L S method can be combined to form a general polynomial fitting method. To fit an m -order polynomial to an M R I phase map, the general method th divides all polynomial coefficients into two parts: low order coefficients and high order 80 Chapter 3 IR Imaging with Non-linear Phase Correction coefficients. The former can be estimated by using the extended A C method, while the latter can be estimated by using the improved W L S method. The balancing point between low and high orders should be optimized to yield a best result rather than being arbitrarily chosen to be m used in the extended A C method, or 1 used in the improved W L S method. However, this more general approach of polynomial fitting may not be necessary for the data studied, since both the two extreme cases — the extended A C method and the improved W L S method — have successfully achieved similar and satisfactory results. 81 Chapter 4 Simplified SPEED with Applications in MRA Chapter 4 Simplified Skipped Phase Encoding and Edge Deghosting (SPEED) with Applications in M R A Phase information in M R I also has great potential use in fast imaging. A successful example of fast imaging using phase information is "Skipped Phase Encoding and Edge Deghosting ( S P E E D ) " recently proposed by X i a n g [3]. S P E E D reduces scan time while suppressing associated aliasing ghosts. Unlike parallel imaging, S P E E D is able to accelerate M R I even with a single coil. S P E E D has been tested and demonstrated to be efficient and robust [3]; however, it is yet to be developed. In my doctoral thesis, I present new development based on the general principle of S P E E D . (1) First, since the current principle of S P E E D has been illustrated only in 2 D , I extend the theory of S P E E D from 2 D to 3 D for more flexible reduction of scan time. (2) Second, S P E E D , as an independent method, is able to achieve improved performance by combining itself with other existing methods. In my doctoral thesis, a combination of Half-Fourier imaging and S P E E D is proposed for further scan-time reduction and therefore is named Half-Fourier-SPEED. 82 Chapter 4 Simplified SPEED with Applications in MRA (3) Third, S P E E D has been demonstrated to reduce scan time considerably with typical M R I data. In this thesis, I simplify S P E E D to further enhance its efficiency by taking advantage of the sparseness of M R angiography ( M R A ) data [4]. (4) Finally, a new parallel imaging method is proposed, which highly accelerates M R I by Skipped Phase Encoding and Edge Deghosting with Array C o i l Enhancement and thus is named " S P E E D - A C E " [5]. S P E E D - A C E combines the principle of S P E E D with that of sensitivity encoding using multiple coils, and thereby can achieve an undersampling factor even greater than the number of receiver coils, which is impossible with conventional parallel imaging methods [5]. To assist in understanding the work in the following chapters, Simplified S P E E D w i l l firstly be discussed because it is a simple version of S P E E D and thus has the most straightforward reconstruction algorithm [4]. Then, the original S P E E D w i l l be reviewed and its two selected extensions — S P E E D with two P E directions in 3 D and HalfFourier-SPEED — w i l l be discussed in Chapter 5. After this, S P E E D - A C E w i l l be explored in Chapter 6. Finally, signal-to-noise ratio ( S N R ) in S P E E D w i l l be studied and a novel approach for S N R improvement w i l l be proposed in Chapter 7. 83 Chapter 4 Simplified SPEED with Applications in MRA 4.1 Introduction Scan time reduction is desirable in M R angiography ( M R A ) , which is often performed in 3D [73]. In the past, the acquisition of M R I data has been greatly speeded up by increasing gradient strength and slew rate. However, the use of this method is limited because of relevant physiological and technological concerns. A s alternatives, many fast imaging strategies have been proposed to accelerate M R A without the need of increasing the performance of gradient. complete Most of them partially sample k-space and reconstruct a image through a non-standard reconstruction. These strategies include undersampled projection reconstruction (PR) [74,75], time-resolved imaging of contrast kinetics ( T R I C K S ) [76], a hybrid of undersampled P R and T R I C K S ( P R - T R I C K S ) [77,78], and parallel imaging including S M A S H [48,79] and S E N S E [49,80,81]. In particular, parallel imaging reduces scan time through skipped sampling of k-space with multiple receiver coils in parallel [48,49,79-81]; and undersampled P R accelerates M R A without intolerable artifacts by reducing the number of projections [74,75,77,78]. Undersampled P R is based on fact that M R A data usually have a very sparse signal distribution. Similarly, based on this characteristic of M R A data, we simplify a fast imaging method named S P E E D [3] to enhance its efficiency. Previously demonstrated to considerably reduce scan time for typical M R I data [3], S P E E D is simplified in this work to accelerate M R A [4]. Both S P E E D and its simplified version (Simplified S P E E D ) partially sample k-space, and separate the overlapping structures resulting from k-space undersampling in a manner similar to that used in motional deghosting [83-84] and 84 Chapter 4 Simplified SPEED with Applications simple spectroscopic imaging [30-33]. in MRA Compared with S P E E D , Simplified S P E E D requires less scan time and uses more straightforward reconstruction [4]. In this chapter, Simplified S P E E D is firstly described. Then it is further demonstrated that Simplified S P E E D can be extended to 3 D for more flexibility in reducing scan time, and that it can be combined with other existing techniques such as Half-Fourier imaging [39] for improved performance. Finally, a generalization of Simplified S P E E D is discussed. 4.2 Method 4.2.1 S i m p l i f i e d S P E E D for Fast M R A In this work, S P E E D is simplified to enhance its efficiency by taking advantage of the sparseness of M R A data [4]. The sparse signal distribution in M R A data is in a close analogy to what happens when several one-dimensional wires are twisted in threedimensional space. G i v e n the sparseness of M R A data, the popular Maximum-IntensityProjection (MIP) algorithm used in M R A considers only the vessel with the brightest signal even when structures overlap [82]. Similarly, Simplified S P E E D accelerates M R A based on the sparseness of M R A data. A s shown in Figure 4.1, Simplified S P E E D partially samples k-space with a skip size of N to obtain two interleaved data sets, and then reconstructs them respectively into two ghosted images: L(r) and I (r), each with 2 ghost overlapping up to N layers. G i v e n the sparse signal distribution of M R A data, Ii(r) and I (r) are modeled with a single-layer structure, and they are subsequently deghosted 2 85 Chapter 4 Simplified SPEED with Applications in MRA by using a least-square-error ( L S E ) algorithm. Finally, a deghosted M R A image I (r) is 0 yielded, and a residual map is output. W i t h the residual map, the quality of the deghosted M R A image can always be monitored. In this way, scan time is reduced by a factor of N/2, corresponding to an acceleration factor of N/2. The algorithm of Simplified S P E E D can be further illustrated with a 2D example as follows, •::.::;z-;::::: S(k)| S (k) 2 V • 1 Figure 4.1: K-space sampling scheme and data flow of Simplified S P E E D with one P E direction. 2D full k-space S (k) is sparsely sampled into two interleaved sparse k-space data sets: S,(k) and S (k), each with the same skip size of N and a different relative shift in P E . After inverse 2DFT, 2 ghosted images —I r) and I (r)— are obtained and subsequently combined into a deghosted image Io(r). 0 2 iv 86 2 Chapter 4 Simplified SPEED with Applications in MRA Firstly, two interleaved data sets are obtained: one is acquired by sampling k-space at every N t h P E step; the other is sampled in the same way but with a different relative shift in P E chosen from 0 to N - l . The sampled data sets are inverse Fourier transformed separately into two ghosted images, each associated with the true image and N - l aliasing ghosts. For simplicity, the true image and N - l aliasing ghosts are all considered as ghosts in this work, and as such results in " N " aliasing ghosts. It is well known that M R A data usually have much sparser signal distribution than typical M R I data. This implies that ghosted M R A images usually have much less ghost overlapping than their typical M R I counterparts. Therefore, differential filtering that is used to reduce ghost overlapping in S P E E D is not used in Simplified S P E E D , and the ghosted M R A images are directly used in the following deghosting. Inspired by the M I P algorithm used in M R A [82], we model the ghosted M R A images with a single dominant layer of ghosts, and write them at each pixel as h=P2\G (4-1) h=PdlG (4.2) H n where G is the dominant layer of ghosts; Pf is a ghost phasor known to have a form of n i2itd (e N )" , where n = 0,1.. . N - l is the order of ghost depending on its location relative to its source, and d is a known relative sampling shift in P E . The unique ghost order parameter n is named as a ghost order index in this work. Although n at each pixel is not 87 Chapter 4 Simplified SPEED with Applications in MPA fully determined, it is known to be an integer chosen from 0 to N - l . Since there are two equations but only one unknown variable G , there must be a L S E solution for each N possible n. The L S E solution GLSE can be analytically derived as follows, 1) First, square error is defined and is given by Square-Error = h = + ~ n d\ G P 2 KJ [ReC/^V ) - R e ( G „ ) ] [Re(I P *)- 2 r + r [Im(/,/>V)- Im(G„ ) ] 2 d 2 2 rf Re(G„ ) ] + [lm(I P *)- n 2 2 d2 + UP?- ~ n d2 Im(G„ ) ] n 2 d 2 + (4.3) 2 where the operators R e ( ) and I m ( ) return the real and imaginary parts of the vector in the brackets, respectively. 2) Second, by considering G as a continuous variable, two partial derivatives of square N error are respectively calculated with respect to Re(G„) and Im(G„). 3) Third, to minimize the square error i n real and imaginary dimensions, the two partial derivatives of square error are set to zero. The real and imaginary parts o f GLSE can therefore be solved as Re(/,/>;, ) + R e ( / P ; 2 and 88 2 ) (4.4) Chapter 4 Simplified SPEED with Applications in MRA Im^/ft j+Im(/ P; 2 2 ) (4.5) 4) Finally, combining Eqs.(4.4) and (4.5) yields the L S E solution 7 W i t h the resolved GLSE, GLSE (4.6) LSE its L S E can be directly obtained by substituting GLSE for G N in E q . (4.3). LSE = h ~GLSEPS\ + 12 h ~GusE^d2 ~hPd\Pd2 j + 2 1 - r P" p n \^d2 d\ l r (4.7) B y using E q . (4.7) at each pixel, the correct value of n can be found by trying all of its N possibilities ranging from 0 to N - l until L S E is minimized. Simultaneously, the correct solution of G can be determined by using E q . (4.6) and the correct value of n. After this, n G is sorted out according to the resolved n and thereby yields N separate deghosted n M R A images. More specifically, G is assorted by assigning its values to the n t h n deghosted image on a pixel-by-pixel basis, resulting in N separate deghosted images with known displacement in P E direction. Each of the separate deghosted images denoted by G„, n = 0,1.. . N - l , can be spatially shifted by — x FOV to its origin location. Finally, the N deghosted images are spatially registered and averaged to yield a final M R A image with 89 Chapter 4 Simplified SPEED with Applications in MRA reduced noise and artifacts. A s shown in Figure 4.2, the above deghosting process can be intuitively illustrated by the following data flow at N =3. In Simplified S P E E D , the skip size of N must be a prime number to avoid "ghost phase degeneracy" [3,83,84], which refers to the situation where ghosts of different orders share the same phase. Besides the resolved n and G„, a residual map is obtained and output for quality control. Each pixel value in the residual map is calculated as the square-root of L S E in E q . (4.7). The residual map can be used to monitor the quality of deghosting. If a single layer is able to adequately model ghosted images, the residual should be close to zero; otherwise, a greater residual w i l l be caused by more than one layer of ghost overlapping. For example, the residual w i l l be intolerably strong i f ghosted M R A images are too complex to be modeled by a single-layer model. T o reduce the residual to a tolerable level or noise level, additional data should be acquired and additional layers should be included into the ghost model to refine the algorithm. This w i l l be discussed later in this work. O f course, a large residue can also be caused by noise and other image artifacts. In summary, given the sparseness of M R A data, Simplified S P E E D models ghosted M R A images with a single-layer structure like M I P does. Based on the single-layer modeling, a deghosted M R A image can be obtained with only two interleaved data sets. In this way, scan time can be reduced by a factor of N / 2 . The algorithm of Simplified S P E E D is summarized by a flowchart shown in Figure 4.3. 90 Chapter 4 Simplified SPEED with Applications in MRA Ur) I,(r) » + » G„(r) GhOSl order index map n(r) Figure 4.2: Data flow of deghosting in Simplified S P E E D with one P E direction at N = 3. A L S E solution G (r) and its associated ghost order map n(r) are firstly solved with two ghosted images: I,(r) and I (r). Then, G is sorted out according to n(r) and thereby yields N = 3 separate deghosted M R A images: G ( r ) , G[(r), and G ( r ) . Finally, they are spatially registered to yield G ^ r ) , G ( r ) , and G ( r ) and averaged into a deghosted image I (r). n 2 n 0 r) 2 r2 0 91 Chapter 4 Simplified SPEED with Applications in MRA Input 2 interleaved datasets S , & S 2 Reconstruct into 2 ghosted images I, & I by direct F T 2 A Y i e l d the L S E solution o f one dominating ghost layer G„ through N trials o f n at each pixel I I Sort out N separate ghosts G„ according to associated n Register the N ghosts G„and average them into a final deghosted image I„ I Output the final image I 0 Figure 4.3: Flowchart of Simplified S P E E D with one P E direction used to achieve an acceleration factor of N / 2 with a single coil. 92 Chapter 4 Simplified SPEED with Applications in MRA 4.2.2 Extension of Simplified SPEED to Two PE Directions Since M R A is often performed in 3 D , the application of Simplified S P E E D needs to be extended from 2 D to 3D to allow further reduction of scan time. Three-dimensional data acquisition is usually accomplished with P E along two orthogonal directions. Therefore, M R A in 3 D can be accelerated by using skipped P E along both of the directions. It can be illustrated with an example as follows: A s shown in Figure 4.4, two relatively shifted and interleaved data sets are sampled respectively with a skip size of M in the horizontal direction and a skip size of N in the vertical direction. The data sets are then inverse Fourier transformed respectively into two ghosted images, each with M x N aliasing ghosts. Based on a single-layer model, the two ghosted images can be written as h=P QlG nm (8) h=P Q?G nm (9) n d n d where G nm is the dominant ghost layer of ghosts, and P Q™ d [2nd ghost phasor (e N is known to have the form of 12K S ) (e n M ) m where m and n are different orders of ghost, depending on their relative location in horizontal direction and vertical direction. The order m must be an integer between 0 and M-l; likewise, the order n must be between 0 and A M . Both 93 Chapter 4 Simplified SPEED with Applications m and n are to be determined. in MRA The integers s and d represent respectively relative shifts in P E along horizontal and vertical directions, chosen from 0 to M-\ or 0 to A -1. 7 Figure 4.4: K-space sampling scheme and data flow of Simplified S P E E D with two P E directions. 2 D cross-section of 3D k-space S (k) is sparsely sampled into two interleaved sparse k-space data sets: S^k) and S (k), each with a skip size of M in horizontal direction and a skip size of N in vertical direction. After inverse 2 D F T , two ghosted images I,(r) and I (r) are obtained and subsequently combined into a deghosted image 0 2 2 Io(r). 94 Chapter 4 Simplified SPEED with Applications in MRA Similar to the deghosting algorithm used in Simplified S P E E D with one P E direction, M x N separate ghosts of different orders can be solved by minimizing L S E through trying all possible values of m and n on a pixel-by-pixel basis. According to the order of each resolved ghost, M x N ghosts are easily aligned by linear registration and are averaged into a single deghosted image with an acceleration factor of ( M x N ) / 2 . Similar to the 2 D imaging situation, M and N must be two different prime numbers in order to avoid ghost phase degeneracy [83,84]. 4.2.3 C o m b i n a t i o n o f S i m p l i f i e d S P E E D and H a l f - F o u r i e r Reconstruction Being an independent method, Simplified S P E E D can be combined with other fast imaging techniques to achieve improved performance. This work demonstrates a fairly straightforward combination of Simplified S P E E D and Half-Fourier imaging [39] (HalfFourier Simplified S P E E D ) as follows: T w o relatively shifted and interleaved data sets are sampled within slightly more than half of k-space. In our practice, data are sampled to span 60% of k-space, corresponding to a typical Half-Fourier coverage. The interleaved data sets are then processed by Simplified S P E E D to yield a final deghosted image. Although Half-Fourier imaging usually requires phase correction to reduce truncation artifacts caused by undersampling, we have found that such phase correction is unnecessary in Simplified S P E E D because truncation artifacts seem almost negligible in the deghosted M R A image. B y combining 95 Chapter 4 Simplified SPEED with Applications in MRA Simplified S P E E D with Half-Fourier imaging, we can further reduce scan time by a factor of nearly 2. 4.2.4 G e n e r a l i z a t i o n o f S i m p l i f i e d S P E E D Although the simple single-layer model is usually applicable in sparse M R A , more layers can be added into the model to refine the algorithm for Simplified S P E E D . For example, in Simplified S P E E D with one P E direction, to handle two layers of ghost overlapping, kspace can be sparsely sampled into three interleaved data sets, each with the same skip size of N steps and a different relative shifts in P E . The three data sets are then reconstructed by direct inverse F T into three ghosted images denoted respectively as I\, h, and 1$. In this way, for each pixel, Eqs. (4.1-4.2) can be extended to I =PdtG +P G (4.10) h=P&G x+P£G (4.11) h=Pd\G +p£G (4.12) n2 l nl d l n2 n n{ nl nl where G„i and G 2 stand for the two dominant layers of ghosts. n In this work, the two dominant ghosts are written into (G„i, G 2) while their orders are denoted as n (nl,n2). Considering the ghost order index (nl, n2) as parameters, we can write Eqs. (4.10-4.12) into a matrix form as follows: 96 Chapter 4 Simplified SPEED with Applications in MRA I =M (4.13) G where / =(I\, h, h) is a known vector, T = G (G„i,G„2) is an unknown vector to be found, r and M is a 3x2 matrix given by >n2 dl pnl dl M = pnl dl pnl pnl d3 (4.14) d3 Eqs. (4.10-4.12) are three equations with two unknowns, indicating an over-determined linear system. Therefore, a L S E solution vector GLSE can be solved for each possible pair of (nf n2) as follows: LSE T = LI (4.15) where L is a 2x3 "pseudo-inverse" matrix given explicitly by L = M M\M r 1 9-QQ 3P f-QP f n d 3^2 d 3Pd f-Q*P f n -QPdi ^di n d -Q*P£ 3Pdi-Qp f n (4.16) d ip f-Q*Pdi d where the vector Q is defined as pnl _ p«l G ~ dl r dl r , pnl + r d l nl p r dl 97 nl pnl p d3 +r d3 r (4.17) Chapter 4 Simplified SPEED with Applications in MRA and the superscripts "t", and " * " denote respectively the operations of the matrix transposed complex conjugate, the matrix inverse, and the complex conjugate. the L S E solution vector ~(G HL )LSE )LSE. GLSE Finally, can be explicitly expressed as, 3P i -QP f -Q> f 3P f -Q*p i -QP f 1 d d ~9-QQ* d d d d -OPS? 3P£ 3P f -Q> f d d The L S E or residual energy for this solution can be given by LSE = energy[I-MG ] (4.19) LSE where the operator energy returns the square sum of all vector elements in the square bracket. E q . (4.15) is a general expression of GLSE, which can be reduced to Eq. (4.6) for a system with two equations and one unknown, and to E q . (4.18) for a system with three equations and two unknowns. B y trying all possible pairs of ( n l , nl), the correct solutions of ( G \ , G i) and ( n l , n2) can N n be simultaneously determined when L S E is minimized at each pixel [3]. Given the resolved ( n l , nl), ( G \ , G,a) can be sorted out, spatially registered, and averaged to yield N a final deghosted M R A image, resulting in an acceleration factor of N/3. Similarly, to handle two layers of ghost overlapping in Simplified S P E E D with two P E directions, Eqs. (4.8-4.9) can be extended to, 98 Chapter 4 Simplified SPEED with Applications \ ~ dl l r 2~ l Us\ n\ml u r r + r f 2 &S2 Lr r 3 ~ d3^S3 nlml 1 d2 &S1 + ^d2^is2 nim\ t in MRA Lr d3 &S3 + r n2m2 Lr (4.20) n2m2 (4.21) n2m2 (4.22) U Lr where I\, h, and h are the three ghosted images obtained from three interleaved and relative shifted data sets, each with a skip size of M in horizontal direction and a skip size of N in vertical direction; G \ \ and G 2 2 are the two dominant layers of ghosts. n m n m In a fashion similar to Simplified S P E E D with one P E direction, a final deghosted M R A image can be obtained through deghosting, resulting in an acceleration factor of ( M x N ) / 3 . So far we have described the extension from a single-layer model to a double-layer model, it is straightforward to extend Eqs.(4.10-4.12) and Eqs.(4.20-4.22) to handle three or more layers of ghost overlapping. Let us consider an extreme situation where k-space is sampled into N interleaved and relatively shifted data sets, each with a skip size of N steps in one P E , corresponding to a full k-space sampling. It is well known that sampling k-space with a skip size of N results in ghost overlapping up to N layers. T o handle the ghost overlapping, Simplified S P E E D can use an N-layer model with N interleaved data sets (i.e. full k-space data) to yield a perfect deghosted image, which is identical to the gold standard reconstructed directly by inverse F T from full k-space data. This implies that the difference between a Simplified S P E E D solution and the gold standard decreases as the number o f layers i n the model increases. The difference w i l l vanish i f an N-layer 99 Chapter 4 Simplified SPEED with Applications in MRA model and N interleaved data sets are used. In other words, Simplified S P E E D solutions are becoming closer and closer to the gold standard i f more and more layers in the ghost model are used, analogous to Taylor expansion. In this work, the multi-layer expansion process is named Ghost Layer Unfolding Expansion ( G L U E ) . In Taylor expansion, a function is expressed as a polynomial with many terms. Similarly, sparse ghosted M R A images can be expressed through a G L U E series by progressively adding more layers of ghosts for increased accuracy. In Taylor expansion, a satisfactory approximation can often be achieved by using the first few terms. Likewise, given the sparseness of M R A data, Simplified S P E E D can be used to achieve a satisfactory deghosted M R A image with a single-layer model or a double-layer model. The L S E deghosting algorithm of Simplified S P E E D ensures that the solutions based on a singlelayer model or a double-layer model must be the most dominant one or two layers of ghosts in the multiple layers of ghost overlapping. Moreover, image quality is monitored by a residual map, and more data can always be added later to improve the quality of the result. Therefore, to optimize its performance, Simplified S P E E D should be used in realtime M R A . In real-time M R A , the solution of Simplified S P E E D should be firstly solved based on a single-layer model with only N / 2 k-space data. Because the image quality is monitored and the reconstruction of Simplified S P E E D algorithm is fast, more data can be added in real time until a satisfactory result is achieved. In a word, Simplified S P E E D allows flexibility in reducing M R A scan time. 100 Chapter 4 Simplified SPEED with Applications in MRA To quantify image quality, the deghosted results can be evaluated with total relative error ( T R E ) defined by Z[M (x,y)-M (x,y)]< 0 TRE- g x,y 2ZM (x,y) (4.23) 0 where Mo(x,y) is the deghosted image, and M (x,y) g is the gold standard reconstructed from full k-space data, both in magnitude form. 4.2.5 E x p e r i m e n t s The fast M R A methods proposed in this work were tested with experimental data from a vasculature phantom, which mimicked blood vessels with a twisted long plastic hose filled with water as shown in Figure 4.5. The vasculature phantom was scanned on a 1.5T clinical scanner ( E D G E , Picker International, Cleveland, O H ) with a spin-echo sequence (matrix 256x256, F O V 30cm, T R 1000 ms, T E 10ms, slice thickness 5 m m , 2 averages, 25 slices). The full k-space data were saved and used to study the feasibilities of the proposed methods. Furthermore, the proposed methods were tested with a phase-contrast (PC) gradient-echo head scan on a volunteer by using a 1.5 T clinical scanner (TwinSpeed, G E Health Care, Milwaukee, WI) (matrix 256x256, F O V 24 cm, T R 34 ms, T E 6.4 ms, slice thickness 5 101 Chapter 4 Simplified SPEED with Applications mm, flip angle 20°, no average, 25 slices). in MRA The P C gradient-echo head scan was conducted with a typical clinical protocol. In the scan, two complete sets of images were sampled and saved: one was velocity-sensitized with V E N C = 8.6 cm/s, in which the velocity encoding direction was perpendicular to the sagittal plane; the other was a reference with the identical imaging parameters but without velocity encoding. To study the feasibility of Simplified S P E E D , both the velocity-sensitized and reference data sets were partially sampled with skipped phase encoding. The complex values of the undersampled velocity-sensitized data and the corresponding reference data were then subtracted on a pixel-by-pixel basis. The resulting data were still aliased, but without unwanted static tissue background. Finally, the resulting data were processed by Simplified S P E E D to yield a deghosted M R A image. The proposed fast M R A algorithms were implemented in C programming language on a personal computer with 2 G B R A M and a clock rate of 3.6 G H z . The reconstruction of a 256x256 image by Simplified S P E E D was typically accomplished within one second. Figure 4.5: The photo of the vasculature phantom (left) and the coverage of the 25 slices. 102 Chapter 4 Simplified SPEED with Applications in MRA 4.3 Results Figures 4.6-4.8 show the results from the vasculature phantom scan. For easy visualization, all images are shown here in magnitude although they are actually complex. A s shown in Figure 4.6, one of the 25 slices is reconstructed by Simplified S P E E D with a P E skip size of N = 7 steps. Figure 4.6 (a) is the gold standard image reconstructed from full k-space data. Figure 4.6 (b) is a noise field in which each pixel value is randomly taken from the background of the real and imaginary parts of the gold standard shown in Figure 4.6 (a). Figure 4.6 (c) is the summation of Figure 4.6 (a) and 4.6 (b). W i t h E q . (4.23), Figure 4.6 (c) is measured to yield a T R E of 0.10, which is normalized square root of error power of noise. Figure 4.6 (d) is one of the ghosted images. In the ghosted image, seven-fold aliasing ghosts of various orders are along the vertical P E direction, and most of them are seen as one dominant layer. Figure 4.6 (e) is the reconstructed image by Simplified S P E E D based on a single-layer model, with an acceleration factor of N / 2 = 3.5. Figure 4.6 (e) is measured with E q . (4.23) to have a T R E of 0.17. Figure 4.6 (f) is the corresponding residual map of Figure 4.6 (e) for quality control. Compared with Figure 4.6 (d), the residual map shows much less pixel intensities, suggesting the adequacy of the single-layer modeling. Figure 4.6 (g) is the reconstructed image by Simplified S P E E D based on a double-layer model, with an acceleration factor of N / 3 = 7/3 ~ 2.3. Figure 4.6 (g) is measured with E q . (4.23) to have a T R E of 0.09. The pixel intensities in its residual map are close to zero as shown in 103 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.6 (h). Compared with Figure 4.6 (e), the image quality as indicated by residual map and T R E is improved in Figure 4.6 (g), but scan time is increased by 50%. T o allow more information to be observed, the contrast in Figure 4.6 (a), (c), (e) and (g) are reduced with a lower window and level setting (-100, 100) as shown in Figure 4.6 (i), (j), (k) and (1), respectively. 104 Chapter 4 Simplified SPEED with Applications in MRA (i) Li) 00 (1) Figure 4.6: 2 D reconstructed images from the phantom scan, (a) is one of the 25 slices reconstructed from full k-space data and is considered as the gold standard, (b) is a noise field where each pixel value is randomly taken from the background of the real and imaginary parts of the gold standard shown in (a), (c) is the summation of (a) and (b). (c) has a T R E of 0.10. (d) is one of the ghosted images, each obtained by sampling kspace with a skip size of N = 7 along vertical direction, (e) is the deghosted image yielded by Simplified S P E E D based on a single-layer model with an acceleration factor of 7/2 (f) is the residual map of (e). (g) is the deghosted image yielded by Simplified S P E E D based on a double-layer model with an acceleration factor of 7/3 (h) is the residual map of (g). Compared with the gold standard (a), both (e) and (g) have the errors comparable to the noise level of the gold standard quantified by T R E . (e) has a T R E of 0.17 while (g) has a T R E of 0.09. T o allow more information to be observed, the contrast in (a), (c), (e) and (g) are reduced with a lower window and level setting (-100, 100) as shown in (i), (j), (k) and (1), respectively. 105 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.7 shows the results of the phantom scan reconstructed by Simplified S P E E D with one P E direction based on a double-layer model. T o handle double-layer ghost overlapping, three interleaved data sets are sampled with a skip size of N =11 along vertical direction in k-space. Figure 4.7 (c) is the deghosted image yielded by Simplified S P E E D based on a double-layer model, with an acceleration factor of N / 3 = 11/3 ~ 3.7. Figure 4.7 (c) is measured to have a T R E of 0.12. Figure 4.7 (d) is the result reconstructed by Half-Fourier Simplified S P E E D based on a double-layer model from 17% k-space data, which are acquired through the same sampling scheme used in obtaining Figure 4.7 (c) but within 60% k-space along vertical direction, corresponding to Half-Fourier coverage. Compared with the gold standard shown in Figure 4.7 (a), the deghosted image shown in Figure 4.7 (d) has the errors comparable to the noise level of the gold standard quantified by T R E , but with a high compression factor close to 6. It implies that we may reduce the scan time from 8.5 minutes (or 2 5 6 x 1 0 0 0 x 2 ms) with full k-space sampling to 1.4 minutes for this phantom scan. Figure 4.7 (d) is measured with E q . (4.23) to have a T R E of 0.18. Figure 4.7 (e-h) are the corresponding M I P images of the 25 slices reconstructed in the same way as Figure 4.7 (a-d). 106 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.7: Results of phantom scan reconstructed from three interleaved data sets with a skip size of N = 11 along vertical direction, (a) is the gold standard, which is the same as Figure 4.6 (a), (b) is ghosted image reconstructed by direct 2 D F T from one of the three skipped k-space data sets, (c) is the deghosted image yielded by Simplified S P E E D based on a double-layer model with an acceleration factor of 11/3 ~ 3.7. (d) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a double-layer model with an acceleration factor close to 6. (e-h) are corresponding M I P images of the 25 slices reconstructed in the same way as (a-d). Both (c) and (d) have the errors comparable to the noise level of the gold standard (a), and their corresponding M I P images (g) and (h) have quite good image qualities, (c) has a T R E of 0.18 while (d) has a T R E of 0.12. W i t h a lower window and level setting (-100, 100), the contrasts i n (a-d) are reduced to allow more information to be observed as shown in (i-1), respectively. 107 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.8 shows the results of the phantom scan reconstructed by Simplified S P E E D with two P E directions based on a double-layer model. In this case, k-space is sampled into three interleaved and relatively shifted data sets, each sampled with a skip size of N = 5 along vertical direction and with a skip size of M = 2 along horizontal direction. Figure 4.8 (a) is one of the three ghosted images reconstructed by inverse F T from the three interleaved data sets. Figure 4.8 (b) is the deghosted image reconstructed by Simplified S P E E D based on a double-layer model, with an acceleration factor of M x N / 3 = 10/3 ~ 3.3. Figure 4.8 (c) is the deghosted image reconstructed by Half-Fourier Simplified S P E E D in a way similar to Figure 4.8 (b) but from three interleaved data sets sampled within 60% k-space coverage along vertical direction, corresponding to only 19% k-space data or an acceleration factor of 5.5. The acceleration factor of 5.5 implies the possible reduction of scan time from 8.5 minutes to 1.5 minutes. Figure 4.8 (d-f) are corresponding M I P images of the 25 slices reconstructed in the same way as Figure 4.8 (a-c). Compared with the gold standard shown i n Figure 4.7 (e), the reconstructed M I P images shown in Figure 4.8 (e) and (f) have the errors comparable to the noise level of the gold standard quantified by T R E . Figure 4.8 (b) and (c) are measured with E q . (4.23) to have T R E values of 0.13 and 0.20, respectively. 108 Chapter 4 Simplified SPEED with Applications (d) in MRA (e) i * > "\ \ • »* *« (f) * V '* •» V * '%*•»» % *• V , * * (g) (h) i . V I * * (i) Figure 4.8: Results of phantom scan reconstructed from three interleaved data sets with a skip size of N = 5 along vertical direction and a skip size of M = 2 along horizontal direction. Figure 4.7 (a) and (e) are respectively the gold standard and the corresponding M I P image, (a) is one of the three ghosted images, (b) is the deghosted image reconstructed by Simplified S P E E D based on a double-layer model with an acceleration factor of 10/3. (c) is the deghosted image yielded by combination of Simplified S P E E D and Half-Fourier imaging based on a double-layer model with an acceleration factor of 5.5. (d-f) are corresponding M I P images of the 25 slices reconstructed in the same way as (a-c). The reconstructed images (b-c) and their corresponding M I P images (e-f) are all have reasonably good image qualities. Compared with the gold standard shown in Figure 4.7 (a), (b) and (c) are respectively measured to have T R E values of 0.13 and 0.20. For easy visualization, the contrasts in (a), (b) and (c) are reduced with a lower window and level setting (-100, 100) as shown in (g), (h) and (i), respectively. 109 Chapter 4 Simplified SPEED with Applications in MRA Figures 4.9-4.12 show the results from the in vivo P C gradient-echo head scan. Shown in Figure 4.9 (a) is one of the 25 slices reconstructed from full k-space data and is considered as the gold standard. Figure 4.9 (b) is what Figure 4.9 (a) becomes after being superimposed by noise points randomly taken from its background, resulting in a T R E of 0.67. (c) Figure 4.9: (a) is the gold standard image with full k-space sampling, (b) is what (a) becomes after being superimposed by noise points randomly taken from its background, resulting in a T R E of 0.67. To reveal more information, the contrast in (a) and (b) are reduced with a lower window and level setting (-25, 25) as shown in (c) and (d), respectively. 110 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.10 illustrates the results of the P C head scan reconstructed by Simplified S P E E D with one P E direction from two or three interleaved data sets, each sampled with a skip size of N = 5 and a different relative shift in P E along the vertical direction within 60% of k-space, corresponding to Half-Fourier coverage. Figure 4.10 (a) is one of 25 slices reconstructed from full k-space data and is considered as the gold standard. Figure 4.10 (a) is identical to Figure 4.9 (a). Figure 4.10 (b) is a ghosted image reconstructed by one of the interleaved k-space data sets sampled with a skip size of N = 5 within 60% of kspace. Figure 4.10 (c) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a single-layer model, with an acceleration factor of 4. Figure 4.10 (c) is measured to have T R E of 0.62, which is comparable to the noise level in Figure 4.9 (b). Figure 4.10 (d) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a double-layer model with an acceleration factor of 2.8. Figure 4.10 (d) is measured to have T R E of 0.51 by using E q . (4.23). Figure 4.10 (e-h) are corresponding axial M I P images of 25 slices reconstructed in a way similar to Figure 4.10 (a-d). Figure 4.10 (i-1) are the M I P images that are similar to Figure 4.10 (e-h) but at a different angle. Compared with the gold standard, the deghosted images and corresponding M I P images have quite good image qualities. T o allow more information to be observed, the contrast in Figure 4.10 (a-d) are reduced with a lower window and level setting (-25, 25) as shown in (m-p), respectively 111 Chapter 4 Simplified SPEED (m) (n) with Applications in MRA fo) (p) Figure 4.10: Results of the in vivo P C head scan reconstructed from various interleaved data sets, each sampled with a skip size of N = 5 and a different relative P E shift along vertical direction within 60% of k-space. (a) is the gold standard, which is identical to Figure 4.9 (a), (b) is the ghosted image reconstructed by direct 2 D F T from one of the interleaved k-space data sets, (c) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a single-layer model with an acceleration factor of 4. (d) is the deghosted image based on a double-layer model with an acceleration factor of 2.8. (c) and (d) are respectively measured to have T R E values of 0.62 and 0.51. (e-h) are corresponding axial M I P images of the 25 slices reconstructed in the same way as (a-d). (i-1) are the M I P images that are similar to (e-h) but at a different view angle. Compared with (a), (c) and (d) are both reasonably good while scan time has been highly reduced, (m-p) are the same as (a-d) but with a lower scale setting (-25, 25). 112 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.11 is similar to Figure 4.10 in that it shows the results reconstructed by Simplified S P E E D with one P E direction, but the former uses a skip size of N = 7 in P E . A s shown in Figure 4.11 (a), although seven-fold aliasing ghosts are seen in the ghosted images, most pixels can still be adequately modeled with a single layer because of the sparse signal distribution. Figure 4.11 (b) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a single-layer model, with an acceleration factor close to 6. It implies that we may reduce the scan time from 3.6 minutes with full k-space sampling to 0.6 minutes for this P C gradient-echo head scan. Figure 4.11 (c) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a double-layer model, with an acceleration factor of nearly 4. W i t h E q . (4.23), Figures 4.11 (b) and 4.11 (c) are respectively measured to have T R E values of 0.68 and 0.59. Figure 4.11 (d-f) are corresponding axial M I P images of the 25 slices reconstructed in a way similar to Figure 4.11 (a-c). Figure 4.11 (g-i) are the M I P images that are similar to Figure 4.11 (df) but at a different view angle. W i t h scan time being significantly reduced, the images reconstructed by Half-Fourier Simplified S P E E D and corresponding M I P images shown in Figure 4.11 are reasonably good compared with the gold standard Figure 4.10. 113 shown i n Chapter 4 Simplified SPEED with Applications in MRA (b) (a) (c) *i • •* * 1 <<i) (0 (e) ^ ,i i I j *i * J *< J *. 4 "j * • ^'4 f 1 (h) (g) 1 i f V ) H »' / (») f- m,. 1 (j) . •»4 (k) Figure 4.11: Results of the in vivo P C head scan reconstructed from various interleaved data sets, each sampled with a skip size of N = 7 and a different relative P E shift along vertical direction within 60% of k-space. Figure 4.10 (a) shows the gold standard; Figure 4.10 (e) and (i) are the corresponding M I P images of the gold standard at two different view angles, (a) is one of the ghosted images, (b) is the deghosted image yielded by Half-Fourier Simplified S P E E D based on a single-layer model with an acceleration factor close to 6. (c) is the deghosted image yielded based on a double-layer model with an acceleration factor close to 4. (b) has a T R E of 0.68 while (c) has a T R E of 0.59. (d-f) are corresponding axial M I P images of the 25 slices reconstructed in the same way as (ac). (g-i) are the M I P images that are similar to (d-f) but at a different view angle. Both (b) and (c) have the errors comparable to the noise level of the gold standard but with high compression factors. For easy visualization, the contrasts in (a-c) are reduced with a lower window and level setting (-25, 25) as shown in (j-1), respectively. 114 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.12 shows the results of the P C head scan reconstructed from various interleaved data sets by Simplified S P E E D with two P E directions, each sampled with a skip size of N = 3 along vertical direction and a skip size of M = 2 along horizontal direction within 60% of k-space in vertical direction. A s shown in Figure 4.12 (a), aliasing ghosts of various orders are seen along both vertical and horizontal directions. Figure 4.12 (b) is the deghosted image yielded by Fourier Simplified S P E E D based on a single-layer model with an acceleration factor of 5, resulting in a reduction of scan time from 3.6 minutes to 0.7 minutes for the P C head scan. Figure 4.12 (b) is measured to yield a T R E of 0.67. Figure 4.12 (c) is the deghosted image yielded by Fourier Simplified S P E E D based on a double-layer model with an acceleration factor of 3.3. B y using E q . (4.23), Figure 4.12 (c) is measured to have a T R E of 0.56. Figure 4.12 (d-f) are corresponding M I P images of the 25 slices reconstructed in a way similar to Figure 4.12 (a-c). Figure 4.12 (g-h) are the M I P images that are similar to Figure 4.12 (d-f) but at a different angle. It can be seen that Simplified S P E E D can compress scan time along two orthogonal directions to achieve rather good results, illustrating the flexibility of the acceleration algorithm. 115 Chapter 4 Simplified SPEED with Applications in MRA Figure 4.12: Results of the in vivo P C head scan reconstructed from various interleaved and relatively shifted data sets, each sampled respectively with a skip size of N = 3 along vertical direction and a skip size of M = 2 along horizontal direction within 60% of kspace in vertical direction. Figure 4.10 shows the gold standard and the corresponding M I P images, (a) is one of the ghosted images, (b) is the deghosted image yielded by HalfFourier Simplified S P E E D based on a single-layer model with an acceleration factor of 5. (c) is the deghosted image yielded based on a double-layer model with an acceleration factor of 3.3. (b) has a T R E of 0.67 while (c) has a T R E of 0.56. (d-f) are corresponding axial M I P images of the 25 slices reconstructed in the same way as (a-c). (g-i) are the M I P images that are similar to (d-f) but at a different view angle. A l l the aliasing ghosts in (a) have been successfully unfolded and the results shown in (b) and (c) are both visually good, (j-1) are the same as (a-c) but with a lower scale setting (-25, 25). 116 Chapter 4 Simplified SPEED with Applications in MRA 4.4 Discussion Based on the sparseness of M R A data, S P E E D is simplified to enhance its efficiency. L i k e S P E E D , Simplified S P E E D samples k-space at selective P E steps with a certain skip size, and suppresses associated aliasing ghosts. However, while S P E E D uses differential operation to generate edge maps, Simplified S P E E D directly models the sparse ghosted images with a single-layer model, and thereby achieves deghosting with only two interleaved data sets. Moreover, Simplified S P E E D is not limited to this single-layer model; it can be extended to include additional layers of ghosts as ghost overlapping increases. In other words, analogous to Taylor expansion, the solution of Simplified S P E E D can be refined through a G L U E series by adding more layers into the model. If all N terms are included in the G L U E series, a Simplified S P E E D solution w i l l converge to a perfect solution. A nice feature of a G L U E series is that it always begins with the most dominant layer and ends with the least important one. It follows that a good approximation can be rapidly reached with the first few dominant layers. In real time M R A , additional data can be progressively acquired to refine the deghosting solution until a satisfactory result is achieved. Since the quality of Simplified S P E E D solution is monitored by its residual, the automatically. above real-time processing can be accomplished In an extremely challenging case where the vascular structure is too complicated to be compressed, all N terms of G L U E series must be included and full kspace data are sampled, resulting in a perfect deghosted image identical to the gold 117 Chapter 4 Simplified SPEED with Applications in MRA standard. Fortunately, this situation rarely occurs because M R A data are generally very sparse, therefore Simplified S P E E D can always be used to reduce scan time in most cases. Generally speaking, Simplified S P E E D can be considered as a scan-time reduction technique for sparse data. Relying on a principle distinctive from those used in other fast M R A methods, Simplified S P E E D is an independent fast M R A approach that works even with one single coil. It is well known that undersampled P R is a non-Cartesian fast M R A method based on conventional P R [74,75]. Compared with undersampled P R , Simplified S P E E D is a Cartesian fast M R A technique. Unlike T R I C K S [76], which uses temporal information, or parallel imaging [48,49,79-81], which uses multiple coils to reduce scan time, Simplified S P E E D uses solely phase information that is available even with a single coil. It can be combined with other existing fast imaging methods to achieve more efficiency. In this work, a combination of Simplified S P E E D and Half-Fourier imaging is proposed and its feasibility is demonstrated. Furthermore, by taking advantage of the significant correlation in k-space and time within dynamic M R A series, Simplified S P E E D may be extended to four-dimensional k-t space to reduce more scan time in a way similar to P R T R I C K S [77,78]. Where multiple coils are available, it is also possible to combine Simplified S P E E D with sensitivity encoding to accelerate M R A . A s a general method, Simplified S P E E D can be extended to accelerate scan time along two or more P E directions for improved performance. In this work, Simplified S P E E D 118 Chapter 4 Simplified SPEED with Applications in MRA with either one or two P E directions has been demonstrated. Compared with Simplified S P E E D with one P E direction, Simplified S P E E D accelerating M R A along two P E directions has more potential in reducing scan-time. Although Simplified S P E E D is demonstrated with only phantom data and in vivo P C M R A data in this work, it also has potential applications in Contrast Enhanced M R A ( C E - M R A ) [73,85,86]. In C E - M R A , both precontrast and postcontrast data sets are sampled for improving vasculature visualization. Once postcontrast data set is obtained, the precontrast data set is subtracted from it. Thus, all the background structures on the precontrast images are subtracted out to yield a contrast-enhanced image, similar to P C MRA. Such a subtraction can be combined with Simplified S P E E D . Unlike the traditional subtraction technique described above, precontrast and postcontrast data sets can be subtracted first and then processed by Simplified S P E E D , resulting in a final contrast-enhanced image. In summary, based on the sparse nature of M R A data, this work successfully simplified S P E E D and has demonstrated its potential to accelerate M R A with a single coil by a factor close to 6. 119 Chapter 5 Accelerating MRI by SPEED and its Extensions Chapter 5 Accelerating MRI by Skipped Phase Encoding and Edge Deghosting (SPEED) and its Extensions 5.1 Introduction In the previous chapter, the fast imaging method of Skipped Phase Encoding and Edge Deghosting ( S P E E D ) [3] was simplified to enhance its efficiency by taking advantage of the sparseness of M R A data. Simplified S P E E D is proposed to compress scan time for sparse data such as M R A , while the original S P E E D can be used to compress scan time for typical M R I data. In this chapter, the original S P E E D w i l l be reviewed and its two selected extensions (i.e. S P E E D with two P E directions i n 3D and Half-Fourier-SPEED) w i l l be studied. 5.2 Methods 5.2.1 Review of the Original SPEED S P E E D is a recently developed fast imaging method [3]. A s shown i n Figure 5.1, S P E E D partially samples k-space at selective P E steps into three sparse interleaved data sets: S,(k), S (k) and S (k). Each data set has the same skip size o f N but a different 2 3 120 Chapter 5 Accelerating MRI by SPEED and its Extensions relative shift in P E . These data sets are then Fourier transformed separately into three ghosted images I^r), I (r) and I (r), each associated with N-fold aliasing ghosts resulted 2 from the N-step 3 skipped undersampling. Obviously, each ghosted image is a superposition of the N aliasing ghosts and as such has ghost overlapping up to N layers. Unlike sensitivity encoding, S P E E D does not use multi-coil sensitivities to unfold the signal superposition. Rather, it high-pass filters the ghosted images into three ghosted edge maps: E,(r), E (r) and E (r), i n which the number of ghost overlapping layers is 2 significantly reduced. 3 B y modeling the three ghosted edge maps with a double-layer structure, a deghosted edge map E (r) can be solved through an analytical algorithm 0 based on least-square-error ( L S E ) minimization, i n the same way as it is done i n Simplified S P E E D . Finally, the deghosted edge map E (r) can be inverse-filtered into a 0 final deghosted image I (r). T o quantify the efficiency of S P E E D , the skip size of N in 0 P E is hereafter referred to as an undersampling factor of N . The principle o f S P E E D can be further illustrated i n the following specific case, in which scan time is reduced by a factor of 3/5. 121 Chapter 5 Accelerating MRI by SPEED and its Extensions Figure 5.1: K-space sampling scheme and data flow o f 2 D S P E E D [3]. 2 D full k-space S (k) is sampled into three interleaved sparse k-space data sets: S[(k), S (k) and S (k), each with a skip size of N and a relative shift in P E . After inverse 2 D F T followed by a differential filtering, S,(k), S (k), and S (k) are firstly turned into three ghosted images: I^r), I (r) and I (r), and then three ghosted edge maps: E,(r), E (r) and E (r). Finally, they are deghosted into a ghost-free edge map Eo(r), which is subsequently inversefiltered to yield a deghosted image In(r). 0 2 2 2 3 3 3 2 122 3 Chapter 5 Accelerating MRI by SPEED and its Extensions Three interleaved data sets (Si, S2, S3) are obtained by sampling k space at every 5 PE th step with three different relative shifts: dl, dl and d3, which are arbitrarily chosen from 0, 1, 2, 3, and 4. The data sets are reconstructed into three ghosted edge maps (Ei, E2, E 3 ) after a standard FT followed by a differential operation along PE direction. The differential operation is just equivalent to a high-pass filtering in k-space. If PE is performed in the y-direction, the y-directional high-pass filtering operator is defined in image space as E(x,y) = I(x,y + l)-I(x,y) (5.1) where I{x,y) stands for a ghosted image, and E{x,y) denotes the edge map of I{x,y). In a typical ghosted edge map, ghost overlapping can simply be modeled with a double-layer structure [3]. Thus, the three ghosted edge maps (Ei, E2, E 3 ) can be mathematically expressed at each pixel as dl E d2 E = d\ n\ P + dl n2 G P = d2 n\ P (-) G 5 + d2 n2 G P (-) G 5 E Z=Pd\ nX+P£ n2 G 2 3 (5-4) G d where Ed stands for the ghosted edge map obtained from the data set with a relative shift of d in PE; G„i and G„2 are the two dominant ghost layers; and Pj is a ghost phasor 1 123 Chapter 5 Accelerating MRI by SPEED and its Extensions ilnd known to have a form of (e )" , where n is the order of ghost depending on its relative 5 location and must be chosen from 0 to 4. Although ghost order index (nl,n2) at each pixel remain to be fully determined, they are known to be two different integers chosen from 0, 1, 2, 3, and 4. Therefore, we can write Eqs. (5.2-5.4) into a matrix form described by E q . (5.5), in a fashion similar to Simplified S P E E D discussed in Chapter 4. E = MG where E = (E \, Eai, E 3) is known vector, G = (G i,G 2) T d (5.5) T d n n is an unknown vector yet to be found, and M is a 3x2 matrix containing parameters (nl,n2) given by >«2 >nl d\ mi M = dl d3 l 1 d\ >n2 dl -.nl dZ (5.6) Similar to E q . (4.14), E q . (5.5) describes an over-determined linear system, in which there are three equations and two unknowns. Therefore, a L S E solution vector GLSE can be solved for each possible pair of (nl, n2) as follows: LSE G ~ L (5.7) E where L is referred to as a "pseudo-inverse" of matrix M and is a 2x3 matrix defined as M (5.8) 124 Chapter 5 Accelerating MRI by SPEED and its Extensions The L S E for this solution GLSE can be expressed as LSE = energy[E - MG LSE ] (5.9) Because of the symmetric positions of (nl, n2) in Eqs.(5.2-5.4), the total number of possible (nl, n2) pairs is / V ( . / V - i ) / 2 = 10 for N = 5. Therefore, through 10 trials, the correct (nl, n2) pair and (Gni, G i) solution can be found through L S E minimization. A s n long as N is a prime number, the L S E minimum w i l l not have the ghost phase degeneracy [83,84], where ghosts of different orders have the same phase. According to the associated (nl, nl) values, the resolved ghosts (G„i, G i) can be sorted out to yield five n separate edge ghosts. The five separate edge ghosts are spatially registered and averaged to yield a single deghosted edge map with reduced noise and artifacts. The above deghosting process is accomplished in the same way as it is done in Simplified S P E E D . Finally, the deghosted edge map is inverse-filtered to yield a final deghosted image with an undersampling factor of 5/3 - 1 . 7 . According to the properties of discrete F T [71], the Fourier relationship between an image I(x,y) and its k-space data S(k , x k ) can be y mathematically expressed as ) (5.10) where N and N are respectively k-space matrix dimensions along the JC and y directions. x y 125 Chapter 5 Accelerating MRI by SPEED and its Extensions W i t h E q . (5.10) at hand, the y-differential operator defined in E q . (5.1) is found to have its k-space counterpart as a multiplicative filter given by ky_ F(k ) y =e" 2 % y -1 (5.11) The inverse-filtering can therefore be interpreted as dividing the Fourier transform of the deghosted edge map by the k-space filter defined in E q . (5.11). However, the operation should be done cautiously because streak artifacts may arise when dividing the numbers close to zero, particularly at k-space data points near k = 0. y The direct division would result in streak artifacts in the final image along y direction [3]. T o avoid such artifacts, extra true data (e.g. a band of 8-32 k-space lines) are symmetrically sampled near the kspace center and then used to replace the inverse-filtered data. In this work, central 32 kspace data lines are fully sampled and used to aid the inverse-filtering. Although additional scan time is spent on sampling extra data, this is still relatively insignificant. Furthermore, this penalty can be reduced by using the shared views already existed in the three interleaved data sets: S i , S2, and S3 [3]. Figure 5.2 shows a flowchart of the S P E E D algorithm used in the example. Similar to Simplified S P E E D , original S P E E D not only provides a deghosting solution, but also outputs a residual defined as the square root of the L S E in Eq.(5.9). It is the residual that is minimized in deghosting and indicates how adequate the double-layer 126 Chapter 5 Accelerating MRI by SPEED and its Extensions model is. Therefore, the quality of the S P E E D solution can be monitored by using the residual. Compared with Simplified S P E E D using sparseness of M R A data, S P E E D accelerates M R I by taking advantage of the sparse nature of image edge. The differences between SPEED and Simplified reconstruction. SPEED are present in both data acquisition and image In data acquisition, S P E E D samples k-space in the same way as Simplified S P E E D except that S P E E D acquires additional central k-space data to aid the inverse filtering. In image reconstruction, both S P E E D and Simplified S P E E D use the same deghosting technique except that S P E E D uses high-pass filtering prior to deghosting and inverse filtering after a deghosted edge map is obtained. S P E E D is efficient and robust; however, S P E E D has been so far demonstrated only in 2D. It is desirable to extend the principle of S P E E D from 2 D to 3D for improved performance. 127 Chapter 5 Accelerating MRI by SPEED and its Extensions Input 3 i n t e r l e a v e d datasets w i t h N =5 R e c o n s t r u c t i n t o 3 g h o s t e d edge m a p s E^E^, & by h i g h - p a s s filtering & F T Y i e l d the L S E s o l u t i o n o f 2 d o m i n a t i n g ghost l a y e r s G „ , & G „ , t h r o u g h 10 trials o f n l & n 2 a t e a c h p i x e l I I S o r t out 5 separate edge ghostsG,, a c c o r d i n g to a s s o c i a t e d n 1 & n 2 R e g i s t e r the 5 edge ghosts G „ a n d a v e r a g e I I t h e m i n t o a f i n a l d e g h o s t e d edge m a p Y i e l d a final d e g h o s t e d i m a g e !„• by inverse filtering ' O u t p u t the f i n a l i m a g e I„> Figure 5.2: Flowchart of the SPEED with one PE direction used to achieve an undersampling factor of 1.7 with a single coil [3]. 128 Chapter 5 Accelerating MRI by SPEED and its Extensions 5.2.2 E x t e n s i o n o f S P E E D to T w o P E directions Three-dimensional M R imaging is a popular area; however, imaging time is often a concern. In 3 D , P E is typically applied along two orthogonal directions. S P E E D can therefore accelerate data acquisition along the two directions with a reduced chance of having edge ghost overlapping, which allows S P E E D to achieve more flexible and more efficient reduction of scan time in 3 D . A s shown in Figure 5.3, S P E E D with two P E directions samples full k-space S (k) at selective P E steps into three interleaved data sets: 0 S^k), S (k), and S (k). Each data set has a skip size of N along one P E direction and has 2 3 another skip size of M along the other P E direction. After an inverse 2 D F T on each dataset, S,(k), S (k), and S (k) are reconstructed separately into three ghosted images: 2 3 Ij(r), I (r), and I (r). They are then turned into three ghosted edge maps — E,(r), E (r) 2 3 2 and E (r) — by using a differential filter. The differential operation can be performed 3 along either P E direction or both. In this work, the differential filtering along one direction is found to be capable of obtaining a sparse image edge, and as such is used i n 3D for simplicity. W i t h an analytical algorithm based on a L S E solution, a deghosted edge map E (r) is obtained and is subsequently inverse-filtered into a final deghosted 0 image I (r). The algorithm of S P E E D with two P E directions can be further illustrated in 0 the following specific example, in which scan time is reduced by an undersampling factor of 2. 129 Chapter 5 Accelerating MRI by SPEED and its Extensions Figure 5.3: K-space sampling scheme and data flow of S P E E D with two P E directions. F u l l k-space S (k) is sampled into three interleaved sparse k-space data sets: S[(k), S (k) and S (k), each with a skip size of N along one P E direction and a different skip size o f M along the other P E direction. After inverse 2 D F T followed by a differential filtering, S,(k), S (k) and S (k) are first turned into three ghosted images: I,(r), I (r) and I (r), and then three ghosted edge maps: E,(r), E (r) and E (r). Finally, they are deghosted into a ghost-free edge map En(r), which is subsequently inverse-filtered into a deghosted image Io(r). 0 2 3 2 3 2 2 3 130 3 Chapter 5 Accelerating MRI by SPEED and its Extensions Three interleaved and relatively shifted data sets ( S i , S , S3) are obtained by sampling k2 space with a skip size of N = 3 in one P E and with a skip size of M = 2 in another P E . Each data set has a unique relative shift pair (d, s) where d is chosen from 1, 2, and 3 and 5 is chosen from 0 and 1. The three data sets are then Fourier transformed into three ghosted images (Ii, I , I3), with ghost overlapping up to six layers. T o reduce the ghost 2 overlapping, the y-differential filter defined in E q . (5.1) is used to turn the ghosted images into three sparse ghosted edge maps ( E i , E , E 3 ) . 2 B y modeling them with a double-layer structure, the three ghosted edge maps can be expressed as follows: P-l = PdlQslGnlml 2 = Pd2Q?2 n\m\ E 3 E where G imi and n G = PdlQ?l n\m\ G G 2m2 n to have a form of (e Pd\Q7\Gn2ml (5.12) + Pd2Q?2 n2m2 (5.13) d3Q?* n2m2 (5.14) + G + P G are the two dominant layers of ghosts; P is a ghost phasor known d / , where n is the order of ghost depending on its relative 3 location and must be chosen from 0,1, and 2; and Q™ is another ghost phasor known to I2K s have a form of (e 2 ), m where m is the order of ghost depending on its relative location and must be chosen from 0 and 1. Since the values of nl, n2 must be chosen from 0 to 2 and the values of ml, ml must be chosen from 0 and 1, all possible pairs of (nl, n2, ml, m2) are known with a limited total 131 Chapter 5 Accelerating number of choices. MRI by SPEED and its Extensions B y trying all possible pairs, six separate edge ghosts of different orders can be solved by minimizing the L S E on a pixel-by-pixel basis by repeating Eqs. (5.5-5.9). The six edge ghosts are then spatially registered and averaged into a single deghosted edge map En. In a fashion similar to S P E E D with one P E direction, the deghosted edge map En can be inverse-filtered to yield a final deghosted image In with an undersampling factor of N x M / 3 = 2. In this work, central 32 k-space data lines i n y direction are fully sampled and used to aid the inverse filtering. 5.2.3 C o m b i n a t i o n o f S P E E D and H a l f - F o u r i e r I m a g i n g A s discussed previously, S P E E D acquires data through skipped sampling of k-space while eliminating the resulting aliasing ghosts through edge deghosting. L i k e Simplified S P E E D , S P E E D is an independent method and is able to achieve improved performance by combining itself with other existing methods. This work combines S P E E D [3] and Half-Fourier imaging [39] into a new fast M R I method named Half-Fourier-SPEED. It is simple to implement and reduces scan time further by a factor of nearly 2 without the need of additional hardware. A s illustrated by the data flow shown in Figure 5.4, Half-Fourier-SPEED can be divided into three major steps. 132 Chapter 5 Accelerating MRI by SPEED and its Extensions (1) First, three interleaved and relatively shifted data sets are sampled with a skip size of N within slightly more than half k-space coverage, resulting in an undersampling factor of nearly 2 N / 3 . In practice, data are usually sampled to span 60% of k-space, corresponding to typical Half-Fourier coverage. The three data sets are then Fourier transformed separately into three ghosted images, each associated with N fold aliasing ghosts and truncation artifacts due to k-space undersampling. (2) Second, to suppress the associated aliasing ghosts, the three ghosted images are processed by S P E E D to yield a deghosted image. Fully sampled central k-space data (e.g. 32 out of 256 lines) are used to aid the inverse filtering and to measure a phase estimate for Half-Fourier reconstruction used in the next step. (3) Finally, to reduce the truncation artifacts resulted from Half-Fourier coverage, the deghosted image is processed by Half-Fourier reconstruction to yield a final image. More specifically, Half-Fourier reconstruction applies a Harming roll-off filter to the deghosted image and then recovers a desired image by phase correction with a phase estimate [39]. The phase estimate is generally obtained by using the phase of the low-resolution image reconstructed from fully sampled central k-space data. For IR M R I , the phase of the lowresolution image cannot be directly used as the phase estimate, but a special phase unwrapping method can be used to obtain the phase estimate [2,72]. In this work, the improved W L S method [2] is used to calculate the phase estimate for IR images. 133 Chapter 5 Accelerating MRI by SPEED and its Extensions S,(k) S (k)J : |s,(k) I,(r)| \l,(r) E,(r) E (r) E,(r) ( E»(r) v Halt-Fourier I roeonstrucuoti : ^HJ I„(r) • Figure 5.4: K-space sampling scheme and data flow of Half-Fourier-SPEED. Slightly more than half k-space S (k) is sampled into three interleaved k-space data sets: S,(k), S (k), and S (k), each with a skip size of N and a relative shift in P E . After inverse 2 D F T followed by a differential filtering, Sj(k), S (k) and S (k) are first turned into three ghosted images: I,(r), I (r), and I (r), and then three ghosted edge maps: E,(r), E (r), and E (r). They are deghosted into a ghost-free edge map Eo(r), which is subsequently inverse-filtered to yield a deghosted image Io(r). Finally, In(r) is processed by HalfFourier reconstruction to yield a final image Ifinai(r)0 2 3 2 2 3 3 2 3 134 Chapter 5 Accelerating MRI by SPEED and its Extensions 5.2.4 G e n e r a l i z a t i o n o f S P E E D Although S P E E D is discussed and demonstrated with a double-layer structure in this chapter, the application of S P E E D is not limited to the double-layer model. Based on Ghost Layer Unfolding Expansion ( G L U E ) introduced in Chapter 4, the sparse ghosted edge maps can be expanded through a G L U E series by adding progressively more layers of ghosts for more accuracy, analogous to the sparse ghosted M R A images in Simplified SPEED. For example, when the double-layer structure is not adequate to model the ghosted edge maps ( E i , E2, E ) , additional interleaved data sets and more layers in the 3 ghost model can be used for better image quality until a satisfactory result with a tolerable residual is achieved. In real-time M R I , the solution of S P E E D should be firstly solved based on a single-layer model. Then, the solution is refined through a G L U E series by adding more layers into the model until S P E E D yields a satisfactory deghosted edge map, analogous to the generalized Simplified S P E E D . 135 Chapter 5 Accelerating MRI by SPEED and its Extensions 5.3 Results Various versions of S P E E D have been tested with in vivo data and demonstrated below in Figures. 5.5 and 5.6 with a spin-echo transverse head scan on a clinical 1.5 T scanner. Figure 5.5 shows the results reconstructed by S P E E D with one P E direction from the head scan. For comparison, Figure 5.5 (a) shows the "gold standard" image reconstructed from full k-space data. Figure 5.5 (b) is the Fourier transform of one of the three interleaved datasets. A s shown in Figure 5.5 (b), aliasing ghosts of various orders are seen along the horizontal P E direction and most overlapping pixels have more than two ghost layers. Figure 5.5 (c) is one of the three ghosted edge maps, which are obtained after a differential filtering is applied to Figure 5.5 (b). Because the number of overlapping ghost layers is effectively reduced by the differential filtering, most ghost overlapping pixels can be adequately modeled with a double-layer structure. Figure 5.5 (d) is the deghosted edge map reconstructed by S P E E D with one P E direction. Figure 5.5 (e) is the final deghosted image obtained after the deghosted edge map is inverse-filtered. Figure 5.5 (e) is measured with E q . (4.23) to yield a T R E of 0.07. The difference between Figures 5.5 (e) and 5.5 (a) is hardly noticeable although the scan time has been reduced by a factor of 0.6. Figure 5.5 (f) is a residual map, which is the square root of the L S E in Eq.(5.8). The pixel intensities in the residual map are much less than those in Figure 5.5 (c) and therefore demonstrate an adequate double-layer modeling. 136 Chapter 5 Accelerating MRI by SPEED and its Extensions Figure 5.5: Results from S P E E D with one P E direction and a standard reconstruction for a spin-echo transverse head scan, (a) is the gold standard with full k-space sampling, (b) is one of the three ghosted images with P E skip size N = 5, having up to five layers of ghost overlapping, (c) is the corresponding edge map after high-pass filtering, with ghost overlapping layers being reduced to be only 2 in most pixels, (d) is the edge map after deghosting. (e) is the final deghosted image by inverse filtering (d). (e) is measured to have a T R E of 0.07. W i t h an undersampling factor of 5/3 ~ 1.7, the final image (e) looks comparable to the gold standard (a), (f) is a residual map for quality control. 137 Chapter 5 Accelerating MRI by SPEED and its Extensions Figure 5.6 shows the results reconstructed by S P E E D with two P E directions from the head scan. For comparison, Figure 5.6 (a) shows the "gold reconstructed from full k-space data by a standard 2 D F T . standard" image Figure 5.6 (b) is one of the three ghosted images obtained with a skip size N = 3 along vertical direction and another skip size M = 2 along horizontal direction. A s shown in Figure 5.6 (b), aliasing ghosts of various orders are seen along both directions and ghost overlapping is up to six layers. Figure 5.6 (c) is the corresponding edge map obtained after a high-pass filtering or simple differential operation is applied to Figure 5.6 (b). A s shown in Figure 5.6 (c), the number of ghost overlapping layers is reduced to be mostly only 2 through the high-pass filtering. Figure 5.6 (d) is the final deghosted image reconstructed by S P E E D with two P E directions with an undersampling factor of 2. The T R E of the deghosted image is measured to be 0.09. Compared with Figure 5.6 (a), the reconstructed image by S P E E D with two P E directions appears comparable to the gold standard image. 138 Chapter 5 Accelerating MRI by SPEED and its Extensions (c) (d) Figure 5.6: Results from S P E E D with two P E directions and a standard reconstruction for a spin-echo transverse head scan, (a) is the gold standard with full k-space sampling, (b) is one of the three ghosted images with a skip size N = 3 along vertical direction and another skip size M = 2 along horizontal direction, (c) is the corresponding edge map. (d) is the final deghosted image, which is comparable to the gold standard (a), (d) is measured to have a T R E of 0.09. In this case, the undersampling factor is 6/3 = 2. 139 Chapter 5 Accelerating MRI by SPEED and its Extensions Half-Fourier-SPEED has been successfully tested with in vivo data and is demonstrated below with a spin-echo sagittal knee scan and an IR coronal head scan from a clinical 1.5 T scanner. Around 40% of a 256x256 k-space matrix is sampled as follows: 64 lines are consecutively sampled near the k-space center. Then, half k-space including the 64 central lines is sparsely sampled into three data sets with a skip size of 5 (N=5) and three different relative shift sizes — 1,2, and 3 — in P E . Figure 5.7 (a-b) show reconstructed images from the sagittal knee scan. Figure 5.7 (b) is the result of Half-Fourier-SPEED reconstructed from about 40% k-space data, and is measured to have a T R E of 0.03. The result shown in Figure 5.7 (b) is comparable to the "gold-standard" reconstructed from 62.5% k-space data by conventional Half-Fourier imaging shown in Figure 5.7 (a). Figure 5.7 (c-d) show reconstructed images from the IR coronal head scan. Figure 5.7 (c) is the "gold-standard" reconstructed by conventional Half-Fourier imaging with 62.5% kspace data, corresponding to Half-Fourier coverage. Figure 5.7 (d) is the reconstructed image by Half-Fourier-SPEED with an undersampling factor of 2.5. The T R E of Figure 5.7 (d) is measured to be 0.05. The difference between Figures 5.7 (c) and 5.7 (d) is hardly noticeable. In both cases, P E is chosen to be along vertical direction. 140 Chapter 5 Accelerating MRI by SPEED and its Extensions Figure 5.7: Results from Half-Fourier-SPEED and conventional Half-Fourier imaging for a spin-echo sagittal knee scan and an IR coronal head scan. Although scan time has been further reduced by a factor of 3/5, the images (b) and (d) reconstructed by Half-FourierS P E E D have no visible difference from the gold standard (a) and (c), each obtained by conventional Half-Fourier reconstruction, (b) and (d) are respectively measured to have T R E values of 0.03 and 0.05, which are both very small. 141 Chapter 5 Accelerating MRI by SPEED and its Extensions 5.4 Summary In the first part of this chapter, the original S P E E D (i.e. S P E E D with one P E direction) is reviewed and its basic principle are illustrated with a specific 2 D example. The second part of the chapter substantially discusses two selected extensions of S P E E D : S P E E D with two P E directions in 3 D and Half-Fourier-SPEED. B y extending the principle of S P E E D from 2 D to 3 D , S P E E D can be used to achieve more flexible and more efficient reduction of scan time. B y combining S P E E D and Half-Fourier imaging, S P E E D is further improved to reduce scan time by another factor of nearly 2. L i k e S P E E D , HalfFourier-SPEED can also be extended to higher dimensions or be combined with other fast imaging methods for potentially reducing scan-time further. Finally, the principle of S P E E D is generalized based on a G L U E series, analogous to generalized Simplified S P E E D . Unlike Simplified S P E E D , which is used only for sparse data, S P E E D can be used to accelerate typical M R I data, which allows S P E E D to have many potential applications. 142 Chapter 6 Highly accelerated MRI with SPEED-ACE Chapter 6 Highly Accelerated MRI by Skipped Phase Encoding and Edge Deghosting with Array Coil Enhancement (SPEED-ACE) 6.1 Introduction A s discussed in Chapter 2, many fast imaging strategies have been proposed to reduce scan time. Some of the strategies physically manipulate spin dynamics to use available magnetization more efficiently [18, 35-38]. Others sparsely sample k-space and reconstruct a complete image through a non-standard reconstruction [3,39-52]. Amongst the strategies is parallel imaging, which reduces scan time by using multiple receiver coils simultaneously to acquire partial k-space data and then uses the information hidden in the partial data to reconstruct a complete image. M a n y parallel imaging methods have been proposed in the past few years [44-52]. Some of them [44,48,52] reconstruct a complete image in k-space, examples of which are S M A S H [48] and G R A P P A [52]. Others [45-47,49-51] do so in image space, examples of which are subencoding [47], S E N S E [49], S P A C E - R E P [50] and P I L S [51]. S E N S E stands for sensitivity encoding [49]. A s its name suggests, it uses coil sensitivity profiles to achieve spatial encoding and to accelerate M R I . The availability of multiple receiver coils on most modern M R I scanners has made S E N S E a popular approach to reduce scan time. 143 Chapter 6 Highly accelerated MRI with SPEED-ACE While parallel imaging reduces scan time by using multiple coils, S P E E D can accelerate M R I with only one single receiver coil. This shows both its efficiency and potential. In this chapter, I propose a new parallel imaging method that combines S P E E D and information available in multiple coils [5]. This method, named "Skipped Phase Encoding and Edge Deghosting with Array C o i l Enhancement ( S P E E D - A C E ) " , partially samples k-space by using multiple receiver coils in parallel and obtains a deghosted image through edge deghosting [5]. S P E E D - A C E extends the principle of S P E E D to improve parallel imaging with multiple coils, and the feasibility of S P E E D - A C E has been demonstrated with in vivo M R I data. 6.2 Methods 6.2.1 R e v i e w o f S E N S E To understand S P E E D - A C E , it is helpful to review the basics of S E N S E . In the Cartesian version of S E N S E , k-space is sampled with a skip size of N in P E with multiple receiver coils, as illustrated in Figure 6.1. The sparse k-space data set obtained from each coil can be directly reconstructed by F T into a sensitivity-weighted and ghosted image. Each ghosted image contains N aliasing ghosts along P E direction. T o yield a ghost-free image, the superposition must be unfolded. This process can be interpreted as separating the contributions from the N aliasing ghosts to the ghosted image at each pixel [49]. In S E N S E , the superposition is unfolded based on the fact that signal in each ghosted image 144 Chapter 6 Highly accelerated MRI with SPEED-ACE is likely weighted by a different coil sensitivity. Typically, to unfold N layers of ghosts, at least N coils are required to allow possible matrix inversion. S E N S E is a very promising method, but has its inherent limit: noise is amplified in an uncontrolled manner as acceleration factors increase [87,88]. Compared with S E N S E , S P E E D reduces scan time with only a single coil based on a different principle. •;,r.:~:;;.7 — S,(k)| I,(r) S (k)| |S (k) 3 3 Ur) Ur) | S.,(k) H u Figure 6.1: K-space sampling scheme and data flow of 2 D S E N S E with four receiver coils. 2 D full k-space is selectively sampled with a skip size of N by using four receiver coils in parallel. The sampled sparse data sets are reconstructed separately into four ghosted images: Ii(r), l2(r), l3(r), and Lt(r). Unfolding signal superposition yields a final deghosted image In(r). 145 Chapter 6 Highly accelerated MRI with SPEED-ACE 6.2.2 R e v i e w o f S P E E D A s illustrated i n Figure 5.1, S P E E D samples k-space at selective P E steps to obtain three interleaved data sets, each with the same skip size of N but a different relative shift in P E . These data sets are then Fourier transformed respectively into three ghosted images, each of which is associated with N-fold aliasing ghosts and as such has ghost overlapping up to N layers. S P E E D then high-pass filters the ghosted images into three ghosted edge maps ( E i , E2, E 3 ) , in which the overlapping ghost layers are significantly less than those in the ghosted images because of the sparse nature of edge maps. With an analytical deghosting algorithm based on a L S E minimization, a deghosted edge map is obtained and subsequently inverse-filtered into a final deghosted image, with an undersampling factor of N / 3 . A few central k-space lines (e.g.32 out of 256) are fully sampled to aid the inverse-filtering. In S P E E D , the skip size N must be a prime number to avoid ghost phase degeneracy [83,84]. S P E E D is an efficient and robust method; however, it is not taking full advantages of multiple coils available on modern M R I scanners: multiple coils can be used in S P E E D to improve S N R , but they are yet to be used to accelerate scan time. This thesis introduces S P E E D - A C E , which uses multiple coils undersampling factor and thereby is a step further than S P E E D . 146 to obtain a higher Chapter 6 Highly accelerated MRI with 6.2.3 SPEED-ACE SPEED-ACE In S P E E D - A C E , the first step is using multiple coils to sample k-space with a skip size of N steps [5]. A s illustrated in Figure 6.2, four sets of data are obtained in this way with four coils in parallel. S P E E D - A C E then Fourier transforms the sampled data into four sensitivity weighted and ghosted images: Ii(r), l2(r), h(r), and Li(r), each associated with N-fold aliasing ghosts resulting from the N-step skipped undersampling. Since coil sensitivity profile usually varies slowly, its spatial differential is often negligible compared with the edges of images. Based on this observation, S P E E D - A C E applies a differential filter to the ghosted images, resulting in a set of sparse ghosted edge maps: Ei(r), E2(r), E3(r), and E4(r). Given that most ghosted edge maps can be adequately modeled with two dominant ghost layers [3], the four ghosted edge maps can be deghosted into a ghost-free edge map Eo(r) through an analytical algorithm based on L S E minimization. Finally, the deghosted edge map Eo(r) is inverse-filtered into a final deghosted image In(r). In this way, S P E E D - A C E achieves an undersampling factor of N , which is three times greater than the undersampling factor of N/3 achieved by S P E E D using a double-layer model and the same skip size of N . It is noteworthy that S P E E D A C E not only uses skipped phase encoding (as S P E E D does) but also takes advantage of the fact that the signal in each ghosted image is likely weighted by a different coil sensitivity. Moreover, while the skip size N has to be a prime number in S P E E D , it can be any arbitrary integer i n S P E E D - A C E , which offers more flexibility in reducing scan time. S P E E D - A C E can be illustrated more specifically with the following example: 147 Chapter 6 Highly accelerated MRI with SPEED-ACE S (k) 4 E»(r) E (r) E,(r) 2 E«(r) B (r) 0 I„(r) Figure 6.2: K-space sampling scheme and data flow of 2 D S P E E D - A C E . B y using four receiver coils in parallel, k-space is selectively sampled with a skip size o f N into four sparse data sets: S^k), S (k), S (k) and S (k). They are then reconstructed separately by direct inverse 2 D F T into four ghosted images: h(r), I (r), I (r), and Li(r), with ghost overlapping up to N layers. After high-pass filtering, four ghosted edge maps — Ei(r), E (r), E3(r), and E4(r) — are obtained and subsequently combined into a deghosted edge map En(r). Finally, Eo(r) is inverse-filtered into a deghosted image In(r). 2 3 4 2 2 148 3 Chapter 6 Highly accelerated MRI with SPEED-ACE If P E is performed in the y-direction, sampling k space at every N P E step with four receiver coils yields four sensitivity-weighted and ghosted images, each with N-fold aliasing ghosts along y-direction. A t N = 5, the application of a y-directional differential filter to the ghosted images would reduce the number of overlapping layers from N = 5 to 2. Because the spatial differential of coil sensitivity profile is usually negligible compared with the image edges (as illustrated intuitively in the example shown in Figure 6.3), the high-pass filtered ghosted maps (i.e. ghosted edge maps) can be described at each pixel as E =Sl P G +SfP G 1 nl 14 n] d + S PfG (6.2) E =SP G + SfpfG (6.3) E = S4 P G nl 3 d n] 1 4d n 2 Hl ] 34 iid n2 E ,d = S?Pj'G 2 where E (6.1) n2 d ! d 2 n2 n2 + S% P G 2 nl (6.4) 2 d n2 stands for the ghosted edge map obtained from the data set sampled by the / coil with a relative shift of d in P E ; G„i and G n2 are the two dominant ghost layers; P is a ghost phasor known to have the form of (e depending th d 5 ) " , where n is the order of ghost on its relative location and must be chosen from 0, 1, 2, 3, and 4; and S" denotes the i t h coil sensitivity profile, shifted by - ^ - x FOV because of ghosting. It is noteworthy that S" is a coil sensitivity factor and P is a S P E E D factor. d 149 Chapter 6 Highly accelerated MRI with SPEED-ACE (a) (b) (c) (d) Figure 6.3: (a) is an M R image from a typical spin-echo scan on a clinical 1.5 T scanner equipped with multiple receiver coils and is denoted as C(x,y). (b) is the corresponding coil sensitivity map for one receiver coil and is denoted as S(x,y), where signal region is surrounded by noise background. S(x,y) was slightly smoothed to improve signal-to-noise ratio [49]. If C(x,y) is sensitivity weighted by S(x,y), the sensitivityweighted image can be written as S(x, y)C(x, y). If the vertical direction is defined as the y-direction, the y-directional spatial differential of S(x,y)C(x,y) consists of two terms: S(x,y+l)C(x,y+l)S(x,y)C(x,y) = C(x,y+l)[S(x,y+l)-S(x,y)] + S(x,y)[C(x,y+l)-C(x,y)]. W i t h the same scale, (c) shows the first term, and (d) shows the second term. Compared with (d), C(x,y+l)[S(x,y+l)-S(x,y)] shown in (c) can be assumed to be negligible. This is because coil sensitivity profile S(x,y) is varying slowly. For easy visualization, all images are shown here in magnitude although they are actually complex. 150 Chapter 6 Highly accelerated MRI with SPEED-ACE Similar to the ghost order index i n Simplified S P E E D and S P E E D , (nl, nl) at each pixel are not fully determined, but are known to be two different integers chosen from 0 to 4. B y considering (nl, n2) as parameters, Eqs. (6.1-6.4) are written into a matrix form as E = (6.5) MG where E =( E\ , E24, E34, E^d) "is a known vector, G = (G i,G 2) 1 is an unknown vector T id n n to be found, and M is a 4x2 matrix containing parameters (nl,n2) given by nn2 r» nl p"l ^1 nnl d r S P nl M = n«2 (6.6) -pnl cnl r>n2 ^4 °4 d r d r Equation (6.5) describes an over-determined linear system, i n which there are four equations and two unknowns. Therefore, a L S E solution vector GLSE can be solved for each possible pair o f (nl, nl) and as such is given by = LE LSE T (6.7) where L is referred to as "pseudo-inverse" o f matrix M and is explicitly calculated to be L= M M F \ M ^ 1 u U U -w -V* l 2 2 - V u. (s? P l l nl d &P l l nl d (s?P )[ nl d [s^pfj (sfpfj (sfpfj 151 {sfp J nl d (sfpfj (6.8) Chapter 6 Highly accelerated MRI with SPEED-ACE where Ui, U2, and Vare defined as u =s? *s + sfs? +sfs? + sfs? l nl { (6.9) l 77 _ U2 — ^ ri«2 ' V = (sfpf + , r*n2 r r « 2 ^2 2 } sfpf + 'K , c^nl r>n2 . r<n2 r*nl ^3 ^4 + (s^Pf } sfpf (6.10) + (sfpf } Sfpf + (sfpf } sfpf (6.11) and the superscripts "t", " - 1 " , and " * " denote respectively the transposed complex conjugate, the matrix inverse and the complex conjugate. The L S E solution vector GLSE can be explicitly expressed as, ( ni )LSE (G 2 )LSE. G n (sfpf} (sfpf} (sfpf} (sfpf} [sfpf} (sfpf} (sfpf} (sfpf} 1 u u -w l -V* 2 (6.12) 2,d E 3,d E 4,d E The L S E or residual energy for this solution can be expressed as LSE = energy[E - The noise variances of (G i)isE n and MG LSE (G 2)LSE n ] (6.13) depend on the "pseudo-inverse" matrix L and can be calculated as 152 Chapter 6 Highly accelerated MRI with cr 2 n l °«2 =( L , , + L 2 = (-^21 2 1 2 + ^22 SPEED-ACE +L 2 1 3 + ^23 +L + 2 1 4 )<7 ^24 •)(To (6.14) 2 0 (6-15) where Ly are elements of the matrix L defined in E q . (6.10), and <r 0 is noise variance of ghosted edge maps. For simplicity, we have assumed that edge maps from all four coils have equal noise variance <r 0 . Ghost order index (nl, n2) represent different orders of ghosts, so nl * n2 and as such the total number of possible pairs of (nl, n2) is N ( N - l ) / 2 =10. A l l the different pairs of (nl, nl) are tried for each pixel in the image, and the best pair is selected by minimizing the L S E in E q . (6.15) for that pixel alone. With the resolved (nl, n2), the correct solutions of (G i, n Gni) can be calculated by using E q . (6.12) on a pixel-by-pixel basis. Then, the edge ghosts (G i, G 2.) can be sorted out according to their associated (nl, n2) values and as n n such can yield five separate edge ghosts. The five separate edge ghosts can therefore be easily registered and averaged into a single deghosted edge map En(r) with reduced noise and artifacts in the same way as it is done in Simplified S P E E D and S P E E D . In addition, the square-root of the minimized L S E is output as a residual for quality control. Finally, S P E E D - A C E inverse-filters the deghosted edge map En(r) to yield a final deghosted image In(r) in the same manner as S P E E D [3]. T o avoid the artifacts in the inverse-filtering, the central k-space data (e.g. 32 out of 256 lines) sampled for measuring 153 Chapter 6 Highly accelerated MRI with SPEED-ACE coil sensitivity profiles are used to replace the inversed data at those k-space points near k = 0, making more efficient use of these data. In this way, only 1/5 k-space data with y full 32 central k-space lines are used to obtain a deghosted image with a high undersampling factor of 5, and scan time is consequently reduced to 30% i f full central kspace sampling is taken into account. The algorithm of S P E E D - A C E used i n this example can be summarized in its flowchart shown i n Figure 6.4. 154 Chapter 6 Highly accelerated MRI with SPEED-ACE Input 4 datasets with N =5 from 4 coils I Reconstruct into 4 sensitivity-weighted and ghosted edge maps E, , E , , E , and E b y high-pass filtering & F T d d 3 d 4 d I I I I I Yield the L S E solution of 2 dominating ghost layers G„,& G„, through 10 trials of n l & n2 at each pixel Sort out 5 separate edge ghosts G„ according to associated n 1 & n2 Register the 5 edge ghosts G„ and average them into a final deghosted edge map E„ Yield a final deghosted image 1„ by inverse filtering Output the final image 1„ Figure 6.4: Flowchart of S P E E D - A C E with single skipped sampling used to achieve an undersampling factor of 5 with four receiver coils. 155 Chapter 6 Highly accelerated MRI with SPEED-ACE 6.2.4 S P E E D - A C E w i t h M u l t i p l e S k i p p e d S a m p l i n g s In the above example, S P E E D - A C E uses one skipped sampling with multiple coils in parallel. This is, however, only one of many possible choices of sampling schemes. The method can be generalized for multiple skipped samplings with multiple receiver coils in parallel to achieve a better result. In the example, four equations — Eqs.(6.1-6.4) — are obtained, each describing a ghosted edge map. A s illustrated by the matrix M defined in E q . (6.6), edge deghosting depends on different coil sensitivity factors S", while S P E E D factors Pj 1 do not help in distinguishing ghost contributions. However, the role of S P E E D factors P£ w i l l become helpful i f four additional data sets are acquired with a relative shift in P E . The four equations — Eqs.(6.1-6.4) — w i l l be extended to eight equations, in which each weighting factor of the same edge ghost G„ w i l l be different from one another because either coil sensitivity factors S" or S P E E D factors P j , or both 1 are different. Thus, both S" and P j are used to distinguish ghost contributions and 1 thereby are both highly involved in edge deghosting. factor P j 1 In particular, the use of S P E E D in edge deghosting ensures good condition of the pseudo-inverse matrix L defined in E q . (6.8). B y repeating Eqs. (6.5-6.15), a deghosted edge map can be solved and be subsequently inverse-filtered into a final deghosted image. In this case, the undersampling factor is N / 2 , which is 50% greater than N / 3 achieved by S P E E D using a double-layer model and the same skip size of N . 156 Chapter 6 Highly accelerated MRI with SPEED-ACE 6.2.5 G e n e r a l i z e d S P E E D - A C E w i t h M u l t i - l a y e r M o d e l s Although the simple double-layer model works well in many cases, S P E E D - A C E can be generalized to include more layers in the model as chances of ghost overlapping increase. In particular, analogous to a Taylor expansion, ghosted edge maps can always be expanded through a G L U E series by progressively adding more layers of ghosts to increase accuracy. If all layers of ghost overlapping are included in the model, S P E E D A C E using single skipped sampling w i l l work in the same way as S E N S E [49]; however, rather than going that far, S P E E D - A C E has to find a balance point, because it faces the same challenge as S E N S E [87,88]: noise is amplified as the number of layers in the model increases. It is well known that noise amplification in S E N S E is quantified by using the "g-factor". Similarly, the behavior of amplified noise in S P E E D - A C E can be quantified by using noise variance of each deghosting solution on a pixel-by-pixel basis. In a Taylor expansion, a satisfactory approximation is often achieved by using the first several terms. Likewise, given the sparseness of edge maps, a satisfactory deghosted image can be achieved with a single-layer model or a double-layer model without significant noise amplification. The L S E deghosting algorithm ensures that the solutions based on a single-layer model or a double-layer model are the most dominant ones in the multiple layers of ghost overlapping. This implies that a good approximation can be rapidly reached with the first few layers. Since image quality is monitored by a residual map and noise variances of L S E solutions, generalized S P E E D - A C E can be used to 157 Chapter 6 Highly accelerated MRI with SPEED-ACE extend the deghosting solutions from a zero-layer model up to an N-layer model at each pixel until the residual is reduced to the level of the amplified noise at the same pixel. If the summation of residual and amplified noise is not tolerable, there is always an option of adding more data to improve the deghosting solution. In real-time M R I , this option provided by S P E E D - A C E enables its users to obtain a satisfactory result with good image quality in considerably less scan time. So that reduction of scan time and image quality can mostly be tailored to the object being imaged. 6.2.6 E x p e r i m e n t s S P E E D - A C E are tested with existing patient data obtained from a clinical spin-echo multi-slice coronal pelvis scan on a 1.5 T scanner ( E D G E , Picker International, Cleveland, O H ) (matrix 256x256, F O V 30 cm, T R 800 ms, T E 24 ms, slice thickness 4 mm, no average, 20 slices) equipped with four receiver coils. In this work, we assume that a desired image from a body coil is formed as root-mean-square ( R M S ) of the four individual magnitude images, each reconstructed from full k-space data of the corresponding receiver coil. To measure the sensitivity profile of each receiver coil, the central 32 k-space lines of the Fourier transform of the desired body image are sampled, and are then reconstructed by F T into one low-resolution body image. Similarly, the central 32 k-space lines of the Fourier transform of each receiver coil image are sampled, and are Fourier transformed into a corresponding low-resolution coil image. Finally, the 158 Chapter 6 Highly accelerated MRI with SPEED-ACE division of each low-resolution coil image by the low-resolution body image yields the sensitivity profile for the corresponding receiver coil. In this work, the desired body image is referred to as the "gold standard" for comparison with the images reconstructed by S P E E D - A C E . The S P E E D - A C E algorithm is implemented in C programming language on a personal computer with 2 G B R A M and a clock rate of 3.6 G H z . The reconstruction of a 256x256 image can be typically accomplished within three seconds. 6.3 Results Figures 6.5-6.9 show the results from the second slice of the pelvis scan, which is 4 m m away from the first slice. A s shown in Figure 8, four data sets are obtained by sampling k-space at every 5 th P E steps with four receiver coils in parallel. Since coil sensitivity profiles vary along the vertical direction but are similar along the horizontal direction, the vertical direction must be used as the P E direction to use spatially varying sensitivity profiles, which is similarly required also by other parallel imaging methods such as SENSE and SMASH. Figure 6.5 shows individual low-resolution coil images, sensitivity-weighted and ghosted coil images, and their edge maps. It can be seen that coil sensitivities i n this particular data set are not very different, which creates difficulties for other parallel imaging methods and also challenges the proposed method as a good test case. 159 Chapter 6 Highly accelerated MRI with SPEED-ACE •... Figure 6.5: Results from the second slice of the pelvis scan. In the top row are individual low-resolution coil images obtained from 32 central k-space lines. In the middle row are individual sensitivity-weighted and ghosted images, each reconstructed from the data set sampled with a skip size of 5. Aliasing ghosts of various orders are seen along vertical P E direction. In the bottom row are corresponding ghosted edge maps of middle row and shows much less ghost overlapping. 160 Chapter 6 Highly accelerated MRI with SPEED-ACE Figure 6.6 shows the results reconstructed by S P E E D - A C E with double skipped samplings, each accomplished by using a skip size of N = 5 steps and four receiver coils in parallel. Figure 6.6 (a) is the gold standard, the desired body image reconstructed from full k-space data. Figure 6.6 (b) is what Figure 6.6 (a) becomes after being superimposed by noise points randomly taken from its background, and is measured with E q . (4.23) to yield a T R E of 0.12, which is normalized square root of error power of noise. The difference between Figures 6.6 (a) and 6.6 (b) is hardly noticeable. Figure 6.6 (c) is one of the eight individual sensitivity-weighted and ghosted images. Figure 6.6 (d) is the R M S of the four individual low-resolution coil images shown in the top row of Figure 6.5, and shows significantly more blurring than the gold standard shown in Figure 6.6 (a). Figure 6.6 (e) is the R M S of coil images obtained by a direct inverse F T of all sampled data, including both fully sampled k-space center and undersampled outer regions of kspace where all the unsampled k-space points are filled with zero. Compared with the gold standard shown in Figure 6.6 (a), Figure 6.6 (e) is not only blurred, but also has strong visible aliasing artifacts. Figure 6.6 (f) is the deghosted image yielded by S P E E D A C E based on a bilayer model, with an undersampling factor of 5/2 = 2.5. Figure 6.6 (f) has a T R E of 0.082. Figure 6.6 (f) appears to be comparable to the gold standard in Figure 6.6 (a), and has a T R E value less than that of Figure 6.6 (b), which represents just the noise level. 161 Chapter 6 Highly accelerated MRI with SPEED-ACE • •91^ (a) (b) ii (c) it (d) - ^Mit.. L . / ' J (e) (f> Figure 6.6: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with double skipped samplings, each accomplished by using a skip size of N = 5 steps and four receiver coils in parallel, (a) is the desired body image considered as the gold standard, (b) is what (a) becomes after being superimposed by noise points randomly taken from the background of the real and imaginary parts of the gold standard shown in (a), resulting in a T R E of 0.12. (c) is one of the eight ghosted images, (d) is the R M S of the four individual low-resolution coil images shown in the top row of Figure 6.5. (e) is the R M S of the coil images obtained by a direct inverse F T of all sampled data, (f) is the reconstructed images by S P E E D - A C E based on a bilayer model with an undersampling factor of 2.5. (f) has a T R E 0.082, and appears to be comparable to the gold standard (a). 162 Chapter 6 Highly accelerated MRI with SPEED-ACE Figure 6.7 shows the results reconstructed by S P E E D - A C E with single skipped sampling. S P E E D - A C E samples k-space with the skip size N = 3 in P E by using four receiver coils in parallel, and reconstructs deghosted images with multi-layer models. Figure 6.7 (a) is one of the four ghosted images. Figure 6.7 (b) is the corresponding ghosted edge map. Figure 6.7 (c) is R M S of the coil images obtained by a direct inverse F T of all sampled data. Figures 6.7 (d) and (e) are the deghosted images yielded by S P E E D - A C E based respectively on a single-layer model and a double-layer one. Figure 6.7 (d) has a T R E 0.09, while Figure 6.7 (e) has a T R E of 0.11. Figures 6.7 (d) and (e). Figure 6.7 (f) is the simple average of W i t h an undersampling factor of N = 3, all the reconstructed S P E E D - A C E images have the errors comparable to the noise level of the gold standard shown in Figure 6.6 (a) quantified by T R E . 163 Chapter 6 Highly accelerated MRI with (d) SPEED-ACE (e) (f) Figure 6.7: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with single skipped sampling. S P E E D - A C E samples k-space with a skip size of N = 3 in P E by using four receiver coils in parallel and reconstructs deghosted images with multiple-layer models. Figure 6.6 (a) is the gold standard reconstructed from full kspace data, (a) is one of the four individual sensitivity-weighted and ghosted images, (b) is corresponding ghosted edge map of (a), (c) is the R M S of coil images obtained by a direct inverse F T of all sampled data, (d) and (e) are the reconstructed images by S P E E D A C E based respectively on a single-layer model and a double-layer one. (d) has a T R E 0.09 while (e) has a T R E of 0.11. (f) is the simple average of (d) and (e). W i t h an undersampling factor of 3, all the images shown in (d-f) are quite good compared with the gold standard shown in Figure 6.6 (a). 164 Chapter 6 Highly accelerated MRI with SPEED-ACE Figure 6.8 shows the results reconstructed by S P E E D - A C E with single skipped sampling. S P E E D - A C E samples k-space at N = 5th P E steps by using four receiver coils in parallel, and reconstructs deghosted images with multi-layer models. Figure 6.8 (a) is one of the four ghosted images with five-fold ghosts shown i n the middle row of Figure 6.5. Figure 6.8 (b) is R M S of coil images obtained by a direct inverse F T of all sampled data, which is significantly blurred. Figures 6.8 (c-f) are the deghosted images yielded by S P E E D A C E based respectively on zero-layer, single-layer, double-layer, and triple-layer models. A s shown in Figures 6.8 (c-f), noise gets increasingly amplified as the number of overlapping layers increases, similar to what occurs in S E N S E [87,88], which solves for all N layers at a high price of noise amplification. Figure 6.8 (g) is the final image reconstructed by generalized S P E E D - A C E with single skipped sampling. Figure 6.8 (g) is measured with E q . (4.23) to have a T R E of 0.12, which is comparable to the noise level quantified by T R E of Figure 6.6 (b). Figures 6.8 (h-k) are automatically generated model maps indicating the number of overlapping layers considered at each pixel for generating Figure 6.8 (g). Figure 6.8 (h) is the zero-layer model map, in which a white pixel has a value of unity representing a zero-layer model applied at that pixel. Likewise, Figures 6.8 (i-k) are respectively single-layer, double-layer, and triple-layer model maps. Figure 6.8 (1) is the corresponding residual map of Figure 6.8 (g), and is used for quality control. A s compared with the four ghosted edge maps shown in the bottom row of Figure 6.5, Figure 6.8 (1) shows much less pixel intensities, indicating the adequacy of the modeling. A key feature of Figure 6.8 (g) is that the solution based on a single-layer model seems 165 Chapter 6 Highly accelerated MRI with SPEED-ACE dominant in the final image shown in Figure 6.8 (g). This is because, as shown in ghosted edge maps in Figure 6.5, the majority of ghost edges are not overlapped with one another, and as such can be adequately modeled with a single-layer model. W h i l e scan time is greatly reduced with an undersampling factor of 5 (or an effective speedup factor of 10/3 when full central k-space sampling is taken into account), the image shown in Figure 6.8 (f) has quite good image quality as compared with the gold standard shown in Figure 6.6 (a). This implies that the scan time can be reduced from 3.4 minutes (or 256 x800ms) to 1 minute for this pelvis scan. The results demonstrate that S P E E D - A C E can achieve an undersampling factor that is even greater than the number of coils, which is considered impossible for other parallel imaging methods such as S E N S E [49] and S M A S H [48]. This is because S P E E D - A C E is proposed by extending the principle of S P E E D to multiple coils, and S P E E D itself is able to accelerate M R I scan considerably even with only a single coil [3]. 166 Chapter 6 Highly accelerated MRI with SPEED-ACE Chapter 6 Highly accelerated MRI with SPEED-ACE Figure 6.8: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with single skipped sampling. S P E E D - A C E samples k-space with a skip size of N = 5 in P E by using four receiver coils and reconstructs deghosted images with multi-layer models. Figure 6.6 (a) is the gold standard reconstructed from full k-space data, (a) is one of the four individual sensitivity-weigh ted and ghosted images, (b) is the R M S of the coil images obtained by direct inverse F T of all sampled data, (c-f) are the reconstructed images by S P E E D - A C E based respectively on zero-layer, single-layer, double-layer, and triple-layer models, (g) is the final image reconstructed by generalized S P E E D - A C E and has a T R E of 0.12. With an undersampling factor of 5, the final image shown in (g) has the errors comparable to the noise level of the gold standard shown in Figure 6.6 (a), (h-k) are automatically generated model maps for (g), which indicate the number of overlapping layers in the model used at each pixel. M o r e specifically, (h) is the zerolayer model map, in which a white pixel has a value of unity representing a zero-layer model applied at that pixel. Likewise, (i-k) are respectively the single-layer, doublelayer, and triple-layer model maps, (i) is the residual map of (g). Figure 6.9 shows the results from the entire 20 slices of the coronal pelvis scan. W i t h a sampling scheme identical to that used in Figure 6.8, the entire 20 slices are reconstructed by generalized S P E E D - A C E with single skipped sampling, and show consistent results ranging from the most superficial slice to the deepest one. Compared with the gold standard shown in the top four rows of Figure 6.9, the images reconstructed by S P E E D A C E all have quite good image qualities with a high undersampling factor of 5. T R E values of these images are measured and illustrated in Table 6.1. 168 The Chapter 6 Highly accelerated MRI with •••• • SPEED-ACE .. . s : fill w 1 %m *\.iT*" I •'• ii _ I — — • _ — <* \ » §i \ £jk *» Figure 6.9: Results from the 20-slice pelvis scan: In the top four rows are the gold standard from the 1 slice to the 2 0 slice, each reconstructed from full k-space data. In the bottom four rows are reconstructed images by generalized S P E E D - A C E with an undersampling factor of 5. st th 169 Chapter 6 Highly accelerated MRI with Slice Distance from No.l Slice SPEED-ACE Total Relative Error Slice TRE Distance from No.l Slice Total Relative Error TRE 0 mm 0.144 No. 11 40 mm 0.139 4 mm 0.117 No. 12 44 mm 0.134 8 mm 0.118 No. 13 48 mm 0.140 12 mm 0.129 No. 14 52 mm 0.136 No. 5 16 mm 0.141 No. 15 56 mm 0.139 No.6 20 mm 0.145 No. 16 60 mm 0.130 No.7 24 mm 0.148 No. 17 64 mm 0.123 No.8 28 mm 0.151 No. 18 68 mm 0.115 No.9 32 mm 0.143 No. 19 72 mm 0.120 36 mm 0.143 No.20 76 mm 0.117 No.l No.2 No.3 No.4 No. 10 Table 6.1: T R E values for the images reconstructed by generalized S P E E D - A C E shown in the bottom four rows of Figure 6.9. 170 Chapter 6 Highly accelerated MRI with SPEED-ACE 6.4 Discussion In summary, S P E E D - A C E undersamples k-space with multiple coils in parallel and reconstructs satisfactory images with a high undersampling factor, which can be greater than the number of coils in certain applications, which is so far impossible for other parallel imaging methods. This part of the chapter addresses further several issues regarding S P E E D - A C E . A s discussed previously, S P E E D - A C E reduces scan time by taking advantage of different coil sensitivities and sparse edge maps. It follows that image quality is affected by two factors: noise amplification and artifacts, depending on the number of layers used in the model of ghosted edge maps. When the layers in the model are not sufficient to handle edge-rich regions, artifacts become significant in those regions and blurring effects w i l l appear in the image reconstructed by S P E E D - A C E . Conversely, as the number of layers in the model increases, the blurring effect rapidly decreases, but noise is significantly amplified, similar to what occurs in S E N S E [87,88], which always effectively includes all N layers of ghosts. Noise is amplified because S P E E D - A C E relies heavily on the coil sensitivity factor S" in distinguishing ghost contributions. To constrain noise amplification, S P E E D - A C E can use multiple skipped samplings at the expense of greater scan time. Alternatively, S P E E D - A C E can refine the deghosting solution on each pixel by progressively including additional layers till a balance can be reached between the degrading effects caused by modeling artifacts and those caused by amplified noise. 171 Chapter 6 Highly accelerated MRI with SPEED-ACE In S P E E D - A C E , image quality is monitored by a residual map similar to the residual in S P E E D , and the noise variances of L S E solutions similar to the "g-factor" in S E N S E . In S P E E D - A C E , residual is used to measure modeling artifacts, and noise variance is used to quantify noise amplification. A s the number of layers in the model increases, the residual is reduced (similar to what occurs in S P E E D [3]), but noise is accordingly amplified (similar to what occurs in S E N S E [87,88]). Therefore, S P E E D - A C E should extend its deghosting solution to a cutoff layer at which the total errors caused by both residue and amplified noise is minimized. This process ensures a reasonably good image with a high undersampling factor. In real-time M R I , additional data can always be acquired to improve image quality till the reconstructed image becomes satisfactory. In this work, S P E E D - A C E is demonstrated with in vivo data obtained with four receiver coils; however, none of them has well localized sensitivity patterns as shown in Figure 6.5. If coils with well localized sensitivity patterns are available, P I L S [51] can be used to accelerate M R I [51], and S P E E D - A C E can be improved by taking advantage of the increased independence of individual coils or the increased difference of coil sensitivities. Alternatively, because each coil image is well localized and thus has a relatively sparse signal distribution, each coil image can be reconstructed separately by using S P E E D with a high undersampling factor. A l l the reconstructed coil images can be then combined into a final image through R M S . In this case, scan time is still significantly reduced, while noise is not greatly amplified because coil sensitivity factors — which would otherwise amplify noise — are excluded from the pseudo-inverse matrix. 172 Chapter 6 Highly accelerated MRI with SPEED-ACE In S P E E D - A C E , fully sampled central k-space data are used to serve two purposes: calculating sensitivity profiles and aiding the inverse-filtering. In practice, we have found that the size of full central k-space sampling has an impact on the final deghosted image. In this work, 32 out of 256 lines are chosen as a reasonable compromise between image quality and scan-time reduction. The size of fully sampled central k-space should be optimized for improved performance, which may be addressed in the future. It has been reported that S P E E D has its simplified version. Similarly, S P E E D - A C E can also be simplified to enhance its efficiency by taking advantage of the sparseness of M R A data. Since M R A data usually have very sparse signal distribution, differential filtering is not used any more and the assumption of slowly varying sensitivity profiles is naturally omitted. In other words, Simplified S P E E D - A C E can achieve improved performance with more straightforward reconstruction. Moreover, to reduce scan time further, S P E E D - A C E and its simplified version can also be combined with other existing fast imaging techniques such as Half-Fourier imaging [39]. In conclusion, S P E E D - A C E is a new parallel imaging method using edge deghosting. It can yield rather satisfactory results with a high undersampling factor, and has many potential 2 D and 3 D applications. 173 Chapter 7 SNR Improvement for SPEED and Its Extensions Chapter 7 SNR Improvement for SPEED and Its Extensions Several fast imaging methods have been proposed in the previous four chapters. In Chapter 4, Simplified S P E E D is proposed to accelerate M R A by taking advantage of the sparseness of M R A data. After this, the original S P E E D is reviewed. Its extension to 3D and its combination with Half-Fourier imaging have been explored in Chapter 5. Both Simplified S P E E D and S P E E D are efficient and robust; however, they are not making full use of multiple coils available on modern M R scanners. T o improve upon this, a new parallel imaging method named S P E E D - A C E is proposed in Chapter 6, which combines the principle of S P E E D with that of sensitivity encoding using multiple coils. Similar to S P E E D , S P E E D - A C E can also be potentially simplified to further accelerate M R A by using its sparse nature. Chapters 4, 5 and 6 focus on new fast M R I methods. As a continuation of this work, noise behaviors in the proposed fast M R I methods w i l l be studied and a novel approach for S N R improvement w i l l be proposed in this chapter. 174 Chapter 7 SNR Improvement for SPEED and Its Extensions 7.1 Introduction Signal-to-noise considerations are very important in M R I . In conventional M R I , it is well known that S N R w i l l decrease as less k-space data are sampled [39,48,49]. For example, an image can be reconstructed from slightly more than half of k-space data by using Half-Fourier imaging; however, the S N R of the image w i l l be reduced by a factor o f - 4 2 , compared with the S N R of the image reconstructed from full k-space data [39]. Consider the k-space matrix in an ordinary 2 D F T - M R I experiment [89]. S N R w i l l be reduced by a factor of nearly VN" i f k-space is partially sampled with a skip size of N . In this work, it has been found that most of the S N R loss caused by k-space undersampling can be recovered in S P E E D . Based on this finding, we propose an approach to improve image S N R for S P E E D and its extensions. This finding and the approach are both demonstrated in this chapter with computer simulations and in vivo data. The following discussions focus on S P E E D because its extensions such as Simplified S P E E D and S P E E D - A C E use similar deghosting algorithms. 175 Chapter 7 SNR Improvement for SPEED and Its Extensions 7.2 Methods 7.2.1 S h e p p - L o g a n p h a n t o m To study noise behaviors in S P E E D , a computer simulation is conducted, in which a digital Shepp-Logan phantom [90] with simulated Gaussian noise [69] is used. The Shepp-Logan phantom used i n this work is originally developed for C T and has been used as a standard test phantom to mimic the anatomy of the brain [90]. In this work, the Shepp-Logan phantom is used in an M R I computer simulation, which can be divided into three steps: (1) First, the Shepp-Logan phantom is used as the real part of an ideal noise-free complex M R image while its imaginary part is set to zero. (2) Second, a set of complex-valued "white"-Gaussian noise is simulated, which is uniformly distributed (i.e. "white") across the entire k-space and shows a Gaussian distribution in histogram [69]. (3) Finally, a simulated noisy M R image is generated by adding the simulated "white"Gaussian noise to the ideal noise-free M R image. In order to quantify image quality, this work uses signal-to-noise ratio ( S N R ) and defines it as the ratio of the root-mean-square ( R M S ) signal and noise across the entire imagespace field-of-view ( F O V ) . The S N R can be mathematically expressed as 176 Chapter 7 SNR Improvement for SPEED and Its Extensions SNR = RMS( Signal) RMS( Noise) (7.1) The computer simulation and S N R are further illustrated below with an example. A s shown in Figure 7.1, a Shepp-Logan phantom is superimposed by Gaussian noise to yield a simulated noisy M R image, where S N R is measured to be 11.5 by using E q . (7.1). The simulated M R images can be used to demonstrate skipped sampling in the next section. Figure 7.1: Simulated M R images with a digital Shepp-Logan phantom: (a) is the ideal noise-free complex M R image; (b) is the noisy M R image imposed with simulated "white"- Gaussian noise. 7.2.2 Sampling with Skipped Phase Encoding A s discussed previously, skipped sampling means sampling k-space with a skip size of N in P E . If k-space is sampled at every N t h P E steps, the sampled data can be reconstructed by direct inverse F T into a full F O V ghosted image, which contains N aliasing ghosts along the P E direction. If the reconstructed image from full k-space data is considered as 177 Chapter 7 SNR Improvement for SPEED and Its Extensions a desired full-FOV image, each aliasing ghost can be treated as the desired image after being shifted by a certain distance and being rotated in phase by a certain angle. Such a repetitive ghost pattern implies that only 1/N of the F O V is independent and the rest of F O V is just redundant. In other words, sampling k-space with a skip size of N imposes a periodic signal modulation and therefore reduces the F O V by a factor of N , resulting in ghost aliasing. This repetitive ghost pattern can be demonstrated in Figure 7.2, which shows that every 1/4 F O V is identical to one another due to a skipped P E sampling of 4. (c) (d) Figure 7.2: Ghosted images obtained with a skip size of N = 4 and their corresponding edge maps, (a) is noise-free ghosted image obtained by sampling k-space with a skip size of N = 4. (b) is the edge map of (a), (c) is simulated noisy ghosted image obtained with N = 4. (d) is the edge map of (c). It can be shown that every 1/4 of F O V is just the same as one another. 178 Chapter 7 SNR Improvement for SPEED and Its Extensions 7.2.3 S N R I m p r o v e d b y D a t a L e a k a g e o f D F T i n S P E E D Ideally, every 1/N F O V of the ghosted image should be identical in both signal and noise. In practice, however, they are found to be slightly different when full F O V size is not an integer multiple of N , because data leakage of discrete F T occurs i f the total number of steps in a normal P E is not divisible by N . Despite being quite small, data leakage is illustrated in Figure 7.3, which shows the noise-free ghosted images and edge maps obtained with a skip size of N = 5. A s shown in Figure 7.3 (e), magnitudes in edge regions between ghosts of different orders are considerably different, although ideally they should be identical; likewise, magnitudes in non-edge regions are non-zero, although ideally they should be zero. Moreover, because data leakage is location-dependent, it changes the periodic ghost pattern and thereby makes each 1/N F O V slightly different from one another. 179 Chapter 7 SNR Improvement for SPEED and Its (a) (b) (c) (d) Extensions (e) Figure 7.3: Ghosted images obtained with a skip size of N = 5 and corresponding edge maps, (a) is the ideal noise-free complex M R image, which is identical to Figure 4 (a), (b) is one of three ghosted images, (c) and (d) are respectively the edge maps of (a) and (b). (e) is profile of column N o . 128 of (d). It can be easily seen from (e) that every a 1/5 F O V is slightly different from one another. 180 Chapter 7 SNR Improvement for SPEED and Its Extensions Based on this fact, S P E E D processes f u l l - F O V data to yield N separate edge ghosts, which are slightly different and as such can be averaged into a final deghosted edge map with an improved S N R . The question arises why a small data leakage can play a big role in S P E E D ; the question can be answered as follows: Rather than unfolding all signal superposition like S E N S E does, S P E E D reduces ghost overlapping by turning ghosted images into sparse edge maps with a differential filter. G i v e n the sparseness of edge maps, each edge map can be modeled with only two dominant layers denoted as (G„i, G i), n associated with unique ghost order index (nl,n2). Once (nl,n2) are determined, (G„],Gn2) can be solved and be sorted out according to (nl,n2), yielding N separate edge ghosts. In S P E E D , solving (nl,n2) based on a double-layer model is quite reliable in edge regions, but not so in other regions where noise is dominating. Therefore, noiselevel non-periodic data leakage can easily affect the solution of (nl,n2) in noise regions while leaving edge regions intact, resulting in a non-periodic (nl,n2) distribution. This implies that all of the N separate edge ghosts are different, especially in noise regions, and therefore should be averaged to improve S N R . Then the averaged edge map is inverse-filtered to yield a final deghosted image, which has S N R comparable to the "gold standard" reconstructed from full k-space data. Such S N R improvement by data leakage of D F T can be demonstrated with computer simulations as follows: A s shown in Figure 7.4, the noise-free M R image is Fourier transformed into k-space, and is sparsely sampled into three interleaved k-space data sets, each with a skip size of 5 and a different relative shift in vertical direction. The three data sets are processed by 181 Chapter 7 SNR Improvement for SPEED and Its Extensions S P E E D to yield five separate edge ghosts. Figure 7.4 (d) is one of the five edge ghosts. Figure 7.4 (e) is the average of the five edge ghosts and is inverse-filtered into the final deghosted image shown in Figure 7.4 (f). Figure 7.4 (f) is measured to have a T R E of 0.007, which is considered as the error caused by artifacts in the noise-free situation. Figure 7.5 shows the deghosted edge maps and images yielded by S P E E D (with noise being taken into account) in the computer simulations. Figure 7.5 (a) is the simulated noisy complex M R image, which is the same as Figure 7.1 (b). Figure 7.5 (b) is the edge map of Figure 7.5 (a) with a S N R of 1.5. Figure 7.5 (c) is one of the five deghosted edge maps. It is noteworthy that Figure 7.5 (c) has a S N R of 1.3, which is greater than theoretically expected value 1.5/yfl ~ 0.7. The greater S N R can be justified by the fact that each deghosted edge map solved by S P E E D has 3/5 noise regions automatically filled with zero due to the double-layer modeling. If this is taken into account, the S N R of 1.3 implies a factor of 1.3xV2/5 =0.8, which is mostly consistent with the theoretically expected value. Figure 7.5 (d) is the average of the five deghosted edge maps with a S N R of 1.5, which is the same as S N R of Figure 7.5 (b). Figure 7.5 (f) is the final deghosted image obtained by inverse-filtering Figure 7.5 (d). The image quality of the final deghosted image is measured with S N R and appears to be comparable to the "gold standard" in Figure 7.5 (a). B y using E q . (4.23), the T R E of Figure 7.5 (f) is measured to be 0.08, which is considered as the total error caused by both artifacts and noise in the noisy situation. Assume artifacts are not significantly affected by the added 182 Chapter 7 SNR Improvement for SPEED and Its Extensions noise as suggested by the stable solutions shown in both Figure 7.4 and Figure 7.5. It can be seen that the error caused by artifacts (quantified by a T R E of 0.007) is much less than noise (quantified by a T R E of -y/(0.08) - ( 0 . 0 0 7 ) 2 2 = 0.079). In this case, we obtain a deghosted image at the cost of mostly negligible artifacts, but with a comparable S N R to that of the gold standard. (d) (e) (£) Figure 7.4: Noise-free deghosted images and edge maps yielded by S P E E D with a skip size of 5. (a) is the ideal noise-free complex M R image, which is the same as Figure 7.1 (a), (b) is the edge map of (a), (c) is one of the three ghosted edge maps, each obtained by sampling k-space with a skip size of 5. (d) is one of the five deghosted edge maps, (e) is the average of the five deghosted edge maps, (f) is the final deghosted image, which is obtained by inverse-filtering (e). (f) has a T R E of 0.007. 183 Chapter (d) 7 SNR Improvement for SPEED and Its (c) Extensions (f) Figure 7.5: Noisy deghosted images and edge maps yielded by S P E E D with a skip size of N = 5. (a) is the simulated noisy complex M R image, which is the same as Figure 7.1 (b) with an S N R of 11.5. (b) is the edge map of (a) with an S N R of 1.5. (c) is one of the five deghosted edge maps, and has an S N R of 1.3. (d) is the average of the five deghosted edge maps, and has an S N R of 1.5, which is the same as S N R of (b). (e) is obtained by inverse-filtering (c) and has an S N R of 9.0. (0 is the final deghosted image, which is obtained by inverse-filtering (d) and has an S N R of 11.5, which is the same as S N R of the original image in (a), (f) is measured to have a T R E of 0.08. 184 Chapter 7 SNR Improvement for SPEED and Its Extensions 1.2.4 P r o p o s e d M e t h o d to I m p r o v e S N R for S P E E D and its E x t e n s i o n s Despite being usually small and even smaller than noise, data leakage plays a significant role in determining (nl, nl) in noise regions, leading to an improved S N R . However, data leakage has negative effects that surely degrade the quality of the final image. Moreover, data leakage occurs only i f full F O V size is not divisible by the skip size of N . Based on the understanding of data leakage and its effects, this work proposes a new method to improve S N R in a similar way to data leakage, but without its "side effects". A s discussed previously, the determination of (nl, nl) is very crucial because (G \, n are solved based on (nl, nl). G i) n If we can randomize (nl, nl) i n noise regions while leaving edge regions intact, noise w i l l be reduced through averaging i n noise regions, and signal w i l l remain unaffected in edge regions. T o do so, a certain amount of random noise is simulated and added to the complex data ( E i , E2, E 3 ) to randomize (nl, nl) distribution like data leakage does. The added random noise must be below noise level so that solutions of (nl, nl) in edge regions are not affected. In practice, we firstly measure the noise level of image background, then add into the complex data ( E i , E2, E 3 ) simulated Gaussian noise comparable to 10% of the measured value, and finally use these data for determining (nl, nl). The (nl, nl) solved through the above process are used for edge deghosting with the original complex data ( E i , E2, E ) before the 10% simulated noise is 3 added, resulting in a final deghosted image with an improved S N R . ! It should be emphasized that the purpose of adding noise is only to randomize (nl, nl) in non-edge 185 Chapter 7 SNR Improvement for SPEED and Its Extensions regions, and that the final L S E solution is still obtained from the original data ( E i , E , E 3 ) 2 once (nl, nl) are determined. The feasibility of the proposed method has been demonstrated in S P E E D as shown in Figure 7.6. This S N R improvement method is especially important for S P E E D - A C E , where an arbitrary skip size can be used and consequently data leakage may not be available. 186 Chapter 7 SNR Improvement for SPEED and Its Extensions <d> (0 (e) Figure 7.6: Noisy deghosted images and edge maps yielded by S N R improved S P E E D with a skip size of N = 5. (a) is the simulated noisy complex M R image, which is the same as Figure 7.1 (b) with an S N R of 11.5. (b) is the edge map of (a) with an S N R of 1.5. (c) is one of the five deghosted edge maps, and has an S N R of 1.3. (d) is the average of the five deghosted edge maps, and has an S N R of 1.7, which is slightly greater than S N R of (b). (e) is obtained by inverse-filtering (c) and has an S N R of 9.0. (f) is the final deghosted image, which is obtained by inverse-filtering (d) and has an S N R of 12.0, which is slightly greater than the S N R of Figure 7.5 (f) obtained without randomizing (nl, nl). In addition, the T R E of (f) is measured to be 0.075, which is slightly less than that of Figure 7.5 (f). 187 Chapter 7 SNR Improvement for SPEED and Its Extensions 7.2.5 E x p e r i m e n t s The S N R improvement method was used to improve S P E E D - A C E and was tested with existing in vivo data obtained from a clinical spin-echo multi-slice coronal pelvis scan on a 1.5 T scanner ( E D G E , Picker International, Cleveland, O H ) (matrix 256x256, F O V 30 cm, T R 800 ms, T E 24 ms, slice thickness 4 mm, no average, 20 slices) equipped with four receiver coils. The algorithm of S P E E D - A C E with S N R improvement was implemented in C programming language on a personal computer with 2 G B R A M and a clock rate of 3.6 G H z . The reconstruction time is nearly identical to that of typical S P E E D - A C E because the computing time of randomizing the ghost order index (nl,n2) is almost negligible . 7.3 Results Figure 7.7 shows the results reconstructed by S P E E D - A C E with single skipped sampling, which is accomplished by using a skip size of N = 4 steps and four receiver coils in parallel. In this case, data leakage does not occur because 256 is divisible by 4. The proposed S N R improvement method should therefore be used to improve image quality. A s shown in Figure 7.7, the quality of the image reconstructed by S P E E D - A C E is significantly enhanced after the proposed S N R improvement method is applied. The improvement has been quantified consistently by the T R E ; Figure 7.7 (b) has a T R E of 0.18, and Figure 7.7 (c) has a T R E of 0.15. 188 Chapter 7 SNR Improvement for SPEED and Its Extensions (a) (b) (c) Figure 7.7: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with single skipped sampling, which is accomplished by using a skip size of N = 4 steps and four receiver coils in parallel, (a) is the gold standard reconstructed from full kspace data, (b) is the image reconstructed by S P E E D - A C E based on a double-layermodel without S N R improvement, with an undersampling factor of 4. (c) is the image reconstructed by S P E E D - A C E based on a double-layer model with S N R improvement, with an undersampling factor of 4. (b) has a T R E of 0.18, and (c) has a T R E of 0.15. Figure 7.8 shows the results reconstructed by S P E E D - A C E with double skipped samplings, each accomplished by using a skip size of N = 8 steps and four receiver coils in parallel. In this case, data leakage does not occur because 256 is divisible by 8. The proposed S N R improvement method should therefore be used to improve image quality. Figure 7.8 (a) is one of the eight sensitivity-weighted and ghosted images. It shows that all aliasing ghosts are greatly overlapped. Figure 7.8 (b) is the corresponding ghosted edge map of Figure 7.8 (a). In Figure 7.8 (b), ghost overlapping is much less than that in Figure 7.8 (a), and most pixels seem to be adequately represented with a single-layer or a double-layer model. Figure 7.8 (c) shows the image reconstructed by S P E E D - A C E based on a double-layer model with an undersampling factor of 8/2 = 4. Figure 7.8 (d) is 189 Chapter 7 SNR Improvement for SPEED and Its Extensions reconstructed in a similar way to Figure 7.8 (c), but with an improved S N R . Compared with the gold standard shown in Figure 7.7 (a), Figure 7.8 (c) has a T R E of 0.14, while Figure 7.8 (d) has a T R E of 0.11. It can be seen that image quality has been improved by using the proposed S N R improvement method, as indicated both by the T R E value and by visual inspection. (a) (b) (c) (d) Figure 7.8: Results from the second slice of the pelvis scan reconstructed by S P E E D A C E with double skipped samplings, each accomplished by using a skip size of N = 8 steps and four receiver coils in parallel. Figure 7.7 (a) is the gold standard with full kspace sampling, (a) is one of the eight individual sensitivity-weighted and ghosted images, (b) is the corresponding ghosted edge map of (a), (c) is the reconstructed images by S P E E D - A C E based on a double-layer model with an undersampling factor of 4. (d) is reconstructed similar to (d) but with an S N R improvement, (c) has a T R E 0.14, while (d) has a T R E of 0.11, which is comparable to the noise level quantified by T R E in Figure 6.6 (b). 190 Chapter 7 SNR Improvement for SPEED and Its Extensions 7.4 Summary A s continuation of the work on fast M R I , noise behaviors in S P E E D are studied and a novel approach is proposed for S N R improvement in this chapter. In the first part of this chapter, data leakage of discrete F T is discussed and its effects on S N R improvement are explored. A s discussed previously, the skip size N used in S P E E D is always a prime number in order to avoid ghost phase degeneracy. Because the most commonly used F O V sizes in M R I are usually not divisible by a prime number other than 2, data leakage usually occurs in original S P E E D . A s illustrated previously, data leakage in S P E E D can improve S N R only to a limited extent, because data leakage does not yield a totally random field and is not even available when full F O V size is divisible by the skip size of N. In this chapter, we propose a practical method, which does not use data leakage but has more potential in achieving optimal S N R improvement. This new method uses a random noise field, and thereby yields a nearly random (nl, n2) distribution in non-edge regions. Moreover, the new method can be used regardless of whether full F O V size is divisible by the skip size of N . W i t h in vivo data, the new method has been demonstrated to successfully improve the images reconstructed by S P E E D - A C E , particularly when data leakage is not available. 191 Chapter 8 Conclusions and Recommendations for Future Work Chapter 8 Conclusions and Recommendations for Future Work This chapter briefly reviews this thesis and proposes future work exploring further the potentials of phase information in M R I . 8.1 Conclusions of the Thesis This thesis presents two major applications of M R I phase information: contrast enhancement in Inversion Recovery (ER) imaging by using an improved statistical nonlinear phase correction method, and scan time reduction through strategic spatial phase encoding. The former is discussed in chapter 3, and the latter is elaborated in Chapters 4, 5, 6, and 7. A s discussed in Chapter 3, polynomial based phase correction is improved by using an extended statistical algorithm, which is more effective for handling non-linear phase terms. The improvement has been successfully demonstrated with 2 D in vivo IR imaging data [2]. The most significant contribution of this work is the introduction of a general npixel-shift R D F , in which signal can be enhanced effectively by a larger pixel shift n while the corresponding noise is not amplified. The unique phase signal behavior of npixel-shift R D F might have other potential applications, which are beyond the scope of phase correction described i n this work. 192 Chapter 8 Conclusions and Recommendations for Future Work From Chapters 4 to 6, a fast imaging method named S P E E D [3] and its extensions are investigated. Based on the principle of S P E E D , a new fast M R A method named Simplified S P E E D is developed, which has been demonstrated with in vivo data to successfully accelerate phase-contrast head scan by a factor close to 6 with a single coil [4]. Moreover, a new parallel imaging method named S P E E D - A C E is developed to highly accelerate M R I [5]. A s demonstrated with clinical data, the method can successfully accelerate pelvis scan by a factor close to 5 with only four receiver coils, which is considered impossible with other parallel imaging methods such as S E N S E [49] and S M A S H [48]. In addition, S P E E D is also extended from 2 D to 3 D for more flexible reduction of scan time, and is combined with Half-Fourier imaging to accelerate M R I further by a factor of nearly 2. The key idea of S P E E D and its extensions is to use partial information to reconstruct a complete image based on the sparse signal distribution of M R A or image edge. In S P E E D and its extensions, k-space is undersampled and S N R is an important consideration. In Chapter 7, I also study noise behaviors and demonstrate that S N R loss caused by undersampling can be mostly recovered in S P E E D . Based on the study, a novel S N R improvement method is proposed for S P E E D , Simplified S P E E D , and S P E E D - A C E . 193 Chapter 8 Conclusions and Recommendations for Future Work 8.2 Recommendations for Future Work Based on the above work, I propose the following research in the areas of flow imaging, water-fat imaging, fast imaging, and S N R improvement in M R I . 8.2.1 Other Applications of Phase Correction with N-Pixel-Shift RDF A s discussed in Chapter 3, we have found that a larger pixel shift can enhance phase signal in n-pixel-shift R D F without amplifying the corresponding noise, allowing for more accurate phase correction. Although it is applied and demonstrated with only IR imaging in this thesis, n-pixel-shift R D F can be used to improve other M R I techniques that require phase correction. Amongst them are flow imaging and water-fat imaging, which are briefly reviewed in Chapter 2. For example, in the single-point Q S water-fat imaging, water and fat images can be separated by using one single Q S acquisition in the absence of phase errors. However, the Q S method is limited because phase errors do occur in reality. It has been reported that the original A C method was used for linear phase correction in the QS method [32]. In the future, the extended A C method based on n-pixel-shift R D F may be used to improve the QS method. Another example is quantitative flow imaging, in which the velocity of blood flow is encoded into phase. A s discussed in Chapter 2, flow velocities can be mapped simply by calculating the phase differences between a velocity-sensitized image and a reference one in P C M R A , in which most of phase errors are removed [12]. However, the remaining phase errors after 194 Chapter 8 Conclusions and Recommendations for Future Work this normal phase correction by subtraction may often be problematic and as such should be removed. It has been recently reported that polynomial-based phase correction can be used to improve Phase-Contrast velocity quantification of pulmonary artery flow [93]. In the future, a novel method based on n-pixel-shift R D F may be proposed to further improve qualitative flow imaging. 8.2.2 H i g h l y A c c e l e r a t e d M R A w i t h S i m p l i f i e d S P E E D - A C E M y future work w i l l also cover fast M R A . A s discussed in Chapter 4, M R A is usually quite time consuming and inefficient, because it is often performed in 3 D and usually acquires less information than typical M R I using the same scan time. T o reduce scan time in M R A , a method named Simplified S P E E D is proposed in Chapter 4. Simplified S P E E D , however, does not take advantage of multiple coils like parallel imaging does. In the same way that S P E E D is extended to S P E E D - A C E , Simplified S P E E D can be developed with Array C o i l Enhancement (Simplified S P E E D - A C E ) to highly accelerate MRA. Simplified S P E E D - A C E can also be considered as a simplified version of S P E E D - A C E that applies specifically to M R A data. Unlike S P E E D - A C E , Simplified S P E E D - A C E uses neither differential filtering nor the assumption of slowly varying sensitivity profiles, and thus can achieve more efficiency in reducing scan time and more accuracy in deghosting. A s a general and independent fast M R A method, Simplified 195 Chapter 8 Conclusions and Recommendations for Future Work S P E E D - A C E can also be extended to 3 D or be combined with other existing techniques such as Half-Fourier imaging for improved performance. 8.2.3 Applications of SPEED in Multiple Acquisitions Another interesting project in my future research is to accelerate multiple acquisitions, which are rather time consuming and are required in many M R I applications such as P C M R A with double acquisitions, three-point water-fat imaging and diffusion imaging. In particular, seven acquisitions are often performed in diffusion tensor imaging [12]. Given their strong similarities in structures, multiple acquisitions generally have a certain amount of redundancy, which is similar to what occurs in dynamic M R I . Dynamic M R I is used to capture a moving object by acquiring a series of images at a high temporal rate [91]. Since the image series generally share most of structures in common, there is a certain amount of redundancy in the data. M a n y fast imaging methods has been proposed to accelerate dynamic M R I by minimizing the redundancy, which include U N F O L D [92], T R I C K S [76], and k-t S E N S E [91]. Similarly, by minimizing redundancy within multiple acquisitions, S P E E D may be used to further reduce scan time. of the idea has been demonstrated with preliminary results. 196 The feasibility Chapter 8 Conclusions and Recommendations for Future Work 8.2.4 S N R Improvement b y U s i n g S P E E D As discussed in Chapter 7, we have found that S N R loss caused by k-space undersampling can be mostly recovered in S P E E D reconstruction. The question arises whether S P E E D reconstruction can be used to improve S N R instead of accelerating acquisition. More specifically, given a full k-space data set, can the S P E E D algorithm help reconstruct an image with a better S N R than that obtained through standard 2 D F T reconstruction? These questions point to an interesting direction for future research. If the future project leads to positive results, it would imply a significant impact on M R I : perhaps most of M R I images should be reconstructed with S P E E D algorithm. This idea has been studied with computer simulations; the preliminary results demonstrate that, depending on the image contents, it is feasible to use the S P E E D algorithm to improve SNR. 197 Bibliography Bibliography [I] C . B . A h n and Z . H . Cho, " A new phase correction method i n N M R imaging based on autocorrelation and histogram analysis," IEEE Trans. Med. Imag., MI-6:32-36, 1987. [2] Z . Chang and Q.S. Xiang, "Non-linear Phase Correction with an Extended Statistical Algorithm," IEEE Trans. Med. Imag., 24: 791-798, 2005. [3] Q.S. Xiang, "Accelerating M R I by skipped phase encoding and edge deghosting ( S P E E D ) , " Magn. Reson. Med., 53:1112-1117, 2005. [4] Z . Chang and Q.S. Xiang, "Accelerating M R A with simplified skipped phase encoding and edge deghosting ( S P E E D ) , " In Proc. 13th Annu. Meet. ISMRM, 2005, p. 2318. [5] Z . Chang and Q.S. Xiang, "Highly accelerated M R I by sensitivity encoding with skipped phase encoding and edge deghosting ( S E N - S P E E D ) , " In Proc. 13th Annu. Meet. ISMRM, 2005, p. 2448. [6] F . Bloch, W . W . Hansen, and M . Packard, "Nuclear Induction," Phys. Rev. 69:127, 1946. [7] F. Bloch, "Nuclear Induction," Phys. Rev. 70: 460-474, 1946. [8] E . M . Purcell, H . C . Torrey, and R . V . Pound, "Resonance absorption by nuclear magnetic moments in a solid," Phys. Rev. 69:37-38, 1946. [9] P . C . Lauterbur, "Image formation by induced local interactions: examples employing nuclear magnetic resonance," Nature 242:190-191, 1973. [10] P. Mansfield, P . K . Grannell, " N M R 'diffraction' i n solids?" J. Phys. C: Solid Phys., 6:LA22-IA21, State 1973. [II] Q.S. Xiang, Class Notes of U B C Phys542 "Nuclear Magnetic Resonance Imaging", 2002. [12] M . A . Bernstein, K . F . K i n g , and X . J . Zhou, Handbook Elsevier/Wiley, 2004 198 of MRI Pulse Sequences, Bibliography [13] David Pickens, "Magnetic Resonance Imaging," In: J. Beutel, H . L . Kundel, R . L . V a n Metter, eds, Handbook of Medical Imaging, Volume 1: Physics and Psychophysics, SPIE, 2000. [14] D . B . Twieg, "The k-trajectory formulation of the N M R imaging process with applications in analysis and synthesis of imaging methods," Med. Phys., 10:610-621, 1983. [15] E . L . Hahn, "Spin echoes," Phys Rev, 80: 580-594, 1950. [16] R . L . D i x o n , K . E . Ekstrand, "The physics of proton N M R , " Med. Phys., 9:807-818, 1982. [17] G . M . Bydder, I.R. Young, " M R imaging: clinical use of the inversion recovery sequence," J. Comput. Assist. Tomogr., 9(4):659-675,1985. [18] P. Mansfield, "Multi-planar image formation using N M R spin-echoes," J. Phys. C: Solid State Phys., 10:L55-L58, 1977. [19] I. R. Young and G . M . Bydder, "Phase imaging," In: D . D . Stark and W . G . Bradley, eds, Magnetic resonance imaging. 2nd ed. St. Louis, M O : M o s b y Year Book, 1992. [20] S. Chavez, Q.S. Xiang, and L . A n , "Understanding phase maps in M R I : A new cutline phase unwrapping method," IEEE Trans. Med. Imag., v o l . 21, pp. 966-977, August 2002. [21] J.R. Singer, " B l o o d flow rates by nuclear magnetic resonance measurements," Science, 130, 1652-1653, 1959. [22] E . L . Hahn, "Detection of sea-water motion by nuclear procession," Journal Geophysical Research, 65, 776-777, 1960. [23] T. Grover and J.R. Singer, " N M R spin-echo flow measurements," J. Appl. 42,938-940,1971. of Phys., [24] P.R. Moran, " A flow velocity zeugmatographic interface for N M R imaging in humans," Magn. Reson. Imag., 1:197-203, 1982. [25] R . K . Harris, Nuclear Magnetic Resonance Spectroscopy, LondomPitman, 1983 [26] T . R . Brown, M . B . Kincaid, K Ugurbil, " N M R chemical shift imaging in three dimensions." In Proc Natl Acad Sci USA, 79:3523-3526, 1982. 199 Bibliography [27] A . A . Maudsley, S . K . H i l a l , W . H . Perman, and H . E . Simon, "Spatially Resolved High Resolution Spectroscopy by 'Four-Dimensional' N M R . " / . Magn. Reson., 51:147152,1983. [28] L L . Pykett and B . R . Rosen, "Nuclear magnetic resonance: i n vivo proton chemical shift imaging." Radiology, 149:197-201, 1983. [29] R . E . Sepponen, J.T. Sipponen, and J.I.Tanttu, " A method for chemical shift imaging: demonstration of bone marrow involvement with proton shift imaging." J. Comput. Assist. Tomogr. 8:585-587, 1984. [30] Q.S. X i a n g and L . A n , "Water-fat imaging with direct phase encoding," / . Magn. Reson. Imag., 7:1002-1015, 1997. [31] W . T . D i x o n , "Simple proton spectroscopic imaging," Radiology, 153:189-194, 1984. [32] C . B . A h n , S . Y . Lee, O. Nalcioglu, and Z . H . Cho, "Spectroscopic imaging by quadrature modulated echo time shifting," Magn. Reson. Imag., 4: 110-111, 1986. [33] L . A n and Q.S. Xiang, "Chemical shift imaging with spectrum modeling," Magn. Reson. Med., 46:126-130, 2001. [34] C C . Lodes, J.P. Felmlee, R . L . Ehman, et al, "Proton chemical shift imaging using double and triple phase contrast acquisition methods," J. Comput. Assist. Tomogr., 13:855-861, 1989. [35] C . B . A h n , J . H . K i m , and Z . H . Cho, "High-speed spiral-scan echo planar N M R imaging," IEEE Trans. Med. Imag., M I - 5 : 2-7, 1986. [36] A . Haase, J . Frahm, D . Matthaei, W . Hanicke, and K . D . Merboldt, " F L A S H imaging: rapid N M R imaging using low flip-angle pulses," / . Magn. Reson., 67: 258-266, 1986. [37] J. Hennig, A . Nauerth, and H . Friedburg, " R A R E imaging: A Fast imaging method for clinical M R , " Magn. Reson. Med., 3: 823-833, 1986. [38] K . Oshio and D . A . Feinberg, " G R A S E M R imaging: A novel fast M R I technique," Magn. Reson. Med., 20:344-349, 1991. [39] P . M . Margosian, G . DeMeester, and H . L i u , "Partial Fourier acquisition in M R I , " In: D . M . Grant and R . K . Harris, eds. Encyclopedia of Nuclear Magnetic Resonance, John W i l e y & Sons, 5: 3463-3467, 1996. 200 Bibliography [40] R . T . Constable and R . M . Henkelman, "Data extrapolation for truncation artefact removal," Magn. Reson. Med., 17: 108-118, 1991. [41] Q.S. X i a n g and R . M . Henkelman, "Dynamic image reconstruction (DIR): M R movies from motion ghosts," J. Magn. Reson. Imag., 2: 679-685, 1992. [42] J.J. van Vaals, M . Brummer, and T . W . Dixon, H . H . Tuithof, H . Engels, R . C . Nelson B . M . Gerety, J . L . Chezmar, and J . A . den Boer, " ' K e y h o l e ' method for accelerating imaging of contrast agent uptake," J. Magn. Reson. Imag., 3: 671-675, 1993. [43] X . H u and T. Parrish, "Reduction of field of view for dynamic imaging," Magn. Reson. Med., 31: 691-694, 1994. [44] J . W . Carlson, " A n algorithm for N M R imaging reconstruction based on multiple R F receiver coils," J. Magn. Reson., 74:376-380, 1987. [45] M . Hutchinson and U . Raff, "Fast M R I data acquisition using multiple detectors," Magn. Reson. Med, 6: 87-91, 1988. [46] D . Kwiat, S. Einav, G . Navon, " A decoupled coil detector array for fast image acquisition in magnetic resonance imaging," Med. Phys., 18: 251-265, 1991. [47] J . B . R a and C . Y . R i m , "Fast imaging using subencoding data sets from multiple detectors," Magn. Reson. Med, 30:142-145, 1993. [48] D . K . Sodickson and W . J . Manning, "Simultaneous acquisition of spatial harmonics ( S M A S H ) : Ultra-fast imaging with radiofrequency coil arrays," Magn. Reson. Med., 38: 591-603,1997. [49] K . P . Pruessmann, W . Weiger, M . B . Scheidegger, P. Boesiger, " S E N S E : Sensitivity encoding for fast M R I , " Magn. Reson. Med, 42: 953-962, 1999. [50] W . E . Kyriakos, L . P . Panych, D . F . Kacher, C F . Westin, S . M . Bao, R . V . Mulkern, and F . A . Jolesz, "Sensitivity profiles from an array of coils for encoding and reconstruction in parallel ( S P A C E RIP)," Magn. Reson. Med, 44: 301-308, 2000. [51] M . A . Griswold, P . M . Jakob, M . Nittka, J.W. Goldfarb,and A . Haase, "Partial parallel imaging with localized sensitivities (PELS)," Magn. Reson. Med., 44: 602-609, 2000. [52] M . A . Griswold, P . M . Jakob, R . M . Heidemann, M . Nittka, V . Jellus, J. Wang, B . Kiefer, and A . Haase, "Generalized autocalibrating partially parallel acquisitions ( G R A P P A ) , " Magn. Reson. Med., 47: 1202-1210, 2002. 201 Bibliography [53] D . C . Ghiglia and M . D . Pritt, Two-Dimensional Phase Algorithms, and Software. 1st ed. N e w Y o r k : Wiley, 1998. Unwrapping: Theory, [54] J. M . Huntley, "Noise-immune phase unwrapping algorithm," Appl. Opt.—Letters the Editor, 28:3268-3270, 1989. to [55] D . J . Bone, "Fourier fringe analysis: The two-dimensional phase unwrapping problem," Appl. Opt., 30: 3627-3632, 1991. [56] C . DeVeuster, P. Slangen, Y . Renotte, L . Berwart, and Y . L i o n , "Disk-growing algorithm for phase-map unwrapping: Application to speckle interferograms," Appl. Opt., 35: 240-247, 1996. [57] M . Takeda and T. A b e , "Phase unwrapping by a maximum cross-amplitude spanning tree algorithm: A comparative study," Opt. Eng. 35, pp. 2345-2351, 1996. [58] L . A n , Q.S. Xiang, and S. Chavez, " A fast implementation of the minimum spanning tree method for phase unwrapping," IEEE Trans. Med. Imag., 19:805-808, 2000. [59] T.J. Flynn, "Two-dimensional phase unwrapping discontinuity," J. Opt. Soc. Amer. A, 14: 2692-2701, 1997. with minimum weighted [60] S. M . Song, S. Napel, N . J. Pelc, and G . H . Glover, "Phase unwrapping of M R phase images using Poisson equation," IEEE Trans. Image Processing, 4: 667-676, 1995. [61] D . C . Ghiglia and L . A . Romero, "Minimum-norm two-dimensional unwrapping,"/. Opt. Soc. Amer. A, 13: 1999-2013, 1996. phase [62] J.R. M a c F a l l , N . J . Pelc, and R . M . V a v r e k , " Correction of spatially dependent phase shifts for partial Fourier imaging," Magn. Reson. Imag.,. 6: 143-155, 1988. [63] G . H . McGibney, M . R . Smith, and S.T. Nichols, "Estimation of phase errors in magnetic resonance imaging," In Proc. 9th Annu. Meet. Soc. Magn. Reson. Med., 1990, p.563. [64] J . Hua and G . C . Hurst, "Noise and artifact comparison for Fourier and polynomial phase correction used with Fourier reconstruction of asymmetric data sets," / . Magn. Reson. Imag., 2: 347-353, 1992. [65] E . Schneider and G . H . Glover, "Rapid in vivo proton shimming," Magn. Med., 18: 335-347, 1991. 202 Reson. Bibliography [66] Z . P . Liang, " A model-based method for phase unwrapping," IEEE Trans. Med. Imag., 15: 893-897, 1996. [67] M . Hedley and D . Rosenfeld, " A new two-dimensional phase unwrapping algorithm for M R I images," Magn. Reson. Med., 24:177-181, 1992. [68] Pratt, Digital Image Processing, New Y o r k : Wiley, 1978, p. 132. [69] R . M . Henkelman, "Measurement of signal intensities in the presence of noise i n M R images," Med. Phys., 12: 232-233, 1985. [70] J . A . Derbyshire, Y . P . D u , and M . Saranathan, "The Fourier conjugate of the A h n algorithm with application for real-time signal procession for M R I , " In Proc. 7th Annu. Meet. ISMRM, Philadelphia, 1999. p. 2002. [71] R . N . Bracewell, The Fourier M c G r a w - H i l l , 1978. transform and its application, New York, N Y , [72] Q.S. Xiang, "Inversion recovery image reconstruction with multiseed regiongrowing spin reversal," J. Magn. Reson. Imag., 6: 775-782, 1996. [73] M . R . Prince, E . Y u c e l , J. Kaufman, D . Harrison, and S. Geller, "Dynamic gadolinium-enhanced three-dimensional abdominal M R arteriography," J. Magn. Reson. Imag., 3:877-881, 1993. [74] D . C . Peters, F . R . Korosec, T . M . Grist, W . F . Block, K . K . Vigen, J.E. Holden, and C . A . Mistretta, "Undersampled projection reconstruction applied to M R angiography," Magn Reson Med., 43:91-101, 2000. [75] J. D u , F . R . Korosec, F . J . Thornton, T . M . Grist, and C . A . Mistretta, " H i g h Resolution Multistation Peripheral M R Angiography Using Undersampled Projection Reconstruction Imaging," Magn. Reson. Med., 52:204-208, 2004. [76] F . R . Korosec, R. Frayne, T . M . Grist, and C . A . Mistretta, "Time-resolved contrast enhanced 3 D M R angiography," Magn. Reson. Med, 36:345-351, 1996. [77] K . K . Vigen, D . C . Peters, T . M . Grist, W . F . Block, and C . A . Mistretta, "Undersampled projection-reconstruction imaging for time-resolved contrast-enhanced imaging," Magn. Reson. Med., 43:170-176, 2000. [78] A . V . Barger, W . F . Block, Y . Toropov, T . M . Grist, and C . A . Mistretta, "Timeresolved contrast-enhanced imaging with isotropic resolution and broad coverage using an undersampled 3 D projection trajectory," Magn. Reson. Med., 48:297-305, 2002. 203 Bibliography [79] D . K . Sodickson, C . A . M c K e n z i e , W . L i , S. Wolff, W . J . Manning, and R . R . Edelman. "Contrast-enhanced 3D M R angiography with simultaneous acquisition of spatial harmonics: a pilot study," Radiology, 217:284-289, 2000. [80] M . Weiger, K . P . Pruessmann, A . Kassner, G . Roditi, T. Lawton, A . Reid, P. Boesiger, "Contrast-enhanced 3D M R A using S E N S E , " J. Magn. Reson. Imag., 12:671— 677,2000 [81] M . E . Huber, S. Kozerke, K . P . Pruessmann, J. Smink, and P. Boesiger, "SensitivityEncoded Coronary M R A at 3T," Magn. Reson. Med., 52: 221-227, 2004. [82] G . A . Laub and W . A . Kaiser, " M R angiography with gradient motion refocusing," J. Comput. Assist. Tomogr., 12:377-382, 1988. [83] Q.S. X i a n g and R . M . Henkelman, "K-space description for the M R imaging of dynamic objects," Magn. Reson. Med., 29: 422-428, 1993. [84] Q.S. X i a n g and R . M . Henkelman, " M o t i o n artifact reduction with three-point ghost phase cancellation," J. Magn. Reson. Imag., 1: 633-642, 1991. [85] D . L . Parker, J.S. Tsuruda, K . C . Goodrich, A . L . Alexander, and H . R . Buswell, "Contrast en-hanced M R angiography of cerebral arteries," Investigative Radiology, 33:560-572, 1998. [86] M . R . Prince, T . M . Grist, and J.F. Debatin, 3D contrast MR angiography. Springer, 1999. Berlin: [87] F. Wiesinger, P. Boesiger, and K . P. Pruessmann, "Electrodynamics and ultimate S N R i n parallel M R imaging," Magn. Reson. Med., 52: 376-390, 2004. [88] M . A . Ohliger, A . K . Grant, and D . K . Sodickson, "Ultimate intrinsic signal-to-noise ratio for parallel M R I : electromagnetic field considerations," Magn. Reson. Med., 50: 1018-1030, 2003. [89] W . A . Edelstein, J . M . S . Hutchison, G . Johnson, and T. Redpath, "Spin warp N M R imaging and applications to human whole-body imaging," Phys. Med. Biol., 25: 751-756, 1980. [90] L . Shepp and B . Logan, "The Fourier reconstruction of a head section," IEEE Trans. Nuc. Sci., NS-21 (6): 21-43, 1974. [91] J . Tsao, P . Boesiger, and K . P . Pruessmann, " k-t B L A S T and k-t S E N S E : dynamic M R I with high frame rate exploiting spatiotemporal correlations," Magn. Reson. Med., 50: 1031-1042, 2003. 204 Bibliography [92] B . Madore, G . H . Glover, and N . J. Pelc, "Unaliasing by Fourier-encoding the overlaps using the temporal dimension ( U N F O L D ) , applied to cardiac imaging and f M R I , " Magn. Reson. Med., 42: 813-828, 1999. [93] J.W. Lankhaar, M . B . M . Hofman, J.T. Marcus, J . J . M . Zwanenburg, T.J.C. Faes, and A . Vonk-Noordegraaf, "Correction of phase offset errors in main pulmonary artery flow quantification," J. Magn. Reson. Imag., 22:73-79, 2005. 205
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Study of phase information in MRI with applications...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Study of phase information in MRI with applications to fast imaging and inversion recovery imaging Chang, Zheng 2005
pdf
Page Metadata
Item Metadata
Title | Study of phase information in MRI with applications to fast imaging and inversion recovery imaging |
Creator |
Chang, Zheng |
Date Issued | 2005 |
Description | In Magnetic Resonance Imaging (MRI), signal is complex valued and can be represented by a magnitude image and a phase map. Although magnitude images are used much more frequently in clinical diagnosis than phase maps, the latter should not be undervalued because a lot of valuable information is encoded only in phase, which has great potential in medical applications. My doctoral thesis focuses on phase information in MRI and its applications in the areas of fast imaging and Inversion Recovery (IR) imaging. IR imaging is one of the most useful techniques in MRI contrast manipulation, but its use is limited because of the presence of phase errors. The thesis proposes a new method to correct phase errors that occur in IR imaging. The method models phase variation with a polynomial, whose coefficients are statistically determined by calculating relative vector rotations of complex fields after being shifted by n pixels. As discovered in this thesis, increasing the pixel shift effectively enhances phase signal without amplifying the corresponding noise, and thereby improves phase correction. The method has been successfully demonstrated with 2D in vivo IR imaging data. Phase information has great potential also in fast imaging in MRI. For example, a recently published fast imaging method named "Skipped Phase Encoding and Edge Deghosting (SPEED)" conducts strategic spatial phase encoding and thereby accelerates MRI. SPEED is promising but is demonstrated so far only in 2D and only with a single coil. This thesis presents new developments based on the principle of SPEED: First, SPEED is extended from 2D to 3D to reduce scan time with more flexibility and efficiency; Second, SPEED is combined with Half-Fourier imaging to accelerate MRI further by a factor of nearly 2; Third, SPEED is simplified based on the sparseness of M R angiography data; Fourth, a new parallel imaging strategy named "SPEED with Array Coil Enhancement (SPEED-ACE)" is proposed, which extends SPEED from a single coil to multiple coils to further accelerate MRI. Finally, signal-to-noise ratio (SNR) in SPEED is studied, and a novel approach for SNR improvement is proposed and demonstrated with computer simulations and in vivo data. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085585 |
URI | http://hdl.handle.net/2429/18231 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2006-129808.pdf [ 26.4MB ]
- Metadata
- JSON: 831-1.0085585.json
- JSON-LD: 831-1.0085585-ld.json
- RDF/XML (Pretty): 831-1.0085585-rdf.xml
- RDF/JSON: 831-1.0085585-rdf.json
- Turtle: 831-1.0085585-turtle.txt
- N-Triples: 831-1.0085585-rdf-ntriples.txt
- Original Record: 831-1.0085585-source.json
- Full Text
- 831-1.0085585-fulltext.txt
- Citation
- 831-1.0085585.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085585/manifest