University of Southern California Dissertations and Theses

New theory and methods for accelerated MRI reconstruction

(USC Thesis Other)

New theory and methods for accelerated MRI reconstruction

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NEW THEORY AND METHODS FOR ACCELERATED MRI RECONSTRUCTION by Rodrigo A. Lobos A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (ELECTRICAL ENGINEERING) December 2022 Copyright 2023 Rodrigo A. Lobos To my sister. ii Acknowledgements This dissertation corresponds to the end of one of the most challenging, enriching, and fascinating experiences that I have had the chance to live, and it has been only possible due to the immense support and love that I have received from outstanding people over the last seven years. I would like to thank the members of my doctoral committee { Professors Richard Leahy, John Wood, Krishna Nayak, and Antonio Ortega { for their guidance, time, and eort. I am truly grateful to my advisor Professor Justin Haldar for his men- toring and constant support. Thank you Justin for inspiring me every day to keep growing as a scientist. I am also thankful for my labmates Tae Hyung Kim, Daeun Kim, Divya Varayanan, Yunsong Liu, Chin-Cheng Chan, and Jiayang Wang. Thank you for reminding me every day how great and rewarding is working on science. I would like to specially thank Eduardo Pavez, Karel Mundnich, and Omid Sani for being by my side during all these years (we made it!). I will be forever grateful to Richard Kittlaus, Hayden Buckley, Alexandra Guardado, Leslie Ramirez, Yecenia De La Rosa, Abdallah AlMutairi, and all the incredible people at Crosst Muse for reinforcing my consistency, discipline, persistence, and iii sense of community during the last years of my PhD. You were all a fundamental pillar in the last part of this process. I would like to acknowledge the generous nancial support provided by the Na- tional Institute of Health and the Ming Hsieh Department of Electrical Engineering for my graduate studies. Finally, I am immensely grateful to my family for their unconditional and innite love over these years. iv Table of Contents Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Main Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Organization of the Document . . . . . . . . . . . . . . . . . . . . . 5 Chapter 2: Navigator-free EPI Ghost Correction with Structured Low-Rank Matrix Models: New Theory and Methods . . . . . . . 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Background and Notation . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Analysis of SLM EPI Ghost Correction . . . . . . . . . . . . . . . . . 19 2.4 Constrained SLM Formulations . . . . . . . . . . . . . . . . . . . . . 25 2.4.1 Formulation using SENSE Constraints . . . . . . . . . . . . . . 25 2.4.2 Formulation using AC-LORAKS Constraints . . . . . . . . . . 28 v 2.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.6 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 48 2.7 Supplementary Information . . . . . . . . . . . . . . . . . . . . . . . 51 2.7.1 Proof of Theorem 1 . . . . . . . . . . . . . . . . . . . . . . . . 51 2.7.2 Proof of Corollary 1 . . . . . . . . . . . . . . . . . . . . . . . . 55 2.7.3 Proof of Corollary 2 . . . . . . . . . . . . . . . . . . . . . . . . 55 2.7.4 Majorize-Minimize Algorithm . . . . . . . . . . . . . . . . . . 56 Chapter 3: Robust Autocalibrated Structured Low-Rank EPI Ghost Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.2.1 Background: Structured Low-Rank EPI ghost correction . . . 62 3.2.2 RAC-LORAKS . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.3.1 Datasets used for Evaluation . . . . . . . . . . . . . . . . . . . 69 3.3.1.1 Gradient-Echo EPI Brain data . . . . . . . . . . . . . 69 3.3.1.2 Diusion-encoded EPI Brain Data . . . . . . . . . . . 71 3.3.1.3 Cardiac EPI Data . . . . . . . . . . . . . . . . . . . . 72 3.3.2 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.3.3 Data Processing . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 vi 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Chapter 4: On the Shape of Convolution Kernels in MRI Recon- struction: Rectangles versus Ellipsoids . . . . . . . . . . . . . . . . 98 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Chapter 5: New Theory and Faster Computations for Subspace- Based Sensitivity Map Estimation in Multichannel MRI . . . . . 114 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.2 Characteristics of the Inverse Problem . . . . . . . . . . . . . . . . . 119 5.3 Nullspace/Linear Predictability Theory for Subspace-Based Sensitiv- ity Map Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.3.1 Summary of Linear Predictability for Multichannel MRI . . . . 121 5.3.2 Linear Predictability and Sensitivity Map Estimation . . . . . 127 5.3.3 Empirical Demonstration of Nullspace Characteristics . . . . . 130 5.3.4 Naive Algorithm for Nullspace-Based Sensitivity Map Estimation132 5.3.5 Memory Requirements of the Nullspace Algorithm and the Di- rect Calculation of G(x) . . . . . . . . . . . . . . . . . . . . . 134 5.3.6 Comparisons of the Nullspace Algorithm and ESPIRiT . . . . 136 5.4 Proposed Computational Methods . . . . . . . . . . . . . . . . . . . 142 vii 5.4.1 Ellipsoidal versus Rectangular Convolution Kernels . . . . . . 143 5.4.2 FFT-Based Calculation of Nullspace Vectors . . . . . . . . . . 144 5.4.3 Smoothness-Based Interpolation . . . . . . . . . . . . . . . . . 146 5.4.4 Power Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . 149 5.4.5 Combined Approach . . . . . . . . . . . . . . . . . . . . . . . . 154 5.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . 157 Chapter 6: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 viii List of Tables 2.1 NRMSE values for the images shown in Fig. 2.8. . . . . . . . 35 2.2 NRMSE values for the images shown in Fig. 2.13. . . . . . . 38 3.1 NRMSEs for the multi-channel simulation results shown in Fig. 3.8. For each acceleration factor, the smallest values are highlighted in bold. 84 3.2 NRMSEs for the multi-channel inverted contrast simulation results shown in Fig. 3.12 . For each acceleration factor, the smallest values are highlighted in bold. . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.3 NRMSEs for the single-channel simulation results shown in Fig. 3.13. For each acceleration factor, the smallest values are highlighted in bold. 89 5.1 Normalized Projection Residual with 3232 calibration data. . . . . 140 5.2 Computation Time (secs.) with 3232 calibration data. . . . . . . . 140 5.3 Memory Usage (GB) with 3232 calibration data. . . . . . . . . . . 141 ix List of Figures 2.1 (a) EPI magnitude (top) and phase (bottom) images corresponding to (left) RO + data and (right) RO data. These alias-free images correspond to real data obtained with the PLACE method [123]. (b) Images corresponding to the same data from (a), except that the odd k-space lines for RO + and the even k-space lines for RO have been multiplied by -1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2 Plots of the singular values for the LORAKS matrices from Eq. (2.3) for the k-space datasets from Figs. 2.1(a) and (b). . . . . . . . . . . 23 2.3 Letting k 1 and k 2 denote the k-space data for the images in Figs. 1(a) and 1(b), respectively, we plot the cost function valueL C k 1 + (1)k 2 as a function of . Setting = 0 yields the cost function value for the images from Fig. 2.1(a), setting = 1 yields the cost function value for the images from Fig. 2.1(b), while setting = 0:5 yields the cost function value for the zero-lled solution. Results are shown for dierent (a) convex and (b-d) nonconvex choices of J(). The rank parameter r = 40 has been used in the nonconvex cases shown in (b-d), both without (a,b) and with (c,d) additional constraints. . . . 24 x 2.4 Comparison of dierent reconstruction techniques using prospectively undersampled (except the R = 5 case, which is retrospectively un- dersampled) in vivo single-shot EPI data at dierent parallel imaging acceleration factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Comparison between (a) DPG and (b) LORAKS with AC-LORAKS constraints for the real single-shot EPI in vivo brain data from Fig. 2.4. Instead of showing coil-combined images, a single representative chan- nel is shown to avoid contamination of the phase characteristics in- duced by coil combination. . . . . . . . . . . . . . . . . . . . . . . . 32 2.6 (top) Magnitude and (bottom) phase images corresponding to recon- struction of unaccelerated (R = 1) single-channel single-shot EPI phantom data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 xi 2.7 Phase dierences between RO + and RO for a representative subset of coils from the prospectively undersampled data shown in Fig. 4. In cases with accelerated data (where no gold standard is available), we show the phase dierence obtained after AC-LORAKS reconstruction. Unlike most other phase images shown in this paper (which do not include background masking and show the entire phase range from to), we have taken steps to make the phase images in this gure easier to visualize. In particular, we have masked the background noise to make it easier to focus on the signal structure. In addition, we have shown a restricted phase range for both the gold standard (black = 0.95 radians, white = 1.02 radians) and the reconstructions using AC-LORAKS (black = 0.89 radians, white = 1.21 radians). . . . 36 2.8 Comparison of dierent reconstruction techniques using retrospec- tively undersampled in vivo single-shot EPI data to simulate dierent parallel imaging acceleration factors. . . . . . . . . . . . . . . . . . . 37 2.9 (top) Magnitude and (bottom) phase images corresponding to recon- struction of unaccelerated (R = 1) single-channel single-shot EPI in vivo human brain data. . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.10 Comparison of dierent reconstruction techniques using prospectively undersampled in vivo double-oblique single-shot EPI data at dierent parallel imaging acceleration factors. . . . . . . . . . . . . . . . . . 39 xii 2.11 Phase dierences between RO + and RO for a representative subset of coils from the prospectively undersampled double-oblique data shown in Fig. 2.10. In cases with accelerated data (where no gold standard is available), we show the phase dierence obtained after reconstruction using LORAKS with AC-LORAKS constraints. Unlike most other phase images shown in this paper (which do not include background masking and show the entire phase range from to), we have taken steps to make the phase images in this gure easier to visualize. In particular, we have masked the background noise to make it easier to focus on the signal structure. In addition, we have shown a restricted phase range for both the gold standard (black = 0.95 radians, white = 1:02 radians) and the reconstructions using AC-LORAKS (black = 0.86 radians, white = 1.18 radians). . . . . . . . . . . . . . . . . 40 2.12 Comparison between (a) DPG and (b) LORAKS with AC-LORAKS constraints for the real single-shot EPI in vivo brain double-oblique data from Fig. 2.10. Instead of showing coil-combined images, a single representative channel is shown to avoid contamination of the phase characteristics induced by coil combination. . . . . . . . . . . . . . . 41 2.13 Comparison of dierent reconstruction techniques using retrospec- tively undersampled in vivo double-oblique data to simulate single- shot EPI experiments at dierent parallel imaging acceleration factors. 42 xiii 2.14 (top-left) Magnitude and (bottom-left) phase images corresponding to reconstruction of unaccelerated (R = 1) in vivo axial multi-channel single-shot EPI data acquired with a small FOV, and with extra simu- lated nonlinear phase dierences between RO + and RO . The images on the right show the phase dierences that were used between RO + and RO , which were set dierently for the ACS data (used for calibra- tion) than they were for the gold standard (used for reconstruction). 43 2.15 Evaluation with multi-shot EPI data for a phantom with simulated respiratory eects. A representative set of four segmented images is shown, extracted from a longer acquisition spanning several minutes. (a) Phase images corresponding to one channel of the data, with im- ages reconstructed without compensating for the mismatches between dierent shots and dierent gradient polarities. (b) Phase images for one channel of the DPG reconstruction, with reconstruction performed using multi-channel data. (c) Phase images for single-channel LO- RAKS reconstruction with AC-LORAKS constraints. (d) Plot show- ing the relative respiratory position across EPI shots (TR=60msec), as measured with an ultrasound transducer coupled to a respiratory phantom [41]. The line-plot peaks show points where the phantom air bag is maximally in ated. The sampling times for the ACS data and for the shots used to generate each of the two-shot images from (a)-(c) are marked as labeled in the legend. . . . . . . . . . . . . . . 44 xiv 2.16 Images reconstructed using LORAKS with AC-LORAKS constraints for accelerated (R = 2) data in both single-channel and multi-channel contexts. (a) Phantom dataset with simulated respiration. (b) In vivo human brain dataset. Phase images from a single channel are shown for zero-lled data reconstructed without compensating for mismatches between the gradient polarities, LORAKS reconstruc- tion with AC-LORAKS constraints from single-channel data, and LORAKS reconstruction with AC-LORAKS constraints from multi- channel data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.17 (top) Magnitude and (bottom) phase images corresponding to re- construction of unaccelerated (R = 1) in vivo double-oblique single- channel single-shot EPI data. . . . . . . . . . . . . . . . . . . . . . . 46 2.18 (top, middle) Magnitude and (bottom) phase images corresponding to reconstruction of unaccelerated (R = 1) in vivo double-oblique single-channel single-shot EPI with strong simulated phase dierences between RO + and RO . The middle row shows a dierent windowing of the magnitude images to more clearly visualize ghost-artifacts. The top-right image shows the simulated phase dierence between RO + and RO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 xv 3.1 Illustration of EPI ghost correction. The top row of this gure shows EPI images obtained from dierent methods, while the bottom row shows the same images with 10 intensity amplication to high- light the ghost characteristics. If EPI data is naively reconstructed without accounting for the systematic dierences between data ac- quired with positive and negative readout gradient polarities (\Uncor- rected"), then strong Nyquist ghosts appear in the image as indicated with arrows. Modern EPI techniques frequently try to eliminate these artifacts using navigator information to estimate the systematic dif- ferences between the data collected with dierent readout polarities. In the navigator-based example we show (\Navigator"), the naviga- tor information was collected using a 3-line EPI acquisition with the phase encoding gradients turned o, and the dierence between posi- tive and negative gradient polarities was modeled using constant and 1D linear phase terms. Although this approach substantially reduces Nyquist ghosts, it is common for some amount of residual ghosting to still be present in the images, particularly in cases where simple 1D phase modeling is inadequate to capture the dierences between the two gradient polarities. We also show an example of our proposed approach (\RAC-LORAKS"), which can account for more compli- cated variations between the dierent gradient polarities, and which is substantially more successful at suppressing Nyquist ghosts in this example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 xvi 3.2 Illustration of the orientation of the double-oblique gradient-echo EPI dataset. The double-oblique slices are shown in red, overlaid on a structural T1-weighted image of the same subject. The double-oblique slice used for the results in Fig. 3.7 is shown with a yellow rectangle. 71 3.3 Illustration of the EPI and ACS datasets used in simulation. The rst and second top rows show coil-combined multi-channel data for the case when the EPI and ACS data have similar and inverted contrast, respectively, while the bottom row shows representative single-channel images. We also show the interpolarity phase dierence for the coil- combined EPI data, as well as the dierence in the interpolarity phase dierence between the coil-combined EPI and ACS data. . . . . . . 75 3.4 ACS data and reconstruction results for in vivo gradient-echo EPI brain data with an axial slice orientation for dierent parallel imag- ing acceleration factors. Note that the rst four acceleration factors (R = 1-4) were acquired from one subject during a single scan ses- sion while the last two acceleration factors (R = 5; 6) were acquired from a dierent subject on a dierent day, which explains the visual discontinuity between these cases. . . . . . . . . . . . . . . . . . . . . 80 3.5 The same results shown in Fig. 3.4 , but with a 5 intensity ampli- cation to highlight the ghost characteristics. . . . . . . . . . . . . . . 81 xvii 3.6 DPG results corresponding to the same data shown in Fig. 3.4 and Fig. 3.5. The same mDPG results shown in Fig. 3.4 and Fig. 3.5 are also reproduced in this gure for reference. Note that the processing steps of DPG cause the image intensities to be mismatched from the intensities of mDPG and the other reconstruction methods, which precludes a quantitative comparison. . . . . . . . . . . . . . . . . . 82 3.7 ACS data and reconstruction results for in vivo gradient-echo EPI brain data with a double-oblique slice orientation for dierent par- allel imaging acceleration factors. For reference, we also show the interpolarity phase dierence as estimated from a coil-combined RAC- LORAKS result. The degree of phase nonlinearity is an indicator of how dicult ghost correction is expected to be. As can be seen, com- plicated 2D nonlinear phase dierences are present in many of these cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.8 Reconstruction results for the rst set of multi-channel simulations (with similar contrast between ACS and EPI data, but with a hy- perintensity added to the EPI data) with dierent parallel imaging acceleration factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.9 mDPG and DPG results corresponding to the same multi-channel simulated data from Fig. 3.8. . . . . . . . . . . . . . . . . . . . . . . 86 xviii 3.10 ESPs for the multi-channel simulation results shown in Fig. 3.8. The vertical axis of each ESP uses a consistent range to enable comparisons between dierent acceleration factors. . . . . . . . . . . . . . . . . . 87 3.11 Reconstruction results for multi-channel simulated data with dier- ent parallel imaging acceleration factors. These simulations are iden- tical to those reported in Fig. 3.8, except that the images used to generate EPI data and the images used to generate ACS data were interchanged. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.12 Reconstruction results for the second set of multi-channel simulations (with inverted contrast between ACS and EPI data) with dierent parallel imaging acceleration factors. . . . . . . . . . . . . . . . . . 89 3.13 Reconstruction results for single-channel simulated data with dierent acceleration factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.14 ACS data and reconstruction results for three representative slices from in vivo diusion brain data (R = 3). A 10 intensity ampli- cation is also shown for each slice to better highlight the ghosting characteristics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 xix 3.15 ACS data and reconstruction results for in vivo diusion EPI brain data for dierent parallel imaging acceleration factors. For improved visualization, zoomed-in versions of these results (corresponding to the spatial region marked with a yellow rectangle in the rst column and rst row) are shown in Fig. 3.16. It should be noted that the subject appears to have slightly moved between scans, so that there is not perfect correspondence between anatomical image features across dierent acceleration factors. . . . . . . . . . . . . . . . . . . . . . . 92 3.16 The same results shown in Fig. 3.15, but zoomed-in to a region of interest for improved visualization. . . . . . . . . . . . . . . . . . . . 93 3.17 ACS data and reconstruction results for unaccelerated in vivo car- diac EPI data. The two rows show the same results, but the second row has 5 intensity amplication to better highlight the ghosting characteristics. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.1 Rectangular kernels (rst row) versus ellipsoidal kernels (second row) for dierent kernel sizes. The center of the kernel is marked in red, while other locations within the support are shown in white. For each conguration, the total number of kernel coecients is indicated in yellow in the right-bottom corner. As the kernel size increases, the ellipsoidal kernels have substantially fewer coecients than the corresponding rectangular kernels. . . . . . . . . . . . . . . . . . . . 101 xx 4.2 The T2-weighted (rst row) and the T1-weighted (second row) datasets with the corresponding sampling patterns used in our experiments. We show results after coil combination (root sum-of-squares). . . . . 105 4.3 NRMSE for dierent reconstruction methods for the T2-weighted dataset as a function of kernel size and kernel shape. . . . . . . . . . . . . . 108 4.4 NRMSE for dierent reconstruction methods for the T1-weighted dataset as a function of kernel size and kernel shape. . . . . . . . . . . . . . 109 4.5 Computation times for dierent reconstruction methods for the T2- weighted dataset as a function of kernel size and kernel shape. . . . 109 4.6 Computation times for dierent reconstruction methods for the T1- weighted dataset as a function of kernel size and kernel shape. . . . . 110 4.7 Memory usage for dierent reconstruction methods for the T2-weighted dataset as a function of kernel size and kernel shape. . . . . . . . . . 110 4.8 Memory usage for dierent reconstruction methods for the T1-weighted dataset as a function of kernel size and kernel shape. . . . . . . . . . 111 4.9 Histograms of the NRMSE values obtained on the validation set by the CNN-based reconstruction for the dierent kernel shapes. . . . . 112 5.1 Depiction of the three datasets we use for illustration and validation. 131 5.2 Normalized singular value curves for the C matrices corresponding to the three datasets from Fig. 5.1. . . . . . . . . . . . . . . . . . . . . 132 xxi 5.3 Spatial maps of the 32 eigenvalues of the G(x) matrices for the Brain TSE dataset. The eigenvalues have been normalized byjj and sorted from largest to smallest, with the largest in the top right and the smallest in the bottom right. . . . . . . . . . . . . . . . . . . . . . . 132 5.4 Magnitude of sensitivity maps for 16 representative coils estimated using the nullspace-based method (rst row), ESPIRiT (second row), and the nullspace-based method combined with the PISCO techniques (third row). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.5 Quantitative comparison between the nullspace-based approach and ESPIRiT for the Brain TSE data in terms of (a) accuracy (as mea- sured through the normalized projection residual) and (b) speed (as measured through the total computation time). . . . . . . . . . . . . 137 5.6 Assessment of the nullspace-based algorithm using lters with a rect- angular (Nullspace) and an ellipsoidal support (Nullspace + S). . . . 144 5.7 Assessment of the nullspace-based algorithm calculating the nullspace vectors in Step (2) with (Nullspace + CO) and without (Nullspace) the FFT-based approach. . . . . . . . . . . . . . . . . . . . . . . . . 147 5.8 Assessment of the nullspace-based algorithm using (Nullspace + I) and not using (Nullspace) the FFT-based interpolation approach in Step (4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 xxii 5.9 Comparison between using the location-wise SVD approach (Nullspace) and the all-at-once PowerIteration-based approach (Nullspace + P) in Step (4) of the nullspace-based algorithm. . . . . . . . . . . . . . . . 152 5.10 Comparison between the nullspace-based algorithm with (Nullspace + PISCO) and without (Nullspace) the PISCO techniques. . . . . . . 153 5.11 Estimated sensitivity maps using the nullspace method with and with- out PISCO for the three considered datasets. We show the magnitude of the sensitivity maps of four representative coils. . . . . . . . . . . 156 xxiii Abstract Magnetic Resonance Imaging (MRI) is a safe, versatile, and noninvasive imaging modality that has revolutionized medicine by providing high-quality images of liv- ing tissue. Even though MRI technology has been constantly evolving since its origins in the 1970's, the time needed to acquire MRI data can still be restric- tively long in real clinical applications. This has limited MRI from reaching its full potential, reason why constant eorts have been made in order to accelerate the data-acquisition process. Most of the adopted approaches to reduce the data- acquisition time involve the individual or combined use of advanced hardware and partial data-acquisition, however, conventional approaches based on these princi- ples can be subject to negative eects in the signal-to-noise ratio, resolution, and/or the presence of undesired artifacts. In this study we propose new computational imaging methods which are able to successfully reconstruct partially acquired data obtained using advanced MRI hardware. Additionally, in this study we propose powerful computational techniques which allow dramatic improvements in compu- tational eciency when performing the reconstruction task. In the rst part of this study we propose reconstruction methods for partially acquired data which are able to automatically identify multiple linear predictability xxiv relationships present in the MRI data. For this purpose we rely on modern struc- tured low-rank modeling theory and advanced optimization techniques. Based on these tools, we make the novel observation that reference data can be used to learn additional linear predictability relationships which considerably improve re- construction performance in challenging scenarios. For instance, we show that linear relationships learned from reference data allow the reconstruction of highly accelerated data acquired using an undersampling scheme with a uniform structure. Notably, we show that these linear predictability relationships can be learned even in cases where the reference data are not pristine. In the second part of this study we propose computational techniques in order to improve the eciency of MRI reconstruction methods. We show that reconstruction methods which are based on the linearly predictable structure of the MRI data can be signicantly enhanced in terms of computational complexity by a simple modication in the model used for prediction. Specically, we show that a missing sample can be predicted using an ellipsoidal neighborhood of samples instead of a rectangular neighborhood of samples, which allows important improvements in computational time and memory usage with negligible eects on performance. Finally, we study how to improve the eciency of parallel imaging sensitivity map estimation methods. An accurate estimation of sensitivity maps is a funda- mental piece in many modern parallel imaging MRI reconstruction methods, which should also be performed eciently in order to reduce the overall reconstruction time. In this study we provide a novel theoretical description of the sensitivity xxv map estimation problem by leveraging on the linearly predictable structure of the MRI data. Then, relying on these theoretical results, we propose a powerful set of computational techniques which allow massive improvements in computational complexity when integrated to sensitivity map estimation methods based on sub- spaces. We show that widely used estimation methods can achieve approximately 100-fold acceleration in computational time and dramatic savings in memory usage without sacricing estimation accuracy. Remarkably, each of the proposed com- putational techniques can also be used individually to improve the eciency of methods in other signal processing applications beyond MRI. xxvi Chapter 1 Introduction 1.1 Motivation Magnetic Resonance Imaging (MRI) is one of the most relevant and exible medical imaging modalities in clinical medicine, which allows high-quality imaging of living tissues in a safe and noninvasive manner. Even though MRI has been constantly pushing the frontiers of diagnostic clinical medicine over the last 50 years, obtain- ing high-quality images in MRI is still subject to restrictive long data-acquisition processes. This limits scanning throughput, aects accessibility, and contributes to patient discomfort. Over the years many eorts have been made in order to decrease the acquisition time, and one popular approach has been to develop fast acquisition protocols given the constant improvements on MRI hardware. A second widely used approach, which does not involve hardware modications, corresponds to acquiring data partially. In this approach, also known as accelerated MRI acqui- sition, only part of the MRI data is acquired, and eects due to the missing data 1 are compensated in post-acquisition stages by relying on advanced optimization algorithms. Over the last decades researchers have been constantly trying to extend the limits of accelerated MRI by increasing the degree of missing data in the acquisi- tion [28, 31, 36, 38, 74, 86, 103, 112, 118, 132]. However, reconstruction methods are not always able to generate a useful image when the number of acquired samples is too small, and one way to surpass this limitation has been to include prior knowl- edge in the reconstruction process. A popular approach corresponds to leveraging autoregressive linear predictability relationships which typically exists in the MRI data [37,74]. Nevertheless, current reconstruction methods based on this approach have not been able to take full advantage of this type of relationships, and many possibilities to obtain further improvements in image-quality remain still open. Another essential aspect in the reconstruction of accelerated MRI data is com- putational complexity. In many clinical applications MRI data can achieve massive dimensions, which challenges the eciency of reconstruction methods in terms of computation time and memory usage. For instance, the potential high-dimensional characteristics of MRI accelerated data might result in long reconstruction pro- cesses, which could negatively aect clinical applications where data need to be analyzed soon afterwards the acquisition. Then, even though many reconstruction methods exhibit powerful characteristics in terms of image-quality performance, their high computational complexity tends to hinder their translation to clinical 2 applications. Therefore, the development of computational techniques to improve reconstruction eciency has been an active area of research in MRI. In this work we make novel contributions in order to solve the aforementioned problems in accelerated MRI reconstruction. These contributions are summarized in the following section. 1.2 Main Contributions • We develop novel reconstruction methods for accelerated MRI data which leverage on multiple linear predictability relationships existent within the data, and between the data and correlated reference datasets. Notably, these lin- ear relationships are learned automatically by relying on structured low-rank modeling theory. We theoretically and empirically show that the integration of linear relationships learned from correlated reference datasets allow the successful reconstruction of data acquired using challenging undersampling schemes, which typically cannot be addressed by conventional methods. Re- markably, we show that the proposed methods are robust to imperfections in the reference data, and that they can be used synergistically with fast acquisi- tion protocols based on advanced MRI hardware. Thus, the proposed methods open new opportunities when aiming at decreasing the total scanning time. • We provide novel computational techniques to improve the eciency of re- construction methods. Specically, we rely on advanced signal processing and linear predictability theory to develop computational techniques which 3 allow important improvements regarding computational time and memory us- age of conventional and state-of-the-art reconstruction methods. Remarkably, the improvements in computational complexity are achieved without compro- mising reconstruction quality. We show that the eciency of reconstruction methods based on linear predictability principles can be signicantly improved by a simple modication of the model used for prediction. We make a sys- tematic comparison between the previous and the proposed predicting model by considering several conventional and state-of-the-art methods. We show that signicant improvements in eciency are obtained in all the cases when adopting the proposed model. • We propose novel theory and fast computational methods in the context of parallel imaging sensitivity map estimation. The accurate and ecient esti- mation of sensitivity maps is a fundamental piece in many conventional and modern MRI reconstruction methods which leverage on parallel imaging con- straints. In this work we provide a novel and intuitive theoretical description of the sensitivity map estimation problem which relies on linear predictability and structured low-rank modeling theory. Then, equipped with these theoret- ical results, we propose a set of computational techniques which allow massive improvements in the computational complexity of sensitivity map estimation methods based on subspaces. Notably, we empirically show that conventional subspace-based sensitivity map estimation methods can achieve approximately 4 a 100-fold acceleration in time when combined with all the proposed compu- tational techniques. Remarkably, this eciency improvement aects the esti- mation accuracy negligibly. Moreover, the proposed computational techniques can be used combined or individually in several reconstruction methods and other important signal processing applications beyond MRI. 1.3 Organization of the Document • In Chapter 2 we describe 1 how we proposed new methods for accelerated MRI reconstruction based on leveraging autoregressive shift-invariant linear pre- dictability relationships existent within the data, and between the data and correlated datasets. Using structured low-rank modeling theory we show how correlated reference datasets allows the reconstruction of accelerated MRI data acquired using undersampling schemes with uniform structure. Based on our theoretical results, we propose two reconstruction methods that leverage on constraints learned from reference data in the spatial and in the Fourier do- main, respectively. We empirically show that using constraints in the Fourier domain of the reference data oers better performance than constraints related to the spatial domain. By using these Fourier-based and additional nonconvex low-rank constraints, we show that the proposed method successfully recon- struct MRI data in challenging data-acquisition scenarios. Specically, we 1 The work presented in this Chapter 2 can also be found in [84]. Copyright belongs to IEEE Transactions on Medical Imaging. 5 show that echo planar imaging (EPI) data can be reconstructed without the presence of undesired artifacts even when the data has been 5-fold accelerated. • In Chapter 3 we develop 2 an enhanced accelerated MRI reconstruction method based on structured low-rank modeling principles, which is able to leverage on linear predictability relationships learned from correlated reference datasets even in cases where they present imperfections. We account for these refer- ence data imperfections by proposing a novel optimization approach for the reconstruction problem. By relying on our proposed optimization techniques, we additionally show that the enhanced method is able to leverage on further constraints learned from the reference data, which allows the reconstruction of data under challenging acceleration regimes. Remarkably, we show that the enhanced method is able to reconstruct EPI data even when the data has been 6-fold accelerated. • In Chapter 4 we study 3 how to improve the eciency of MRI reconstruction methods based on linear predictability principles. We theoretically and em- pirically show the eects of adopting a predicting model which assumes an ellipsoidal shape for the set of samples involved in the prediction of missing 2 The work presented in this Chapter 3 can also be found in [82]. Copyright belongs to Magnetic Resonance in Medicine. 3 The work presented in this Chapter 4 can also be found in [81]. Copyright belongs to Magnetic Resonance in Medicine. 6 samples, instead of a rectangular shape as implemented in conventional ap- proaches. We conclude that the ellipsoidal approach oers signicant improve- ments in computational time and memory usage for several MRI reconstruction methods while preserving reconstruction quality. • In Chapter 5 we provide a novel theoretical framework to describe the parallel imaging MRI sensitivity map estimation problem by relying on linear pre- dictability and structured low-rank modeling theory. We then provide a set of powerful computational techniques which allow massive improvements in eciency while maintaining performance when integrated to subspace-based sensitivity map estimation methods. We show that approximately a 100-fold reduction in computational time is achieved when the proposed techniques are used. Finally, we provide our conclusions in Chapter 6. 7 Chapter 2 Navigator-free EPI Ghost Correction with Structured Low-Rank Matrix Models: New Theory and Methods 2.1 Introduction Echo-planar imaging (EPI) [95] is currently one of the fastest MRI pulse sequences and one of the most popular sequences for functional, diusion, and perfusion imaging. EPI uses a train of gradient echoes to measure multiple lines of k-space from a single excitation, but is prone to artifacts because it employs a long readout, uses rapidly-switching high-amplitude gradients, and measures alternating lines of k-space with dierent gradient polarities [5]. In conventional single-shot EPI, even and odd lines of k-space are acquired with alternating gradient polarities. In practice, hardware imperfections, eddy currents, eld inhomogeneity, concomitant elds, system delays, and similar phenomena can introduce signal phase errors between k-space lines acquired with dierent readout gradient polarities. If these phase errors are not correctly compensated, a Nyquist 8 (or N=2) ghost artifact is observed corresponding to an aliased image that is po- sitioned a half eld-of-view (FOV) away from the true spatial position along the phase-encoding direction. In multi-shot EPI, full k-space coverage is achieved by using multiple excitations, where a dierent segment of k-space is acquired using EPI for each shot. Multi-shot EPI is used to reduce the EPI echo train, which subsequently reduces distortion and spin-dephasing eects from local eld inho- mogeneity. Images can then be reconstructed by interleaving the multi-shot data together. As in the single-shot case, the mismatch between dierent gradient polar- ities also leads to Nyquist ghost artifacts for multi-shot data. However, multi-shot data may also exhibit additional ghost artifacts if there happen to be inconsisten- cies between each shot { including system drift, subject motion, and, particularly in the case of gradient-recalled EPI, respiration. Many approaches to ghost correction have been proposed over the years, which we group into two main categories. The rst category contains simple model- based approaches such as Refs. [12,15,21,46,51,111,125,127], which assume a low- dimensional model to describe a systematic phase mismatch between even and odd lines. The parameters of the mismatch model are often estimated using separate navigator data, and can then be used to correct the mismatch in the measured EPI data. While these methods are widely used and can work well when the mismatch model is accurate, phenomena such as eddy currents can lead to more complicated data mismatches that are not fully captured by simple models. 9 This paper focuses on the second category of methods, which includes those described in Refs. [19, 42, 44, 123, 124]. These methods rely on a more exible model in which the data samples of each gradient polarity/shot are assumed to be coming from dierent but highly-correlated images. For example, it is often assumed that the images corresponding to dierent gradient polarities or shots have the same image magnitudes but dierent image phases. 1 This is similar to how parallel imaging methods like SENSE [103] and GRAPPA [28] assume that the dierent channels of an array receiver coil acquire images that are dierent (i.e., modulated by dierent coil sensitivity proles) but highly correlated. As a result, it is not surprising that many recent ghost correction approaches can be viewed as adaptations of previous parallel imaging methods to the ghost correction context. An example is the dual-polarity GRAPPA (DPG) method [42], which treats dierent polarities as if they were dierent virtual coils, and uses a dual GRAPPA kernel (with the GRAPPA weights divided into two halves corresponding to the two dierent gradient polarities) to synthesize a ghost-free fully-sampled image. Even though methods from the second category have been shown to have state-of-the-art performance in many challenging scenarios, they can still suer from artifacts in certain cases. For example, DPG can fail to successfully correct ghost artifacts if there are mismatches between the measured EPI data and the autocalibration signal (ACS) used to train the dual GRAPPA kernel. This type 1 This assumption is a simplication of the true imaging physics, and may not fully account for any time- dependent eects that evolve dynamically during data acquisition (e.g., due to eddy currents). More detailed modeling has been considered in some previous work, e.g., [12], though the simplied image-domain model described above underlies much of the recent literature. 10 of mismatch can occur because of changes in the measured data as a function of time, e.g., due to respiration [41]. In addition, most of the methods in the second category rely on the use of multichannel data, and are not easily applicable to single-channel ghost-correction. Recently, novel image reconstruction methods have been proposed that enable calibrationless single-channel and multi-channel image reconstruction from under- sampled k-space data using structured low-rank matrix (SLM) completion ap- proaches [32,33,35,38,56,66,97,110]. SLM approaches are based on the assumption that there exist linear dependencies in k-space due to limited image support, smooth image phase variations, parallel imaging constraints, and/or transform-domain im- age sparsity. While such constraints have been used before in EPI ghost correction, e.g., Refs. [12,58], the SLM formulation of these constraints is distinct from classical approaches. SLM approaches were not originally developed for EPI data, but have very recently been adapted to such contexts [66, 71, 88, 91, 93]. These approaches have demonstrated to yield state-of-the-art performance in highly-accelerated EPI image reconstruction [66] and the ability to perform navigator-free EPI ghost cor- rection [71,88,91,93]. In this paper, we analyze theoretical aspects of navigator-free EPI ghost correc- tion using SLM approaches and obtain new insights that have major implications for ghost correction performance. Specically, we prove that the SLM completion prob- lem associated with ghost correction either has a non-unique solution or a unique solution that is undesirable. Based on this result, we observe that constraints are 11 needed to ensure the performance of SLM-based ghost correction, and investigate two approaches that achieve substantially improved results. A preliminary account of portions of this work was previously given in Ref. [80]. In our rst approach, we combine ideas from the LORAKS [32, 35, 38, 66] SLM framework with coil sensitivity maps within the SENSE framework, as has pre- viously been done for EPI reconstruction [66] and EPI ghost correction [91, 93]. Compared to MUSSELS [91,93] (a similar SENSE-based ghost correction method), our new approach makes use of a nonconvex regularization function from earlier LORAKS work [32,33,35,38,66] which yields improved results, both in theory and in practice, than the convex approach used in MUSSELS. Additionally, MUSSELS uses one of the simpler forms of SLM construction (named the C-matrix in the terminology of LORAKS [32, 35]) that incorporates support and parallel imaging constraints, but does not leverage image phase constraints. Our new SENSE-based approach takes advantage of a more advanced SLM construction (named the S- matrix in the terminology of LORAKS [32,35]). Our second approach combines LORAKS with k-space domain parallel imaging linear predictability constraints, like those used in GRAPPA [28], SPIRiT [86], and PRUNO [132]. In our implementation, these constraints are imposed within the broader framework of autocalibrated LORAKS (AC-LORAKS) [36]. To the best of our knowledge, this is the rst time that this type of information has been combined with structured low-rank matrix completion methods in the context of EPI ghost correction. This second new approach not only works for multi-channel 12 data as expected, but remarkably, we observe it also works for ghost correction of single-channel data in some cases. This paper is organized as follows. Section 2.2 reviews SLM approaches and denes the notation used in the rest of the paper. Section 2.3 presents our novel theoretical analysis of unconstrained SLM methods for EPI ghost correction. Sec- tion 2.4 describes new constrained SLM formulations that we propose to overcome the theoretical limitations of unconstrained approaches. Section 2.5 presents a systematic evaluation of these approaches with respect to current state-of-the-art approaches. Finally, discussion and conclusions are presented in Sec. 2.6. 2.2 Background and Notation While many SLM descriptions have appeared in the literature, our description of SLM will focus on the perspectives and terminology from the LORAKS framework. For simplicity, we only present a high-level review of LORAKS for the 2D case, and refer interested readers to Refs. [32,33,35,38] for more general descriptions and additional details. LORAKS is a exible constrained reconstruction framework that uses SLM mod- eling to unify and jointly impose several dierent classical and widely used MRI reconstruction constraints: limited image support constraints, smooth phase con- straints, parallel imaging constraints, and sparsity constraints [32, 33, 35, 38]. LO- RAKS is based on the observation that if any of these constraints is appropriate 13 for a given image, then the Nyquist-sampled k-space data for that image will pos- sess shift-invariant linear prediction relationships. These relationships mean that missing/corrupted data samples can be extrapolated or imputed as a weighted linear combination of neighboring points in k-space. Such linear prediction rela- tionships imply that the k-space data will lie in a low-dimensional subspace, and SLM approaches can implicitly learn and impose this subspace structure directly from undersampled/low-quality data. Importantly, while the data will possess linear prediction relationships and lie in a low-dimensional subspace under assumptions about the image support, phase, etc., the LORAKS approach is agnostic to the original source of these relationships. Instead, the approach attempts to identify and utilize all of the linear prediction relationships that are present in the k-space data, regardless of their source. This means that, while LORAKS reconstruction may be easier when support, phase, and parallel imaging constraints are applicable simultaneously, the LORAKS approach can still function when one or more of these constraints is inapplicable, as long as there are sucient sources of linear predictability in the data. The basic premise of the LORAKS support constraint [35] is that, if there are large regions of the FOV in which the true image is identically zero and if s(k x ;k y ) represents the Fourier transform of the true image, then there exist innitely many k-space functions f(k x ;k y ) such that s(k x ;k y )f(k x ;k y ) 0, where denotes the standard convolution operation. If we let k denote the vector of samples of s(k x ;k y ) on the Cartesian Nyquist grid for the FOV and let f represent the samples 14 of f(k x ;k y ) on the same Cartesian grid, then the convolution relationship can be expressed in matrix-vector form asP C (k)f 0, where the operatorP C (k) forms a Toeplitz-structured convolution matrix (called the LORAKS C-matrix [35]) out of the entries of k. Since there are many such vectors f that satisfy this relationship, we observe that the LORAKS C-matrix will be approximately low-rank. The basic premise of the LORAKS phase constraint [32,35] is that, if the image has smoothly-varying phase and the image has limited support, then there are innitely many functions h(k x ;k y ) such that s(k x ;k y )h(k x ;k y ) s(k x ;k y ) h(k x ;k y ) 0, where s(k x ;k y ) and h(k x ;k y ) are respectively the complex conjugates ofs(k x ;k y ) andh(k x ;k y ). Similar to the previous case, this convolution relationship can be expressed in a matrix-vector form as P S (k)h 0, where the operator P S (k) combines a Toeplitz-structured convolution matrix with a Hankel-structured convolution matrix (resulting in what we call the LORAKS S-matrix [32, 35]) out of the entries of k, and h is the vector of Nyquist samples ofh(k x ;k y ) and h(k x ;k y ). Since there are many such vectors h that satisfy this relationship, we observe that the LORAKS S-matrix will also be approximately low-rank. These low-rank matrix constructions are easily generalized to the context of parallel imaging. Specically, assume that data is acquired from N c channels, and let k n denote the vector of k-space samples from the nth channel, and let k tot 15 denote the vector containing the k-space samples from all channels. It has been shown that the concatenated matrix C P (k tot ) = P C (k 1 ) P C (k 2 ) P C (k Nc ) (2.1) will generally have low rank [38,110,132], and that the concatenated matrix S P (k tot ) = P S (k 1 ) P S (k 2 ) P S (k Nc ) (2.2) will also generally have low rank [38]. Note that Eqs. (2.1) and (2.2) reduce to the standard single-channel case whenN c = 1, so we will use these expressions for both the single-channel and the multi-channel cases. The observation that these matrices have approximately low rank is valuable, be- cause these low-rank characteristics can be exploited to improve image reconstruc- tion quality. Specically, by enforcing one or more of these low-rank constraints during image reconstruction, it becomes possible to reconstruct high-quality images from highly accelerated and/or unconventionally sampled k-space data. The preceding paragraphs described SLM approaches for single-channel and multi-channel image reconstruction for general contexts, but without specializa- tion to ghost correction for EPI. However, as described in the introduction, there is a straightforward analogy between parallel imaging and Nyquist ghost correction. For the sake of simplicity and without loss of generality, we will describe the SLM 16 matrix construction for this case in the context of single-shot imaging with posi- tive and negative readout polarities (denoted RO + and RO , respectively), noting that the extension to multi-shot imaging is trivial (obtained by concatenating to- gether the SLM matrices for each shot as if the dierent shots were coming from dierent receiver coils in a parallel imaging experiment). Let k + tot and k tot repre- sent hypothetical vectors of Nyquist-sampled Cartesian k-space data for the two dierent readout gradient polarities from either a single-channel or multi-channel experiment. Under the assumption that we can treat dierent readout polarities in the same way we treat dierent receiver coils in parallel imaging, we expect the matrices C P (k + tot ) C P (k tot ) (2.3) and S P (k + tot ) S P (k tot ) (2.4) to be approximately low-rank. 2 However, due to the form of single-shot EPI imag- ing, we only measure a subset of the phase encoding lines of k + tot and k tot . Specif- ically, let the measured data for the RO + and RO be respectively denoted as d + tot and d tot , respectively, with d + tot = A + k + tot and d tot = A k tot , where A + and A are simple subsampling matrices that extract the measured entries of k + tot and 2 As noted in the introduction, the widely-used assumption that dierent readout polarities are dierent modula- tions of some original image (which forms the basis for the analogy between parallel imaging and ghost correction) is a simplication of the true imaging physics. However, it should be noted that Eqs. (2.3) and (2.4) may still possess low-rank structure even if this assumption is violated. In particular, these LORAKS matrices will have low-rank structure as long as there exist shift-invariant linear prediction relationships in k-space, and the LORAKS approach is agnostic to the original source of such relationships. We believe that deriving the existence of linear prediction relationships in the presence of more realistic models (e.g., accounting for eddy currents) may be feasible, although it is beyond the scope of the present work. Nevertheless, the results we show later with real data seem to imply that linear predictability assumptions are reasonably applicable in the real scenarios we have examined. 17 k tot (i.e., A + and A are formed by concatenating the rows of the identity matrix corresponding to the k-space sampling masks for each polarity). With these denitions, we can now describe previous SLM-based EPI ghost cor- rection methods with a consistent language. For example, the earliest and simplest such approach, ALOHA [71], can be viewed as a special case of the optimization problem n ^ k + tot ; ^ k tot o = argmin fk + tot ;k tot g J C P (k + tot ) C P (k tot ) ; (2.5) subject to the additional data-consistency constraints that A + ^ k + tot = d + tot and A ^ k tot = d tot . Here, J () is a cost function that depends only on the singular values of its matrix argument, and promotes low-rank solutions. In the sequel, we will use the notation L C (k tot ) =J C P (k + tot ) C P (k tot ) ; (2.6) where k tot concatenates k + tot and k tot . Similarly, we will also use L S () to denote the function with the same form as Eq. (2.6), but switching from the LORAKS C-matrix to the LORAKS S-matrix by replacing all instances of C P with S P . A popular choice for J () in the general low-rank matrix completion literature (and the choice made by Ref. [71]) is the nuclear norm, which is a convex function that is known to encourage minimum-rank solutions [105]. The nuclear norm of a matrix G is dened as kGk = rank(G) X i=1 i (G); (2.7) 18 where i (G) is theith singular value of G. Another potential choice ofJ (), which was proposed in the original LORAKS work [35] but which has not been used in the EPI ghost correction work by other groups, is dened by J r (G) = rank(G) X i=r+1 ( i (G)) 2 ; (2.8) where r is a user-selected parameter. This cost function is nonconvex, and J r (G) will equal zero whenever rank(G) r. However, if rank(G) > r, then J r (G) will be nonzero, and equal to the squared Frobenius norm error that is incurred when G is optimally approximated by a rank-r matrix. As a result, this cost function will encourage the reconstructed image to have a LORAKS matrix that is approximately rank-r or lower. The following subsection provides a novel theoretical analysis of the optimization problem from Eq. (2.5), which reveals that it has several undesirable characteristics. 2.3 Analysis of SLM EPI Ghost Correction For our analysis, we assume a typical setup in which the nominal fully-sampled k-space dataset has equally-spaced consecutive phase encoding positions. For the sake of brevity, we will assume fully-sampled single-channel EPI imaging 3 in which d + tot corresponds to the full set of measured even phase encoding positions, while d tot corresponds to the full set of measured odd phase encoding positions. 3 Generalized theoretical results for the case of parallel imaging with uniformly undersampled phase encoding can also be derived using the same principles we used for the single-channel fully-sampled case. We have elected not to show these derivations because they are intellectually straightforward extensions of the single-channel fully-sampled case, but require a lot of additional notation to describe. 19 Notice that the form of A + implies that A H + A + is a diagonal projection matrix, and that multiplying any vector of k-space samples by A H + A + is equivalent to preserving the values of the even phase encoding lines while setting the values of the odd phase encoding lines to zero. Similarly, A H A is a diagonal projection matrix, and multiplying any vector of k-space samples by A H A is equivalent to preserving the values of the odd phase encoding lines while setting the values of the even phase encoding lines to zero. Additionally, we have that A H + A + = IA H A , where I is the identity matrix. Using these facts together with the vector space concepts of orthogonal com- plements and direct sums [85], we know that iff ^ k + tot ; ^ k tot g obeys the data delity constraint from Eq. (2.5), then there exist corresponding vectors y and z such that we can write ^ k + tot = A H + d + tot + A H A y ^ k tot = A H d tot + A H + A + z: (2.9) We have the following theoretical results: Theorem 1. Given the context described above and arbitrary vectors y and z, the singular values of the matrix C P ( ^ k + tot ) C P ( ^ k tot ) (2.10) 20 are identical to the singular values of the matrix C P ( ~ k + tot ) C P ( ~ k tot ) ; (2.11) with ~ k + tot = A H + d + tot A H A y ~ k tot = A H d tot A H + A + z: (2.12) Note that ~ k + tot and ~ k tot from Eq. (2.12) are identical to the vectors ^ k + tot and ^ k tot appearing in Eq. (2.9), except that the estimates of the unmeasured data samples have been multiplied by -1. The proof of this theorem is sketched in the Supplementary Information sec- tion. Some basic intuition for this result is that we can multiply our estimates for the unmeasured k-space lines by -1 without impacting delity with the measured data. Due to uniform subsampling of each gradient polarity by a factor of 2, this multiplication procedure is equivalent to applying linear phase in k-space, which corresponds to a spatial shift of the image by half the FOV along the phase en- coding dimension (and a 180 constant phase oset for the RO polarity). This shifting procedure has no eect on the image support or on the correlations that exist between the dierent coils, and thus has no impact on the singular values or the rank of the LORAKS matrix. This means that if we have one solution to Eq. (2.5), then it is easy for us to construct another solution to Eq. (2.5), and this optimization problem will generally not have a unique useful solution. 21 The following corollaries formalize some of these statements and provide addi- tional useful insight. Corollary 1. Equation (2.5) either has the unique solution f ^ k + tot ; ^ k tot g = fA H + d + tot ; A H d tot g which corresponds to zero-lling of the measured data, or it has at least two distinct optimal solutions that share exactly the same cost function value. Corollary 2. If the cost function J() is chosen to be convex (e.g., the nuclear norm), then the zero-lled solution is always an optimal solution of Eq. (2.5). If Eq. (2.5) has more than one optimal solution, then it has innitely many optimal solutions. Corollary 3. Theorem 1 and Corollaries 1 and 2 are still true if we replace C P () in Eqs. (2.5), (2.10), and (2.11) with S P (). Corollary 1 is proven in the Supplementary Information section, and implies that the optimization problem of Eq. (2.5) either has a trivial undesirable solution corresponding to zero-lling of the measured data, or it is an ill-posed optimization problem that does not possess a unique solution. While some of the solutions to Eq. (2.5) may be desirable, there are no guarantees that the algorithm we use to minimize Eq. (2.5) will yield one of these desirable solutions. Corollary 2 is also proven in the Supplementary Information section, and suggests that the use of convex cost functions to impose LORAKS constraints is likely to be suboptimal relative to the use of nonconvex cost functions. It should be noted that, unlike recent Nyquist ghost correction methods [71, 91, 93] which have made use of the 22 (a) Original Data (b) Modulated Data Figure 2.1: (a) EPI magnitude (top) and phase (bottom) images corresponding to (left) RO + data and (right) RO data. These alias-free images correspond to real data obtained with the PLACE method [123]. (b) Images corresponding to the same data from (a), except that the odd k-space lines for RO + and the even k-space lines for RO have been multiplied by -1. 0 20 40 60 80 100 0 2 4 6 8 ×10 −2 SingularValueIndex SingularValues DatafromFig. 2.1(a) DatafromFig. 2.1(b) Figure 2.2: Plots of the singular values for the LORAKS matrices from Eq. (2.3) for the k-space datasets from Figs. 2.1(a) and (b). convex nuclear norm, the early structured low-rank matrix completion methods for MRI all made use of nonconvex cost functions [32, 35, 38, 110]. These nonconvex options are likely to be better for this problem setting. Corollary 3 is stated without proof (but can be proved using an approach that is similar to our proof of Theorem 1), and indicates that the deciencies of Eq. (2.5) are not alleviated by switching from the LORAKS C-matrix to the LORAKS S-matrix. 23 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.1 1.2 1.3 1.4 1.5 α (a) Convex (Eq. (5) with (7)) −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 ×10 −3 α (b) Nonconvex (Eq. (5) with (8)) −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 1 2 3 ×10 −4 α (c) Noncvx. w/SENSE (Eq. (13)) −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 2 4 6 ×10 −2 α (d)Noncvx.w/AC-LORAKS(Eq.(14)) Figure 2.3: Lettingk 1 andk 2 denote the k-space data for the images in Figs. 1(a) and 1(b), respectively, we plot the cost function value L C k 1 + (1)k 2 as a function of . Setting = 0 yields the cost function value for the images from Fig. 2.1(a), setting = 1 yields the cost function value for the images from Fig. 2.1(b), while setting = 0:5 yields the cost function value for the zero-lled solution. Results are shown for dierent (a) convex and (b-d) nonconvex choices of J(). The rank parameter r = 40 has been used in the nonconvex cases shown in (b-d), both without (a,b) and with (c,d) additional constraints. Practical illustrations of these theoretical results are shown in Figs. 2.1{2.3. Fig- ure 2.1 shows two dierent sets of EPI images that are both perfectly consistent with standard fully-sampled EPI data. The dierence between the two datasets is the same as the dierence between Eqs. (2.9) and (2.12). As expected, this k-space phase dierence leads to shifting of the images for both RO + and RO by half the FOV, as well as adding a constant phase oset for the RO image. Figure 2.2 shows a plot of the singular values of the LORAKS C-matrices corresponding to these two datasets. As expected from Theorem 1, the singular values are identical in both cases. Figure 2.3(a,b) illustrates the dierence in behavior between convex and nonconvex cost functionsJ(). Figure 2.3(a) shows that the zero-lled solution is a minimum of Eq. (2.5) in the convex case, as expected from Corollary 2, and that there are many dierent images with very similar cost function values. Notably, the images from Fig. 2.1(a) and (b) are not optimal solutions in this convex case, even though they both have high quality and appear to be devoid of ghost artifacts. 24 Figure 2.3(b) shows that the cost function has more desirable behavior in the non- convex case (e.g., the zero-lled solution is no longer a minimum of Eq. (2.5), and there are sharp local minima in the vicinity of the images from Fig. 2.1), although the solution to Eq. (2.5) is still not unique in this case as we should expect based on Corollary 1. While it may be possible to get a useful result from solving Eq. (2.5), it should be noted that in the presence of multiple global minimizers, it is dicult to ensure that an optimization algorithm will always converge to a desirable minimum. Incor- porating additional constraints on the solution is a straightforward way to reduce the ambiguity associated with Eq. (2.5), and we describe two practical approaches for this in the next section. 2.4 Constrained SLM Formulations 2.4.1 Formulation using SENSE Constraints A natural approach to imposing additional constraints on SLM reconstruction is to use coil sensitivity map information within the SENSE framework [103], assuming that coil sensitivity proles are available and that data is acquired using a multi- channel receiver array. This style of approach has been used previously for both EPI reconstruction (assuming ghosts have been precorrected using navigator data) [66] and for navigator-free EPI ghost correction [91,93]. Our proposed approach can be viewed as a combination of these two previous formulations. 25 In this work, we propose to use the following formulation for navigator-free EPI ghost correction using SENSE: f^ + ; ^ g = arg min f + ; g kE + + d + tot k 2 2 (2.13) +kE d tot k 2 2 +L S (k tot ); (2.14) subject to the constraints that E + = k + tot , that E = k tot , and that k tot is the concatenation of k + tot and k tot . In this formulation, we are using SENSE to reconstruct one image for RO + ( + ) and another image for RO ( ), and the only coupling that occurs between the two comes from the SLM regularization term. We have also used E + , E , and E to denote the standard SENSE matrices (mapping from the image domain to k-space) corresponding to RO + subsampling, RO subsampling, and full Nyquist sampling, respectively. These three matrices all use exactly the same sensitivity proles, and dier only in the associated k-space sampling patterns. In addition, is a regularization parameter, and we suggest the use of the nonconvex regularizer from Eq. (2.8) for L S (). A major dierence between this proposed approach and our previous SENSE- LORAKS work [66] is the separation of the RO + and RO datasets, which enables navigator-free ghost correction. The main dierence between this proposed ap- proach and MUSSELS [91,93] is that MUSSELS used the LORAKS C-matrix and nuclear norm regularization, while we advocate use of the LORAKS S-matrix with nonconvex regularization. Another major dierence from MUSSELS is that, in 26 the multi-channel case, we use the multi-channel LORAKS matrices (concatenat- ing 2N c SLMs) of Eqs. (2.3) and (2.4) [66], instead of concatenating only 2 SLMs formed from the SENSE reconstructions of each gradient polarity. Equation (2.14) has been written assuming single-shot data. In the multi-shot case with phase inconsistencies between dierent shots, we generalize Eq. (2.14) by reconstructing a separate image for each polarity and each shot, with a separate SENSE encoding matrix and data delity term for each. Note that separating the data from dierent shots increases the eective acceleration factor for each data consistency term and is also associated with additional computational complexity. There are many ways to solve the optimization problem in Eq. (2.14). In this paper, we use the simple majorize-minimize approach from Ref. [66], as described in the Supplementary Information section. As shown in Fig. 2.3(c), incorporating SENSE constraints changes the shape of the cost function and can help to resolve the uniqueness issues associated with the inverse problem. In particular, we observe that in this case, the extra infor- mation provided by SENSE constraints reduces the ambiguity between solutions, and causes the local minimum associated with = 0 to be substantially preferred over the local minimum associated with = 1. However, it appears that in this case, the nonconvex cost function may still have multiple local minima, in which case careful initialization may be necessary to ensure that the iterative approach converges to a good local minimum. 27 2.4.2 Formulation using AC-LORAKS Constraints Another natural approach is to use k-space constraints like those of GRAPPA [28], SPIRiT [86], and PRUNO [132]. In this work, we use a formulation based on AC-LORAKS [36] (with strong similarities to PRUNO [132]). Specically, in the single-shot case, we solve n ^ k + tot ; ^ k tot o = arg min fk + tot ;k tot g C P (k + tot )N 2 F (2.15) + C P (k tot )N 2 F +L S (k tot ); (2.16) subject to the constraints that d + tot = A + k + tot , d tot = A k tot , and that k tot is the concatenation of k + tot and k tot . In this expression,kk F denotes the Frobenius norm, and the matrix N is an estimate of the approximate right nullspace of a LORAKS C-matrix formed from ACS data acquired in a standard parallel imaging calibration pre-scan. The rst two terms of Eq. (2.16) are similar to the rst two terms of Eq. (2.14), in the sense that they impose support and parallel imaging constraints derived from some form of prescan, but do not make any assumptions about the relationship between k + tot and k tot or the relationship between the image- domain phase characteristics of the calibration data and the image-domain phase characteristics of the EPI data being reconstructed. Similarly, the third terms in Eqs. (2.14) and (2.16) are the only terms that introduce coupling between k + tot and k tot , and the only terms that use the LORAKS S-matrix to introduce constraints on the image phase. The use of phase constraints is useful both for partial Fourier 28 EPI acquisition and for stabilizing the reconstruction of symmetrically-acquired EPI data [66]. As before, we suggest using the nonconvex regularizer from Eq. (2.8) for L S (). As shown in Fig. 2.3(d), incorporating AC-LORAKS constraints also changes the shape of the cost function. Just like with SENSE constraints, we observe that the additional information provided by the AC-LORAKS constraints leads to less ambiguity, and a clear preference towards the desired solution at = 0. However, the AC-LORAKS approach appears to be even more benecial than the SENSE approach in this case, since we no longer observe a local minimum associated with = 1 in this case. As a result, we may expect that the AC-LORAKS approach is less sensitive to local minima and yields better solutions than the SENSE approach. Similar to before, Eq. (2.16) is written for the single-shot case, but the general- ization to multi-shot EPI is straightforward by separating and jointly reconstructing images for each polarity and shot. And similar to Eq. (2.14), there are also many ways to solve the optimization problem in Eq. (2.16). These two optimization prob- lems have very similar structure (i.e., the rst two terms are least-squares penalties, while the third term encourages low-rank matrix structure), and as a consequence, we have used a minor modication of the algorithm we used for solving Eq. (2.14) to also solve Eq. (2.16). An interesting feature of our proposed AC-LORAKS formulation is that it can potentially work with single-channel data [36]. This is a major advantage over our proposed SENSE formulation, which is not expected to produce good results unless 29 a multiple-channel receiver array is used. It should be noted that while many human MRI experiments use parallel imaging, single-channel reconstruction is still highly relevant in a variety of situations, including animal studies and studies of human anatomy that make use of specialized receiver coil technology (e.g., single-channel prostate coils). . SENSE LORAKS (unconstrained) LORAKS (SENSE-Convex) MUSSELS LORAKS (SENSE) LORAKS (AC-LORAKS) . . R = 1 . R = 2 . R = 3 . R = 4 . R = 5 Figure 2.4: Comparison of dierent reconstruction techniques using prospectively undersampled (except the R = 5 case, which is retrospectively undersampled) in vivo single-shot EPI data at dierent parallel imaging acceleration factors. 30 2.5 Results This section describes evaluations of our new LORAKS-based EPI ghost correction methods using navigator-free gradient-echo EPI data acquired from phantoms and in vivo human brains. The phantom data shown in Fig. 2.6 were acquired using a gradient echo EPI sequence on a Siemens 3T Tim Trio using a standard product 12- channel receiver array (compressed in hardware down to 4 channels). The imaging parameters were: FOV = 24 mm 24 mm; matrix size 6464; slice thickness = 5 mm; TR = 1.46 sec; TE = 38 msec. The phantom data used in Figs. 2.15 and 2.16(a) were acquired using a gradient echo EPI sequence on a Siemens 3T Tim Trio using a standard product 12-channel receiver array. The imaging parameters were: FOV = 24 mm 24 mm; matrix size 9696; slice thickness = 5 mm; TR = 60 msec; TE = 28 msec. The in vivo data were acquired using a gradient echo EPI sequence on a Siemens (Erlangen, Germany) 3T Prisma Fit using a standard product 32- channel receiver array. The imaging parameters were: FOV = 220 mm 220 mm; matrix size 128128; slice thickness = 3 mm; TR = 2.08 sec; TE = 47 msec; acceleration factor R2f1; 2; 3; 4; 5g. For each dataset, ACS data (used both for estimating nullspaces and for estimating sensitivity maps) was acquired using the same approach as previously used in DPG [42]. Most of the reconstructions we show in this section estimate separate images for each gradient polarity and each shot, and some also estimate separate images for each coil. While various approaches exist for combining together multiple images from dierent coils/polarities/shots for visualization, for simplicity and consistency we have combined the multiple 31 R = 1 R = 2 R = 3 R = 4 (a) DPG R = 1 R = 2 R = 3 R = 4 (b) LORAKS (AC-LORAKS) Figure 2.5: Comparison between (a) DPG and (b) LORAKS with AC-LORAKS constraints for the real single-shot EPI in vivo brain data from Fig. 2.4. Instead of showing coil-combined images, a sin- gle representative channel is shown to avoid contamination of the phase characteristics induced by coil combination. images into a single image using principal component analysis, which is a standard method for parallel imaging coil compression/combination [14,48]. Unless otherwise specied, our LORAKS-based results also always use the nonconvex regularization penalty from Eq. (2.8), with the rank threshold r chosen based on the singular values of the LORAKS matrix formed from ACS data. Specically, r was chosen as the point at which the plot of the singular values appears to atten out, which is a standard approach to matrix rank estimation in the presence of noise. For methods that use regularization parameters, was initially set to a small value ( = 10 3 ), and if necessary based on visual assessment of image ghost artifacts, was gradually increased until good reconstructions were observed. LORAKS-based reconstructions were performed based on adaptations of publicly-available code [31]. Figure 2.4 shows a comparison of dierent parallel imaging reconstruction and EPI ghost correction methods for in vivo single-shot EPI data. A gold standard image with fully-sampled RO + data and fully-sampled RO images was obtained 32 using PLACE [123] with a 32-channel receiver coil and an 128128 acquisition ma- trix. We also acquired standard fully-sampled EPI (acceleration factorR = 1, with each gradient polarity undersampled by a factor of two) and prospectively acceler- ated EPI acquisitions for a range of acceleration factors (R = 2, 3, 4). Additionally, the acceleration factor ofR = 5 was simulated by retrospectively undersampling the PLACE data. Reconstructions were performed using unconstrained LORAKS as in Eq. (2.5) but using the LORAKS S-matrix, LORAKS with SENSE constraints as in Eq. (2.14) in both convex (Eq. (2.7)) and nonconvex (Eq. (2.8)) variations, and LORAKS with AC-LORAKS constraints as in Eq. (2.16). For SENSE-based recon- struction, sensitivity proles were estimated using ESPIRiT [118]. For comparison, we also performed MUSSELS reconstruction [91, 93] and independent SENSE re- construction of each gradient polarity [44] without LORAKS-based regularization (equivalent to setting = 0 in Eq. (2.14)). The gure shows that SENSE without SLM regularization works well for low- acceleration factors, though faces challenges at high acceleration factors. This be- havior is expected as an EPI acceleration factor ofR = 5 is an eective acceleration factor of R = 10 for each readout polarity, which is a very challenging case for SENSE reconstruction. We also observe that unconstrained LORAKS reconstruc- tion has severe problems, as should be expected based on our theoretical analysis of Eq. (2.5). The results also show that the convex SLM approaches (MUSSELS and LORAKS with SENSE and convex regularization) are eective at low acceleration factors, but start demonstrating artifacts as the acceleration factor increases. On 33 Uncorrected ALOHA LORAKS (unconstrained) LORAKS (SENSE) DPG LORAKS (AC-LORAKS) Figure 2.6: (top) Magnitude and (bottom) phase images corresponding to reconstruction of unaccelerated (R = 1) single-channel single-shot EPI phantom data. the other hand, both of our proposed new formulations are substantially more suc- cessful, achieving high quality reconstruction results even at very high acceleration factors. This result is consistent with our theoretical expectations from Corollaries 1-3 that nonconvex cost functions can lead to a better-posed reconstruction prob- lem. At the highest acceleration factors, LORAKS with AC-LORAKS constraints was more eective than LORAKS with SENSE constraints, which displayed unre- solved aliasing artifacts. The data shown in Fig. 2.4 was also reconstructed using the state-of-the-art DPG method [42] using the same ACS data, and a comparison against LORAKS with AC- LORAKS constraints is shown in Fig. 2.5. While DPG generally works well, a close examination of the reconstructed magnitude and phase images demonstrates that DPG still has small residual ghost artifacts that are not present in the LORAKS- based reconstruction. These artifacts are particularly visible in the phase images, since the phase is highly sensitive to ghosting in regions of the image where the magnitude is small. A deeper examination of the data leads us to believe that the 34 Table 2.1: NRMSE values for the images shown in Fig. 2.8. R SENSE LORAKS LORAKS MUSSELS LORAKS LORAKS (unconstr) (SENSE-Conv) (SENSE) (AC) 1 0.17 0.70 0.17 0.17 0.09 0.03 2 0.23 0.86 0.20 0.20 0.10 0.05 3 0.38 1.04 0.27 0.26 0.11 0.06 4 0.53 0.94 0.42 0.49 0.15 0.07 5 0.65 1.00 0.62 0.67 0.25 0.12 ghost artifacts we see for DPG are the result of systematic changes (between the relative phases of the dierent gradient polarities) that have occurred in part due to the length of time that passed between the collection of the ACS data and the acquisition of the accelerated EPI data that is being reconstructed. The data used in Figs. 2.4 and 2.5 was prospectively sampled, but it is hard to quantify accuracy for prospectively sampled data because the image phase charac- teristics (due to eddy currents, etc.) can vary as a function of sequence parameters like the acceleration factor, and because our subject is living and breathing (which leads to variations over time). An illustration of the phase dierences between RO + and RO for dierent acceleration factors is shown in Fig. 2.7. To enable a quantitative comparison of dierent approaches, we also performed reconstruc- tions of retrospectively undersampled versions of the fully sampled PLACE data. Results are shown in supplementary Fig. 2.8, with normalized root-mean-squared error (NRMSE) shown in Table 2.1. The retrospective results are consistent with our prospective results, and the quantitative NRMSE results are consistent with our qualitative evaluations. The data shown in Figs. 2.4 and 2.5 was acquired with an axial slice orientation, and therefore had less phase mismatch between RO + and RO than it could have if 35 . Coil 6 Coil 12 Coil 18 Coil 24 Coil 27 Coil 30 Gold Standard R = 2 (AC-LORAKS) R = 3 (AC-LORAKS) R = 4 (AC-LORAKS) Figure 2.7: Phase dierences between RO + and RO for a representative subset of coils from the prospec- tively undersampled data shown in Fig. 4. In cases with accelerated data (where no gold standard is available), we show the phase dierence obtained after AC-LORAKS reconstruction. Unlike most other phase images shown in this paper (which do not include background masking and show the entire phase range from to ), we have taken steps to make the phase images in this gure easier to visualize. In particular, we have masked the background noise to make it easier to focus on the signal structure. In addition, we have shown a restricted phase range for both the gold standard (black = 0.95 radians, white = 1.02 radians) and the reconstructions using AC-LORAKS (black = 0.89 radians, white = 1.21 radians). we had used a less-conventional slice orientation. In addition, the slice we used was far from sources of eld inhomogeneity, and therefore had a relatively smooth phase prole. To demonstrate a more complicated case, we have also performed ghost correction of EPI data acquired with a double-oblique slice orientation. Double- oblique orientation is non-traditional for EPI imaging of the brain, but is known to give rise to 2D nonlinear phase-dierences between gradient polarities due, e.g., to concomitant eld eects [22, 27, 42, 106, 125]. In addition, the double-oblique slice we selected passes close to air-tissue interfaces and has a substantially less-smooth phase prole. We have performed reconstructions of prospectively undersampled 36 Gold Standard SENSE LORAKS (un- constrained) LORAKS (SENSE- Convex) MUSSELS LORAKS (SENSE) LORAKS (AC-LORAKS) . R = 1 . R = 2 . R = 3 . R = 4 . R = 5 Figure 2.8: Comparison of dierent reconstruction techniques using retrospectively undersampled in vivo single-shot EPI data to simulate dierent parallel imaging acceleration factors. double-oblique data (similar to Figs. 2.4 and 2.5), and the results are shown in Figs. 2.10 - 2.12. We have also performed reconstructions of retrospectively under- sampled double-oblique data (similar to Fig. 2.8 and Table 2.1), and the results are shown in Fig. 2.13 and Table 2.2. As can be seen from the phase dierence maps in Fig. 2.11, we have more signicant nonlinear 2D phase mismatches between RO + and RO in this case. Our results with double-oblique data are consistent with our 37 Uncorrected ALOHA LORAKS (unconstrained) LORAKS (SENSE) DPG LORAKS (AC-LORAKS) Figure 2.9: (top) Magnitude and (bottom) phase images corresponding to reconstruction of unaccelerated (R = 1) single-channel single-shot EPI in vivo human brain data. Table 2.2: NRMSE values for the images shown in Fig. 2.13. R SENSE LORAKS LORAKS MUSSELS LORAKS LORAKS (unconstr) (SENSE-Conv) (SENSE) (AC) 1 0.06 0.70 0.09 0.09 0.07 0.03 2 0.07 0.86 0.13 0.15 0.08 0.07 3 0.19 0.91 0.20 0.25 0.11 0.10 4 0.46 0.99 0.35 0.40 0.18 0.12 5 0.61 0.95 0.65 0.74 0.29 0.11 results from axial data, suggesting that our proposed approaches can still work well in more challenging scenarios. The results described above were all based on data acquired with a somewhat loose FOV, which may be benecial for the use of support constraints. However, there are also imaging scenarios of interest in which the FOV is much tighter. To test performance in the presence of a tight FOV, we acquired similar data to that shown in Figs. 2.4 and 2.5, but with half the FOV along the phase encoding dimension. This causes our gold standard image to demonstrate aliasing. To make the reconstruction problem even more challenging, we also simulated an additional 2D nonlinear phase dierence between RO + and RO , and this phase dierence was chosen dierently for the ACS data and the gold standard used for simulation. 38 . SENSE LORAKS (un- constrained) LORAKS (SENSE- Convex) MUSSELS LORAKS (SENSE) LORAKS (AC-LORAKS) . . R = 1 . R = 2 . R = 3 . R = 4 . R = 5 Figure 2.10: Comparison of dierent reconstruction techniques using prospectively undersampled in vivo double-oblique single-shot EPI data at dierent parallel imaging acceleration factors. 39 . Coil 1 Coil 9 Coil 18 Coil 22 Coil 26 Coil 30 Gold Standard R = 2 (AC-LORAKS) R = 3 (AC-LORAKS) R = 4 (AC-LORAKS) R = 5 (AC-LORAKS) Figure 2.11: Phase dierences between RO + and RO for a representative subset of coils from the prospectively undersampled double-oblique data shown in Fig. 2.10. In cases with accelerated data (where no gold standard is available), we show the phase dierence obtained after reconstruction using LORAKS with AC-LORAKS constraints. Unlike most other phase images shown in this paper (which do not include background masking and show the entire phase range from to ), we have taken steps to make the phase images in this gure easier to visualize. In particular, we have masked the background noise to make it easier to focus on the signal structure. In addition, we have shown a restricted phase range for both the gold standard (black = 0.95 radians, white = 1:02 radians) and the reconstructions using AC-LORAKS (black = 0.86 radians, white = 1.18 radians). Since it is dicult to apply conventional SENSE-based methods in the presence of a tight FOV (because of the presence of aliasing in the gold standard), we only performed reconstructions using methods that do not use sensitivity maps, i.e., DPG and LORAKS with AC-LORAKS constraints. The results are reported in Fig. 2.14. As can be seen, DPG does not perform very well in this very complicated scenario, particularly because of the phase mismatch between the ACS data and 40 R = 1 R = 2 R = 3 R = 4 R = 5 (a) DPG R = 1 R = 2 R = 3 R = 4 R = 5 (b) LORAKS (AC-LORAKS) Figure 2.12: Comparison between (a) DPG and (b) LORAKS with AC-LORAKS constraints for the real single-shot EPIin vivo brain double-oblique data from Fig. 2.10. Instead of showing coil-combined images, a single representative channel is shown to avoid contamination of the phase characteristics induced by coil combination. the data being reconstructed. On the other hand, AC-LORAKS is substantially more successful. In addition to navigator-free multi-channel settings, the proposed methods were also evaluated in navigator-free single-channel settings, which are expected to be substantially more challenging. Single-channel datasets were obtained by isolating the information from a single coil in multi-channel acquisitions. Reconstructions were performed using ALOHA [71] (Eq. (2.5) with Eq. (2.7)), unconstrained LO- RAKS (Eq. (2.5) with Eq. (2.8) but using the LORAKS S-matrix), LORAKS with SENSE constraints (Eq. (2.14)), DPG [42], and LORAKS with AC-LORAKS con- straints (Eq. (2.16)). Note that, since sensitivity-map estimation is not feasible in 41 Gold Standard SENSE LORAKS (un- constrained) LORAKS (SENSE- Convex) MUSSELS LORAKS (SENSE) LORAKS (AC-LORAKS) . R = 1 . R = 2 . R = 3 . R = 4 . R = 5 Figure 2.13: Comparison of dierent reconstruction techniques using retrospectively undersampled in vivo double-oblique data to simulate single-shot EPI experiments at dierent parallel imaging acceleration factors. the single-channel setting, our SENSE-based results used a binary support mask (that has value 1 inside the support of the image and value 0 everywhere else) in place of a coil sensitivity map. This has the eect of imposing prior knowledge of the image support on the reconstructed image. Note also that DPG was not originally designed to be used with single-channel data, although the formulation can still be applied to the single-channel case. Unaccelerated (R = 1) single-channel single- shot EPI results are shown for a phantom dataset (6464 acquisition matrix) in 42 Gold Standard Uncorrected DPG LORAKS (AC-LORAKS) Phase Dierence (ACS) Phase Dierence (Gold Standard) Figure 2.14: (top-left) Magnitude and (bottom-left) phase images corresponding to reconstruction of unaccelerated (R = 1) in vivo axial multi-channel single-shot EPI data acquired with a small FOV, and with extra simulated nonlinear phase dierences between RO + and RO . The images on the right show the phase dierences that were used between RO + and RO , which were set dierently for the ACS data (used for calibration) than they were for the gold standard (used for reconstruction). Fig. 2.6 and for one channel of the previous in vivo human brain dataset in Fig. 2.9. The results are consistent in both cases. Images obtained without compensating the mismatch between RO + and RO have obvious ghost artifacts, and these arti- facts are not solved (and are potentially even amplied) when using unconstrained SLM approaches (i.e., ALOHA or LORAKS without constraints). The LORAKS reconstruction with \SENSE" constraints (i.e., support constraints) helps to elim- inate some of the ghost artifacts that appeared outside the support of the original object, although residual aliasing artifacts are still observed within the support of the object. These artifacts are most visible in the phantom image, though close inspection also reveals the appearance of aliasing artifacts in the brain image. We also observe that DPG is unsuccessful in this single-channel case, which we believe is due both to the diculties of the single-channel problem as well as systematic dierences between the ACS data and the data being reconstructed. On the other hand, the LORAKS results with AC-LORAKS constraints are substantially more successful than any of the previous methods. 43 Image 1 Image 2 Image 3 Image 4 (a) Uncorrected (b) DPG (c) LORAKS (AC-LORAKS) 10 12 14 16 18 0.1 0.2 0.3 0.4 0.5 0.6 Time(sec) Relativerespiratoryposition(unitless) ACS1stshot ACS2ndshot Image1 Image2 Image3 Image4 1 (d) Simulated Respiration Figure 2.15: Evaluation with multi-shot EPI data for a phantom with simulated respiratory eects. A representative set of four segmented images is shown, extracted from a longer acquisition spanning several minutes. (a) Phase images corresponding to one channel of the data, with images reconstructed without compensating for the mismatches between dierent shots and dierent gradient polarities. (b) Phase images for one channel of the DPG reconstruction, with reconstruction performed using multi-channel data. (c) Phase images for single-channel LORAKS reconstruction with AC-LORAKS constraints. (d) Plot showing the relative respiratory position across EPI shots (TR=60msec), as measured with an ultrasound transducer coupled to a respiratory phantom [41]. The line-plot peaks show points where the phantom air bag is maximally in ated. The sampling times for the ACS data and for the shots used to generate each of the two-shot images from (a)-(c) are marked as labeled in the legend. We also applied the dierent reconstruction approaches to single-channel double- oblique data, for the same reasons as in the multi-channel case. Results are shown in Fig. 2.17, and are consistent with the results observed with axial data. One potential criticism of the single-channel results we have shown is that, at least for these datasets, the phase mismatch between RO + and RO is not very severe, which causes the uncorrected reconstructions to have relatively low levels of ghost artifact. To demonstrate the performance in a more severe case, we performed reconstruction of this dataset again with an exaggerated phase mismatch between RO + and RO . The results of this simulation are shown in Fig. 2.18. As can be seen, this case demonstrates much stronger ghost artifacts without correction. 44 Zero-filled Uncorrected LORAKS (single-channel) LORAKS (multi-channel) (a) Phantom dataset with simulated respiration (b) In vivo human brain dataset Figure 2.16: Images reconstructed using LORAKS with AC-LORAKS constraints for accelerated (R = 2) data in both single-channel and multi-channel contexts. (a) Phantom dataset with simulated respiration. (b) In vivo human brain dataset. Phase images from a single channel are shown for zero-lled data recon- structed without compensating for mismatches between the gradient polarities, LORAKS reconstruction with AC-LORAKS constraints from single-channel data, and LORAKS reconstruction with AC-LORAKS constraints from multi-channel data. However, similar to the previous cases, LORAKS with AC-LORAKS constraints is still the most successful at reducing ghost artifacts. Reconstructions were also performed using multi-shot data. Figure 2.15 shows results using fully sampled R = 1 two-shot data (128128 acquisition matrix, 12-channel receiver coil) at dierent time points for a phantom with physically- simulated respiration eects [41]. Note that with R = 1 and two-shots, each gradi- ent polarity for each shot has an eective undersampling factor of 4. LORAKS with AC-LORAKS constraints is compared against DPG for segmented EPI [40], and to make reconstruction even more challenging for LORAKS, LORAKS reconstruction was performed from single-channel data while DPG was provided with the full set of multi-channel data. Due to simulated respiration, there are mismatches between 45 Uncorrected ALOHA LORAKS (un- constrained) LORAKS (SENSE) DPG LORAKS (AC-LORAKS) Figure 2.17: (top) Magnitude and (bottom) phase images corresponding to reconstruction of unaccelerated (R = 1) in vivo double-oblique single-channel single-shot EPI data. Uncorrected ALOHA LORAKS (un- constrained) LORAKS (SENSE) DPG LORAKS (AC-LORAKS) Phase Dierence Figure 2.18: (top, middle) Magnitude and (bottom) phase images corresponding to reconstruction of unaccelerated (R = 1) in vivo double-oblique single-channel single-shot EPI with strong simulated phase dierences between RO + and RO . The middle row shows a dierent windowing of the magnitude images to more clearly visualize ghost-artifacts. The top-right image shows the simulated phase dierence between RO + and RO . 46 the measured ACS data and the EPI data being reconstructed. DPG is not ro- bust to these mismatches, and displays residual ghost artifacts with time-varying characteristics. On the other hand, LORAKS with AC-LORAKS demonstrates robustness against the time-varying changes in this dataset, even despite the chal- lenging single-channel multi-shot nature of this reconstruction problem. For further insight, Fig. 2.16(a) shows an even more challenging case where LORAKS with AC-LORAKS constraints is used to reconstruct a single-shot of the previous multi-shot dataset. Note that this corresponds to EPI with an acceleration factor of R = 2 (i.e., an eective acceleration factor of R = 4 for each gradient po- larity). Remarkably, we observe that single-channel LORAKS with AC-LORAKS constraints is still successful in this very dicult scenario, with similar quality to that obtained using multi-channel LORAKS with AC-LORAKS constraints. How- ever, it is also important for us to point out that single-channel reconstruction with R = 2 is not always successful, as illustrated in Fig. 2.16(b) using previously de- scribed in vivo human brain data. In this brain case, we observe that single-channel LORAKS reconstruction is unsuccessful at correctly reconstructing the image, while the multi-channel case (also shown in Fig. 2.5) yields accurate results as described previously. We suspect that the dierence between the phantom result and the in vivo result is explained by dierences in the size of the FOV relative to the size of the object. Specically, the size of the FOV is more than twice the size of the phantom, while the same is not true for the in vivo data. We expect that 47 LORAKS-based single-channel image reconstruction will be harder for tight FOVs than it is for larger FOVs that contain a larger amount of empty space. 2.6 Discussion and Conclusions This paper derived novel theoretical results for EPI ghost correction based on struc- tured low-rank matrix completion approaches. Key theoretical results include the observation that the corresponding matrix completion problem is ill-posed in the absence of additional constraints, and that convex formulations have undesirable characteristics that are somewhat mitigated by the use of nonconvex formulations. These theoretical results led to two novel problem formulations that use additional constraints and nonconvex regularization to avoid the problems associated with ill- posedness. Our results showed that these new approaches are both eective relative to state-of-the-art ghost correction techniques like DPG, and are capable of han- dling nonlinear 2D phase mismatches between RO + and RO . It was also observed that the proposed variation that uses AC-LORAKS constraints appears to be more eective than the proposed variation that uses SENSE constraints in challenging scenarios with single-channel data or highly-accelerated multi-channel data. In addition, we were surprised to observe that the variation using AC-LORAKS con- straints can even be successful when applied to undersampled single-channel data. We believe that these approaches will prove valuable across a range of dierent applications, but especially those in which navigator-based ghost correction meth- ods are ineective, in which the need for navigation places undesirable constraints 48 on the minimum achievable echo time or repetition time, in cases where single- channel imaging is unavoidable, or in cases where imaging conditions are likely to change over time over the duration of a long experiment (e.g., functional imaging or diusion imaging). This paper focused on simple proof-of-principle demonstrations of the proposed approaches, and there are opportunities for a variety of improvements. For ex- ample, the choices of r and were made manually using heuristic approaches in this paper, though it is likely that these same decision processes could be made automatically using standard techniques for rank estimation (e.g., [17]) or ghost correction parameter tuning (e.g., [111]). Additionally, extensions to the case of simultaneous multi-slice imaging are possible within the LORAKS framework [61], and are likely to be practically useful in a range of experiments given the advantages and the modern popularity of combining this approach with EPI [101]. While we derived our new theory and methods in the context of EPI ghost correction, we believe that this paper also has more general consequences. For ex- ample, we believe that it is straightforward to generalize our theoretical results to show that unconstrained LORAKS-based reconstruction will be ill-posed for any application that uses uniformly-undersampled Cartesian k-space trajectories (also see similar comments in Ref. [66]). We also believe that the combination of LO- RAKS with additional constraints will always be benecial when the constraints are accurate, and so encourage the use of additional constraints whenever the LO- RAKS reconstruction problem is ill-posed. In addition, we believe that our novel 49 LORAKS formulation with AC-LORAKS constraints is an important innovation that is likely to be useful in other applications, similar to how LORAKS with SENSE constraints has already proven to be useful in other settings [63, 66]. For example, it may be benecial to use AC-LORAKS constraints to augment existing SLM approaches that have recently been proposed for gradient delay estimation in non-Cartesian MRI [55,92]. Finally, while recent work has proposed the use of con- vex LORAKS-based formulations, our empirical experience since we rst started exploring LORAKS several years ago [34,35] has consistently been that nonconvex formulations are substantially more powerful than convex ones. 50 2.7 Supplementary Information 2.7.1 Proof of Theorem 1 Assuming the notation of Eq. (2.5), let B denote the zero-lled LORAKS matrix B = C P (A H + d + tot ) C P (A H d tot ) (2.17) and let D denote the matrix corresponding to the unmeasured data samples D = C P (A H A y) C P (A H + A + z) : (2.18) Due to the way the LORAKS C-matrix is constructed, if the entry in the mth column and nth row of the matrix B is nonzero, then the corresponding entry of the matrix D is required to be zero and vice versa. Note also that the matrix from Eq. (2.10) can be written as B+D, while the matrix from Eq. (2.11) can be written as B D. To avoid additional tedious notation, we will assume in our proof sketch that the rows and columns of the matrix B have been permuted in such a way that samples from even and odd lines in k-space are never adjacent to one another in the matrix, which is always possible based on the convolutional structure of the LORAKS C-matrix. This allows the B matrix to be written in a \checkerboard" form 51 B = 2 6 6 6 6 6 6 6 6 6 4 b + 11 0 b + 13 0 b + 15 0 b 12 0 b 14 0 0 b + 22 0 b + 24 0 b 21 0 b 23 0 b 25 b + 31 0 b + 33 0 b + 35 0 b 32 0 b 34 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 7 7 7 7 7 7 7 7 7 5 ; (2.19) where b + ij and b ij are the nonzero entries of the B matrix corresponding to pos- itive and negative readout polarities, respectively. Using the same permutation scheme, the matrix D can similarly be written in the corresponding complemen- tary \checkerboard" form: D = 2 6 6 6 6 6 6 6 6 6 4 0 d + 12 0 d + 14 0 d 11 0 d 13 0 d 15 d + 21 0 d + 23 0 d + 25 0 d 22 0 d 24 0 0 d + 32 0 d + 34 0 d 31 0 d 33 0 d 35 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 7 7 7 7 7 7 7 7 7 5 ; (2.20) 52 Consider the diagonal matrix Q 1 which has the same number of columns as B, and whose diagonal entries alternate in sign in a way that follows the non-zero pattern of the rst row of B: Q 1 = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 0 0 0 0 0 0 0 0 0 01 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 01 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 01 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 01 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 01 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 : (2.21) 53 Similarly, consider the diagonal matrix Q 2 which has the same number of rows as B, and whose diagonal entries alternate in sign in a way that follows the non-zero pattern of the rst column of B: Q 2 = 2 6 6 6 6 6 6 6 6 6 4 1 0 0 0 1 0 0 0 1 . . . . . . . . . . . . 3 7 7 7 7 7 7 7 7 7 5 : (2.22) We make the following observations: • The matrices Q 1 and Q 2 are unitary and satisfy Q 1 1 = Q 1 and Q 1 2 = Q 2 . • The matrix B is structured in such a way that BQ 1 = Q 2 B. • The matrix D is structured in such a way that DQ 1 =Q 2 D. • The matrix Q 2 (B + D)Q 1 simplies according to Q 2 (B + D)Q 1 = Q 2 Q 2 (B D) = B D: (2.23) From Eq. (2.23), we can infer that if we write the singular value decomposition of B + D as B + D = UV H , then the matrix B D can be written as ~ U ~ V H , where ~ U = Q 2 U and ~ V = Q 1 V. Since ~ U and ~ V are matrices with orthonormal columns, we must have that ~ U ~ V H is a valid singular value decomposition of 54 B D. Thus, we can conclude that B D and B + D have identical singular values. This completes the proof of the theorem. 2.7.2 Proof of Corollary 1 Assume that the vectors y and z are chosen such that the expressions in Eq. (2.9) represent an optimal solution to Eq. (2.5). Theorem 1 then tells us that the vectors from Eq. (2.12) represent another optimal solution to Eq. (2.5). These two solutions are identical to one another if and only if A H A y = 0 and A H + A + z = 0, in which case both solutions are equal to the zero-lled solution. In this case, the zero- lled solution is clearly an optimal solution to Eq. (2.5), and must be the unique optimal solution if Eq. (2.5) only has a single solution. If either A H A y6= 0 or A H + A + z6= 0, then Eqs. (2.9) and (2.12) represent two distinct solutions to Eq. (2.5), indicating that Eq. (2.5) does not have a unique solution. 2.7.3 Proof of Corollary 2 Based on the proof of Corollary 1, we know that if Eq. (2.5) has a solution that is not equal to the zero-lled measured data, then we can use Eqs. (2.9) and (2.12) to obtain a pair of two distinct solutions to Eq. (2.5). Let k 1 and k 2 denote these two solutions, and notice that by the denition of an optimal solution, we must have that L C (k 1 ) = L C (k 2 ) and that L C (k 1 ) L C (k) for all possible candidate solutions k. 55 Corollary 2 is easily proven based on the denition of a convex function. Specif- ically, if L C (y) is convex, then it must satisfy [85] L C (y 1 + (1)y 2 )L C (y 1 ) + (1)L C (y 2 ); (2.24) for every possible pair of vectors y 1 and y 2 and for every real-valued scalar between 0 and 1. Setting y 1 = k 1 and y 2 = k 2 in Eq. (2.24) leads to L C k 1 + (1)k 2 L C k 1 +(1)L C k 2 =L C k 1 ; (2.25) Combining Eq. (2.25) with the previous observation that L C (k 1 ) L C (k) for all possible candidate solutions k implies that L C k 1 + (1)k 2 =L C k 1 . As a result, k 1 + (1)k 2 must also be an optimal solution of Eq. (2.5) for every possible choice of 0 1, and we have successfully proven that there exist an innite number of solutions. Specically, an innite set of solutions is given byfk k 1 + (1)k 2 j2 [0; 1]g. Additionally, the zero-lled solution is obtained as one of these solutions, corresponding to the specic choice of = 0:5. 2.7.4 Majorize-Minimize Algorithm For the sake of completeness, this section provides an operational description of the majorize-minimize algorithm we use to solve Eq. (2.14). The key observation [31,35] 56 that enables this algorithm is that the function J r (G) from Eq. (2.8) is majorized at a point ^ G (j) by the function g G; ^ G =kG ^ G (j) k 2 F : (2.26) Using this majorant, Eq. (13) can be optimized by iteratively solving easy-to- optimize surrogate problems. Specically, using ^ (j) to denote the estimate of f + ; g at the jth iteration, and given some initial guess ^ (0) , we minimize Eq. (2.14) by iteratively solving the surrogate problems ^ (j+1) = arg min f + ; g kE + + d + tot k 2 2 +kE d tot k 2 2 +k S P (E + + ) S P (E ) G (j) r k 2 F (2.27) forj=0; 1; 2;::: until convergence, where G (j) r is the optimal rank-r approximation (obtained by truncating the singular value decomposition) of the matrix formed by numerically evaluating S P (E + + ) S P (E ) at the pointf + ; g = ^ (j) . These surrogate problems have the form of simple linear least squares problems, and can be solved using standard iterative algorithms like the conjugate gradient algorithm. See Refs. [31,35] for a more detailed description and justication of this algorithmic approach. 57 Chapter 3 Robust Autocalibrated Structured Low-Rank EPI Ghost Correction 3.1 Introduction Echo-planar imaging (EPI) is a widely-used high-speed MRI acquisition strategy [113], but is subject to several undesirable artifacts [5]. Nyquist ghosts are one of the most common EPI artifacts, and occur because of systematic dierences between the interleaved lines of k-space that are acquired with dierent readout gradient polarities, and/or because of systematic dierences between interleaved lines of k-space data that are acquired with dierent shots in a multi-shot acquisition. Despite substantial eorts over several decades to solve this problem [5, 12, 15, 19, 21,42,44,46,51,71,77,84,89,93,111,123{125,127], the widely-deployed modern ghost correction schemes are still prone to incomplete ghost suppression, as illustrated in Fig. 3.1. 58 Uncorrected Navigator RAC-LORAKS Figure 3.1: Illustration of EPI ghost correction. The top row of this gure shows EPI images obtained from dierent methods, while the bottom row shows the same images with 10 intensity amplication to highlight the ghost characteristics. If EPI data is naively reconstructed without accounting for the systematic dierences between data acquired with positive and negative readout gradient polarities (\Un- corrected"), then strong Nyquist ghosts appear in the image as indicated with arrows. Modern EPI techniques frequently try to eliminate these artifacts using navigator information to estimate the sys- tematic dierences between the data collected with dierent readout polarities. In the navigator-based example we show (\Navigator"), the navigator information was collected using a 3-line EPI acquisition with the phase encoding gradients turned o, and the dierence between positive and negative gradient polarities was modeled using constant and 1D linear phase terms. Although this approach substantially reduces Nyquist ghosts, it is common for some amount of residual ghosting to still be present in the images, particularly in cases where simple 1D phase modeling is inadequate to capture the dierences between the two gradient polarities. We also show an example of our proposed approach (\RAC-LORAKS"), which can account for more complicated variations between the dierent gradient polarities, and which is substantially more successful at suppressing Nyquist ghosts in this example. Recently, structured low-rank matrix methods for ghost correction [71,77,84,89, 93] have received increasing attention for their ability to provide excellent ghost- suppression performance without the need for additional \navigator" information (i.e., reference scans collected alongside each EPI readout that allow estimation of the systematic inconsistencies between dierent gradient polarities or dierent shots). These methods can suppress ghosts better than navigator-based methods, and eliminate the need to acquire navigators for each EPI readout. Among dif- ferent structured low-rank matrix approaches, a ghost correction method based on Autocalibrated LORAKS (AC-LORAKS) [36] was previously demonstrated to yield high-quality results across a range of dierent scenarios [84]. To eliminate 59 a fundamental ambiguity in structured low-rank matrix recovery from uniformly undersampled EPI data [84], AC-LORAKS makes use of parallel imaging subspace information estimated from autocalibration (ACS) data acquired in a pre-scan. This ACS-based approach is similar to standard autocalibrated parallel imaging methods like GRAPPA [28], SPIRiT [86], and PRUNO [132]. While the AC-LORAKS approach to ghost correction generally works well when the ACS data is pristine and well-matched to the EPI data to be reconstructed, there are many situations where experimental conditions (e.g., subject motion, eddy currents, etc.) can lead to artifacts within the ACS data or mismatches between the ACS and EPI data. The performance of the AC-LORAKS ghost correction procedure degrades in the presence of these ACS artifacts and mismatches. Note that this kind of issue is not unique to AC-LORAKS or to ghost correction, and imperfect ACS/calibration data is a longstanding and commonly-reported problem for all calibration-based image reconstruction methods [4,42,45,109,114,130]. For AC-LORAKS ghost correction, imperfect ACS data can be especially troublesome in contexts where the prescan would be done once before acquiring a long sequence of multiple EPI images (e.g., in BOLD fMRI or diusion MRI applications), and then used to reconstruct each image in the sequence. In this paper, we propose an extension of AC-LORAKS for EPI ghost correction that is more robust to imperfections in the ACS data. The new method, called Ro- bust Autocalibrated LORAKS (RAC-LORAKS), has two major dierences from the previous AC-LORAKS approach. First, RAC-LORAKS does not completely 60 trust the subspace information learned from the ACS data, but rather uses a novel structured low-rank matrix formulation that learns subspace information jointly from both the (imperfect) ACS data and the EPI data being reconstructed. To the best of our knowledge, no previous methods have used this kind of approach to address the longstanding issue of imperfect ACS data. And second, RAC-LORAKS uses the ACS data to provide additional complementary information for the recon- struction of the EPI data within a multi-contrast joint reconstruction framework [9]. Preliminary accounts of the rst strategy were previously reported in recent con- ferences [79, 83], although we have not previously reported the combination with the second strategy. 3.2 Theory Due to space constraints, our descriptions in this paper will assume that the reader is already familiar with the basic physics of EPI. Readers interested in a more detailed explanation are referred to classic references [5, 113]. For simplicity, our description of EPI ghost correction will generally assume that we are correcting ghosts associated with the dierences between data acquired with dierent readout gradient polarities in a single-shot EPI experiment. However, since the ghost model for bipolar gradients is nearly identical to the ghost model for multi-shot acquisition, the same approach is easily adapted mutatis mutandis to multi-shot acquisition with an arbitrary number of shots [84]. 61 3.2.1 Background: Structured Low-Rank EPI ghost correction Structured low-rank matrix methods for EPI ghost correction [71,77,84,89,93] can be viewed as an extension of structured low-rank matrix methods for conventional MR image reconstruction [35, 38, 56, 75, 97, 110, 118, 132], and are based on the same underlying theoretical principles. In particular, it has been shown that when MRI images have limited support, smooth phase variations, multi-channel corre- lations, or transform-domain sparsity, then the MRI k-space data will be linearly predictable [37], which means that convolutional Hankel- or Toeplitz-structured matrices formed from the k-space data will possess low-rank characteristics. This observation means that MRI reconstruction can be reformulated as structured low- rank matrix recovery. Importantly, these structured low-rank matrix recovery meth- ods can even be successful in calibrationless scenarios where ACS data or other prior information about the spatial support, phase, or multi-channel sensitivity proles is not available [35,38,110]. Structured low-rank EPI ghost correction methods combine these principles with the earlier observation that EPI data acquired from dierent gradient polarities or dierent shots can be treated as coming from dierent eective \channels" in a parallel imaging experiment, where the systematic dierences between dierent po- larities or shots lead to dierent phase or magnitude modulations of the underlying EPI image [19, 42, 44, 124]. Since structured low-rank methods for conventional image reconstruction automatically account for the unknown sensitivity maps that modulate the underlying image in a parallel imaging experiment, it is reasonable 62 to apply these same types of methods to handle the unknown polarity- or shot- dependent modulations that manifest in EPI ghost correction. For the sake of brevity, we will focus the remainder of our review on the AC- LORAKS method for EPI ghost correction [84], since our proposed RAC-LORAKS method is a generalization of AC-LORAKS. The AC-LORAKS method for EPI ghost correction is based on solving the following regularized optimization problem subject to exact data consistency constraints: n ^ k + ; ^ k o = arg min fk + ;k g P C (k + )N 2 F + P C (k )N 2 F +J r ([P S (k + )P S (k )]): (3.1) In this expression, k + and k respectively represent the ideal fully-sampled multi- channel Cartesian k-space data for the positive and negative readout gradient polar- ities;P C () is the LORAKS operator that maps the k-space data into a structured matrix that should possess low-rank if the multi-channel image possess limited support and/or interchannel parallel imaging correlations;P S () is the LORAKS operator that maps the k-space data into a structured matrix that should possess low-rank if the multi-channel image possess limited support, smooth phase, and/or interchannel parallel imaging correlations; the matrix N comprises an orthonormal (i.e., N H N = I) basis for the nullspace of the matrix 2 6 6 4 P C (k acs +) P C (k acs ) 3 7 7 5 , where k + acs and k acs respectively represent the ACS data for the positive and negative readout gra- dient polarities; is a regularization parameter; J r () is a regularization penalty that promotes low-rank characteristics; andkk F denotes the Frobenius norm. Due 63 to space constraints, this paper will not provide a detailed recipe for implementing the LORAKS operatorsP C () andP S (), and simply note that our implementa- tions for this paper are identical to those that are described in detail in earlier LORAKS papers [35, 38]. There are theoretical benets to choosing a nonconvex low-rank regularization penalty [84], and the previous AC-LORAKS approach for ghost correction [84] used the nonconvex function proposed in the original LORAKS paper [35] dened by J r (X) = min Y kX Yk 2 F s.t. rank(Y)r; (3.2) where r is a user-selected rank parameter, X is a matrix representing the point at which we are evaluating the function J r (X), and Y is an optimization variable of the same size as X. This penalty encourages matrices that have accurate rank-r approximations. The rst two terms appearing on the right hand side of Eq. (3.1) respectively im- pose limited support and parallel imaging constraints on the reconstructions of the positive and negative readout gradient polarities. The constraints that are used in these terms are implicit in the low-rank characteristics of the structured LORAKS matrices, as captured by the nullspace matrix N. The nullspace matrix is learned in advance from the ACS data, and as a result, there is an implicit assumption that the support and parallel imaging constraints that were valid for the ACS data are also valid for the EPI data to be reconstructed. Note that if the third term were 64 removed from Eq. (3.1), then these rst two terms would reduce to performing sep- arate PRUNO [132] or conventional AC-LORAKS [36] reconstructions of the data from each polarity. Acquiring ACS/calibration data is relatively fast and easy, and is already a standard part of most modern parallel imaging protocols, so is not very burdensome on the acquisition. Using ACS data can also be important in this context, since it has been mathematically proven that structured low-rank ma- trix methods for ghost correction suer from fundamental ambiguities unless some form of side information is available [84]. While other options exist for removing ambiguity (e.g., using SENSE-like [103] image-domain constraints [84, 93]), it was previously observed that the AC-LORAKS approach (i.e., using GRAPPA-like [28] Fourier-domain constraints) oered better performance [84]. The third term of Eq. (3.1) couples the reconstruction of the two polarities to- gether, allowing the reconstruction of one polarity to benet from information from the other polarity, while also introducing phase constraints to allow the reconstruc- tion to benet from k-space conjugate symmetry characteristics. In particular, the third term implicitly and automatically imposes the following constraints when- ever they are compatible with the measured data: limited image support, smooth phase, interchannel parallel imaging correlations, and interpolarity correlations. Notably, these constraints are all imposed implicitly through the nullspace of a structured matrix, and if a given constraint is not compatible with the measured data, then that constraint will automatically not be imposed by the reconstruction procedure [37]. 65 The ACS data for AC-LORAKS ghost correction has typically been acquired using the same process used by Dual Polarity GRAPPA (DPG) [42, 43, 76, 100]. In particular, assuming a parallel imaging acceleration factor of R, DPG employs a 2R-shot EPI prescan. The data from dierent shots and dierent gradient po- larities is then rearranged and interleaved to form one fully-sampled ACS dataset comprised only of data acquired with a positive readout gradient polarity (k + acs ) and another fully-sampled ACS dataset comprised only of data acquired with a negative readout gradient polarity (k acs ). Since this ACS acquisition strategy is based on a multi-shot approach, it therefore may be prone to ghosting artifacts due to shot-to- shot variations. In addition, since the ACS data is often acquired only once at the beginning of a long multi-image EPI scan (e.g., in BOLD fMRI or diusion MRI experiments), the ACS data acquired at the beginning of the experiment may grad- ually become mismatched with data acquired at later time points due to scanner drift, subject motion, etc. As noted previously, the ghost correction performance of AC-LORAKS can be substantially degraded when there are mismatches between the ACS data and the EPI data to be reconstructed. Although a pre-processing procedure has been previously developed to correct for shot-to-shot variations in the ACS data for DPG [42], this approach is not sucient for the present context. In particular, this approach undesirably modies the magnitude and phase charac- teristics of the ACS data in ways that are not well-suited for AC-LORAKS, and only addresses ACS artifacts without accounting for mismatches that may exist between the ACS data and the EPI data. 66 3.2.2 RAC-LORAKS Our proposed RAC-LORAKS method is based on solving the following optimization problem n ^ k + ; ^ k ; ^ N o = arg min fk + ;k ;Ng P C (k + )N 2 F + P C (k )N 2 F (3.3) + P C (k + acs )N 2 F + P C (k acs )N 2 F +J r P S (k + ) P S (k ) P S (k + acs ) P S (k acs ) subject to exact data consistency constraints on k + and k and subject to or- thonormality constraints on N such that N H N = I. This optimization problem involves four user-selected parameters: the regularization parameters and , the rank parameterr, and the number of columnsp of the matrix N (which determines the dimension of the approximate nullspace). Equation (3.3) has two main dierences from Eq. (3.1). First, instead of choosing a predetermined value of the approximate nullspace matrix N that depends only on the ACS data, N is now an optimization variable that depends on both the ACS data and the EPI data to be reconstructed. This allows the reconstruction to be more robust against possible imperfections in the ACS data. The extent to which the ACS data is trusted is controlled by the user-selected parameter. In the limit as !1, the approximate nullspace matrix N will converge to the xed matrix from Eq. (3.1). 67 The second dierence is that the nal term of Eq. (3.3) now includes structured matrices formed from the ACS data, in addition to the previous structured matrices formed from the EPI data to be reconstructed. By concatenating the ACS data in this way, we are essentially treating the ACS data in the same way that we would treat additional channels in a parallel imaging experiment. Although the ACS data may not have the same contrast as the EPI data to be reconstructed, it has previously been shown that treating multi-contrast information like additional channels in a parallel imaging experiment often leads to improved reconstruction performance [9]. While this improvement has been justied empirically, some level of theoretical justication for this approach can be obtained by modeling dierent image contrasts as dierent modulations of some latent image [37]. Algorithmically, Eq. (3.3) can be minimized using existing algorithms for LO- RAKS optimization [35, 38, 62, 84]. In particular, it is not hard to show that the solution to Eq. (3.3) can be equivalently obtained by solving: n ^ k + ; ^ k o = arg min fk + ;k g J (Cp) 0 B B B B B B B B B @ 2 6 6 6 6 6 6 6 6 6 4 P C (k + ) P C (k ) p P C (k + acs ) p P C (k acs ) 3 7 7 7 7 7 7 7 7 7 5 1 C C C C C C C C C A +J r P S (k + ) P S (k ) P S (k + acs ) P S (k acs ) ; (3.4) where J (Cp) () is the same as J r () but replacing the rank parameter r with the rank parameter (Cp), whereC is the number of columns of the LORAKS matrix 68 formed byP C (). Equation (3.4) is convenient because it takes the same form as previous LORAKS optimization problems involving multiple J r () terms [35]. For this paper, we use a multiplicative half-quadratic majorize-minimize algorithm to minimize this objective function [62], which takes advantage of FFT-based matrix multiplications to improve computational complexity [98]. The RAC-LORAKS solution is obtained through the optimization of a non- convex cost function. As such, the algorithm has the potential to converge to an undesirable local minimum. For the results shown in this paper, we initialize RAC- LORAKS using a naive initialization with minimal processing cost as explained in the next section. Other choices could potentially result in even higher performance, but are not considered here. 3.3 Methods 3.3.1 Datasets used for Evaluation As described below, we evaluated the characteristics of RAC-LORAKS using data from several dierent contexts. All in vivo data were acquired under IRB-approved written informed consent. 3.3.1.1 Gradient-Echo EPI Brain data In one set of experiments, we acquired in vivo human brain data using a gradient- echo EPI sequence with parameters that are somewhat similar to a BOLD fMRI 69 experiment. Data was acquired on a Siemens 3T Prisma Fit scanner using a stan- dard 32-channel receiver array. The data was acquired using FOV = 220 mm 220 mm; matrix size = 128 128; slice thickness = 3 mm; and TR = 2 sec. In one subject, data was acquired without acceleration (R = 1) with TE = 47 msec. From this same subject, data was also acquired for parallel imaging acceleration factors of R = 2; 3; 4 with TE = 35 msec. In a second subject, data was acquired for parallel imaging acceleration factors of R = 5; 6 with TE = 35 msec. In all cases, fully-sampled ACS data was acquired using the same interleaved 2R-shot EPI prescan as used for DPG [42]. The previous datasets were acquired with a conventional axial slice orientation. However, because Nyquist ghost problems tend to be more extreme with oblique acquisitions due to concomitant elds that can produce substantial nonlinear 2D phase dierences between positive and negative readout polarities [22, 27, 42, 106, 125], we also acquired an additional dataset with a double-oblique slice orientation from a third subject to test performance in a more challenging scenario. The slice orientation in this case is nonstandard and likely dicult to interpret for many readers, so we have depicted its position in Fig. 3.2. For this case the data was acquired with TR = 2.08 sec and TE = 35 msec for parallel imaging acceleration factors of R = 1; 2; 3; 4; 5; 6. 70 Sagittal Coronal Axial Figure 3.2: Illustration of the orientation of the double-oblique gradient-echo EPI dataset. The double- oblique slices are shown in red, overlaid on a structural T1-weighted image of the same subject. The double-oblique slice used for the results in Fig. 3.7 is shown with a yellow rectangle. 3.3.1.2 Diusion-encoded EPI Brain Data In another set of experiments, we acquired in vivo human brain data using a diusion-encoded spin-echo EPI sequence. Diusion EPI data might be consid- ered more challenging than the previous gradient-echo EPI data, due to the fact that diusion MRI data usually suers from random image-to-image phase varia- tions, and can also have lower SNR than gradient-echo EPI. In addition, the rapid switching of strong diusion gradients can introduce substantial additional eddy current eects that can cause systematic dierences between the ACS data and the diusion EPI data if they are acquired with dierent diusion gradient settings [70]. A rst diusion dataset was acquired on a Siemens 3T Prisma Fit scanner using a standard 32-channel receiver array. For the sake of computational complexity, this data was subsequently reduced to 16 channels using standard coil-compression techniques. The data was acquired using FOV = 220 mm 220 mm; matrix size = 220 220; slice thickness = 5 mm; TR = 2.8 sec; TE = 63 msec; b-values of 0 sec/mm 2 and 1000 sec/mm 2 ; 6 diusion encoding directions; parallel imaging 71 acceleration factor R = 3; and 6/8ths partial Fourier sampling. ACS data was acquired using the same interleaved 2R-shot EPI prescan as used for DPG [42], except that the data was acquired with lower resolution along the phase encoding dimension (i.e., we only acquired 45 phase-encoding lines for the ACS data). Due to the random phase variations associated with diusion encoding gradients, the ACS data was acquired without diusion weighting, which means that the ACS data has very dierent contrast characteristics from the EPI data. To show results across a broader range of acceleration factors, a second set of acquisitions was performed withR = 2; 3; 4; 5. Other parameters were identical to the previous case, except for matrix size = 110 110; slice thickness = 2 mm; TR = 11.4 sec; TE = 73 msec; and fully-sampled ACS data. 3.3.1.3 Cardiac EPI Data In a third set of experiments, we acquired in vivo human cardiac data during diastole using a spin-echo EPI sequence with parameters that are typical for a myocardial arterial spin labeling experiment [54]. Data was acquired on a GE 3T Signa HDx scanner with an 8-channel cardiac coil. The acquisition used FOV = 280 mm 140 mm; matrix size = 128 64; slice thickness = 10 mm; TR = 55 msec; TE = 32.9 msec; velocity cuto = 5 cm/s; no parallel imaging acceleration (R = 1); and 5/8ths partial Fourier sampling. ACS data was acquired using the same interleaved 2R- shot EPI prescan as used for DPG [42], but with 5/8ths partial Fourier sampling. 72 Data was acquired with a double-oblique slice orientation to achieve a mid-short axis view. 3.3.2 Simulations In addition to in vivo data, the dierent methods were also evaluated using simu- lations where a gold standard was present. To form a gold standard with realistic EPI characteristics, we took two in vivo gradient-echo EPI brain datasets (as de- scribed in Section 3.3.1.1) with axial slice orientation andR = 1 from the same scan session, and reconstructed them both using SENSE. Each gradient polarity was re- constructed separately, providing a realistic representation of typical interpolarity image dierences. This procedure provides two sets of fully-sampled multi-channel dual-polarity gold standard images. One of these sets was used for ACS data, while the other was undersampled (including parallel imaging acceleration, along with interleaving the data from positive and negative gradient polarities) to sim- ulate EPI data. These datasets were acquired roughly 5 minutes apart, allowing time for mismatches to evolve. Since ghost correction is frequently more dicult for EPI datasets with 2D nonlinear phase dierences between the two polarities, we applied an additional 2D nonlinear phase pattern to make the problem more challenging. This additional phase dierence was designed to be roughly 3-larger than we observed in the real data from Fig. 3.7 . In a rst set of simulations, to mimic the situation where a localized image feature is dierent between the ACS data and EPI data (e.g., as may happen in a dynamic 73 experiment), we added a Gaussian-shaped additive image hyperintensity to the EPI data that we did not add to the ACS data. The hyperintensity was designed to follow both the coil sensitivity maps (obtained by applying ESPIRiT [118]) and the phase characteristics of the original data. We also performed simulations with these two datasets interchanged, i.e., with the hyperintensity in the ACS data but not in the EPI data. In another set of simulations, to mimic the situation where the ACS data and EPI data have very dierent contrasts, we inverted the magnitude image for the ACS data to create an image with dierent contrast [7], while still following the coil sensitivity maps and the phase characteristics of the original data. In addition to performing multi-channel simulations, we also performed a sim- ulation in a very challenging single-channel setting. For this, single-channel data was obtained by a linear combination of the multi-channel data [13]. Single-channel ghost correction is a dicult setting where only a few previous methods have had any success [71,84]. This case is hard because even with unaccelerated data (R = 1), each polarity has an eective acceleration factor of R = 2 when the data for each readout gradient polarity is separated, and it can be dicult to reconstruct R = 2 data without multi-channel information. The fully-sampled ACS and EPI datasets used for all three simulations are il- lustrated in Fig. 3.3 . 74 . EPI data ACS data EPI Phase diff. ” (EPI Phase diff.) - (ACS Phase diff.) ” Multi-Channel Multi-Channel (Inverted Contrast) ” ” Single Channel ” ” Figure 3.3: Illustration of the EPI and ACS datasets used in simulation. The rst and second top rows show coil-combined multi-channel data for the case when the EPI and ACS data have similar and inverted contrast, respectively, while the bottom row shows representative single-channel images. We also show the interpolarity phase dierence for the coil-combined EPI data, as well as the dierence in the interpolarity phase dierence between the coil-combined EPI and ACS data. 3.3.3 Data Processing RAC-LORAKS was applied to perform reconstruction and ghost correction on these datasets. For comparison against existing methods, the datasets were also recon- structed using the previous AC-LORAKS ghost-correction method [84], DPG [42], and MUSSELS [93]. For some of the datasets we consider, the ACS data may be incomplete due to low-resolution ACS acquisition (i.e., the rst brain EPI diusion data) or partial Fourier ACS acquisition (i.e., the cardiac EPI data). In such cases, we modify RAC-LORAKS to consider the fully sampled ACS data vectors k + acs and k acs as additional variables to be optimized in Eq. (3.4), subject to ACS data consistency constraints. 75 For RAC-LORAKS and AC-LORAKS, the regularization parameters and were selected manually based on subjective visual inspection reconstruction quality and ghost-reduction performance for in vivo data, and to minimize quantitative er- ror measures for simulated data. The rank-related parametersp andr were selected based on the singular value characteristics of LORAKS matrices formed from the ACS data. The rank parameters were set based on the points at which the singular value curves begin to atten out, which is a common rank estimation technique for noisy matrices. This decision was made manually (based on visual inspection) for the results shown later in the paper, although fully automatic approaches would also be viable. DPG is a ghost correction method that treats dierent gradient polarities like dierent coils in a parallel imaging experiment, and uses a dual GRAPPA kernel estimated from ACS data for image reconstruction [42]. In order to use DPG for the initialization of RAC-LORAKS, we have adapted DPG to output two sets of images (with calibration based on the raw uncorrected multi-channel ACS data), one for the original k + data and one for the original k data. This is dierent than the original DPG implementation, which applies ACS pre-processing to try and correct for errors in the ACS data, and then directly fuses information from the two polarities together into a single virtual \hybrid" output [42]. This hybrid output can have dierent magnitude and phase characteristics than the original k + and k data, so is not useful as an initialization for RAC-LORAKS. We refer to our adapted version as modied DPG (mDPG) from now on. In some cases, we also compare 76 against the original version of DPG (including the original ACS pre-processing procedure to correct for shot-to-shot variations in the ACS data [42]), although note that such comparisons are necessarily qualitative, since the magnitude and phase characteristics of the hybrid output DPG images do not match the images generated using other methods. MUSSELS is a structured low-rank matrix recovery method that uses SENSE- type parallel imaging constraints together with nuclear norm regularization to im- pose low-rank constraints [93]. While MUSSELS was originally developed for multi- shot EPI ghost correction, it can apply equally well to the ghost correction problem associated with dierent gradient polarities. Sensitivity maps for MUSSELS were estimated by applying ESPIRiT [118] to the same ACS data used for the other methods. The regularization parameter for MUSSELS was selected manually based on subjective visual inspection of reconstruction quality and ghost-reduction per- formance in the case of in vivo data, or to minimize quantitative error measures for simulated data. Note that DPG and MUSSELS were both developed for the multi-channel set- ting. We can adapt DPG to the single-channel setting in straightforward ways [84], and we apply this adaptation to the single-channel simulated data. We did not adapt MUSSELS to the single-channel case. Note that the SENSE-based con- straints used by MUSSELS would reduce to a simple spatial-domain support con- straint in the single-channel case, which is not strong enough to yield good perfor- mance results. 77 For all methods, results were visualized by using a standard square-root sum-of- squares technique to combine the images from dierent coils and dierent gradient polarities into a single image. Results from in vivo experiments were evaluated qualitatively, since a gold standard reference was not available in these cases. Simulation results were evaluated quantitatively using the normalized root mean- squared error (NRMSE): NRMSE, r ^ k + k + gold 2 2 + ^ k k gold 2 2 q k + gold 2 2 + k gold 2 2 ; (3.5) where k + gold and k gold are respectively the gold standard values for the positive and negative gradient polarities. We also plotted Fourier Error Spectrum Plots (ESPs) to gain further insight into how the errors were distributed across dierent spatial resolutions scales [65]. An ESP is designed to reveal the spectral characteristics of the error, and for example, can discriminate between methods that make more errors in the low-resolution features of an image versus methods that make more errors in high-resolution features. 3.4 Results Figure 3.4 shows ACS data and reconstruction results from the in vivo gradient- echo EPI brain data with an axial slice orientation. The ACS data in this case does not have strong artifacts, although close inspection does reveal that ACS ghost artifacts are present. This can be further appreciated in Fig. 3.5 where the same 78 images are shown with amplied image intensity to highlight ghost characteristics in the image background. As can be seen, all ghost correction methods work well at smaller acceleration factors, although performance begins to degrade at larger acceleration factors. We observe that, compared to other methods, the visual qual- ity of the MUSSELS reconstruction seems to degrade most rapidly as a function of acceleration factor, which is consistent with previous observations [84]. The mDPG method had qualitatively better performance than MUSSELS in this case. How- ever, a close inspection of the images reveals that the mDPG results are not entirely ghost-free even for the unaccelerated (R = 1) case. This may be expected due to the artifacts and mismatches that are present in the ACS data. Although mDPG does not attempt to correct the ACS artifacts, it should be noted that the original DPG method does try to correct them through pre-processing. Results showing the qualitative performance of the original DPG method are shown in Fig. 3.6, where we observe that the ghost artifacts still exist, though as expected, are less prominent than were observed for mDPG. In spite of the ACS artifacts, the AC- LORAKS reconstruction still has good performance at low acceleration factors and does a good job of suppressing ghosts in the background regions of the image at all acceleration factors, although exhibits substantial degradation in image qual- ity at the highest acceleration factors (with artifacts similar to those observed for highly-accelerated parallel imaging reconstructions). However, the RAC-LORAKS reconstruction appears to have much higher quality than the other methods, even at very high acceleration factors like R = 6. (Note that when R = 6, the eective 79 acceleration factor is R = 12 when each readout gradient polarity is considered separately. This leads to a highly ill-posed inverse problem). R = 1 ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 R = 2 R = 3 R = 4 R = 5 R = 6 Figure 3.4: ACS data and reconstruction results for in vivo gradient-echo EPI brain data with an axial slice orientation for dierent parallel imaging acceleration factors. Note that the rst four acceleration factors (R = 1-4) were acquired from one subject during a single scan session while the last two acceleration factors (R = 5; 6) were acquired from a dierent subject on a dierent day, which explains the visual discontinuity between these cases. Figure 3.7 shows results from the in vivo gradient-echo EPI brain data with a double-oblique slice orientation. This case is more challenging than the previous one due to the complicated nonlinear 2D spatial phase dierences we observed between data acquired with positive and negative polarities (as visualized in the last 80 R = 1 ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 R = 2 R = 3 R = 4 R = 5 R = 6 Figure 3.5: The same results shown in Fig. 3.4 , but with a 5 intensity amplication to highlight the ghost characteristics. column of Fig. 3.7), the proximity to air-tissue interfaces that result in substantial magnetic eld inhomogeneity eects, as well as more substantial ghosting artifacts present in the ACS data. Note that the ACS data corresponding to the R = 5 case is particularly corrupted, which can be attributed to the unpredictable shot- to-shot variations that frequently occur in these kinds of multi-shot acquisitions. Despite the more extreme scenario, the dierent ghost reconstruction methods have 81 hola R = 1 R = 2 R = 3 R = 4 R = 5 R = 6 mDPG mDPG (5) DPG DPG (5) Figure 3.6: DPG results corresponding to the same data shown in Fig. 3.4 and Fig. 3.5. The same mDPG results shown in Fig. 3.4 and Fig. 3.5 are also reproduced in this gure for reference. Note that the processing steps of DPG cause the image intensities to be mismatched from the intensities of mDPG and the other reconstruction methods, which precludes a quantitative comparison. similar characteristics to those observed in the previous case, with RAC-LORAKS appearing to demonstrate the cleanest overall results. Figure 3.8 shows reconstruction results from the rst set of multi-channel simu- lations (with similar contrast between ACS and EPI data, but with a hyperintensity added to the EPI data). Quantitative NRMSE values are reported in Table 3.1 with corresponding ESPs shown in Fig. 3.10. Qualitatively, the results from Fig. 3.8 have similar characteristics to the results observed with in vivo data. Notably, RAC-LORAKS is able to consistently reconstruct a high-quality image that bears close resemblance to the gold standard image, while methods like mDPG and AC- LORAKS have artifacts due to the small mismatches between the ACS and EPI 82 R = 1 ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS Phase Diff. ” R = 1 R = 2 ” R = 3 ” R = 4 ” R = 5 ” R = 6 ” Figure 3.7: ACS data and reconstruction results for in vivo gradient-echo EPI brain data with a double- oblique slice orientation for dierent parallel imaging acceleration factors. For reference, we also show the interpolarity phase dierence as estimated from a coil-combined RAC-LORAKS result. The degree of phase nonlinearity is an indicator of how dicult ghost correction is expected to be. As can be seen, complicated 2D nonlinear phase dierences are present in many of these cases. data. The visual assessment of reconstruction quality matches well with the quan- titative NRMSE assessment shown in Table 3.1. AC-LORAKS and RAC-LORAKS have a similar performance at R = 1 and 2, with RAC-LORAKS having the best performance at high acceleration factors. Reconstructions were also performed using the original DPG formulation as shown in Fig. 3.9. In this case, DPG has similar ghost artifacts to mDPG, which 83 Table 3.1: NRMSEs for the multi-channel simulation results shown in Fig. 3.8. For each acceleration factor, the smallest values are highlighted in bold. MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 0.059 0.024 0.016 0.020 R = 2 0.104 0.045 0.035 0.042 R = 3 0.271 0.083 0.056 0.055 R = 4 0.572 0.127 0.132 0.064 R = 5 0.741 0.161 0.269 0.085 is expected because there are no artifacts in the ACS data, while there is a prob- lematic mismatch between the ACS data and the EPI data that neither DPG nor mDPG address. Notably, for both DPG and mDPG, we observe aliasing artifacts that seem to be associated with the hyperintensity that was present in the EPI data but was not in the ACS data. The ESP plots in Fig. 3.10 enable a more nuanced analysis. These results suggest that RAC-LORAKS has good (i.e., among the best, even if it is not always the best) performance at all spatial frequencies, meaning that it is good at reconstructing image features across the whole range of resolution scales. Figure 3.11 shows a similar simulation result to that shown in Fig. 3.8, with the main dierence being that the previous EPI images (with the hyperintensity) were used as ACS data and the previous ACS images (without the hyperintensity) were used to generate EPI data. Consistent with the previous case, we observe good per- formance for RAC-LORAKS, and do not observe the features of the hyperintensity being erroneously transferred into the RAC-LORAKS reconstruction results. Figure. 3.12 and Table 3.2 show simulation results for the case where the ACS data has an even more substantial contrast dierence (i.e. inverted contrast) with respect to the EPI data. For this case we observe a degradation in performance for 84 hola MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 R = 2 R = 3 R = 4 R = 5 Figure 3.8: Reconstruction results for the rst set of multi-channel simulations (with similar contrast between ACS and EPI data, but with a hyperintensity added to the EPI data) with dierent parallel imaging acceleration factors. all methods compared to the previous cases, although RAC-LORAKS still showed the best overall qualitative and quantitative performance. This result suggests that RAC-LORAKS may have better performance when the contrast is similar between the ACS and EPI data, although can still provide benets when the contrast dif- ference is substantial. Table 3.2: NRMSEs for the multi-channel inverted contrast simulation results shown in Fig. 3.12 . For each acceleration factor, the smallest values are highlighted in bold. MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 0.098 0.177 0.027 0.025 R = 2 0.131 0.182 0.065 0.045 R = 3 0.280 0.190 0.141 0.077 R = 4 0.532 0.213 0.236 0.117 R = 5 0.754 0.309 0.373 0.160 85 hola R = 1 R = 2 R = 3 R = 4 R = 5 mDPG mDPG (5) DPG DPG (5) Figure 3.9: mDPG and DPG results corresponding to the same multi-channel simulated data from Fig. 3.8. Figure 3.13 shows reconstruction results from the single-channel simulation, with quantitative NRMSE values reported in Table 3.3. While previous work [84] reported that mDPG and AC-LORAKS can be reasonably successful for single- channel data withR = 1 when the ACS data is pristine, our new results demonstrate that this performance can be sensitive to the quality of the ACS data. In partic- ular, we observe strong ghost artifacts for both of these methods, even though we do observe that the AC-LORAKS reconstruction has successfully suppressed ghost artifacts in the image background (outside of the support of the true image). In contrast, RAC-LORAKS is substantially more successful forR = 1. Notably, RAC- LORAKS also performed well for the even more challengingR = 2 case, unlike the 86 R = 1 R = 2 0 20 40 60 80 0 0.5 1 1.5 SpatialFrequency(NyquistUnits) Relativeerror RAC-LORAKS AC-LORAKS MUSSELS mDPG 1 0 20 40 60 80 0 0.5 1 1.5 SpatialFrequency(NyquistUnits) Relativeerror RAC-LORAKS AC-LORAKS MUSSELS mDPG 1 R = 3 R = 4 0 20 40 60 80 0 0.5 1 1.5 SpatialFrequency(NyquistUnits) Relativeerror RAC-LORAKS AC-LORAKS MUSSELS mDPG 1 0 20 40 60 80 0 0.5 1 1.5 SpatialFrequency(NyquistUnits) Relativeerror RAC-LORAKS AC-LORAKS MUSSELS mDPG 1 R = 5 0 20 40 60 80 0 0.5 1 1.5 SpatialFrequency(NyquistUnits) Relativeerror RAC-LORAKS AC-LORAKS MUSSELS mDPG 1 Figure 3.10: ESPs for the multi-channel simulation results shown in Fig. 3.8. The vertical axis of each ESP uses a consistent range to enable comparisons between dierent acceleration factors. other methods. For reference, note that even with high-quality ACS data, the pre- vious AC-LORAKS method did not yield good results with similar single-channel R = 2 data [84]. Figure 3.14 shows reconstruction results from the rst set of in vivo diusion EPI brain data, including a 10 intensity amplication to highlight the ghost char- acteristics. As can be seen, the ACS data has ghost artifacts in all cases, and both MUSSELS and mDPG reconstructions also exhibit unsuppressed ghosting artifacts. 87 hola MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 R = 2 R = 3 R = 4 R = 5 Figure 3.11: Reconstruction results for multi-channel simulated data with dierent parallel imaging ac- celeration factors. These simulations are identical to those reported in Fig. 3.8, except that the images used to generate EPI data and the images used to generate ACS data were interchanged. On the other hand, both AC-LORAKS and RAC-LORAKS are relatively ghost- free in this example and have only minor dierences from one another (it might be argued that the RAC-LORAKS result has a slightly less-noisy appearance than the AC-LORAKS result, but if so, this dierence is very subtle). While this result does not demonstrate an obvious advantage for RAC-LORAKS over AC-LORAKS, it should be observed that this diusion result is at least consistent with the previous gradient-echo EPI data results, in which we also did not observe a substantial dier- ence between RAC-LORAKS and AC-LORAKS whenR = 3. In addition, this case involves a very substantial contrast dierence between the ACS data and the EPI 88 hola MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 1 R = 2 R = 3 R = 4 R = 5 Figure 3.12: Reconstruction results for the second set of multi-channel simulations (with inverted contrast between ACS and EPI data) with dierent parallel imaging acceleration factors. Table 3.3: NRMSEs for the single-channel simulation results shown in Fig. 3.13. For each acceleration factor, the smallest values are highlighted in bold. mDPG AC-LORAKS RAC-LORAKS R = 1 0.170 0.584 0.046 R = 2 0.494 0.809 0.073 data. This dierence does not appear to have adversely aected the performance characteristics of these methods in substantial ways. Figure 3.15 shows reconstruction results from the second set of in vivo diusion EPI brain acquisitions (with dierent acceleration factors), with zoom-ins shown in Fig. 3.16 for improved visibility. Consistent with the results shown for the gradient-echo EPI data in Fig. 3.4, we observe that all methods perform well for low acceleration factors. As the acceleration factor increases, the performance of each 89 R = 1 mDPG AC-LORAKS RAC-LORAKS R = 1 R = 2 Figure 3.13: Reconstruction results for single-channel simulated data with dierent acceleration factors. method degrades, with RAC-LORAKS showing a lower qualitative degradation in comparison to the other methods at the very high acceleration factors R = 4; 5. Note that at high acceleration factors (e.g., R = 4; 5) the reconstruction quality for RAC-LORAKS is not quite as good as for the gradient-echo EPI dataset shown in Fig. 3.4. We believe that this should be expected, since as mentioned before, diusion EPI data can be considered more challenging than the gradient-echo EPI data due to SNR issues, eddy current eects, motion-induced phase eects, and contrast mismatches between the ACS and EPI datasets. Finally, Fig. 3.17 shows results from the in vivo cardiac EPI data. While this data was not accelerated (R = 1), this case is challenging because of the double- oblique slice orientation as well as the substantial artifacts present in the ACS data resulting from cardiac motion-induced shot-to-shot variations. In addition, this case can also be challenging for SENSE-based methods (like MUSSELS), due to the use of a small FOV with aliasing. When aliasing is present within the FOV, it violates the standard SENSE modeling assumption of one sensitivity map value per spatial location, which generally leads to artifacts if not properly accounted 90 R = 1 ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS Slice 1 Slice 1 (10) Slice 2 Slice 2 (10) Slice 3 Slice 3 (10) Figure 3.14: ACS data and reconstruction results for three representative slices from in vivo diusion brain data (R = 3). A 10 intensity amplication is also shown for each slice to better highlight the ghosting characteristics. for. The results demonstrate that both MUSSELS and mDPG have substantial residual ghosting artifacts, which might not be surprising given the high degree of corruption that is present in the ACS data. On the other hand, both AC-LORAKS and RAC-LORAKS are more successful at suppressing the ghosts. Without a gold standard reference, it is hard to establish denitively whether AC-LORAKS or RAC-LORAKS is better in this example, although we believe that the RAC- LORAKS result demonstrates slightly less ghosting than AC-LORAKS, particularly 91 hola ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 2 R = 3 R = 4 R = 5 Figure 3.15: ACS data and reconstruction results for in vivo diusion EPI brain data for dierent parallel imaging acceleration factors. For improved visualization, zoomed-in versions of these results (correspond- ing to the spatial region marked with a yellow rectangle in the rst column and rst row) are shown in Fig. 3.16. It should be noted that the subject appears to have slightly moved between scans, so that there is not perfect correspondence between anatomical image features across dierent acceleration factors. on the left side of the image where the ACS data and mDPG both have particularly strong ghost artifacts. 3.5 Discussion The results in the previous section demonstrated that, in the presence of imperfect ACS data, RAC-LORAKS frequently oers similar or better performance to the previous AC-LORAKS ghost correction method that it generalizes, while both of these methods perform substantially better than methods like MUSSELS or DPG. We also observed that RAC-LORAKS appears to have the biggest advantage over AC-LORAKS in scenarios where the parallel imaging acceleration factor was high. 92 hola ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS R = 2 R = 3 R = 4 R = 5 Figure 3.16: The same results shown in Fig. 3.15, but zoomed-in to a region of interest for improved visualization. For these cases, we observed that RAC-LORAKS was able to mitigate ghost arti- facts both inside and outside the support of the original image, while AC-LORAKS was only able to mitigate ghost artifacts outside the support but not inside. This advantage for RAC-LORAKS is likely the result of its improved robustness to ACS errors combined with the multi-contrast linear predictability constraints which help to make the reconstruction problem less ill-posed. However, it should be noted that RAC-LORAKS has one more regularization parameter than AC-LORAKS (i.e., , which controls the level of trust placed in the information from the ACS data). In our experience, manual tuning of this parameter is not hard (i.e., we always started from the small value = 10 3 , and frequently did not have to modify this value to 93 ACS data MUSSELS mDPG AC-LORAKS RAC-LORAKS Figure 3.17: ACS data and reconstruction results for unaccelerated in vivo cardiac EPI data. The two rows show the same results, but the second row has 5 intensity amplication to better highlight the ghosting characteristics. achieve satisfying results). The method would be easier to use if the selection of were automated. Both RAC-LORAKS and AC-LORAKS also depend on the choice of rank pa- rameters, and as described previously, the results shown in this work made a heuris- tic choice based on the empirical rank characteristics of the ACS data. Even though the low-rank characteristics of the structured matrices might vary between the ACS data and the acquired EPI data due to systematic phenomena (e.g., thermal noise, subject motion, respiration, artifacts in the ACS data, etc.), we have not observed major problems associated with inappropriate rank selection in our empirical re- sults. This might be expected, based on the observation that LORAKS recon- struction results are frequently not very sensitive to small variations in the rank 94 parameter [35, 38]. Nevertheless, the development of improved automatic RAC- LORAKS parameter selection methods would be an interesting topic for future work. Although RAC-LORAKS oers good performance, it should be noted that our current implementation of RAC-LORAKS can be more computationally expensive than existing methods. For example, for the results shown with R = 1 in Fig. 1, RAC-LORAKS used 45 min of reconstruction time, while MUSSELS, mDPG, and AC-LORAKS respectively used 15 min, 2 min, and 100 min. All meth- ods were implemented in MATLAB on a standard desktop computer with an Intel Xeon E5-1603 2.8 GHz quad core CPU processor and 32GB of RAM. While this relatively long computation time may be a concern, it should be noted that we are reporting the results of a simple proof-of-principle implementation, and we did not spend much time to optimize the computational eciency of this approach. We believe that major improvements may be possible by leveraging better computa- tional hardware, smarter algorithms, and more ecient implementations. Given the reconstruction performance oered by RAC-LORAKS, we believe that improving its computational performance is a promising topic for future research. However, RAC-LORAKS is notably faster than AC-LORAKS, and it appears that this speed dierence results from the fact that RAC-LORAKS has consistently faster conver- gence than AC-LORAKS in this setting. The reason for this faster convergence is unclear at this stage, although we believe that a detailed analysis of convergence characteristics is beyond the scope of the present paper. 95 While this paper focused on EPI ghost correction for standard single-slice exci- tation, we believe that the extension of these ideas to simultaneous multi-slice EPI acquisitions (similar to Refs. [8,77,89,90]) is a very promising research direction. Finally, although the techniques we developed in this work were described and evaluated in the context of EPI ghost correction, we believe that the overall ap- proach is likely to be useful across a wide range of parallel imaging applications, particularly those for which the measured ACS data is not adequate to resolve all of the reconstruction artifacts. Specically, we believe that the key principles em- ployed by RAC-LORAKS (i.e., using structured low-rank matrix methods to avoid placing complete trust in the accuracy of ACS data, and leveraging ACS data to provide additional information in a multi-contrast framework) are both novel ideas that are applicable to arbitrary image reconstructions involving ACS data, and are not exclusive to ghost correction settings. In addition, we are encouraged by the high-quality reconstruction results that RAC-LORAKS produces even in very highly-accelerated scans. These results suggest to us that there may be value in exploring the usefulness of RAC-LORAKS to other parallel imaging experiments in future work. 3.6 Conclusions This work proposed and evaluated RAC-LORAKS, a new structured low-rank ma- trix method for EPI ghost correction that integrates multiple constraints (including parallel imaging constraints, support constraints, phase constraints, and inter-image 96 linear predictability constraints) to not only mitigate artifacts resulting from im- perfect ACS data and Nyquist ghosts, but also accounting for partial Fourier ac- quisition and reducing parallel imaging artifacts and noise in an integrated fashion. RAC-LORAKS uses ACS data and k-space domain linear predictive modeling to stabilize the solution of the ill-posed inverse problem, and was observed to oer advantages relative to state-of-the-art ghost correction methods like AC-LORAKS, DPG, and MUSSELS. 97 Chapter 4 On the Shape of Convolution Kernels in MRI Reconstruction: Rectangles versus Ellipsoids 4.1 Introduction Constrained MRI reconstruction methods that rely on shift-invariant convolution models have existed for decades. One of the earliest approaches (which still remains quite popular) is to assume that a missing sample of k-space data can be predicted as a linear shift-invariant combination of neighboring samples [37,74]. In particular, it is common to assume that ^ [k] = X z2 w[z][k z] (4.1) for an appropriately-chosen shift-invariant interpolation kernel w[k], where ^ [k] is the missing sample to be interpolated, the interpolation kernel w[k] has a shape dictated by the support set , and the neighboring samples correspond to [k z] for z2 . In many cases, the interpolation kernel w[z] is a quantity that must 98 be learned from some kind of training data. Methods that rely on these kinds of concepts include popular approaches like SMASH [112], GRAPPA [28], SPIRiT [86], and structured low-rank matrix reconstruction methods [35, 56, 97, 110], as well as related subspace methods like ESPIRiT [118]. Another group of methods adopts the same basic structure, but moves from shift-invariant linear prediction to shift- invariant nonlinear prediction relationships to improve performance [1, 20, 64, 87]. And in recent years, due to the growing excitement about machine learning and deep learning in MRI, a huge number of dierent convolutional neural network (CNN) methods have been proposed [67,72]. Each of these methods can be implemented using shift-invariant ltering with appropriate convolution kernels. Most of the existing methods use a rectangular kernel shape, in which the support set of the kernel has a rectangular geometry. For example, the support set of an 1111 rectangular kernel in 2D could be dened as =fk = (n x k x ;n y k y ) : max(jn x j;jn y j) 5g, where k x and k y are the Nyquist sampling intervals along each dimension. In contrast, an ellipsoidal kernel shape would have a support set with an ellipsoidal geometry, and could be dened for an 1111 kernel in 2D as =fk = (n x k x ;n y k y ) : p jn x j 2 +jn y j 2 5g. Al- though we [35,38,64] and other researchers [128] have sometimes used ellipsoidally- shaped , we are not aware of any previous systematic comparisons of rectangular and ellipsoidal kernel shapes. 99 In this work, we perform a detailed evaluation of rectangular versus ellipsoidal kernel shapes, and the results we obtained (to be described later) suggest that el- lipsoidal kernels can indeed oer several advantages over more-common rectangular kernels. A preliminary account of portions of this work was previously presented at a recent conference [78]. 4.2 Theory Many of the earliest convolutional reconstruction methods were designed for 1D reconstruction problems, in which case there is no distinction between a rectangular kernel and an ellipsoidal kernel. For more recent higher-dimensional methods, it is easy to access rectangle-shaped subcomponents of a higher-dimensional array in most programming languages (e.g., in Python, for a 2D array denoted by B, the code B[1:3,1:3] will yield access to a 2D rectangular subcomponent), while accessing an ellipsoid-shaped subcomponent of an array is generally more involved to code. This is likely a major contributing factor to the widespread modern use of rectangular kernels. However, there are some theoretical reasons why an ellipsoidal kernel might be preferred over a rectangular kernel. One reason is that an ellipsoidal kernel can be viewed as a rectangular kernel with the corners removed, which implies that for ellipsoidal and rectangular kernels with \matched size" (i.e., the principal axes of the ellipse have the same lengths as the sides of the corresponding rectangle), an ellipsoidal kernel has smaller area/volume and fewer coecients. In particular, in 100 2D, the area of an ellipse is more than 20% smaller than the area of the rectangle that inscribes it, while in 3D, the volume of an ellipsoid is nearly 50% smaller than the volume of the inscribing hyperrectangle [6]. This dierence in the number of coecients is illustrated for the 2D case in Fig. 4.1. Practically, this means that ellipsoidal kernels are associated with fewer degrees-of-freedom (i.e., fewer coecientsw[z] that need to be learned in the linear case) than rectangular kernels of the same basic size. 3×3 5×5 7×7 9×9 11×11 9 25 49 81 121 5 13 29 49 81 Figure 4.1: Rectangular kernels (rst row) versus ellipsoidal kernels (second row) for dierent kernel sizes. The center of the kernel is marked in red, while other locations within the support are shown in white. For each conguration, the total number of kernel coecients is indicated in yellow in the right-bottom corner. As the kernel size increases, the ellipsoidal kernels have substantially fewer coecients than the corresponding rectangular kernels. Choosing the number of parameters of a reconstruction model represents a clas- sical trade-o. A reconstruction model with too few parameters may not have enough representational capabilities to accurately capture the important features of the desired image to be reconstructed. On the other hand, a model with too many parameters can be prone to overtting, can be sensitive to noise, and can 101 require more training data than models with fewer parameters. As a result, it is important to choose a kernel shape that balances these factors appropriately. Another factor to consider is that reconstruction methods that use convolution in k-space can always be interpreted as using multiplication in the image domain [3,35, 37,86,118,122]. In that sense, the shape of the support set implicitly determines the characteristics of the spatial-domain function corresponding to the convolution kernel w[k]. Interestingly, it is well-established that functions with rectangular k- space support are expected to have highly-anisotropic spatial resolution (with much higher spatial resolution along diagonal lines), while functions with ellipsoidal k- space support are expected to have more isotropic spatial resolution characteristics [6]. This has led some authors to suggest that, in the context of data sampling to achieve a certain target resolution, it suces to acquire an ellipsoidal region of k-space and acquiring the corners of k-space is often inecient [6]. By the same logic, it is reasonable to hypothesize that the corners of rectangular convolution kernels may not be very important in practice, and that ellipsoidal convolution kernels may be able to achieve similar capabilities with better eciency. 4.3 Methods We carried out a systematic comparison between the rectangular and the ellipsoidal kernel shapes, by assessing the performance and eciency of seven representative k- space convolutional reconstruction methods: GRAPPA, SPIRiT, ESPIRiT, SAKE, LORAKS , AC-LORAKS, and a CNN-based approach. 102 In what follows we give a very brief explanation of the selected methods. GRAPPA [28] is a noniterative reconstruction method that interpolates a missing point as a shift-invariant linear combination of neighboring samples (with the set of neighbors dened by the kernel shape) using kernel weights previously estimated from au- tocalibration (ACS) data. SPIRiT [86] is an iterative reconstruction method that uses the constraint that every point in k-space (regardless of whether that point was sampled or missing in the original acquisition) can be predicted as a shift-invariant linear combination of neighboring samples. As before, the set of neighbors is de- ned by the kernel shape, and the kernel weights are estimated using ACS data. The ESPIRiT [118], SAKE [110], LORAKS [35], and AC-LORAKS [36] methods are all based on structured matrices, in which convolution-structured (i.e., Hankel or Toeplitz) matrices are formed from k-space data such that multiplying the struc- tured matrix with a vector is equivalent to a convolution with the k-space data. In all of these methods, the kernel shape in uences the way that the structured matrices are constructed. Each method uses structured matrices in slightly dier- ent ways. ESPIRiT is a sensitivity-map estimation method that forms a structured matrix from ACS data, and obtains sensitivity maps from the eigenvectors of this matrix. SAKE is a reconstruction method that recovers missing k-space samples by enforcing a constraint that a structured matrix formed from the reconstructed k- space data is expected to have low-rank characteristics due to parallel imaging and image support constraints. Similarly, LORAKS uses similar low-rank matrix con- straints to recover missing k-space samples, but constructs the structured matrix 103 and enforces the low-rank constraints dierently than SAKE does. For the results shown in this work, we utilized the version of LORAKS based on the \S-matrix", which simultaneously imposes parallel imaging, image support, and smooth-phase constraints. AC-LORAKS is an autocalibrated version of LORAKS, in which the nullspace of the LORAKS matrix is estimated from ACS data. This allows image reconstruction to be reformulated as a simple least-squares problem instead of a more complicated low-rank matrix recovery problem, greatly improving computa- tion speed. The CNN-based method we considered is the U-Net approach that has been used as a benchmark in the FastMRI challenges [131]. This approach uses a CNN to remove the artifacts from coil-combined images formed through simple zero-lled reconstructions, where the convolution lters are trained using thousands of reference datasets. It should be noted that the U-Net is based on image-domain convolution instead of k-space convolution (which all of the other methods are based on). We assessed the performance of GRAPPA, SPIRiT, SAKE, LORAKS, and AC- LORAKS by reconstructing retrospectively-undersampled data using the ten dif- ferent kernel supports shown in Fig. 4.1, comprising ve rectangular (square) ker- nels with dierent side lengths and ve ellipsoidal (circular) kernels with dierent radii. Two datasets were used: 12-channel T2-weighted brain data acquired with a 256 187 matrix size; and 4-channel T1-weighted brain data acquired with a 258 256 matrix size. These datasets are shown in Fig. 4.2. Retrospective un- dersampling was performed for several dierent sampling patterns as shown in 104 Fig. 4.2. The sampling patterns always include 24 fully-sampled lines at the center of k-space to be used as ACS data. Outside of this region, we used conventional uniform undersampling with acceleration factors ofR = 2; 3, and 4. In order to test both kernel shapes with other undersampling schemes, the T2-weighted dataset was accelerated using a partial Fourier acquisition scheme with the same number of k- space lines used for the T1-weighted dataset. The lines outside the ACS region were distributed proportionally such that only 6/8 of the k-space includes fully-sampled lines. Reconstruction results were evaluated using normalized root-mean-squared error (NRMSE), and we also evaluated the amount of computation time and the amount of RAM used in each case. These results were obtained on a computer with an Intel Xeon E5-1603 2.8 GHz quad core CPU processor and 32GB of RAM. ’ ’ R = 2 R = 3 R = 4 ’ ’ R = 2 R = 3 R = 4 Figure 4.2: The T2-weighted (rst row) and the T1-weighted (second row) datasets with the corresponding sampling patterns used in our experiments. We show results after coil combination (root sum-of-squares). 105 ESPIRiT was evaluated using the same two datasets, using the same ten kernel supports and the same 24 lines of ACS data described previously. The quality of the estimated sensitivity maps was assessed by calculating the size of the residual error (quantied as NRMSE) after projecting the fully-sampled data onto the subspace spanned by the estimated sensitivity maps [118]. The U-Net was evaluated with knee data from the 2019 FastMRI challenge [131]. We compared the original U-Net (which used 3 3 rectangular kernels) against a modied U-Net (which used 3 3 ellipsoidal kernels, which was achieved using 3 3 rectangular kernels but forcing the values in the corners to be zero). Both the original and modied U-Nets were trained for 27 epochs from the same random initialization using 3; 474 single-channel training examples. We used the undersampling scheme used in [131]. For validation, we used 7; 135 single-channel datasets that were not part of the training set. 4.4 Results Figure 4.3 shows the NRMSE results for the T2-weighted dataset obtained for GRAPPA, SPIRiT, ESPIRiT, SAKE, LORAKS, and AC-LORAKS, when using the kernel shapes and kernel sizes shown in Fig. 4.1, with the corresponding results for the T1-weighted dataset shown in Supporting Information Fig. S1. For both datasets it can be observed that, in most cases, the NRMSE achieved by both kernel shapes was quite similar, which can likely be attributed to the fact that both shapes have practically the same base spatial resolution. It should also be noted 106 that for GRAPPA, the 33 rectangular kernel has signicantly better performance than the 3 3 ellipsoidal kernel. This is expected since, for this small kernel size, the ellipsoidal kernel only uses two neighboring k-space samples (one on each side) to interpolate each missing sample, which is much smaller than the six neighboring k-space samples used by the rectangular kernel. For SPIRiT we observed a similar NRMSE for both kernel shapes comparable to the values obtained for the other methods, however, for the T1-weighted dataset, both kernel shapes exhibited a substantial decrease in performance when the kernel size was 11 11. The same phenomenon was observed for the T2-weighted dataset but only for the rectangular shape. This poor performance can be attributed to the fact that these large kernels have a large number of parameters that were dicult to estimate reliably based on the relatively small number of ACS lines. Figure 4.5 and Supporting Information Fig. S2 respectively show corresponding computation times for the T2-weighted and T1-weighted datasets. We observe that the elliptical kernels shapes constitently yielded faster computation times than the rectangular kernel shapes, with generally bigger dierences for larger kernel sizes. This matches the theoretical expectations described previously. Supporting Information Figs. S3 and S4 respectively show the corresponding memory usage for the T2-weighted and T1-weighted datasets. Analogous results are shown in Supporting Information Fig. 4 for the T1-weighted dataset. In all cases, the ellipsoidal kernels used less memory than the rectangular kernels, as expected based on previous theoretical arguments. 107 3×3 5×5 7×7 9×9 11×11 0.1 0.15 0.2 0.25 KernelSize NRMSE GRAPPA R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0.1 0.15 0.2 0.25 KernelSize NRMSE SPIRiT R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0.1 0.15 0.2 0.25 KernelSize NRMSE ESPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0.1 0.15 0.2 0.25 KernelSize NRMSE SAKE R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0.1 0.15 0.2 0.25 KernelSize NRMSE LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0.1 0.15 0.2 0.25 KernelSize NRMSE AC-LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 Figure 4.3: NRMSE for dierent reconstruction methods for the T2-weighted dataset as a function of kernel size and kernel shape. Figure 4.9 shows the histogram of NRMSE values obtained from the validation set of images for CNN-based reconstruction. As can be seen, there is negligible dierence in NRMSE performance for the two dierent kernel shapes, despite the fact that the network based on ellipsoidal kernels had substantially fewer parameters (i.e., 410 6 ) than the network based on rectangular kernels (i.e., 710 6 ). This dierence could be potentially leveraged to reduce the complexity of the training process, although for simplicity, we utilized a naive implementation that did not (and would not be expected to) demonstrate computational benets. 108 3×3 5×5 7×7 9×9 11×11 0 5⋅10 −2 0.1 0.15 0.2 0.25 KernelSize NRMSE GRAPPA R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5⋅10 −2 0.1 0.15 0.2 0.25 KernelSize NRMSE SPIRiT R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5⋅10 −2 0.1 0.15 0.2 0.25 KernelSize NRMSE ESPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5⋅10 −2 0.1 0.15 0.2 0.25 KernelSize NRMSE SAKE R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5⋅10 −2 0.1 0.15 0.2 0.25 KernelSize NRMSE LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5⋅10 −2 0.1 0.15 0.2 0.25 KernelSize NRMSE AC-LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 Figure 4.4: NRMSE for dierent reconstruction methods for the T1-weighted dataset as a function of kernel size and kernel shape. 3×3 5×5 7×7 9×9 11×11 0 10 20 30 40 KernelSize Time(sec) GRAPPA R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 0 10 20 30 40 KernelSize Time(sec) SPIRiT R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5 10 15 KernelSize Time(sec) ESPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 1,000 2,000 3,000 KernelSize Time(sec) SAKE R= 2 Elliptical R= 2 Rectangular R= 3 Elliptical R= 3 Rectangular R= 4 Elliptical R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 200 400 600 KernelSize Time(sec) LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5 10 15 20 KernelSize Time(sec) AC-LORAKS R= 2 Elliptical R= 2 Rectangular R= 3 Elliptical R= 3 Rectangular R= 4 Elliptical R= 4 Rectangular 1 Figure 4.5: Computation times for dierent reconstruction methods for the T2-weighted dataset as a function of kernel size and kernel shape. 109 3×3 5×5 7×7 9×9 11×11 0 1 2 3 4 5 KernelSize Time(sec) GRAPPA R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5 10 15 KernelSize Time(sec) SPIRiT R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 2 4 6 KernelSize Time(sec) ESPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 500 1,000 1,500 KernelSize Time(sec) SAKE R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 20 40 60 80 100 120 KernelSize Time(sec) LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 5 10 15 KernelSize Time(sec) AC-LORAKS R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 Figure 4.6: Computation times for dierent reconstruction methods for the T1-weighted dataset as a function of kernel size and kernel shape. 3×3 5×5 7×7 9×9 11×11 0 0.1 0.2 0.3 0.4 0.5 KernelSize Megabytes GRAPPA R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 0.1 0.2 0.3 KernelSize Megabytes SPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 1 2 3 4 5 KernelSize Megabytes ESPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 200 400 600 800 1,000 KernelSize Megabytes SAKE Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 500 1,000 1,500 2,000 KernelSize Megabytes LORAKS Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 500 1,000 1,500 2,000 KernelSize Megabytes AC-LORAKS Ellipsoidal Rectangular 1 Figure 4.7: Memory usage for dierent reconstruction methods for the T2-weighted dataset as a function of kernel size and kernel shape. 110 3×3 5×5 7×7 9×9 11×11 0 0.5 1 1.5 2 ×10 −2 KernelSize Megabytes GRAPPA R= 2 Ellipsoidal R= 2 Rectangular R= 3 Ellipsoidal R= 3 Rectangular R= 4 Ellipsoidal R= 4 Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 1 2 3 ×10 −2 KernelSize Megabytes SPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0.5 1 1.5 KernelSize Megabytes ESPIRiT Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 100 200 300 400 500 KernelSize Megabytes SAKE Elliptical Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 200 400 600 800 1,000 KernelSize Megabytes LORAKS Ellipsoidal Rectangular 1 3×3 5×5 7×7 9×9 11×11 0 200 400 600 800 1,000 KernelSize Megabytes AC-LORAKS Ellipsoidal Rectangular 1 Figure 4.8: Memory usage for dierent reconstruction methods for the T1-weighted dataset as a function of kernel size and kernel shape. 111 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 0 100 200 300 400 500 NRMSE #Images NRMSEHistograms Ellipsoidal Rectangular 1 Figure 4.9: Histograms of the NRMSE values obtained on the validation set by the CNN-based recon- struction for the dierent kernel shapes. 112 Discussion and Conclusions To the best of our knowledge, this work represents the rst systematic evaluation of kernel shapes for convolution-based MRI reconstruction methods. Although ellipsoidal kernels are not very popular in the modern literature, our results suggest that ellipsoidal kernels often oer advantages (in computation time, memory usage, and the number of model parameters) over more-common rectangular kernels, with largely similar NRMSE performance. We expect that these new insights may be valuable for improving the eciency of MRI reconstruction in the future. In our evaluations, we have focused only on square-shaped and circle-shaped kernels, without taking full advantage of the degrees of freedom oered by rectangles and ellipses (with dierent sizes along dierent axes). We anticipate that additional improvements might be obtained by leveraging this untapped degree of freedom. It should also be noted that the ellipsoidal kernel shape is not necessarily optimal, and that other shapes (perhaps learned directly from empirical data) might achieve even better results. Exploration of such ideas is a potentially interesting direction for further research, though is beyond the scope of the present work. 113 Chapter 5 New Theory and Faster Computations for Subspace-Based Sensitivity Map Estimation in Multichannel MRI 5.1 Introduction Most modern magnetic resonance imaging (MRI) experiments are performed using multichannel phased-array receiver coils [10,18,25,28,49,50,59,68,69,86,102{104, 107,112,129,132]. For acquisitions involvingQ receiver channels, the measurements obtained in a multichannel MRI experiment are usually modeled using the Fourier transform as d q (k m ) = Z c q (x)(x)e i2k T m x dx + qm (5.1) for m = 1;:::;M and q = 1;:::;Q. In this expression, M represents the total number of measured k-space positions k m 2 R D , where D is the dimension of the image; 1 d q (k m ) represents the complex-valued k-space data measured from the 1 For simplicity, we have adopted the use of \spatial" notation throughout this paper, which is consistent with typical use cases in 2D or 3D imaging. In these cases, we would have eitherD = 2 withx = [x;y] T andk = [kx;ky ] T orD = 3 withx = [x;y;z] T andk = [kx;ky;kz ] T . However, it should be noted that our approach is also compatible with other scenarios such as spatiotemporal imaging with time-varying sensitivities, which, for example, might use D = 3 with x = [x;y;f] T and k = [kx;ky;t] T . 114 qth receiver channel at the mth position in k-space; (x) is the complex-valued underlying MR image (re ecting the state of the excited magnetization at the time of data acquisition) as a function of the spatial position x2R D ;c q (x) is a complex- valued sensitivity map describing the sensitivity of the qth receiver channel to excited magnetization at spatial position x; and qm represents complex-valued measurement noise. We will also use q (x),c q (x)(x) to denote theqth sensitivity- weighted image. This paper is focused on the problem of estimating the sensitivity proles c q (x) from a set of calibration measurements (i.e., a set of k-space data samples that are sampled at the Nyquist rate). While prior knowledge of the coil sensitivity maps is not always required, there are many tasks that can benet from estimated sensitivity maps, including coil combination (in which the underlying MRI image (x) is estimated from the multichannel images q (x) [16,107,121]) and accelerated image reconstruction (in which the underlying MRI image (x) is estimated from multichannel k-space data d q (k m ) that is sampled below the Nyquist rate [10, 25, 69,102,103,112,129]). Over the years, many dierent sensitivity map estimation methods have been proposed [2, 16, 29, 45, 57, 96, 99, 103, 108, 112, 117, 118, 121, 130]. Among these, subspace-based approaches [29, 96, 108, 118] have proven to be particularly power- ful and popular, with the subspace-based ESPIRiT method [118] rising to become the most widely-used sensitivity map estimation method in the modern literature because of its simplicity (with few tuning parameters) and excellent performance. 115 The mathematical principles underlying subspace-based sensitivity map estima- tion are nontrivial and are not always easy to understand. In the rst part of this work, we present a novel theoretical derivation of subspace-based sensitivity map estimation that is based on a nullspace/linear predictability perspective [37], and makes use of concepts from the recent literature on structured low-rank matrix modeling [35{38,53,56,97,110]. This approach is distinct from and complementary to existing theoretical explanations, and we personally nd it to be more intuitive than the alternatives. In the second part of this work, we introduce several new computational ap- proaches that can substantially reduce the time and memory required for subspace- based sensitivity map estimation. This can be important because existing subspace- based methods can be computationally demanding, with substantial memory re- quirements and slow computation times. While slow computations can already be problematic when considering individual datasets, they can be especially limiting in modern machine learning contexts where it may be of interest to calculate sen- sitivity maps for every image within a database containing hundreds or thousands of datasets. Our computational contributions include the following: • We leverage our previous observation that using ellipsoidal convolution kernels instead of rectangular convolution kernels can lead to substantial reductions in memory usage and computational complexity, without sacricing accuracy [81]. 116 • We introduce a fast FFT-based approach that allows us to rapidly calculate the nullspace vectors associated with a multichannel convolution-structured matrix, without ever having to directly construct the actual convolution- structured matrix (which would normally require substantial memory and time to assemble). This approach can be viewed as a multichannel gener- alization of an approach that was previously developed for the single-channel case [98]. In addition to being useful for subspace-based sensitivity map es- timation, this type of approach may also be useful for the broad class of structured low-rank modeling methods for multichannel MRI image recon- struction [35{38,53,56,97,110]. Note that a similar FFT-based approach was proposed for multichannel MRI reconstruction in Ref. [133]. Compared to that approach, our implementation uses substantially fewer FFTs, which is enabled by making use of substantial redundancies that are present in the convolution-structured matrix. • Noting that sensitivity maps arise from Maxwell's equations and must be smooth, we observe that estimating sensitivity maps at the nominal desired spatial resolution can be computationally wasteful. Instead, we propose a method that rst estimates the sensitivity maps on a low-resolution spatial grid, and then uses a fast FFT-based spatial interpolation procedure to obtain sensitivity maps at the desired spatial resolution. 117 • Subspace-based sensitivity map estimation methods often require assembling a large number of small-scale matrices (e.g., one matrix for every spatial loca- tion in the image) which in total require substantial memory. We propose a memory-ecient approach based on signal processing techniques to calculate these matrices. • Subspace-based sensitivity map estimation methods and coil combination meth- ods often require computing a large number of small-scale singular value de- compositions (SVDs) (e.g., one SVD for every spatial location in the im- age) [118,121]. We introduce an ecient iterative approach (based on classical power iteration [26]) to calculate many small-scale partial SVDs simultane- ously. While this approach is directly useful for subspace-based sensitivity map estimation, we expect that it may also have some utility for a range of locally low-rank modeling and denoising methods that also require calculat- ing partial SVDs over local image patches for every spatial location in the image [24,94,115,116,120]. Combined together, these algorithmic improvements | which we collectively call PISCO (Power iteration over simultaneous eciently-assembled patches, Interpo- lation, ellipSoidal kernels, and FFT-based COnvolution) | enable a substantial (roughly 100-fold in some cases) computational acceleration for subspace-based sensitivity map estimation. 118 5.2 Characteristics of the Inverse Problem Before describing our novel theory and methods, we will rst review some important characteristics of the sensitivity map estimation problem. As mentioned in the previous section, the problem we consider in this paper is the estimation of the sensitivity mapsc q (x) forq = 1;:::;Q from a set of Nyquist- sampled calibration measurements. We will assume that these calibration measure- ments consist of a collection of M k-space samples from each of the Q channels, for a total ofMQ observations. Our goal will be to calculate the values of the sen- sitivity maps on a discrete voxel grid of N spatial locations x n , for n = 1;:::;N, resulting in a total of NQ unknowns. Most practical scenarios will involve sub- stantially more spatial locations than k-space samples (i.e., NQ MQ), such that the inverse problem is underdetermined (with fewer measurements than un- knowns) unless additional assumptions are made about the characteristics ofc q (x). To overcome this issue, the literature on sensitivity map estimation typically makes the assumption that the sensitivity maps c q (x) are spatially smooth, which can be justied based on physical principles (e.g., Maxwell's equations). A more pernicious issue is that, while our goal is to estimate the unknown sen- sitivity maps c q (x), the image (x) that appears in Eq. (5.1) is also generally an unknown. This means that Eq. (5.1) can be viewed as a bilinear inverse problem, which is problematic because bilinear inverse problems have well-known scaling ambiguities. In particular, given one set of estimated sensitivity maps c q (x) for q = 1;:::;Q and a corresponding estimated image (x) that together provide a 119 good t to the measured data, it is straightforward to observe that another equally- good t to the measured data can be obtained by combining a rescaled version of the sensitivity maps (x)c q (x) for q = 1;:::;Q together with a rescaled version of the image 1 (x) (x). The complex-valued scaling function (x) must be nonzero for every point x within the eld of view (FOV), but is otherwise completely un- constrained, and can have arbitrary spatial variations. As such, unless additional strong assumptions are made, there will be innitely many ways of choosing (x) that all lead to sensitivity map estimates that provide equally-good ts to the mea- sured data. Practically, this means that both the magnitude and phase of the estimated sensitivity maps are inherently ambiguous, and we can only hope for sensitivity map estimates that are unique up to these ambiguities. Because of the inherent ambiguity, it is common for sensitivity map estimation methods to search for one set of estimated sensitivity maps c q (x) that t the data well, and subsequently rescale them with a heuristically-chosen scaling function (x) that can be designed to endow the scaled sensitivity maps (and the corre- sponding scaled image) with properties that are deemed to be useful. Note also that the sensitivity map valuesc q (x) will have no impact on measured data at spatial locations where (x) = 0. As a result, the values of the sensitivity maps are completely unidentiable at spatial locations outside the support of(x), leading to additional ambiguity. 120 5.3 Nullspace/Linear Predictability Theory for Subspace- Based Sensitivity Map Estimation In this section, we present a novel derivation of subspace-based sensitivity map estimation from a nullspace/linear predictability perspective [37]. We will start with a quick summary of the relevant linear predictability relationships that are known to exist for multichannel MRI, before making novel links to sensitivity map estimation. 5.3.1 Summary of Linear Predictability for Multichannel MRI Dene the spatial support of(x) as the set of points =fx2R D :j(x)j> 0g. In what follows, we will assume (without loss of generality) that the spatial coordinate system has been normalized such that the support is completely contained within a hypercube with sides of length one, dened as ,fx2R D :kxk 1 < 1 2 g. We will refer to as the FOV. This assumption allows us to equivalently represent the continuous images (x) and q (x) using innite Fourier series representations as (x) = 1 (x) X n2Z D s[n]e i2n T x (5.2) and q (x) = 1 (x) X n2Z D s q [n]e i2n T x (5.3) forq = 1;:::;Q, wheres[n] ands q [n] respectively represent samples of the Fourier transforms of (x) and q (x) on the rectilinear Nyquist grid (taking k = 1 along 121 each dimension, as enabled by our assumptions about ), and 1 (x) is the indicator function for the FOV, with 1 (x) = 8 > > < > > : 1; x2 0; else: (5.4) A basic assumption of the linear predictability formalism [37] is that the multi- channel collection of sequencesfs q [n]g Q q=1 will be autoregressive, satisfying multiple multichannel shift-invariant linear predictability relationships [23, 28, 29, 35{38,53, 56, 86, 96, 108, 110, 118, 132]. Specically, there should exist multiple dierent Q- channel nite impulse response (FIR) lters h r [n;q] that each satisfy Q X q=1 X m2 h r [m;q]s q [n m] 0 for8n2Z D (5.5) for r = 1;:::;R, where R is the number of dierent lters and Z D is a nite index set describing the support of each FIR lter. There are many dierent situations in which the existence of such ltersh r [n;q] has been theoretically proven [37]. We review a few of these below: 2 • Limited Image Support [23,35,37]. In many situations, the actual spatial support of the original image(x) will be smaller than the FOV , such that the complement of in (denoted by { , n ), is a measurable nonempty set. In such cases, it is possible to dene nonzero functions f(x) that have 2 Note that shift-invariant linear predictability relationships can also be proven for images that possess smooth phase [35, 37, 47] and sparsity characteristics [37, 53, 56, 73, 75, 97], although these relationships are not directly relevant for sensitivity map estimation and we do not discuss them further. 122 their supports contained within { . Since(x) andf(x) have nonoverlapping support, we must have f(x)(x) = 0 for8x2 . If we let ~ f[n] denote the Nyquist-sampled Fourier series representation off(x), then taking the Fourier transform of f(x)(x) and applying the Fourier convolution theorem yields that X m2Z D ~ f[m]s[n m] = 0 for8n2Z D : (5.6) If { is large enough, then it becomes possible to choose relatively smooth functions f(x) that satisfy the conditions specied above. If we assume that these functions are smooth enough that they satisfy ~ f[n] 0 for8n = 2 (i.e., the functionsf(x) are approximately bandlimited to the frequency range dened by ), 3 then we can approximate ~ f[n] as an FIR lter, and obtain X m2 ~ f[m]s[n m] 0 for8n2Z D : (5.7) Note also that if we have f(x)(x) 0 for8x2 , then the relationship between (x) and q (x) implies we must also have that f(x) q (x) 0 for 8x2 for q = 1;:::;Q. This implies that the multichannel FIR lter h r [n;q] = 8 > > < > > : q ~ f[n]; n2 0; else (5.8) 3 While we can often obtain very accurate approximations of this form, this must always be an approximation, and will never be exact. Specically, since f(x) is assumed to have nite support, then ~ f[n] cannot have exactly- limited support by the Fourier uncertainty principle [11]. While this paper describes this approximation informally to simplify the exposition, more formal treatments exist in the literature [37]. 123 will satisfy Eq. (5.5), for any choices of q 2 C for q = 1;:::;Q. It is easy to obtain many dierent h r [n;q] lters from this relationship because of the exibility oered by choosing dierent functionsf(x) and dierent coecients q . Using similar arguments, additional channel-specic lters can also be obtained in common situations where some sensitivity-weighted images q (x) have smaller support than the underlying image (x). • Smooth Sensitivity Maps [37, 53, 96, 108]. As mentioned previously, it is often assumed that the sensitivity maps c q (x) are spatially smooth. Linear predictability relationships can be derived if we make the assumption that they are smooth enough that they can be treated as approximately bandlimited to the frequency range dened by , such that c q (x) 1 (x) X n2 ~ c q [n]e i2n T x (5.9) for appropriate Fourier coecients ~ c q [n]. In this case, it can be observed that c i (x) j (x)c j (x) i (x) = 0 for8x2 for any choices of i and j, which is sometimes called the \cross-relation" in the multichannel blind deconvolution literature [30,39,126]. Taking the Fourier transform of the cross-relation and applying the Fourier convolution theorem leads to X m2 ~ c i [m]s j [n m] X m2 ~ c j [m]s i [n m] 0 for8n2Z D : (5.10) 124 This implies that the multichannel FIR lter h r [n;q] = 8 > > > > > < > > > > > : ~ c j [n]; n2 ;q =i ~ c i [n]; n2 ;q =j 0; else (5.11) will satisfy Eq. (5.5) for any choices ofi andj withi6=j. As before, it is easy to obtain many dierent h r [n;q] lters using dierent choices of i and j. • Compressible Receiver Arrays [14,48,60]. It is frequently observed that the dierent channels in a multichannel receiver array possess linear depen- dencies, such that there exist multiple distinct choices of linear combination weights w q that will each satisfy Q X q=1 w q q (x) 0 for8x2 : (5.12) This observation is often used for dimensionality reduction (\coil compres- sion") [14, 48, 60]. However, this observation also implies that a relationship in the form of Eq. (5.5) will hold if h r [n;q] is chosen as h r [n;q] = 8 > > < > > : w q ; n = 0 0; else: (5.13) Since there often exist multiple distinct choices of w q that satisfy Eq. (5.12), there can also be many dierent lters h r [n;q] that arise from this kind of relationship. 125 The preceding arguments all support the existence of multiple lters h r [n;q] that satisfy Eq. (5.5). In addition, if we were given full information about the support and the sensitivity mapsc q (x), then we would be able to calculate lters h r [n;q] directly from Eqs. (5.8), (5.11), and (5.13). However, in practice, we do not have information about or c q (x) (these are quantities that we are interested in estimating!), and we must instead resort to other methods to obtain the FIR lters h r [n;q]. Towards this end, an important observation from the literature on structured low-rank matrix modeling [35{38, 53, 56, 97, 110] is that the convolution appearing in Eq. (5.5) can be represented in matrix form as Ch 0, where the matrix C is given by C, C 1 C 2 ::: C Q ; (5.14) the matrix C q 2C Pjj is a convolution-structured matrix (i.e., a Hankel or Toeplitz matrix) formed from Fourier samples s q [n], with each row of C q corresponding to one value of n from Eq. (5.5);P represents the number of dierent n values that are used in the construction of each C q matrix; and h2C Qjj is the vector of h r [n;q] samples. See Ref. [38] for a detailed description of the matrix construction we use. Even though the lter coecients h may be a priori unknown, all ltersh r [n;q] that satisfy Eq. (5.5) must appear as approximate nullspace vectors of an appropriately- constructed C matrix. Importantly, it is possible to directly form such a C matrix 126 from Nyquist-sampled calibration data, and its approximate nullspace can be iden- tied by applying the SVD. It is therefore possible to learn a set ofR linearly inde- pendent ltersh r [n;q] (withR equal to the dimension of the approximate nullspace of C) in a data-driven way, without requiring any additional prior information [36]. Note that this nullspace-based method of computation will not result in lters that directly match the forms given in Eqs. (5.8), (5.11), or (5.13). Instead, the obtained nullspace vectors are likely to be linear combinations of such lters. In addition, the full set of approximate nullspace vectors for an appropriately-constructed C matrix is expected to form a basis for the set of all h r [n;q] lters that are consistent with the FIR lter support that was used in the construction of C. 5.3.2 Linear Predictability and Sensitivity Map Estimation In what follows, we assume that we have access toR ltersh r [n;q] that each satisfy Eq. (5.5), which can, e.g., be obtained from the nullspace of a structured calibration matrix as described in the previous subsection. Observe that the linear prediction relationship from Eq. (5.5) can be rewritten in the image domain as Q X q=1 h r (x;q) q (x) 0 for8x2 ; (5.15) where h r (x;q), X n2 h r [n;q]e i2n T x : (5.16) 127 Using the relationship between(x) and q (x), Eq. (5.15) can be further simplied to (x) Q X q=1 h r (x;q)c q (x) 0 for8x2 : (5.17) Sincej(x)j> 0 for8x2 (by the denition of ), this implies that the summation Q X q=1 h r (x;q)c q (x) 0 for8x2 ; (5.18) demonstrating that the sensitivity maps c q (x) and the spatial-domain representa- tion of the ltersh r (x;q) are strongly related to each other, and that must hold for all x2 . On the other hand, the summation in Eq. (5.18) could result in arbitrary nonzero values for x2 { , since the fact that (x) = 0 at these spatial locations will still ensure that Eq. (5.17) is satised regardless of the properties of h r (x;q) and c q (x). Since we haveR ltersh r [n;q] that each satisfy Eq. (5.5), the relationship from Eq. (5.18) can be expressed simultaneously for all R lters in matrix form as H(x)c(x) 0 for8x2 ; (5.19) where H(x)2 C RQ has entries [H(x)] rq = h r (x;q), and c(x)2 C Q is the vector of c q (x) values. Importantly, this demonstrates that the sensitivity map values at location x2 will be one of the approximate nullspace vectors of the corresponding matrix H(x). 128 Ideally, the matrix H(x) will only have a single approximate nullspace vector for each spatial position x2 , since that would allow us to estimate c q (x) for x2 by calculating a basis for the approximate nullspace of H(x) (which can, e.g., be done using the SVD). Note that singular vectors associated with distinct singular values are always unique up to scaling ambiguities, which would allow us to obtain unique sensitivity map valuesc q (x) up to the same type of scaling ambiguities that are inherent to sensitivity map estimation (see previous discussion in Sec. 5.2). Importantly, the ability to estimate unique sensitivity mapsc q (x) (up to scaling ambiguities) for x2 is predicated on the matrix H(x) being approximately rank Q1. A necessary condition for this to be true is thatR>Q. However, havingR> Q is not sucient to ensure that H(x) possesses the desired rank characteristics. For example, lters associated with limited image support (Eq. (5.8)) will satisfy h r (x;q) 0 for8x 2 , and will therefore not contribute signicantly to the rank of H(x). In practice however, we have observed empirically 4 that if we use the entire collection of lters from the approximate nullspace of C, then the H(x) matrices always have an approximate rank ofQ1 for x2 under the assumptions made in this paper. 5 Interestingly, when x2 { (i.e., the situation where the sensitivity maps are unidentiable) and we use the entire collection of lters from 4 Our experience is based on identifying the approximate nullspace by applying a heuristic thresholding rule to the SVD. Specically, ifn denotes thenth singular value ofC (with singular values indexed in order from largest to smallest), then our thresholding rule includes the nth right singular vector within the approximate nullspace whenevern < 0:05 1 . The threshold we chose of 0:05 1 was determined heuristically and is unlikely to be optimal, although worked well for the datasets considered in this work. 5 Note that dierent rank characteristics may be observed under dierent assumptions. For example, if the image support is larger than the FOV and if aliasing occurs within the FOV, then we empirically observe multiple sets of nullspace vectors at the spatial locations where the image support aliases onto itself. This is expected since such aliased voxels will be a linear combination of multiple sets of sensitivity-weighted images with distinct sensitivity maps. This enables the computation of multiple sensitivity maps, similar to what has been described previously in work on ESPIRiT [118], although we do not belabor this point as this scenario is not our primary interest. 129 the approximate nullspace of C, we empirically observe that the matrix H(x) always has rank Q rather than Q 1, which provides a simple mechanism for identifying both and { . In what follows, we will use G(x)2 C QQ to denote the positive semidenite matrix formed from H(x) as G(x), (H(x)) H H(x). The G(x) matrix is useful because it can be much smaller than H(x) in the common case that RQ, while the two matrices share many characteristics. Specically, the singular values and right singular vectors of H(x) have a one-to-one correspondence with the sorted squared eigenvalues and eigenvectors of G(x). Empirical demonstrations of the nullspace characteristics of H(x) are presented in the next subsection. 5.3.3 Empirical Demonstration of Nullspace Characteristics Throughout this paper, we will use the following three multichannel MRI datasets to illustrate theoretical principles and to validate our proposed methods: • Brain MPRAGE. We used an MPRAGE sequence on a 3T scanner to ac- quire a T1-weighted human brain image using a 32-channel receiver array, as depicted in Fig. 5.1(a). • Brain TSE. We used a turbo spin-echo sequence on a 3T scanner to acquire a T2-weighted human brain image using a 32-channel receiver array, as depicted in Fig. 5.1(b). 130 (a) Brain MPRAGE (b) Brain TSE (c) Knee TSE Figure 5.1: Depiction of the three datasets we use for illustration and validation. • Knee TSE. We used a 15-channel knee dataset from the fastMRI database [131], as depicted in Fig. 5.1(c). This image was acquired using a fat-suppressed turbo spin-echo sequence on a 3T scanner, producing a proton-density weighted image. To illustrate the low-rank characteristics of the C matrix as described in Sec. 5.3.1, Fig. 5.2 plots the singular values of C for all three datasets. In each case, the C matrix was calculated assuming a 7 7 rectangular FIR lter support (i.e., =fn2 Z 2 :knk 1 3g), and the central 2424 k-space samples from the Nyquist grid. As expected from theory, the singular values decay rapidly in each case, and all three C matrices are approximately low-rank with a substantial ap- proximate nullspace. To illustrate the characteristics of the H(x) and G(x) matrices described in Sec. 5.3.2, Fig. 5.3 shows spatial maps of the Q = 32 eigenvalues of G(x) for the Brain TSE dataset. As expected, the smallest eigenvalue is approximately zero within the support of the image (i.e., H(x) and G(x) are approximately rank 131 0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 1 n σ n 1 (a) Brain MPRAGE 0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 1 n σ n 1 (b) Brain TSE 0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 1 n σ n 1 (c) Knee TSE Figure 5.2: Normalized singular value curves for theC matrices corresponding to the three datasets from Fig. 5.1. Figure 5.3: Spatial maps of the 32 eigenvalues of the G(x) matrices for the Brain TSE dataset. The eigenvalues have been normalized byjj and sorted from largest to smallest, with the largest in the top right and the smallest in the bottom right. =Q 1 for x2 ), while the matrices have full rank (i.e., rank =Q) outside the image support. 5.3.4 Naive Algorithm for Nullspace-Based Sensitivity Map Estimation Based on the theoretical framework described in the previous subsections, a naive nullspace-based algorithm for calulating sensitivity maps would proceed as follows: 1. Choose an FIR lter support set , and construct the corresponding C matrix from calibration data. 2. Calculate a set of R approximate nullspace vectors of the C matrix using the SVD. (In our implementation, we will compute the SVD of C H C rather than the SVD of C, since the two matrices have the same nullspace while C H C is generally much smaller.) 132 Nullspace ESPIRiT Nullspace + PISCO Figure 5.4: Magnitude of sensitivity maps for 16 representative coils estimated using the nullspace-based method (rst row), ESPIRiT (second row), and the nullspace-based method combined with the PISCO techniques (third row). 3. Using the nullspace vectors from Step (2), calculate the H(x) matrices for each spatial position x, which can be simplied into the G(x) matrices as described in the previous subsection. 4. Calculate the eigenvalue decomposition of G(x) for each spatial location x, and set c q (x) equal to the qth entry of the eigenvector associated with the smallest eigenvalue of G(x). 5. (Optional) If desired, choose a scaling function(x) that has useful properties, and use it to normalize the sensitivity maps. 6. (Optional) If desired, apply thresholding to the spatial map of the smallest eigenvalue of G(x) to obtain a binary mask for the image support . As we will demonstrate in the sequel, this naive algorithm works well although can require substantial computation and memory usage, while a modied algorithm that uses our proposed PISCO approach leads to substantial computational accelerations and memory savings. 133 5.3.5 Memory Requirements of the Nullspace Algorithm and the Direct Calculation of G(x) One of the limiting factors in the previously described nullspace-based algorithm occurs in Step (3), where it is necessary to calculate the R spatial-domain mul- tichannel lters h r (x;q) in order to construct the H(x) matrices for each spatial location x. This can be computationally demanding in terms of memory usage when R is big, which is commonly observed in practice. To overcome this limita- tion, we propose a computational procedure to directly calculate the matrices G(x) for each spatial location x without constructing H(x) rst. This procedure corre- spond to the rst of the proposed PISCO techniques. Before describing how this task is carried out, we introduce the matrix V2C QjjR which columns correspond to the approximate nullspace vectors of C or, equivalently, to the Fourier-domain multichannel lters satisfying Eq. (5.5). Specically, V has the form V, 2 6 6 6 6 6 6 6 6 6 4 h (1) 1 h (1) 2 ::: h (1) R h (2) 1 h (2) 2 ::: h (2) R . . . . . . ::: . . . h (Q) 1 h (Q) 2 ::: h (Q) R 3 7 7 7 7 7 7 7 7 7 5 ; (5.20) where h (q) r 2 C jj is the vectorized version of h r [n;q]; n2 . Now we describe our ecient procedure to calculate G(x). The entries of this matrix are given by [G(x)] pq = P R r=1 h r (x;p)h r (x;q), where h r corresponds to the complex conjugate 134 version of h r . This expression can also be calculated using the vectorized Fourier- domain lters as follows [G(x)] pq = f H (x) R X r=1 h (q) r (h (p) r ) H ! f(x); (5.21) where f(x), [e i2n T 1 x e i2n T 2 x e i2n T jj x ] T is the vector of complex sinusoids associated to the frequency componentsfn s g jj s=1 within . Equation (5.21) is im- portant since the matrix W qp , P R r=1 h (q) r (h (p) r ) H 2C jjjj can be easily extracted from the projection matrix VV H , which is not expensive to calculate. Then, if we denote thesth column of W qp by w (s) qp , we have that W qp = P jj s=1 w (s) qp e H s , where e s corresponds to the canonical vector with only zeros except for the sth entry which is equal to one. This allows us to write [G(x)] pq = f H (x) 0 @ jj X s=1 w (s) qp e H s 1 A f(x) = jj X s=1 w (s) q (x;p)e i2n T s x ; (5.22) where w (s) q (x;p) corresponds to the spatial-domain lter associated to w (s) qp . Alter- natively, we can perform the previous calculation in the Fourier domain. Denoting by w (s) q [n;p] the Fourier representation of w (s) q (x;p) (which nonzero components are the entries of w (s) qp ) and byF 1 the inverse Fourier transform operator, we have that [G(x)] pq corresponds toF 1 P jj s=1 w (s) q [n + n s ;p] at location x. Therefore, in order to calculate the matrices G(x) for each location x, we only need to store Q spatial-domain multichannel lters (i.e.,w (s) q (x;p) for8q2f1;:::;Qg), which is 135 usually substantially smaller than R. It should also be noted that a more ecient calculation of VV H is possible using the row-space vectors of C, which are usually considerably less than the approximate nullspace vectors. In fact, we have that VV H = I Qjj UU H , where U is the matrix which columns correspond to the row-space vectors, and I Qjj corresponds to theQjjQjj identity matrix. In the following we consider this approach to calculate the projection matrix VV H , given that calculating UU H is faster than calculating VV H . In the following experiments the proposed nullspace-based algorithm is imple- mented considering the previously described approach to calculate the G(x) ma- trices in Step (3). Otherwise, the memory requirements become restrictive making the nullspace-based algorithm unfeasible. 5.3.6 Comparisons of the Nullspace Algorithm and ESPIRiT Before moving on to describe the rest of our novel PISCO computational accelera- tion techniques, we will rst spend a little time to describe how PISCO compares to the popular ESPIRiT approach for sensitivity map estimation [118], which is also subspace-based. While the two approaches are derived using very dierent argu- ments, the resulting approaches end up having quite similar structure. Specically, in our nullspace-based approach, we use the nullspace vectors of the C matrix to identify matrices G(x) that satisfy G(x)c(x) 0, where c(x)2 C Q is the vector of c q (x) sensitivity map values. This causes us to estimate sensitivity maps using 136 24×24 32×32 64×64 128×128 256×256 0 0.1 0.2 0.3 0.4 0.5 SizeofCalibrationData NormalizedProjectionResidual Nullspace ESPIRiT 1 Figure 5.5: Quantitative comparison between the nullspace-based approach and ESPIRiT for the Brain TSE data in terms of (a) accuracy (as measured through the normalized projection residual) and (b) speed (as measured through the total computation time). the eigenvector associated with the smallest eigenvalue of G(x). In contrast, ES- PIRiT extracts the rowspace vectors of the C matrix, and uses these to construct matrices B(x)2C QQ that are expected to satisfy B(x)c(x) c(x). As a result, the sensitivity maps in ESPIRiT can be estimated using the eigenvector associated with the largest eigenvalue of B(x). Interestingly, it can be derived (with a sub- stantial amount of tedious-but-straightforward derivations that are too involved to reproduce here) that the B(x) matrix used by ESPIRiT will satisfy B(x) = I Q G(x) jj : (5.23) In addition to theoretical similarities between ESPIRiT and the nullspace-based approach, the two approaches also empirically produce sensitivity map estimates that are nearly identical to one another. Results are illustrated in Fig. 5.4, where we show representative sensitivity maps estimated from the Brain TSE dataset. 137 For both methods we chose an FIR lter support of dimensions 7 7 and a cali- bration data set of dimensions 24 24 to calculate C. We calculated the singular vectors of C H C corresponding to singular values which were greater than the 5% of the greatest singular value. These singular vectors are supposed to approximate a basis for the rowspace of C H C, and can be used to calculate the matrices B(x) for each spatial location. Analogously, we calculated the singular vectors of C H C corresponding to singular values smaller than the 5% of the greatest singular value (i.e., approximate nullspace vectors), and use them to directly calculate the matri- ces G(x) (i.e., using the method described in Sect. 5.3.5). Then, we estimated the sensitivity maps at each spatial location by calculating the rst eigenvector of B(x) in ESPIRiT, and the last eigenvector of G(x) in the nullspace-based algorithm. Next, the estimated sensitivity maps for both methods were normalized in magni- tude using a sum-of-squares approach, and referenced in phase with respect to the phase of the sensitivity map estimate corresponding to the rst channel. Finally, both sets of sensitivity maps were masked by identifying locations where the small- est eigenvalue of G(x) was smaller than the 8% of the greatest eigenvalue. From the results in Fig. 5.4 it can be observed that the visual dierences between the two approaches are negligible. Figure 5.4 also shows that the results are negligibly dierent from those obtained when combining the nullspace-based approach with PISCO (which leads to dramatic improvements in computation speed and memory usage), although we will save that discussion for later in the paper after we have discussed the remaining details of PISCO. 138 To evaluate the quality of a sensitivity map estimate in a more quantitative way, we have also computed a \normalized projection residual" to quantify how well the estimated sensitivity maps match the fully-sampled data. Specically, we rst take the fully-sampled multichannel data and project it onto the subspace spanned by the sensitivity maps. If the sensitivity map estimation procedure worked well, then the projection onto the sensitivity maps should capture all of the true signal, with only noise remaining in the residual after projection [118]. Our \normalized projection residual" metric is obtained by taking the ` 2 -norm of the dierence between the original multichannel data and the projected data, and normalizing by the ` 2 -norm of the original multichannel data. A smaller value of this metric is indicative of better sensitivity map estimation performance. Figure 5.5(a) shows plots of this metric for the Brain TSE data, for both ESPIRiT and the nullspace- based approach, using a 77 FIR lter support and dierent calibration data sizes. As can be seen, both methods have very similar residual values. Table 5.1 also shows the normalized projection residual for all three datasets for both ESPIRiT and the nullspace-based approach, this time using 3232 calibration data while still using a 77 FIR lter support. As before, the dierences between ESPIRiT and the nullspace-based approach are negligible. (This table also shows PISCO results, although we will defer a detailed discussion of these results until after we have described PISCO). 139 Table 5.1: Normalized Projection Residual with 3232 calibration data. Dataset ESPIRiT ESPIRiT Nullspace Nullspace +PISCO +PISCO Brain MPRAGE 0.136 0.140 0.136 0.140 Brain TSE 0.081 0.087 0.081 0.087 Knee TSE 0.292 0.293 0.292 0.293 Table 5.2: Computation Time (secs.) with 3232 calibration data. Dataset ESPIRiT ESPIRiT Nullspace Nullspace +PISCO +PISCO Brain MPRAGE 44.3 1.0 NA 1.0 Brain TSE 99.2 1.1 NA 1.1 Knee TSE 66.2 0.3 NA 0.3 To evaluate the speed of dierent algorithmic approaches, we measured the com- putation time employed by both approaches using a MATLAB in-house implemen- tation on a computer with an Intel Xeon E5-1603 2.8GHz quad core CPU and 32GB of RAM. Table 5.2 shows the computation time for all three datasets for ES- PIRiT, the nullspace-based approach, and the combination of these two algorithms with the PISCO techniques. We have not included the time corresponding to the nullspace-based approach since its implementation makes use of the procedure to directly calculate the matrices G(x) described in Sect. 5.3.5. Therefore, comparing its computation time against ESPIRiT would be an unfair comparison given that this type of techniques is not used in the original ESPIRiT description. However, as we will discuss later in the paper (and as is already evident from Table 5.2), there is no computational speed dierence between ESPIRiT and the nullspace-based method when these are combined with all the PISCO techniques. We also evaluated the memory usage for the nullspace-based algorithm and ES- PIRiT, by measuring the RAM memory associated to the most expensive variables 140 Table 5.3: Memory Usage (GB) with 3232 calibration data. Dataset ESPIRiT ESPIRiT Nullspace Nullspace +PISCO +PISCO Brain MPRAGE 1.01 0.55 25.33 0.55 Brain TSE 1.71 1.44 68.20 1.44 Knee TSE 1.91 0.43 18.89 0.43 needed in each algorithm. These correspond to the spatial-domain multichannel lters needed to naively calculate the matrices G(x) and B(x) in the nullspace- based algorithm and ESPIRiT, respectively 6 . We also added the RAM memory associated to the convolutional matrix C, which is also considerable. These results are reported in Table 5.3 from where it can be observed how massive the memory requirements are for the nullspace-based algorithm in comparison to ESPIRiT (if a naive implementation is used). However, we can also observe that the memory as- sociated to the considered variables is reduced drastically for both algorithms when combined with the PISCO techniques, making them similar in terms of memory usage. The detailed explanation of these memory savings will be given after our description of the rest of the PISCO techniques in the next section. Overall, these results demonstrate strong theoretical, empirical, and computa- tional similarities between ESPIRiT and the nullspace-based algorithm. Of course, another contribution of this paper is to introduce PISCO, which will enable sub- stantial improvements in computational complexity for the nullspace method, ES- PIRiT, and a number of related techniques. The remaining details of PISCO and a numerical evaluation of its performance will be presented in the next section. 6 It should be noted that in this assessment we are not considering the memory savings due to the procedure described in Sect. 5.3.5. This type of memory savings are only considered in the reported PISCO results. 141 5.4 Proposed Computational Methods Now we describe the remaining of our computational contributions included in PISCO, noticing that one of them has already been described in Sect. 5.3.5. As we explain next, each of these computational methods provides important improve- ments either in memory usage and/or computational complexity when incorporated into the naive algorithm described in Section 5.3.4. Notably, these eciency im- provements do not compromise estimation quality and, given the structural sim- ilarities between ESPIRiT and the nullspace-based method, it is also straightfor- ward to incorporate PISCO into ESPIRiT. In the following, besides providing a theoretical description of the computational methods in PISCO, we empirically as- sess eciency improvements when including these methods individually into the nullspace-based algorithm. In all the cases we use the Brain TSE dataset and con- sider dierent sizes for the calibration data as in Fig. 5.5. It should be mentioned that one of the PISCO techniques was introduced in Sect. 5.3.5, which makes the nullspace-based algorithm feasible in terms of memory usage. In the following ex- periments, the reported results regarding accuracy and computation time of our baseline nullspace-based method were obtained using the procedure described in Sect. 5.3.5 for the calculation of the G(x) matrices. However, unless it is specied, the reported results regarding memory usage were not necessarily obtained using the procedure described in Sect. 5.3.5. 142 5.4.1 Ellipsoidal versus Rectangular Convolution Kernels In previous work we have shown that considering an ellipsoidal support for the FIR lters in the construction of C oers signicant eciency improvements in comparison to using a rectangular support [81]. Specically, if we consider the ellipsoidal neighborhood e =fn2 Z 2 :knk ` 2 g in Eq. (5.7), instead of the rectangular neighborhood r =fn2Z 2 :knk 1 g, the size of C is considerably reduced given thatj e j<j r j. Then, memory usage and computation time should be reduced when using the ellipsoidal support instead of the rectangular support, given that C has smaller dimensions and therefore it is faster to construct. In Fig. 5.6 we assess the nullspace-based algorithm when using both types of support. In order to have a fair comparison we considered an ellipsoidal support e with = 3, and the rectangular support r with size (2 + 1) (2 + 1). We can see from Fig. 5.6(a) that both types of support obtained similar accuracy, however, there are considerable gains in memory usage when using the ellipsoidal support as we see in Fig. 5.6(b). In this case we are reporting the RAM memory associated to the spatial-domain multichannel lters needed to calculate the G(x) matrices in addition to the RAM memory associated to C. The computation time for both types of support was similar in this case, given that the dierence in size for both supports is not big enough. We should expect bigger dierences in memory usage and computational complexity for higher values of as shown in [81]. We have decided to only consider = 3 for the results in Fig. 5.6 and in the following sections, given that the estimation quality varied negligibly for higher values of . 143 24×24 32×32 64×64 128×128 256×256 0 0.1 0.2 0.3 0.4 0.5 SizeofCalibrationData NormalizedProjectionResidual Nullspace Nullspace+S 1 (a) 24×24 32×32 64×64 128×128 256×256 10 5 10 6 SizeofCalibrationData Megabytes(log) 1 (b) Figure 5.6: Assessment of the nullspace-based algorithm using lters with a rectangular (Nullspace) and an ellipsoidal support (Nullspace + S). 5.4.2 FFT-Based Calculation of Nullspace Vectors One important step in the nullspace-based method and ESPIRiT is the calculation of nullspace and rowspace vectors of C, respectively. As indicated in the Step (2) of the nullspace-based algorithm in Section 5.3.4, it is convenient to work with C H C instead of C, since C H C has the same nullspace vectors and usually much smaller dimensions. Even though the small dimensions of C H C permit a fast calculation of the SVD, we still need to calculate C in prior, which is a demanding task in terms of memory usage and computation time. Following previous work [98, 133], we propose to calculate C H C directly without calculating C rst. In [98] the authors proposed an FFT-based approach to calculate C H C directly in the single- channel case (i.e., Q = 1), which was extended to the multichannel case in [133]. However, the multichannel approach proposed in [133] can be demanding in terms of computation time. In the following we make the observation that it is possible 144 to rapidly calculate C H C in the multichannel case by leveraging on its structure and signal processing techniques. In the proposed approach the idea is to calculate each column of C H C separately. Let e i 2C Qjj be the ith canonical vector. Then, if we write e i using a stack of Q subvectorsfe (i) q g Q q=1 of sizejj, we have that all but theq(i),bi=jjc+1 subvector has a nonzero components. Then, the ith column of C H C would be given by C H Ce i = 2 6 6 6 6 6 4 C H 1 . . . C H Q 3 7 7 7 7 7 5 C q(i) e (i) q(i) ; (5.24) which reduces to calculatingQ products of the form C H p C q f wherep;q2f1;:::;Qg and f2C jj is a canonical vector. Given the convolutional structure of the matrices fC q g Q q=1 , it is possible to approximate these products eciently using the FFT [62, 98]. It can be theoretically shown that this procedure to calculate approximations of the columns of C H p C q is equivalent to performing a sliding-window that extracts a series of patches of sizejj from F 1 (F(Z(~ p ))F(Z(~ q ))); where each patch is centered at a dierent frequency component in . In the last expression ~ p ,F( p );Z(~ p ) correspond to a zero-padded version of ~ p ; and denotes the element-wise product. It should be noted that this sliding-window approach has been previously proposed in [98] for the single-channel. To the best 145 of our knowledge this is the rst time that this approach has been extended to the multichannel case. In Fig. 5.7 we show how the performance and memory usage of the nullspace- based algorithm changes when calculating C H C directly instead of calculating C rst. In this case we only report the RAM memory associated to C and C H C given that all the other variables remain the same. We observe that there are dramatic improvements in memory usage for big calibration data sizes, which correspond to the cases when the dimensions of C H C are considerably smaller than the dimensions of C. It should be noted that there are cases where the memory usage is worse for the proposed approach. These cases correspond to calibration data sizes where the dimensions of C H C are higher than the dimensions of C. However, in these cases the memory consumption by both matrices is not considerable. Notably, the improvements provided by the proposed approach negligibly aect estimation accuracy and computation time (not reported). 5.4.3 Smoothness-Based Interpolation In the nullspace-based algorithm proposed in Sect. 5.3.4 Step (4) corresponds to a location-wise estimation where the vector c(x) is estimated for every location x using the eigenvalue decomposition of G(x). This can be expensive in terms of computation time given that usually the number of locations is big. In this work we make the observation that it is possible to reduce the number of locations where the previously described process is performed, and estimate c(x) for the rest of locations 146 24×24 32×32 64×64 128×128 256×256 0 0.1 0.2 0.3 0.4 0.5 SizeofCalibrationData NormalizedProjectionResidual Nullspace Nullspace+CO 1 (a) 24×24 32×32 64×64 128×128 256×256 10 1 10 2 10 3 10 4 SizeofCalibrationData Megabytes(log) 1 (b) Figure 5.7: Assessment of the nullspace-based algorithm calculating the nullspace vectors in Step (2) with (Nullspace + CO) and without (Nullspace) the FFT-based approach. using a fast interpolation procedure. In fact, sensitivity maps arise from Maxwell's equations and therefore they possess smooth spatial characteristics. Based on this principle, we propose to estimate c(x) for a low-resolution spatial grid using the eigenvalue decomposition of G(x), and then use an FFT-based interpolation pro- cedure to obtain c(x) for the rest of locations considered at the desired spatial resolution. However, for this approach to be successful, it is necessary to obtain smooth sensitivity map estimates for the low-resolution spatial grid. This is not usually the case when conventional computational methods are used to calculate the eigenvalue decomposition of G(x), since the phase of the resulting sensitiv- ity map estimates commonly exhibits a random structure. To solve this issue, we reference the phase of the low-resolution spatial grid sensitivity map estimates to the phase of a low-resolution image obtained from the calibration data. To obtain 147 this low-resolution image we propose to use the low-resolution spatial grid sensi- tivity map estimates (with random phase) and a ltered version of the calibration data to obtain a SENSE-based reconstructed image [103]. We have observed that ltering the calibration data with a Gaussian window produces a smooth SENSE- reconstructed image with adequate phase characteristics to be used as reference. However, the Gaussian window might not be necessarily optimal, and other win- dows could achieve even better results. In Fig. 5.8 we show how the performance and computation time of the nullspace-based algorithm change when using the pro- posed FFT-based interpolation approach in Step (4). We can see that there are negligible variations in accuracy and dramatic improvements in computation time. 24×24 32×32 64×64 128×128 256×256 0 0.1 0.2 0.3 0.4 0.5 SizeofCalibrationData NormalizedProjectionResidual Nullspace Nullspace+I 1 (a) 24×24 32×32 64×64 128×128 256×256 0 20 40 60 80 100 120 SizeofCalibrationData Time(sec) 1 (b) Figure 5.8: Assessment of the nullspace-based algorithm using (Nullspace + I) and not using (Nullspace) the FFT-based interpolation approach in Step (4). 148 5.4.4 Power Iteration One of the most expensive steps in the nullspace-algorithm given in Sec. 5.3.4 corresponds to Step (4), where c(x) is estimated by calculating the eigenvector corresponding to the smallest eigenvalue of G(x). One possible approach to do this is to calculate the SVD of each G(x) for every location x, and then extract the last eigenvector of each matrix. However, this process can be computationally expensive when the number of locations is big. We make the observation that, given that the nullspace of G(x) tends to be unidimensional, it would be enough to only calculate one vector in this subspace instead of the whole SVD. For this purpose, we propose to use an iterative algorithm based on Power Iteration [26]. Remarkably, the steps of this algorithm can be easily implemented in parallel for all the locations, which highly improves the computation time. Before introducing the proposed algorithm, we review how to use Power Iteration to estimate c(x) for a particular location x, and then we show how we can perform the same procedure for multiple locations at the same time. The power iteration procedure for a particular location is summarized in Algorithm 1, where corresponds to an upper bound for Algorithm 1 Power Iteration procedure to estimate c(x) for a particular location x Input: G(x) Output: c(x) Initialization: Initialize c(x) randomly for k 1 to K do c(x) (I Q G(x))c(x) c(x) c(x)=kc(x)k 2 end for the magnitude of the eigenvalues of G(x). Even though this process is faster than 149 SVD, repeating it for all the considered locations can still be slow. We propose a computational method to execute the steps in Algorithm 1 for all the locations at the same time. For this purpose we rst dene the vectorsfg pq g Q p;q=1 , where g pq is a vector which contains the (p;q) entry of every matrix G(x). Then, we dene the vectorsf~ g pq g Q p;q=1 such that ~ g pq = 8 > > < > > : max 1 g pq ; if p =q g pq ; if p6=q where max corresponds to an upper bound for the eigenvalues of all the matrices G(x) and 1 is a vector with all entries equal to one. Given the structure of G(x), a simple choice for this bound corresponds to max =jj. Then, we also dene the vectorsfc q g Q q=1 where c q correspond to the vectorized version of c q (x) considering every location x. It should be noted that estimating c(x) for every location, as indicated in Step (4) of the nullspace-based algorithm, is equivalent to estimating the vectorsfc q g Q q=1 . Using these denitions, we propose Algorithm 2, which is equivalent to execute the steps in Algorithm 1 for all the locations in an all-at-once fashion. Notably, Algorithm 2 can be eciently implemented using programming languages like MATLAB or Python, where vectorization techniques can be used to rapidly calculate summations and reduce the number of for-loops. In our in-house implementation we only consider the outer for-loop associated to the variablek. The inner for-loop and summations are all performed using vectorization techniques. 150 Algorithm 2 All-at-once PowerIteration-based Algorithm Input:f~ g pq g Q p;q=1 Output: Estimated sensitivity map vectorsfc q g Q q=1 Initialization: Initializefc q g Q q=1 randomly for k 1 to K do c p P Q q=1 ~ g pq c q for all p2f1:::;Qg Compute the vector with entries [] m = 1= v u u t Q X p=1 j[c p ] m j 2 c p c p for all p2f1:::;Qg end for A fast convergence of the previous algorithms depends on having a big dierence in magnitude between the last two eigenvalues of the matrices G(x). In practice we have observed that this is usually the case (e.g., see Fig. 5.3), therefore, only a few iterations are needed for convergence. However, there are cases where the last two eigenvalues are close to each other, which corresponds to situations where more than one sensitivity map is required to describe a specic location. It should also be mentioned that there are other more advanced algorithms besides Power Iteration which can also be used to compute only one eigenvector of the matrices G(x), such as the Lanczos algorithm, which theoretically has better convergence that Power Iteration. However, we have decided to use Power Iteration since its steps can be easily implemented in all-at-once fashion as we previously described, while algorithms like Lanczos include normalization steps which are more involved to implement for multiple matrices at the same time. In Fig. 5.9 we show how the performance and computation time of the nullspace-based algorithm change when 151 using the proposed PowerIteration-based approach instead of a location-wise SVD decomposition approach in Step (4). We can see that there are negligible variations in accuracy and signicant improvements in computation time. For each case we selected K = 10 as the number of iterations in Algorithm 2. The same number of iterations is used in the experiments in the following subsection where all the PISCO techniques are used at the same time. 24×24 32×32 64×64 128×128 256×256 0 0.1 0.2 0.3 0.4 0.5 SizeofCalibrationData NormalizedProjectionResidual Nullspace Nullspace+P 1 (a) 24×24 32×32 64×64 128×128 256×256 0 20 40 60 80 100 120 SizeofCalibrationData Time(sec) 1 (b) Figure 5.9: Comparison between using the location-wise SVD approach (Nullspace) and the all-at-once PowerIteration-based approach (Nullspace + P) in Step (4) of the nullspace-based algorithm. 152 24×24 32×32 64×64 128×128 256×256 0 0.1 0.2 0.3 0.4 0.5 SizeofCalibrationData NormalizedProjectionResidual Nullspace Nullspace+PISCO 1 (a) 24×24 32×32 64×64 128×128 256×256 0 20 40 60 80 100 120 SizeofCalibrationData Time(sec) 1 (b) 24×24 32×32 64×64 128×128 256×256 10 3 10 4 10 5 SizeofCalibrationData Megabytes(log) 1 (c) Figure 5.10: Comparison between the nullspace-based algorithm with (Nullspace + PISCO) and without (Nullspace) the PISCO techniques. 153 5.4.5 Combined Approach Now that we have introduced each of the PISCO computational techniques, we empirically show the eciency improvements that are achieved when they are all integrated at the same time to the nullspace-based estimation algorithm. In Fig. 5.10 we compare the nullspace-based algorithm with its enhanced version after in- tegrating all the PISCO techniques. Remarkably, it can be observed that massive improvements in computation time and memory usage are achieved while maintain- ing estimation accuracy. For the memory usage analysis of the naive nullspace-based algorithm we have considered the RAM memory associated to the spatial-domain multichannel lters needed to calculate the G(x) matrices, in addition to the RAM memory associated to C. As expected, this memory usage is huge for each of the considered calibration data sizes. Notably, the memory usage is reduced drastically when using PISCO. In this case, we are considering the RAM memory associated to the G(x) matrices (which are calculated eciently using the method in Sect. 5.3.5), in addition to the RAM memory associated to C H C (which is calculated eciently using the method in Sect. 5.4.2). In Fig. 5.11 we show the estimated sensitivity maps for four representative channels of the three considered datasets when using the nullspace-based algorithm and its enhanced version using PISCO. For both algorithms we used a calibration data size equal to 32 32. For the three datasets it can be observed that there are strong similarities between the sensitivity maps estimated with each of the methods. In addition, we report estimation accu- racy, computation times, and memory usage in Tables I, II, and III, respectively. 154 For the sake of completeness, we also report the results of integrating the PISCO techniques into ESPIRiT, which is straightforward given the strong similarities in structure between the nullspace-based algorithm and ESPIRiT. From these tables it can be observed that the nullspace-based algorithm and ESPIRiT provide simi- lar results when including PISCO, which is expected given the theoretical analysis provided in Sect. 5.3.6. We obtained that in all the cases PISCO oers massive improvements in computation time and memory usage without compromising esti- mation accuracy. 155 ” c 1 c 2 c 3 c 4 Nullspace Nullspace + PISCO Nullspace Nullspace + PISCO Nullspace Nullspace + PISCO Figure 5.11: Estimated sensitivity maps using the nullspace method with and without PISCO for the three considered datasets. We show the magnitude of the sensitivity maps of four representative coils. 156 5.5 Discussion and Conclusions In this work we have proposed a novel theoretical framework for the subspace- based sensitivity map estimation problem. In comparison to the existent liter- ature, our results oer an alternative and complementary theoretical description that in our opinion provides a more intuitive mathematical explanation. Using this framework, we have proposed a sensitivity map estimation method that relies on nullspace and structured low-rank matrix modeling theory. We showed that the proposed nullspace-based method possesses strong theoretical and structural simi- larities with ESPIRiT. In the second part of this work we proposed a set of power- ful computational methods that allow an ecient implementation of the proposed nullspace-based method, and that can be straightforwardly integrated into existent methods like ESPIRiT. Remarkably, we showed that our proposed computational techniques produce massive improvements in memory usage and computation time in comparison to naive implementations of subspace-based methods, while main- taining estimation accuracy. As an example, roughly a 100-fold acceleration can be achieved when integrating our techniques to ESPIRiT. For the results shown in this paper we used an in-house MATLAB implementa- tion for all the considered methods. It should be mentioned that even bigger gains in performance could be obtained with more ecient implementations and other programming languages like C. As an example, the BART software package [119] provides a C-based implementation of ESPIRiT which is expected to be much faster 157 than MATLAB implementations. Notably, using our PISCO-based MATLAB im- plementation we were able to estimate sensitivity maps much faster than the C- based implementation of ESPIRiT provided in BART. For instance, using the Brain TSE dataset and a calibration data of size 24 24, our PISCO-based MATLAB implementation completed in 1 second while the BART C-based implementation completed in 19 seconds. It should be mentioned that both implementations ob- tained similar estimation accuracy. Therefore, substantial eciency improvements should be expected if the computational methods of this paper are implemented using more powerful programming languages. Moreover, the automatic selection of the dierent parameters in the proposed computational methods should also con- tribute to improving the obtained results. For instance, approaches based on the SURE estimator, as the one proposed in [52], could be combined with the proposed methodology. Finally, it should be highlighted that the computational methods developed in this work can also be applied beyond the sensitivity map estimation context. For instance, many accelerated MRI reconstruction methods are based on the construc- tion of convolutional structured matrices as the ones considered in this work (e.g., SPIRiT [86], PRUNO [132], SAKE [110], LORAKS [35], ALOHA [56]). Therefore, these methods could be beneted from the the proposed computational method described in Sect. 5.4.2, which eciently calculates this type of matrices using FFT-based convolutions. 158 Chapter 6 Conclusion During the last ve decades MRI has been constantly evolving in terms of per- formance and eciency, providing researchers and clinicians with an outstanding imaging technology capable of generating high-quality images of the living tissue in a safe and noninvasive manner. Even though MRI is an omnipresent resource in modern medicine, there are still several challenges that have limited MRI from reaching its full potential. Particularly, the time needed for the data-acquisition process can be restrictively long, which aects scanning throughput, patient com- fort, and accessibility to MRI technology. To address this problem, widely used approaches have consisted of accelerating the data-acquisition process by consider- ing advanced hardware and/or the acquisition of only a fraction of the data. In the case of partial acquisition, the accelerated data need to be eciently reconstructed before the generation of the corresponding images. In this work we made signif- icant contributions in this area. Firstly, we provided powerful algorithms able to successfully reconstruct accelerated MRI data in challenging acquisition scenarios. Secondly, we provided computational techniques which allow substantial eciency 159 improvements for reconstruction algorithms without compromising image-quality performance. A more detailed description of our main contributions is given in the following: 1. We have theoretically and empirically studied reconstruction algorithms based on structured low-rank modeling, which are capable of leveraging on con- straints learned from reference datasets in order to improve reconstruction performance. These algorithms account for the missing data samples by auto- matically identifying linear predictability relationships embedded in the data, and complement the reconstruction task with in prior knowledge learned from correlated datasets. We empirically showed that the use of reference data con- straints in the Fourier domain oers performance advantages over constraints in the spatial domain. We showed that, in the case of accelerated data acquired using an uniform undersampling scheme, the use of a nonconvex cost-function leads to better reconstruction performance than the use of a convex cost- function. Notably, we show that these developments allow the reconstruction of highly accelerate data in important MRI applications such as Echo Planar Imaging acquisition. 2. We have developed a robust reconstruction method based on structured low- rank modeling which is able to account for imperfections in the reference data. Additionally, this reconstruction method is able to learn linear predictabil- ity relationships between the reference data and the accelerated data, which produces important improvements in reconstruction performance. We have 160 empirically shown that this method is able to successfully reconstruct data in high acceleration regimes. 3. By making a systematic comparison we showed that the eciency of MRI reconstruction methods based on linear predictability principles can be sub- stantially improved by a simple modication in the underlying model used for prediction. When predicting missing data samples, we showed that using an el- lipsoidal neighborhood of samples instead of using a rectangular neighborhood of samples has a negligible eect on reconstruction performance. However, the computational time and memory usage are substantially improved. We showed this eect empirically by implementing both neighborhood shapes in several conventional and state-of-the-art reconstruction methods. 4. We provided a novel theoretical description of the parallel imaging MRI sen- sitivity map estimation problem based on linear predictability and structured low-rank modeling theory. These theoretical results provide an intuitive math- ematical perspective to study subspace-based sensitivity map estimation meth- ods. In addition, by relying on advanced signal processing techniques, we de- veloped and tested powerful computational techniques which allow massive eciency improvements when integrated to conventional sensitivity map es- timation methods. For instance, we showed that approximately a 100-fold acceleration in computational time can be achieved without compromising 161 estimation accuracy. Therefore, the overall reconstruction time of MRI re- construction methods that estimate sensitivity maps in prior should also be substantially reduced. Even though MRI can be considered as one the most remarkable technologies in medicine, there are still many areas where MRI can be improved, as the time needed for the data-acquisition. In this work we have developed theoretical and compu- tational tools that have a direct impact on reducing the acquisition time, and we have empirically shown how they allow substantial improvements in performance and eciency. However, there are still exciting research lines that we have not explored in this work. For example, in the rst part of this work, we have focused the analysis of our proposed reconstruction methods on performance by following a proof-of-principles approach, and we have not necessarily provided the most e- cient computational implementation. Therefore, the exploration of a more ecient implementation appears as a compelling research line with the potential of allowing the translation of our developed technology to clinical settings. Consequently, given that our proposed methods oer substantial advantages over conventional methods, we expect that their use in clinical applications will open new exciting research op- portunities for the MRI community. There are also promising research lines related to the methods developed in the second part of this work. We have focused the assessment of the proposed computational techniques in specic contexts related to MRI reconstruction, however, these techniques can be straightforwardly translated to other important medical imaging applications such as computed tomography. 162 Moreover, the proposed computational techniques can also be used to improve the eciency of computational imaging applications in areas beyond medical imaging. 163 Bibliography [1] Mehmet Ak cakaya, Steen Moeller, Sebastian Weing artner, and K^ amil U gurbil. Scan-specic robust articial-neural-networks for k-space interpo- lation (RAKI) reconstruction: Database-free deep learning for fast imaging. Magn. Reson. Med., 81(1):439{453, 2019. [2] Michael J Allison, Sathish Ramani, and Jerey A Fessler. Accelerated regular- ized estimation of MR coil sensitivities using augmented Lagrangian methods. IEEE Trans. Med. Imag., 32(3):556{564, 2012. [3] J. A. Bankson and S. M. Wright. Generalized partially parallel imaging with spatial lters. In Proc. Int. Soc. Magn. Reson. Med., page 1795, 2001. [4] Corey A Baron and Christian Beaulieu. Motion robust GRAPPA for echo- planar imaging. Magn. Reson. Med., 75(3):1166{1174, 2016. [5] M. A. Bernstein, K. F. King, and X. J. Zhou. Handbook of MRI Pulse Se- quences. Elsevier Academic Press, Burlington, 2004. [6] Matt A Bernstein, Sean B Fain, and Stephen J Riederer. Eect of windowing and zero-lled reconstruction of MRI data on spatial resolution and acquisi- tion strategy. Magn. Reson. Med., 14(3):270{280, 2001. [7] Chitresh Bhushan, Justin P Haldar, Soyoung Choi, Anand A Joshi, David W Shattuck, and Richard M Leahy. Co-registration and distortion correction of diusion and anatomical images based on inverse contrast normalization. Neuroimage, 115:269{280, 2015. [8] B. Bilgic, C. Liao, M. K. Manhard, Q. Tian, I. Chatnuntawech, S. S. Iyer, S. .F Cauley, T. Feiweier, S. Giri, andS. Y. Huang Y. Hu, J. R. Polimeni, L. L. Wald, and K. Setsompop. Robust high-quality multi-shot EPI with low-rank prior and machine learning. In Proc. Int. Soc. Magn. Reson. Med., page 1250, 2019. [9] Berkin Bilgic, Tae Hyung Kim, Congyu Liao, Mary Kate Manhard, Lawrence L Wald, Justin P Haldar, and Kawin Setsompop. Improving paral- lel imaging by jointly reconstructing multi-contrast data. Magn. Reson. Med., 80(2):619{632, 2018. 164 [10] Martin Blaimer, Felix Breuer, Matthias Mueller, Robin M Heidemann, Mark A Griswold, and Peter M Jakob. SMASH, SENSE, PILS, GRAPPA: how to choose the optimal method. Topics in Magnetic Resonance Imaging, 15(4):223{236, 2004. [11] R. N. Bracewell. The Fourier Transform and its Applications. McGraw-Hill, Boston, 2000. [12] H. Bruder, H. Fischer, H.-E. Reinfelder, and F. Schmitt. Image reconstruction for echo planar imaging with nonequidistant k-space sampling. Magn. Reson. Med., 23:311{323, 1992. [13] M. Buehrer, P. Boesiger, and S. Kozerke. Virtual body coil calibration for phased-array imaging. In Proc. Int. Soc. Magn. Reson. Med., page 760, 2009. [14] M. Buehrer, K. P. Pruessmann, P. Boesiger, and S. Kozerke. Array com- pression for MRI with large coil arrays. Magn. Reson. Med., 57:1131{1139, 2007. [15] Michael H. Buonocore and Lisheng Gao. Ghost artifact reduction for echo planar imaging using image phase correction. Magn. Reson. Med., 38:89{100, 1997. [16] M. Bydder, D. J. Larkman, and J. V. Hajnal. Combination of signals from array coils using image-based estimation of coil sensitivity proles. Magn. Reson. Med., 47:539{548, 2002. [17] E. J. Cand es, C. A. Sing-Long, and J. D. Trzasko. Unbiased risk estimates for singular value thresholding and spectral estimators. IEEE Trans. Signal Process., 61:4643{4657, 2013. [18] J. W. Carlson. An algorithm for NMR imaging reconstruction based on mul- tiple RF receiver coils. J. Magn. Reson., 74:376{380, 1987. [19] H.-C. Chang and N.-K. Chen. Joint correction of Nyquist artifact and mi- nuscule motion-induced aliasing artifact in interleaved diusion weighted EPI data using a composite two-dimensional phase correction procedure. Magn. Reson. Imag., 34:974{979, 2016. [20] Y. Chang, D. Liang, and L. Ying. Nonlinear GRAPPA: A kernel approach to parallel MRI reconstruction. Magn. Reson. Med., 68:730{740, 2011. [21] N.-k. Chen and Alice M. Wyrwicz. Removal of EPI Nyquist ghost artifacts with two-dimensional phase correction. Magn. Reson. Med., 51:1247{1253, 2004. [22] Nan-Kuei Chen, Alexandru V. Avram, and Allen W. Song. Two-dimensional phase cycled reconstruction for inherent correction of echo-planar imaging Nyquist artifacts. Magn. Reson. Med., 66:1057{1066, 2011. 165 [23] K. F. Cheung and R. J. Marks II. Imaging sampling below the Nyquist density without aliasing. J. Opt. Soc. Am. A, 7:92{105, 1990. [24] L. Cordero-Grande, D. Christiaens, J. Hutter, A. N. Price, and J. V. Haj- nal. Complex diusion-weighted image estimation via matrix recovery under general noise models. NeuroImage, 200:391{404, 2019. [25] Anagha Deshmane, Vikas Gulani, Mark A. Griswold, and Nicole Seiberlich. Parallel mr imaging. J. Magn. Reson., 36(1):55{72, 2012. [26] G. Golub and C. van Loan. Matrix Computations. The Johns Hopkins Uni- versity Press, London, third edition, 1996. [27] Stuart M. Grieve, Andrew M. Blamire, and Peter Styles. Elimination of Nyquist ghosting caused by read-out to phase-encode gradient cross-terms in EPI. Magn. Reson. Med., 47:337{343, 2002. [28] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, and A. Haase. Generalized autocalibrating partially parallel acqui- sitions (GRAPPA). Magn. Reson. Med., 47:1202{1210, 2002. [29] Derya Gol Gungor, Rizwan Ahmad, and Lee C Potter. A subspace approach to blind coil sensitivity estimation in parallel MRI. Journal of Cardiovascular Magnetic Resonance, 16(1):1{2, 2014. [30] M. _ I G urelli and C. L. Nikias. EVAM: An eigenvector-based algorithm for multichannel blind deconvolution of input colored signals. IEEE Trans. Signal Process., 43:134{149, 1995. [31] J. P. Haldar. Low-rank modeling of local k-space neighborhoods (LORAKS): Implementation and examples for reproducible research. Technical Report USC-SIPI-414, University of Southern California, Los Angeles, CA, April 2014. http://sipi.usc.edu/reports/abstracts.php?rid=sipi-414. [32] J. P. Haldar. Low-rank modeling of local k-space neighborhoods: From phase and support constraints to structured sparsity. In Wavelets and Sparsity XVI, Proc. SPIE 9597, page 959710, 2015. [33] J. P. Haldar and T. H. Kim. Computational imaging with LORAKS: Recon- structing linearly predictable signals using low-rank matrix regularization. In Proc. Asilomar Conf. Sig. Sys. Comp., 2017. [34] Justin P. Haldar. Calibrationless partial Fourier reconstruction of MR images with slowly-varying phase: A rank-decient matrix recovery approach. In ISMRM Workshop on Data Sampling & Image Reconstruction, Sedona, 2013. [35] Justin P Haldar. Low-Rank Modeling of Local-Space Neighborhoods (LO- RAKS) for Constrained MRI. IEEE Trans. Med. Imag., 33:668{681, 2014. 166 [36] Justin P Haldar. Autocalibrated LORAKS for fast constrained MRI recon- struction. In Proc. IEEE Int. Symp. Biomed. Imag., pages 910{913, 2015. [37] Justin P Haldar and Kawin Setsompop. Linear predictability in magnetic res- onance imaging reconstruction: Leveraging shift-invariant Fourier structure for faster and better imaging. IEEE Signal Process. Mag., 37:69{82, 2020. [38] Justin P. Haldar and Jingwei Zhuo. P-LORAKS: Low-rank modeling of lo- cal k-space neighborhoods with parallel imaging data. Magn. Reson. Med., 75:1499{1514, 2015. [39] G. Harikumar and Y. Bresler. Perfect blind restoration of images blurred by multiple lters: theory and ecient algorithms. IEEE Trans. Image Process., 8:202{219, 1999. [40] W. S. Hoge and J. R. Polimeni. Artifact correction in accelerated-segmented EPI data via dual-polarity GRAPPA. In Proc. Int. Soc. Magn. Reson. Med., page 449, 2017. [41] W. S. Hoge, F. R. Preiswerk, J. R. Polimeni, S. Yengul, P. A. Ciris, and B. Madore. Ultrasound monitoring of a respiratory phantom for the develop- ment and validation of segmented EPI reconstruction methods. In Proc. Int. Soc. Magn. Reson. Med., page 1307, 2017. [42] W Scott Hoge and Jonathan R Polimeni. Dual-polarity GRAPPA for simulta- neous reconstruction and ghost correction of echo planar imaging data. Magn. Reson. Med., 76:32{44, 2016. [43] W Scott Hoge, Kawin Setsompop, and Jonathan R Polimeni. Dual-polarity slice-GRAPPA for concurrent ghost correction and slice separation in simul- taneous multi-slice EPI. Magn. Reson. Med., 80(4):1364{1375, 2018. [44] W Scott Hoge, Huan Tan, and Robert A Kraft. Robust EPI Nyquist ghost elimination via spatial and temporal encoding. Magn. Reson. Med., 64:1781{ 1791, 2010. [45] H Christian M Holme, Sebastian Rosenzweig, Frank Ong, Robin N Wilke, Michael Lustig, and Martin Uecker. ENLIVE: An ecient nonlinear method for calibrationless and robust parallel imaging. Sci. Rep., 9(1):1{13, 2019. [46] Xiaoping Hu and Tuong Huu Le. Artifact reduction in EPI with phase- encoded reference scan. Magn. Reson. Med., 36:166{171, 1996. [47] F. Huang, W. Lin, and Y. Li. Partial Fourier reconstruction through data tting and convolution in k-space. Magn. Reson. Med., 62:1261{1269, 2009. [48] F. Huang, S. Vijayakumar, Y. Li, S. Hertel, and G. R. Duensing. A software channel compression technique for faster reconstruction with many channels. Magn. Reson. Imag., 26:133{141, 2008. 167 [49] M. Hutchinson and U. Ra. Fast MRI data acquisition using multiple detec- tors. Magn. Reson. Med., 6:87{91, 1988. [50] J. S. Hyde, A. Jesmanowicz, W. Froncisz, J. B. Kneeland, and T. M. Grist. Parallel image acquisition from noninteracting local coils. J. Magn. Reson., 70:512{517, 1986. [51] Julianna D. Ianni, E. Brian Welch, and William A. Grissom. Ghost reduc- tion in echo-planar imaging by joint reconstruction of images and line-to-line delays and phase errors. Magn. Reson. Med., 79:3114{3121, 2018. [52] Siddharth Iyer, Frank Ong, Kawin Setsompop, Mariya Doneva, and Michael Lustig. SURE-based automatic parameter selection for ESPIRiT calibration. Magn. Reson. Med., 84(6):3423{3437, 2020. [53] M. Jacob, M. P. Mani, and J. C. Ye. Structured low-rank algorithms: Theory, magnetic resonance applications, and links to machine learning. IEEE Signal Process. Mag., 37:54{68, 2020. [54] A. Javed and K. S. Nayak. Single-shot EPI for ASL-CMR. Magn. Reson. Med., 84:738{750, 2020. [55] W. Jiang, P. E. Z. Larson, and M. Lustig. Simultaneous estimation of auto- calibration data and gradient delays in non-Cartesian parallel MRI using low- rank constraints. In Proc. Int. Soc. Magn. Reson. Med., page 939, 2016. [56] K. H. Jin, D. Lee, and J. C. Ye. A general framework for compressed sens- ing and parallel MRI using annihilating lter based low-rank Hankel matrix. IEEE Trans. Comput. Imaging, 2:480{495, 2016. [57] Stephen L Keeling and Roland Bammer. A variational approach to magnetic resonance coil sensitivity estimation. Applied Mathematics and Computation, 158(2):359{388, 2004. [58] P. Kellman and E. R. McVeigh. Phased array ghost elimination. NMR Biomed., 19:352{361, 2006. [59] J. R. Kelton, R. L. Magin, and S. M. Wright. An algorithm for rapid image acquisition using multiple receiver coils. In Proc. Int. Soc. Magn. Reson. Med., page 1172, Amsterdam, 1989. [60] D. Kim, S. F. Cauley, K. S. Nayak, R. M. Leahy, and J. P. Haldar. Region- optimized virtual (ROVir) coils: Localization and/or suppression of spatial regions using sensor-domain beamforming. Magn. Reson. Med., 86:197{212, 2021. [61] T. H. Kim and J. P. Haldar. SMS-LORAKS: Calibrationless simultaneous multislice MRI using low-rank matrix modeling. In Proc. IEEE Int. Symp. Biomed. Imag., pages 323{326, 2015. 168 [62] T. H. Kim and J. P. Haldar. LORAKS software version 2.0: Faster implemen- tation and enhanced capabilities. Technical Report USC-SIPI-443, University of Southern California, Los Angeles, CA, May 2018. [63] Tae Hyung Kim, Berkin Bilgic, Daniel Polak, Kawin Setsompop, and Justin P. Haldar. Wave-LORAKS for faster wave-CAIPI MRI. In Proc. Int. Soc. Magn. Reson. Med., page 1037, 2017. [64] Tae Hyung Kim, Pratyush Garg, and Justin P. Haldar. LORAKI: Autocali- brated recurrent neural networks for autoregressive reconstruction in k-space. Preprint, 2019. arXiv:1904.09390. [65] Tae Hyung Kim and Justin P Haldar. The Fourier radial error spectrum plot: A more nuanced quantitative evaluation of image reconstruction quality. In Proc. IEEE Int. Symp. Biomed. Imag., pages 61{64, 2018. [66] Tae Hyung Kim, Kawin Setsompop, and Justin P. Haldar. LORAKS makes better SENSE: Phase-constrained partial Fourier SENSE reconstruction with- out phase calibration. Magn. Reson. Med., 77:1021{1035, 2016. [67] F. Knoll, K. Hammernik, C. Zhang, S. Moeller, T. Pock, D. K. Sodickson, and M. Akcakaya. Deep-learning methods for parallel magnetic resonance image reconstruction: A survey of the current approaches, trends, and issues. IEEE Signal Process. Mag., 37:128{140, 2020. [68] D. Kwiat, S. Einav, and G. Navon. A decoupled coil detector array for fast image acquisition in magnetic-resonance-imaging. Med. Phys., 18:251{265, 1991. [69] W. E. Kyriakos, L. P. Panych, D. F. Kacher, C.-F. Westin, S. M. Bao, R. V. Mulkern, and F. A. Jolesz. Sensitivity proles from an array of coils for encoding and reconstruction in parallel (SPACE-RIP). Magn. Reson. Med., 44:301{308, 2000. [70] Denis Le Bihan, Cyril Poupon, Alexis Amadon, and Franck Lethimonnier. Artifacts and pitfalls in diusion MRI. J. Magn. Reson. Imag., 24(3):478{ 488, 2006. [71] Juyoung Lee, Kyong Hwan Jin, and Jong Chul Ye. Reference-free EPI Nyquist ghost correction using annihilating lter-based low rank Hankel matrix for K- space interpolation. Magn. Reson. Med., 76:1775{1789, 2016. [72] D. Liang, J. Cheng, Z. Ke, and L. Ying. Deep magnetic resonance image reconstruction: Inverse problems meet neural networks. IEEE Signal Process. Mag., 37:141{151, 2020. [73] Z.-P. Liang. Constrained image reconstruction from incomplete and noisy data: A new parametric approach. PhD thesis, Case Western Reserve Uni- versity, Cleveland, OH, USA, 1989. 169 [74] Z.-P. Liang, F. Boada, T. Constable, E. M. Haacke, P. C. Lauterbur, and M. R. Smith. Constrained reconstruction methods in MR imaging. Rev. Magn. Reson. Med., 4:67{185, 1992. [75] Z.-P. Liang, E. M. Haacke, and C. W. Thomas. High-resolution inversion of nite Fourier transform data through a localised polynomial approximation. Inverse Probl., 5:831{847, 1989. [76] Congyu Liao, Mary Kate Manhard, Berkin Bilgic, Qiyuan Tian, Qiuyun Fan, Sohyun Han, Fuyixue Wang, Daniel Joseph Park, Thomas Witzel, Jianhui Zhong, et al. Phase-matched virtual coil reconstruction for highly accelerated diusion echo-planar imaging. NeuroImage, 194:291{302, 2019. [77] Y. Liu, M. Lyu, M. Barth, Z. Yi, A. T. Leong, F. Chen, Y. Feng, and E. X. Wu. PEC-GRAPPA reconstruction of simultaneous multislice EPI with slice- dependent 2D Nyquist ghost correction. Magn. Reson. Med., 81:1924{1934, 2019. [78] R. A. Lobos and J. P. Haldar. Improving the performance of accelerated image reconstruction in k-space: The importance of kernel shape. In Proc. Int. Soc. Magn. Reson. Med., page 2407, 2019. [79] R. A. Lobos, A. Javed, K. S. Nayak, W. S. Hoge, and J. P. Haldar. Robust autocalibrated LORAKS for improved EPI ghost correction with structured low-rank matrix models. In Proc. Int. Soc. Magn. Reson. Med., page 3533, 2018. [80] R. A. Lobos, T. H. Kim, W. S. Hoge, and J. P. Haldar. Navigator-free EPI ghost correction using low-rank matrix modeling: Theoretical insights and practical improvements. In Proc. Int. Soc. Magn. Reson. Med., page 0449, 2017. [81] Rodrigo A Lobos and Justin P Haldar. On the shape of convolution kernels in MRI reconstruction: Rectangles versus ellipsoids. Magn. Reson. Med., 87(6):2989{2996, 2022. [82] Rodrigo A Lobos, W Scott Hoge, Ahsan Javed, Congyu Liao, Kawin Setsom- pop, Krishna S Nayak, and Justin P Haldar. Robust autocalibrated structured low-rank EPI ghost correction. Magn. Reson. Med., 85(6):3403{3419, 2021. [83] Rodrigo A Lobos, Ahsan Javed, Krishna S Nayak, W Scott Hoge, and Justin P Haldar. Robust autocalibrated LORAKS for EPI ghost correction. In Proc. IEEE Int. Symp. Biomed. Imag., pages 663{666, 2018. [84] Rodrigo A Lobos, Tae Hyung Kim, W Scott Hoge, and Justin P Haldar. Navigator-free EPI ghost correction with structured low-rank matrix models: New theory and methods. IEEE Trans. Med. Imag., 37(11):2390{2402, 2018. 170 [85] D. G. Luenberger. Optimization by Vector Space Methods. John Wiley & Sons, New York, 1969. [86] M. Lustig and J. M. Pauly. SPIRiT: Iterative self-consistent parallel imaging reconstruction from arbitary k-space. Magn. Reson. Med., 65:457{471, 2010. [87] J. Lyu, U. Nakarmi, D. Liang, J. Sheng, and L. Ying. KerNL: Kernel-based nonlinear approach to parallel MRI reconstruction. IEEE Trans. Med. Imag., 38:312{321, 2019. [88] M. Lyu, M. Barth, V. B. Xie, Y. Liu, Y. Feng, and E. X. Wu. Robust 2D Nyquist ghost correction for simultaneous multislice (SMS) EPI using phase error correction SENSE and virtual coil SAKE. In Proc. Int. Soc. Magn. Reson. Med., page 515, 2017. [89] M. Lyu, M. Barth, V. B. Xie, Y. Liu, X. Ma, Y. Feng, and E. X. Wu. Robust SENSE reconstruction of simultaneous multislice EPI with low-rank enhanced coil sensitivity calibration and slice-dependent 2D Nyquist ghost correction. Magn. Reson. Med., 80:1376{1390, 2018. [90] M. Mani, M. Jacob, G. McKinnon, B. Yang, B. Rutt, A. Kerr, and V. Mag- notta. SMS-MUSSELS: A navigator-free reconstruction for slice-accelerated multi-shot diusion imaging. In Proc. Int. Soc. Magn. Reson. Med., page 0233, 2019. [91] M. Mani, V. Magnotta, D. Kelley, and M. Jacob. Comprehensive reconstruc- tion of multi-shot multi-channel diusion data using MUSSELS. In Proc. IEEE Eng. Med. Bio. Conf., pages 1107{1110, 2016. [92] M. Mani, S. Poddar, V. Magnotta, and M. Jacob. Trajectory correction of radial data using MUSSELS. In Proc. Int. Soc. Magn. Reson. Med., page 1197, 2017. [93] Merry Mani, Mathews Jacob, Douglas Kelley, and Vincent Magnotta. Multi- shot sensitivity-encoded diusion data recovery using structured low-rank ma- trix completion (MUSSELS). Magn. Reson. Med., 78:494{507, 2017. [94] J. V. Manj on, P. Coup e, L. Concha, A. Buades, D. L. Collins, and M. Robles. Diusion weighted image denoising using overcomplete local PCA. PloS one, 8:e73021, 2013. [95] Peter Manseld. Multi-planar image formation using NMR spin echoes. J. Phys. C: Solid State Phys., 10:L55, 1977. [96] Robert L Morrison, Mathews Jacob, and Minh N Do. Multichannel esti- mation of coil sensitivities in parallel MRI. In 2007 4th IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pages 117{120. IEEE, 2007. 171 [97] G. Ongie and M. Jacob. O-the-grid recovery of piecewise constant images from few Fourier samples. SIAM J. Imaging Sci., 9:1004{1041, 2016. [98] G. Ongie and M. Jacob. A fast algorithm for convolutional structured low- rank matrix recovery. IEEE Trans. Comput. Imaging, 3:535{550, 2017. [99] Xi Peng, Bradley P Sutton, Fan Lam, and Zhi-Pei Liang. DeepSENSE: Learn- ing coil sensitivity functions for SENSE reconstruction using deep learning. Magn. Reson. Med., 87(4):1894{1902, 2022. [100] Jonathan R Polimeni, Himanshu Bhat, Thomas Witzel, Thomas Benner, Thorsten Feiweier, Souheil J Inati, Ville Renvall, Keith Heberlein, and Lawrence L Wald. Reducing sensitivity losses due to respiration and motion in accelerated echo planar imaging by reordering the autocalibration data acquisition. Magn. Reson. Med., 75(2):665{679, 2016. [101] B. A. Poser and K. Setsompop. Pulse sequences and parallel imaging for high spatiotemporal resolution MRI at ultra-high eld. NeuroImage, 2017. Early View. [102] Klaas P. Pruessmann, Markus Weiger, Peter B ornert, and Peter Boesiger. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn. Reson. Med., 46:638{651, 2001. [103] Klaas P. Pruessmann, Markus Weiger, Markus B. Scheidegger, and Peter Boesiger. SENSE: Sensitivity encoding for fast MRI. Magn. Reson. Med., 42:952{962, 1999. [104] J. B. Ra and C. Y. Rim. Fast imaging using subencoding data sets from multiple detectors. Magn. Reson. Med., 30:142{145, 1993. [105] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM Rev., 52:471{ 501, 2010. [106] Scott B. Reeder, Ergin Atalar, Anthony Z. Faranesh, and Elliot R. McVeigh. Referenceless interleaved echo-planar imaging. Magn. Reson. Med., 41:87{94, 1999. [107] P. B. Roemer, W. A. Edelstein, C. E. Hayes, S. P. Souza, and O. M. Mueller. The NMR phased array. Magn. Reson. Med., 16:192{225, 1990. [108] H. She, R.-R. Chen, D. Liang, Y. Chang, and L. Ying. Image reconstruction from phased-array data based on multichannel blind deconvolution. Magn. Reson. Imag., 33:1106{1113, 2015. [109] D Sheltraw, B Inglis, V Deshpande, and M Trumpis. Simultaneous reduc- tion of two common autocalibration errors GRAPPA EPI time series data. Preprint, 2012. arXiv:1208.0972. 172 [110] P. J. Shin, P. E. Z. Larson, M. A. Ohliger, M. Elad, J. M. Pauly, D. B. Vi- gneron, and M. Lustig. Calibrationless parallel imaging reconstruction based on structured low-rank matrix completion. Magn. Reson. Med., 72:959{970, 2014. [111] S Skare, DB Clayton, R Newbould, M Moseley, and R Bammer. A fast and robust minimum entropy based non-interactive Nyquist ghost correction algorithm. In Proc. Int. Soc. Magn. Reson. Med., page 2349, 2006. [112] D. K. Sodickson and W. J. Manning. Simultaneous acquisition of spatial harmonics (SMASH): Fast imaging with radiofrequency coil arrays. Magn. Reson. Med., 38:591{603, 1997. [113] M. K. Stelhing, R. Turner, and P. Manseld. Echo-planar imaging: Magnetic resonance imaging in a fraction of a second. Science, 254:43{50, 1991. [114] S Lalith Talagala, Joelle E Sarlls, Siyuan Liu, and Souheil J Inati. Im- provement of temporal signal-to-noise ratio of GRAPPA accelerated echo planar imaging using a FLASH based calibration scan. Magn. Reson. Med., 75(6):2362{2371, 2016. [115] J. Trzasko and A. Manduca. Local versus global low-rank promotion in dy- namic MRI series reconstruction. In Proc. Int. Soc. Magn. Reson. Med., page 4371, 2011. [116] J. D. Trzasko and A. Manduca. CLEAR: Calibration-free parallel imaging using locally low-rank encouraging reconstruction. In Proc. Int. Soc. Magn. Reson. Med., page 517, 2012. [117] M. Uecker, T. Hohage, K. T. Block, and J.. Frahm. Image reconstruction by regularized nonlinear inversion { joint estimation of coil sensitivites and content. Magn. Reson. Med., 60:674{682, 2008. [118] M. Uecker, P. Lai, M. J. Murphy, P. Virtue, M. Elad, J. M. Pauly, S. S. Vasanawala, and M. Lustig. ESPIRiT { an eigenvalue approach to autocali- brating parallel MRI: Where SENSE meets GRAPPA. Magn. Reson. Med., 71:990{1001, 2014. [119] Martin Uecker, Frank Ong, Jonathan I Tamir, Dara Bahri, Patrick Virtue, Joseph Y Cheng, Tao Zhang, and Michael Lustig. Berkeley advanced recon- struction toolbox. In Proc. Int. Soc. Magn. Reson. Med., page 2486, 2015. [120] J. Veraart, D. S. Novikov, D. Christiaens, B. Ades-Aron, J. Sijbers, and E. Fieremans. Denoising of diusion MRI using random matrix theory. Neu- roImage, 142:394{406, 2016. [121] David O Walsh, Arthur F Gmitro, and Michael W Marcellin. Adaptive recon- struction of phased array MR imagery. Magn. Reson. Med., 43(5):682{690, 2000. 173 [122] J. Wang, B. Zhang, K. Zhong, and Y. Zhuo. Image domain based fast GRAPPA reconstruction and relative SNR degradation factor. In Proc. Int. Soc. Magn. Reson. Med., page 2428, 2005. [123] Qing-San Xiang and Frank Q Ye. Correction for geometric distortion and N/2 ghosting in EPI by phase labeling for additional coordinate encoding (PLACE). Magn. Reson. Med., 57:731{741, 2007. [124] Victor B. Xie, Mengye Lyu, Yilong Liu, Yanqiu Feng, and Ed X. Wu. Ro- bust EPI Nyquist ghost removal by incorporating phase error correction with sensitivity encoding (PEC-SENSE). Magn. Reson. Med., 79:943{951, 2018. [125] Dan Xu, Kevin F. King, Yuval Zur, and R. Scott Hinks. Robust 2D phase correction for echo planar imaging under a tight eld-of-view. Magn. Reson. Med., 64:1800{1813, 2010. [126] G. Xu, H. Liu, L. Tong, and T. Kailath. A least-squares approach to blind channel identication. IEEE Trans. Signal Process., 43:2982{2993, 1995. [127] Uten Yarach, Myung-Ho In, Itthi Chatnuntawech, Berkin Bilgic, Frank Go- denschweger, Hendrik Mattern, Alessandro Sciarra, and Oliver Speck. Model- based iterative reconstruction for single-shot EPI at 7T. Magn. Reson. Med., 78:2250{2264, 2017. [128] Ernest N Yeh, Charles A McKenzie, Michael A Ohliger, and Daniel K Sod- ickson. Parallel magnetic resonance imaging with adaptive radius in k-space (PARS): Constrained image reconstruction using k-space locality in radiofre- quency coil encoded data. Magn. Reson. Med., 53(6):1383{1392, 2005. [129] Leslie Ying and Zhi-Pei Liang. Parallel MRI using phased array coils. IEEE Signal Process. Mag., 27(4):90{98, 2010. [130] Leslie Ying and Jinhua Sheng. Joint image reconstruction and sensitivity estimation in SENSE (JSENSE). Magn. Reson. Med., 57:1196{1202, 2007. [131] Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, et al. fastmri: An open dataset and benchmarks for accelerated mri. arXiv preprint arXiv:1811.08839, 2018. [132] J. Zhang, C. Liu, and M. E. Moseley. Parallel reconstruction using null oper- ations. Magn. Reson. Med., 66:1241{1253, 2011. [133] Shen Zhao, Lee C Potter, and Rizwan Ahmad. High-dimensional fast con- volutional framework (HICU) for calibrationless MRI. Magn. Reson. Med., 86(3):1212{1225, 2021. 174

Abstract (if available)

Abstract Magnetic Resonance Imaging (MRI) is a safe, versatile, and noninvasive imaging modality that has revolutionized medicine by providing high-quality images of living tissue. Even though MRI technology has been constantly evolving since its origins in the 1970's, the time needed to acquire MRI data can still be restrictively long in real clinical applications. This has limited MRI from reaching its full potential, reason why constant efforts have been made in order to accelerate the data-acquisition process. Most of the adopted approaches to reduce the data-acquisition time involve the individual or combined use of advanced hardware and partial data-acquisition, however, conventional approaches based on these principles can be subject to negative effects in the signal-to-noise ratio, resolution, and/or the presence of undesired artifacts. In this study we propose new computational imaging methods which are able to successfully reconstruct partially acquired data obtained using advanced MRI hardware. Additionally, in this study we propose powerful computational techniques which allow dramatic improvements in computational efficiency when performing the reconstruction task.

In the first part of this study we propose reconstruction methods for partially acquired data which are able to automatically identify multiple linear predictability relationships present in the MRI data. For this purpose we rely on modern structured low-rank modeling theory and advanced optimization techniques. Based on these tools, we make the novel observation that reference data can be used to learn additional linear predictability relationships which considerably improve reconstruction performance in challenging scenarios. For instance, we show that linear relationships learned from reference data allow the reconstruction of highly accelerated data acquired using an undersampling scheme with a uniform structure. Notably, we show that these linear predictability relationships can be learned even in cases where the reference data are not pristine.

In the second part of this study we propose computational techniques in order to improve the efficiency of MRI reconstruction methods. We show that reconstruction methods which are based on the linearly predictable structure of the MRI data can be significantly enhanced in terms of computational complexity by a simple modification in the model used for prediction. Specifically, we show that a missing sample can be predicted using an ellipsoidal neighborhood of samples instead of a rectangular neighborhood of samples, which allows important improvements in computational time and memory usage with negligible effects on performance.

Finally, we study how to improve the efficiency of parallel imaging sensitivity map estimation methods. An accurate estimation of sensitivity maps is a fundamental piece in many modern parallel imaging MRI reconstruction methods, which should also be performed efficiently in order to reduce the overall reconstruction time. In this study we provide a novel theoretical description of the sensitivity map estimation problem by leveraging on the linearly predictable structure of the MRI data. Then, relying on these theoretical results, we propose a powerful set of computational techniques which allow massive improvements in computational complexity when integrated to sensitivity map estimation methods based on subspaces. We show that widely used estimation methods can achieve approximately 100-fold acceleration in computational time and dramatic savings in memory usage without sacrificing estimation accuracy. Remarkably, each of the proposed computational techniques can also be used individually to improve the efficiency of methods in other signal processing applications beyond MRI.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Shift-invariant autoregressive reconstruction for MRI

PDF

Optimization methods and algorithms for constrained magnetic resonance imaging

PDF

Improved brain dynamic contrast enhanced MRI using model-based reconstruction

PDF

Efficient and accurate 3D FISP-MRF at 0.55 Tesla

PDF

Methods for improving reliability and consistency in diffusion MRI analysis

PDF

Functional connectivity analysis and network identification in the human brain

PDF

Technology for improved 3D dynamic MRI

PDF

Model based PET image reconstruction and kinetic parameter estimation

PDF

Correction, coregistration and connectivity analysis of multi-contrast brain MRI

PDF

Applications of graph theory to brain connectivity analysis

PDF

High-dimensional magnetic resonance imaging of microstructure

PDF

Diffusion MRI of the human brain: signal modeling and quantitative analysis

PDF

Novel theoretical characterization and optimization of experimental efficiency for diffusion MRI

PDF

Improving sensitivity and spatial coverage of myocardial arterial spin labeling

PDF

Seeing sleep: real-time MRI methods for the evaluation of sleep apnea

PDF

Point-based representations for 3D perception and reconstruction

PDF

Direct wholebody Patlak and Logan image estimation from listmode PET data

PDF

Functional real-time MRI of the upper airway

PDF

Measuring functional connectivity of the brain

PDF

Autoregression and structured low-rank modeling of sinograms

Asset Metadata

Creator Lobos, Rodrigo Alejandro (author)

Core Title New theory and methods for accelerated MRI reconstruction

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Electrical Engineering

Degree Conferral Date 2022-12

Publication Date 06/14/2023

Defense Date 09/28/2022

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

Tag biomedical imaging,computational imaging,magnetic resonance imaging,OAI-PMH Harvest,signal processing

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Haldar, Justin (committee chair), Leahy, Richard (committee member), Wood, John (committee member)

Creator Email r.lobos.morales@gmail.com,rlobos@usc.edu

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

Unique identifier UC112620924

Identifier etd-LobosRodri-11376.pdf (filename)

Legacy Identifier etd-LobosRodri-11376

Document Type Dissertation

Format theses (aat)

Rights Lobos, Rodrigo Alejandro

Internet Media Type application/pdf

Type texts

Source 20221214-usctheses-batch-997 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu