{'fraction_non_alphanumeric': 0.030455606041999265, 'fraction_numerical': 0.014531908796921691, 'mean_word_length': 4.824988078206962, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
MRI brain tumor segmentation and uncertainty estimation using 3D-UNet architectures

Automation of brain tumor segmentation in 3D magnetic resonance images (MRIs) is key to assess the diagnostic and treatment of the disease. In recent years, convolutional neural networks (CNNs) have shown improved results in the task. However, high memory consumption is still a problem in 3D-CNNs. Moreover, most methods do not include uncertainty information, which is especially critical in medical diagnosis. This work studies 3D encoder-decoder architectures trained with patch-based techniques to reduce memory consumption and decrease the effect of unbalanced data. The different trained models are then used to create an ensemble that leverages the properties of each model, thus increasing the performance. We also introduce voxel-wise uncertainty information, both epistemic and aleatoric using test-time dropout (TTD) and data-augmentation (TTA) respectively. In addition, a hybrid approach is proposed that helps increase the accuracy of the segmentation. The model and uncertainty estimation measurements proposed in this work have been used in the BraTS'20 Challenge for task 1 and 3 regarding tumor segmentation and uncertainty estimation.


Introduction
Brain tumors are categorized into primary, brain originated; and secondary, tumors that have spread from elsewhere and are known as brain metastasis tumors. Among malignant primary tumors, gliomas are the most common in adults, representing 81% of brain tumors [7]. The World Health Organization (WHO) categorizes gliomas into grades I-IV which can be simplified into two types (1) "low grade gliomas" (LGG), grades I-II, which are less common and are characterized by low blood concentration and slow growth and (2) "high grade gliomas" (HGG), grades III-IV, which have a faster growth rate and aggressiveness. This work has been partially supported by the project MALEGRA TEC2016-75976-R financed by the Spanish Ministerio de Economía y Competitividad.
The extend of the disease is composed of four heterogeneous histological sub-regions, i.e. the peritumoral edematous/invaded tissue, the necrotic core (fluid-filled), the enhancing and non-enhancing tumor (solid) core. Each region is described by varying intensity profiles across MRI modalities (T1-weighted, postcontrast T1-weighted, T2-weighted, and Fluid-Attenuated Inversion Recovery-FLAIR), which reflect the diverse tumor biological properties and are commonly used to assess the diagnosis, treatment and evaluation of the disease. These MRI modalities facilitate tumor analysis, but at the expense of performing manual delineation of the tumor regions which is a challenging and time-consuming process. For this reason, automatic mechanisms for region tumor segmentation have appeared in the last decade thanks to the advancement of deep learning models in computer vision tasks. Despite these recent advances, the segmentation of brain tumors in multimodal MRI scans is still a challenging task in medical image analysis due to the highly heterogeneous appearance and shape of the problem.
The Brain Tumor Segmentation (BraTS) [1][2][3][4][5] challenge started in 2012 with a focus on evaluating state-of-the-art methods for glioma segmentation in multimodal MRI scans. BraTS 2020 training dataset includes 369 cases (293 HGG and 76 LGG), each with four 3D MRI modalities rigidly aligned, re-sampled to 1mm 3 isotropic resolution and skull-stripped with size 240x240x155. Each provides a manual segmentation approved by experienced neuro-radiologists. Training annotations comprise the enhancing tumor (ET, label 4), the peritumoral edema (ED, label 2), and the necrotic and non-enhancing tumor core (NCR/NET, label 1). The nested sub-regions considered for evaluation are: whole tumor WT (label 1, 2, 4), tumor core TC (label 1, 4) and enhancing tumor ET (label 4). The validation set includes 125 cases, with unknown grade nor ground truth annotation. The test set is composed of 166 cases.
The goal of this work is to develop a 3D convolutional neural network (CNN) for brain tumor segmentation from 3D MRIs and provide an uncertainty measure to assess the confidence on the model predictions. The proposed methods are used to participate in BraTS'20 Challenge for tasks 1 and 3, respectively. In task 1, we explore the use of two well-known 3D-CNN for medical imaging -V-Net [6] and 3D-UNet [27]-and apply some modifications to their baselines. With both networks, the usage of sampling techniques is necessary due memory limitations as well as data augmentation to prevent over-fitting. For task 3, the work provides voxel-wise uncertainty measures computed at test time, with global and per sub-region information. Uncertainty is estimated using both epistemic and aleatoric [8] uncertainties using test-time dropout (TTD) [9] and data augmentation, respectively.

Semantic Segmentation
Brain tumor segmentation methods include generative and discriminative approaches. Generative methods try to incorporate prior knowledge and model probabilistic distributions whereas discriminative methods extract features from image representations. This latter approach has thrived in recent years thanks to the advancement in CNNs, as demonstrated in the winners of the previous BraTS. The biggest break through in this area was introduced by DeepMedic [10] a 3D CNN that exploits multi-scale features using parallel pathways and incorporates a fully connected conditional random field (CRF) to remove false positives. [11] compares the performances of three 3D CNN architectures showing the importance of the multi-resolution connections to obtain fine details in the segmentation of tumor sub-regions. More recently, EMMA [12] creates and ensemble at inference time which reduces overfitting but at high computational cost, and [13] proposes a cascade of two CNNs, where the first network produces raw tumor masks and the second network is trained on the vecinity of the tumor to predict tumor regions. BraTS 2018 winner [14] proposed an asymmetrically encoder-decoder architecture with a variational autoencoder to reconstruct the image during training, which is used as a regularizer. Isensee, F [15] uses a regular 3D-U-Net optimized on the evaluation metrics and co-trained with external data. BraTS 2019 winners [16] use a two-stage cascade U-Net trained end-toend. Finally, [17] applies several tricks in three categories: data processing, model devising and optimization process to boost the model performance.

Uncertainty
Uncertainty information of segmentation results is important, specially in medical imaging, to guide the clinical decisions and help understand the reliability of the provided segmentation, hence being able to identify more challenging cases which may require expert review. Segmentation models for brain tumor MRIs tend to label voxels with less confidence in the surrounding tissues of the segmentation targets [19], thus indicating regions that may have been miss-segmented.
Last year's BraTS challenge already started introducing uncertainty measurements. [18] computes epistemic uncertainty using TTD. They obtain a posterior distribution generated after running several epochs for each image at test-time. Then, mean and variance are used to evaluate the model uncertainty. A different approach is proposed by Wang G [19], who uses TTD and data augmentation to estimate the voxel-wise uncertainty by computing the entropy instead of the variance. Finally, [20] proposes to incorporate uncertainty measures during training as they define a loss function that models label noise and uncertainty.

Dataset Statistics
The biggest complexity for brain tumor segmentation is derived from the class imbalance. The tumor regions account for a 5-15% of the brain tissue and each tumor region is an even smaller portion. Fig. 1 provides a graphical representation of the distribution per each tumor class: ET, NCR, ED; without healthy tissue. It can be seen, that ED is more probable than ET and NCR and that there is high variability between subjects in the NCR label. Another complexity is the difference between glioma grades as LGG patients are characterized by low blood concentration which is translated to low appearance of ET voxels and higher number of voxels for NCR and NET regions.

Data Pre-processing and Augmentation
MRI intensity values are not standardized as the data is obtained from different institutions, scanners and protocols. Therefore we normalize each modality of each patient independently to have zero mean and unit std based on non-zero voxels only, which represent the brain region.
We also apply data augmentation techniques to prevent over-fitting by trying to disrupt minimally the data. For this, we apply Random Flip (for all 3 axes) with a 50% probability, Random 90 • Rotation on two axis with a 50% probability, Random Intensity Shift between (−0.1..0.1 of data std) and Random Intensity Scale on all input channels at range (0.9..1.1).

Sampling Strategy
3D-CNNs are computationally expensive and in many cases, the input data cannot be fed directly to the network. Patch-wise training helps to free memory resources so more images can be fed in one batch. However, there is a trade-off between patch size and batch size. Bigger batches will have a more accurate representation of the data but will require smaller patches (due to memory constraints) that provide local information but lack contextual knowledge.
Another key aspect to consider when selecting the patching strategy is to maintain the class distribution. Losing this distribution can generate a biased model, i.e. if the model only sees small patches with tumor it will likely missclassify healthy tissue.
In this work, we have used two approaches depending on the patch size.
-Binary Distribution: Small patches, equal or lower than 64 3 are randomly selected with a 50% probability of being centred on healthy tissue and 50% probability on tumor [10]. -Random Tumor Distribution: Bigger sizes, 112 3 or 128 3 , are selected randomly but always centred in tumor region, as the patches will contain more healthy tissue and background information.

Loss
The Dice score coefficient (DSC) is a measure of overlap widely used to assess segmentation performance when ground truth is available. Proposed in Milletari et al. [6] as a loss function for binary classification, it can be written as: where N is the number of voxels, p i and g i correspond to the predicted and ground truth labels per voxel respectively, and is added to avoid zero division.
Many variations of the dice loss have been proposed in the literature. For instance, the Generalized Dice Loss (GDL) [26] which is based on the generalized dice score (GDS) [28] for multiple class evaluation. Its goal is to correct the correlation between region size and dice score, by weighting the contribution of each label with the inverse of its volume. It is described as: where L represents the number of classes and w l the weight given to each class. We use the GDL variant as it is more suited for unbalanced segmentation problems.

Network Architecture
This work proposes three networks, variations of V-Net [6] and 3U-Net [27] architectures, for brain tumor segmentation and creates an ensemble to mitigate the bias in each independent model.
The different models are trained using the ADAM optimizer, with start learning rate of 1e − 4, decreased by a factor of 5 whenever the validation loss has not improved in the past 30 epochs and regularized with a l2 weight decay of 1e − 5. They all use the GDL loss.

V-Net
The V-Net implementation has been adapted to use four output channels (Non-Tumor, ED, NCR/NET, ET) and uses Instance Normalization [21] in contrast to Batch Normalization, which normalizes across each channel for each training example instead of the whole batch. Also, as proposed in [15], we have increased the number of feature maps to 32 at the highest resolution, instead of 16 as proposed by the original implementation. Figure 2 shows the network architecture with an input patch size of 64x64x64.
The network has been trained using a patch size of 96 3 and the random tumor distribution strategy (see 3.3). The maximum batch size due to memory constraints is 2. Fig. 2: V-Net [6] architecture with instance normalization, PreLU non-linearities, 32 feature channels at the highest resolution. Feature dimensionality is denoted at each block. The network outputs the segmentation and the softmax prediction.

3D-UNet
We use the original implementation with some minor modifications. Batch Normalization is changed for Group Normalization and, as in V-Net, we use 32 feature maps at the highest resolution.
The network architecture is divided into symmetric Encoder and Decoder parts. The Encoder is composed of two convolutional blocks -with 3DConv + ReLu + GroupNorm structure. The downsampling is performed with 2 3 Max-Pooling and the corresponding upsampling is performed with interpolation. All convolutional layers have kernel size 3 3 , except for the last one that has 1x1x1 kernel and 4 feature maps as output. In this case, we use ReLu non-linearity and the skip-connections are joined with a concatenation step. The network outputs a four-channel segmentation map with the training labels as well as a softmax. The detailed architecture can be seen in Figure 3.
The Basic 3D-UNet is trained with a patch size of 112 3 and a batch size of 2.  This network is trained following two different strategies. The first one, 3D-UNet-residual uses a patch size of 112 3 and a batch size of 2 for the whole training, whereas 3D-UNet-residual-multiscale varies the sampling strategy so the network sees local and global information. For that, the first half of the training uses a patch size of 128 3 with a batch size of 1. Then, the patch size is reduced to 112 3 and the batch increased to 2.

Post-Processing
In order to correct the appearance of false positives in the form of small and separated connected components, this work uses a post-processing step that keeps the two biggest connected components if their proportion is bigger than some threshold -obtained by analysing the training set. With this process, small connected components that may be false positives are removed but big enough components are kept as some of the subjects may have several tumors.
Moreover, one of the biggest difficulties of this challenge is to provide an accurate segmentation of the smallest sub-region, ET, which is particularly difficult to segment in LGG patients, as almost 40% have no enhancing tumor in the training set. In the evaluation step, BraTS awards a Dice score of 1 if a label is absent in both the ground truth and the prediction. Conversely, only a single false positive voxel in a patient where no enhancing tumor is present in the ground truth will result in a Dice score of 0. Therefore, some previous works [15,16] propose to replace enhancing tumor voxels for necrosis if the total number of enhancing voxels is smaller than some threshold, which is found for each experiment independently. However, we were not able to find a threshold that improved the performance as it helped for some subjects but made some other results worse.

Uncertainty
This year's BraTS includes a third task to evaluate the model uncertainty and reward methods with predictions that are: (a) confident when correct and (b) uncertain when incorrect. In this work, we model the voxel-wise uncertainty of our method at test time, using test time dropout (TTD) and test-time data augmentation (TTA) for epistemic and aleatoric uncertainty respectively.
We compute epistemic uncertainty as proposed in Gal et.al [23], who uses dropout as a Bayesian Approximation in order to simplify the task. Therefore, the idea is to use dropout both at training and testing time. The paper suggests to repeat the prediction a few hundred times with random dropout. Then, the final prediction is the average of all estimations and the uncertainty is modelled by computing the variance of the predictions. In this work, we perform B = 20 iterations and use dropout with a 50% probability to zero out a channel. The uncertainty map is estimated with the variance for each sub-region independently.
be the vector that represents the i-th voxel's predicted labels, the voxel-wise uncertainty map, for each evaluation region, is obtained as the variance: Uncertainty can also be estimated with the entropy, as [19] showed. However, the entropy will provide a global measure instead of map for each sub-region. In this case, the voxel-wise uncertainty is calculated as: wherep i m is the frequency of the m-th unique value in Y i and X represents the input image.
To model aleatoric uncertainty we apply the same augmentation techniques from the training step plus random Gaussian noise, in order to add modifications not previously seen by the network. The final prediction and uncertainty maps are computed following the same strategies as in the epistemic uncertainty.
All that begin said, we hope to evaluate the model's behaviour w.r.t to input and model variability by defining the several experiments:

Results
The code 1 has been implemented in Pytorch [24] and trained on the GPI 2 servers, based on 2 Intel(R) Xeon(R) @ 2.40GHz CPUs using 16GB RAM and a 12GB NVIDIA GPU, using BraTS 2020 training dataset. We report results on training, validation and test datasets. All results, prediction and uncertainty maps, are uploaded to the CBICA's Image Processing Portal (IPP) for evaluation of Dice score, Hausdorff distance (95th percentile), sensitivity and specificity per each class. Specific uncertainty evaluation metrics are the ratio of filtered TN (FTN) and the ratio of filtered TP (FTP).

Segmentation
The principal metrics to evaluate the segmentation performance are the Dice Score, which is an overlap measure for pairwise comparison of segmentation mask X and ground truth Y: and the Hausdorff distance, which is the maximum distance of a set to the nearest point in the other set, defined as: where sup represents the supremum and inf the infimum. In order to have more robust results and to avoid issues with noisy segmentation, the evaluation scheme uses the 95th percentile. Tables 1 and 2 show Dice and Hausdorff Distance (95th percentil) scores for training and validation sets respectively.   The model used with the test set is the Residual 3D-UNet-multiscale with post-processing. Table 3 shows the results in the training, validation and test sets for comparison. All the proposed models are greatly penalized when no ET is present on the ground truth. In addition, the V-Net suffers more from false positives and 3D-UNet based models from false negatives. The excess of false positives may be caused due to the usage of small patches instead of using the whole volume, which provokes a variation in the proportion of healthy tissue against tumor regions. On the other hand, 3D-UNet models use bigger patch sizes and pooling layers instead of strided convolutions which may be the cause of having a larger number of false negatives. Increasing the patch size helps reduce false positives but it misses local information, which is reflected in label miss-classification on the region's boundaries. Figure 5 shows a visual comparison of the models with a representation on the explained behaviours.

Uncertainty
BraTS requires to upload three uncertainty maps, one for each subregion (WT, TC, ET) together with the prediction map. Values must be normalized between 0-100 such that "0" represents the most certain prediction and "100" represents the most uncertain. The metrics used are the FTP ratio defined as F T P = (T P 100 − T P T )/T P 100 , where T represents the threshold used to filter the more uncertain values. The ratio of filtered true negatives (FTN) is calculated in a similar manner. The integrated score will be calculated as follows: From this point forward all experiments are performed on the model Residual 3D-UNet-multiscale, as it is the one with more balanced results across the different regions. Table 4 shows the results for the epistemic, aleatoric and hybrid uncertainties when computed with entropy or variance. As a general overview, we can see that the AUC-Dice, which is computed by averaging the segmentation results for several thresholds that filter uncertain predictions, improves 2 to 3 points w.r.t the results obtained in the segmentation task (W T : 0.8172, T C : 0.7664, ET : 0.7071). Although the metrics are not the same, it indicates that the model is more certain on the TP and less certain on FP and FN. Moreover, the AUC-Dice is higher when using entropy as the uncertainty measure.
Our results show that the model is more uncertain in LGG patients, particularly on epistemic uncertainty; meaning the model requires more data to achieve a more confident prediction. If we compare the behaviour between the uncertainty types, we see that (1) aleatoric focuses on the region boundaries, with small variations (2) epistemic improves results on the ET region but filters more TP and TN and(2) the hybrid approach achieves the best Dice-AUC results when using entropy as the uncertainty measurement.  We participate in the challenge using test time augmentation (TTA) when uncertain values are computed using variance as it achieves the highest integrated scores per sub-region on the validation set. Table 5 shows the obtained results in both validation and test sets. The achieved integrated scores for validation are 0.93, 0.91 and 0.89 and for test 0.93, 0.93, 0.91 for WT, TC and ET respectively. We see a two point improvement on the ET and TC sub-regions for the test set.

Discussion and Conclusions
This work proposes a set of models based on two 3D-CNNs specialized in medical imaging, V-Net and 3D-UNet. As each of the trained models performs better in a particular tumor region, we define an ensemble of those models in order to increase the performance. Moreover, we analyze the implication of uncertainty estimation on the predicted segmentation in order to understand the reliability of the provided segmentation and identify challenging cases, but also as a means of improving the model accuracy by filtering uncertain voxels that should refer to wrong predictions. We use the Residual 3D-UNet-multiscale as our model to participate at the BraTS'20 challenge. The best results in the validation set are obtained when creating an ensemble of the proposed models, as we can leverage the biases of each model, but are still far from the current state the art. These results may be caused by a bad training strategy where the sampling technique does not reflect the correct label distribution, thus providing more false detections. This is reflected more in the ET region as all models predict more tumor voxels of this label, which is greatly penalized when the ground truth does not contain it. In order to improve results, future work should try to provide a better representation of the labels, not just increase the patch size, but maybe let the network see both local and more global information.
Another potential problem is the model's simplicity. Although previous works achieve good results using a 3D-UNet, i.e. [15], adding more complexity to the network may help boost the performance. Therefore a possible line of work would be to extend the proposed models into a cascaded network, where each nested evaluation region -WT, TC and ET-is learnt as a binary problem. Also, LGG subjects usually achieve lower accuracy on the prediction. In order to improve the results, we could research other post processing techniques and design them specifically to target each one of the glioma grades, as they may be differentiated by the sub-region distribution.
For uncertainty estimation, the work evaluates the usage of aleatoric, epistemic and a hybrid approach using the entropy as a global measure and variance to evaluate uncertainty on each evaluation region. In the provided results, it has been seen that using uncertainty information actually helps improve the accuracy of the network, achieving the best Dice Score (AUC, estimated from filtering uncertain voxels) when using the hybrid approach and entropy as the uncertainty measure. Our method achieves a score of 0.93, 0.93, 0.91 for WT, TC and ET respectively on the test set.