Advanced Image Reconstruction

Advanced Image Reconstruction#

This section introduces a more general formulation of image reconstruction in MRI. This formulation can then be used to describe a variety of basic and advanced image reconstruction methods, several of which are described.

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.

Fourier Transform Reconstruction#

The definition of \(\mathbf{E}\) is critical, and this can be derived from the MRI signal equation. If we only consider the effects of magnetic field gradients, then

\[s(t) = \int m(\vec{r}) \exp(-i2\pi \vec{k}(t) \cdot \vec{r}) \ d\vec{r}\]

Our measurements, \(y\), taken from \(s(t)\) at discrete time points. Discretizing this relationship for Cartesian sampled data:

\[\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 repetition/TR \(m\) at sample \(n\).

We want to reconstruct the following image on a grid:

\[\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.

In this most simple case, the encoding matrix, \(\mathbf{E}\), includes just a discrete Fourier Transform matrix, \(\mathbf{F}\), describing the fact that our data is the Fourier Transform of the transverse magnetization, evaluated at k-space locations, as shown in the above signal equation. 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 \]

This can equivalently be computed using the Fast Fourier Transform (FFT) algorithm when the data is fully sampled on a Cartesian grid.

Parallel Imaging#

For parallel imaging (PI), we take advantage and include in our model the coil sensitivity profiles, \(\mathbf{C}_q\), for each RF coil. These are included in the encoding matrix along with a Fourier Transform encoding matrix, \(\mathbf{F}\). For performing undersampled reconstruction, a k-space sub-sampling operator, \(\mathbf{S}\), is also added to the encoding matrix, acknowledging the fact we know the sampling criteria are not satisfied, and allowing \(\mathbf{F}\) and \(\mathbf{F^H}\) to be computed with the FFT. With these factors, the measurements from each RF coil, \(\mathbf{y}_q\) are:

\[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 diagonal matrix with entries of 1 or 0, describing whether an expect grid location was sampled.

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}\]

Model-based Methods#

Model-based parallel imaging methods use some simulation of the MRI system and include the methods of SENSE, ESPIRiT, and NLINV. These reconstruct a single underlying image by solving an optimization problem. The popular SENSE methods solves

\[\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}\]

Linear Predictability Methods#

Linear predictability parallel imaging methods assume that each k-space sample is a linear combination of other k-space samples and inlcude the methods of AUTO-SMASH, GRAPPA, and SPIRiT. These utilize a calibration kernel, computed from the data itself and captured in the matrix \(\mathbf{G}\), and can 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 simultaneously enforces data consistency (first term) as well as self-consistency of the calibration (second 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.

Compressed Sensing Methods#

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.

Popular sparsifying transforms, \(\mathbf{W}\), include total variation (TV), total generalized variation (TGV), and wavelets, where MRI data is typically sparse.

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.

Deep Learning Reconstructions#

Deep learning (DL) 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 deep learning techniques. This means that the deep 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 deep neural network, \(g_{DNN}()\), 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

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

Reference#

For a recent, comprehensive reference on these 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

Contents