Synopsis

Keywords: Image acquisition: Fast imaging, Image acquisition: Reconstruction

Magnetic Resonance Imaging (MRI) is a non-invasive medical imaging technique that has emerged as a pivotal clinical diagnostic tool over the last decades. Yet, its extended scanning times often compromise patient comfort and attainable image resolution. In this course, I will review standard acceleration techniques to shorten MRI scans: First, I will discuss parallel imaging methods based on multicoil acquisition, deterministic under-sampling in k-space and linear image reconstruction. Second, I will expose how Compressed Sensing, which relies on incoherent under-sampling, sparsity and non-linear image reconstruction, has been instantiated in MRI, notably to reach higher acceleration regimes.

Background

As the nuclear magnetic resonance (NMR) phenomenon is transient and of short lifespan, collecting MRI data is a sequential and lengthy procedure that requires multiple repetitions of the same sequence (RF pulse delivery through the transmit coil, spatial encoding using time-varying gradient waveforms and effective readout through the receptive coil using an analog-to-digital converter). Due to the spatial encoding performed along gradient waveforms, MRI data represent the Fourier samples of the spatial Fourier transform of the organ to be imaged. The sequential aspect of MRI data acquisition in the Fourier space (also called k-space) means that the sampling pattern is segmented in time over multiple shots, the latter being either straight lines in case of Cartesian sampling or more complex smooth curves in case of off-the-grid or non-Cartesian sampling. In this context, Parallel Imaging (PI) and Compressed Sensing (CS) are two main acceleration approaches for MRI data acquisition that consist in under-sampling the k-space, thus retaining fewer shots during k-space coverage. Alternative acceleration principles based on multiplexing like simultaneous multi-slice acquisition exist but they won’t be addressed in this course.

Parallel Imaging

Parallel imaging methods in 2D MRI first rely on uniform and deterministic under-sampling in k-space along the phase encoding direction, as the latter is responsible for scan duration. This principle directly extends to 3D MRI by under-sampling along the phase and partition directions, thereby allowing for higher acceleration factors in 3D MRI. To compensate for the violation of the Shannon-Nyquist criterion and the resulting aliasing artifacts (see Fig. 1), a phased array of RF coils is used to collect multiple data sets in k-space of the same object, which are modulated in space by spatially varying coil sensitivities. Various PI reconstruction methods have emerged in the late 90s and early 2000s, with different strategies to perform aliasing-free image reconstruction. Some proceed in the image domain (e.g. SENSE) while other in the k-space (e.g. GRAPPA).
SENSE method. The SENSitivity Encoding or SENSE method (Pruessmann et al, 1999) was first introduced in the literature and proceeds in the image domain (see Fig. 2). This means that each coil-specific k-space data is first inverse Fourier transformed. SENSE additionally assumes that coil sensitivity maps are known. This assumption is tenable as the coil sensitivity maps are spatially smooth, hence they can be extracted from the central lines in k-space (low frequencies), which often remain fully sampled for calibration purpose. Then, the full FOV image is reconstructed by solving in the least squares sense a sequence of small-dimension and independent linear systems that involve the aliased coil-specific images and coil sensitivities. As such, SENSE allows for an embarrassingly parallel implementation. The image reconstruction problem remains well posed (i.e. the solution is unique and stable) as long as the number of coils L is larger than the under-sampling factor R.
GRAPPA method. Alternatively, the Generalized Autocalibrating Partially Parallel Acquisition or GRAPPA technique (Griswold et al, 2002) is a k-space based reconstruction method that estimates kernel weights based on autocalibrated lines in the central portion of k-space. These weights relate the source k-space points from all coils to a target k-space point in one specific coil. For an under-sampling factor of R, (R-1) different kernels must be specified to fill the missing k-space data (or lines in Cartesian sampling) that have been skipped during under-sampling. Once all the lines for a particular coil are reconstructed, an inverse Fourier transform is applied to generate the uncombined MR image for that coil. Once this process has been repeated for each coil of the array, the full set of uncombined MR images can be obtained. These images are eventually combined in a single one using the square root of the sum of squares reconstruction.

Compressed Sensing

In the last 15 years, the application of Compressed Sensing (CS) theory (Donoho, 2006; Candès et al, 2006) to MRI (Lustig et al, 2007) has received considerable interest and led to major improvements in terms of accelerating data acquisition without degrading image quality in low acceleration regimes. In contrast to parallel imaging, CS for MRI relies on variable density sampling (VDS) to reach partially incoherent under-sampling. VDS here means that the center of k-space is more densely sampled than its periphery (see Fig. 3). Although VDS has been originally evidenced empirically (Lustig et al, 2007), it’s a key concept to theoretically minimize the number of measurements while preserving good image quality (Chauffert et al, SIAM 2014). VDS can be implemented in a suboptimal way (1D) along the phase encoding direction in 2D Cartesian acquisition, or more efficiently along complex trajectories (e.g. spirals, cones, etc) in non-Cartesian acquisitions. For an extensive survey of non-Cartesian trajectories in 2D and 3D imaging and basic non-Cartesian MR image reconstruction, the participants are invited to install the MRI-NUFFT Python package [https://github.com/mind-inria/mri-nufft].
However, depending on the selected target sampling density, analytical curves may no longer meet the hardware constraints on the gradients (maximum amplitude and slew rate). To overcome this limitation, the SPARKLING (Spreading Projection Algorithm for Rapid K-space sampLING) methodology has been developed (Lazarus et al, 2019; Chaithya GR et al, 2022). From a CS perspective, SPARKLING is optimal as it implements variable density sampling while being locally uniform (see Fig. 4).
Additionally, CS theory claims that exact image reconstruction can be achieved in noise-free scenarios even from an amount of k-space data far below the Shannon-Nyquist bound as long as the underlying MR image is sparse or compressible in an appropriate transform domain and sparsity is enforced during image reconstruction. Most MR images are compressible in hand-crafted transforms like wavelet basis (see Fig. 3), tight frames (e.g. curvelets) or in the image gradient domain. Then, adding an appropriate l1-norm regularization term to a quadratic data consistency cost function enforces sparsity and makes the MR image reconstruction problem well-posed (see Fig. 5). The resulting cost function is convex but nonsmooth. As such its global minimizer can be computed iteratively by any modern proximal gradient method like FISTA (Beck and Teboulle, 2009), POGM (Kim and Fessler, 2018) or any primal-dual method (Condat, 2013) depending on the analysis vs synthesis prior formulation that has been adopted in the definition of the cost function. Several open source packages in the MRI community implement these ideas, notably the BART toolbox [https://mrirecon.github.io/bart/], pysap-mri [https://github.com/CEA-COSMIC/pysap-mri] and SigPy [https://github.com/mikgroup/sigpy-mri-tutorial] in Python, the MIRT toolbox [https://github.com/JeffFessler/mirt] in Julia, etc.

Acknowledgements

This work was supported in part by 1) ANR (grant ANR-20-CE19-0027), 2) High Performance Computing (HPC) Resources of Institut du Développement et des Ressources en Informatique Scientifique (IDRIS) through Grand Equipement National de Calcul Intensif (GENCI) under Grant 2021-AD011011153 and 3) by the Cross-Disciplinary Program on Numerical Simulation of French Alternative Energies and Atomic Energy Commission for the project SILICOSMIC.