Accelerated Imaging Methods#

MRI scans can be accelerated by applying methods that reduce the number of k-space samples required, also known as undersampling. These require a corresponding advanced image reconstruction method to allow for the reduced sampling. The techniques covered here include parallel imaging, compressed sensing, and deep learning methods.

Learning Goals#

  1. Describe how images are formed

    • Describe how accelerated imaging methods work

  2. Understand the constraints and tradeoffs in MRI

    • Identify accelerated imaging methods

    • Understand when accelerated imaging methods can be applied

  3. Manipulate MRI sequence parameters to improve performance

    • Define which parameters are modified when accelerated imaging methods are applied

  4. Manipulate and analyze MRI data

    • Reconstruct an image from undersampled raw data

General Formulation of MRI Reconstruction#

For advanced image reconstruction methods, it is helpful to reformulate the reconstruction problem as a linear system and using standard mathematical notation describing linear systems as:

\[ y = \mathbf{E}x + n \]

where \(y\) is the acquired data, \(\mathbf{E}\) is the encoding matrix, \(x\) is the spatial distribution of the transverse magnetization (e.g. image), and \(n\) is noise. In this formulation, the image is vectorized. For example, 2D FT sampled data and the corresponding 2D image would be converted to:

\[\begin{split} y = \left[ \begin{array}{c} s_1(t_1) \\ s_1(t_2) \\ \vdots \\ s_1(t_{Nx}) \\ s_2(t_1) \\ \vdots \\ s_{Ny}(t_{Nx})\\ \end{array} \right]\end{split}\]

where \(s_m(t_n)\) is the data from the \(m\)th TR at the \(n\)th sample.

\[\begin{split} x = \left[ \begin{array}{c} m(x_1,y_1) \\ m(x_2,y_1) \\ \vdots \\ m(x_{Nx},y_1) \\ m(x_{1},y_2) \\ \vdots \\ m(x_{Nx},y_{Ny}) \end{array} \right]\end{split}\]

where \(x_i\) and \(y_i\) are equally spaced locations in image space.

The encoding matrix, \(\mathbf{E}\), at a minimum includes a discrete Fourier Transform matrix, \(\mathbf{F}\), representing our data is the Fourier Transform of the transverse magnetization, evaluated at k-space locations. For this matrix, the entries are defined by k-space location being sampled, \(\vec{k}_i\) and the spatial locations, \(\vec{r}_j\):

\[ \mathbf{f}_{i,j} = \exp(-i 2 \pi \vec{k}_i \cdot \vec{r}_j)\]

In this case, the measurement is the discrete Fourier Transform of the image

\[ y = \mathbf{F}x + n \]

and image reconstruction is performed by inverse discrete Fourier Transform, which for fully sampled 2D FT imaging is well defined by matrix inversion:

\[ \hat{x} = \mathbf{F^{H}} y \]

Parallel Imaging#

For parallel imaging (PI), we need to consider the coil sensitivity profiles, \(\mathbf{C}_q\), for each RF coil into encoding matrix along with a Fourier Transform encoding matrix, \(\mathbf{F}\), as well as a k-space sub-sampling operator, \(\mathbf{S}\), for the measurements from each RF coil, \(\mathbf{y}_q\):

\[ y_q = \mathbf{E}_q x + n_q = \mathbf{S} \mathbf{F} \mathbf{C_q} x + n_q \]

The coil sensitivity profiles \(\mathbf{C}_q\) is a diagonal matrix with entries corresponding to the coil sensitivity profile at each location in \(\mathbf{x}\).

The k-space sub-sampling operator \(\mathbf{S}\), is a diagnoal matrix with entries of 1 or 0, describing whether an expect grid location was sampled. It is convenient to use when describing image reconstruction as filling in points in k-space.

Here we can return to our original formulation by creating augmented matrices, concatenating along the coil elements dimension for example as:

\[\mathbf{y} = [y_1 \ y_2 \ldots y_N]\]
\[\mathbf{E} = [\mathbf{E}_1 \ \mathbf{E}_2 \ldots \mathbf{E}_N]\]
\[\mathbf{n} = [n_1 \ n_2 \ldots n_N]\]

resulting in:

\[ \mathbf{y} = \mathbf{Ex} + \mathbf{n} \]

Image-space Methods#

Image-space parallel imaging methods (e.g. SENSE) can be formulated as the following optimization problem

\[\mathbf{\hat{x}}_{PI} = \arg \min_\mathbf{x} \frac{1}{2} \| \mathbf{y} - \mathbf{Ex} \|^2_2\]

This can be solved directly by a pseudo-inverse, but requires measurements of the coil sensitivity profiles to calculate \(\mathbf{E}\):

\[\mathbf{\hat{x}}_{PI} = (\mathbf{E}^H\mathbf{E})^{-1} \mathbf{E}^H \mathbf{y} \]

The k-space sampling patterns used for these methods typically use regular undersampling, meaning there is a consistent pattern of acquired and skipped k-space lines.

Auto-calibrated K-space Methods#

K-space parallel imaging methods (e.g. GRAPPA) utilize a calibration kernel, computed from the data itself and captured in the matrix \(\mathbf{G}\). These can also be generally formulated as the following optimization problem

\[\mathbf{\hat{x}}_{PI} = \arg \min_\mathbf{x} \| \mathbf{y} - \mathbf{Ex} \|^2_2 + \lambda \| (\mathbf{G} - \mathbf{I}) \mathbf{x} \|^2_2 \]

Where \(\lambda\) is a tuning or regularization parameter and \(\mathbf{I}\) is the identity matrix. This simultaneous enforces data consistency (first term) as well as self-consistency of the calibration (sescond term).

The calibration kernel, \(\mathbf{G}\), can either be directly measured from the data as in GRAPPA, or iteratively estimated while solving the optimization, as in SPIRiT.

The k-space sampling patterns used for these methods typically use regular undersampling, meaning there is a consistent pattern of acquired and skipped k-space lines, combined with full sampling in the center of k-space.

SNR in Parallel Imaging#

There is an SNR penalty when using these methods that varies in severity depending on the conditioning of the undersampled reconstruction matrix. Typically, more RF coil elements in the direction of the undersampling leads to a more well-conditioned reconstruction and lower SNR penalty. Conversely, fewer RF coil elements in the direction of the undersampling leads to a more ill-conditioned reconstruction matrix and higher SNR penatly. Also, regions will little difference betweeen RF coils (typically in the center of the body), also tend to have larger SNR penalties.

This is characterized by the “g-factor”, where \(g(\vec{r}) \geq 1\) characterizes a spatially-varying SNR loss that is dependent on the coil loading and geometry, k-space sampling pattern, and parallel imaging reconstruction method. When appyling an acceleration factor of \(R\), the total readout time reduces the SNR by \(\sqrt{R}\) as well, leading to the SNR relationship:

\[SNR_{PI} = \frac{SNR_{full}}{g(\vec{r})\sqrt{R}}\]

Artifacts with Parallel Imaging#

Parallel Imaging artifacts typically appear as aliasing. The undersampled data will have aliasing when not using a parallel imaging reconstruction method, and sometimes the aliasing is not completely removed. There are also aliasing like artifacts when there is a mismatch between the data and measured coil sensitivity profiles.

Compressed Sensing Methods#

Compressed Sensing (CS) theory says that an image that is compressible in some domain can then be reconstructed from a subset of data samples. In other words, we can further accelerate our data acquisition.

Compressed Sensing is formulated as the following optimization problem, specifically using the \(\ell_1\)-norm (\(\|x\|_1 = \sum_{i=1}^n |x_i|\)) that promotes sparsity in the solution:

\[\mathbf{\hat{x}}_{CS} = \arg \min_\mathbf{x} \frac{1}{2} \| \mathbf{y} - \mathbf{Ex} \|^2_2 + \lambda_{CS} \| \mathbf{Wx} \|_1 \]

which includes a data consistency term where the data multiplied by the encoding matrix must match the reconstructed image, and a regularization term that enforces that the image is sparse using the \(\ell_1\) norm in some domain through the sparsifying transform, \(\mathbf{W}\). There is a regularization factor, \(\lambda_{CS}\), that must be chosen to balance the data consistency and sparsity terms.

These methods specifically require k-space sampling patterns with pseudo-random undersampling, which is illustratic by the diagram below of the undersampled reconstruction procedure:

Compressed Sensing Illustration

Caption (Figure from https://doi.org/10.1109/MSP.2007.914728 ): Heuristic procedure for reconstruction from undersampled data. A sparse signal (a) is 8-fold undersampled in its 1-D k -space domain (b). Equispaced undersampling results in signal aliasing (d) preventing recovery. Pseudo-random undersampling results in incoherent interference (c). Some strong signal components stick above the interference level, are detected and recovered by thresholding (e) and (f). The interference of these components is computed (g) and subtracted (h), thus lowering the total interference level and enabling recovery of weaker components.

Popular sparsifying transforms include total variation (TV), total generalized variation (TGV), and wavelets, where MRI data is typically more sparse. The pseudo-random undersampling typically uses a variable density that preferentially has more samples near the center of k-space.

This regularization term, \(\| \mathbf{Wx} \|_1\) can also be combined with iterative parallel imaging reconstruction methods as well, for example in the \(\ell_1\)-ESPIRiT method.

SNR in Compressed Sensing#

It is difficult to define SNR when using compressed sensing reconstruction methods, as they inherently perform some denoising when constraining the reconstruction based on sparsity. The apparent SNR can also vary significantly depending on the choice of the regularization factor, \(\lambda_{CS}\).

Artifacts with Compressed Sensing#

The most common artifacts from compressed sensing come from overfitting to the sparsity penalty. This result in either an over-smoothed, unnatural appearance or, in the case of a wavelet sparsifying transform, blocking artifacts.

Deep Learning Reconstructions#

Deep learning (DL) can be used to reconstructed undersampled k-space into image data from a subset of data samples as well. Conceptually, these methods can be trained to learn how to incorporate coil sensitivity information like parallel imaging and typical image sparsity patterns like compressed sensing. Since they learn from real-world data, they can learn information that is the most relevant to MRI data and thus have been shown to support higher acceleration factors. For example, they can learn what patterns within the image and/or k-space domain and most relevant and most common in MRI.

These methods can often be considered general solutions to a regularized least-squares objective

\[\mathbf{\hat{x}}_{reg} = \arg \min_\mathbf{x} \frac{1}{2} \| \mathbf{y} - \mathbf{Ex} \|^2_2 + \mathcal{R}(\mathbf{x}) \]

where \(\mathcal{R}(\cdot)\) can be a parallel imaging or compressed sensing regularizer as above, but can also be implicitly implemented via machine learning techniques. This means that the machine learning reconstruction is attempting to constrain or regularize the reconstruction.

The most straightforward approach using deep learning is to train a CNN to invert the encoding operation in image space, e.g. after applying the inverse Fourier Transform to the data:

\[\mathbf{\hat{x}}_{DL} = g_{DNN}(\mathbf{E^{H} y})\]

The learned deeep neural network, \(g_{NN}()\), includes many weights that must be trained prior to being able to perform the reconstruction.

However, so-called un-rolled networks have become more popular for MRI reconstruction, as they can better incorporate parallel imaging information, data consistency, and also are often more compact. They are designed to mimic iterations of classical algorithms to solve optimization problems, and each iteration has been “un-rolled” into a series of cascaded neural networks. Some popular methods include MoDL and Variational Networks.

Variational Networks (VarNets) for undersampled MRI reconstruction are formulated as the above regularlized least-squares objective, but models the architecture after a single gradient update in iterative gradient descent (used to solve optimization problems) and uses a CNN (e.g. UNet) within the iterations

Variational Network for MRI Reconstruction

Caption (from https://arxiv.org/abs/2004.06688): The E2E-VN model takes under-sampled k-space as input and applies several cascades, followed by IFT and RSS. Bottom: The DC module brings intermediate k-space closer to measured values, the Refinement module maps multi-coil k-space to one image, applies a U-Net and maps back to k-space, and the SME estimates sensitivity maps used in the refinement step.

The Model-based Deep Learning (MoDL) method uses the formulation where the regularizer is a CNN that estimates the noise and aliasing patterns

\[ \mathcal{R}_{MoDL}(\mathbf{x}) = \lambda \| \mathcal{N}_w(\mathbf{x}) \|^2\]

and the estimate depends on a set of learned parameters, \(w\), in the neural network, \(\mathcal{N}_\mathbf{w}(\cdot)\). It’s structure is overall similar to VarNets, as both are unrolled architectures, but with slightly different formulations.

Deep Learning Requirements#

  1. Data

    • Generally a lot of data (\(10^4\) in many radiology applications, but also depends a lot on the question being answered and data quality)

    • Fully sampled data typically required for supervised learning of undersampled MRI reconstruction

    • Data augmentation – realistic modifications of existing data to increase training data size

    • For reconstruction, fastMRI provides thousands of high quality training data: http://fastmri.med.nyu.edu, https://doi.org/10.48550/arXiv.1811.08839

  2. Architecture

    • U-net, ResNet popular for working in image-domain

    • Un-rolled networks popular for image reconstruction.

    • Generative models possible

    • Transformers

  3. Computation

    • GPUs are essential for training of deep networks

  4. Training

    • Send input data forward through network, compare to the expected output (e.g. fully sampled image or data) update the network weights through back-propagation

    • Choose loss function, such as l1 or l2

    • Choose batch sizes, Learning rate, etc

SNR in DL Reconstructions#

It is difficult to define SNR when using machine learning reconstruction methods, as they inherently perform some denoising when constraining the reconstruction.

Artifacts with DL Reconstructions#

A challenge with machine learning reconstructions is that artifacts can be difficult to identify. They can result in an over-smoothed appearance, but often a well-trained network will retain realistic texture and noise appearance features. These methods can also overfit to the learned anatomy and data, meaning features might be erased or hallucinated, although using physics-based techniques provide some assurances of data-consistency.

Reference#

For a recent, comprehensive reference on thes methods I reccommend:

Hammernik K, Küstner T, Yaman B, Huang Z, Rueckert D, Knoll F, Akçakaya M. Physics-Driven Deep Learning for Computational Magnetic Resonance Imaging: Combining physics and machine learning for improved medical imaging. IEEE Signal Process Mag. 2023 Jan;40(1):98-114. doi: 10.1109/msp.2022.3215288. Epub 2023 Jan 2. PMID: 37304755; PMCID: PMC10249732.

https://arxiv.org/pdf/2203.12215.pdf