{'fraction_non_alphanumeric': 0.04250992306561474, 'fraction_numerical': 0.013622776498260401, 'mean_word_length': 4.604146642866951, 'pattern_counts': {'":': 0, '<': 1, '<?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': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perturb-and-MPM: Quantifying Segmentation Uncertainty in Dense Multi-Label CRFs

This paper proposes a novel approach for uncertainty quantification in dense Conditional Random Fields (CRFs). The presented approach, called Perturb-and-MPM, enables efficient, approximate sampling from dense multi-label CRFs via random perturbations. An analytic error analysis was performed which identified the main cause of approximation error as well as showed that the error is bounded. Spatial uncertainty maps can be derived from the Perturb-and-MPM model, which can be used to visualize uncertainty in image segmentation results. The method is validated on synthetic and clinical Magnetic Resonance Imaging data. The effectiveness of the approach is demonstrated on the challenging problem of segmenting the tumor core in glioblastoma. We found that areas of high uncertainty correspond well to wrongly segmented image regions. Furthermore, we demonstrate the potential use of uncertainty maps to refine imaging biomarkers in the case of extent of resection and residual tumor volume in brain tumor patients.


Introduction
Image segmentation plays a pivotal role in the analysis of medical imaging data. Information extracted via segmentation (e.g. volume/position) can be used for diagnosis, treatment planning and monitoring. Graphical models offer a way to describe structure in medical images and embed it in a probabilistic framework. Markov or Conditional Random Fields (CRFs) have become one of the most widely-used graphical models in image understanding (Wang et al., 2013). Although CRFs are probabilistic graphical models (they model a Gibbs distribution), the final image segmentation corresponds usually to the most probable (MAP) labeling (=hard labeling). Sampling from the underlying Gibbs distribution would require the use of Markov Chain Monte Carlo (MCMC) methods, which are computationally very expensive and thus prohibitive for clinical applications. Efficient sampling would enable quantification of segmentation uncertainty, which could better support the decision making process of clinicians employing the segmentations. More specifically, image regions that are segmented (automatically) with high uncertainty can be e.g. reassessed by a human expert or excluded from subsequent analysis.
The most popular CRF model corresponds to a structured grid, where voxels are represented by nodes in a graph. For neighboring nodes pairwise potentials are defined which aim to induce smoothness of the image segmentation. The use of a Potts prior in grid-structured pairwise CRFs is known to result in excessive smoothing (also known as shrinking bias (Kohli et al., 2009;Meier et al., 2016b)) of the object boundary. Dense CRFs establish pairwise potentials between all nodes in a graph allowing for long-range interactions among voxels. Consequently, the shrinking bias of grid-structured pairwise CRFs is reduced allowing for more detailed segmentation results. In two key papers, Krähenbühl et al. proposed an efficient inference (Krähenbühl and Koltun, 2011) and learning algorithm (Krähenbühl and Koltun, 2013) for dense CRFs, making their application in medical image volumes feasible. As a consequence, dense CRFs are becoming increasingly popular for a wide range of segmentation problems in medical image analysis including brain tumor segmentation (Kamnitsas et al., 2016), lung segmentation (Gao et al., 2016), retinal vessel segmentation (Orlando and Blaschko, 2014;Orlando et al., 2016) or liver segmentation (Christ et al., 2016). However, the assessment of segmentation uncertainty in dense CRFs has not been possible so far.
The need for uncertainty quantification in segmentation has recently been addressed by Lê et al. (2016), who proposed a Gaussian Process for sampling candidate segmentation boundaries close to a reference segmentation. Although applicable to multiple label maps simultaneously, their approach works on contours and is thus primarily designed for binary segmentation problems. In addition, the method is also decoupled from the segmentation method that generated the reference segmentation in the first place, and it is driven by the fuzzyness of image boundaries. Related to our work are also Bayesian approaches that rely on parameter sampling using MCMC techniques to quantify segmentation uncertainty (e.g. (Iglesias et al., 2013;Le Folgoc et al., 2017)). In contrast to such methods, we do not aim to sample model parameter values but aim to rather perturb a fixed parameter value in a principled fashion. Recently, Papandreou et al. proposed random MAP perturbations (Perturb-and-MAP (Papandreou and Yuille, 2011)) for CRFs that effectively allow drawing samples from the underlying Gibbs distribution. Alberts et al. (2016b) employed Perturb-and-MAP within a grid-structured CRF for segmenting brain tumors and quantifying the uncertainty of volume estimates. In contrast to (Alberts et al., 2016b;Lê et al., 2016), we focus on multi-label segmentation. More importantly, Perturband-MAP requires an exact MAP estimation, which is not feasible for dense CRFs in medical images due to their size. Moreover, in multi-label segmentation problems the MAP solution can usually only be approximated. Despite these obstacles, we propose a novel perturbation model, referred hereafter as Perturb-and-MPM, for dense multi-label CRFs and investigate its feasibility for quantifying segmentation uncertainty.
Recently, Kim et al. (Kim and Ermon, 2016) proposed a method that enables approximate sampling from CRFs for which MAP inference can be reformulated as an Integer Linear Program (ILP). Similarly to our approach, they also rely on random perturbations in order to realize the sampling. However, their approach requires to solve the respective ILP, which in case of a dense multi-label CRF would rely on a very large number of constraints and thus result in a running time prohibitive for clinical use.
The contribution of this paper is a novel approach to quantify segmentation uncertainty in dense multi-label CRFs, which is computationally efficient and easy to implement. The method is validated on synthetic data and on Magnetic Resonance (MR) Imaging datasets of brain tumor patients. A clinical application to the challenging problem of estimating the extent of resection and residual tumor volume in postoperative images is presented and serves as an example on how uncertainty estimation can improve the extraction of volumetric imaging biomarkers.

Preliminaries
This section introduces the notation and concepts this work is based on.

Notation
In the following, the input image is denoted by I. The set of voxels in I is denoted by V and the total number of voxels by N . The available label values are contained in the set L = {1, . . . , m}. A labeling of I is referred to by X = {x i } i∈V with x i being a scalar value that indicates the segmented image region (e.g. tissue compartment), i.e. x i ∈ L.

Conditional Random Field
In a CRF model, voxels of an image are represented by nodes in a graph. Hence, we interpret the image I as an undirected graph G consisting of nodes and edges. Every node in G is associated with a random variable x i . The pair (I, X) is a CRF, if for any given I the distribution P (X|I) factorizes according to G. The conditional distribution corresponds to a Gibbs distribution P (X|I) = 1 Z(I) exp (−E(X|I)), where Z(I) is the partition function. The most probable (MAP) labeling for a given image I can then be estimated by minimizing the energy, i.e. (1) The energy of a pairwise CRF corresponds to a sum of unary potentials ψ u (=data-likelihood) and pairwise potentials ψ p (=prior). In contrast to gridstructured CRFs, pairwise potentials in dense CRFs are defined between all pairs of voxels, i.e. G is a complete graph.

Dense CRF and MPM-inference
Krähenbühl et al. (Krähenbühl and Koltun, 2011) proposed an efficient approximate Maximum Posterior Marginal (MPM) inference scheme based on a mean field approximation and cross bilateral filtering techniques. In short, a mean field with the approximate distribution Q(X) = i Q i (x i ) is introduced. The mean field approximation algorithm aims at minimizing the KL-Divergence between the approximate distribution Q and the true distribution P . The algorithm then iteratively performs updates of the marginals Q i (x i ) via an efficient filtering scheme and estimates the MPM solution via x M P M = arg max l Q i (x i = l) for all voxels in the image (where l ∈ L). The update equation for a marginal Q i is given by: where for notational convenience we define the energy as and θ refers to both ψ u and ψ p . The marginal probabilities of the original Gibbs distribution are approximated by:

Perturb-and-MAP
The main idea of Perturb-and-MAP (Papandreou and Yuille, 2011) is to locally perturb potentials of the energy function of the CRF by adding independent and identically distributed (i.i.d.) random noise γ. Subsequently, the MAP labeling of the perturbed energy can be estimated as: If the perturbation density follows a Gumbel distribution, it can be shown that the Perturb-and-MAP model approximates the Gibbs distribution of the corresponding random field (Tomczak, 2016). More specifically, we consider γ(X) to be drawn i.i.d. from a Gumbel distribution with zero mean and cumulative distribution function F (t) = exp(− exp(−(t + c))), with c being the Euler constant. Under perturbation of all potentials (cf. (Papandreou and Yuille, 2011), Proposition 3), the MAP labeling of the perturbed energy coincides with the corresponding sample of the Gibbs distribution. In this case, the distribution of the perturbed MAP solutions (see Equation (4)) is equivalent to the original Gibbs distribution: We denote the LHS of Equation (5) as Perturb-and-MAP model. This relationship is deeply connected to the so-called "Gumbel-max-trick" (Yellott, 1977) which is summarized in the following Lemma (Papandreou and Yuille, 2011): Lemma 1. Let {θ 1 , . . . , θ m }, with θ j ∈ R. We additively perturb them byθ j = θ j + γ j with γ j i.i.d. Gumbel samples. Then the probability thatθ j attains the minimum value is P (arg min{θ 1 , . . . ,θ m } = j) = exp(−θ j )/ m j =1 exp(−θ j ).

Estimation of marginal probability distributions
We are interested in the marginal probability of a voxel i being labeled with a particular label l, i.e. the event x i = l with l ∈ L. The corresponding marginal probability distribution can be expressed as an expectation (1 {·} denotes the indicator function): The expectation can be estimated via sample approximations. Consider S = x (1) i , · · · ,x (m) i to be |S| samples of possible labels for voxel i. The marginal probability for x i = l can then simply be estimated by This corresponds to the approximation of P (x i ) with the relative frequencies of its events, i.e. via the respective histogram.

Error measures
In the following, we rely on the voxel-wise hamming loss for quantifying differences between two image labelingsX and X: Moreover, for measuring the similarity between two probability distributions, we introduce the (voxel-wise) total variation distance:

Perturb-and-MPM
In this section, we introduce the main contribution of our work, which is a perturbation-based sampling approach for dense multi-label CRFs. We start by employing a random field over the random variables X = {x 1 , . . . , x N } conditioned on an image I, where every voxel i = 1, . . . , N is associated with a random variable x i . The random variables take values from the previously defined label set L = {1, 2, . . . , m}. As presented in Section 2, the dense CRF defines a Gibbs distribution P (X|I) with energy E (X|I) composed of unary potentials ψ u (x i ) and pairwise potentials ψ p (x i , x j ), which are collectively referred to as θ.

Perturbation Model
In the following, we adapt the perturbation strategy proposed by Papandreou et al. (Papandreou and Yuille, 2011), to enable sampling from dense multi-label CRFs. The adaptation considers three steps: First, the potentials θ of the initial energy function are perturbed with i.i.d. Gumbel noise γ. Second, the mean field approximation is performed. Third, the perturbed (voxel-wise) The Perturb-and-MPM model f M P M can be constructed by iterating the previous three steps and aggregating the MPM solutions. This enables the estimation of the marginal distribution f M P M (x i ) via the histogram of the aggregated MPM solutions (see Section 2.5). The maximum number of samples (i.e. aggregated MPM solutions) is denoted by T . Note that f M P M (x i = l) is used as a shorthand for P (x M P M,i = l) and they can be used interchangeably. The main idea is visualized in Figure 1 and the procedure is summarized in Algorithm 1. Perturb potentials of E(X|I; θ) →θ (t) := θ + γ (t)

Distribution of MPM-solutions
A natural question that arises from the previous perturbation model is how the perturbed MPM solutions are distributed. This aspect is crucial since our estimation of segmentation uncertainty will be based on this distribution. Given Figure 1: Conceptual overview of Perturb-and-MPM. The input to Perturb-and-MPM is an image represented by a dense CRF and its energy function E(X|I), respectively. Sampling is realized by first perturbing the potentials θ. Second, by computing the mean field approximation Q(X;θ) and third by estimating the MPM solution of the perturbed mean field. After iterating these three steps, the sampled labelmapsX L, one out of m label values can be assigned to a variable x i . We proceed by reparametrizing Equation (3) to explicitly take all possible label assignments into account (cf. fully expanded potential table of (Papandreou and Yuille, 2011) Now, we can state the following proposition: Proposition 1. If we perturb each entry of the fully expanded potential vector θ with i.i.d. Gumbel noise samples γ j , j = 1, · · · , m, then the marginals of the Perturb-and-MPM and mean field model coincide, i.e., Proof. The proof is analogous to the proof of Proposition 3 in (Papandreou and Yuille, 2011) and makes use of the reparametrization in Equation (10) and Lemma 1.
The consequence of Proposition 1 is that the Perturb-and-MPM model is equivalent to the original mean field approximation. However, Proposition 1 does not hold for the Algorithm 1, because the mean field approximation Q i (x) will be different in every iteration (i.e. for every perturbation, cf. step 4&5 in Algorithm 1). This observation plays a central role in the proposed Perturband-MPM framework. It allows the generation of marginal distributions that are possibly closer to the marginals of the original Gibbs distribution P (x i ) than mean field marginals Q i (x i ). From Equation (5), we know that the distribution of perturbed MAP-solutions is equivalent to the original Gibbs distribution, i.e. P (X M AP = X) = P (X) withX M AP = arg max X −E(X;θ). If we assume that our MPM solution coincides with the MAP, i.e. arg max X Q(X;θ) = arg max X −E(X;θ) holds, then we can state that P (X M P M = X) = P (X M AP = X) = P (X). This requirement of coinciding maxima, illustrated in Figure 2, is much less restrictive than requiring a perfect mean field approximation (i.e. KL-Divergence=0). Mean field approximations have shown to suffer from bad local optima (Hoffman and Blei, 2015). Thus, we anticipate that f M P M (X) may potentially yield marginals that are closer to the true distribution than the mean field approximation.

Deriving error bounds of Perturb-and-MPM
In the previous section, we demonstrated that the Perturb-and-MPM model can generate marginals different from the mean field marginals Q i (x i ) and hence possibly better approximate the marginals of the original Gibbs distribution than a mean field approximation could do. The aim of this section is to perform an error analysis on the derived Perturb-and-MPM model. We are interested in obtaining bounds on the approximation error of Perturb-and-MPM, as this serves as a basis for deriving additional bounds on the error of the final uncertainty quantification. The Perturb-and-MAP and Perturb-and-MPM solutions are given byX M AP = arg max X −E(X;θ) andX M P M = arg max X Q(X;θ), respectively. In order to effectively tackle the analysis, we differentiate between two different situations: i)X M AP =X M P M and ii)X M AP =X M P M . The corresponding loss function is the voxel-wise hamming loss as defined in Equation the identity P (X) = P (X M AP = X) holds. We have f M P M (X) = P (X M P M = X) withX M P M = arg max X Q(X;θ). If we assume that our MPM solution coincides with the MAP (arg max X Q(X;θ) = arg max X −E(X;θ)), we can then state that f M P M (X) = P (X M AP = X) = P (X). In other words, the Perturb-and-MPM model coincides with the original Gibbs distribution. The model f M P M (X) is estimated based on a finite number of samples |S|. We are thus interested i) in knowing if the respective approximation error is bounded and ii) in how many samples are necessary to achieve this bound.
We define the marginals of the Perturb-and-MPM model to be As presented in Section 2.5, the expectation can be estimated via sample approximation: Note that these voxel-wise marginals are computed independently from each other as the mean over a set of samples S. We denote the previous sample approximation byf M P M (x i = l). This sum approximates the true but unknown Perturb-and-MPM marginal f M P M . Hence, we refer to it as the empirical Perturb-and-MPM model. Furthermore, 1 {x M P M,i = l} is a Bernoulli random variable and thus the estimatef M P M (x i = l) defines |S| independent Bernoulli trials, each with success probability f M P M (x i = l). Consequently, we can bound the absolute error between the empirical and the true but unknown Perturb-and-MPM marginal via the Hoeffding bound: If we choose δ such that 2 exp(−2 |S| 2 ) ≤ δ, we can say that with a probability of at least 1 − δ it holds that: In this situation (givenX M AP =X M P M holds) the approximation error of the empirical Perturb-and-MPM model with respect to the original Gibbs distribution is bounded: The first equality follows from the previous reasoning that ifX M  Note that this holds for any voxel i in the image. For fixed ( , δ) the required sample size to obtain a reliable Perturb-and-MPM model is then given by: The previous bound in Equation (16) is limited since it holds only for a fixed label l. Thus, we generalize it in order to know how likely it is that any of the empirical marginalsf M P M (x i ) significantly deviates from P (x i ). This can be described by the probability We proceed as follows: Based on this, we can refine the required sample size in Equation (22) and obtain Therefore, the required number of samples |S| to obtain a reliable (in the sense of ( , δ)) Perturb-and-MPM model, grows logarithmically with the number of labels |L|.
3.3.2. SituationX M AP =X M P M Second, we are looking at the situation whereX M AP =X M P M . In this situation, the previously derived approximation error (Equation (21)) is not guaranteed to shrink to zero anymore but to a residual error. We can characterize this error by investigating the loss that is induced through the different labelingsX M AP =X M P M . We can state the following Proposition: with · 1 being the (voxel-wise) total variation distance.
Proof. We assume for the Perturb-and-MAP (f M AP ) and Perturb-and-MPM (f M P M ) model that for both models the same number of samples |S| was generated. Furthermore, we assume a perturbation of all potentials such that P (X) = f M AP (X) holds.
This shows that the approximation error is indeed bounded. For < |V | the bound is tighter than the trivial 0 ≤ P (X) −f M P M (X) 1 ≤ |V |. As a consequence, the voxel-wise error can be bounded in the same fashion 0 ≤ P (x i ) −f M P M (x i ) 1 ≤ i with i = 1 being the worst case.
In summary, we found for the situation whereX M AP =X M P M holds: • The approximation error P (X) −f M P M (X) 1 is bounded and will go to zero with an increasing number of samples |S|.
• The required number of samples for a reliable Perturb-and-MPM model f M P M (X) scales logarithmically with the number of labels |L|, i.e. |S| ∝ log(|L|).
For the situation whereX M AP =X M P M holds, we found: • Given the hamming loss function (X M AP ,X M P M ), the approximation error P (X) − f M P M (X) 1 can be bounded.

Uncertainty Quantification
Based on Algorithm 1, we can effectively generate the Perturb-and-MPM model f M P M . The uncertainty of the segmentation model for a possible labeling X can then be derived from the MPM-marginals f M P M (x i ). The marginals reflect voxel-wise histograms. A natural choice to quantify uncertainty is via the Shannon-entropy contained in these histograms. The entropy for a discrete probability distribution P over possible labels l ∈ L is given by: H P (x i ) = − l∈L P (x i = l) log 2 P (x i = l). The entropy is maximal for uniformly distributed labels. We can now use the previously derived bound in Equation (31) and results from information theory (Equation (11) in Sason et al. (Sason, 2013)) to provide an upper bound on the difference in entropy of the Perturband-MPM model and the original Gibbs distribution: with h denoting the binary entropy function. For the worst case scenario of i = 1 the bound reduces to: As an example, if we set |L| = 4, we obtain the value 1.585 for the RHS in Equation (43). Hence, this bound is tighter than the trivial bound |H P (x i ) − H f M P M (x i )| ≤ log 2 (|L|), which in this case is log 2 (|L| = 4) = 2.

Implementation of Perturb-and-MPM
Sampling from the zero mean, unit variance Gumbel distribution can be realized via drawing samples u from the standard uniform distribution U (0, 1) and transforming them according to − log(− log(u)). Note that the number of perturbations for perturbing all potentials of the initial energy function grows exponentially with the number of nodes N in the CRF. Recent studies (Hazan et al., 2013;Papandreou and Yuille, 2011) have empirically shown that it is sufficient to perform low-order Gumbel perturbations rather than perturbing all realizations of both unary and pairwise potentials. Hence, we employ an Order-1 perturbation which perturbs each of the unary potentials with Gumbel noise.

Validation on synthetic data
In order to empirically validate the proposed Perturb-and-MPM scheme, we employed a dense CRF defined over grids with different number of nodes (N ∈ {6, 9, 12}) which can take on binary labels, i.e. x i ∈ {0, 1}. The unary potentials were defined to be ψ u = −log(p u ), where p u (x i = 0) ∼ U [0, 1) and p u (x i = 1) = 1 − p u (x i = 0). The pairwise potentials contained a Potts prior and a gaussian kernel over the node positions. For every grid size, we estimated f M P M with a varying number of samples (|S| ∈ 10, 50, 10 2 , 10 3 , 10 4 , 10 5 , 10 6 . We performed a comparison of unperturbed mean field marginals Q(x i ) with marginals of f M P M and exact marginals (estimated with an exact sampler of the R package "CRF" (Wu, 2014)) over 20 random initializations of ψ u (a larger number of initializations did not change the result). Similarity of the unperturbed mean field Q(x i ) and f M P M marginals with exact marginals was quantified using the 1 -norm. The corresponding error plot is shown in Figure 3. For an increasing number of samples |S| the mean (over the 20 initializations) of the approximation error converges to a residual error and from |S| = 50 on appears to be consistently lower than the mean approximation error of the unperturbed mean field (clearly visible for N = 9, N = 12).

Assessing multi-label segmentation uncertainty in glioblastoma patient images
For evaluating the clinical utility of our uncertainty estimation, we looked at the task of multi-label segmentation of glioblastoma from multi-sequential Magnetic Resonance Imaging (MRI) data. In a first application, we look at the most challenging subtask in glioblastoma segmentation, which is the segmentation of the tumor core, including the identification of necrosis, non-enhancing and enhancing tumor. The differentiation of edema from non-enhancing and enhancing tumor is connected to a high uncertainty for human raters, which in turn is reflected in the rather poor performance of supervised learning-based segmentation methods for this task (Menze et al., 2015). In this application, we aim to quantify the quality of the uncertainty prediction. In a second application, we are computing the extent of resection and residual tumor volume from the segmentation results of pre-and postoperative images. The extent of resection as well as the residual tumor volume are important clinically-relevant indicators of success of the neurosurgical tumor resection, as well as prognostic biomarkers in glioblastoma patients (Grabowski et al., 2014). The aim is to study the use of uncertainty quantification in order to improve their extraction. The dataset used for the first application contains preoperative images of 14 patients, while the second application contains pre-as well as immediate postoperative images (acquired within 72 hours after surgery) of 19 patients. All cases are newly diagnosed and histologically confirmed glioblastoma. For the second dataset of 19 patients, nine patients underwent subtotal extirpations/partial resection of enhancing tumor (PRET) and ten patients underwent complete resections of enhancing tumor (CRET).
For both datasets the Magnetic Resonance protocol encompassed T 1 -weighted, T 1 -weighted gadolinium-enhanced, T 2 -weighted and FLAIR-weighted sequences. The FLAIR-sequence was an axial acquisition with an anisotropic voxel size of 1mm × 1mm × 3mm, whereas the remaining sequences were sagittal acquisitions with an isotropic voxel size of 1mm × 1mm × 1mm. The field of view for all sequences was 256 × 256mm 2 . Manual segmentation of all tumor tissues for the first dataset was performed by one expert rater (Neuroradiologist with more than six years of experience). For the second dataset, segmentation of the enhancing tumor was performed by four expert raters. These manual segmentations were used as groundtruth data for the estimation of the extent of resection. More details on the MR acquisition, the segmentation protocol as well as the inter-observer variability of the individual raters can be found in a related clinical study (Meier et al., 2017).
For performing the experiments, we employed a recently proposed and clinicallyevaluated segmentation method ( (Meier et al., 2016b,a)). The method relies on a decision forest classifier for definition of unary potentials in a dense CRF. Details on the feature extraction, decision forest and CRF can be found in (Meier et al., 2014(Meier et al., , 2016b. Parameters of the CRF were learned by minimizing the Intersection-Over-Union loss according to (Krähenbühl and Koltun, 2013). We remark that the training dataset of this model is an independent one from the datasets used for evaluation here, and encompassed 54 patient cases, described in more detail in (Meier et al., 2016a).

Visualization of segmentation uncertainty
The segmentation uncertainty is estimated as explained in Section 3.4 and can be overlayed on the corresponding patient image in form of a heatmap, where voxel-wise uncertainty values are color-coded. In Figure 4, the impact of an increasing number of samples to estimate f M P M and subsequently the uncertainty is visualized. With an increasing number of samples the segmentation uncertainty estimation appears spatially more homogenous. We observed that the uncertainty estimation starts to saturate from |S| = 100 on (implications of these results are discussed in detail in Section 5).

4.2.2.
Quantifying the quality of the segmentation uncertainty in multi-label tumor core segmentation This experiment concerned the segmentation of the tumor core, which involves the differentiation of necrosis, enhancing and non-enhancing tumor from the surrounding edema. The aim of the experiment was to quantify the quality of the uncertainty estimation. In order to do so, we performed a comparison between the predicted segmentation, uncertainty and the manual ground truth data. Four situations were of interest: • Voxels that are misclassified and uncertain (=True Positives).
• Voxels that are misclassified and certain (=False Negative).
• Voxels that are correctly classified and uncertain (=False Positive).
• Voxels that are correctly classified and certain (=True Negative).
A total of |S| = 200 samples were generated in order to estimate segmentation uncertainty. For simplicity, we defined voxels to be uncertain if they showed a non-zero value in the uncertainty estimation. Sensitivity and specificity were assessed accordingly on the preoperative dataset containing 14 patients. The mean and standard deviation of the sensitivity and specificity (mean±sd) for the tumor core is (0.45±0.23) and (0.76±0.13). The values for the enhancing tumor are (0.53 ± 0.1) and (0.88 ± 0.14), respectively. By introducing a threshold value for exclusion of voxels with small uncertainty values, we obtain better specificity with the cost of having a reduced sensitivity. As an example a threshold for the uncertainty of 0.1 yields for the tumor core a sensitivity and specificity of (0.39 ± 0.21) and (0.81 ± 0.12), respectively, and for the enhancing tumor (0.44 ± 0.09) and (0.93 ± 0.11). Likewise, larger threshold values resulted in an increased specificity but lowered sensitivity.

Extraction of prognostic biomarkers in glioblastoma patients
In order to investigate the clinical utility of the uncertainty estimation, we looked at the challenging problem of estimating the extent of resection (EOR) and residual tumor volume (RTV). The extent of resection is defined by the preoperative enhancing tumor volume V pre and the residual enhancing tumor volume V post after surgery. It corresponds to the ratio: An EOR of 1 (or 100%) reflects a complete resection of the enhancing tumor. Pre-and postoperative images of the second dataset containing 19 patients were segmented automatically by our approach (Meier et al., 2016b). The automatically estimated volumes of V pre and V post (=RTV) were used to compute EOR. The automatically estimated EOR and RTV were compared to the groundtruth data of four different raters. In a second step, we neglect all voxels that show a non-zero value in the uncertainty estimation for computing the EOR and RTV, respectively. |S| = 200 samples were generated to estimate segmentation uncertainty. The "corrected" EOR and RTV were then again compared to the groundtruth data. The resulting absolute error in EOR and RTV for both corrected and original estimation are presented as boxplots in Figures 5 and 6. We performed a statistical analysis under a significance level of α = 0.05. The change in error between the original and corrected automatic estimates were assessed using a paired, exact Wilcoxon signed-rank test (R package "exac-tRankTests", Bonferroni correction was applied for multiple comparisons). We found a statistically significant decrease in error of the corrected automatic estimates for the RTV including the comparisons to rater 1 (p = 5.34 × 10 −5 ), rater 2 (p = 6.45 × 10 −4 ) and rater 4 (p = 1.64 × 10 −4 ). An exemplary uncertainty quantification of the segmentation for a postoperative case with a low EOR is shown in Figure 7. Similarly, exemplary uncertainty maps for two different slices of a postoperative patient case with completely resected enhancing tumor (=high EOR) are shown in Figure 8.

Discussion & Conclusion
We have introduced a perturbation-based approach for performing efficient, approximate sampling in dense multi-label CRFs. The methodology enables the approximation of segmentation uncertainty contained in image labelings from dense CRFs. In medical image analysis, image segmentation is used to provide spatial information (volume/position) on anatomical regions of interest. In brain tumor image analysis, segmentation is used to partition glioblastoma into different tumor compartments. The utility of the spatial information of these tumor compartments is manifold. Possible areas of applications are longitudinal tumor volumetry (Alberts et al., 2016b;Meier et al., 2016a), planning and assessment of neurosurgical interventions (Meier et al., 2016b;Porz et al., 2016) and radiomic analyses (Cui et al., 2015;Gutman et al., 2015;Rios Velazquez et al., 2015). We strongly believe that an estimate of segmentation uncertainty will help in any clinical decision-making process that is based on information from image segmentations. Uncertainty estimation could be used for various purposes such as the exclusion of uncertain labeled voxels from further analysis, or the guidance of an expert rater in the correction of automatically generated segmentations. The former situation would be of great importance for the emerging field of radiomics in which imaging biomarkers are discovered in a high-throughput setting based on features extracted from image segmentations (Gillies et al., 2016;Yip and Aerts, 2016).

Methodological aspects
The proposed Perturb-and-MPM methodology enables approximate sampling from dense multi-label CRFs. We showed that the approximation error of Perturb-and-MPM is bounded and converges with an increasing number of samples to a residual error. For tractable problems, this residual error was shown to be consistently lower than the error of a mean field approximation. This empirically confirms our intuition in Section 3.2 ( Figure 2). Furthermore, our theoretical analysis suggested that the residual error is driven by the error occurring in the situationX M AP =X M P M since forX M AP =X M P M the approximation error vanishes with an increasing number of samples. The segmentation uncertainty is quantified based on the Shannon entropy contained in the voxel-wise marginal distributions estimated by Perturb-and-MPM. As a consequence, the error in uncertainty quantification can also be bounded. Moreover, we emphasize that the estimated uncertainty is directly linked to the probabilistic model used for generating the segmentation itself. The computational overhead of the method is modest (sampling one segmentation corresponds to one inference step). Perturb-and-MPM can essentially be applied to any dense CRF that relies on a mean field approximation, e.g. (Kamnitsas et al., 2016).

Comparison to other methods
Methods used for sampling image segmentations can be grouped in two categories: Methods that perform sampling over the complete image grid, and methods that perform sampling from a parametric representation of the region of interest. The first group of methods relies conventionally on MCMC-based sampling approaches. They exhibit several disadvantages such as burn-in phase and typically bad scalability, which leads to a slow computation time (Iglesias et al., 2013). The computation time can be improved by employing parallelization techniques (Byrd et al., 2010). Additionally, MCMC-based sampling approaches usually rely on a number of sampling hyperparameters (e.g. step sizes, etc.) which need to be fine-tuned. In contrast, our method does neither require a burn-in phase nor does it require the tuning of any hyperparameters. The second group of methods typically employs curves as parametric representation for segmentation boundaries (Fan et al., 2007;Lê et al., 2016). The parametric representation enables a more efficient sampling than in conventional MCMC-based approaches. However, these approaches do not handle sampling of multi-label segmentations naturally. In contrast, our method enables efficient sampling of multi-label segmentations over the complete image grid (corresponding to the CRF). Kohli and Torr (2008) proposed a method based on min-marginals for the computation of labeling uncertainty in pairwise, grid-structured CRFs. Parisot et al. (2014) used this technique within a combined segmentation-registration framework in order to segment binary tumor masks of low-grade glioma. Alberts et al. (2016a) used it in context of longitudinal tumor volumetry of gliomas. The estimation of min-marginals relies on the computation of st-cuts, which in case of multi-label problems is only possible for a very limited class of energy functions (Ishikawa, 2003). More importantly, their estimation is not possible for non-convex priors such as e.g. the Potts prior commonly used in image segmentation.

Experimental aspects
We investigated the quality and potential radiological utility of the uncertainty estimation for the multi-label segmentation of the tumor core in glioblastoma patients. We chose this particular segmentation task since recent benchmarks (Menze et al., 2015) highlighted that the segmentation of the tumor core appears to be the most challenging segmentation task in glioblastoma for current state-of-the-art algorithms. In addition, the individual compartments of the tumor core have shown to be associated with molecular characteristics Naeini et al., 2013) and patient survival (Pope et al., 2005;Rao et al., 2016;Rios Velazquez et al., 2015), emphasizing the importance of their accurate definition via automatic segmentation. Quality was defined in terms of sensitivity and specificity of the uncertainty estimation with regards to segmentation errors. We found that uncertain regions correspond well to wrongly segmented regions. However, we observed that uncertain voxels are mostly situated at the border/periphery of such a region, whereas central areas show less or no uncertainty (cf. Figure 4 and 8). Intuitively, it does make sense that an algorithm is most uncertain in delineating the interface of different tissues. This behavior can also be seen in other perturbation-based approaches (Alberts et al., 2016b;Papandreou and Yuille, 2014). In general, the apparent uncertainty is dependent on the marginal distributions of the model. Thus, one can observe that the spatial extent of segmentation uncertainty does not change anymore after a number of samples (cf. Figure 4). If the model at hand has very peaky distributions (i.e. is overconfident), the corresponding uncertainty map will be less informative. As a consequence, wrongly labeled image regions may not be reflected in the corresponding uncertainty map. This observation brings interesting aspects for future developments, where the interplay between Perturb-and-MPM, to assess uncertainty, and the design as well as training of the underlying segmentation model, is to be considered.

Clinical aspects
The clinical utility of Perturb-and-MPM was demonstrated for the estimation of the EOR and RTV in glioblastoma patients. The EOR and RTV are important prognostic factors and hence their estimation is essential for patient prognosis (Chaichana et al., 2014;Grabowski et al., 2014;Stummer et al., 2008). We assessed the ability of Perturb-and-MPM to indicate uncertain voxels to be excluded from the computation of EOR and RTV, respectively. This resulted in a consistent improvement in the estimation of both parameters with respect to ground truth data of four expert raters. Consequently, we would expect also an improvement in association of such automatically extracted and corrected imaging biomarkers with patient survival. Similarly, we expect that these methods can contribute for a better assessment of disease progression and response to therapy, where a precise longitudinal assessment of tumor volumes is very important for clinical decision making Meier et al., 2016a).

Conclusion
Perturb-and-MPM is an easy-to-implement approach for performing approximate sampling from dense multi-label CRFs. It allows the generation of spatially-resolved uncertainty maps for image segmentations. This opens up opportunities to create time-effective human-machine interfaces necessary to monitor and correct results of automated segmentation. More importantly, it can have a substantial impact in clinical scenarios such as radiotherapy and neurosurgery, where an accurate delineation of the tumor is needed in order to precisely target the tumor while preserving as much healthy tissues as possible. In turn, uncertainty estimates of the segmentation can be used as data quality measures in high-throughput radiomics studies relying on these segmentation results to extract imaging biomarkers of diagnosis, assessment of prognosis, and prediction of therapy response.

Funding
This project has received funding from the European Unions Seventh Framework Programme for research, technological development and demonstration under grant agreement N o 600841 and the Swiss National Science Foundation project 205321 169607.