Deep evidential fusion with uncertainty quantification and reliability learning for multimodal medical image segmentation

The information provided by a mass function $m$ can also be represented by a belief function $Bel$ or a plausibility function $Pl$ from $2^{\Theta}$ to $[0,1]$ defined, respectively, as:

and

for all $A\subseteq\Theta$, where $\bar{A}$ denotes the complement of $A$. The quantity $Bel(A)$ can be interpreted as a degree of support for $A$, while $Pl(A)$ is a measure of lack of support against $A$. The contour function $pl$ associated to $m$ is the function that maps each element $\theta$ of $\Theta$ to its plausibility, i.e.,

As shown below, this function can be easily computed when combining several pieces of evidence; it plays an important role in decision-making.

In DST, the beliefs about a certain question are established by aggregating independent pieces of evidence represented by belief functions over the same frame of discernment [69]. Given two mass functions $m_{1}$ and $m_{2}$ derived from two independent items of evidence, the mass function $m_{1}\oplus m_{2}$ representing the pooled evidence is defined as

This operation is called Dempster’s rule of combination. It is commutative and associative. The combined mass function $m_{1}\oplus m_{2}$ is called the orthogonal sum of $m_{1}$ and $m_{2}$. Mass functions $m_{1}$ and $m_{2}$ can be combined if and only if $\kappa<1$. Let $pl_{1}$, $pl_{2}$ and $pl_{1}\oplus pl_{2}$ denote the contour functions associated with, respectively, $m_{1}$, $m_{2}$ and $m_{1}\oplus m_{2}$. The following equation holds:

The complexity of calculating the combined contour function using (2) is linear in the cardinality of $\Theta$, whereas computing the combined mass function using (1) has, in the worst-case, exponential complexity.

Given a mass function $m$ and a nonempty subset $A$ of $\Theta$ such that $Pl(A)>0$, the conditional mass function $m(\cdot|A)$ is defined as the orthogonal sum of $m$ and the mass function $m_{A}$ such that $m(A)=1$. Conversely, given a conditional mass function $m_{0}$ given $A$ (expressing one’s beliefs in a context where it is only known that the truth lies in $A$), its conditional embedding [73] is the least precise mass function $m$ on $\Theta$ such that $m(\cdot|A)=m_{0}$; it is obtained by transferring each mass $m_{0}(C)$ to $C\cup\overline{A}$, for all $C\subseteq A$. Conditional embedding is a form of “deconditioning”, i.e., it performs the inverse of conditioning.

Once a mass function representing the combined evidence has been computed, it is often used to make a decision. Decision-making methods in DST are reviewed in [18]. Here, we will use the simplest method [11], which consists in computing a probability distribution on $\Theta$ by normalizing the plausibilities of the singletons,

Once probabilities have been computed, a decision can be made by maximizing the expected utility. We note that this method fits well with Dempster’s rule, as the plausibility of the singletons can be easily computed from (2) without computing the whole combined mass function.

In [58], the authors generalize the discounting operation using the notion of contextual discounting, which makes it possible to account for richer metaknowledge about the reliability of a source in different contexts, i.e., conditionally on different hypotheses regarding the variable of interest. In the corresponding refined model, $m(\cdot|R)$ and $m(\cdot|\neg R)$ are defined as before, but our beliefs about the reliability of the source are now defined by $K$ coefficients $\beta_{1},\ldots,\beta_{K}$, one for each state in $\Theta$. More specifically, we have $K$ conditional mass functions defined by $m_{\mathcal{R}}(\{R\}|\theta_{k})=\beta_{k}$ and $m_{\mathcal{R}}(\mathcal{R}|\theta_{k})=1-\beta_{k}$, for $k=1,\ldots,K$. In this model, $\beta_{k}$ is, thus, the degree of belief that the source of information is reliable, given that the true state is $\theta_{k}$. As shown in [58], combining the conditional embeddings of $m(\cdot|R)$ and $m_{\mathcal{R}}(\cdot|\theta_{k})$ for $k=1,\ldots,K$ by Dempster’s rule yields the following discounted mass function,

for all $A\subseteq\Theta$, where ${\boldsymbol{\beta}}=(\beta_{1},\ldots,\beta_{K})$ is the vector of all reliability coefficients, and a product of terms is equal to 1 if the index set is empty. In many applications, we actually do not need to compute the whole mass function (5): we can compute only the associated contour function ${}^{{\boldsymbol{\beta}}}pl$, which is all we need for decision-making. As shown in [58], this contour function is equal to

It can be computed in linear time with respect to the size of $\Theta$, instead of exponential time for ${}^{{\boldsymbol{\beta}}}m$. An evidential $k$ nearest neighbor rule based on the contextual discounting operation was introduced in [22].

The proposed framework was tested on two multimodal medical image datasets.

The PET-CT lymphoma dataset contains 3D images from 173 patients who were diagnosed with large B-cell lymphomas and underwent PET-CT examination³³3The study was approved as a retrospective study by the Henri Becquerel Center Institutional Review Board.. For lymphoma segmentation, PET imaging helps identify active tumor sites by highlighting areas of increased metabolic activity. In contrast, CT imaging provides anatomical information about the size, shape, location, and surrounding structures of lymphoma tumors. While PET image makes it possible to obtain functional information about the tumor and surrounding tissues, CT images provide complementary anatomical details allowing for more accurate segmentation. The lymphomas in mask images were delineated manually by experts and considered as ground truth. Figure 4 shows an example of PET and CT images of a patient with lymphomas. The PET and CT images and the corresponding mask images have different sizes and spatial resolutions due to the use of different imaging machines and operations. For CT images, the size varies from $267\times 512\times 512$ to $478\times 512\times 512$. For PET images, the size varies from $276\times 144\times 144$ to $407\times 256\times 256$.

The multi-MRI brain tumor dataset was made available for the BraTS2021 challenge [5]. The original BraTS2021 dataset comprises training, validation, and test sets with, respectively, 1251, 219, and 570 cases. There are four modalities: FLAIR, T1Gd, T1, and T2 with $240\times 240\times 155$ voxels. Figure 5 shows examples of four-modality MRI slices for one patient. The appearance of brain tumors varies in different modalities [5]. T1Gd MRI images are obtained following the administration of a gadolinium-based contrast agent that enhances areas with disrupted blood-brain barrier such as tumor regions, making tumors appear hyperintense (bright) and improving the visibility of tumor margins. FLAIR MRI images suppress the signal from cerebrospinal fluid (CSF), highlighting pathological changes while suppressing the CSF signal. T2 MRI images are sensitive to tissue water content and provide good contrast between soft tissues. Tumors with increased water content often appear hyperintense (bright) on T2 images. T1 MRI images are crucial for identifying tumor location and structural details by their excellent anatomical detail. Annotations of scans comprise gadolinium (GD)-enhancing tumor (ET), necrotic and non-enhancing tumor core (NRC/NET), and peritumoral edema (ED). The task of the BraTS2021 challenge was to segment the images into three overlapping regions: ET, tumor core (TC, the union of ET and NRC/NET), and whole tumor (WT, the union of ET, NRC/NET, and ED). In this work, we evaluated the segmentation performances with respect to these three overlapping regions to allow a fair comparison with other state-of-the-art methods. Additionally, we also compared the results with respect to the three original non-overlapping tumor regions to highlight the impact of contextual discounting on subregion segmentation.

For the PET-CT dataset, we first normalized the PET, CT and mask images: (1) for PET images, we applied a random intensity shift and scale to each channel with a shift value of 0 and scale value of 0.1; (2) for CT images, the shift and scale values were set to 1000 and 1/2000; (3) for mask images, the intensity value was normalized into the $[0,1]$ interval by replacing the outside value by $1$. We then resized the PET and CT images to $256\times 256\times 128$ by linear interpolation and mask images to $256\times 256\times 128$ by nearest neighbor interpolation. Lastly, CT and PET images were registered using B-spline interpolation. Following [37], we randomly divided the 173 scans into subsets of size 138, 17, and 18 for, respectively, training, validation, and test. The training process was then repeated five times to test the stability of our framework, with different data used exactly once as the validation and test data.

For the BraTS2021 dataset, we used the same pre-processing operation as in [60]. We first performed a min-max scaling operation and clipped intensity values to standardize all volumes; we then cropped/padded the volumes to a fixed size of $128\times 128\times 128$ by removing the unnecessary background (the cropping/padding operation was only applied to training data). No data augmentation technique was applied, and no additional data was used in this study. Since the ground truth labels are unavailable for the validation and test sets, we trained and tested our framework with the training set. Following [60], we randomly divided the 1251 training scans into subsets of 834, 208, and 209 cases for training, validation, and testing, respectively. The process was repeated five times to test the stability of our framework. All the preprocessing methods mentioned in this paper can be found in the SimpleITK [53] toolkit.

All the compared methods used the same dataset composition and pre-processing operations. They were implemented in Python with the PyTorch-based medical image framework MONAI⁴⁴4More details about how to use those models can be found in MONAI core tutorials https://monai.io/started.html##monaicore..

At the FE stage, the number of filters in UNet was set to $(8,16,32,64,128)$ with kernel size equal to five and convolutional strides equal to $(2,2,2,2)$ for layers from left to right. For nnUNet used in Section 4.3, the kernel size was set to $(3,(1,1,3),3,3)$ and the upsample kernel size was set to $(2,2,1)$ with strides $((1,1,1),2,2,1)$. For nnFormer used in Section 4.3, the crop size was set to $(128,128,128)$ with embedding dimension set to $96$ and the number of heads was set to $(3,6,12,24)$. The number of extracted features was $H=2$ for the PET-CT lymphoma dataset and $H=4$ for the multi-MRI BraTS2021 dataset.

To train our fusion framework, we proceeded in three steps. First, FE modules (i.e., UNet, nnUNet, or nnFormer) were pre-trained independently for each modality during 50 epochs. Then, the weights of the FE modules were fixed, and the parameters of the EM and MMEF modules were optimized. Finally, the whole framework was fine-tuned for a few epochs. The initial values of parameters $\alpha_{i}$ and $\gamma_{i}$ in the EM modules were set to 0.5 and 0.01, and the membership degrees $u_{ik}$ were initialized randomly by drawing uniform random numbers, and normalizing. We used, $I=10$ prototypes for the PET-CT lymphoma dataset, and $I=20$ prototypes for the more complex multi-MRI BraTS2021 dataset. These prototypes were randomly initialized from a normal distribution with zero mean and an identity covariance matrix. Details about the initialization of the EM module can be found in [37]. The reliability coefficients $\beta^{t}_{k}$ in the MMEF module were initialized at 0.5.

For both datasets, we used the Adam optimization algorithm with an early stopping strategy: training was stopped when there was no improvement in performance on the validation set during ten epochs. The initial learning rate was set to $0.01$. The batch size was set to 4. For all the compared methods, the model with the best performance on the validation set was saved as the final model for testing⁵⁵5The code is available at https://github.com/iWeisskohl/Deep-evidential-fusion..

Although many authors have shown that segmentation performance can be improved by merging multimodal medical images into deep neural networks [63, 2], the reliability of information sources and the quality of uncertainty quantification have rarely been investigated. Here, the former issue will be addressed by analyzing the reliability coefficients $\beta_{k}^{t}$ defined in Section 3.1.3. To assess the quality of uncertainty quantification, we will use three metrics: the Brier score [9], the negative log-likelihood (NLL), and Expected Calibration Error (ECE) [31]. These metrics provide a robust evaluation framework for the uncertainty of the segmentation results, with smaller values indicating better performance. Their definitions are recalled below.

The Brier Score and NLL are defined, respectively, as

and

where $G_{n}$ is the ground truth of voxel $n$, $P_{n}$ is the predicted probability of voxel $n$, and $N$ is the number of voxels.

The ECE measures the correspondence between predicted probabilities and ground truth. The output normalized plausibilities of the model are first discretized into equally spaced bins $E_{b}$, $b\in[1,B]$ ($B=10$ in this paper). The accuracy of bin $E_{b}$ is defined as

where $S_{n}$ is the predicted class label for voxel $n$ and $\boldsymbol{1}(\cdot)$ is the indicator function. The average confidence of bin $E_{b}$ is defined as

The ECE is the weighted average of the difference in accuracy and confidence of the bins:

A model is perfectly calibrated when $\textsf{acc}(E_{b})=\textsf{conf}(E_{b})$ for all $b\in\{1,...,B\}$, in which case $\textsf{ECE}=0$.

Since our dataset has imbalanced foreground and background proportions, we only considered voxels belonging to the foreground or tumor region to calculate the above three indices. For the PET-CT lymphoma dataset, focusing only on the tumor region is not easy since the lymphomas are scattered throughout the whole body. Thus, we focused on the foreground region for this dataset. For the BraTS2021 dataset, we followed the suggestion from [66] to focus on the tumor region for the reliability evaluation. For each patient in the test set, we defined a bounding box covering the foreground or tumor region and calculated the corresponding values in this bounding box. For all segmentation performance criteria, the reported results were obtained by calculating the criteria for each test 3D scan and then averaging over the patients.

In addition to evaluating segmentation reliability, we also measured segmentation accuracy using the Dice score. In a segmentation task, the Dice score measures the volume of the overlapping region of the predicted object and the ground truth object as

where $TP$, $FP$, and $FN$ denote, respectively, the numbers of true positive, false positive, and false negative voxels.

The results concerning uncertainty estimation are reported in Table 1. Our model (MMEF-UNet) was compared to

1. 1.

   UNet with a softmax decision layer (the baseline);

2. 2.

   UNet with Monte-Carlo (MC) dropout [29] and deep ensemble [45], two popular techniques for improving the uncertainty quantification capabilities of probabilistic deep neural networks;

3. 3.

   ENN-UNet, composed of UNnet as the FE module and the EM module in place of the softmax layer; this is the architecture studied in [37];

4. 4.

   RBF-UNet, an alternative model composed of UNnet and a radial-basis function (RBF) module in place of the softmax layer; as shown in [37], this model makes it possible to compute output belief functions that are similar to those computed by ENN-UNet.

We can remark that approaches 1 to 4 above implement pixel-level fusion, whereas our approach is based on decision-level fusion. As for uncertainty quantification, UNet, UNet-MC and UNet-Ensemble are probabilistic methods. UNet only computes point estimates of class probabilities without taking into account second-order uncertainty. UNet-MC applies dropout during both training and inference, sampling multiple forward passes to estimate uncertainty by averaging the predictions. UNet-Ensemble quantifies uncertainty by averaging the predictions obtained from multiple independently-trained models. In contrast, ENN-UNet and RBF-UNet are evidential methods: they both calculate belief functions to represent segmentation evidence and uncertainty under the DST framework. For UNet-MC, the dropout rate was set to 0.2 and the number of samples was set to five; we averaged the five output probabilities at each voxel as the final output of the model. For UNet-ensembles, the number of samples was set to five; the five output probabilities were then averaged at each voxel as the final output of the model. The settings of ENN-UNet and RBF-ENN are the same as those reported in [37].

From Table 1, we can see that Monte-Carlo dropout and deep ensembles do not significantly improve the segmentation reliability as compared to the baseline UNet model, as shown, e.g., by the higher NLL values. In contrast, the addition of the EM module to the FE module, as implemented in ENN-UNet, brings a significant improvement, particularly according to NLL; the RBF-UNet model yields similar results. The decision-fusion framework MMEF-UNet brings an additional improvement according to all three criteria (ECE, Brier score, NLL) and outperforms the other models: specifically, we observe decreases of 1.1%, 0.9%, and 13% in ECE, Brier score, and NLL, respectively, as compared to UNet. We can conclude that, compared to the baseline model, both the EM and MMEF modules contribute to a higher segmentation reliability.

These findings are, to some extent, confirmed by Figure 6, which shows the calibration plots (also known as reliability diagrams) for the compared methods on the lymphoma dataset. Calibration plots are graphical representations showing how well the probabilistic predictions of a segmentation model are calibrated, i.e., how well confidence matches accuracy. In the left graph of Figure 6, we can see that the curve corresponding to UNet-Ensemble is closer to the diagonal than those of UNet and UNet-MC, which indicates better calibration. Looking at the right graph in Figure 6, we can see that the three DST-based models, ENN-UNet, RBF-UNet, and MMEF-UNet, have better calibration performance than the probabilistic ones, as shown by their calibration curves closer to the diagonal. Among them, MMEF-UNet shows the best calibration performance as ENN-UNet is slightly overconfident, while RBF-UNet is slightly underconfident.

The segmentation accuracy was measured by the Dice score, as shown in Table 1. Compared with the baseline model UNet, our proposal MMEF-UNet significantly increases segmentation performance, as shown by a 4.1% increase in the Dice score. Compared with the two DST-based deep evidential segmentation methods, MMEF-UNet has a higher Dice score (although the difference with ENN-UNet is not statistically significant). Figure 7 shows an example of visualized segmentation results obtained by UNet, ENN-UNet, RBF-UNet, and MMEF-UNet. We can see that UNet and RBF-UNet are more conservative (they correctly detect only a subset of the tumor voxels), while ENN-UNet is more radical (some of the voxels that do not belong to tumors are predicted as tumors). In contrast, the tumor regions predicted by MMEF-UNet better overlap the ground-truth tumor region, especially for the isolated lymphomas, which is also reflected by the promising Dice score value. These conclusions are consistent with the calibration trends displayed in Figure 6.

Table 2 reports the learned reliability coefficients. We can see that they are higher for the PET modality. This is consistent with domain knowledge, as mentioned in Section 4.1: PET images provide functional information about tumor activity and make it possible to identify active tumor sites, whereas CT images essentially provide detailed anatomical information (e.g., size, shape, and location) about lymph nodes and surrounding tissues and are used as a complement to PET images. This is also confirmed by the results presented in Table 3, showing that the performance of UNet and ENN-UNet with either CT alone or PET alone, the latter configuration yielding better results.

For the BraTS2021 dataset, we tested the segmentation performance of our fusion framework with UNet as well as two alternative FE modules: nnUNet and nnFormer. The nnUNet model was reported to have the best performance in the BraTS2021 challenge [54] and nnFormer is now one of the state-of-the-art brain tumor segmentation models. The complete frameworks with nnUNet and nnFormer as a feature extractor are referred to, respectively, as MMEF-nnUNet and MMEF-nnFormer. We compared our results with three baseline models: UNet, nnUNet, and nnFormer, and three Monte Carlo-based uncertainty segmentation models: UNet-MC, nnUNet-MC, and nnFormer-MC. Since the results obtained in Section 4.2, as well as those reported in [37] have shown that ENN-UNet and RBF-UNet yield similar results, here we only compared the performance of the ENN-based models, i.e., ENN-UNet, ENN-nnUNet and ENN-nnFormer. Moreover, we did not test the performance of deep ensemble models because applying them to larger-scale datasets exceeds our computation resources.

As with the lymphoma dataset, we used the ECE, Brier score, and NLL metrics to assess segmentation uncertainty. The results with UNet, nnUNet and nnFormer in the FE module are presented, respectively, in Tables 4, 5 and 6. We can see that our fusion model consistently outperforms the baseline models with all three FE models and across all uncertainty evaluation metrics, although the differences are more significant when UNet is used as a feature extractor. Indeed, the fusion mechanism can be expected to have a smaller impact when information sources are more informative. Overall, MMEF-nnUNet achieves the highest segmentation reliability with the lowest ECE, Brier score, and NLL values, and MMEF-nnFormer yields the second-best results.

Segmentation accuracy was evaluated by the Dice score for the three overlapping regions, ET, TC, and WT, as well as by the mean Dice score. The results with UNet, nnUNet and nnFormer as feature extractors are reported, respectively, in Tables 7, 8 and 9. Again, we can see that our fusion strategy improves segmentation accuracy for all three FE models. Overall, the highest segmentation accuracy was achieved by MMEF-nnFormer, with an increase of 1.5 % in the mean Dice score compared with the second-best method, ENN-nnFormer.

We also report the Dice score for the segmentation of the three original tumor regions: ED, ET, and NRC/NET in Table 10. As we can see, the baseline nnFormer shows good performance for segmenting ED and ET, while it does not perform as well for segmenting NRC/NET. Indeed, the lack of clear contrast, the similar signal intensities to normal brain tissue, the infiltrative growth patterns, and the need for multi-modal data make the segmentation of NRC/NET inherently more challenging compared to ED and ET. When the MMEF-nnFormer approach was applied, the Dice scores for the ED, ET, and NRC/NET improved by 0.6%, 1.6%, and 6.5%, respectively. The substantial improvement in NRC/NET segmentation is particularly encouraging, as it demonstrates the effectiveness of the proposed fusion method for delineating fuzzy tumor boundaries and solving challenging segmentation tasks.

Figures 8 and 9 show two segmentation cases when using nnFormer as the feature extractor. Figure 8 shows an easy segmentation case where only one tumor type is present. Both the Flair and T1Gd images exhibit good segmentation performance with only a few mislabeled voxels. It is surprising to see that concatenating multimodal medical images as the input for nnFormer resulted in worse outcomes, with the most mislabeled voxels. This might be due to the hard fusion strategy of nnFormer, i.e., image concatenation, which cannot mitigate the impact of noisy information. Consequently, the fused results are sometimes not as good as those from single-modality inputs. The proposed MMEF-nnFormer approach achieves the best performance, with fewer mislabeled voxels compared to other methods. Figure 9 illustrates a challenging segmentation scenario involving a tumor with ED, ET, and NRC/NET components. We can remark that the FLAIR image alone provides sufficient information to accurately segment ED, which is consistent with domain knowledge. Overall, the MMEF-nnFormer model yields the best results in this case. This example illustrates the ability of our method to improve segmentation accuracy by appropriately weighting and combining information from different modalities.

We first recall some clinical domain knowledge of MRI images in segmenting brain tumors:

1. 1.

   T1Gd images are particularly useful for delineating tumor boundaries by making tumor regions hyperintense (bright);

2. 2.

   FLAIR images help delineate tumor boundaries, assess tumor infiltration into surrounding brain tissue, and are particularly sensitive to peritumoral edema, which appears hyperintense (bright) on FLAIR sequences;

3. 3.

   T2 images help delineate tumor extent, identify peritumoral edema, and assess the relationship between the tumor and surrounding brain structures;

4. 4.

   Tumors typically appear hypointense (dark) on T1 images, while the contrast between the tumor and surrounding normal brain tissue may not always be sufficient for accurate segmentation.

Figure 10 shows the learned reliability coefficients $\beta_{k}^{t}$ estimated by MMEF-nnFormer, for the four modalities and the three tumor classes. It can be seen that the evidence from the T1Gd modality is reliable when the true class is ED, ET, or NRC/NET, with all the reliability values greater than 0.9. In contrast, the evidence from the FLAIR modality is more reliable for the ED class with a high-reliability coefficient of 0.879 against, respectively, 0.26 and 0.39 for ET and NRC/NET. The evidence from the T2 modality shows similar reliability in segmenting the three classes with a reliability coefficient of around 0.5. The evidence from the T1 modality is the least reliable one, compared with the other three MRI modalities. These results are consistent with domain knowledge about these modalities as reported in [5] and recalled at the beginning of this section, i.e., T1Gd images are useful for delineating tumor boundaries, FLAIR images are sensitive to ED, and T1 images are not sufficient for accurate tumor segmentation. This transparency and explainability of the decision-making process can be expected to enhance end-users’ trust and can be seen as significant advantages of the proposed multimodal evidence fusion approach, as opposed to the “black box” nature of conventional deep learning segmentation models.

The main advantage of our framework is its ability to model and learn the reliability of each image modality, which can be crucial when dealing with diverse, potentially noisy, or low-quality data. While multimodal medical image segmentation tasks are the focus of this paper, the proposed deep evidential fusion framework can be applied to a broader range of challenging medical tasks involving heterogeneous data sources. For instance, in medical tasks such as diagnosing dementia or Alzheimer’s disease, various heterogeneous medical data are available [10]. These data can include lower-quality brain MRI images due to brain degeneration, textual data on disease history and progression, time-series data on blood-brain-barrier integrity, cerebrovascular information, and other relevant physiological measures. Traditional models struggle to effectively address this heterogeneous data within a single neural network [13], and recent work also proposed to address data heterogeneity with model ensembles and hard decision fusion [78]. Our deep evidential fusion framework could be well-suited to analyze such heterogeneous medical tasks. By learning the reliability coefficients for each of the modalities, our model can effectively combine the evidence from heterogeneous sources to reach a more informed and explainable diagnostic decision.

Beyond medical image processing, our approach could be applied to multimodal data fusion in other domains, such as reviewed in [44] and [8]. As examples of potential application domains where heterogeneous data need to be processed to make decisions, we can mention remote sensing and earth observations, in which light detection and ranging (LiDAR), synthetic aperture radar (SAR), and hyperspectral images need to be combined for, e.g., improved classification of objects. As noted in [44], SAR and LiDAR use different electromagnetic frequencies and thus interact differently with materials and surfaces. It would thus be beneficial to apply different discounting (reliability) coefficients to these sensor data depending on the nature of the objects of interest. This conjecture needs, of course, to be validated experimentally, which goes beyond the scope of this paper.

Although the operations of DST have, in the worst case, exponential complexity, the mass functions computed in the EM module have only $K$ focal sets, where $K$ is the number of classes, and the contextual discounting operation computed in the MMEF is applied to the contour function, as explained in Section 3.1.3. Consequently, the number of operations performed in the EM and MMEF modules is only linear in the number of classes. More precisely, as shown in [16], each forward and backward propagation for one voxel and one modality in the EM module has complexity $O(I(H+K))$, where $I$ is the number of prototypes, $H$ is the number of features extracted by the FE module, and $K$ is the number of classes. In the MMEF module, the discounting of the $T$ mass functions for each voxel using (9) and their combination using (10) can be performed in $O(KT)$ operations, and the backward pass (gradient calculation) requires the same computational effort. Overall, the complexity of our model is, thus, similar to that performed in standard neural network architectures based on weighted sums. In terms of computing times, pre-training each of the FE modules with the nnFormer architecture took approximately one hour on our machine⁶⁶6All models were trained on an NVIDIA A100-SXM4 graphics card with 40 GB GPU memory. for the BraTS2021 dataset, and training the whole system end-to-end took 2.3 hours. The total training time (6.4 hours) is slightly less than that of nnFormer with the four modalities (7.8 hours). As far as state-of-the-art uncertainty quantification techniques are concerned, Monte Carlo dropout does not significantly impact training time, while the deep ensemble method is notoriously time-consuming because it implies training several models. Overall, our framework based on DST and decision fusion is at least as efficient as alternative uncertainty quantification approaches.

Our approach is based on combining high-level information extracted from each modality by the FE and EM modules in the form of mass functions. It, thus, has all the advantages and limitations of decision-level fusion approaches. On the plus side, it is highly modular and can still provide sensible results when only some of the modalities are available. This advantage is not crucial in multimodal image segmentation applications because all modalities are usually available, but it can matter in other potential applications such as remote sensing, as mentioned above. Another advantage of decision fusion is that the fusion process is simple and transparent, as already discussed in Sections 4.2 and 4.3. On the minus side, decision-level fusion is, at least in principle, suboptimal because it does not consider all input data globally: we can always construct a classification task in which a single classifier trained with a set of features will perform better than a combination of classifiers trained with each of the features. The good performances of our approach reported in Sections 4.2 and 4.3 show that this potential suboptimality is not an issue in the considered medical image segmentation applications, but it could be in other applications. Another limitation of our approach is that, to keep computations simple, we do not combine the whole discounted mass functions in the MMEF module, but only the contour functions. As a result, the output at each voxel is not a full mass function (with $2^{K}-1$ focal sets), which prevents us from harnessing the full power of DST, such as some of the decision rules reviewed in [18]. This and other limitations will be addressed in future work.