1 Introduction↩︎

Accelerated MRI reconstruction has substantially improved in recent years through the use of machine learning methods. Currently, diffusion models (DMs) [1], [2] achieve state-of-the-art MRI reconstruction scores [3]–[5]. An often neglected area in the literature is the inherent uncertainty of the reconstruction process. Accelerated MRI reconstruction is an ill-posed problem and infinite plausible reconstructions exist for an undersampled MRI image. However, in the majority of recent work only a single point estimate (a single reconstruction) of the solution space is provided. This common practice is problematic in safety-critical applications like medical imaging, as it does not accurately reflect the inherent uncertainty and no estimation of the confidence of the solution can be provided.

A common approach to quantify the reconstruction uncertainty is to reconstruct a large number of samples and analyze their pixel-wise variance. This kind of uncertainty prediction has been investigated with probabilistic networks [6], [7], Monte Carlo dropout [8] ensembling [9] or DMs [4], [5], [10].

Recently however, [11] has shown that Repeated Reconstruction (RR) in natural image reconstruction with DMs leads to similar images which display a lack of meaningful diversity. Even though the diffusion model is trained to model the posterior distribution accurately, the heavy tail of the true posterior distribution is often largely ignored in practice and plausible yet unlikely solutions with different semantic meaning are not found. This is a problem for medical imaging, where certain pathologies might be rare and unlikely. However, it is important to recover these cases precisely and include them as possible solutions when reconstructing multiple samples for uncertainty prediction.

The definition of meaningful diversity can depend on the downstream task and might be different for e.g. pathology detection (a pathology exists or not) or semantic segmentation (the segmented tissue might be smaller or larger). In this work we focus on diversity of semantic segmentation, as it is an important task in clinical practice and the results are more easily interpretable for non-experts.

In this paper, instead of random sampling, we propose to reconstruct two MRI-images per segmentation class corresponding to the upper and lower bound segmentation volume. To this end, we introduce a novel method called Segmentation Guided Reconstruction (SGR). SGR builds on DM-based MRI reconstruction. However, rather than standard sampling, we guide the reconstruction process with the gradient of segmentation losses that either penalize small segmentation volumes or large segmentation volumes. The resulting lower and upper bound reconstructions and segmentations are an intuitive way of understanding the segmentation uncertainty which arises due to the ill-posed nature of the reconstruction problem. See Fig. 1 for a visual explanation of the motivation. In summary, our contributions are:

• We introduce a novel method Segmentation-Guided Reconstruction (SGR) to guide the reconstruction process of undersampled MRI data, leading to meaningfully diverse reconstructions.

• We show that the diverse reconstructions lead to diverse segmentations that correspond to an upper and lower bound on the segmentation. We introduce the concept of ‘uncertainty boundary’ for the volume between both bounds.

• We show that our method captures the inherent uncertainty more faithfully than the commonly used Repeated Reconstruction (RR), potentially leading to safer MRI reconstruction methods

2 Related Work↩︎

Diffusion based inverse problem solvers have seen large success in MRI reconstruction [4], [5], [10], [12], achieving state-of-the-art reconstruction scores [3]–[5]. In contrast to supervised learning-based methods [13], [14], trained DMs can be used for several acceleration factors. Recent reconstruction methods [3] allow for more than \(80 \times\) faster inference time, reducing the computational burden for RR. With RR it is possible to reconstruct several distinct images, allowing to estimate the reconstruction uncertainty by analyzing the variance of the voxel values. However, the distinct images do not necessarily display meaningful diversity [11].

While the reconstruction and segmentation problem are mostly tackled separately, some work exists aiming to improve segmentation results by combining reconstruction and segmentation in an end-to-end fashion [15]–[19]. Wu et al. [18] show that fine-tuning a pre-trained reconstruction network with a task-specific network head leads to improved segmentation results. Similarly, the winning method of the K2S challenge [19] also first trains a reconstruction network before fine-tuning the network for segmentation. These approaches, however, lead to degraded reconstructions. Moreover, all of these methods lack an analysis of the effects of the inherent reconstruction uncertainty on segmentation.

The literature on reconstruction uncertainty propagation to downstream tasks is relatively sparse and most methods utilize repeated sampling of probabilistic networks that might not offer meaningful diversity [11]. Feiner et al. [20] propose to utilize Monte Carlo sampling to estimate the reconstruction uncertainty that can then be propagated to classification tasks. Cheung et al. [7] use repeated sampling of CT reconstructions to obtain a distribution of several downstream metrics, which they use to calculate statistically valid prediction sets using conformal prediction. Fischer et al. [6] are the first to consider the propagation of the aleatoric uncertainty due to the reconstruction process in the context of segmentation. By repeated sampling of a probabilistic reconstruction network and subsequent segmentation, they were able to quantify the variance of the segmentation at pixel-level. However, data-consistency was not enforced and no upper or lower bound of the segmentation has been estimated.

3 Method↩︎

3.1 Background diffusion models and MRI-reconstruction↩︎

In the inverse problem of MRI reconstruction, we aim to reconstruct a plausible image \(\boldsymbol{x}_0\) that is coherent with the measurement data \(\boldsymbol{y}\) when passing through a forward operator (i.e. Fourier Transform + Masking) \(\boldsymbol{A}\), such that the following equation holds: \[\boldsymbol{y} = \boldsymbol{A} \boldsymbol{x}_0^{(i)} \label{eq:inverse95problem} \;\text{.}\tag{1}\]

The problem is typically ill-posed and infinite \(\boldsymbol{x}_0^{(i)}\) exist that are possible solutions for Eq. 1 . In recent years DMs have shown to perform exceptionally well on the MRI-reconstruction tasks, being able to sample a set of distinct solutions \(\boldsymbol{x}_0^{(i)}\) from the posterior.

DMs generate samples from the data distribution \(p_{data}(\boldsymbol{x}_0)\) by reversing a diffusion process from timestep \(t=0\) to \(t=T\) described by a conditional density

\[q(\boldsymbol{x}_t|\boldsymbol{x}_0) = \mathcal{N}(\boldsymbol{x}_t|\sqrt{\bar{\alpha_t}} \boldsymbol{x}_0, (1-\bar{\alpha_t})\boldsymbol{I})\] with a pre-defined defined \(\bar{\alpha_t}\) schedule.

A DM can then be trained by predicting the noise through epsilon matching: \[\min\limits_{\theta}\mathbb{E}_{\boldsymbol{x}_t \sim q(\boldsymbol{x}_t|\boldsymbol{x}_0), \boldsymbol{x}_0 \sim p_{data}(\boldsymbol{x}_0), \boldsymbol{\epsilon} \sim \mathcal{N}(0,\boldsymbol{I})} (||\boldsymbol{\epsilon}_{\theta}^{(t)}(\boldsymbol{x}_t) - \boldsymbol{\epsilon}||_2^2) \;\text{.}\] The learned \(\boldsymbol{\epsilon}\)-matching function can then be used to predict the clean image \(\hat{\boldsymbol{x}}_{0|t}\) at every step \(t\) of the diffusion process using Tweedie denoising [21]:

\[\hat{\boldsymbol{x}}_{0|t} = (\boldsymbol{x}_t - \sqrt{1-\bar{\alpha_t}} \boldsymbol{\epsilon}_{\theta}^{(t)}(\boldsymbol{x}_t)) / \sqrt{\bar{\alpha_t}} \text{.}\]

Recently, Chung et al. [3] introduced a DM based reconstruction method that utilizes conjugate gradient (CG) optimization: During the CG-optimization process, starting from \(\hat{\boldsymbol{x}}_{0|t}\) we search a \(\hat{\boldsymbol{x}}_{0|t}'\) that minimizes the distance \(\boldsymbol{A}^*\boldsymbol{A} \hat{\boldsymbol{x}}_{0|t}' - \boldsymbol{A}^*\boldsymbol{y}\) with \(\boldsymbol{A}^*\) being the backward operator mapping from the measurement space back to the image space:

\[\label{eq:conjugate95gradient} \hat{\boldsymbol{x}}_{0|t}' = CG(\boldsymbol{A}^{*}\boldsymbol{A}, \boldsymbol{A}^*\boldsymbol{y}, \hat{\boldsymbol{x}}_{0|t}) \text{.}\tag{2}\] The update rule using the DDIM method [22] is then given as follows: \[\boldsymbol{x}_{t-1} = \sqrt{\bar{\alpha}_{t-1}} \hat{\boldsymbol{x}}_{0|t}' + \sqrt{1-\bar{\alpha}_{t-1} - \eta^2 \tilde{\beta}_t^2} \boldsymbol{\epsilon}_{\theta}^{(t)}(x_t) + \eta \tilde{\beta}_t \boldsymbol{\epsilon} \text{.} \label{eq:ddim95update}\tag{3}\] We modify the algorithm from [3] slightly to achieve full data-consistency. For the last diffusion step \(t=1 \rightarrow t=0\), we replace Eq. 2 and make sure that the reconstructed image \(\boldsymbol{x}_0\) is fully data-consistent using \(\boldsymbol{y}\) and Eq. 1 .

3.2 Segmentation Guidance for diverse sampling↩︎

In order to generate a set of meaningfully diverse solutions \(\boldsymbol{x}_0^{(i)}\), we reconstruct two images per segmentation class (\(\boldsymbol{x}_0^{\uparrow}\) and \(\boldsymbol{x}_0^{\downarrow}\)) with corresponding upper and a lower bound segmentations (\(S^{\uparrow}\) and \(S^{\downarrow}\)) as well as segmentation volumes (\(V_{seg}^{\uparrow}\) and \(V_{seg}^{\downarrow}\)). The volume between \(S^{\downarrow}\) and \(S^{\uparrow}\) is the uncertainty boundary \(S_{unc}\). We guide the reconstruction process of the diffusion model using segmentation guidance. To achieve this, we propose two loss-functions \(\mathcal{L}^{\uparrow}\) and \(\mathcal{L}^{\downarrow}\) that are suitable for segmentation-based model guidance.

Let \(p_{\phi}^c(\boldsymbol{x}_{0,jk})\) be the predicted probability of pixel \((j,k)\) in the image \(\boldsymbol{x}_0\) belonging to the segmentation class \(c \in \mathcal{C}\), and let \(p_{\phi}\) be a segmentation network. We abbreviate the term as \(\boldsymbol{p}_{jk}^{(c)}\). The volume of \(c\) can be calculated by summing the number of pixels that belong to \(c\) multiplied by voxel-volume \(V_{voxel}\):

\[V_c(\boldsymbol{x}_0) = \sum_{j,k} (\mathbb{I}_c(\boldsymbol{p}_{jk}^{(c)})) \cdot V_{voxel} \text{,}\] where \(\mathbb{I}_c\) is an indicator function

\[\mathbb{I}_c = \begin{cases} 1 & \text{if } \boldsymbol{p}_{jk}^{(c)} > \boldsymbol{p}_{jk}^{(c')} \text{ } \forall c' \in \mathcal{C}, c' \neq c \\ 0 & \text{otherwise} \end{cases} \text{.}\]

During optimization, we consider only one class at a time, such that we can simplify the objective to a binary segmentation class. We build our guidance-loss based on the widely-used Binary-Cross-Entropy loss:

\[\mathcal{L}_{BCE}(\boldsymbol{y}, \boldsymbol{p}, c) = - \sum_{j,k} \Bigl( \boldsymbol{y}_{jk} \log(\boldsymbol{p}_{jk}^{(c)}) + (1-\boldsymbol{y}_{jk}) \log(1-\boldsymbol{p}_{jk}^{(c)}) \Bigr) \text{.} \label{eq:binary95ce}\tag{4}\]

We can maximize the optimization objective \(\mathcal{L}^{\uparrow}\) to find the upper bound segmentation, as well as \(\mathcal{L}^{\downarrow}\) to find the lower bound segmentation:

\[\mathcal{L}^{\uparrow} = \mathcal{L}_{BCE}(\boldsymbol{0}, \boldsymbol{p}, c) = - \sum_{j,k} \log(1-\boldsymbol{p}_{jk}^{(c)}) \text{,}\]

\[\mathcal{L}^{\downarrow} = \mathcal{L}_{BCE}(\boldsymbol{1}, \boldsymbol{p}, c) = - \sum_{j,k} \log(\boldsymbol{p}_{jk}^{(c)}) \text{.}\]

The binary cross-entropy loss is unbounded, meaning that pixel probabilities are further optimized, even if the target class has already been reached. Following Croce et al. [23] we can avoid this behavior by using a masked cross-entropy loss, which excludes pixels already segmented as the target class from the optimization process. The optimization objectives can therefore be written as follows:

\[\mathcal{L}_{m}^{\downarrow} = \mathbb{I}_c \mathcal{L}^{\downarrow}; \quad \mathcal{L}_{m}^{\uparrow} = (1-\mathbb{I}_c) \mathcal{L}^{\uparrow} \text{.}\]

The derived optimization objective can be used during the reverse diffusion process. We therefore calculate the gradient of \(\mathcal{L}_m^{\uparrow }\) (or \(\mathcal{L}_m^{\downarrow}\), respectively) with respect to \(x_t\) at every step \(t\), to find the perturbation necessary to increase the optimization loss. We add the gradient \(\nabla_{x_t} \mathcal{L}_{m}^{\uparrow}\) (or \(\nabla_{x_t} \mathcal{L}_{m}^{\downarrow}\), respectively) during the DDIM update in Eq. 3 :

\[\boldsymbol{x}_{t-1} = \sqrt{\bar{\alpha}_{t-1}} \hat{\boldsymbol{x}}_{0|t}' + \sqrt{1-\bar{\alpha}_{t-1} - \eta^2 \tilde{\beta}_t^2} (\boldsymbol{\epsilon}_{\theta}^{(t)}(x_t)+ \gamma \nabla_{\boldsymbol{x}_t}\mathcal{L}_{m}^{\uparrow / \downarrow}) + \eta \tilde{\beta}_t \boldsymbol{\epsilon} \text{,}\] where the parameter \(\gamma\) makes sure that the \(l_2\) norm of \(\nabla_{x_t} \mathcal{L}_{m}^{\uparrow / \downarrow}\) stays small with respect to \(\boldsymbol{\epsilon}_{\theta}\). We set \(\gamma\) to the following value:

\[\gamma = \begin{cases} b \cdot ||\boldsymbol{\epsilon}_{\theta}^{(t)}||_2 / ||\nabla_{\boldsymbol{x}_t} \mathcal{L}_{m}^{\uparrow / \downarrow}||_2 & \text{if ||\nabla_{\boldsymbol{x}_t} \mathcal{L}_{m}^{\uparrow / \downarrow}||_2 > b ||\boldsymbol{\epsilon}_{\theta}^{(t)}||_2} \\ 1 & \text{otherwise} \text{,} \end{cases}\] where \(b\) is a small constant that we set to \(b=0.005\) in our experiments. An ablation study on the effects of the parameter value can be found in the Appendix. A summary of our method is shown in Fig. 2.

For the baseline RR, we sample 16 images per slice and select \(x_0^{\uparrow}\) and \(x_0^{\downarrow}\) for every segmentation class as the reconstruction leading to the highest segmentation volume (\(V_{seg}^{\uparrow}\)) and lowest segmentation volume (\(V_{seg}^{\downarrow}\)), respectively.

3.3 Data and models↩︎

We train and test our method on the publicly available SKM-TEA [24] dataset and focus only on the \(E1\)-echoes. As the spacing of the SKM-TEA data is very fine-grained (0.3125 mm in the axial (\(=\)sampling) direction), neighboring images tend to be similar. We therefore save computational resources during testing by only considering every 8th slice in the relevant regions, corresponding to a spacing of 2.5mm between the slices. In total we analyze \(869\) slices for each acceleration factor and each method.

All experiments are based on the same diffusion model. We train the diffusion model on the fully-sampled images of the SKM-TEA data. We normalize the data by dividing by the 0.99 quantile of every slice. The diffusion model is based on the one introduced in [25]. We train the model for 1 million steps using a batch-size of 2 on a 3080TI GPU. For segmentation, we train a U-Net on the fully-sampled MRI-data in axial direction. Due to the strong class-imbalance on the axial imaging plane, where only a minority of slices show elements of tibial cartilage or meniscus as these parts tend to be relatively flat, we perform a foreground oversampling of 40%. To reduce overfitting, we apply mirroring and Gaussian Noise augmentations during training.

4 Experiments and results↩︎

We perform experiments on 4 acceleration factors: \(4\times\), \(8\times\), \(12\times\), and \(16\times\), analyzing the behavior of the uncertainty boundaries of SGR and RR. We assess the quality of the generated segmentations \(S^{\downarrow}\) and \(S^{\uparrow}\) by calculating the precision for \(S^{\downarrow}\), recall for \(S^{\uparrow}\) and the uncertainty volume \(V_{unc}\). We focus on precision on \(S^{\downarrow}\), as we want to avoid overestimating the segmentation and creating false-positives, and measure recall for \(S^{\uparrow}\), because we want to avoid segmentations containing false-negatives. An example image where RR overestimates \(S^{\downarrow}\) while our method predicts accurate boundaries is shown in Fig. 3. A similar observation can be found in Fig. 5, where our method sustains high precision and recall values, while these scores degrade with higher acceleration for the baseline method RR. Reconstruction quality as measured by SSIM and PSNR remains very high for SGR compared to RR, as shown in Tab. 1.

In Fig. 4, we see that the uncertainty boundary increases with higher acceleration for our method, while it stays approximately constant for the baseline method RR. This result is also shown in Fig. 5, where the ratio \(V_{unc}/V_{seg}^{\uparrow}\) increases faster for SGR compared to RR. Given the relatively low precision and recall scores for RR and the higher scores for our method, one can conclude that RR underestimates the uncertainty boundaries, while this problem is largely reduced for the uncertainty boundaries estimated by SGR.

5 Conclusion and Discussion↩︎

We have introduced a novel method for segmentation-guided MRI reconstruction and have shown that it can generate meaningfully diverse reconstructions that lead to upper and lower bound segmentations. These segmentation can offer an intuitive way of understanding the segmentation uncertainty that is caused by the ill-posed nature of the reconstruction problem. As our method also offers more reliable uncertainty boundaries compared to the standard method of repeated reconstruction, we believe that it can be helpful in clinical practice, whenever a trustworthy uncertainty estimation is required.

In future work, we want to analyze if robust segmentation networks introduced in [23] can help to reduce \(V_{unc}\) while maintaining high precision and recall scores. In this work, however, we have focused on more commonly used neural networks, to show the easy applicability of our method for existing pipelines.

6 Supplementary Material↩︎

6.1 Ablation study for the weighting hyperparameter \(b\)↩︎

References↩︎