Benoît Audelan
email: benoit.audelan@inria.fr

Hervé Delingette

Unsupervised quality control of segmentations based on a smoothness and intensity probabilistic model

Keywords: Unsupervised quality control, image segmentation, Bayesian learning, spatial regularization

Monitoring the quality of image segmentation is key to many clinical applications. This quality assessment can be carried out by a human expert when the number of cases is limited. However, it becomes onerous when dealing with large image databases, so partial automation of this process is preferable.

works have proposed both supervised and unsupervised methods for the automated control of image segmentations. The former assume the availability of a subset of trusted segmented images on which supervised learning is performed, while the latter does not. In this paper, we introduce a novel unsupervised approach for quality assessment of segmented images based on a generic probabilistic model. Quality estimates are produced by comparing each segmentation with the output of a probabilistic segmentation model that relies on intensity and smoothness assumptions. Ranking cases with respect to these two assumptions allows the most challenging cases in a dataset to be detected. Furthermore, unlike prior work, our approach enables possible segmentation errors to be localized within an image. The proposed generic probabilistic segmentation method combines intensity mixture distributions with spatial regularization prior models whose parameters are estimated with variational Bayesian techniques. We introduce a novel smoothness prior based on the penalization of the derivatives of label maps which allows an automatic estimation of its hyperparameter in a fully data-driven way. Extensive evaluation of quality control on medi-

Introduction

Semantic segmentation of an image is the process of associating a label to every pixel in an image. This task is particularly important in a medical context since it impacts downstream algorithms using image segmentations as input, but also the decisions that clinicians may make about the patient. For instance, in radiotherapy planning, the delineations of tumor lesions directly influence the extent of the dose delivered around the tumor. Also, obtaining reliable image segmentations is mandatory to use image derived biomarkers in a clinical setting [START_REF] Keshavan | Mindcontrol: A web application for brain segmentation quality control[END_REF]. Finally, the development of supervised learning for image segmentation requires the accumulation of potentially large sets of manually or semi-manually segmented image databases that need to be quality controlled. Such segmentations are prone to inter-rater variability [START_REF] Visser | Inter-rater agreement in glioma segmentations on longitudinal MRI[END_REF] in addition to plain errors. It is therefore of great importance to automatically detect possible failed segmentation cases, whether those segmentations are generated by an algorithm or a human rater. The challenge is to perform this monitoring in the absence of ground truth segmentations.

In prior work, evaluation methods can be categorized either as supervised or as unsupervised, depending on whether a reference segmentation is required or not [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF]. A first set of supervised methods is based on a classifier which accepts or rejects the proposed segmentation based on combined features.

For instance, in Hui [START_REF] Zhang | Meta-Evaluation of Image Segmentation Using Machine Learning[END_REF], decision trees based on handcrafted features depending on the image (texture, color space...) and on the geometry of the segmented region (perimeter, compactness...) are combined in a single classifier. In [START_REF] Shamir | Evaluation of head segmentation quality for treatment planning of tumor treating fields in brain tumors[END_REF], a decision tree predicts the Dice score of head segmentations with an application to the treatment of brain tumors. In [START_REF] Xu | Automated Quality Control for Segmentation of Myocardial Perfusion SPECT[END_REF], a framework to detect failures in cardiac segmentation based on a shape parameter and an intensity feature has been proposed. The number of features taken into account is increased in [START_REF] Kohlberger | Evaluating Segmentation Error without Ground Truth[END_REF], where the model decision relies on 42 shape and appearance features. They are combined in an SVM classifier regressing the Dice coefficient between the given segmentation and the unknown ground truth. While in [START_REF] Xu | Automated Quality Control for Segmentation of Myocardial Perfusion SPECT[END_REF] the features were specific to cardiac segmentation, the approach taken in [START_REF] Kohlberger | Evaluating Segmentation Error without Ground Truth[END_REF] is more generic and was trained on segmentations of 8 different organs.

Reverse Classification Accuracy (RCA) has also been proposed for quality control assessment in [START_REF] Valindria | Reverse Classification Accuracy: Predicting Segmentation Performance in the Absence of Ground Truth[END_REF]. Assuming the availability of a set of trusted images with ground truth, the proposed segmentation on a new image is compared to the predicted one based on those reference images, which can result in rejection if discrepancies are too large. This approach was tested on larger databases in [START_REF] Robinson | Automated quality control in image segmentation: application to the UK Biobank cardiovascular magnetic resonance imaging study[END_REF] where the authors showed the ability of the method to highlight poor quality segmentations but pointed out the relatively long computation times as a bottleneck.

Another family of supervised approaches uses deep learning to estimate the quality of a segmentation. For instance, in [START_REF] Robinson | Real-Time Prediction of Segmentation Quality[END_REF], a neural network is trained to predict the Dice coefficient of cardiac segmentations. The Jaccard index (intersection over union) is predicted by neural networks in [START_REF] Arbelle | QANet -Quality Assurance Network for Image Segmentation[END_REF][START_REF] Huang | QualityNet: Segmentation quality evaluation with deep convolutional networks[END_REF][START_REF] Shi | Segmentation quality evaluation based on multi-scale convolutional neural networks[END_REF] where the original image and the proposed segmentation mask are provided as input. Some authors have proposed exploiting the uncertainty of segmentations in order to assess their quality, within a deep learning framework. Uncertainty quantification also adds some interpretability to the quality assessment as it provides information about the location of possible errors. Bayesian QuickNat proposed by [START_REF] Roy | Bayesian QuickNAT: Model uncertainty in deep whole-brain segmentation for structure-wise quality control[END_REF] uses Monte Carlo dropout at test time to generate several segmentation samples.

The average over the samples gives the final segmentation map while variability across the different samples gives an estimate of the uncertainty of the segmentation. The authors show a good correlation between the measured uncertainty and the Dice coefficient between the segmentation and the unknown ground truth. Other methods to evaluate the uncertainty were explored in [START_REF] Jungo | Assessing Reliability and Challenges of Uncertainty Estimations for Medical Image Segmentation[END_REF] and the results suggest that none is superior to the others. Finally in DeVries and Taylor (2018), a first network outputs a segmentation map and an uncertainty map at the pixel level, which are then taken as inputs by a second network which regresses a quality score at the image level.

A limitation shared by these methods is their supervised design, meaning that they require the extraction of a subset of segmented data that is considered to be "ground truth". This trusted subset is used by the models to learn how a "good" segmentation looks. The resulting decision rules making a new segmentation acceptable or not may thus be biased by the composition of the trusted set, which must be large enough for training a deep-learning-based framework.

Further, access to large annotated datasets remains an issue in many domains including medical imaging. Finally, supervised methods often lack generality as their performance depends on the type of images and segmented structures in the training set.

In contrast, unsupervised approaches do not rely on a subset of trusted images but rather on assumptions about the appearance and shape of the foreground and background regions [START_REF] Rosenberger | Unsupervised and supervised image segmentation evaluation[END_REF][START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF].

These assumptions are then translated into a set of segmentation metrics. For instance, common hypothesis is that a "good" segmentation exhibits high levels of intra-region homogeneity and inter-region heterogeneity [START_REF] Johnson | Unsupervised image segmentation evaluation and refinement using a multi-scale approach[END_REF], and several handcrafted features have been proposed to measure them [START_REF] Chabrier | Unsupervised Performance Evaluation of Image Segmentation[END_REF][START_REF] Gao | A Novel Unsupervised Segmentation Quality Evaluation Method for Remote Sensing Images[END_REF][START_REF] Johnson | Unsupervised image segmentation evaluation and refinement using a multi-scale approach[END_REF][START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF]. The main limitation of these approaches is that it is difficult to design discriminative indices and to find a proper way to combine them. Moreover, as mentioned by [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF], most of those metrics assume a single underlying intensity distribution, typically Gaussian, in both foreground and background regions which is overly simplistic and sensitive to outliers.

Last but not least, interpretability is a desirable property, as knowing the problematic regions could facilitate the segmentation curation. However, it is often an issue since many of the previous methods, supervised or not, are black boxes outputting a simple score, which does not help to understand why a segmentation has failed.

In this paper, we propose a novel unsupervised approach for automated quality assessment of image segmentations. It is based on the comparison between a proposed segmentation S produced by an algorithm or a human rater and the segmentation M given by a generic probabilistic segmentation model. The generic model is based on two simple intensity and smoothness assumptions, the underlying hypothesis being that explainable segmentations correspond to clearly visible boundaries in the image well captured by M . On the contrary, segmentations far from M are categorized as difficult or challenging as they would require priors other than intensity and smoothness to be explained. The quality assessment of a set of segmented images is then performed by studying how the distance between the proposed segmentation S and the modelled segmentation M varies within the dataset. Segmentations that are lying on the tails of this distance distribution are considered to be atypical and are candidates for manual verification. We show the effectiveness of this approach to extract suspicious segmentations on various public datasets ranging from photographic images for object detection and segmentation (COCO dataset) to lung and brain medical images (LIDC and BRATS datasets). We also show that this approach can be used in some cases to predict the performance of segmentation algorithms.

Our main contributions are twofold:

• Instead of relying on an arbitrary subset of selected segmentations as a training set, we propose an unsupervised approach based on intensity and smoothness hypotheses without any prior knowledge of the structure to be segmented. It removes the bias related to the selection of the reference images and allows the quality of segmentations to be assessed when few or even no other segmentations are available from a database. Our method differs from previous unsupervised segmentation quality indices with a more complex and robust approach to modeling the intensity of the different regions in the image. In addition, it allows a combination of the key factors defining a "good" segmentation (i.e., the intra-region homogeneity and the inter-region heterogeneity) in a data-driven way. Last but not least, our method is visually interpretable. For instance, when dealing with 3D medical images, it allows automatic retrieval of the slices with suspicious segmentations. Finally, the result can be useful to guide the manual correction of poorly segmented cases.

• We provide different spatial regularization strategies to enforce the spatial continuity. In particular, we introduce a novel prior, denoted by FDSP (Finite Difference Spatial Prior), based on the penalization of the squared norm of the derivatives of the prior label map, which allows an adaptative learning of the hyperparameter. It is compared to the classical Markov random field (MRF) and another spatial prior based on a weighted combination of spatially smooth kernels introduced in an earlier work of the authors [START_REF] Audelan | Unsupervised Quality Control of Image Segmentation Based on Bayesian Learning[END_REF], which will be denoted by GLSP (Generalized Linear Spatial Prior) throughout the paper. This paper expands [START_REF] Audelan | Unsupervised Quality Control of Image Segmentation Based on Bayesian Learning[END_REF] by proposing a different spatial regularization strategy for which the hyperparameter can be estimated.

In addition, the novel regularization prior can be entirely inferred with a variational Bayes method (no Laplace approximation needed) and leads to much faster computations. We also provide more extensive experiments on different datasets and added a qualitative and quantitative comparison with unsupervised segmentation quality control indices proposed in prior works. The code with the different regularization strategies is available in this repository: https://gitlab.inria.fr/epione/unsegqc.

The rest of the paper is organized as follows. Section 2 presents the general framework of our unsupervised quality control assessment. In section 3, we present our appearance model and the spatial priors. Section 4 describes the model inference depending on the regularization. Finally, we show in section 5 the relevance of our approach for segmentation quality control on several datasets. Supervised segmentation quality control methods require the existence of a trusted subset of data from which quality assessment is learned. Instead, we follow an unsupervised approach (see Fig. 1) based on a probabilistic segmentation model relying only on two simple smoothness and intensity assumptions.

Unsupervised quality control workflow

Its great advantage is that it is agnostic with respect to the structure to be segmented and therefore can be run automatically even in the absence of ground truth.

Input segmentation

The input of the proposed method is a binary segmentation S on an image I into foreground and background regions for which we would like a quality estimate. There are no restrictions regarding the origin of S as it can have been created by an algorithm or a human rater. Note that this is in contrast to several other methods that require the input segmentation to have been generated by a specific algorithm, like the uncertainty-based methods in deep learning [START_REF] Devries | Leveraging Uncertainty Estimates for Predicting Segmentation Quality[END_REF][START_REF] Jungo | Assessing Reliability and Challenges of Uncertainty Estimations for Medical Image Segmentation[END_REF][START_REF] Roy | Bayesian QuickNAT: Model uncertainty in deep whole-brain segmentation for structure-wise quality control[END_REF].

Probabilistic model

Given the segmentation S, we produce a smooth contour or surface M close to S which is mostly aligned with visible contours in the image. We stress that the objective is not to build a surrogate ground truth, but instead to use M only as a comparison tool.

Intensity assumption. The first hypothesis of our approach is that intensity distribution variations in the image can help to understand segmentations.

Given the segmentation S, two intensity models are built for the foreground and background regions.

Spatial smoothness assumption. The second hypothesis relies on the generally accepted assumption that two neighbouring voxels share a higher probability of belonging to the same label region. This is classically enforced by the use of discrete priors such as MRF. In [START_REF] Audelan | Unsupervised Quality Control of Image Segmentation Based on Bayesian Learning[END_REF], we proposed a regularization strategy based on a combination of spatially smooth kernels (GLSP). In addition to these two possibilities, we introduce in this paper a novel way to take into account the spatial organization of the voxels, which we call the Finite Difference Spatial Prior (FDSP). This approach allows full tractability of the hyperparameter in an efficient manner which is not possible for the MRF and GLSP formulations.

The two assumptions are combined into a probabilistic model that outputs a new segmentation map M . By construction, M is typically a smooth contour which is mostly aligned with the visible intensity boundaries in the image.

Again, M should only be seen as a representation used to benchmark the input segmentation S.

To measure the adequacy of S with respect to M , we employ the average asymmetric surface error (ASE) defined as

E S = d(S, M ) = 1 ∂S x∈∂S min y∈∂M d(x, y)
where ∂ denotes the segmentation surface. We discard the metric d(M, S) as being uninformative since M is not a surrogate ground truth. An alternative measure to the ASE used in this paper is the Dice score computed between the segmentations M and S.

Detection of challenging cases

Segmentations S close to M are identified as being explained by the model.

In that case, the two intensity and spatial smoothness assumptions upon which the probabilistic model is based are sufficient to understand the contours. However, segmentations S far from M are classified as unexplained or challenging.

Typically, contours crossing large regions of uniform intensity distribution would be identified as unexplained by our model. It is important to note that having an unexplained segmentation does not imply that this segmentation is wrong.

It simply means that other priors besides those of smoothness and intensity are required to understand its boundaries.

Use case

We believe our approach is particularly interesting when dealing with a whole set of segmentations. For instance, say we are given a set of images with corresponding annotations. The comparison of adequacies between S and M for all images allows the detection of atypical cases which behave differently from the majority of the distribution, and for which a visual inspection might be worthy. On the contrary, applying the method on a single image is not the ideal use case as the analysis of the result is difficult without any comparison with similar images.

Our approach is unsupervised, generic, and based on few simple assumptions.

However this comes with intrinsic limitations. For instance, any irrelevant contour following visible intensity boundaries will not be considered as suspicious.

This limitation is common to all previously proposed unsupervised methods.

More generally, the proposed method is not intended to return all erroneous segmentations inside a dataset (which is expected from a supervised approach) but instead to extract some suspicious cases when limited information is available.

Method

In this section, we review the details of our probabilistic model. We consider a binary image segmentation problem for isolating a single structure from an image I made of N voxels in a grid of dimension D (D = 2, 3) having intensity

I n ∈ R v , n = 1, . . . , N , where v ≥ 1 (v = 1,
3 and 4 in practice). We introduce for each voxel a binary hidden random variable Z n ∈ {0, 1} with Z n = 1 if voxel n belongs to the structure of interest.

Mixtures of multivariate Student's t-distributions

Appearance models of the foreground and background regions of S are defined respectively by the two image likelihoods p(I n |Z n = 1, θ 1 I ) and p(I n |Z n = 0, θ 0 I ) where θ 0 I , θ 1 I are parameters governing those models. In this paper, we consider generic parametric appearance models as variational mixtures of multivariate Student's t-distributions [START_REF] Archambeau | Robust Bayesian clustering[END_REF]. The Student's t generalizes the Gaussian distribution with heavy tails and leads to robust mean and covariance estimates. The number of components in the mixture is automatically estimated by using a sparsity-inducing Dirichlet prior over the mixture proportions which automatically prunes the components with a small number of samples. Finally, we introduce the appearance probability ratio r n defined as:

r n (I, θ 0 I , θ 1 I ) p(I n |Z n = 1, θ 1 I ) p(I n |Z n = 0, θ 0 I ) + p(I n |Z n = 1, θ 1 I ) , (1)
which is the posterior label probability with a non-informative prior (p(Z n = 1) = 0.5).

Spatial smoothness prior

The spatial smoothness prior allows the spatial organization between voxels to be taken into account and a certain degree of continuity to be enforced.

To this end, different strategies can be employed. In this paper, we propose to compare one discrete prior (MRF) with two continuous priors (GLSP and FDSP), the third one being novel to the best of our knowledge.

MRF prior

The classical MRF formulation relies on labels of neighbouring voxels. In a binary segmentation problem, a natural way to enforce spatial smoothness is the Ising model. Assuming β to be the hyperparameter of the MRF, the label prior probability is given by:

p(Z|β) = 1 T (β) exp    β 2 N i=1 j∈δi Z i Z j    , (2)
where δ i are the neighbouring voxels of i and T (β) is the partition function. In practice, we consider 4-and 6-connectivity neighborhoods for 2D and 3D images, respectively. The value of β represents the strength of association between neighbouring voxels: β = 0 corresponds to a model with no spatial prior, while large positive values encourage neighbouring voxels to have the same label. The Ising model may be replaced by an image contrast sensitive prior as performed for instance in the GrabCut algorithm [START_REF] Rother | GrabCut": Interactive Foreground Extraction Using Iterated Graph Cuts[END_REF].

The computation of the partition function T (β), needed for an automatic estimation of the model's hyperparameter β, requires considering all possible configurations of the MRF which is not computationally tractable for large lattices. Therefore, β has to be fixed by the user.

Generalized Linear Spatial Prior

In Audelan and Delingette (2019), we proposed a continuous label prior denoted by Generalized Linear Spatial Prior (GLSP) to enforce the spatial continuity. The prior is defined through a generalized linear model of spatially smooth functions. More precisely, the prior probability p(Z n = 1) is defined as a Bernoulli distribution whose parameter is a spatially random function specified as a generalized linear model:

p(Z n = 1|W ) = σ L l=1 Φ l (x n )w l , (3)
where x n ∈ R D is the voxel position in an image of dimension D and the link function σ(u) is the sigmoid function σ(f ) = 1/(1 + exp(-f )). The basis {Φ l (x)} are L functions of space, typically radial basis functions (for instance, Gaussian functions) defined on a regular grid, and w l ∈ W are weights considered as random variables. Thus the prior probabilities of two geometrically close voxels are related to each other through the smoothness of the function

f (x n ) = L l=1 Φ l (x n )w l = Φ T n W , writing Φ T n = [Φ 1 (x n ), • • • , Φ L (x n )].
The smoothness of the label prior σ (f (x n )) depends on the choice of the L basis functions {Φ l (x)} which are commonly uniformly spread over the image domain. The key parameters are the spacing between the basis centers, the standard deviations (or radii) r of the Gaussian functions and the position of the origin basis. Together, they influence the amount of smoothing brought by the label prior, large spacing and standard deviations leading to smoother prior probability maps.

To obtain a robust description, the weight vector W = [w 1 , . . . , w L ] T is fitted with a zero mean Gaussian prior parameterized by the diagonal precision matrix αI L : p(W ) = N (0, α -1 I L ). Finally, a non-informative prior is chosen for α, p(α) ∝ 1. In contrast to the MRF formulation, a Bayesian inference of the hyperparameter is possible here, as shown in section 4.2.

Finite Difference Spatial Prior

As a third regularization strategy, we introduce in this paper the Finite Difference Spatial Prior (FDSP). The prior probability p(Z n = 1) is again defined as a Bernoulli distribution whose parameter belongs to a spatially smooth random field:

p(Z n = 1|W ) = σ (w n ) , (4)
where σ(u) is once more the sigmoid function. The smoothness of the label field is caused by a prior applied to the vector W = [w 1 , . . . , w n ] T penalizing the squared norm of its derivatives of order p:

p(W |α) = 1 T (α) exp -α N n=1 ∆ p (w n ) 2 , (5)
where ∆ p (w n ) is the p order central finite difference operator at w n and T (α)

is the normalization factor. The quantity ∆ p (w n ) is a tensor of order p approximating the p-order derivatives of the scalar field defined by w n . Since

the function h(x) = ∆ p (x)
2 is 2-homogenous, we know that the normalization factor has the form T (α) = cα -N/2 where c is constant independent of α [START_REF] Pereyra | Maximum-aposteriori estimation with unknown regularisation parameters[END_REF]. One can easily show that p(W |α) is a zero mean Gaussian distribution whose precision matrix consists of difference operators. The value of the parameter α controls the amount of the spatial regularization applied to the weights W .

In this paper, we consider only first order derivatives (p = 1) corresponding to the discretization of the Dirichlet energy. In that case Eq. 5 is written:

p(W |α) = cα N 2 exp - α 4 N n=1 D d=1 (w δ d (n+1) -w δ d (n-1) ) 2 , ( 6
)
where δ d (n + i) represents the neighbor of index i of voxel n in the dimension d.

The graphical models of the different segmentation frameworks are shown in Figure 2.

Implementation

The second step is to compute a variational approximation of the posterior

P (Z n |I n ) which
I n Z n β θ k I N (a) I n Z n w l α θ k I x n N L (b) I n Z n w n α θ k I N (c)

Probabilistic inference

MRF regularization

A classical way to maximize the log likelihood log p(I) with an MRF prior is to use variational inference with a mean field approximation [START_REF] Ambroise | Convergence of an EM-type algorithm for spatial clustering[END_REF][START_REF] Roche | On the convergence of EM-like algorithms for image segmentation using Markov random fields[END_REF]. The label posterior distribution q(Z) is assumed to factorize as i q i (Z i ), which leads to the following fixed-point equation for voxel i at iteration m + 1:

q m+1 ip = r i exp{β j∈δi q m jp } 1 k=0 r k i (1 -r i ) 1-k exp{β j∈δi q m jk } , (7)
where p ∈ {0, 1}, q ik represents q i (Z i = k), r i is the appearance probability ratio for voxel i and β is fixed by the user.

GLSP regularization

A type-II maximum likelihood approach is used to estimate the model parameters. A Gaussian approximation for the weights posterior distribution is found by computing a Laplace approximation through iterative reweighted least squares. The parameter α is then updated by maximizing the marginal likelihood. We refer to the original paper for more details [START_REF] Audelan | Unsupervised Quality Control of Image Segmentation Based on Bayesian Learning[END_REF].

The algorithm requires several steps of covariance matrix inversion which can be prohibitive for large images. The problem was addressed in Audelan and Delingette (2019) by splitting the narrow band into smaller overlapping patches that were then merged. In this paper, the code was improved and the decomposition into patches is no longer needed.

FDSP regularization

We propose a variational inference scheme to estimate prior and hyperprior parameters U = {Z, W, α}. Variational inference approximates the true posterior p(U |I) by a chosen family of distributions q(U ). Maximizing the data log likelihood log p(I) implies minimizing the Kullbach-Leibler divergence between q(U ) and p(U |I) or equivalently maximizing the lower bound L(q):

log p(I) = U q(U ) log p(I, U ) q(U ) dU L(q) + KL [q(U )||p(U |I)] . (8)
We assume that the approximation of the posterior can be factorized as q(U ) = q Z (Z)q W (W )q α (α). The lower bound can thus be re-written as:

log p(I) L(q) = Z α W q Z (Z)q W (W )q α (α) log p(I|Z)p(Z|W )p(W |α) q Z (Z)q W (W )q α (α) dW dα . (9)
We can further expand the factors defining the joint probability:

p(I|Z) = n r Zn n (1 -r n ) 1-Zn .
The spatial prior p(Z n |W ) can be likewise written as [σ(w n )] Zn [σ(-w n )] 1-Zn and the weights prior p(W |α) is given by ( 6) for first order derivatives.

However, the right hand side of ( 9) is intractable because the spatial prior does not belong to the exponential family (due to the sigmoid function). As an alternative to the Laplace approximation, we use a local variational bound as introduced in [START_REF] Jaakkola | Bayesian parameter estimation via variational methods[END_REF] in the context of logistic regression. In this case, we replace the sigmoid function with a well-chosen lower bound:

σ(x) g(x, ξ) = σ(ξ) exp (x -ξ)/2 -λ(ξ)(x 2 -ξ 2 )
. ξ is a variational parameter and λ(ξ) = tanh(ξ/2)/(4ξ). The spatial prior p(Z|W ) can thus be approximated by

F (Z, W, ξ) = n [g(w n , ξ n )] Zn [g(-w n , ξ n )] 1-Zn . This
approximation leads to a new lower bound J (q) on the lower bound L(q):

log p(I) L(q) J (q) = Z α W q Z (Z)q W (W )q α (α) log p(I|Z)F (Z, W, ξ)p(W |α) q Z (Z)q W (W )q α (α) dW dα . (10)
This new lower bound J (q) is now tractable and the optima q * for each of the variational posteriors can be derived by variational calculus (See Appendix B for details of the derivations). q * Z (Z) is therefore given by q * Z (Z) = n η Zn n1 η 1-Zn n0 with η nk = ρ nk / k ρ nk for k ∈ {0, 1} and:

ρ nk = r k n (1 -r n ) 1-k σ(ξ n ) exp (-1) 1-k E[w n ] 2 - ξ n 2 -λ(ξ n )(E[w 2 n ] -ξ 2 n ) . (11)
By further assuming that q W (W ) = n q wn (w n ), the variational optimization for q wn (w n ) yields a normal distribution of the form q * wn (w n ) = N (µ wn , Σ wn ). A fixed-point equation is found for updating the mean. For first order derivatives, we have:

Σ wn = 2λ(ξ n ) + 2 d α 2 -1 , (12)
µ wn = Σ wn η n1 - 1 2 + α 2 d µ w δ d (n+2) + µ w δ d (n-2) . (13)
The variational posterior q α (α) is assumed to be a Dirac distribution which leads to the following update:

α -1 = 1 2N n D d=1 E (w δ d (n+1) -w δ d (n-1) ) 2 . ( 14
)
Finally, following [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF], maximizing (10) with respect to ξ n gives an update formula of the form:

ξ 2 n = E[w 2 n ] . (15)
To compute ( 11), ( 14) and ( 15), we need the expectations

E[w n ], E (w δ d (n+1) -w δ d (n-1) ) 2
and E[w 2 n ] with respect to the variational distribution q wn . They can be easily evaluated to give

E[w n ] = µ wn , E[w 2 n ] = Σ wn +µ 2 wn and E (w δ d (n+1) -w δ d (n-1) ) 2 = µ 2 w δ d (n+1) + µ 2 w δ d (n-1) -2µ w δ d (n+1) µ w δ d (n-1) + Σ w δ d (n+1) + Σ w δ d (n-1
) . After convergence, the variational distribution q Z (Z) gives an approxima-tion to the posterior label probability p(Z n = 1|I, W ), which combines prior and intensity likelihoods. Finally, the maximum a posteriori estimate of the segmented structure is obtained as the isosurface p(Z n = 1|I, W ) = 0.5. This approach has some advantages in comparison with the first two. First, it allows an automatic estimation of all its parameters. For the MRF, the user needs to fix β and for the GLSP, the layout of the basis functions and their radii are also user-defined. Moreover, a lower bound (Fig. 3e) on the marginal likelihood can be computed in this case, which can be used to monitor the convergence and is helpful to compare segmentation results. The computation of the lower bound is given in Appendix C.

Results

Datasets

The proposed method was evaluated on four publicly available datasets: the BRATS 2017 training and validation datasets [START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BRATS)[END_REF], the LIDC dataset [START_REF] Armato | The lung image database consortium (LIDC) and image database resource initiative (IDRI): A completed reference database of lung nodules on CT scans[END_REF], the training data from the MSSEG challenge [START_REF] Commowick | Objective Evaluation of Multiple Sclerosis Lesion Segmentation using a Data Management and Processing Infrastructure[END_REF] and finally the COCO 2017 validation dataset [START_REF] Lin | Microsoft COCO: Common Objects in Context[END_REF]. The LIDC dataset comprises 1018 pulmonary CT scans with 0.6 mm to 5.0 mm slice thickness. The in-plane pixel size ranges from 0.461 mm to 0.977 mm. Each scan was reviewed by 4 radiologists who annotated lesions of sizes ranging from 3 mm to 30 mm. Annotations include localization and manual delineations of the nodules. Up to 4 segmentations can be available for the same nodule, depending on the number of radiologists who considered the lesion to be a nodule. In this paper, all scans were first re-sampled to 1 mm 3 resolution as pre-processing step, and we restrict the analysis to nodules of diameter above 20 mm, i.e. 309 segmentations.

The MSSEG training dataset contains MR data from 15 multiple sclerosis (MS) patients. Manual delineations of lesions were performed on the FLAIR sequence by seven experts.

Finally, COCO is a large-scale object detection and segmentation dataset of real world images. The 2017 validation set contains 5000 images with 80 object categories, ground truth object classification, object localization and segmentation. To annotate such a large number of images, the authors resorted to a crowd-sourcing annotation pipeline.

Unsupervised indices

As discussed in section 1, different indices have been proposed in prior works for unsupervised segmentation evaluation. We selected 4 of them in order to provide a qualitative and quantitative comparison with our approach. They all involve the computation of 2 metrics, the former measuring the intra-region uniformity while the latter gives an estimate of the inter-region disparity.

Three out of the four indices are taken from [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF]: Zeb, η and F RC . The last one was introduced in Johnson and Xie (2011) and is denoted by GS in this paper. Formula are given in Appendix A.

Setting hyperparameters

Width of the narrow band

As noted in section 3.3, the analysis is restricted to a narrow band alongside the input segmentation's contour. The width of this narrow band controls the extent of the region taken into account for learning the appearance models for both background and foreground and fitting the regularization model.

We assessed the sensitivity of the results to this hyperparameter on the BRATS training set and LIDC dataset. We applied the algorithm for several narrow band widths using FDSP as a spatial prior. Different ASE values were obtained for each segmentation depending on the narrow band setting. We then While the sets from the BRATS training set are rather stable, those from the LIDC dataset show some variability. An explanation of the sensitivity of LIDC segmentations to the narrow band width can be found in Fig. 5. If the narrow band is too wide, the high intensity differences between the pleura and the lung parenchyma lead appearance models of nodules close to the pleura to leak outside the lung. In brief, the sensitivity of the algorithm with respect to the narrow band's width varies from case to case. As the computation time is not a bottleneck for FDSP regularization, we propose in practice to perform the analysis with different width settings and then choose the one leading to the most stable and reasonable results.

We would also like to underline that previously published unsupervised indices are likewise sensitive to the width of the narrow band. An example is shown in Fig. 4 for the indicator GS. In order to provide a fair comparison between approaches, the computations of the selected unsupervised indices are always performed on the same narrow band as the one used for our method.

Other hyperparameters

Among the parameters that need to be defined by the user is the number of components for the mixtures of multivariate Student's t-distributions. It is fixed to 7 in all our experiments. This parameter is not so sensitive as unnecessary components will be pruned by the Dirichlet prior and removed from the model.

The number of remaining parameters depends on the chosen spatial prior.

For an MRF prior, the user needs to provide a value for β, which controls the strength of the regularization. We tested 3 values for this hyperparameter throughout our experiments: 0.2, 1 and 3. For a GLSP regularization, the user has to define a dictionary of basis functions whose key parameters are the step between each basis function and their radii. They likewise control the amount of regularization. In this paper, we set the step to 6 vx and the radius to 17 vx, except for the LIDC dataset for which the step was set to 4 vx and the radius to 12 vx. Finally, for the regularization using an FDSP prior, no further parameter needs to be set by the user as the model's hyperparameters are all learnt automatically, which is a great advantage in comparison with the first two approaches. In the case of segmentations produced by human raters, possibly with the help of interactive annotation tools, it is very useful to be able to rank segmentations, highlight potentially difficult segmentations and track possible errors in large databases.

Qualitative analysis

In this section we present some results from two datasets of medical images, Furthermore, one can see that contours in the right tail samples from BRATS are more irregular and that intensity variations in some regions are very weak making their accuracy questionable. Those contours were probably extracted through thresholding instead of being manually drawn as was permitted in the annotation process [START_REF] Jakab | Segmenting Brain Tumors with the Slicer 3D Software[END_REF][START_REF] Menze | The Multimodal Brain Tumor Image Segmentation Benchmark (BRATS)[END_REF]. Similarly, some contours in the right tail samples from LIDC cross regions of uniform intensity and therefore require other priors like shape to be explained. Yet, the contours are far from obvious in some areas in comparison with the left tail samples. Therefore, our approach fulfills its role of extracting challenging, possibly suspicious, cases within a dataset.

We present in Fig. 9 a qualitative comparison between the spatial priors proposed for our approach and the unsupervised indices presented in section 5.2. First, we note that our approaches and the unsupervised indices yield different sets of suspicious segmentations, as the IoU score is always less than 0.4. Furthermore, the unsupervised indices lead to inconsistent results on both datasets which make their performance highly unreliable on medical images.

One possible explanation is that those methods were designed for 2D color images with large contrast and may not scale well to 3D medical images.

If we now compare the different regularization strategies proposed for our approach, we observe that the level of regularization has some impact. Indeed, there is a significant variability of results with the value of β for the MRF prior.

This observation supports using the last regularization strategy proposed, the FDSP prior, as in this case all hyperparameters are learnt in a data-driven way.

Quantitative analysis

In order to perform a quantitative comparison of the different methods, we need to have a grading of the quality of all segmentations. As they are easier to obtain for real world images than medical images, we propose to conduct this quantitative assessment on the COCO dataset which contains real world pictures with a large variability among them, with grayscale or color images, segmented structures of variable sizes and large ranges of noise level. For each distribution, we extracted segmentations from the right tail corresponding to 20% of the total distribution. If the number of extracted segmentations was larger than 40, only the 40 cases with largest ASE/score were retained. All segmentations were pooled leading to a total of 703 delineations.

Quality grading process

Six raters were then recruited in order to grade each segmented image as good or poor through a custom application presenting the segmentations in a random order. Raters were asked to repeat the annotation twice in order to estimate the intra-rater variability. The intra-rater variability was found to be slightly lower than the inter-rater variability, with a mean rate of identical responses of 83% for the former and of 73% for the latter.

Performance comparison

The objective is to compare the segmentation quality among the right tails of the distributions given by the different approaches. These tails correspond to segmentations with the largest ASE or index score. The percentage of cases rated as poor by the raters strongly depends on the size of the set of segmentations extracted from the right tail of each distribution. Yet, it is useful to compare two quality control algorithms since a better algorithm is expected to have a greater proportion of segmentations annotated as poor by the raters than a worse one.

Each segmentation was assigned to a quality category, good or poor, after taking the mean across the raters' responses. Proportions of poor segmentations per approach were then derived for each object category. Distributions of these proportions over the 7 object categories are shown in Fig. 11. Two observations can be made. First, our approaches show competitive results as they lead to higher mean and median proportions of poor segmentations than the unsupervised indices. Second, no regularization strategy seems better than the others.

In particular, variations of the value of β do not affect the results very much for the MRF.

Medians Means

Figure 11: Distribution of the proportion of poor segmentations over the 7 object categories for each approach. The mean over the raters is taken as the final label for each segmentation.

taken into account, depending on the level of agreement among raters' responses.

Fisher's exact test is used to assess the difference between the proportion for a

given approach and the one obtained with an FDSP prior. Three unsupervised indices η, Zeb and GS, are found to give significantly different results than our approach with FDSP regularization, regardless of the threshold. More generally, our approaches always lead to a higher percentage of poor segmentations than the indices. Again, all regularization strategies seem to be appropriate. The results seem to be stable with respect to the level of regularization enforced by β.

Comparison with inter-rater variability

Assessing the quality of segmentations inside a medical imaging dataset is difficult without any expert knowledge. However, some datasets provide several segmentations of the same image produced by different experts. For instance, up to four segmentations are available for each nodule in the LIDC dataset and MS lesions in the MSSEG training dataset were delineated by seven radiologists.

The inter-rater variability measures the level of agreement between the experts.

It is reasonable to assume that images for which there is a low level of agreement Figure 12: Proportion of poor segmentations after pooling all object categories. Only segmentations with a raters' agreement above a given threshold are retained in the computation of the proportion. Results found to be significantly different from the ones given by the FDSP prior with Fisher's exact test and p-value 0.05 are marked with star symbols .

between the experts are more challenging than others. Therefore we study in this section the relationship between inter-rater variability and the score produced by our unsupervised model.

Tabs. 1 and 2 show the correlation coefficient between the inter-rater variability and the Dice score or ASE produced by our model on the LIDC and MSSEG datasets, respectively. We also compare with the four unsupervised indices selected earlier. The inter-rater variability was quantified in three manners: by computing the average Dice score between all pairs of experts, the average pairwise Hausdorff distance (HD) and the average 95% percentile of the pairwise Hausdorff distance (95% HD). It is compared to the average score computed on the different raters' segmentations for each unsupervised method. For the LIDC dataset, we discarded all nodules annotated by a single radiologist, leaving a total of 87 nodules.

Better correlations are achieved on the MSSEG dataset than on the LIDC dataset. Furthermore, our approach differs significantly from the other unsu- We further analyse the link with inter-rater variability by showing some examples from both datasets on Fig. 13. The first row presents results on the MSSEG dataset, where the correlation is quite good (0.81). Case A has a high inter-rater variability and is labelled as challenging by our model (low average Dice between the inputs and the model). Indeed, only three raters out of seven considered that some lesions were visible on the slice presented in Fig. 13b. Moreover, the low intensity contrast does not help to understand the segmentations. On the other hand, case B is better explained by the model with a good agreement between the experts, as shown in Fig. 13c.

The bottom row shows poorer results on the LIDC dataset. Two contradictory cases are highlighted. The first one, case C, has a low inter-rater variability but is predicted as challenging by our model (low average Dice score between S and M ). The two radiologists are indeed giving close contours (Fig. 13e) but it is also clear that the case is challenging according to the assumptions of our model. In the image regions highlighted by the arrows, the contours are indeed crossing areas of uniform intensity distribution, which make them more difficult to understand. On the other hand, case D is a typical case illustrating the limitations of our model (Fig. 13f). Raters disagree about the extent of the nodule, but all segmentations correspond to visible boundaries and match the assumptions of our model. One possible explanation for the poorer correlation obtained on the LIDC dataset is that the annotations were made in two stages, the second stage allowing radiologists to see the annotations made by the other experts in the first stage. This may have led to a decrease in inter-rater variability.

This analysis shows that in some cases, the inter-rater variability may not be a good surrogate of the difficulty of a segmentation. Raters may provide similar segmentations despite the fact that they are not close to visible boundaries (Case C in Fig. 13e) in the image. In that case, a low inter-rater variability is associated with a difficult segmentation.

Results interpretability

The previous section demonstrated how well our approach performed in extracting suspicious segmentations from a dataset in an unsupervised manner.

However, it also differs from approaches proposed in the literature regarding the output of the algorithm. For instance, the unsupervised indices output only a scalar score as a ratio of 2 metrics measuring the intra-region homogeneity and the inter-region dissimilarity. In our case, the output of the algorithm is a new segmentation used as a comparison tool. Although this segmentation must not be seen as a surrogate ground truth, it can help to visually understand why a segmentation is considered atypical, that is, has a large ASE, which is not possible with the indices.

Voxels lying on the input segmentation border can thus be colored depending on their distance to the model segmentation contour, as shown in Fig. 14. When dealing with 3D medical images with a large number of slices, it is useful to be able to retrieve quickly the most problematic regions according to the model.

Identifying the most suspicious slices is not possible with approaches outputting a simple score. Last but not least, the model segmentation could also be used as a guide for the correction of poor cases. regions have different mixtures of Student's t-distributions.

Discussion

The proposed unsupervised quality control method was shown to efficiently and automatically isolate challenging or atypical segmentation cases from a whole dataset. It was shown to outperform four previously introduced segmentation quality indices on the COCO dataset. Furthermore, those four indices do not provide stable results on the LIDC and BRATS medical datasets. The proposed algorithm does not produce a classification between good or poor segmentations but rather a ranking between cases within a dataset.

The genericity of the algorithm allows it to work on any type of object category or image (2D RGB or 3D grayscale images). We demonstrated the ability of the method to handle a wide range of segmentations, from small structures (lung nodules) to large brain tumor delineations. Yet, the approach is not suited for very tiny objects since a reasonable size is required to have a reliable estimation of the intensity parameters. Also, the spatial prior is likely to wipe out the segmentation if its area is really too small. Furthermore, the genericity of the algorithm may also be considered as a limitation when focusing on a specific structure of interest. For instance, if we aim at segmenting objects from the car category on the COCO dataset, a contour perfectly following intensity boundaries but around another object category would not be identified as atypical. To this end, one would need to also monitor several specific features of that structure such as its color, size or shape, which amounts to performing a supervised quality control as in [START_REF] Xu | Automated Quality Control for Segmentation of Myocardial Perfusion SPECT[END_REF]. This limitation is shared by all unsupervised quality control methods.

Another limitation is the difficulty to distinguish boundaries in areas with low intensity contrast. Our method is based on mixtures of Student's t-distributions, which is already a far more general assumption than some previous unsupervised approaches that hypothesize a unique Gaussian distribution in each region [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF]. Furthermore, our Bayesian formulation integrates intensity and smoothness assumptions into a single probabilistic model, as opposed to previous unsupervised methods, which require weighting of the heterogeneity and homogeneity terms.

Different spatial regularization strategies are proposed and tested in this paper. Quantitative assessment on COCO seems to indicate that all approaches lead to similar results. However, the FDSP prior based on derivative penalization does not require any hyperparameters to be set while keeping the computation time low, supporting its use in preference to the others.

Finally, compared to learning-based approaches such as [START_REF] Kohlberger | Evaluating Segmentation Error without Ground Truth[END_REF] or [START_REF] Robinson | Real-Time Prediction of Segmentation Quality[END_REF] and also to previous unsupervised indices which only output a score, our method provides an explanation for the mismatches between the posterior probabilities M and the input segmentation S. This is a major advantage considering the growing importance of providing interpretable models.

Conclusion

Image segmentation is an important task in medical image analysis and computer vision. Quality control assessment of segmentations is therefore crucial,

Figure 1 :

 1 Figure 1: Unsupervised segmentation quality control workflow.

  involves solving an inference problem depending on the choice of spatial prior. After convergence of the probabilistic model, a new segmentation M is generated by thresholding the posterior p(Z n |I n ) at the level 0.5.

Figure 2 :

 2 Figure 2: Graphical model of the framework with a discrete MRF prior (2a), a GLSP prior (2b) or a FDSP prior (2c).

Figure 3 :

 3 Figure 3: Quality control workflow on a glioblastoma segmentation from the BRATS 2017 dataset with FDSP regularization. (3a) Original image. (3b) Input segmentation S. (3c) Narrow band along the ground truth boundary with foreground (in red) and background (in blue) regions. (3d) Appearance ratio, label spatial prior and label posterior. (3e) Evolution of the lower bound J (q).

  The BRATS 2017 datasets consist of multisequence preoperative MR images of patients diagnosed with malignant brain tumors. It includes 285 patients for the training dataset and 46 for the validation set. Four MR sequences are available for each patient: T1-weighted, post-contrast (gadolinium) T1-weighted, T2weighted and FLAIR. All the images have been pre-processed: skull-stripped, registered to the same anatomical template and re-sampled to 1 mm 3 resolution. Ground truth segmentations of the brain tumors are provided only for the training set.

  analysed the stability of the sets made by the 40 segmentations with the largest ASE values by computing pairwise intersection over union (IoU) coefficients. A value of 1.0 indicates that the 40 images are the same for a pair of narrow band widths. The outcome is shown in Fig. 4.

Figure 4 :

 4 Figure 4: Sensitivity analysis of the narrow band width for the BRATS and LIDC datasets with our approach using FDSP as a spatial prior (first row) or with the unsupervised indicator GS (second row). Matrices show IoU scores computed between the sets made of the 40 segmentations with largest ASE (first row) or largest GS score (second row).

Figure 5 :

 5 Figure 5: Example of a nodule segmentation from LIDC where the result of the quality assessment is different depending on the narrow band width. If too large, the appearance model of the foreground leaks inside the pleura leading to an irrelevant result.

Figure 6 :

 6 Figure6: ASE distributions for the analysis of ground truth segmentations from the BRATS (6a) and LIDC datasets (6b). Samples from the left tail are identified as explained by the model while samples from the right tail are classified as challenging.

  whole brain tumor segmentations from the BRATS 2017 training set and pulmonary nodule segmentations from the LIDC dataset. On average, one minute is required to complete the quality control workflow for a 3D image from the BRATS dataset using an MRF or FDSP regularization. The inference time increases to 4 minutes for a model with a GLSP prior. The computation time of course also depends on the size of the segmented structure and on the extent of the narrow band. Computation of the ASE for each segmentation allows the distribution for the whole dataset to be drawn. Histograms obtained with FDSP regularization are shown in Fig. 6. They present a similar shape, with a short left tail, a single peak and a heavier right tail. Cases in the right tail isolated from the rest of the distribution are atypical and possibly include errors. Samples from the left and right tail are shown in Figs. 7 and 8, respectively. For both datasets, cases with larger ASE are clearly more challenging than the cases taken from the left tail.

Figure 7 :

 7 Figure 7: Segmentations with the smallest ASE taken from the left tail of the distributions. Cases are ranked according to their ASE value (Largest values to the right) and slices with largest ground truth area are shown. The width of the narrow band is 30 vx for BRATS and 10 vx for LIDC.

Figure 8 :Figure 9 :

 89 Figure 8: Segmentations with the largest ASE taken from the right tail of the distributions. Cases are ranked according to their ASE value (Largest values to the right) and slices with largest ground truth area are shown. The width of the narrow band is 30 vx for BRATS and 10 vx for LIDC.

Figure 10 :

 10 Figure 10: Examples of ground truth segmentations from the COCO dataset representing the 7 selected object categories.

Fig. 12

 12 Fig. 12 is obtained after pooling all object categories. To assess the robustness of the results, different thresholds are used to select the segmentations

Figure 13 :

 13 Figure 13: Correlation between the inter-rater variability and the difficulty of a segmentation as predicted by our model on the MSSEG dataset (top row) and LIDC dataset (bottom row).

Figure 14 :Figure 16 :

 1416 Figure 14: Interpretability of the result given by our approach on a brain tumor segmentation from BRATS. (14a) Ground truth segmentation and label posterior given by the model with FDSP regularization. (14b) Coloring of voxels lying on the ground truth border depending on their distance to the ouput of the probabilistic model.

E

  but the trend towards the generation of large databases makes any human-based monitoring onerous if not impossible. This paper introduces a new framework for generic quality control assessment which relies on a simple and unsupervised model. It has the advantage of not requiring a priori any knowledge aboutE[log p(I|Z)] = n η n1 log r n + η n0 log(1r n ) . (C.2) E[log F (Z, W , ξ)] = n η n1 E[w n ] + log σ(ξ n ) [log q Z (Z)] = n η n1 log η n1 + η n0 log η n0 .

Table 1 :

 1 Values of the correlation coefficient between the inter-rater variability and the average score given by different methods on the LIDC dataset. The width of the narrow band is 10 vx and FDSP was used as a spatial prior for our model.

	pervised indices with much larger correlation values. The others (except Zeb)
	exhibit indeed coefficients close to zero.

Table 2 :

 2 Values of the correlation coefficient between the inter-rater variability and the average score given by different methods on the MSSEG dataset. The width of the narrow band is 20 vx and FDSP was used as a spatial prior for our model.

Acknowledgments

This work was partially funded by the French government, through the UCA JEDI and 3IA Côte d'Azur "Investments in the Future" projects managed by the National Research Agency (ANR) with the reference numbers ANR-15-IDEX-01 and ANR-19-P3IA-0002 and supported by the Inria Sophia Antipolis -Méditerranée "NEF" computation cluster.

Surrogate segmentation performance

In this section we investigate if metrics estimated by our segmentation quality assessment algorithm can be correlated with the overall segmentation performance of an algorithm. In particular, we consider the segmentations generated However, we do not find that this approach always predicts the performance of segmentation algorithms well. For instance, we have noticed poor predictions for most categories in the COCO dataset. This can be explained by the fact that good performance predictions can only be obtained when the segmented structure follows the model assumptions, that is, the background and foreground the segmented objects nor a subset of trusted images to be extracted. This is especially suited to the monitoring of manually created segmentations, where potential errors can be found, as shown by our results. Its application to segmentations generated by algorithms is also of great interest and in some cases can be used as a surrogate for segmentation performance.

The proposed generic segmentation model produces contours of variable smoothness that are mostly aligned with visible boundaries in the image. Three regularization strategies were proposed in this paper and produced similar results. However, the prior based on derivative penalization has the great advantage of allowing an automated estimation of all hyperparameters with variational Bayesian inference, which is not possible within the classical MRF framework.

Extensive testing has been performed on different datasets containing various types of images and segmented structures, showing the ability of the method to isolate atypical cases and therefore to perform quality control assessment.

Comparison with unsupervised indices from the literature proved our approach to be effective and competitive. Coping with multiple foreground labels may be an interesting extension to process multiple regions of interest jointly rather than sequentially. Finally, an interactive use of the proposed algorithm during the manual delineation of structures in images is an exciting perspective to help reduce the inter-rater variability in the context of crowdsourcing.

Appendix A. Unsupervised indices

We give in this section the formula used to compute the unsupervised indices.

We denote by R the number of regions inside an image (typically 2 here, for the foreground and background regions). R j denotes the set of voxels in region j and |R j | is the number of voxels in region j. Each indicator requires the computation of an intra-region uniformity metric IU and an inter-region disparity metric ID.

Appendix A.1. Zeb [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF])

where W (s) is the neighborhood of voxel s and:

where b(R j ) is the set of pixels on the border of R j .

The final indicator is given by:

Appendix A.2. F RC [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF])

where:

where W (R j ) is the set of neighboring regions of R j and:

The final indicator is given by: [START_REF] Zhang | Image segmentation evaluation: A survey of unsupervised methods[END_REF] The background is denoted here by b, while f denotes the foreground.

where N b and N f are the number of voxels in the background and foreground, respectively, and e 2 (R j ) is defined as previously. .12) where ÎRj = 1 ν ν i=1

Îi Rj . The final indicator is given by:

Appendix A.4. GS [START_REF] Johnson | Unsupervised image segmentation evaluation and refinement using a multi-scale approach[END_REF])

where V j is the variance of region j.

The inter-region disparity metric used is the Global Moran's I, defined as:

where w ii = 0, w ij = 1 if R i and R j are neighbors and 0 otherwise. y i is the mean intensity value of region R i and y is the mean intensity value of the image.

The final indicator is given by: GS = IU + ID . (A.16)

Appendix B. FDSP prior -variational inference

We present in this section the derivation of the variational update formula (11), ( 12), ( 13), ( 14) and ( 15). The likelihood of the model p(I, Z, W, α) factorizes as p(I|Z)p(Z|W )p(W |α)p(α).

Recalling that p(I|Z) = n r Zn n (1r n ) 1-Zn , we have:

The prior p(Z|W

to give:

2) and (B.3) and taking the exponential, we have q * Z (Z) ∝ n ρ 1-Zn n0 ρ Zn n1 where the expressions of ρ n0 and ρ n1 are given by ( 11). With the normalization constraint, we finally obtain q * Z (Z) = n η Zn n1 η 1-Zn n0 with

With the expression of p(W |α) given in (5) and assuming that q W (W ) = n q wn (w n ), we obtain:

By identifying the quadratic and linear terms in w n , we obtain the formula for Σ wn and µ wn given in ( 12) and ( 13).

Assuming q * α (α) to be a Dirac distribution, we take the derivative of (B.6) with respect to α which leads to the update formula given in ( 14).

Appendix B.4. Update of q * ξn (ξ n )

Taking the derivative with respect to ξ n and setting it equal to zero gives λ (ξ n )(E[w 2 n ]ξ 2 n ) = 0. As λ (ξ n ) 0, we finally obtain the formula reported in (15).

Appendix C. FDSP prior -lower bound

The lower bound on the log-likelihood is used as a stopping criterion. To compute J (q), we need to evaluate the right hand side of (10):

The values of the different expectations can be computed and are reported below.


Direct Measurement of Ionization Charges in Single-phase Liquid Xenon Compton Telescope for 3γ Medical Imaging

Yajing Xing, Mounir Abaline, Stéphane Acounis, Nicolas Beaupère, Jean-Luc Beney, Julien Bert, Stéphane Bouvier, Clotilde Canot, Thomas Carlier, Hervé Chanal, Michel Cherel, Jean-Pierre Cussonneau, Sara Diglio, Debora Giovagnoli, Jérôme Idier, Franc ¸oise Kraeber-Bodéré, Patrick Le Ray, Frédéric Lefèvre, Samuel Manen, Julien Masbou, Eric Morteau, Nicolas Pillet, Didier Roy, Laurent Royer, Jean-Sébastien Stutzmann, Richard Vandaele, Dimitris Visvikis, Yuwei Zhu and Dominique Thers.

Abstract-We report the study of direct measurement of ionization charges in an innovative liquid xenon Compton camera of small animal imaging, named XEMIS2. It was combined with a novel 3γ medical imaging modality, showing a "TOF-like" PET performance, targeted to reduce the administered activity diametrically while preserving the image quality in oncology diagnosis. For the Compton cones reconstruction, the sequential low energy electronic recoils from Compton scattering of third γray (∼ 1 MeV) are measured directly from ionization signals under a high electric field of 2 kV/cm. The Geant4/NEST simulations studies reveal that the intrinsic energy and spatial resolutions of charges carriers measuring in liquid xenon are dominated respectively by the electron-ion recombination and the abstruse trajectory of recoil electrons. Besides, a novel ASIC front-end electronic, XTRACT, and data acquisition chain were developed, dedicated to optimizing the accuracy and efficiency of ionization signals measurement in XEMIS2. It allows a continuous read-out with ultra-low charges threshold and negligible deadtime per individual channel. The performing characteristics were calibrated in the prototype XEMIS1 under the XEMIS2 operation condition, showing a good homogeneous charge linearity response. Meanwhile, the data processing and analysis were developed, where the charge measurement has been optimized through Monte Carlo simulations.

I. INTRODUCTION

T HE personalized medicine as the heart of healthcare priorities of the 21 st century requires molecular imag- ing having better sensitivity and resolution to reduce the administered dose and to shorten the exam time. There are three current evaluations of the Positron Emission Tomography (PET) imaging detectors. The first type is the total-body PET. It reduces the exam time dramatically with the same dose by increasing the field-of-view (FOV), which is highly developed in US [START_REF] Badawi | First human imaging studies with the EXPLORER total-body PET scanner[END_REF]. Another approach much explored in the EU is the time-of-flight PET (TOF-PET), whose aim is to decrease the effective length of a line of response (LOR) through accurate time measurement [START_REF] Lecoq | Pushing the Limits in Time-of-Flight PET Imaging[END_REF]. The third modality aimed at reducing the parallax effects through a depth of interaction on the whole FOV [START_REF] Yoshida | Development of a whole-body dual ring openpet for inbeam pet[END_REF]. All those sophisticated detector technologies are based on solid scintillators.

Our team proposed an innovative alternative liquid XEnon Medical Imaging System (XEMIS) [START_REF] Grignon | Nuclear medical imaging using β+ γ coincidences from Sc-44 radio-nuclide with liquid xenon as detector medium[END_REF]. Liquid xenon (LXe), as an emerging monolithic radiation detection medium, has been widely used for particle physics and dark matter detection. It possesses lots of advantageous features, such as high stopping power, ethical sensibility, and simultaneous scintillation light and ionization charge production. Meanwhile, its efficient scalability affords the possibility to develop a total-body imaging detector. Moreover, a "TOF-like" PET performance could be achieve through the application of an innovative 3γ imaging modality [START_REF] Cussonneau | 3γ Medical Imaging with a Liquid Xenon Compton Camera and 4 4Sc Radionuclide[END_REF]. It is based on a specific (β + , γ) radio-isotope 44 Sc [START_REF] Huclier-Markai | Promising Scandium Radionuclides for Nuclear Medicine: A Review on the Production and Chemistry up to In Vivo Proofs of Concept[END_REF] and a LXe Compton telescope, measuring with good precision the intersection of LOR and Compton cone to realize a direct 3D location of radio-tracers.

The first R&D prototype XEMIS1 [START_REF] Gallego Manzano | XEMIS: A liquid xenon detector for medical imaging[END_REF] has been accomplished and showed promising results. A larger Compton camera XEMIS2 [START_REF] Manzano | XEMIS2: A liquid xenon detector for small animal medical imaging[END_REF] for small animal imaging is developed to demonstrate the advantage of liquid xenon (LXe) as a detection medium for 3γ imaging. The objective is to obtain a highquality image of small animals with only 20 kBq activity for 20 minutes. The sequential low energy electronic recoils from Compton scattering of third γ -ray (∼ 1 MeV) are measured directly from ionization signals under a high electric field of 2 kV/cm. The accurate measurement of ionization charges carriers has thus a vital impact on the Compton cones reconstruction.

II. MOTIVATION OF CHARGE MEASUREMENT

The design of the XEMIS2 Compton camera is shown as Fig. 1. It is composed of two identical single-phase monolithic liquid xenon time projection chambers (LXe-TPC) sharing the same central cathode. Each LXe-TPC has a 12 cm drift length with the radius varies from 7 cm to 19 cm. The periphery of the active zone is covered by 64 one-inch PhotoMultiplier tubes to detect the prompting scintillation light. Two-segmented anodes, possessing a total of twenty thousand pixels with a pitch of 3.1 mm, are equipped on both sides to collect the ionization charge carriers. In XEMIS2, the scintillation signals mainly provide the γrays interaction time and the spatial pre-localization of γ-ray interactions to reduce the occupancy of the TPC [START_REF] Zhu | XEMIS2: A liquid xenon Compton camera to image small animals[END_REF]. The deposited energy and position of each interaction vertex are measured by ionization signals. The ionization charges carriers are directly collected by the segmented anodes which are shielded with Frisch grids, in our single-phase LXe-TPC, without through a gaseous phase electroluminescence amplification.

The accurate measurement of ionization charges carriers plays a vital role in the Compton cones reconstruction, which provides the possibility of Compton interactions sequential analysis and discrimination. Compared with the high cost of PMT/SiPM, the direct ionization signal measurement also has a tremendous economic advantage for the medical imaging application.

III. STUDY OF ENERGY AND SPATIAL RESOLUTION IN LXE

The performing energy and spatial resolutions of the XEMIS2 are mainly defined by the intrinsic limits generated by the physic process in LXe and the instrumental measurement method. For example, the spatial resolution of X-Y position reconstruction from the instrument is about 100 µm, which is dominated by the pixel size. In this part, we discussed only the intrinsic limits given by LXe.

The precision of charge carriers measurement in LXe-TPC is intrinsically limited by the electrons capture from the electronegative impurities and the electron-ion recombination, leading to attenuation and statistical fluctuation of charge carriers collection.

With our scintillation signal measurement strategy [START_REF] Zhu | Scintillation Signal in XEMIS2, a Liquid Xenon Compton Camera with 3γ Imaging Technique[END_REF], the ionization/scintillation yield anticorrelation will not be applied to correct the recombination influence effect. An intrinsic energy resolution has thus been carried out. In order to optimize the resolution of ionization signals measurement, a high electric field of 2 kV/cm will be applied in XEMIS2 to increase the ionization yield. Besides, this fluctuation of electrons collection from the electron-ion recombination in LXe can be estimated with NEST [START_REF] Szydagis | NEST: A Comprehensive Model for Scintillation Yield in Liquid Xenon[END_REF]. A Geant4/NEST simulation has been done to evaluate this limit, shown in Fig. 2. To eliminate the impact of charges attenuation, a very high purity level of xenon is needed since the electronegative impurities dissolved in LXe may absorb the electrons during the drift process. The LXe contamination in the XEMIS1 prototype limited to a ppb level through the purification system [START_REF] Xing | XEMIS: Liquid Xenon Compton Camera for 3γ Imaging[END_REF]. The experimental results demonstrated that the attenuation influence in ultrahigh purity LXe could be neglected before the recombination effect. In XEMIS2, the same/better purity level is expected to be achieved through the new purification-cryogenic subsystem ReStoX [START_REF] Virone | Gravity assisted recovery of liquid xenon at large mass flow rates[END_REF] equipped a gaseous purification system, realizing a re-circulation closed-loop. The ReStoX and XEMIS2 cryostat are installed in the small animal medical imaging center CIMA at Nantes University Hospital. The current statue is presented in Fig. 3. The intrinsic spatial resolution in LXe is dominated by the distribution of ionized electrons cloud produced from electronic recoil. This effect has also been considered and studied by Geant4/NEST simulation, which presents a spatial measurement influence with the energy of electron recoil (about 100 µm at 511 keV).

IV. IONIZATION SIGNALS MEASUREMENT SYSTEM OF XEMIS2

To achieve this ultra-low activity imaging of small animals, we need to consider three critical issues concerning ionization measurement: resolution optimization, event efficiency maximization and high-frequency signals readout with processable data flow due to low threshold triggering. The published performance calibration of XEMIS1 presented promising results with a roughly 100 µm position resolution along the drift coordinate (Z-position) and a good energy resolution of 5% (9% FWHM) for 511 keV photoelectric events under an electric field of 1kV/cm [START_REF] Gallego Manzano | XEMIS: A liquid xenon detector for medical imaging[END_REF].

To extend the same/better performance from 64 pixels to more than 20000 pixels under a higher electric field of 2 kV/cm has still several technical hurdles to overcome. Moreover, a considerable quantity of data caused by the ultra-low threshold and self-triggering circuit needs to be managed. For the sake of continuous read-out with negligible dead-time during 20 minutes, the first version of new developed higher capability Data Acquisition (DAQ) chain dedicated to XEMIS2 is accomplished.

A. Front-end electronics and data acquisition chain

To optimize the event efficiency of XEMIS2 with the negligible dead time, we chose a low threshold self-triggered data acquisition mode for ionization measurement. Besides, in light of the massive number of pixels, the idea of signal information extraction has been brought out to reduce the readout data flow. Based on the above considerations, a novel self-triggered scalable signal readout array with ultra-low charges threshold has been developed in the Subatech laboratory with collaboration.

It consists of four function units: IDeF-X HD-LXe [START_REF] Lemaire | Development of a Readout Electronic for the Measurement of Ionization in Liquid Xenon Compton Telescope Containing Micro-patterns[END_REF], XTRACT, XEMIS Process Unit (XPU) and XEMIS Data Concentrator (XDC). IDeF-x has been validated in XEMIS1, and all the other parts are novel developments, especially XTRACT. It is a multi-channel ASIC (×32) with analog and common digital blocks. Each channel includes a Constant Fraction Discrimination (CFD) block which performs the selftriggering function with a limited time walk and an improved timing resolution. XTRACT memorizes only pixel address, amplitude and CFD time of each over-threshold ionization signal. Then, each XPU could connect a maximum of 16 XTRACTs realizing data collection and analog-to-digital conversion. The XDC consists of a Field-Programmable Gate Array (FPGA) component, whose fundamental function is to output the binary data and to synchronize ionization and scintillation DAQ chain. A real-time data filter is also proposed to be implemented in further development. XEMIS2 equipped 640 IDeF-Xs and XTRACTs, 64 XPUs and two FPGAs to read out all ionization signals of 20000 pixels with negligible dead time. This electronic chain is currently tested in XEMIS1 under the XEMIS2 operation condition for precise characterization before integration in XEMIS2.

B. Performing characterization of charge measurement

The performing characteristics of the recent version of the novel ASIC/DAQ chain were calibrated in the prototype XEMIS1 under the XEMIS2 operation condition. The experimental results have shown that all tested 64 channels of two XTRACTs had a good charge linearity response and quite homogeneous behavior, as presented in Fig. 4. The electronic noise at the DAQ level of each channel is quite uniform and low (about 150 electrons) at the first attempt. Besides, the electronic time resolution and the CFD time walk have also been measured. The DAQ intrinsic time measurement RMS is less than 10 ns, which is limited by FPGA sampling at 100 MHz. Besides, results showed that walk correction on the small measured charge is needed. All 800 XTRACTs have been fabricated and the current bench test at room temperature shows a high rate of qualified to 99%.

V. DATA PROCESSING AND CHARGE MEASUREMENT

OPTIMIZATION

An ultra-low threshold used for self-triggering has been applied to optimize the charge measurement. Hence, the raw data contents abundant over-threshold noise signals during the acquisition. To realize the final interaction vertex reconstruction, we developed a dedicated data processing. After binary data, the electronic calibration results were firstly applied for the sake of pedestal subtraction, gain application and CFD time walk correction. Then, an algorithm has been programmed on soft to realize a data filter in time, reducing the raw data by a factor of 80 (most of the time-isolated noise signals with few ultra-low under-threshold physic signals). We attempt to achieve a real-time data filter in the further R&D step of XDC. Finally, the practiced spatial Clustering is reconstructed with a time window of 2 µs, defined by time resolution and time window studies.

A specific detector response Monte Carlo simulation has been done to improve the charge measurement accuracy, concerning the clustering method optimization and Compton sequence reconstruction. The results showed that the ionization signal measurement resolution could be degraded by the mix of two interactions, due to the ionized electrons cloud size, pixel size and electronic integration time. A further finer study is under going. Besides, we found that the threshold effect of clustering resulted in a charge loss, presenting a dependence on charge collected number and z position. A simulation result is shown in Fig. 5, who has been appliqued during the data analysis.

VI. CONCLUSION

The accurate Compton scattering reconstruction is a core issue for the success of 3γ imaging in XEMIS2. To estimate the ultimate performance, we explored the intrinsic performances of LXe-TPC with Geant4/NEST simulations. Besides, a novel low threshold self-triggered ionization signals measurement system dedicated to XEMIS2 has been currently accomplished and tested in XEMIS1 who presents a good linearity charges response at LXe temperature. All 800 XTRACTs have been fabricated and the first bench test at room temperature shows a high rate of qualified to 99%. Meanwhile, the associated data processing has also been developed. Meanwhile, the data processing and analysis were developed, where charge measurement has been optimized through Monte Carlo simulations. A new approach to Compton scattering angular calibration is underway for further study.

Fig. 1 .

 1 Fig. 1. Design of the XEMIS2 camera for small animal imaging.

Fig. 2 .

 2 Fig. 2. Intrinsic contribution estimated using NEST comparing with the detector contribution estimated with the simulation of the electronic response. The recombination effect dominates the measurement resolution of large ionization signals (carried charges bigger than 5000 electrons).

Fig. 3 .

 3 Fig. 3. Installation of the cryogenic subsystem ReStoX and XEMIS2 cryostat in the small animal medical imaging center CIMA at Nantes University Hospital.

Fig. 4 .

 4 Fig. 4. Experimental charge linearity response of 64 channels of two XTRACTs under XEMIS2 operation condition.

Fig. 5 .

 5 Fig. 5. Mean charge loss with clustering charge and z position estimated by simulation under the conditions: considering the effects of electron-ion recombination and charge diffusion in LXe; simulating electronics response; neglecting the recoil electron trajectory; under an homogeneous electric field of 2 kV/cm; neglected attenuation influence; equivalent noise charge per pixel equal to 150 electrons.

  Manuscript received December 13, 2019. The research presented in this paper has been funded by the E.U., the Region Pays de la Loire in France and by grants from the French National Agency for Research.

	Yajing Xing is with SUBATECH, IMT Atlantique, CNRS/IN2P3, Uni-
	versité de Nantes, 44307 Nantes, France (corresponding author e-mail: ya-
	jing.xing@subatech.in2p3.fr).
	Stéphane Acounis, Nicolas Beaupère, Jean-Luc Beney, Stéphane Bouvier,
	Clotilde Canot, Jean-Pierre Cussonneau, Sara Diglio, Patrick Le Ray, Frédéric
	Lefèvre, Julien Masbou, Eric Morteau, Jean-Sébastien Stutzmann, Yuwei Zhu
	and Dominique Thers are with SUBATECH, IMT Atlantique, CNRS/IN2P3,
	Université de Nantes, 44307 Nantes, France.
	Mounir Abaline, Hervé Chanal, Samuel Manen, Nicolas Pillet, Laurent
	Royer, Richard Vandaele are with LPC Clermont-Ferrand, 24 Avenue des
	Landais, Clermont-Ferrand France
	Julien Bert, Debora Giovagnoli and Dimitris Visvikis are with INSERM,
	UMR1101, LaTIM, CHRU Morvan, 29600 Brest, France.
	Thomas Carlier and Franc ¸oise Kraeber-Bodéré are with Centre Hospitalier
	Universitaire de Nantes, 44093 Nantes, France.
	Michel Cherel is with INSERM U892 équipe 13, 44000 Nantes, France.
	Jérôme Idier is with LS2N, Ecole Centrale de Nantes, CNRS/Inp, Université
	de Nantes, 44307 Nantes, France.

978-1-7281-4164-0/19/$31.00 c 2019 IEEE