Before discussing undersampling and the associated aliasing problem, let us first discuss the simple case of reconstructing an MR image, $\textbf{x}\in\mathbb{C}^{N_{1}\times N_{2}}$, from a fully sampled k-space measurement, $\textbf{y}^{full}\in\mathbb{C}^{N_{1}\times N_{2}}$, using the Fourier transform $\mathscr{F}(\cdot)$:

where $\eta\sim\mathcal{N}(0,\Sigma)$ is the associated measurement noise typically assumed to have a Gaussian distribution (Virtue and Lustig, 2017) when the k-space measurement is obtained from a single receiver coil.

Modern MR scanners support parallel acquisition using an array of overlapping receiver coils $n$ modulated by their sensitivities $S_{i}$. So Eqn. 1 changes to: $\textbf{y}_{i}^{full}=\mathscr{F}S_{i}\textbf{x}+\eta$, where $i=\{1,2,\cdots,n\}$ is the number of receiver coils. We use $\textbf{y}_{i}$ for the undersampled k-space measurement coming from the $i^{th}$ receiver coil. To speed up the data acquisition process, multiple lines of k-space data (for cartesian sampling) are skipped using a binary sampling mask $\mathcal{M}\in\mathbb{C}^{N_{1}\times N_{2}}$ that selects a subset of k-space lines from $\textbf{y}^{full}$ in the phase encoding direction:

An example of $\textbf{y}^{full}$, y, $\mathcal{M}$ is shown in Fig 1.

To estimate the MR image x from the measurement, a data fidelity loss function is typically used to ensure that the estimated data is as close to the measurement as possible. A typical loss function is the squared loss, which is minimized to estimate x:

We borrow this particular formulation from (Sriram et al., 2020a; Zheng et al., 2019). This squared loss function is quite convenient if we wish to compute the error gradient during optimization.

However, the formulation in Eqn. 3 is under-determined if data is undersampled and does not have a unique solution. Consequently, a regularizer $\mathcal{R}(\textbf{x})$ is typically added to solve such an ill-conditioned cost function: ¹¹1The regularization term, $\mathcal{R}(\textbf{x})$ is related to the prior, $p(\textbf{x})$, of a maximum a priori (MAP) extimation of x, i.e. $\hat{\textbf{x}}=\arg\min_{\textbf{x}}(-\log p(\textbf{y}|\textbf{x})-\log p(\textbf{x}))$. In fact, in Ravishankar et al. (2019) the authors loosely define $\beta\mathcal{R}(\textbf{x})=-\log p(\textbf{x})$, which promotes some desirable image properties such as spatial smoothness, sparsity in image space, edge preservation, etc. with a view to get a unique solution.

Please note that each $\mathcal{R}_{i}(\textbf{x})$ is a separate regularizer, while the $\lambda_{i}$s are hyperparameters that control the properties of the reconstructed image $\hat{\textbf{x}}$ while avoiding over-fitting. Eqn. 3 along with the regularization term can be optimized using various methods, such as (i) the Morozov formulation, $\hat{\textbf{x}}=\min\{\mathcal{R}(\textbf{x});\text{such that},||A\textbf{x}-\textbf{y}||\leq\delta\}$; (ii) the Ivanov formulation, i.e. $\hat{\textbf{x}}=\min\{||A\textbf{x}-\textbf{y}||;\text{such that},\mathcal{R}(\textbf{x})\leq\epsilon\}$; or (iii) the Tikonov formulation, $\hat{\textbf{x}}=\min\{||A\textbf{x}-\textbf{y}||+\lambda\mathcal{R}(\textbf{x})\}$, discussed in (Oneto et al., 2016).

In general, the Tikonov formulation can be designed using a physics based, sparsity promoting, dictionary learning, or a deep learning based model. But there are several factors that can cause loss in data quality (especially small anatomical details) such as inaccurate modeling of the system noise, complexity, generalizability etc. To overcome these limitations, it is essential to develop inverse mapping methods that not only provide good data fidelity but also generalize well to unseen and unexpected data. In the next section, we shall describe how DL methods can be used as priors or regularizers for MR reconstruction.

We begin our discussion by considering DL methods with learn-able parameters $\theta$. The learn-able parameters $\theta$ can be trained using some notion of a learning rule that we shall discuss in Sec. 3. A DL training process helps us to find a function $G_{DL}(\cdot)$ that acts as a regularizer to Eqn. 4 with an overarching concept of an inverse mapping, i.e; (please note that, we shall follow (Zheng et al., 2019) to develop the optimization formulation)

and z is a latent variable capturing the statistical regularity of the data samples, while c is a conditional random variable that depends on a number of factors such as: undersampling of the k-space (Shaul et al., 2020; Oksuz et al., 2019a; Shitrit and Raviv, 2017), the resolution of the image (Yang et al., 2017; Yuan et al., 2020), or the type of DL network used (Lee et al., 2019). Based on the nature of the learning, there are two types of DL methods known as generative models, and non-generative models. We shall start with a basic understanding of DL methods to a more in-depth study of different architectures in Secs. 4 and 5.

In generative modeling, the random variable z is typically sampled from a noisy Gaussian distribution, $\textbf{z}\sim\mathcal{N}(0,I)$ with or without the presence of the conditional random variable, i.e.;

There are various ways to learn the parameters of Eqn. 6. For instance, the Generative Adversarial Network (GAN) (Goodfellow et al., 2014) allows us to learn the generator function $G_{GEN}(\cdot)$ using an interplay between two modules, while the Variational Autoencoders (VAEs) (Kingma and Welling, 2013) learns $G_{GEN}(\cdot)$ by optimizing the evidence lower bound (ELBO), or by incorporating a prior in a Bayesian Learning setup as described in Section 4.2. It is shown in the literature that a generative model can efficiently de-alias an MR image that had undergone a $4\times$ or $8\times$ undersampling in k-space (Zbontar et al., 2018).

The non-generative models on the other hand do not learn the underlying latent representation, but instead learn a mapping from the measurement space to the image space. Hence, the random variable z is not required. The cost function for a non-generative model is given by:

The function $G_{NGE}(\cdot)$ is a non-generative mapping function that could be a Convolutional Neural Network (CNN) (Zheng et al., 2019; Akçakaya et al., 2019; Sriram et al., 2020b), a Long Short Term Memory (LSTM) (Hochreiter and Schmidhuber, 1997), or any other similar deep learning model. The non-generative models show a significant improvement in image reconstruction quality compared to classical methods. We shall describe the generative and the non-generative modeling based approaches in detail in Secs. 4 and 5 respectively. Below, we give a brief overview of the classical or the non-DL based approaches for MR image reconstruction.

In the literature, several approaches can be found that perform an inverse mapping to reconstruct the MR image from k-space data. Starting from analytic methods (Fessler and Sutton, 2003; Laurette et al., 1996), there are several works that provide ways to do MR reconstruction, such as the physics based image reconstruction methods (Roeloffs et al., 2016; Tran-Gia et al., 2016; Maier et al., 2019; Tran-Gia et al., 2013; Hilbert et al., 2018; Sumpf et al., 2011; Ben-Eliezer et al., 2016; Zimmermann et al., 2017; Schneider et al., 2020), the sparsity promoting compressed sensing methods (Feng and Bresler, 1996; Bresler and Feng, 1996; Candès et al., 2006), and low-rank based approaches (Haldar, 2013). All these methods can be roughly categorized into two categories, i.e. (i) GRAPPA-like methods: where prior assumptions are imposed on the k-space; and (ii) SENSE-like methods: where an image is reconstructed from the k-space while jointly unaliasing (or dealiasing) the image using sparsity promoting and/or edge preserving image regularization terms.

A few k-space methods estimate the missing measurement lines by learning kernels from an already measured set of k-space lines from the center of the k-space (i.e., the auto-calibration or ACS lines). These k-space based methods include methods such as SMASH (Sodickson, 2000), VDAUTOSMASH (Heidemann et al., 2000), GRAPPA and its variations (Bouman and Sauer, 1993; Park et al., 2005; Seiberlich et al., 2008). The k-t GRAPPA (Huang et al., 2005) takes advantage of the correlations in the k-t space and interpolates the missing data. On the other hand, sparsity promoting low rank based methods are based on the assumption that, when the image reconstruction follows a set of constraints (such as sparsity, smoothness, parallel imaging, etc.), the resultant k-space should follow a structure that has low rank. The low rank assumption has been shown to be quite successful in dynamic MRI (Liang, 2007), functional MRI (Singh et al., 2015), and diffusion MRI (Hu et al., 2019). In this paper we give an overview of the low-rank matrix approaches (Haldar, 2013; Jin et al., 2016; Lee et al., 2016; Ongie and Jacob, 2016; Haldar and Zhuo, 2016; Haldar and Kim, 2017) in Sec. 2.1. While in k-t SLR (Lingala et al., 2011), a spatio-temporal total variation norm is used to recover the dynamic signal matrix.

The image space based reconstruction methods, such as, the model based image reconstruction algorithms, incorporate the underlying physics of the imaging system and leverage image priors such as neighborhood information (e.g. total-variation based sparsity, or edge preserving assumptions) during image reconstruction. Another class of works investigated the use of compressed sensing (CS) in MR reconstruction after its huge success in signal processing (Feng and Bresler, 1996; Bresler and Feng, 1996; Candès et al., 2006). Compressed sensing requires incoherent sampling and sparsity in the transform domain (Fourier, Wavelet, Ridgelet or any other basis) for nonlinear image reconstruction. We also describe dictionary learning based approaches that are a spacial case of compressed sensing using an overcomplete dictionary. The methods described in (Gleichman and Eldar, 2011; Ravishankar and Bresler, 2016; Lingala and Jacob, 2013; Rathi et al., 2011; Michailovich et al., 2011) show various ways to estimate the image and the dictionary from limited measurements.

The main contributions of this paper are:

• •

  We give a holistic overview of MR reconstruction methods, that includes a family of classical k-space based image reconstruction methods as well as the latest developments using deep learning methods.

• •

  We provide a discussion of the basic DL tools such as activation functions, loss functions, and network architecture, and provide a systematic insight into generative modeling and non-generative modeling based MR image reconstruction methods and discuss the advantages and limitations of each method.

• •

  We compare eleven methods that includes classical, non-DL and DL methods on fastMRI dataset and provide qualitative and quantitative results in Sec. 7.

• •

  We conclude the paper with a discussion on the open issues for the adoption of deep learning methods for MR reconstruction and the potential directions for improving the current state-of-the-art methods.

Classical k-space based reconstruction methods are largely based on the premise that the missing k-space lines can be interpolated (or extrapolated) based on a weighted combination of all acquired k-space measurement lines. For example, in the SMASH (Sodickson, 2000) method, the missing k-space lines are estimated using the spatial harmonics of order $m$. The k-space signal can then be written as:

where, $w_{i}^{m}$ are the spatial harmonics of order $m$, $\Delta k_{2}=\frac{2\pi}{FOV}$ is the minimum k-space interval ($FOV$ stands for field-of-view), and $\textbf{y}(k_{1},k_{2})$ is the k-space measurement of an image x, $(k_{1},k_{2})$ are the co-ordinates in k-space along the phase encoding (PE) and frequency encoding (FE) directions, and $j$ represents the imaginary number. From Eqn. 8, one can note that the $m^{th}$-line of k-space can be generated using $m$-number of spatial harmonics and hence we can estimate convolution kernels to approximate the missing k-space lines from the acquired k-space lines. So, in SMASH a k-space measurement x, also known as a composite signal in common parlance, is basically a linear combination of $n$ number of component signals (k-space measurements) coming from $n$ receiver coils modulated by their sensitivities $S_{i}$, i.e.

We borrow the mathematical notation from (Sodickson, 2000) and represent the composite signal in Eqn. 9 as follows:

However, SMASH requires the exact estimation of the sensitivity of the receiver coils to accurately solve the reconstruction problem.

To address this limitation, AUTO-SMASH (Jakob et al., 1998) assumed the existence of a fully sampled block of k-space lines called autocalibration lines (ACS) at the center of the k-space (the low frequency region) and relaxed the requirement of the exact estimation of receiver coil sensitivities. The AUTO-SMASH formulation can be written as:

We note that the AUTO-SMASH paper theoretically proved that it can learn a kernel that is, $\sum_{i=1}^{n}w_{i}^{m}\textbf{y}_{i}(k_{1},k_{2})=\sum_{i=1}^{n}w_{i}^{m}\textbf{y}_{i}^{ACS}(k_{1},k_{2}+m\Delta k_{2})$, that is, it can learn a linear shift invariant convolutional kernel to interpolate missing k-space lines from the knowledge of the fully sampled k-space lines of ACS region. The variable density AUTO-SMASH (VD-AUTO-SMASH) (Heidemann et al., 2000) further improved the reconstruction process by acquiring multiple ACS lines in the center of the k-space. The composite signal x is estimated by adding each individual component $\textbf{y}_{i}$ by deriving linear weights $w^{m}_{i}$ and thereby estimating the missing k-space lines. The more popular, generalized auto calibrating partially parallel acquisitions (GRAPPA) (Bouman and Sauer, 1993) method uses this flavour of VD-AUTO-SMASH, i.e. the shift-invariant linear interpolation relationships in k-space, to learn the coefficients of a convolutional kernel from the ACS lines. The missing k-space are estimated as a linear combination of observed k-space points coming from all receiver coils. The weights of the convolution kernel are estimated as follows: a portion of the k-space lines in the ACS region are artificially masked to get a simulated set of acquired k-space points $\textbf{y}^{ACS_{1}}$ and missing k-space points $\textbf{y}^{ACS_{2}}$. Using the acquired k-space lines $\textbf{y}^{ACS_{1}}$, we can estimate the weights of the GRAPPA convolution kernel $K$ by minimizing the following cost function:

where $\circledast$ represents the convolution operation. The GRAPPA method has shown very good results for uniform undersampling, and is the method used in product sequences by Siemens and GE scanners. There are also recent methods (Xu et al., 2018; Chang et al., 2012) that show ways to learn non-linear weights of a GRAPPA kernel.

The GRAPPA method regresses k-space lines from a learned kernel without assuming any specific image reconstruction constraint such as sparsity, limited support, or smooth phase as discussed in (Kim et al., 2019). On the other hand, low-rank based methods assume an association between the reconstructed image and the k-space structure, thus implying that the convolution-structured Hankel or Toeplitz matrices leveraged from the k-space measurements must show a distinct null-space vector association with the kernel. As a result, any low-rank recovery algorithm can be used for image reconstruction. The simultaneous autocalibrating and k-space estimation (SAKE) (Shin et al., 2014) algorithm used the block Hankel form of the local neighborhood in k-space across all coils for image reconstruction. Instead of using correlations across multiple coils, the low-rank matrix modeling of local k-space neighborhoods (LORAKS) (Haldar, 2013) utilized the image phase constraint and finite image support (in image space) to produce very good image reconstruction quality. The LORAKS method does not require any explicit calibration of k-space samples and can work well even if some of the constraints such as sparsity, limited support, and smooth phase are not strictly satisfied. The AC-LORAKS (Haldar, 2015) improved the performance of LORAKS by assuming access to the ACS measurements, i.e.:

where $\mathcal{P}(\cdot)$ is a mapping function that transforms the k-space measurement to a structured low-rank matrix, and the matrix $N$ is the null space matrix. The mapping $\mathcal{P}(\cdot)$ basically takes care of the constraints such as sparsity, limited support, and smooth phase. In the PRUNO (Zhang et al., 2011) method, the mapping $\mathcal{P}(\cdot)$ only imposes limited support and parallel imaging constraints. On the other hand, the number of nullspace vectors in $N$ is set to 1 in the SPIRiT method (Lustig and Pauly, 2010). The ALOHA method (Lee et al., 2016) uses the weighted k-space along with transform-domain sparsity of the image. Different from them, the method of (Otazo et al., 2015) uses a spatio-temporal regularization.

These methods directly estimate the image from k-space by imposing prior knowledge about the properties of the image (e.g., spatial smoothness). Leveraging image prior through linear interpolation works well in practice but largely suffers from sub-optimal solutions and as a result the practical cone beam algorithm (Laurette et al., 1996) was introduced that improves image quality in such a scenario. The sensitivity encoding (SENSE) method (Pruessmann et al., 1999) is an image unfolding method that unfolds the periodic repetitions from the knowledge of the coil. In SENSE, the signal in a pixel location $(i,j)$ is a weighted sum of coil sensitivities, i.e.;

where $N_{2}$ is the height of image $\textbf{x}\in R^{N_{1}\times N_{2}}$, $N_{c}$ is the number of coils, and $S_{k}$ is the coil sensitivity of the $k^{th}$ coil. The $I_{k}$ is the k^th coil image that has aliased pixels at a certain position, and $i$ is a particular row and $j=\{1,\cdots,N_{2}\}$ is the column index counting from the top of the image to the bottom. The $S$ is the sensitivity matrix that assembles the corresponding sensitivity values of the coils at the locations of the involved pixels in the full FOV image x. The coil images $I_{k}$, the sensitivity matrix $S$, and the image x in Eqn. 14 can be re-written as;

By knowing the complex sensitivities at the corresponding positions, we can compute the generalized inverse of the sensitivity matrix:

Please note that, $I$ represents the complex coil image values at the chosen pixel and has length $N_{c}$. In k-t SENSE and k-t BLAST (Tsao et al., 2003) the information about the spatio-temporal support is obtained from the training dataset that helps to reduce aliasing.

The physics based methods allow statistical modeling instead of simple geometric modeling present in classical methods and reconstruct the MR images using the underlying physics of the imaging system (Roeloffs et al., 2016; Tran-Gia et al., 2016; Maier et al., 2019; Tran-Gia et al., 2013; Hilbert et al., 2018; Sumpf et al., 2011; Ben-Eliezer et al., 2016; Zimmermann et al., 2017; Schneider et al., 2020). These types of methods sometimes use very simplistic anatomical knowledge based priors (Chen et al., 1991; Gindi et al., 1993; Cao and Levin, 1997) or “pixel neighborhood” (Szeliski, 2010) information via a Markov Random Field based regularization (Sacco, 1990; Besag, 1986).

A potential function based regularization takes the form $\mathcal{R}(\textbf{x})=\sum_{i=2}^{n}\psi(\textbf{x}_{i}-\textbf{x}_{i-1})$, where the function, $\psi(\cdot)$, could be a hyperbolic, Gaussian (Bouman and Sauer, 1993) or any edge-preserving function (Thibault et al., 2007). The Total Variation (TV) could also be thought of as one such potential function. The (Rasch et al., 2018) shows a variational approach for the reconstruction of subsampled dynamic MR data, which combines smooth, temporal regularization with spatial total variation regularization.

Different from Total Variation (TV) approaches, (Bostan et al., 2012) proposed a stochastic modeling approach that is based on the solution of a stochastic differential equation (SDE) driven by non-Gaussian noise. Such stochastic modeling approaches promote the use of nonquadratic regularization functionals by tying them to some generative, continuous-domain signal model.

The Compressed Sensing (CS) based methods impose sparsity in the image domain by modifying Eqn. 2 to the following:

where $\Gamma$ is an operator that makes x sparse. The $l_{1}$ norm is used to promote sparsity in the transform or image domain. The $l_{1}$ norm minimization can be pursued using a basis pursuit or greedy algorithm (Boyd et al., 2004). However, use of non-convex quasi-norms (Chartrand, 2007; Zhao and Hu, 2008; Chartrand and Staneva, 2008; Saab et al., 2008) show an increase in robustness to noise and image non-sparsity. The structured sparsity theory (Boyer et al., 2019) shows that only $\mathcal{O}(M+M\log N)$ measurements are sufficient to reconstruct MR images when $M$-sparse data with size $N$ are given. The kt-SPARSE approach of (Lustig et al., 2006) uses a spatio-temporal regularization for high SNR reconstruction.

Iterative sparsity based methods (Ravishankar and Bresler, 2012; Liu et al., 2015, 2016) assume that the image can be expressed as a linear combination of the columns (atoms) from a dictionary $\Gamma$ such that $\textbf{x}=\Gamma^{T}\textbf{h}$ and h is the coefficient vector. Hence Eqn. 4 becomes:

The SOUP-DIL method (Bruckstein et al., 2009) uses an exact block coordinate descent scheme for optimization while a few methods (Chen et al., 2008; Lauzier et al., 2012; Chen et al., 2008) assume to have a prior image $\textbf{x}_{0}$, i.e. $\mathcal{R}(\textbf{x}-\textbf{x}_{0})$, to optimize Eqn. 4. The method in (Caballero et al., 2014) optimally sparsifies the spatio-temporal data by training an overcomplete basis of atoms.

The method in (Hongyi Gu, 2021) shows a DL based approach to leverage wavelets for reconstruction.

The transform based methods are a generalization of the CS approach that assumes a sparse approximation of the image along with a regularization of the transform itself, i.e., $\Gamma\textbf{x}=\textbf{h}+\epsilon$, where h is the sparse representation of x and $\epsilon$ is the modeling error. The method in (Ravishankar and Bresler, 2015) proposed a regularization as follows:

where $Q(\Gamma)=-\log|det\Gamma|+0.5||\Gamma||^{2}$ is the transform regularizer. In this context, the STROLLR method (Wen et al., 2018) used a global and a local regularizer.

In general, Eqn. 5 is a non-convex function and cannot be optimized directly with gradient descent update rules. The unrolled optimization algorithm procedure decouples the data consistency term and the regularization term by leveraging variable splitting in Eqn 4 as follows:

where the regularization is decoupled using a quadratic penalty on x and an auxiliary random variable z. Eqn 20 is optimized via alternate minimization of

and the data consistency term:

where the $\textbf{h}_{i},\textbf{h}_{i-1}$ are the intermediate variables at iteration $i$. The alternating direction method of multiplier networks (ADMM net) introduce a set of intermediate variables, $\textbf{h}_{1},\textbf{h}_{2},\cdots,\textbf{h}_{n}$, and eventually we have a set of dictionaries, $\Gamma_{1},\Gamma_{2},\cdots,\Gamma_{n}$, such that, $\textbf{h}_{i}=\Gamma_{i}\textbf{x}$, collectively promote sparsity. The basic ADMM net update (Yang et al., 2018) is as follows:

where $g(\cdot)$ can be any sparsity promoting operator and $\beta$ is called a multiplier. The iterative shrinkage thresholding algorithm (ISTA) solves this CS optimization problem as follows:

Later in this paper, we shall show how ISTA and ADMM can be organically used within the modern DL techniques in Sec. 5.3.

Perceptron: The journey of DL started in the year 1943 when Pitts and McCulloch (McCulloch and Pitts, 1943) gave the mathematical model of a biological neuron. This mathematical model is based on the “all or none” behavioral dogma of a biological neuron. Soon after, Rosenblatt provided the perceptron learning algorithm (Rosenblatt, 1957) which is a mathematical model based on the behaviour of a neuron. The perceptron resembles the structure of a neuron with dendrites, axons and a cell body. The basic perceptron is a binary classification algorithm of the following form:

where $\textbf{x}_{i}$’s are the components of an image vector x, $w_{i}$’s are the corresponding weights that determine the slope of the classification line, and $b$ is the bias term. This setup collectively resembles the “all or none” working principle of a neuron. However, in the famous book of Minsky and Papert (Minsky and Papert, 2017) called “Perpetron” it was shown that the perceptron can’t classify non-separable points such as an exclusive-OR (XOR) function.

The activation function, $\psi(\cdot)$, operates on a node or a layer of a neural network and provides a boolean output, probabilistic output or output within a range. The step activation function was proposed by McCulloch & Pitts in (McCulloch and Pitts, 1943) with the following form: $\psi(\textbf{x})=1,\text{if},\textbf{x}\geq 0.5$ and 0 otherwise. Several initial works also used the hyperbolic tangent function, $\psi(\textbf{x})=\tanh(\textbf{x})$, as an activation function that provides value within the range $[-1,+1]$. The sigmoid activation function, $\psi(\textbf{x})=\frac{1}{1+e^{-x}}$, is a very common choice and provides only positive values within the range $[0,1]$. However, one major disadvantage of the sigmoid activation function is that its derivative, $\psi^{{}^{\prime}}(\textbf{x})=\psi(\textbf{x})(1-\psi(\textbf{x}))$, quickly saturates to zero which leads to the vanishing gradient problem. This problem was addressed by adding a Rectified Linear Unit (ReLu) to the network (Brownlee, 2019), $\psi(\textbf{x})=\max(0,\textbf{x})$ with the derivative $\psi^{{}^{\prime}}(\textbf{x})=1,if,\textbf{x}>0$ or 0 elsewhere.