Training Physics-Driven Deep Learning Reconstruction without Raw Data Access for Equitable Fast MRI

MRI raw data is acquired in the frequency domain of the image, referred to as k-space. For fast MRI, data is acquired in the sub-Nyquist regime by undersampling the acquisition in k-space. In this case, the forward acquisition model relating the image $\bf{x}\in{\mathbb{C}}^{n}$ to these raw MRI data (or k-space) measurements is given as:

$$\mathbf{y}_{\Omega}=\mathbf{E}_{\Omega}\bf{x}+\bf{n},$$

(1)

where $\mathbf{y}_{\Omega}$ denotes the acquired k-space data corresponding to the undersampling pattern $\Omega$ with $|\Omega|=m<n$. $\mathbf{E}_{\Omega}$ denotes the multi-coil encoding operator that includes information from $n_{c}$ receiver coils, each of which are sensitive to a different part of the image (Hamilton et al., 2017). When the acceleration rate $R=n/m$ is less than $n_{c}$, this system of equations is over-determined due to the redundancies among the receiver coils. Parallel imaging uses these redundancies to solve the maximum likelihood estimation problem under i.i.d. Gaussian noise (Pruessmann et al., 1999):

$${\bf x}_{\textrm{PI}}=\arg\min_{\bf x}\|\mathbf{y}_{\Omega}-\mathbf{E}_{\Omega}\mathbf{x}\|_{2}^{2}=({\bf E}_{\Omega}^{H}{\bf E}_{\Omega})^{-1}{\bf E}_{\Omega}^{H}{\bf y}_{\Omega}.$$

(2)

Numerically, this can be solved directly for certain undersampling patterns (Pruessmann et al., 1999) or more broadly iteratively using conjugate gradient (CG) (Pruessmann et al., 2001). Using the equivalence of multiplication in image domain and convolutions in k-space (Uecker et al., 2014), it can also be solved as an interpolation problem in k-space (Griswold et al., 2002). Parallel imaging remains the most clinically used acceleration method for MRI, with some MR systems using the image-based reconstruction, while others utilizing the equivalent k-space interpolation.

In modern computational MRI, additional regularization is often incorporated into the objective function (Hammernik et al., 2023):

$$\arg\min_{\bf x}\|\mathbf{y}_{\Omega}-\mathbf{E}_{\Omega}\mathbf{x}\|_{2}^{2}+\mathcal{R}(\mathbf{x}),$$

(3)

where $\mathcal{R}(\cdot)$ denotes a regularizer. For instance, compressed sensing (CS) uses the idea that images should be compressible in an appropriate transform domain (Lustig et al., 2007), and uses $\mathcal{R}(\mathbf{x})=\tau||{\bf Wx}||_{1}$, where $\tau$ is the regularization weight, ${\bf W}$ is a linear sparsifying transform such as a discrete wavelet transform (DWT) and $||\cdot||_{1}$ is the $\ell_{1}$ norm.

Among different PD-DL methods (Ahmad et al., 2020; Gilton et al., 2021; Knoll et al., 2020a), unrolled networks (Monga et al., 2021) remain the highest performer in reconstruction challenges, as reported a year ago (Hammernik et al., 2023; Muckley et al., 2021). These methods unroll iterative algorithms for solving the regularized least squares objective in (3) (Fessler, 2020), such as proximal gradient descent (Schlemper et al., 2018) or variable splitting with quadratic penalty (VS-QP) (Aggarwal et al., 2019), over a fixed number of steps. VS-QP transforms (3) into 2 sub-problems:

	$$\mathbf{z}^{(i)}=\arg\min_{\bf{z}}\|\mathbf{x}^{(i-1)}-\mathbf{z}\|_{2}^{2}+\mathcal{R}(\mathbf{z}),$$		(4a)
	$$\mathbf{x}^{(i)}=\arg\min_{\bf x}\|\mathbf{y}_{\Omega}-\mathbf{E}_{\Omega}\mathbf{x}\|_{2}^{2}+\mu\|\mathbf{x}-\mathbf{z}^{(i)}\|_{2}^{2},$$		(4b)

where (4a) is the proximal operator for the regularization, implicitly solved using neural networks, while (4b) accounts for data fidelity and has a closed form solution:

$$\mathbf{x}^{(i)}=\left(\mathbf{E}^{H}_{\Omega}\mathbf{E}_{\Omega}+\mu\mathbf{I}\right)^{-1}(\mathbf{E}^{H}_{\Omega}\mathbf{y}_{\Omega}+\mu\mathbf{z}^{(i)}),$$

(5)

which can be solved by CG (Aggarwal et al., 2019). Unrolled networks are conventionally trained using supervised learning over a database, where the reference raw k-space measurements are first retrospectively undersampled to form ${\bf y}_{\Omega}$. Subsequently, the network is trained to map to the original full reference k-space or the corresponding reference image (Hammernik et al., 2018; Aggarwal et al., 2019) by minimizing:

$$\min_{\bm{\theta}}\mathbb{E}\left[\mathcal{L}\left(\mathbf{y}_{\text{ref}},\mathbf{E}_{\text{full}}(f(\mathbf{y}_{\Omega},\mathbf{E}_{\Omega};\bm{\theta}))\right)\right]$$

(6)

where $\bm{\theta}$ are the network parameters, $f(\mathbf{y}_{\Omega},\mathbf{E}_{\Omega};\bm{\theta})$ denotes the network output for inputs ${\bf y}_{\Omega}$ and ${\bf E}_{\Omega}$, $\mathbf{E}_{\text{full}}$ is the fully-sampled encoding operator, $\mathbf{y}_{\text{ref}}$ is the fully-sampled reference k-space data, and ${\cal L}(\cdot,\cdot)$ is a loss function.

Obtaining fully-sampled reference data in MRI can be infeasible due to prolonged scan durations, organ movement in acquisitions such as real-time cardiac imaging or myocardial perfusion (Rajiah et al., 2023), or signal decay in acquisitions like diffusion MRI with EPI (Uğurbil et al., 2013). To enable training of PD-DL networks without fully sampled raw MRI data, a variety of unsupervised learning methodologies have emerged (Akçakaya et al., 2022), including self-supervised learning techniques (Yaman et al., 2020; Chen et al., 2021) and generative modeling approaches (Jalal et al., 2021; Chung & Ye, 2022; Chung et al., 2023).

Self-supervised methods use a masking approach to generate supervisory labels from the undersampled data (Yaman et al., 2020; Millard & Chiew, 2023; Hu et al., 2024). A pioneering method in this field, self-supervision via data undersampling (SSDU) (Yaman et al., 2020, 2022a), involves partitioning the acquired measurement $\Omega$ into two disjoint subsets ($\Omega=\Lambda\cup\Theta$) to train the network in a self-supervised manner:

$$\min_{\bm{\theta}}\mathbb{E}\left[\mathcal{L}\left(\mathbf{y}_{\Lambda},\mathbf{E}_{\Lambda}(f(\mathbf{y}_{\Theta},\mathbf{E}_{\Theta};\bm{\theta}))\right)\right]$$

(7)

Even though these self-supervision based approaches demonstrate exceptional performance across various tasks, they lack the ability to train the model without access to undersampled raw data, as they cannot operate solely using images that are exported from the scanner.

Conversely, generative methods learn the prior distribution of the given dataset, which is then leveraged in conjunction with a log-likelihood data term during the testing phase. Although recent methods based on diffusion/score-based models have shown substantial promise, these methods require large amounts of high-quality images either reconstructed from raw data (Jalal et al., 2021; Luo et al., 2023) or as DICOMs (Chung & Ye, 2022), as well as computational resources to perform the training, both of which may not be feasible in the setups we are focused on.

We conducted a thorough evaluation of our method, assessing its performance through both qualitative and quantitative analyses, and focused on uniform/equidistant patterns which produces coherent artifacts that are more difficult to remove compared to the incoherent artifacts from random undersampling (Knoll et al., 2019). We further note that CUPID demonstrates robust performance across a wide range of $\lambda$ values, provided that $\lambda$ is chosen within a reasonable range. An ablation study on the choice of $\lambda$ is included in Sec. 4.5.

Retrospective Undersampling Setup. In our retrospective studies, we used fully-sampled multi-coil knee and brain MRI data from the fastMRI database (Knoll et al., 2020b). Knee dataset included fully-sampled coronal proton density-weighted (coronal PD) and PD with fat suppression (coronal PD-FS) data. For brain MRI, axial FLAIR (ax-FLAIR) dataset with matrix size of $320\times 320$ is used. The knee and brain MRI datasets comprised data collected from 15 and 20 receiver coils, respectively. Both datasets were retrospectively undersampled using a uniform/equidistant pattern at $R=4$. 24 lines of auto-calibration signal (ACS) from center of the raw k-space data were kept. DICOM images to train our proposed model were reconstructed using parallel imaging (CG-SENSE), solving $\mathbf{x}_{\text{PI}}=(\mathbf{E}_{\Omega}^{H}\mathbf{E}_{\Omega})^{-1}\mathbf{E}_{\Omega}^{H}\mathbf{y}_{\Omega}$. For each dataset, models were trained using $300$ slices, and testing was performed using $380$ slices for knee MRI and $100$ slices for brain MRI, from distinct subjects.

Prospective Undersampling Setup. A multi-echo 3D GRE sequence on a 7T Siemens Magnetom MRI scanner was acquired. In this experiment, we replicate the practical pipeline for CUPID, where data is acquired at the desired high acceleration rate, and reconstructed to ${\bf x}_{\textrm{PI}}$ with noise and aliasing artifacts, using parallel imaging. The corresponding DICOM images are exported and used for fine-tuning the PD-DL model with CUPID. To this end, the brain dataset, with matrix size $=288\times 288$ and in-plane resolution $0.7\times 0.7\textrm{mm}^{2}$, was acquired with prospective undersampling $R=9$ (in $k_{y}$ only), which is the desired target acceleration, much higher than the clinical protocol at $R=3$. Low-resolution images were acquired in the same orientation for sensitivity estimation (Krueger et al., 2023). Training and reconstruction with CUPID was done in a zero-shot subject-specific manner.

More details about the implementation of the PD-DL models are provided in Appendix A.

We compared our method against several database training methods that have access to raw k-space data, including supervised PD-DL (Hammernik et al., 2018; Aggarwal et al., 2019; Knoll et al., 2020a), self-supervision via data undersampling (SSDU) (Yaman et al., 2020), and equivariant imaging (EI) (Chen et al., 2021). All PD-DL methods utilized the same unrolled network and components (Appendix A) to ensure that only the training process differed for fair comparisons.

In addition, we compared our approach with methods that can operate without raw data access as long as ${\bf E}_{\Omega}$ is known at test time. These include compressed sensing (CS) (Lustig et al., 2007), and ScoreMRI (Chung & Ye, 2022). The latter trains a time-dependent score function using denoising score matching on a large dataset of reference fully-sampled images, and uses this score function during inference to sample from the conditional distribution given the measurements. Note both CS and ScoreMRI use $\mathbf{E}_{\Omega}^{H}\mathbf{y}_{\Omega}$ for data fidelity during inference. This can be accessed by multiplying $\mathbf{x}_{\text{PI}}$ with $\mathbf{E}_{\Omega}^{H}\mathbf{E}_{\Omega}$. Note that ${\bf E}_{\Omega}$ includes information about the undersampling pattern $\Omega$, which is completely known from the acquisition parameters, and coil sensitivities, which can be estimated from separate calibration scans in DICOM format (Krueger et al., 2023). A similar observation applies to the data fidelity in (5) for unrolled networks, thus they can be used for inference using only $\mathbf{x}_{\text{PI}}$ and ${\bf E}_{\Omega}$. We emphasize that what sets CUPID apart from other PD-DL strategies is that it is the only one that can train the unrolled network without using ${\bf y}_{\Omega}$. Thus, without loss of generality, ${\bf E}_{\Omega}$ is known both at training and testing for all methods. Finally, we also used CG-SENSE, which was used to generate the original ${\bf x}_{\textrm{PI}}$ as the clinical baseline comparison.

In the zero-shot setup, we compared our zero-shot results with zero-shot SSDU (ZS-SSDU) (Yaman et al., 2022b) as well as ScoreMRI, compressed sensing (CS), and our baseline method, CG-SENSE - all of which are compatible with zero-shot inference. All quantitative evaluations used structural similarity index (SSIM) and peak signal-to-noise ratio (PSNR).

Database Results. Representative results in Fig. 3 show that baseline CG-SENSE, CS, EI and ScoreMRI reconstructions exhibit residual artifacts. In contrast, CUPID successfully eliminates these artifacts from the CG-SENSE image using a well-trained PD-DL network, achieving state-of-the-art reconstruction quality comparable to supervised PD-DL and SSDU, despite only having access to ${\bf x}_{\textrm{PI}}$ for training, and not to raw k-space data unlike these other methods. We observe that parallel imaging reconstruction is not clinically usable at higher acceleration rates, but it is improved using a PD-DL reconstruction trained with CUPID. We further note that mild blurring was observed in some slices for database-training only. This is expected since we are no longer benefiting from redundancies from across multiple coils due to having no access to multi-coil raw k-space data, unlike the comparison methods. Quantitative results presented in Tab. 1 validate the visual observations, demonstrating that CUPID consistently outperforms CG-SENSE, CS, EI, and ScoreMRI across multiple datasets. Moreover, CUPID maintains performance comparable to that of supervised PD-DL and SSDU. Additional qualitative results on ax-FLAIR database and quantitative results including standard error of the mean are given in Appendix C and Appendix D, respectively.

Zero-Shot Learning Results. Fig. 4 shows results from zero-shot reconstructions. Both CG-SENSE and CS suffer from noise amplification and persistent residual artifacts, with CG-SENSE displaying a more pronounced degradation in quality. On the other hand, ScoreMRI exhibits blurring while still displaying residual artifacts. We note that uniform undersampling is used in these datasets, which is consistent with clinical parallel imaging acquisitions, but which have not been previously reported with diffusion models in existing works. Once again, CUPID demonstrates superior artifact and noise reduction over these methods and closely matches the quality of ZS-SSDU, despite not having access to raw data and an explicit self-validation mechanism to prevent overfitting as in the latter.

As discussed in Sec. 4.1, brain data is acquired at the target acceleration rate, reconstructed via parallel imaging and exported in DICOM format to perform zero-shot fine-tuning. Fig. 5 shows reconstruction results for the vendor parallel imaging reconstruction, as well as CS, ScoreMRI and CUPID. CS reduces the noise in the parallel imaging reconstruction, but leads to blurring due to over-regularization. In contrast, ScoreMRI struggles to reconstruct accurately at this high acceleration rate, suggesting generalizability issues for the pre-trained score function to high-resolution imaging at a different field strength, not represented in the training database, and potential difficulties with uniform undersampling. Furthermore, the public implementation of ScoreMRI uses hard constraints for data fidelity, which leads to more pronounced artifacts due to potential phase mismatch with the coils generated from separate calibration data. Our proposed CUPID method effectively mitigates both artifacts and noise in the DICOM image (shown in zoomed insets) without requiring any raw k-space data, attesting to the effectiveness of CUPID in real-world scenarios. Note minor residual artifacts remain since the target acceleration $R=9$ in 1-dimension is very high. We note that ZS-SSDU cannot be applied here due to the unavailability of raw data. We further note that the vendor-provided DICOM was generated using k-space interpolation (Griswold et al., 2002) instead of the image domain formulation in (2). Due to their equivalence, this did not cause any issues for CUPID, as expected.

We carried out two ablation studies to explore key factors influencing the performance of our algorithm. The first study explored the effect of $\lambda$ parameter to the final reconstruction, by training 5 distinct PD-DL networks using $\lambda\in\{0,50,100,200,\infty\}$. We note that using $\lambda=0$ corresponds to using only the compressibility term ($\mathcal{L}_{\text{comp}}$ in (8)), whereas using $\lambda\to\infty$ translates to using solely the parallel imaging fidelity term ($\mathcal{L}_{\text{pif}}$ in (9)). Fig. 6 shows the corresponding reconstruction results for each case. As outlined in Sec. 3, only using $\mathcal{L}_{\text{comp}}$ leads to overly-smooth reconstructions due to network forcing the wavelet coefficients towards zero without maintaining consistency with the data. On the other hand, solely using $\mathcal{L}_{\text{pif}}$ results in DIP-like reconstructions (Ulyanov et al., 2018), where the network overfits the data without any regularization, resulting in noise amplification. CUPID with $\lambda\in\{50,100,200\}$ integrates both loss terms to attain high-fidelity reconstructions. Thus, we conclude that CUPID demonstrates robust performance across a wide range of $\lambda$ values. Our second ablation study focuses on the effect of the number of perturbation patterns used in the training, and is provided in Appendix B.1.

Filtering on Routine Clinical Images. MR scans may include filtering operations applied by some vendors that affect the assumption $\mathbf{\hat{x}}_{\text{PI}}=(\mathbf{E}_{\Omega}^{H}\mathbf{E}_{\Omega})^{-1}\mathbf{E}_{\Omega}^{H}\mathbf{y}_{\Omega}$. This was discussed extensively in Shimron et al. (2022), in the context of using retrospective undersampling of DICOM images to train DL reconstruction, especially highlighting the use of zero-padding, which improves the display resolution compared to the acquisition resolution. It was shown that training of models from retrospective undersampling of databases of DICOM images for PD-DL training using zero-padding may lead to biases and inaccuracies. Conversely, our approach is physics-driven in nature, and the sampling pattern $\Omega$ naturally accounts for the zero-padding operation. However, our method is not immune to other types of filtering/processing, such as implicit intensity correction (Han et al., 2001) or deidentification methods (Van Essen et al., 2013), in which case the filtered ${\bf x}_{\textrm{PI}}$ would need to be treated as the parallel imaging solution corresponding to a filtered version of ${\bf y}$.

Resources for Fine-Tuning of PD-DL Reconstruction. While our method is aimed to improve equitable access to fast MRI in low-resource settings, we do acknowledge that such low-resource MRI centers may lack the necessary hardware to fine-tune PD-DL models. We note that, similar to what has been shown in Yaman et al. (2022b), transfer learning along with zero-shot fine-tuning may be beneficial in decreasing training time and resources. Furthermore, since only DICOM images are needed for CUPID, standard anonymization techniques can be used (Van Essen et al., 2013), and data can be transferred to potential offsite locations with more computational resources. Since no raw data storage or transfer is necessary, the burden for transmission of DICOMs is low, and privacy can be preserved thorough anonymization, though specifics would need to be implemented with proper protocols.

Magnitude-only DICOMs. In some cases, MR vendors only give access to the magnitude of $\mathbf{x}_{\text{PI}}$ due to limitations in the reconstruction pipelines, for instance when using partial Fourier imaging in some vendors. Additionally, if a reference calibration scan is not prospectively acquired at exam time, coil sensitivities, and thus ${\bf E}_{\Omega}$ would be unknown as well, which is the case for most existing DICOM databases. In such cases, the phase of ${\bf x}_{\text{PI}}$ and coil sensitivities need to be estimated jointly. There has been work on the latter in the standard PD-DL setup with access to raw k-space data (Hu et al., 2024; Arvinte et al., 2021). However, this is a modification on the PD-DL objective function, and not the learning procedure considered here. Thus, it is not the focus of our study and it will be investigated in future studies.

Reinforcement Learning for Optimal Perturbations. Finally, a promising direction for future work may involve using reinforcement learning to simultaneously learn the optimal perturbations. This approach may enable the model to adapt and optimize perturbation strategies so that the sample mean estimate in (9) can be better approximated with a few perturbation examples.