INTRODUCTION

Myocardial parameter mapping is an important technique in cardiovascular MR imaging (CMR) (1). There is a growing emphasis on characterizing multiple relaxation times by acquiring multi-contrast (MC) images (2,3). 3D-CMR offers comprehensive information about the entire heart compared to 2D imaging with limited coverage, but it suffers from the long scan time. Deep learning (DL) has emerged as a successful tool in MR reconstruction. However, most DL-based reconstruction methods rely on supervised learning that requires extensive fully sampled datasets for training. Acquiring fully sampled k-space data is challenging for 3D-MC-CMR since the scan time is extremely lengthy. This study endeavors to develop a self-supervised reconstruction method, referred to as SSJDM, for 3D-MC-CMR. It employs a score-based model with a Bayesian convolutional neural network, obviating the need for fully sampled training data.

METHODS

MRI reconstruction model: Let $$$\boldsymbol Y$$$ and $$$ \boldsymbol X$$$ denote the undersampled k-space data and the corresponding MR image, respectively. The reconstruction model to recover $$$\boldsymbol X$$$ is:
\begin{equation} \begin{gathered}\underset{\boldsymbol{X}}{\operatorname{\arg \min}} \frac{1}{2}\left \| \boldsymbol{AX}-\boldsymbol{Y} \right\|_{F}^{2} +\lambda R(\boldsymbol{X}) \end{gathered} \end{equation}
where $$$\boldsymbol A$$$ represents the encoding operator given by $$$\boldsymbol A=\boldsymbol{MFS}$$$, $$$\boldsymbol M$$$ is the undersampling operator, $$$ \boldsymbol F$$$ is the Fourier transform operator, and $$$\boldsymbol S$$$ denotes the coil sensitivity map, $$$R(\boldsymbol{X})$$$ denotes a combination of regularizations. From a Bayesian perspective, $$$R(\boldsymbol{X})$$$ can be considered as the prior model of the data, denoted as $$$p(\boldsymbol{X})$$$ (4).

Self-supervised Bayesian reconstruction network (BCNN): A secondary sub-sampling mask is applied to the measurements, resulting in $$${\boldsymbol{Y}^{\prime}} = {\boldsymbol{M}^{\prime}} \boldsymbol{Y}$$$. A BCNN can be trained using the sub-sampled data $$${\boldsymbol{Y}^{\prime}}$$$ and the zero-filling image $$${\boldsymbol{X}^{\prime}}$$$ of $$$\boldsymbol{Y}$$$with:
\begin{equation}{\boldsymbol{X}^{\prime}} = \boldsymbol{f}_{\boldsymbol{\theta}}( {\boldsymbol{Y}^{\prime}} ) +\boldsymbol{n}_1 \end{equation}
where $$$\boldsymbol{f_{\theta}}$$$ represents a mapping parameterized by $$$\boldsymbol{\theta}$$$ which follows a distribution $$$p(\boldsymbol{\theta})$$$. Then $$$\boldsymbol{f_{\theta}}$$$ is utilized to convert the measured $$$\boldsymbol Y$$$ to its corresponding image $$$\boldsymbol X$$$:
\begin{equation}\boldsymbol{X}= \boldsymbol{f}_{\boldsymbol{\theta}}(\boldsymbol{Y}) +\boldsymbol{n}_2 \end{equation}
where $$$\boldsymbol{n}_1$$$ and $$$\boldsymbol{n}_2$$$ represent Gaussian noise with scales $$$\gamma_1$$$ and $$$\gamma_2$$$. Assume that $$$p(\boldsymbol{n}_1)$$$ $$$\sim$$$ $$$\prod_{t } \exp \left(\frac{-{n}_t^2}{2 {\gamma_1}^2}\right)$$$, and the prior $$$p(\boldsymbol{\theta})$$$ $$$ \sim$$$ $$$\prod_s \exp \left(\frac{-{\theta}_s^2}{2 {\bar{\sigma}}^2}\right)$$$, an approximate distribution $$$q(\boldsymbol{\theta} )$$$ can be obtained by minimizing the following the Kullback-Leibler divergence:
\begin{equation}\begin{aligned}& \min _{\boldsymbol{\mu_{\theta}}, \boldsymbol{\sigma_{\theta}}} \frac{1}{2{\gamma_1}^2}\mathbb{E}_{q(\boldsymbol{\theta} \mid \boldsymbol{\mu_{\theta}}, \boldsymbol{\sigma_{\theta}})} \sum_{i=1}^N \left\|\boldsymbol{f_{\theta}}({\boldsymbol{X}^{\prime}}_i) -\boldsymbol{Y}_i\right\|_2^2\\&+\frac{1}{2\bar{\sigma}^2} \left(\|\boldsymbol{\mu_{\theta}}\|_2^2+\|\boldsymbol{\sigma_{\theta}}\|_2^2\right)- \sum_s \log \frac{\sigma_{s}}{\bar \sigma} +\text {const.}\end{aligned}\end{equation}
MC image reconstruction via score-based model: Utilizing the above estimated $$$q(\boldsymbol{\theta} )$$$, the score matching technique (5) is employed to approximate $$$\nabla \log_{\boldsymbol X} p(\boldsymbol{X})$$$ using a score-matching network parameterized by $$$\boldsymbol{\phi}$$$ in a way that $$$\boldsymbol{s_\phi}(\boldsymbol{X};\boldsymbol \theta)=\nabla \log q(\boldsymbol{X};\boldsymbol \theta)$$$ by training the joint score function with:
\begin{equation}\frac{1}{2 L} \sum_{i=1}^L \mathbb{E}_{p(\boldsymbol{Y}) q(\boldsymbol{\theta})} \mathbb{E}_{p_{\varepsilon_i}(\tilde{\boldsymbol{X}} | \boldsymbol{Y},\boldsymbol{\theta})}\left [\left\| \varepsilon_i \boldsymbol{s}_{\boldsymbol{\phi}}(\tilde{\boldsymbol{X}}, \varepsilon_i)+\frac{\tilde{\boldsymbol{X}}-\boldsymbol{f}_{\boldsymbol{\theta}}(\boldsymbol{Y})}{\varepsilon_i}\right\| ^2\right]\end{equation}
After the score function is estimated, the Langevin MCMC sampling can be applied to reconstruct MC images via the following formula:
\begin{equation}\begin{aligned}\tilde{\boldsymbol{X}}_{i+1} & =\tilde{\boldsymbol{X}}_i+\frac{\eta_i}{2} \nabla \log p(\boldsymbol{\tilde{X}}_i | \boldsymbol{Y})+\sqrt{\eta_i} \boldsymbol{z}_i \\& =\tilde{\boldsymbol{X}}_i+\frac{\eta_i}{2}(\boldsymbol{s}_\phi(\tilde{{\boldsymbol X}}_i, \varepsilon_i)+\frac{\boldsymbol{A}^H(\boldsymbol{A} \tilde {\boldsymbol{X}}_i-\boldsymbol{Y})}{{\gamma_2}^2+\varepsilon_i^2})+\sqrt{\eta_i} \boldsymbol{z}_i\end{aligned}\end{equation}
where $$$\eta_i$$$ serves as the step size.

Figure 1 shows the proposed SSJDM framework. Figure 2 shows the structure of BCNN. The NCSNv2 framework in (6) was used as the score matching network. We utilized a 3D simultaneous cardiac T1 and T1ρ mapping sequence (Figure 1(c)) on a 3T MR scanner (uMR 790, United Imaging Healthcare, Shanghai, China). The 3D-MC-CMR dataset using this sequence was prospectively undersampled with an acceleration rate of R = 6. The training set consisted of 5796 slices from 46 volunteers, while the test set included 1512 slices from 12 volunteers. To evaluate SSJDM, the test data were retrospectively undersampled using a 2D Poisson-disc sampling pattern with R = 11 and 14. The SSJDM was used to reconstruct the MC images, and T1 and T1ρ maps were estimated using the dictionary matching method (7). The results were compared with those obtained using traditional regularized reconstruction (PROST (8)), a self-supervised DL method (SSDU (9)), and BCNN.

RESULTS and DISCUSSION

Figure 3 shows the reconstructed MC images using the four methods at R = 11. Blurring artifacts of PROST become noticeable in scenarios with higher acceleration rates. BCNN preserves more detailed information than SSDU (indicated by the yellow arrows). Images of SSJDM still exhibit sharp boundaries and high texture fidelity. At an even higher acceleration rate of R = 14 (shown in Figure 4), the image quality of SSJDM degrades a little, while other methods exhibit severe blurring artifacts. Figure 5 shows the T1 and T1ρ maps estimated from reconstructions using different methods at R = 11. Similar conclusions can be drawn.

CONCLUSION

This study shows the feasibility of score-based models for reconstructing 3D-MC-CMR images from highly undersampled data in simultaneous whole-heart T1 and T1ρ mapping.

Acknowledgements

This work is supported in part by the National Natural Science Foundation of China under grant no. 62201561, National Key R&D Program of China under grant no. 2021YFF0501402, and the Guangdong Basic and Applied Basic Research Foundation under grant no. 2021A1515110540.