A Deep Cascade of Convolutional Neural Networks for MR Image Reconstruction

TL;DR

This paper introduces a deep cascade of convolutional neural networks for MRI reconstruction from undersampled data, significantly improving speed and accuracy over traditional compressed sensing methods, enabling real-time imaging.

Contribution

The paper presents a novel deep learning framework that outperforms existing methods in MRI reconstruction speed and quality, especially for highly undersampled data.

Findings

Outperforms compressed sensing methods in reconstruction error and perceptual quality.

Reconstructs images in 23 ms, enabling real-time applications.

Produces approximately half the error of previous state-of-the-art methods.

Abstract

The acquisition of Magnetic Resonance Imaging (MRI) is inherently slow. Inspired by recent advances in deep learning, we propose a framework for reconstructing MR images from undersampled data using a deep cascade of convolutional neural networks to accelerate the data acquisition process. We show that for Cartesian undersampling of 2D cardiac MR images, the proposed method outperforms the state-of-the-art compressed sensing approaches, such as dictionary learning-based MRI (DLMRI) reconstruction, in terms of reconstruction error, perceptual quality and reconstruction speed for both 3-fold and 6-fold undersampling. Compared to DLMRI, the error produced by the method proposed is approximately twice as small, allowing to preserve anatomical structures more faithfully. Using our method, each image can be reconstructed in 23 ms, which is fast enough to enable real-time applications.

Figures

Tables

Table 1. Table 1: DLMRI vs. CNN across 10 scans

	3-fold	6-fold
Models	MSE (SD) $\times {10}^{-3}$	MSE (SD) $\times {10}^{-3}$
DLMRI	2.12 (1.27)	6.31 (2.95)
CNN	0.89 (0.46)	3.42 (1.65)

Equations

$$\mathbf{y}={\mathbf{F}}_u\mathbf{x}$$

$$\underset{\mathbf{x}}{\text{min.}}$$

$$\underset{\mathbf{x},\mathbf{D},\{{\mathbf{\gamma }}_i\}}{\text{min.}}$$

$$\underset{\mathbf{x},\mathbf{\theta }}{\text{min.}}$$

$$\mathcal{L}(\mathbf{\theta })=\sum _{({\mathbf{x}}_u,{\mathbf{x}}_t)\in \mathcal{D}}\ell \left({\mathbf{x}}_t,{\mathbf{x}}_{\text{cnn}}\right)$$

$${\hat{\mathbf{x}}}_{\text{rec}}(k)=\{\begin{array}{ll}{\hat{\mathbf{x}}}_{\text{cnn}}(k)& \text{if }k\notin \Omega \\ \frac{{\hat{\mathbf{x}}}_{\text{cnn}}(k)+\lambda {\hat{\mathbf{x}}}_u(k)}{1+\lambda }& \text{if }k\in \Omega \end{array}$$

$$f_{dc}(\hat{\mathbf{x}},{\hat{\mathbf{x}}}_u,\lambda )=\mathbf{\Lambda }\hat{\mathbf{x}}+\frac{\lambda }{1+\lambda }{\hat{\mathbf{x}}}_u$$

$${\mathbf{\Lambda }}_{kk}=\{\begin{array}{ll}1& \text{if }k\notin \Omega \\ \frac{1}{1+\lambda }& \text{if }k\in \Omega \end{array}$$

$$f_L(\mathbf{x},{\hat{\mathbf{x}}}_u,\lambda )={\mathbf{F}}^{-1}\mathbf{\Lambda Fx}+\frac{\lambda }{1+\lambda }{\mathbf{F}}^{-1}{\hat{\mathbf{x}}}_u$$

$$\frac{\partial f_L}{\partial {\mathbf{x}}^T}={\mathbf{F}}^{-1}\mathbf{\Lambda F}$$

Full text

1 Introduction

In many clinical scenarios, medical imaging is an indispensable diagnostic and research tool. One such important modality is Magnetic Resonance Imaging (MRI), which is non-invasive and offers excellent resolution with various contrast mechanisms to reveal different properties of the underlying anatomy. However, MRI is associated with a slow acquisition process. This is because data samples of an MR image are acquired sequentially in $k$-space and the speed at which $k$-space can be traversed is limited by underlying MR physics. A long data acquisition procedures impose significant demands on patients, making the tool expensive and less accessible. One possible approach to accelerate the acquisition process is to undersample $k$-space, which in theory provides an acceleration rate proportional to a reduction factor of a number of k-space traversals required. However, undersampling in $k$-space violates the Nyquist-Shannon theorem and generates aliasing artefacts when the image is reconstructed. The main challenge in this case is to find an algorithm that takes into account the undersampling undergone and can compensate missing data with a-priori knowledge on the image to be reconstructed.

Using Compressed Sensing (CS), images can be reconstructed from sub-Nyquist sampling, assuming the following: firstly, the acquired images must be compressible, i.e. they have a sparse representation in some transform domain. Secondly, one must ensure incoherence between the sampling and sparsity domains to guarantee that the reconstruction problem has a unique solution and that this solution is attainable. In practice, this can be achieved with random sub-sampling of $k$-space, which translates aliasing patterns in the image domain into patterns that can be regarded as correlated noise. Under such assumptions, images can then be reconstructed through nonlinear optimization or iterative algorithms. The class of methods which applies CS to the MR reconstruction problem is termed CS-MRI [13]. A natural extension of these has been to enable more flexible representations with adaptive sparse modelling, where one attempts to obtain the optimal representation from data directly. This can be done by exploiting, for example, dictionary learning (DL) [17].

To achieve more aggressive undersampling, several strategies can be considered. One way is to further exploit the inherent redundancy of the MR data. For example, in dynamic imaging, one can make use of spatio-temporal redundancies [2], [9], [16]. Similarly, when imaging a full 3D volume, one exploit redundancy from adjacent slices [6]. An alternative approach is to exploit sources of explicit redundancy of the data and solve an overdetermined system. This is the fundamental assumption underlying parallel imaging [24]. Similarly, one can make use of multi-contrast information [7] or the redundancy generated by multiple filter responses of the image [14]. These explicit redundancies can also be used to complement the sparse modelling of inherent redundancies [8], [12].

Recently, deep learning has been successful at tackling many computer vision problems. Deep neural network architectures, in particular convolutional neural network (CNN), are becoming the state-of-the-art technique for various imaging problems including image classification [4], object localisation [18] and image segmentation [19]. Deep architectures are capable of extracting features from data to build increasingly abstract representations that are useful for the end-goal being considered, replacing the traditional approach of carefully hand-crafting features and algorithms. For example, it has already been demonstrated that CNNs outperform sparsity-based methods in super-resolution [3], not only for its quality but also in terms of the reconstruction speed [21]. One of the contributions of our work is to explore the application of CNNs in undersampled MR reconstruction and investigate whether they can exploit data redundancy through learned representations. In fact, CNNs have already been applied to compressed sensing from random Gaussian measurements [11]. Despite the popularity of CNNs, there has only been preliminary research on CNN-based MR image reconstruction [22], [25], hence the applicability of CNNs to this problem is yet to be qualitatively and quantitatively assessed in detail.

In this work we consider reconstructing 2D static images with Cartesian sampling using CNNs. Similar to the formulations in CS-MRI, we view the reconstruction problem as a de-aliasing problem in the image domain. However, reconstructing an undersampled MR image is challenging because the images typically have low signal-to-noise ratio, yet often high-quality reconstructions are needed for clinical applications. To resolve this issue, we propose a very deep network architecture which forms a cascade of CNNs. Our cascade network closely simulates the iterative reconstruction of DL-based methods, however, our approach allows end-to-end optimisation of the reconstruction algorithm. We show that under the Cartesian undersampling scheme, our CNN approach is capable of producing high-quality reconstructions of 2D cardiac MR images, outperforming DL-based MRI reconstruction (DLMRI) [17]. Moreover, using the proposed method, each images can be reconstructed in about 23 $ms$, which enables the real-time applications.

2 Problem Formulation

Let $\mathbf{x}\in\mathds{C}^{N}$ represent a complex-valued MR image composed of $\sqrt{N}\times\sqrt{N}$ pixels stacked as a column vector. Our problem is to reconstruct $\mathbf{x}$ from $\mathbf{y}\in\mathds{C}^{M}$, the measurements in $k$-space, such that:

[TABLE]

Here $\mathbf{F}_{u}\in\mathds{C}^{M\times N}$ is an undersampled Fourier encoding matrix. For undersampled $k$-space measurements ($M<<N$), the system of equations (1) is underdetermined and hence the inversion process is ill-defined. In order to reconstruct $\mathbf{x}$, one must exploit a-priori knowledge of its properties, which can be done by formulating an unconstrained optimisation problem:

[TABLE]

$\mathcal{R}$ expresses regularisation terms on $\mathbf{x}$ and $\lambda$ allows the adjustment of data fidelity based on the noise level of the acquired measurements $\mathbf{y}$. For CS-based methods, the regularisation terms $\mathcal{R}$ typically involve $\ell_{0}$ or $\ell_{1}$ norms in the sparsifying domain of $\mathbf{x}$. Our formulation is inspired by DL-based reconstruction approaches, in which the problem is formulated as:

[TABLE]

Here $\mathbf{R}_{i}$ is an operator which extracts an image patch at $i$, $\boldsymbol{\gamma}_{i}$ is the corresponding sparse code with respect to a dictionary $\mathbf{D}$. In this approach, the regularisation terms enforce $\mathbf{x}$ to be approximated by the reconstructions from the sparse code of patches. By taking the same approach, for our CNN formulation, we enforce $\mathbf{x}$ to be well-approximated by the CNN reconstruction:

[TABLE]

Here $f_{\text{cnn}}$ is the forward mapping of the CNN parameterised by $\boldsymbol{\theta}$, which takes in the zero-filled reconstruction $\mathbf{x}_{u}=\mathbf{F}^{H}_{u}\mathbf{y}$ and directly produces a reconstruction as an output. Since $\mathbf{x}_{u}$ is heavily affected by aliasing from sub-Nyquist sampling, the CNN reconstruction can therefore be seen as solving de-aliasing problem in the image domain.

The approach of eq. (LABEL:eq:cnn_rec), however, is limited in the sense that the CNN reconstruction and the data fidelity are two independent terms. In particular, since the CNN operates in the image domain, it is trained to reconstruct the image without a-priori information of the acquired data in $k$-space. However, if we already know some of the $k$-space values, then the CNN should be discouraged from modifying them. Therefore, by incorporating the data fidelity in the learning stage, the CNN should be able to achieve better reconstruction. This means that the output of the CNN is now conditioned on $\Omega$, an index set indicating which $k$-space measurements have been sampled in $\mathbf{y}$. Then, our final reconstruction is given simply by the output, $\mathbf{x}_{\text{cnn}}=f_{\text{cnn}}(\mathbf{x}_{u}|\boldsymbol{\theta},\lambda,\Omega)$. Given training data $\mathcal{D}$ of input-target pairs $(\mathbf{x}_{u},\mathbf{x}_{t})$, we can train the CNN to produce an output that attempts to accurately reconstruct the fully-sampled data by minimising an objective function:

[TABLE]

where $\ell$ is a loss function. In this work, we consider an element-wise squared loss, which is given by $\ell\left(\mathbf{x}_{t},\mathbf{x}_{\text{cnn}}\right)=\|\mathbf{x}_{t}-\mathbf{x}_{\text{cnn}}\|_{2}^{2}$.

3 Data Consistency Layer

In order to incorporate the data fidelity in the network architecture, we first note the following: for a fixed $\boldsymbol{\theta}$, eq. (LABEL:eq:cnn_rec) has a closed-form solution in $k$-space, given as in [17]:

[TABLE]

where $\hat{\mathbf{x}}_{\text{cnn}}=\mathbf{F}f_{\text{cnn}}(\mathbf{x}_{u}|\boldsymbol{\theta})$, $\hat{\mathbf{x}}_{u}=\mathbf{F}\mathbf{x}_{u}$ and $\mathbf{F}$ is the Fourier encoding matrix. The final image is reconstructed by applying the inverse of the encoding matrix $\mathbf{x}_{\text{rec}}=\mathbf{F}^{-1}\hat{\mathbf{x}}_{\text{rec}}$. In the noiseless setting (i.e. $\lambda\to\infty$), we simply replace the $i$th predicted coefficient by the original coefficient if it has been sampled. For this reason, this operation is called data consistency step in $k$-space (DC).

Since the DC step has a simple expression, we can in fact treat it as a layer operation of the network, which we denote as DC layer. When defining a layer of a network, the rules for forward and backward passes must be specified. This is because CNN training can effectively be performed through stochastic gradient descent, where one updates the network parameters $\boldsymbol{\theta}$ to minimise the objective function $\mathcal{L}$ by descending along the direction given by the derivative $\partial\mathcal{L}/\partial\boldsymbol{\theta}^{T}$. For this, it is necessary to define the gradients of each network layer relative to the network’s output. In practice, one uses an efficient algorithm called backpropagation [20], where the final gradient is given by the product of all the Jacobians of the layers contributing to the output. Hence, in general, it suffices to specify a layer operation $f_{L}$ for the forward pass and derive the Jacobian of the layer with respect to the layer input $\partial f_{L}/\partial\boldsymbol{\mathbf{x}}^{T}$ for the backward pass.

Forward pass

The data consistency in $k$-space can be simply decomposed into three operations: Fourier transform, data consistency and inverse Fourier transform. In our case, we take our Fourier transform to be a two-dimensional (2D) discrete Fourier transform (DFT) of the 2D image representation of $\mathbf{x}$, which is written as $\hat{\mathbf{x}}=\mathbf{F}\mathbf{x}$ in matrix form. The inverse transformation is defined analogously, where $\mathbf{x}=\mathbf{F}^{-1}\hat{\mathbf{x}}$. The data consistency $f_{dc}$ performs the element-wise operation defined in eq. (6). We can write it in matrix form as:

[TABLE]

Here $\bm{\Lambda}$ is a diagonal matrix of the form:

[TABLE]

Combining the three operations defined above, we can obtain the forward pass of the layer performing data consistency in $k$-space:

[TABLE]

Backward pass

In general, one requires Wirtinger calculus to derive a gradient in complex domain[1], however, in our case, the derivation greatly simplifies due to the linearity of the DFT matrix and the data consistency operation. The Jacobian of the DC layer with respect to the layer input $\mathbf{x}$ is therefore given by:

[TABLE]

There are several points that deserve further explanation: firstly, unlike many other applications where CNNs process real-valued data, MR images are complex-valued and the network needs to account for this. One possibility would be to design the network to perform complex-valued operations. A simpler approach, however, is to accommodate the complex nature of the data with real-valued operation in a dimensional space twice as large (i.e. we replace $\mathds{C}^{N}$ by $\mathds{R}^{2N}$). In the latter case, the derivations above still hold due to the fundamental assumption in Wirtinger calculus. Secondly, even though the DC layer does not have any additional parameters to be optimised, it allows end-to-end training of CNN, hence benefiting our final reconstruction.

4 Cascading Network

For CS-based methods, in particular for DLMRI, the optimisation problem is solved using a coordinate-descent type algorithm, alternating between the de-aliasing step and the data consistency step until convergence. In contrast, with CNNs, we are performing one step de-aliasing and the same network cannot be used to de-alias iteratively. While CNNs may be powerful enough to learn one step reconstruction, such network could indicate signs of overfitting, unless we have vast amounts of training data. In addition, training such networks may require a long time as well as careful fine-tuning steps. It is therefore best to be able to use CNNs for iterative reconstruction approaches.

A simple solution is to train a second CNN which learns to reconstruct from the output of the first CNN. In fact, we can concatenate a new CNN on the output of the previous CNN to build extremely deep networks which iterate between intermediate de-aliasing and the data consistency reconstruction. We term this a cascading network. In fact, one can essentially view this as unfolding the optimisation process of DLMRI. If each CNN expresses the dictionary learning reconstruction step, then the cascading CNN can be seen as a direct extension of DLMRI, where the whole reconstruction pipeline can be optimised from training.

5 Architecture and Implementation

Incorporating all the new elements mentioned above, we can devise our cascading network architecture. Our CNN takes in a two-channeled image $\mathds{R}^{\sqrt{n}\times\sqrt{n}\times 2}$, where each channel stores real and imaginary parts of the undersampled image. Based on literature, we used the following network architecture for CNN, illustrated in Figure 1: it has $n_{d}-1$ convolution layers $C_{i}$, which are all followed by Rectifier Linear Units (ReLU) as a choice of nonlinearity. For each of them, we used a kernel size $k=3$ [23] and the number of filters were set to $n_{f}=64$. The network is followed by another convolution layer $C_{\text{rec}}$ with kernel size $k=3$ and $n_{f}=2$, which projects the extracted representation back to image domain. We also used residual connection [5], which sums the output of the CNN module with its input. Finally, we form a cascading network by using the DC layers interleaved with the CNN reconstruction modules. For our experiment, we chose $n_{d}=5$ and $n_{c}=5$. We found that our choice of hyperparameters work sufficiently well, however, by no means were they optimised. Hence the result is likely to be improved by changing the architecture and varying the parameters such as kernel size and stride [19], [26].

As mentioned, pixel-wise squared error was used as our objective function. The minibatch size was set to 10, however, for the deeper models with large number of cascades, the minibatch size was reduced to fit the model on a single GPU memory. We initialised our network weights using He initialisation [5]. Adam optimiser [10] was used to train all models, with the parameters $\alpha=10^{-4},\beta_{1}=0.9$ and $\beta_{2}=0.999$. We also added $\ell_{2}$ weight decay of $10^{-7}$. The network was implemented in Python using Theano and Lasagne libraries.

6 Experimental Results

6.1 Setup

Dataset

Our method was evaluated using the cardiac MR dataset used in [2], consisting of 10 fully sampled short-axis cardiac cine MR scans. Each scan contains a single slice SSFP acquisition with 30 temporal frames with a $320\times 320$ mm field of view and 10 mm slice thickness. The raw data consists of 32-channel data with sampling matrix size $192\times 190$, which was zero-filled to the matrix size $256\times 256$. The data was combined into a single complex-valued image using SENSE [15] with no undersampling, retrospective gating and the coil sensitivity maps normalised to a body coil image. The images were then retrospectively undersampled using Cartesian undersampling masks, where we fully sampled along the frequency-encoding direction but undersample in the phase-encoding direction. The strategy was adopted from [9]: for each frame we acquired eight lowest spatial frequencies. The sampling probability of other frequencies along the phase-encoding direction was determined by a zero-mean Gaussian distribution. The acceleration rates are stated with respect to the matrix size of the raw data. Note that similarly to previous studies, [17], [2], since the raw data was combined prior to the simulation, the coil sensitivities were not directly addressed in our reconstruction. This is set for future investigation, where we plan to incorporate the explicit redundancy created by parallel imaging into our model.

Although the dataset is a dynamic sequence, we restrict our experiments to the 2D case only. Therefore, each time frame was treated as an independent image, yielding a total of 300 images. We found that applying rigid transformations as a data augmentation was crucial, as without it, the network quickly overfitted the training data. Moreover, for a fixed undersampling rate, we generated an undersampling mask on-the-fly to allow the network to learn diverse patterns of aliasing artefact.

Metric

We evaluated our method by reconstructing undersampled images from 3-fold and 6-fold acceleration rates. We used mean squared error (MSE) as our quantitative measure. During our experiment, we noticed that even for the same undersampling rate, different undersampling masks yield considerable differences in the reconstruction’s signal-to-noise. To take this into consideration for fair comparison, we assigned an arbitrary but fixed undersampling mask for each image in test data. Apart from the quantitative measure, we also inspected the visual aspect of the reconstructed images for qualitative assessment.

Models

For CNN, we selected the hyperparameters described above. To ensure a fair comparison, we reported the aggregated test result from 2-fold cross-validation (i.e. train on five subjects and test on the other five). For each iteration of the cross validation, the network was initialised using He initialisation, trained end-to-end. For 6-fold undersampling, we initialised the network using the parameters obtained from the trained models from 3-fold acceleration and fine-tuned using Adam optimiser. Each network converged within 3 days on GeForce TITAN X.

We compared our method to DLMRI, a representative of the state-of-the-art CS-based methods. For DLMRI, we simplified the implementation of DLTG from [2], with patch size $6\times 6$. We switched off any de-aliasing along the temporal axis. Since DLMRI is quite time consuming, in order to obtain the results within a reasonable amount of time, we trained a joint dictionary for all time frames within the subject and reconstructed them in parallel. Note that we did not observe any decrease in performance from this approach. For each subject, we ran 400 iterations and obtained the final reconstruction.

6.2 Results

The means of the reconstruction errors across 10 subjects are summarised in table 1. For both 3-fold and 6-fold acceleration, one can see that CNN consistently outperformed DLMRI, and that the standard deviation of the error made by CNN was smaller. The reconstruction from 3-fold acceleration can be found in Figure 2. It can be seen that the CNN approach produced a smaller overall error. The CNN reconstruction produced a more homogeneous reconstruction. On the other hand, DLMRI gave a blocky reconstruction. In some cases, both CNN and DLMRI suffered from small losses of important anatomical structures in their reconstructions (orange), but CNN was able to recover more details (red). The reconstructions from 6-fold acceleration is in Figure 3. Although both methods suffered from significant loss of structures (orange), CNN was still capable of better preserving the texture than DLMRI (red). On the other hand, DLMRI created extremely block-like artefacts due to over-smoothing. 6x undersampling for these images typically approaches the limit of sparsity-based methods, however, CNN was able to predict some anatomical details which was not possible by DLMRI. This could be due to the fact that CNN has more free parameters to tune with, allowing the network to learn complex but more accurate transformations of data.

Comparison of Reconstruction Speed

While training CNN is time consuming, once it is trained, the inference can be done extremely quickly on a GPU. Reconstructing each slice took $23\pm 0.1$ milliseconds on GeForce GTX 1080, which enables real-time applications. To produce the above results, DLMRI took about $6.1\pm 1.3$ hours per subject on CPU. Even though we do not have a GPU implementation of DLMRI, it is expected to take longer than 23ms because DLMRI requires dozens of iterations of dictionary learning and sparse coding steps. Using a fixed, pre-trained dictionary could remove this bottleneck, in exchange of lowering the reconstruction capacity.

7 Discussion and Conclusion

In this work, we evaluated the applicability of CNNs for the MR image reconstruction problem. From the experiment, we have shown that using the network with interleaved data consistency stages, we can obtain a model which can reconstruct images sufficiently well. The CS framework offers mathematical guarantee for the signal recovery, which makes the approach appealing in theory as well as in practice even though the required sparsity cannot be genuinely achieved in medical imaging. However, even though this is not the case for CNNs, we have empirically shown that a CNN-based approach can outperform DL-based MR reconstruction. In addition, at very aggressive undersampling rates, the CNN method was capable of reconstructing most of the anatomical structures more accurately, while CS-based methods do not guarantee such behaviour.

The limitation of this work is that the data was first reconstructed by SENSE, which was then used to simulated the acquisition process. It is, however, more practical to consider images with sensitivity map of the surface coils, which allows the model to be used for parallel imaging reconstruction directly. In fact, a better approach is to exploit the redundancy of the coil sensitivity maps and combine directly into our model, which will be addressed in our future work.

In this work, we were able to show that the network can be trained using arbitrary Cartesian undersampling masks of the fixed sampling rate rather than selecting a fixed number of undersampling masks for training and testing. This suggests that the network was capable of learning a generic strategy to de-alias the images. A further investigation should consider how tolerant the network is for different undersampling rates. Furthermore, it is interesting to consider other sampling patterns such as radial and spiral trajectories. As these trajectories provide different properties of aliasing artefacts, a further validation is appropriate to determine the flexibility of our approach.

Finally, Although CNNs can only learn local representations which should not affect global structure, it remains to be determined how the CNN approach operates when there is a pathology present in images, or other more variable content. We have performed a two-fold cross-validation to ensure that the network can handle unseen data acquired through the same acquisition protocol. Generalisation properties must be evaluated carefully on a larger dataset, however, CNNs are flexible in a way such that one can incorporate application specific priors to its objective to allocate more importance on preserving any features of interest in the reconstruction, provided that such expert knowledge is available at training time. For example, analysis of cardiac images in clinical settings often employs segmentation and/or registration. Multi-task learning is a promising approach to further improve the utility of CNN-based MR reconstructions.

Acknowledgment

The work was partially funded by EPSRC Programme Grant (EP/P001009/1).