SelfCoLearn: Self-supervised collaborative
learning for accelerating dynamic MR imaging

The goal of dynamic MR imaging is to estimate dynamic MR image sequences $\mathbf{x}\in\mathbb{C}^{N}$ from undersampled measurements $\mathbf{y}\in\mathbb{C}^{M}(M\ll N)$ in k-space. $N=N_{h}N_{W}T$ is a vector. $N_{h}$ and $N_{W}$ are height and width of the frame respectively. $T$ represents the number of frames in each sequence. Thus, the imaging model can be described as follows:

where $\mathrm{e}\in\mathbb{C}^{M}$ is noise and $\mathbf{A}=\mathbf{PF}$ is an undersampled Fourier encoding operator, $\mathbf{F}$ is 2D Fourier transform to each frame in dynamic image sequence and $\mathbf{P}$ is the undersampled mask for each frame. In general, the reconstruction problem is formulated as an unconstrained optimization problem to explore the prior knowledge:

where $\mathcal{R}(\mathbf{x})$ represents a prior regularization item on $\mathbf{x}$, and $\lambda$ is the weight of the regularization. $\frac{1}{2}\|\mathbf{A}\mathbf{x}-\mathrm{y}\|_{2}^{2}$ is the data fidelity item, which ensures the reconstruction result to be consistent with the original undersampled measurements.

For fully-supervised deep learning methods, it typically uses a CNN $f_{CNN}\left(\mathbf{y}\mid\theta\right)$, as a regularization term $\mathcal{R}(\mathbf{x})$, by learning the mapping between undersampled/ corrupted data and their corresponding fully sampled data with parameters $\theta$. Its mathematical description can be given as:

where $S$ is the number of dynamic image sequences in the training set. $\mathbf{x}_{i}^{ref}$ is the fully sampled data of the subject data $i$. $\mathcal{L}(\cdot)$ denotes the loss function between the predicted reconstruction output and the fully sampled reference data, which typically adopts the $l_{1}-$norm or $l_{2}-$norm.

In our work, we propose a simple but effective self-supervised training framework for dynamic MR imaging, whose paradigm is shown in Fig. 1. Our framework trains two independent reconstruction networks simultaneously, which have different inputs and different weight parameters. The backbone network can adopt either data-driven network architecture or the iterative un-rolled network, such as CRNN [26], k-t NEXT [24] and SLR-Net [29] et al. Based on the consistency between two network predictions, the network provides complementary information for the to-be-reconstructed dynamic MR images in its peer partner, which is an additional regularization compared to the existing unsupervised methods [32]. The two networks will finally realize consistent reconstruction in the training process. Specifically, given a raw undersampled k-space data sequence $\Omega=\left\{\mathbf{y}_{\Omega}^{t}\right\}_{t=1}^{T}$, we reundersampled the original k-space data $\mathbf{y}_{\Omega}^{t}$ to construct a partial data points sequence $\left\{\mathbf{y}_{\mathrm{u}}^{t}\right\}_{t=1}^{T}$:

where $t$ is the index of sequence, $u$ denotes the index of the two training sequences and $P_{\mathrm{u}}^{t}$ is the undersampled mask for frame $t$. In order to make full use of all data points in $\mathbf{y}_{\Omega}^{t}$ to learn representation, and ensure each network can provide complementary information for the to-be-reconstructed dynamic MR images in its peer network, these training sequences are generated adhere to the following data augmented principles: (1) The union of data points in two training sequences must be equal to the data $\mathbf{y}_{\Omega}^{t}$, i.e., $\mathbf{y}_{\Omega}^{t}=\mathbf{y}_{\Theta}^{t}\cup\mathbf{y}_{\Lambda}^{t}$. (2) The data points in two training sequences should be different, i.e., $\mathbf{y}_{\Theta}^{t}\neq\mathbf{y}_{\Lambda}^{t}$. (3) The training sequences should include most of the low frequency data points and part of the high frequency data points. Following these principles, the two training sequences contain similar data points in the low frequency region, and different points in the high frequency region. Noted that the operation of data reundersampling is only necessary during training, the reconstructed images can be inferred from the test data directly.

1) Data-driven dynamic MR imaging: In the data-driven settings, the common practice is to decouple Eq. 2 into a regularization term and a data fidelity term via utilizing the variable splitting technique [25, 26]. By introducing an auxiliary variable $\mathbf{z}=\mathbf{x}$, Eq.2 can be re-formulated as the penalty function:

where $\mu$ is a penalty parameter. Eq.5 can then be solved iteratively via alternating minimization over $\mathbf{z}$ and $\mathbf{x}$:

where $n\in\left\{0,1,2,\ldots,N-1\right\}$ is the $n$th iteration, $\mathbf{x}^{0}$ is the zero-filling image transformed from original undersampled measurement, $\mathbf{z}^{n}$ denotes the intermediate reconstruction sequence, and $\mathbf{x}^{n}$ denotes the final reconstruction sequence at each iteration. In Eq.7, the operation on the intermediate reconstruction sequence $\mathbf{z}^{n}$ is a data consistency step, which uses the original sampled k-space data points to replace the corresponding data points in the reconstructed k-space data [25]. The iterative optimization process in Eq.6 and Eq.7 is unrolled into a neural network.

CRNN [26] is a typical data-driven method that integrates data consistency in k-space. A single iteration of CRNN can be illustrated as the following process:

where $\mathbf{x}_{rnn}^{(n)}$ is the intermediate reconstruction sequence analogous to $\mathbf{z}^{n}$ in Eq.6, and $\mathbf{x}_{rec}^{(n)}$ denotes the final predicted result at each iteration analogous to $\mathbf{x}^{n}$ in Eq.7. The regularization subproblem in Eq.6 is solved by using a convolutional recurrent neural network. The data consistency subproblem in Eq.7 is treated as a data consistency network layer. The unrolled architecture of CRNN is shown in Fig. 2. More details of CRNN layers can be found in [26].

2) Un-rolled dynamic MR imaging: Another widely-used strategy is un-rolled dynamic MR imaging, which constructs CNNs according to the iterations of traditional optimization algorithms. Different optimization algorithms lead to different network architectures. SLR-Net, which formulates sparse and low rank priors as regularized terms in an optimization algorithm [29], is a typical example of un-rolled method. In SLR-Net, by introducing an auxiliary variable $\mathbf{M}$, Eq.2 can be decoupled as the fidelity term, sparse regularization term and the low rank regularization term:

where $D$ is a sparse transform in a certain sparse domain. $M=R\mathbf{x}$ is a matrix (with size ($N_{\mathbf{h}}\times N_{\mathbf{w}}$, $T$)), in which each column corresponds to one frame in dynamic MR image sequence. $R$ is a reshaping operator. $\|M\|_{*}$ is the nuclear norm. Previous works have proven that nuclear norm minimization is effective in low-rank matrix recovery [37]. More details of the iterative process in SLR-Net can be found in [29].

We have designed a co-training loss to promote accurate dynamic MR image reconstruction in a self-supervised manner. The core idea of the co-training loss is to enforce the consistency not only between the reconstruction results and the original undersampled k-space data, but also between two network predictions. Compared with existing methods with single network, the consistency between two network predictions is an additional regularization, which guides the dual-network to learn more correct information. Specifically, the co-training loss in SelfCoLearn, including an undersampled consistency loss term and a contrastive consistency loss term, is calculated to optimize the whole framework.

Let $f_{SelfCoLearn}\left(\mathbf{y}_{\Omega}^{t}\right)$ denotes SelfCoLearn, $\mathbf{y}_{\Omega}^{t}$ is the original undersampled k-space data. During training, two training sequences $\mathbf{y}_{\Theta}^{t}$ and $\mathbf{y}_{\Lambda}^{t}$ are generated from $\mathbf{y}_{\Omega}^{t}$ following the data augmented principles in section II-B as follows:

where $\mathbf{P}_{\Theta}^{t}$ and $\mathbf{P}_{\Lambda}^{t}$ are reundersampled mask for $\mathbf{y}_{\Omega}^{t}$. The undersampled consistency loss is mainly referred to the actually sampled k-space points in $\mathbf{y}_{\Omega}^{t}$, which ensures the corresponding sampled points in network prediction consistent with actually sampled k-space points in $\mathbf{y}_{\Omega}^{t}$. The corresponding sampled points $\mathrm{y}_{\Theta\rightarrow\Omega}^{t}$ and $\mathrm{y}_{\Lambda\rightarrow\Omega}^{t}$ in these two network predictions can be written as:

where k-space data $f\left(\mathbf{y}_{\Theta}^{t}\right)$ and $f\left(\mathbf{y}_{\Lambda}^{t}\right)$ are transformed from the predicted image sequences of two networks, respectively. $\mathbf{P}^{t}$ is the undersampled mask, which is applied to generate the original undersampled k-space data $\mathbf{y}_{\Omega}^{t}$ from the fully sampled data. Thus, the Undersampled Consistency loss term is used to calculate the mean-square-error between the actually sampled k-space points in $\mathbf{y}_{\Omega}^{t}$ and that in each network prediction as follows:

In the ideal case, when different reundersampled k-space data from the same data are feed into two networks, the network predictions shall approximate the fully sampled reference data with the network optimizations. However, when the fully sampled reference data are absent, these two networks may produce different predicted results only with the undersampled consistency loss, and result in different reconstruction performances. As mentioned above, a contrastive consistency loss is defined to compute the mean-square-error between two network predictions with different reundersampling inputs from the same data. Specially, our contrastive consistency loss term is mainly referred to the corresponding points in network predictions to unsampled k-space points in $\mathbf{y}_{\Omega}^{t}$. These corresponding points $\bar{\mathbf{y}}_{\Theta\rightarrow\Omega}^{t}$ and $\bar{\mathbf{y}}_{\Lambda\rightarrow\Omega}^{t}$ in two network predictions $f\left(\mathbf{y}_{\Theta}^{t}\right)$ and $f\left(\mathbf{y}_{\Lambda}^{t}\right)$ can be written as:

therefore, the Contrastive Consistency loss term is formulated as:

Combining the two loss terms, our final co-training loss function can be written as:

where $\gamma$ is used to balance the weight parameter of undersampled consistency loss and contrastive consistency loss. During testing phase, undersampled data sequence is set as input of collaborative network-1 to obtain the final reconstruction result.

1) Dataset: T1-weighted FLASH sequence is utilized to collect fully sampled cardiac data from 101 volunteers on a 3T scanner. All in vivo experiments have been approved by the Institutional Review Board (IRB) of Shenzhen Institutes of Advanced Technology, and written informed consent have been obtained from all volunteers. Each scan acquires a single slice from the volunteer with 25 temporal frames. The following parameters were used for the FLASH sequence: FOV 330×330 mm, acquisition matrix 192×192, slice thickness = 6 mm, TR/TE = 50 ms/3 ms, and 24 receiving coils. The raw multi-coil data of each frame was combined by the adaptive coil combine method [38] to produce a single-channel complex-valued reconstruction image. Then, the complex-valued images were transformed to k-space data, which simulate a fully sampled single-coil data acquisition. Training of neural networks requires a large amount of data. To this end, we use data augmentation to enlarge the data set. Specifically, we translated the original images along x, y, and t directions, and the translation step size is 128×128×14, and the stride along the three directions are 12, 12, and 3, respectively. Finally, our dataset consists of 6214 complex-valued cardiac MR data sequences with size 128×128×14. 5950 cardiac MR data sequences were randomly selected as the training dataset, 50 cardiac sequences were used as validation dataset, and the remaining sequences were used for testing.

2) Reunderampling K-space Data Augmentation: In the proposed method, the fully sampled data are only used to generate the original undersampled k-space data $\mathbf{y}_{\Omega}^{t}$ with a 2D random retrospective undersampled mask $\mathbf{P}^{t}$. Following the principles of training data augmentation in section II-B, $\mathbf{y}_{\Omega}^{t}$ is augmented to two training sequences $\mathbf{y}_{\Theta}^{t}$ and $\mathbf{y}_{\Lambda}^{t}$ with two 2D random reundersampled masks $\mathbf{P}_{\Theta}^{t}$ and $\mathbf{P}_{\Lambda}^{t}$. $\mathbf{P}_{\Theta}^{t}$ with 2-fold acceleration is used for collaborative network-1, and $\mathbf{P}_{\Lambda}^{t}$, which is combines the complementary set of $\mathbf{P}_{\Theta}^{t}$ with some low-frequency data points of $\mathbf{P}^{t}$, is used for collaborative network-2.

3) Evaluation Metrics: To evaluate the reconstruction performance, the mean-square-error (MSE), peak-signal-to-noise ratio (PSNR), and structural similarity index (SSIM) [39] are calculated as follows:

where $\mathrm{Rec}$ is the reconstructed image sequence, and $\mathrm{Ref}$ represents the reference image sequence. $\mathrm{MAX}_{Ref}$ is the maximum possible value in the image. $\mu_{Ref}$ and $\mu_{Rec}$ are the averaged intensity values of the corresponding images. $\sigma_{Ref}$ and $\sigma_{Rec}$ are the variances. $c_{1}$ and $c_{2}$ are adjustable constants. $\sigma_{Ref,Rec}$ is the covariance. The SSIM index is a multiplicative combination of the luminance term, the contrast term, and the structural term (details shown in [39]).

4) Model Configuration and Implementation Details: The proposed framework is flexible to be integrated with both data-driven and iterative un-rolled networks. Most of our experiments adopt CRNN as the backbone network. In detail, the network is composed of a bidirectional CRNN layer, three CRNN layers, a 2D CNN layer, a residual connection sums output with input and a DC layer. The nonlinear activation function utilized is the Rectified linear unit (ReLU). For the bidirectional CRNN and CRNN layer, the number of convolution filters is set to $n_{f}=64$ with a kernel size of $k=3$. The 2D CNN contains one convolutional layer with $k=3$ and $n_{f}=2$. We use $stride=1$ and the padding is set to half of the filter size (rounded down). The DC layer is followed by the 2D CNN layer, which forces the sampled data points in predicted k-space data to be consistent with that in the inputs.

For model training, the number of iteration step is set to $N=5$. The batch size is set to 1. All training data and test data are normalized to [0,1]. The SelfCoLearn framework with CRNN and k-t NEXT is implemented in PyTorch 1.8.1, and that with SLR-Net is implemented in Tensorflow2.2.0. The experiments are performed on a GPU server with an Nvidia Titan Xp Graphics Processing Unit (GPU, 12GB memory). The model is trained by Adam optimizer [40] with parameters $\beta_{1}=0.5$ and $\beta_{2}=0.999$. The weight parameter $\gamma$ in co-training loss is set to 0.01. The learning rate is set to $10^{-4}$. It takes 42 hours to train SelfCoLearn with CRNN for 40 epochs and each cardiac MR data sequence takes roughly 0.5 second to get the reconstructed result.

To demonstrate the superiority of the proposed SelfCoLearn, we compared it with two self-supervised methods, SS-DCCNN and SS-CRNN, at different acceleration factors. It is worth noting that the state-of-the-art self-supervised method SSDU [34] was developed for static MR imaging. Literature [32] adopted a similar self-supervised training manner as SSDU for dynamic MR imaging. It evaluated several backbone architectures for dynamic MR imaging including DCCNN and CRNN, whereas SSDU adopted ResNet as backbone network. We choose two self-supervised learning methods SS-DCCNN and SS-CRNN [32] for comparison. In this experiment, the proposed SelfCoLearn selects the CRNN as the backbone network. For fair comparison, all methods are carefully tuned to obtain the best performance on the current dataset.

Fig. 3 plots the reconstruction results of different methods at 8-fold acceleration, 12-fold acceleration and 14-fold acceleration, respectively. The first three rows show diastolic reconstruction at the tenth frame of image sequence, and the following three rows show systolic reconstruction at the fifth frame of image sequence. The first row shows, from left to right, the fully sampled reference image and the reconstruction results of respective methods (display range [0, 1]). The second row shows the enlarged views of the heart regions. The third row shows the corresponding error maps (display range [0, 0.2]). The y-t views at $\mathbf{x}=40$ are shown in the seventh row. The corresponding error maps of y-t views are shown in the last row. From the visualization results, the proposed SelfCoLearn generates better reconstruction results than the two self-supervised methods, SS-DCCNN and SS-CRNN, at all acceleration factors. The reconstruction images of SelfCoLearn show finer structural details and more precise heart borders with fewer artifacts, while the SS-DCCNN and SS-CRNN exhibit artifacts within the heart border and chambers. The error maps of SelfCoLearn also indicate minor reconstruction errors, especially at the borders of heart chambers.

The quantitative results of these self-supervised methods are listed in Table 1. Similar conclusions can be drawn that the proposed SelfCoLearn achieves better quantitative performance than the two existing self-supervised learning methods. Therefore, our collaborative learning strategy can effectively capture essential and inherent representations directly from undersampled k-space data.

Fig. 4 give the box plots displaying the median and interquartile range (25th-75th percentile) of the reconstruction results of different self-supervised methods across all test cardiac cine data at 8-fold acceleration, 12-fold acceleration and 14-fold acceleration, respectively. For all dynamic cine sequences, SelfCoLearn outperforms the two self-supervised learning methods (SS-DCCNN and SS-CRNN) at all three acceleration factors.

We further compare our SelfCoLearn with different supervised methods, including supervised U-Net and supervised CRNN [26], at different acceleration factors. Fig. 5 plots the reconstruction results of different methods at 8-fold acceleration, 12-fold acceleration and 14-fold acceleration, respectively. From the visualization results, the proposed SelfCoLearn restores more precise anatomical details of heart regions than supervised U-Net at all acceleration factors, while supervised U-Net fails to recover some details in the heart chambers. The error maps of SelfCoLearn also indicate minor reconstruction errors than those of supervised U-Net.

In addition, SelfCoLearn generates comparable reconstruction results to those of supervised CRNN at low acceleration factors. From the enlarged heart region plots, the images reconstructed by SelfCoLearn are as clear as those generated by supervised CRNN. At higher acceleration factors, such as 14-fold acceleration, the reconstructed images of SelfCoLearn become slightly blur. Nevertheless, most of the structural details in the heart regions are still successfully restored by SelfCoLearn. Quantitative results at all acceleration factors (Table 2) also show the promising results of SelfCoLearn. Therefore, it can be concluded that the proposed SelfCoLearn can achieve comparable reconstruction performance with fully-supervised methods via self-supervised collaborative learning for accelerating dynamic MR imaging.

In this section, we explore the reconstruction results of the proposed self-supervised learning strategy with different backbone networks for dynamic MR imaging. Experiments are conducted using SLR-Net [29], k-t NEXT [24], and CRNN [26] at 8-fold acceleration. The reconstruction results with different backbone networks can be found in Fig. 6 and Table 3. Compared with SS-CRNN [8], the proposed self-supervised learning strategy can achieve better results regardless of the utilized backbone network. Among the three different backbone networks, SLR-Net generates worse results than k-t NEXT and CRNN. The reason for this phenomenon may be that SLR-Net needs to learn a singular value threshold, and the absence of fully sampled reference data causes the learned singular value threshold is suboptimal. However, the proposed self-supervised learning strategy with SLR-Net still obtain acceptable reconstruction results. The qualitative results in Fig. 6. clearly show that the proposed SelfCoLearn can better restore the structural details and achieve clearer reconstructed MR images (especially in the heart regions around the red and yellow arrows) than SS-CRNN. The quantitative results also indicate more accurate reconstructions achieved by the proposed SelfCoLearn. These results confirm that our proposed self-supervised learning framework is flexible regarding the adopted backbone network, and it can achieve promising reconstruction results with both data-driven and iterative un-rolled networks for dynamic MR imaging.

In this section, we investigate the utility of the designed the co-training loss function. The backbone network in these experiments adopts CRNN. Different training strategies at 8-fold acceleration are utilized: Strategy B-I – a single reconstruction network is trained in self-supervised manner. Only the loss function between the output $f\left(\mathbf{y}_{\Theta}^{t}\right)$ of network and $\mathbf{y}_{\Lambda}^{t}$ is used to train the network. This training strategy is similar to that of SSDU. Strategy B-II – a strategy similar to B-I but the loss function here is calculated between the output $f\left(\mathbf{y}_{\Theta}^{t}\right)$ of network and original undersampled k-space data $\mathbf{y}_{\Omega}^{t}$. SelfCoLearn – two networks are trained collaboratively with $\mathcal{L}_{UC}$ and $\mathcal{L}_{CC}$, and the two collaborative networks adopt the same backbone network as that in strategy B-I. Reconstruction images of methods utilizing the different training strategies are plotted in Fig. 7. Quantitative results are listed in Table 4. From both qualitative and quantitative results, we can observe that SelfCoLearn (training two networks collaboratively with both loss terms) achieves the best performance (especially in the heart regions around the red and yellow arrows). In particular, the contrastive consistency loss term results in a large reconstruction performance improvement. For example, PSNR is improved from 31.04dB (Strategy B-II) to 37.27 dB (SelfCoLearn).

In this section, we inspect the effects of loss functions. The backbone network in these experiments adopts CRNN. Reconstruction results at 8-fold acceleration are given in Fig. 8 and Table 5. Three strategies utilizing different loss function settings are investigated. In Strategy C-I, two networks are trained collaboratively with $\mathcal{L}_{UC}$ and $\mathcal{L}_{CC}$. in which $\mathcal{L}_{UC}$ is calculated in the x-t domain, and $\mathcal{L}_{CC}$ is calculated in the k-space domain. In Strategy C-II, both $\mathcal{L}_{UC}$ and $\mathcal{L}_{CC}$ are calculated in the x-t domain. In Strategy C-III, both $\mathcal{L}_{UC}$ and $\mathcal{L}_{CC}$ are calculated in the k-space domain. From both qualitative and quantitative results, we can observe that the influence of utilizing different loss function settings on the reconstruction performance is insignificant. All the other experiments in this work adopt the loss function setting of strategy C-III.