Deep Unfolding Convolutional Dictionary Model for Multi-Contrast MRI Super-resolution and Reconstruction

1 Introduction

Magnetic Resonance Imaging (MRI) is a noninvasive and non-ionizing medical imaging technique widely used for medical diagnosis, clinical analysis, and staging of disease [Song et al., 2019a]. However, due to the physics and physiological constraints, the MRI scanning procedure is always time-consuming. The scanning time of each patient can take up to ten minutes long, which seriously affects the patient experience and leads to high costs [Feng et al., 2022]. To accelerate the MRI procedure, a commonly-used approach is to undersample the k-space measurements and then use the reconstruction algorithm to restore its fully-sampled imaging. The recent mainstream methods for accelerating MRI include MRI reconstruction and MRI super-resolution (SR). MRI reconstruction aims to remove the aliasing artifacts caused by k-space undersampling. MRI SR aims to reconstruct highresolution (HR) MR images from low-resolution (LR) images. Both two methods can effectively improve the quality of MR images and accelerate MRI without upgrading hardware facilities.

∗Corresponding author

In the past few years, numerous MRI SR [Zhao et al., 2019; Wang et al., 2021; Zhang et al., 2021b] and reconstruction [Wang et al., 2016; Zheng et al., 2019; Nitski et al., 2020] methods have been proposed. However, they only use single-contrast MR images, ignoring complementary information frommulti-contrast images. In clinical practice, scanners usually produce exquisite images of body parts with different contrasts, such as T1-weighted, T2-weighted, and proton density weighted images (PD). These contrasts capture the information with the same underlying anatomy and thus exhibit similar structures [Lyu et al., 2020]. Furthermore, since the repetition time and echo time of different contrasts are different. For example, the scanning time of T1 and PD are significantly shorter than those of T2, so it is promising to leverage an HR (or fully sampled) reference image with a shorter acquisition time to reconstruct the modality with a longer scanning time from an LR (or undersampled) image. We name them guided-SR and guided-reconstruction.

Recently, deep convolutional neural network (DCNN) [Liu and Lei, 2021] has become a mainstream approach for solving MRI guided-SR problems. [Lyu et al., 2020] proposed to fuse the multi-contrast information in the high-level feature space. [Feng et al., 2021] used a multi-stage integration network to explore the rich relations among different contrasts. [Li et al., 2022; Fang et al., 2022] designed a transformerbased framework to fully explore the long-range dependencies in both reference and target images. For solving the multi-contrast reconstruction problem, [Xiang et al., 2018] fused multi-contrast MR acquisitions through deep networks to speed up the reconstruction of a certain target image. [Sun et al., 2019] proposed a deep information-sharing network to jointly reconstruct multi-contrast MRI images. [Dar et al., 2020] employed an auxiliary contrast as prior information into the generator of a generative adversarial network (GAN) to accelerate MR imaging. [Zhang et al., 2021a] applied the temporal feature fusion block based on ADMM to achieve multi-modal MRI reconstruction. [Feng et al., 2022] proposed a multi-modal transformer (MTrans) to transfer multiscale features from the target modality to the auxiliary one.

Although these methods have got promising results, their network architectures are black-boxes, making it difficult to understand the underlying mechanism and lack certain physical interpretations. As is well known, the model-based method [Song et al., 2019b] is explainable and easy to incorporate domain knowledge. Most recently, numerous welldesigned model-driven DCNN methods have been proposed to solve various image restoration problems, including deraining [Wang et al., 2020], debluring [Liu et al., 2019], image fusion [Liu et al., 2021] and so on. For MR images, [Yang et al., 2018; Aggarwal et al., 2018; Sriram et al., 2020] proposed model-based reconstruction methods based on the compressed sensing theory. However, they are designed for single-contrast MRI reconstruction.

In this paper, we propose a novel multi-contrast convolutional dictionary (MC-CDic) model to explicitly model the correlations of multi-contrast images for guided-SR and reconstruction tasks. MC-CDic model is designed by the following two observations. First, the multi-contrast images are sampled under the same anatomical structure, they thus share some common features, e.g., edges and texture information. Second, since the two images are generated by different scan settings, they also have their unique contrast information and some unaligned inconsistent structures, as shown in Figure 1. To better reconstruct the target image, we hope to transfer the common texture information from the reference image to the target image while avoiding the interference of inconsistent information. Taking these into consideration, we first extract the common features and unique features from the multicontrast images and then fuse these features to generate the desired high-quality images.

To extract the features from both images, firstly, we adopt the convolutional sparse coding (CSC) technique to model the multi-contrast MR images. Specifically, the multi-contrast images are decomposed into two types of features (i.e., common features and unique features) using two types of dictionaries (i.e., common filters and unique filters). Secondly, we employ implicit regularizations on the decomposed features and use the proximal gradient algorithm to optimize them. Thirdly, we unfold the iterative optimization process into a deep network. Especially, the unfolding process strictly follows the rules of the solution and the proximal operators are replaced by learnable residual blocks. Thus all the parameters can be learned in an end-to-end manner. Fourthly, we further employ multi-scale dictionaries to extract the multiscale sparse representations in each iteration. In this way, the proposed network is not only interpretable but also powerful. Our contributions can be summarized as follows.

• We propose a deep unfolding network based on convolutional dictionary learning for multi-contrast MRI SR and reconstruction. Different from existing methods that design unsophisticated fusion rules, we develop priors and constraints plugged into the model according to the characteristics of multi-contrast MR images.
• We build an observation model for multi-contrast MRI using the CSC technique and use the proximal gradient algorithm to solve the model. Based on the iterative solutions, we unroll the iteration steps into a multi-scale convolutional dictionary network and the proximal operators are replaced by learnable deep residual networks.
• The proposed model combines the good interpretability of the dictionary model with powerful deep neural networks. We test it on guided MRI SR and reconstruction tasks. Extensive experiments on the benchmark datasets demonstrate the effectiveness of the proposed model.

2 MC-CDic Model for MR Images

2.1 Model Formulation

Let ˆx1 ∈ Rm×n and ˆx2 ∈ Rm×n denote two classes of HR (or fully-sampled) MR images with different contrasts (e.g. PD and T2). The corresponding LR (or undersampled) images of ˆx2 is x2. For the guided-SR task, our MC-CDic model aims to restore the HR image ˆx2 from LR image x2 with the aid of the reference HR image ˆx1. Since the two contrast images are sampled under the same anatomical structure, they share some common features. Since they are smapled under different scan setting, they also have their unique features. Based on these observations, inspired by [Deng and Dragotti, 2020], we model the multi-contrast images as:

where ⊗ denotes the convolutional operation. {dck }Kk=1 and {hck }Kk=1 are the common filters of ˆx1 and x2, {duk }Kk=1 and {hvk }Kk=1 are the unique filters of ˆx1 and x2, {ck}Kk=1 represent the common features between ˆx1 and x2, {u}Kk=1 and {v}Kk=1 represent the unique features of ˆx1 and x2. K denotes the channel numbers of the feature representations. Based on the decomposed common features and unique features, the reconstructed HR image can be represented as:

where {qv k }Kk=1 and {qck }Kk=1 are the common filters and the unique filters for reconstructing the HR images x′2 .

For the above modeling, the main problem is to extract the common features and the unique features from the input images. Assuming all the filters are known, they can be solved by optimizing the following problem:

where {C,U, V } ∈ Rm×n×K are the stacked features response to {c, u, v}Kk=1 , λc, λu, λv are the trade-off parameters, ψc(·), ψu(·), ψv(·) are the regularizers to deliver the prior information to the feature maps.

2.2 Optimization Algorithm

To solve this problem, we alternately update each variable with other variables fixed. It leads to the following three subproblems:

Next, we optimize these three subproblems separately using the proximal gradient method. The details are as follows:

Updating U. The unique feature of ˆx1, U can be updated by solving the quadratic approximation [Beck and Teboulle, 2009] of the subproblem (5) as:

where U t-1 is the updating result of the last iteration, ηu is the stepsize parameter, and

where ⊗⊤ denotes the transposed convolution, Du ∈ R1×n×n×K is a 4-D tensor stacked by {duk }Kk=1 . Eq. (10) can be re-writed as:

where Dc ∈ R1×n×n×K is a 4-D tensor stacked by {dck }Kk=1 . Corresponding to general regularization term [Donoho, 1995], the solution of Eq. (8) is:

where proxηuλu (·) is the proximal operator dependent on the regularization term ψu(·) with respect to U. Instead of choosing a fixed regularizer in the model, the form of the proximal operator can be automatically learned from training data.

Updating V . Similarly, the quadratic approximation of the problem (6) with respect to V is:

where

Hv ∈ R1×n×n×K is a 4-D tensor stacked by {hvk }Kk=1 . Eq. 15 can be re-writed as:

where Hc ∈ R1×n×n×K is a 4-D tensor stacked by {hck }Kk=1 .

The solution of Eq. (13) can be represented as:

where proxηvλv (·) is the proximal operator correlated to the regularization term ψv(·) with respect to V .

Updating C. We first denote ˜x1 = ˆx1 - ∑k duk ⊗ uk and ˜x2 = x2 - ∑k hvk ⊗ vk, Eq. (7) can be re-written as follows:

Then the first two terms in Eq. (18) can be further combined and this yields the following optimization problem:

where M is the concatenate of ˜x1 and ˜x2, lck is the concatenation of dck and hck . Thus, Eq. (19) can be solved by the same way as U and V . Specifically, the quadratic approximation of this problem with respect to C is:

where

where Lc ∈ Rn×n×K×2 is a 4-D tensor stacked by {lck }Kk=1 . Eq. (22) can be re-writed as:

The solution of Eq.(20) can be represented as:

where proxηcλc (·) is the proximal operator correlated to the regularization term ψc(·) with respect to C.

3 Deep UnfoldingMC-CDic Model

Before unfolding the above model, we first formulate some relaxations on the dictionaries. Following earlier works [Simon and Elad, 2019; Lecouat et al., 2020b; Lecouat et al., 2020a; Liu et al., 2022], we untie the convolutional dictionaries from their transposed convolutional dictionaries in Eq. (11), (16) and (23) for better performance. To be specific, we initialize the convolution dictionary Du⊗ identically to Du⊗⊤ but are allowed to evolve independently during training. For convenience, we use ̂Du⊗ to replace Du⊗⊤. Similarly, Hv⊗⊤ and Lc⊗⊤ are also replaced by ̂Hv⊗ and ̂Lc⊗ in our deep unfloding model.

3.1 Network Design

(1) Feature Initialization Module. We use three convolutional dictionaries to initialize the common features and the unique features from the multi-contrast images ˆx1 and x2. Specifically, we first use the dictionaries ̂Du and ̂Hv to initialize the unique feature U 0 and V 0 respectively. Then we concatenate ˆx1 and x2, and use ̂Lc to initialize the common feature C 0 .

(2) MC-CDic Iteration Module. As shown in the figure, our model consists T MC-CDic blocks, representing T iterations of the algorithm for solving Eq. (4). Specifically, we take the input images ˆx1 and x2, and the previous outputs U t , V t and C t as inputs, and outputs the updated U t+1 , V t+1 and C t+1 in each iteration block. The iteration blocks are unfolded following the updating rules of Eq. (12), (17) and (24). To be specific, the convolutional dictionary operations can be easily implemented by the commonly-used convolutional layers in normal networks. The key issue of unrolling the algorithm is how to represent the three proximal operators proxηuλu (·), proxηvλv (·) and proxηcλc (·). Inspired by recent unfolding networks [Wang et al., 2020; Wang et al., 2022], we employ the deep residual network (ResNet) to replace the proximal operators. Note that the T MC-CDic iterations share the same parameters in our model.

Thus, at the tth stage, the network for updating U can be built as:

where proxNetu(·) denotes the ResNet with several resblocks. Similarly, the network for updating V at the tth stage can be built as:

where proxNetv(·) denotes the ResNet with several resblocks.

Based on the updated unique features U t and V t , the network for updating C at the tth stage can be built as:

where proxNetc(·) denotes the ResNet with several resblocks, M t is the concatenated features of ¯xt1 and ¯xt2 . According to the updating rules of U, V and C in Eq. (25), (26) and (27), we design the multi-contrast convolutional dictionary (MC-CDic) block. Its overall structure is shown in Figure 2. All the parameters involved can be automatically learned from training data in an end-to-end manner.

(3) Reconstruction Module. After T MC-CDic iterations, we get the final common features C T and the final unique feature U T and V T . Then the reconstructed HR image x′2 can be obtained by:

where P c = {pck }Kk=1 and P u = {puk }Kk=1 are the common filters and unique filters for reconstructing the HR images.

3.2 Multi-Scale Dictionaries

The importance of multi-scale representations has been proven in numerous computer vision tasks. However, most existing CDic models only use the single-scale dictionary to extract single-scale representations, which limits their performance to a large extent. Inspired by [Liu et al., 2022], we employ the multi-scale dictionaries in our MC-CDic model to perceive the multi-scale representations. To be specific, we employ a stripped-down version of the image synthesis process of U-Net by removing all non-linearities, batch normalization, and additive biases. We regard the decoder part of the UNet as a normal CDic layer and regard the encoder part of the UNet as its transposed CDic layer. As shown in Figure 3 (a) and (b), we present the operations of U t = ˆDu ⊗ s and s = Du ⊗ U t . The proximal network for multi-scale dictionaries is shown in Figure 3 (c).

3.3 Loss Function

We employ the simple L1 loss to supervise the multi-contrast images simultaneously. It can be formulated as:

where α is a trade-off parameter to balance the importance of different contrasts. In this paper, we set α = 0.1. Note that the loss on the reference image is to ensure that the reference image can be efficiently decomposed.

4 Experiments

4.1 Datasets

Following [Lyu et al., 2020; Fang et al., 2022], we employ two publicly available multi-modal MR image datasets, IXI and BraTS2018 [Menze et al., 2015] to validate the effectiveness of our model. IXI dataset contains 576 multi-contrast MR volumes and we employ the PD as the reference modal to guide the reconstruction of the target T2 modal. BraTS2018 dataset contains 285 multi-contrast MR volumes and we employ the T1 as the reference modal to guide the reconstruction of the target T2 modal. We splitted the two datasets patientwisely into a ratio of 7:1:2 for training/validation/testing. For each volume, we only select the middle 100 slices for our experiments. We use cartesian sampling mask and the k-space center cropping [Lyu et al., 2020] mask to get the undersampled and LR images, respectively. The acceleration ratio for reconstruction and the scale factor for SR are set at 4.

4.2 Implementation Details

Training Details. Our proposed MC-CDic model is implemented in PyTorch with an NVIDIA RTX3090 GPU. We adopt the Adam optimizer with a batch size of 6. The learning rate is set to 1 × 10-4 and the models are trained 50 epochs. We employ the peak-to-noise-ratio (PSNR), structural similarity index (SSIM), and root mean squared error (RMSE) to evaluate the model performance. The higher PSNR and SSIM, the lower RMSE represent the better result.

Model Details. In our large model MC-CDic-L, the iteration stage T is set to 4. The levels of the multi-scale dictionaries are set to 3. From the first level to the third level, the channel numbers of the features are set to 64, 96, and 128. The numbers of the resblocks for the multi-scale proximal network are set to 5, 3, and 1. In our small model MC-CDicS, we use the single-scale dictionary and other settings are the same with MC-CDic-L.

4.3 Comparison with State-of-the-arts

Evaluations on Guided-SR. Guided-SR aims to reconstruct the HR image from its LR image with the aid of another reference HR image. We compare our methods with various single image SR approaches, including VDSR [Kim et al., 2016], CSN [Zhao et al., 2019], RCAN [Zhang et al., 2018] and SwinIR [Liang et al., 2021], and various multicontrast guided-SR methods, including MCSR [Lyu et al., 2020], CUNet [Deng and Dragotti, 2020], MASA [Lu et al., 2021], MINet [Feng et al., 2021] and MCMRSR [Li et al., 2022]. The qualitative results of these methods are shown in Table 1. Among these methods, our MC-CDic-L model gets the best performance. It outperforms MINet by 0.92dB and 1.29dB in PSNR on the IXI and BrainTS datasets, respectively. Our MC-CDic-S only has 1.9M parameters, but it outperforms most existing SOTA methods. We randomly select two images from the two test sets and visualize the error maps of different methods in Figure 4. As one can see, the results of our methods have fewer errors than other methods. In Figure 5, we visualize the decomposed components of our MC-CDic model on the guided-SR task. The common feature represents the learned complementary information from the reference HR image to the target LR image. By decomposing the multi-contrast images, our method can only focus on the details (red arrows) related to the target image and avoid the interference of inconsistent information (purple arrows).

Evaluations on Guided Reconstruction. Guided reconstruction aims to reconstruct a clear image from its undersampled image with the aid of another auxiliary fullysampled image. Here we compare our MC-CDic model with single-contrast reconstruction methods, including zero filling (ZF), UNet [Ronneberger et al., 2015], MUSC [Liu et al., 2022] and Restormer [Zamir et al., 2022], and multi-contrast guided-reconstruction methods, including MDUNet [Xiang et al., 2018], MTrans [Feng et al., 2022] and Restormer∗ [Zamir et al., 2022]. Note that Restormer∗ employs multicontrast images as input to restore the target contrast images. Quantitative results are shown in Table 2. From the table, we can find that our MC-CDic-L achieves the best performance on the two testing sets. Restormer∗ outperforms our MC-CDic-S, but it has far more parameters than our model. Visual comparisons are shown in Figure 6, our method can restore more anatomical details than other methods. It fully demonstrates the superiority of our method.

4.4 Ablation Analysis

The Effect of the Number of Stages T. To show how the number of iteration stages T affects the reconstruction performance, we have compared the proposed MC-CDic method with different stages. Figure 7 shows the average PSNR and SSIM performance for the guided-SR task and guidedreconstruction task on the IXI validation set when T ∈ [1, 5]. It can be observed that more stages lead to improved reconstruction performance. Specifically, from T = 1 to T = 4, the psnr improves 0.6dB and 1.4dB in the guided-SR and guided-reconstruction tasks, respectively. To balance the performance and the computational complexity, we manually set T = 4 in our final model.

The Effect of the Multi-Scale Dictionaries. We compare the performance of the proposed MC-CDic model under different levels of multi-scale dictionaries in Figure 8. From the figure, we can find that the performance increases with the feature levels. From the trend lines, we can also see that multi-scale features are more important for the reconstruction task. Specifically, three levels of feature improve the PSNR by about 0.5dB and 0.9dB compared with the single-level features on the SR and reconstruction tasks, respectively. More feature levels may achieve better performance, however, the number of parameters and computations will also increase. To balance the performance and the computational complexity, we set the feature level as 1 in our MC-CDic-S and set the feature level as 3 in our MC-CDic-L.

The Effect of the Proximal Network. In this paper, we employ a proximal network to learn the prior information in a data-driven manner. To validate its effectiveness, we compare the performance of the proposed MC-CDic model with and without the proximal network. Note that when we remove the proximal network, we use the soft-thresholding operator to constrain the learned features. The performance of these two models is shown in Table 3. For the guided-SR task, the proximal network improves the performance by 0.43dB PSNR on the IXI-T2 dataset. For the guided-reconstruction task, the proximal network improves the performance by 1.06dB PSNR on the IXI-T2 dataset. It fully demonstrates the effectiveness of the proposed proximal network.

5 Concluding Remarks

We proposed a multi-scale convolutional dictionary model (MC-CDic) for multi-contrast MRI reconstruction and SR. Different from existing DCNN-based methods that manually designed fusion rules to fuse the multi-contrast information, our MC-CDic model was constructed under the guidance of the optimization algorithm with a well-designed data fidelity term. A proximal gradient algorithm was employed to solve the model and the iterative solutions were unfolded into a deep network. Especially, each module in the network could one-to-one correspond to an operation of a solution, thus our model was interpretable. We have tested our MC-CDic model on multi-contrast MRI reconstruction and SR tasks, experimental results demonstrate the superior performance of our method over the existing SOTA methods. In the future, we will explore more applications of our model on multi-modal tasks, such as RGB-guided depth map completion and RGBguided depth map SR.

Acknowledgements

This work was supported by the National Key R&D Program of China (2022ZD0161800), the NSFC-RGC (61961160734), the Shanghai Rising-Star Program (21QA1402500), and the Shanghai Municipal Commission of Economy and Informatization Program (2021-GZLRGZN-01028).