Wang Zihao
email: zihao.wang@inria.fr

Demarcy Thomas

Vandersteen Clair

Gnansia Dan

Raffaelli Charles

Guevara Nicolas

Delingette Hervé

Bayesian Logistic Shape Model Inference : application to cochlea image segmentation

Keywords: Bayesian Inference, Image Segmentation, Shape Modeling

Incorporating shape information is essential for the delineation of many organs and anatomical structures in medical images. While previous work has mainly focused on parametric spatial transformations applied on reference template shapes, in this paper, we address the Bayesian inference of parametric shape models for segmenting medical images with the objective to provide interpretable results. The proposed framework defines a likelihood appearance probability and a prior label probability based on a generic shape function through a logistic function. A reference length parameter defined in the sigmoid controls the trade-off between shape and appearance information. The inference of shape parameters is performed within an Expectation-Maximisation approach where a Gauss-Newton optimization stage allows to provide an approximation of the posterior probability of shape parameters. This framework is applied to the segmentation of cochlea structures from clinical CT images constrained by a 10 parameter shape model. It is evaluated on three different datasets, one of which includes more than 200 patient images. The results show performances comparable to supervised methods and better than previously proposed unsupervised ones. It also enables an analysis of parameter distributions and the quantification of segmentation uncertainty including the effect of the shape model.

Introduction

Several anatomical structures have a typical shape, such that a medical expert can easily recognize them from their three-dimensional representation. This is for instance the case of basal ganglia within the brain (Ashburner and Friston, 2005a), but also of abdominal structures, such as the liver or kidneys. Another emblematic example is the cochlea which is a small organ within the inner ear having a remarkable spiraling configuration where mechanical waves are transformed into electrical stimulation of the auditory nerve. The cochlea shape is complex as it completes around two and a half turns with its centerline closely resembling a logarithmic spiral helix [START_REF] Cohen | Improved and simplified methods for specifying positions of the electrode bands of a cochlear implant array[END_REF][START_REF] Baker | Tracking, modelling and registration of anatomical objects: the human cochlea[END_REF]. Its segmentation from CT images of the temporal bone is challenging since those images have low resolution with respect to the anatomy of the cochlea: the cochlea dimension is about 8.5x7x4.5 mm 3 while the typical CT voxel size is larger than 0.2 mm which is weakly visible for the fine structures of the chambers.In addition, the cochlea is filled with fluids that can be found in the vestibular system and other neighbouring structures, with similar appearance in CT images.

Supervised learning (e.g. Deep Learning) is an effective way to perform image segmentation or processing in many cases. Specifically, in inner ear CT imaging analysis, many works achieved impressive results [START_REF] Lv | Automatic segmentation of temporal bone structures from clinical conventional ct using a cnn approach[END_REF][START_REF] Raabid | Automatic segmentation of inner ear on ct-scan using auto-context convolutional neural network[END_REF][START_REF] Heutink | Multi-scale deep learning framework for cochlea localization, segmentation and analysis on clinical ultra-highresolution ct images[END_REF][START_REF] Wang | Deep learning based metal artifacts reduction in post-operative cochlear implant ct imaging[END_REF][START_REF] Li | Micro-ct synthesis and inner ear super resolution via generative adversarial networks and bayesian inference[END_REF][START_REF] Alshazly | Ensembles of deep learning models and transfer learning for ear recognition[END_REF][START_REF] Zhang | Few-shot learning for ear recognition[END_REF]Wang et al., 2020b). However, supervised learning methods have also many limitations. First, creating dataset annotations is time consuming, possibly preventing the creation of massive training datasets. In the cochlea case, a well trained ENT surgeon would need at least ten minutes to segment each 3D cochlea volume. Second, due to the potential overfitting related to the limited training set, the output of such supervised algorithm is likely to fall outside the shape space of the structure of interest.

Shape-based image segmentation can overcome the above limitations since the optimization of the model can be done in an unsupervised or weakly supervised way. Besides, the recovered shape parameters make a natural compact representation that is useful for shape analysis and even clinical applications. In this paper, we consider shapes that are either defined as an explicit S(θ S ) ∈ R d or implicit S(θ S , x) = 0 parametric shape models where θ S is a set of shape parameters and x ∈ R d , is any point in space (d = 2, 3).

Those parametric shape models serve to guide the delineation of such anatomical structures by constraining the shape space of the segmented object. We can roughly split the shape-based image segmentation methods into two sets of methods. A first set optimizes the shape parameters θ S by minimizing the sum of a regularizing term E R (θ S ) and an image term E I (S(θ S ), I, θ I ) : θS = arg min θ S E I (S(θ S ), I, θ I ) + E R (θ S ) where θ I is a set of image parameters that may also be optimized. This iconic shape fitting principle is typically used in the classical active shape model [START_REF] Cootes | Active Shape Models-Their Training and Application[END_REF][START_REF] Heimann | A shape-guided deformable model with evolutionary algorithm initialization for 3D soft tissue segmentation[END_REF] and their extensions [START_REF] Cremers | Shape Statistics in Kernel Space for Variational Image Segmentation[END_REF]. Various generic image terms may be considered for instance as those explored in [START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF]. A second set of methods uses the shape model S(θ S ) as a shape prior instead of a shape space. Several shape constraints have been introduced within several image segmentation frameworks including level-sets [START_REF] Chan | Level set based shape prior segmentation[END_REF][START_REF] Cremers | A variational framework for image segmentation combining motion estimation and shape regularization[END_REF], free-form deformation space (Rueckert et al., 2003b) or implicit template deformation [START_REF] Prevost | Incorporating shape variability in image segmentation via implicit template deformation[END_REF]. While those methods have greater shape flexibility for delineating structures, it is often difficult to set the coefficients weighting the shape constraint with other image terms. Those two sets of shape based segmentation methods are expressed as energy minimization problems, thus only allowing to have point estimates of shape parameters and not their posterior probabilities.

Another common shape representation consists in specifying a parametric spatial transformation T (θ D ) : R d → R d acting on a template shape S(θ 0 ) ∈ R d leading to an indirect shape parameterization : S(θ D ) = T (θ D ) • S(θ 0 ). This formulation of shape modeling based on a deformable template leads to solving a joint segmentation and registration problem. More precisely, several authors (Ashburner and Friston, 2005b;Pohl et al., 2006a) defined generative image and shape models and performed statistical variational inference to optimize their parameters and hyperparameters. Priors on the deformation space based for instance on minimal elastic energy [START_REF] Van Leemput | Encoding probabilistic brain atlases using bayesian inference[END_REF], were applied on triangular or tetrahedral mesh templates. Other shape priors were defined as restricted Boltzmann machines [START_REF] Agn | A modality-adaptive method for segmenting brain tumors and organs-at-risk in radiation therapy planning[END_REF] or as shape-odds [START_REF] Elhabian | Shapeodds: Variational bayesian learning of generative shape models[END_REF]. In most cases, optimal shape parameters (e.g. mesh vertex positions) are obtained as maximum a posteriori but not their posterior probability. Uncertainty quantification of image registration algorithms has been tackled in some research papers [START_REF] Simpson | Probabilistic inference of regularisation in non-rigid registration[END_REF][START_REF] Wang | Efficient laplace approximation for bayesian registration uncertainty quantification[END_REF][START_REF] Folgoc | Quantifying registration uncertainty with sparse bayesian modelling[END_REF] based on a low dimensional representation of deformation space and Laplace approximation.

In this paper, we propose a novel Bayesian framework for shape constrained image segmentation based on parametric shape models (instead of parametric spatial transformations) where the output segmentation is driven by a shape model but without restricting it to a low dimensional space. The proposed approach is generic as it is suitable for any explicit S(θ S ) and implicit S(θ S , x) = 0 parametric shape models associated with any appearance models representing the intensity distributions inside background or foreground regions. It is based on a logistic shape prior defined as the sigmoid of a shape function (e.g. signed distance map) defined over the image domain. Inferences of shape and intensity parameters are performed by maximizing the joint image and shape parameters probability p(θ S , θ I , I) with an Expectation-Maximization algorithm. We show that this optimization boils down to having the posterior label distribution as close as possible (in terms of Kullback Leibler divergence) from both the likelihood and shape prior distributions. A Gauss-Newton optimization method is introduced to optimize the shape parameters leading to closed form updates similarly to iterative reweighted least squares schemes. It outputs the most probable shape and imaging parameters but also an approximation of the posterior shape parameter probability which is essential for estimating the segmentation uncertainty.

This framework is applied to the problem of cochlea segmentation on CT images based on a parametric shape model with 10 parameters, and an imaging model defined as a mixture of Student's t-distributions. It results in the reconstruction of cochlea structures in 2 small datasets consisting of paired CT and µCT post-mortem images and one large dataset of nearly 200 patient CT images. We showed that the proposed framework leads to state of the art reconstruction performances as well as the recovery of consistent shape parameter distributions and the estimation of segmentation uncertainty.

The main contributions of this article are:

-A novel framework for image segmentation that combines probabilistic appearance and shape models. It is generically defined for parametric shape functions rather than parametric space transformations. The trade-off between the appearance and shape models is governed by an interpretable parameter : the reference length.

-A Gauss-Newton optimization method of the shape parameters which also produces a posterior approximation of those shape parameters.

-A method for uncertainty quantification of image segmentation which takes into account the shape uncertainty.

-A segmentation method of the cochlea in clinical CT images which provides state-of-the-art results and interpretable shape parameters.

We present below the framework of the logistic shape model (section 2), the shape and intensity models used specifically for cochlea segmentation (section 3), and the segmentation results on 3 clinical and pre-clinical datasets (section 4).

Method

Shape-based Generative Probabilistic Model

We consider an observed image I consisting of N voxels I n ∈ R, n = 1, . . . , N, for which we seek to solve a binary segmentation problem guided by a shape model. That model is defined either as in a parametric form as S(θ S ) ∈ R d , d = 2, 3 or in an implicit form as S(θ S , x) = 0. In the case of parametric shape models, one can define an associated implicit function SDM(S(θ S ), x) = 0 as the signed distance map defined at point x. Therefore, we propose to unify notations for both parametric and implicit cases by stating the existence of a shape function S(θ S , x) ∈ R whose zero level defines a shape and whose sign indicates if a point is inside (positive) or outside (negative). Note that with this hypothesis, a shape corresponds to a (smooth) manifold of co-dimension 1 without borders, thus defining a partition of the image into inside and outside regions.

A binary label variable Z n ∈ {0, 1} is defined at each voxel specifying if voxel n belongs to the background or foreground regions. A probabilistic intensity distribution model is defined for each region p(I n |Z n = k, θ k I ), k = 0, 1 controlled by the intensity parameter array θ k I . The arrays for background (k = 0) and foreground (k = 1) are concatenated into the intensity parameter array θ I . This appearance model can be either supervised , e.g. a trained convolutional neural network, or unsupervised, e.g. a Gaussian mixture model. In the remainder, we assume the latter case and therefore we define mechanisms to optimize the appearance parameters θ I . In the supervised case, the steps involving the update of θ I should be ignored.

We enforce a spatial correlation between the label of each voxel by specifying their a priori dependence on the shape model S(θ S , x). More precisely, we define a prior probability for voxel n to belong to the foreground region as follows:

p(Z n = 1|θ S ) = σ S(θ S , x n ) l ref (1) p(Z n = 0|θ S ) = 1 -p(Z n = 1|θ S ) = σ - S(θ S , x n ) l ref
where σ(x) is the sigmoid (or logistic) function, x n is the position of voxel n and l ref is a reference length. With that definition, the prior probability will be close to 1 inside the object, close to 0 outside and equal to 0.5 on the shape boundary. We call this formulation of the label prior, the logistic shape model as it combines shape information into a probability distribution through a logistic function. This definition of the shape prior is related to several prior work in the literature such as probabilistic atlases and LogOdds maps (Pohl et al., 2006b), continuous STAPLE [START_REF] Commowick | A continuous staple for scalar, vector, and tensor images: An application to dti analysis[END_REF], a nd label fusion [START_REF] Sabuncu | A generative model for image segmentation based on label fusion[END_REF].

The quantity l ref is a characteristic length which controls the slope of the prior probability next to the object boundary. This parameter also influences the trade-off between intensity and shape information in the segmentation process as discussed in section 2.5. The shape parameters θ S are themselves regarded as random variables with a multivariate Gaussian prior controlled by hyper-parameters α: p(θ S |α). The intensity parameters may also optionally be considered as random variables with hyperparameter β as p(θ I |β). The shape based generative model is summarized in Fig. 3:(a).

Logistic Shape Model Framework

With the proposed generative model, given an image, the objective is to infer the most probable values of the intensity θI and shape parameters θS which will lead to the estimation of the posterior label probabilities given by :

p(Z n = 1|I n , θ I , θ S ) = p(I n |Z n = 1, θ 1 I )p(Z n = 1|θ S ) 1 k=0 p(I n |Z n = k, θ 1 I )p(Z n = k|θ S ) (2)
That posterior probability is clearly a compromise between shape information stored in the prior p(Z n = 1|θ S ) and appearance information stored into the likelihood p(I n |Z n = 1, θ 1 I ). The segmented region of interest (SROI) then corresponds to voxels for which p(Z n = 1|I n , θ I , θ S ) ≥ 1 2 . In addition, the logistic shape model framework recovers the most likely shape parameters θS that corresponds to the segmented shape instance (SSI) which is the best fit of the shape model in that image. Finally, we will show that we can approximate the posterior shape parameter p(θ S |I) in order to capture the uncertainty in the shape parameter estimation.

The optimization of the intensity and shape parameters is done by maximizing the the log-joint intensity and parameters probability :

( θS , θI ) = arg max θ S ,θ I log p(I, θ S , θ I ) = arg max θ S ,θ I L(θ S , θ I ) L(θ S , θ I ) = log p(I|θ S , θ I ) + log p(θ S ) + log p(θ I ) = N n=1 log         1 k=0 p(I n |Z n = k, θ I ) p(Z n = k|θ S )         + log p(θ S ) + log p(θ I ) (3)
In the log-joint probability L(θ S , θ I ) we have marginalized out the hidden label variables Z n and used the conditional independence of variables I n given θ S .

Expectation-Maximization Inference

The direct optimization of L(θ S , θ I ) can be done by any optimization toolbox but it is difficult due to the possible encountered overflows/underflows caused by the logsum-exp expressions.

This is why we propose to follow the Expectation-Maximization (EM) algorithm which relaxes that optimization problem into several optimizations over simpler problems. We proceed by introducing N variables u n that are surrogates for the posterior label probability p

(Z n = 1|I n , θ S , θ I ) such that u n ∈ [0, 1]. Writing U = {u n }, we intro- duce a new augmented criterion L * (θ S , θ I , U) = log p(I, θ S , θ I )-D KL (U||p(Z|I, θ S , θ I ))
by adding the negative Kullback-Leibler divergence between u n and the posterior label p(Z n |I n , θ S , θ I ).

Maximizing (θ S , θ I , U) over the augmented criterion L * (θ S , θ I , U) leads to the same optima in (θ S , θ I ) than the maximization of L(θ S , θ I ) but with simpler expressions:

L * (θ S , θ I , U) = N n=1 1 k=0 u k n log(p(I n |θ S , θ I )p(Z k n = k|I n , θ S , θ I )) - N n=1 1 k=0 u k n log u k n + log p(θ S ) + log p(θ I ) = Q(U, θ S , θ I ) + N n=1 H(u n ) + log p(θ S ) + log p(θ I )
where

Q(U, θ S , θ I ) = E U (log p(I, Z|θ S , θ I )) is the conditional expectation of the com- plete marginal log-likelihood (a.k.a. evidence) and H(u n ) is the entropy of variable u n . The quantity Q(U, θ S , θ I ) is a lower bound of the log-likelihood since H(u n ) > 0.
The maximization of the augmented criterion L * (θ S , θ I , U) is performed by the successive maximization over the U, θ I and θ S variables. The E-step corresponds to the maximization of L * (θ S , θ I , U) with respect to U which sets the surrogate variable U to the posterior label probability

u n = p(Z n = 1|I n , θ S , θ I ).
The MI-step optimizes the log-joint probability with respect to the appearance variables θ I , which is equivalent to the maximization of L I = -D KL (U||p(I|Z, θ I )) + log p(θ I |β). When the appearance parameters are independent between classes, then log p(θ I |β) = K k=0 log p(θ k I |β k ) and the MI-step splits into 2 independent maximization

over θ k I , k = 0, 1 of L k I = -N n=1 D KL (u k n ||p(I n |Z n = e k , θ k I )) + log p(θ k I |β k ).
For certain well chosen intensity models such as Gaussian mixture models, this optimization leads to closed-form updates of θ I .

Finally, we perform the MS-step corresponding to the maximization over shape variables θ S which is equivalent to the maximization of L S :

L S = -D KL (U||p(Z|θ S )) + log p(θ S |α)
We can see that the EM algorithm preserves an interesting symmetry between shape and appearance information. Indeed, the iterative application of the E, MS and MI steps makes the posterior labels distribution U as close (in terms of KL divergence) as possible from the likelihood p(I|Z, θ I ) and shape prior p(Z|θ S ) that the minimization of D KL (U||p(Z|θ S )) + D KL (U||p(I|Z, θ I )). At convergence, the posterior distribution is therefore clearly a compromise between shape and appearance information.

Optimization of shape parameters p(θ S |I)

The functional L S is a non trivial function of the parameters θ S as it combines 2 non-linear functions : the sigmoid σ() and the shape function S(θ S , x n ):

L S = - N n=1 u n log σ S(θ S , x n ) l ref + (1 -u n ) log σ - S(θ S , x n ) l ref + log p(θ S |α) + cst (4)
The functional gradient ∇ θ S L S cannot be written in closed form since it requires the computation of the gradient of the scaled shape function at each voxel :

d n = ∇ θ S S(θ S ,x n ) l ref ∈ R |θ S | .
Those gradient vectors may be computationally costly to compute, for instance when the shape function is based on a signed distance map of parametric shape models S(θ S , x) = SDM(S(θ S ), x). In that case, the d n values are computed by a costly finite difference approximation except for translation and rotation parameters for which they can be computed efficiently (see Appendix A). After combining all d n terms in a gradient matrix d ∈ R |θ S |×N , the functional gradient can be simplified as

∇ θ S L S = -d(u-µ)+ ∇ θ S log p(θ S |α) where u = (u 1 . . . u N ) T ∈ R N and µ = σ S(θ i S ,x 1 ) l ref . . . σ S(θ i S ,x n ) l ref T ∈ R N .
Thus, a first approach for optimizing the shape parameters is to use any quasi-Newton optimization method such as the BFGS algorithm ( similarly to [START_REF] Demarcy | Segmentation and study of anatomical variability of the cochlea from medical images[END_REF]), since it only requires the computation of the functional gradient and iteratively estimates the Hessian matrix. Yet, this generic optimization was found to be fairly time consuming and sometimes unstable.

Instead, we propose to adopt a Gauss-Newton optimization approach where we approximate the Hessian matrix by ignoring the term involving second order derivatives. More precisely, the Hessian of the functional is computed as

H = ∇ 2 θ S L S = -∇ θ S d ⊗ (u -µ) -d ⊗ ∇ θ S µ + ∇ 2 θ S log p(θ S |α).
After dropping the first term, we get the following approximate Hessian H ≈ H = -d⊗∇ θ S µ+∇ 2 θ S log p(θ S |α). When inserting the expression of the gradient of the prior, we get :

H = d Diag(µ • (1 -µ)) d T + ∇ 2 θ S log p(θ S |α)
where • is the element-wise product between two vectors. This approximate Hessian matrix is positive definite by construction and is then used to perform several Newtons steps.

The sketch of the MS step is shown as algorithm 1 where the shape parameter prior p(δθ i S |α) is arbitrarily chosen as a zero mean Gaussian distribution with covariance Σ 0 θ S . It consists of two intertwined loops, the innermost performing iteratively the Newton updates and updating the mean, gradient and Hessian values. The outer loop updates the shape function gradient which is potentially a costly step. In line 15 of the algorithm, the U variable is updated in an E-step in order to speed-up the convergence of the overall EM algorithm. Since the parameter range is bounded, we perform in practice a truncated Newton step as proposed in [START_REF] Nash | Newton-type minimization via the lanczos method[END_REF].

This Gauss-Newton approach was inspired by the iterative re-weighted least squares algorithm [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] developed for solving logistic regression (LR) problems. Indeed the first term of L S is similar to the log likelihood of LR after replacing u n with a binary variable and linearizing the shape function. The proposed approach is also related to the Fisher scoring algorithm (see [START_REF] Sourati | Intelligent labeling based on fisher information for medical image segmentation using deep learning[END_REF] as an example in medical image analysis) when the point-wise Hessian matrix of the log likelihood is replaced by its expectation thus leading to more stable evaluation. In this particular case, the approximate Hessian is not the expectation of the Hessian since the first term of L S is the expectation of the log-prior with respect to binary variable U instead of Z.

Finally, the proposed algorithm also outputs a Laplace approximation of the shape parameter posterior p(θ S |I) as a Gaussian distribution where the mean is the optimized shape parameter θ S and the covariance is the inverse approximate Hessian matrix Σ θ S = ( H) -1 .

The overall optimization finally consists in iterating a series of outer loop, each loop consisting in optimizing the shape parameters as in Alg. 1 then followed by a series of MI-steps until the relative change of intensity parameters is less than a threshold. The stopping criterion for the outer loop is the relative change of foreground intensity parameters as it is the most impactful parameter.

Influence of the characteristic length l k

ref Based on Eq.2.5 and Eq.1, it is easy to see that for infinitely small value of the characteristic length l ref → 0, then the label prior becomes more and more sharp p(Z n = 1|θ S ) → δ SDM(S(θ S ),x)>0 and the label posterior becomes equal to the label posterior : p(Z n = 1|θ S , θ I , I n ) -→ p(Z n = 1|θ S ). Conversely, for infinitely large value of the characteristic length l ref -→ ∞, the label prior becomes uninformative p(Z n = 1|θ S ) -→ 1 2 and the label posterior converges towards the appearance driven label posterior :

p(Z n = 1|θ S , θ I , I n ) -→ p(I n |Z n = 1, θ 1 I )/(p(I n |Z n = 0, θ 0 I ) + p(I n |Z n = 1, θ 1 I ))
. Therefore the characteristic length controls the relative influence of the shape and appearance information in the probability of assigning a label.

Since it is scaling the signed distance function, l ref can be interpreted as controlling how far the resulting shape given by p(Z n = e 1 |θ S , θ I , I n ) = 0.5 is allowed to deviate from the reference shape given by S(θ S ). More precisely, assuming a uniform distribution of the appearance label probability between 0 and 1, one can compute the expectation of the posterior probability for a voxel located as a distance d n from the Algorithm 1: MS step to compute p(θ S |I)

1 i ← 0; 2 u n = p(Z n = 1|I, θ S ); // E-Step, Update U 3 repeat 4 V ← S(θ i S ,x n ) l ref ∈ R N ; // Shape function 5 d ← ∇ θ S S(θ i S ,x n ) l ref ∈ R |θ S |×N ; // Shape function gradient 6 δθ 0 S ← 0, t ← 0; 7 repeat 8 µ ← σ(V + d T δθ t S ); // Current Prior probability 9 g ← -d(u -µ) -(Σ 0 θ S ) -1 δθ t S ; // Functional Gradient 10 H ← d Diag(µ • (1 -µ)) d T -(Σ 0 θ S ) -1 ; //
u n = p(Z n = 1|I, θ S ); // E-Step, Update U 16 θ i+1 S ← θ i S + δθ t+1 S , i ← i + 1; // end inner loop 17 until δθ t+1 S / θ t+1 S < ; 18 θ S ← θ i S , Σ θ S = Σ ;
// Gaussian posterior reference shape :

E(p(Z n = 1|θ S , θ I , I n )) = 1 0 tS (∆ n ) tS (∆ n ) + (1 -t)(1 -S (∆ n )) dt = 1 -∆ n e -∆ n -e -∆ n (e -∆ n -1) 2 , ∆ n = d n l ref
Based on the graph of Fig. 1, a voxel located at least at 4l ref inside the boundary of the reference shape S(θ S ) (p < -4) will have in average at least 95% probability to be classified as belonging in the object.

Application to Cochlea Shape Recovery

Cochlea shape model

We use a parametric cochlea shape model which is controlled by a set of 4 deformable shape parameters θ S D : {a, b, α, ϕ} as shown in Fig: (2). Those 4 parameters control the deformation of the centerline of the cochlea represented as a generalized cylinder and is detailed in Appendix B. In addition to those 4 deformable parameters, we consider the 6 pose parameters θ S R consisting of rotation {rx, ry, rz} (parameterizing a rotation vector) and translation {tx, ty, tz} values. Therefore, the total number of shape parameters is 10, controlling the rigid and non rigid (deformable) motion: To fit in our framework, signed distance map SDM(S(θ S ), x) from the cochlea triangular mesh surface must be created. This can be performed for instance by using VTK functions [START_REF] Maurer | A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions[END_REF] but that distance map generation may take several seconds on large volumetric images. This is why we have developed a convolution neural network, noted as DLSDM, which outputs an approximation of the signed distance map from the set of deformable shape parameters in few milliseconds on CPU (Wang et al., 2020a). manner. Many simple generative models such as Gaussian mixture models (GMM) with spatial corrections (Pohl et al., 2006a;Ashburner and Friston, 2005b) have been proposed in the literature to describe tissue intensity distributions. For the cochlea segmentation in CT images, we propose an unsupervised approach based on mixture of mixtures of Student's t-distributions, i.e. each background and foreground regions are described as mixtures of Student's t-distributions. Those t-distributions are generalized Gaussian distributions with heavy tails and lead to more robust estimations than GMM since they are less sensitive to extreme intensity values [START_REF] Peel | Robust mixture modelling using the t distribution[END_REF]. In this context, the probability of observing intensity I n knowing the label Z n is parameterized as :

θ S = θ S D ∪ θ S
I n Z n θ S α θ I β l ref N 2.5e+3 2000
p(I n |Z n = k, θ I ) = M k m=1 π k m t(I n |µ k m , σ k m , ν k m ) , (5)
where M k corresponds to the number of mixture components for the class k and mixture coefficients π k m are positive and sum to one M k m=1 π k m = 1. The mean parameter µ k i , standard deviation coefficient σ k i and degrees of freedom ν k i are parameters of the Student's t-distribution defined as :

t(I n |µ, σ, ν) = Γ ν+1 2 Γ ν 2 1 √ πνσ 1 + (I n -µ) 2 σ 2 ν -( ν+1 2 ) , (6)
where Γ(•) is the gamma function. To write the likelihood of this Student's t-distribution mixture of mixtures, we introduce a new categorical variable τ nkm which is a binary 1of-Mk encoding such that τ nkm = 1 if voxel n belongs to the m-th component of region k, and M k m=1 τ nkm = 1. The likelihood then writes as:

p(I n |Z n , τ n ) = 1 k=0 M k m=1 t(I n |µ k m , σ k m , ν k m ) τ nkm Z nk
The inference is performed with closed-form updates of all parameters [START_REF] Peel | Robust mixture modelling using the t distribution[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] after writing the Student's t-distribution as a Gaussian scale mixture. The total number of parameters to estimate is then |θ I | = 4(M 0 + M 1 ). For the cochlea segmentation problem, we assume that the cochlea region mainly consists of two components (M 1 = 2) : the fluid (perilymph and endolymph) component centered around 0 HU and the bony walls centered around 500 HU. For the cochlea background, we consider 4 components (M 0 = 4) centered around 0 HU (fluid), 2000 HU (bony labyrinth), -1000 HU (air due to pneumatization) and 600 HU (temporal bone). The corresponding initial distribution of intensity in the background and foreground regions are shown in Fig. 4 and the exact initialization values are provided in Appendix C.

Results

Synthetic Images

We provide a 2D synthetic example to illustrate the influence of the reference length l ref in the proposed segmentation algorithm. We consider the segmentation of an ellipse with Gaussian intensity distribution on both background and foreground (see Fig. 5 (Top Left)) by using a circle prior shape S(θ S , x) = x -C 2 -R 2 . It illustrates the frequent case where the parametric model used as a prior is far simpler than the shape visible in the image. The intensity model consists of two Gaussian distributions initialized with mean and variance offsets and the circle is parameterized by its center coordinates and radius. The trade-off between imaging information (leading to an ellipse) and prior shape (leading to a circle) is controlled by the l ref parameter. The log-likelihood as a function of l ref exhibits a single maximum for l ref = l opt = 0.021 (Fig. 5 (Middle)) corresponding to the white circle in Fig. 5 (Left) and to the posterior label distribution in Fig. 5 (Right). The resulting segmentation is the isocontour p(Z n = 1) = 0.5, displayed as a yellow curve in Fig. 5 (Left), which closely matches the elliptic shape except at its flat part (see arrow). This optimal value of l ref corresponds to a configuration where the area of the circle is roughly equal to the area of the ellipse. A value of l ref < l opt leads to isocontours p(Z n = e 1 ) = 0.5 that fit more closely the ellipse whereas l ref > l opt leads to isocontours that look more like a circle.

Inner Ear Datasets

The evaluation of the proposed approach is studied on 3 different datasets.

Dataset #1. includes spiral CT temporal bone images of 210 patients from the radiology department of Nice University Hospital of size 512 × 512 × 178 corresponding to a voxel size of 0.185mm, 0.185mm, 0.25mm. They have then been registered to a reference image via an automatic pyramidal blocking-matching (APBM) algorithm [START_REF] Ourselin | Block Matching : A General Framework to Improve Robustness of Rigid Registration of Medical Images[END_REF] from the software MedInria [START_REF] Toussaint | MedINRIA: Medical Image Navigation and Research Tool by INRIA[END_REF] followed by an image reformatting around the cochlea to the dimension (60, 50, 50) with isotropic voxel size of 0.2mm. The relatively robust registration provides a rough alignment of the cochlea visible in the input image with a cochlea reference frame. From that dataset, 5 CT images were manually segmented by an ENT surgeon (see section 4.4). Dataset #2. includes 9 cadaveric cochlea spiral CT images acquired at the face and neck institute at Nice University Hospital with the same size and voxel spacing as dataset #1. In addition to CT images, high resolution X-ray microtomography (a.k.a. µCT ) images with dimension of (1035, 800, 1095) and isotropic voxel spacing of 0.02479mm were acquired each subject. The 9 µCT and spiral CT images have been registered together as shown in Fig 6 and reformatted around the cochlea to the same physical size as for dataset #1 (i.e. 12mm, 10mm, 10mm). The cochlea and its two scala have been segmented on both CT and µCT images by an ENT surgeon with a semi-interactive tool [START_REF] Criminisi | GeoS: Geodesic Image Segmentation[END_REF]. The high resolution µCT masks serve as ground truth information for the location of the cochlea.

Dataset #3. is a human bony labyrinth dataset [START_REF] Wimmer | Human bony labyrinth dataset: Co-registered ct and microct images, surface models and anatomical landmarks[END_REF] which includes 22 bony labyrinth CT images and their corresponding µCT images having isometric voxel size respectively of 0.1562mm and 0.0607mm. Those images were preprocessed and reformatted as for dataset #2 and also contains manually segmented cochlea masks.

CT μCT 1mm 1mm

Figure 6: A visual comparison of imaging resolution between the µCT and conventional CT for cochlea imaging.

4.3. Quantitative evaluation of segmentation on post-mortem µCT /CT datasets #2 and #3 Baseline Approach: We have implemented a 3D atlas based segmentation approach and applied it on dataset #2 and #3 to get a baseline accuracy in terms of Dice score. To this end, we randomly select one image from each dataset as template image and for each input image we perform a multiscale demons deformable registration Logistic Shape Model Inference: In all cases, the deformable shape parameters θ S D were initialized as (a = 4, b = 0.15, α = 0.6, ϕ = 0.2) and the pose parameters were set to zero. The l ref value was set between 0.1 and 0.3 (see section 4.3) and the stopping condition is ∆θ θ < 0.1, thus stopping when parameter updates are less than 10% of the parameter values.

Computational efficiency: We analyze the computational cost of several alternative formulations of our algorithm. More precisely, in Table 1 we compare the computational time of three different implementations of our approach that differ by the choice of the quasi-Newton optimization method in the MS-step (BFGS vs Gauss-Newton) and by the algorithm used for generating signed distance maps ( VTK based vs deep learning based). The various algorithms was applied on the 9 images of dataset #2 and ran on a Dell Precision 7520 computer. It is clear that the Gauss-Newton method described in Algorithm 1 is far more efficient since it uses a much better approximation of the Hessian matrix than in the generic quasi Newton approach. Furthermore, as expected, the trained deep learning method leads to a speedup factor greater than 3.

Influence of the reference length To assess the influence of the hyper parameter reference length: l re f , we analyse the variation of the final Dice score for various reference lengths based on one image of dataset #2. The results are shown in Table 3. We see that the reference length within the range [0.05mm, 0.25mm] has a relatively small influence on the Dice score. To minimize the time of computation, we do not optimize the reference length through a greedy search but simply set its value to 0.1 for dataset #1 and #2 and 0.3 for dataset #3 for shape fitting. To compute the final hard segmentation we use a fixed reference length of 0.25.

Robustness analysis To study the robustness of the method, we randomly initialize the cochlea shape parameters by performing a random uniform sampling within their defined value range. Based on 10 initial random samples, we computed the average This clearly shows the increased robustness with respect to initial shape values obtained by the Gauss-Newton optimization of the MS-step.

Evaluation on CT and µCT images: Datasets #2 and #3 include both CT and µCT images of the same subject that have been registered to each other. Furthermore the cochlea was manually or semi-automatically segmented by an expert on both modalities such that we can use those two binary maps to evaluate the accuracy of the algorithm applied on the CT image. The cochlea binary map from high resolution µCT images have been downsampled and represent a more reliable ground truth than the manual segmentation performed on the CT images.

The proposed algorithm using Gauss-Newton optimization and deep-learning generation of signed distance maps was applied on the 9+22 CT images of the two datasets. Fig. 3 (Right) shows the segmented cochlea in red, the associated shape model, and its evolution during the MS step. Clearly, we see that the resulting segmentation is strongly constrained by the shape model.

In Table 2, we provide two metrics between pairs of binary masks : the Dice score and the 95% and 100% symmetric Hausdorff distance (HD) (computed as the average of two distances). Furthermore, we compare the segmentations produced by the posterior label probability (SROI for p(Z n |I n , θ S , θ I ) = 0.5) and the ones produced by the shape model only (SSI for p(Z n |θ S ) = 0.5) with both manual segmentations obtained on CT and µCT images. To measure the uncertainty in the manual CT segmentation, we also evaluate the metrics between both CT and µCT manual mask images.

The logistic shape model framework produces good segmentation results on both datasets (Dice scores of 0.81 and 0.91) and even slightly outperforms the manual CT segmentation on dataset#2 (0.81 vs 0.7) which is far more challenging dataset #3. The segmented shape instances produced by the shape model are not as accurate as the SROI for the cochlea segmentation (lower Dice score and larger HD). This confirms that the parametric geometric cochlea model is a simplified representation of the cochlea anatomy. Finally, the metrics between the 2 manual segmentations on dataset #2 (DSC of 0.7 with a 95% HD of 0.5mm) shows the difficulty of performing a manual segmentation of the cochlea due to its limited size and contrast.

Semi-quantitative analysis of segmentation on clinical dataset #1

We ran the segmentation framework (with DLSDM) on the 210 CT of dataset #1 on a Dell 6145 and 6420 CPU clusters.

Unsupervised quality control and semi-quantitative evaluation As manual segmentations of the 210 images are not available, we propose instead an original approach to estimate our algorithm's performance while minimizing the manual annotation effort. First, we apply the unsupervised quality control algorithm of Audelan et al. (Au-delan and Delingette, 2019) on the whole dataset in order to sort the 210 segmentations according to their hypothesized performance. More precisely, this quality control algorithm computes for each image segmentation, an average distance between the segmentation provided by our algorithm and a segmentation produced by a simple generic probabilistic method. We can then generate an histogram of such average surface error (ASE) in Fig. 7. Segmentations having a low ASE correspond to those having good intensity contrast across their boundaries while those on the right tail of the distributions are considered as more challenging and suspect of including segmentation errors.

The histogram exhibits a bell shape with few outliers on its right and left tails. Furthermore, we have manually checked that this unsupervised quality control algorithms worked well on this dataset with visually better segmentations localized on the left side of the histogram. To estimate the relation between the ASE and the Dice score, we picked 5 images in order to sample the histogram at different levels of ASE corresponding to images #213, #210, #53, #264 and #143 (see Fig. 7) in ascending order of ASE. Those 5 CT images were manually segmented by an ENT surgeon and the Dice scores of the segmentation produced by our algorithm was reported in Table 4. We see that the Dice score decreases as the ASE increases which indicates that the ASE may be a proper surrogate for the segmentation performance. The cochlea in image #143 was indeed found to be an outlier in terms of shape probably due to a patient malformation. Inspired by [START_REF] Audelan | Unsupervised Quality Control of Image Segmentation based on Bayesian Learning[END_REF], we can make the hypothesis that ASE a good proxy for the Dice score as there is a monotonic relation between ASE and Dice. On this basis, we can extrapolate that the median Dice score over the whole dataset is probably above 0.82. Yet, a more thorough study with far more manual segmentations is necessary to be less speculative about the actual performance on clinical CT data. Parameter analysis The application of the algorithm on dataset #1 resulted in the estimation of 10*210 shape parameters with 210 covariance matrices Σ θ S . In Fig. 8(a) the histograms of the 4 deformable shape parameters are displayed in green. Interestingly, the a and α parameters exhibit a bimodal distribution for which a simple explanation may be provided. Indeed, the left highest mode is probably correspond- ing to straight centerline profiles whereas the rightmost mode may be associated with the "rollercoaster" longitudinal profiles [START_REF] Avci | Variations in microanatomy of the human cochlea[END_REF]. In Fig. 8(b) the average 10 × 10 covariance computed as the log-Euclidean mean [START_REF] Arsigny | Geometric means in a novel vector space structure on symmetric positive-definite matrices[END_REF] is displayed showing potential correlations between pose and deformable parameters. Shape parameter b corresponding the exponent of the logarithmic center line curve has a particularly low variance such that it can be well estimated from the data. Conversely the phase parameter ϕ has much greater variance and is harder to estimate in average. The extraction of the eigenvectors of that covariance matrix confirms that most parameters are independent with other except for parameter a which is a bit correlated with the ϕ parameter and the translation term ty. The relative independence of the parameters for shape fitting indicates that we do not have an overparametrization of the cochla shape. Uncertainty Segmentation analysis The estimated covariance matrix Σ θ S can be used for studying the uncertainty of the output segmentation. We sampled 100 times the multivariate Gaussian approximate posterior distribution of the parameters p(θ S |I) ≈ N(θ S ; Σ ) and generated accordingly 100 random posterior labels p(Z|I, θ S , θ I ) that are then averaged to estimate with Monte-Carlo sampling the marginal posterior p(Z|I, θ

I ) = R |θ S | p(Z|I, θ S , θ I )p(θ S |I) dθ S .
In Fig 9 we show several slices of the resulting probability maps with the 0.5 level curve together with the posterior probability p(Z|I, θ S , θ I ) obtained with the most likely shape parameter θ S . We see that we have a much larger uncertainty in the resulting segmentation when accounting for the uncertainty in the shape parameters than without them. This is a far better approximation of the true uncertainty p(Z|I) than the posterior label probability p(Z|I, θ S , θ I ).

Comparison with the state-of-the-art

We consider below the prior work on cochlea segmentation evaluated on clinical CT images while discarding the literature on the segmentation of µCT images [START_REF] Kjer | Semiautomatic anatomical measurements on microCT 3D surface models[END_REF][START_REF] Ruiz Pujadas | Random walks with shape prior for cochlea segmentation in ex vivo µCT[END_REF][START_REF] Pujadas | Cochlea segmentation using iterated random walks with shape prior[END_REF] or of the scala tympani and vestibuli located inside the cochlea [START_REF] Noble | Statistical shape model segmentation and frequency mapping of cochlear implant stimulation targets in CT[END_REF][START_REF] Noble | Image-Guidance enables new methods for customizing cochlear implant stimulation strategies[END_REF]. Table 5 summarises the relevant publications on cochlea segmentation that are split into unsupervised and supervised methods. The former approaches are mostly based on cochlear shape fitting based on template image registration [START_REF] Baker | Model-image registration of parametric shape models: fitting a shell to the cochlea[END_REF], parametric shape model [START_REF] Baker | Tracking, modelling and registration of anatomical objects: the human cochlea[END_REF]. The supervised methods are based on statistical deformation models (Ruiz [START_REF] Ruiz Pujadas | Random walks with statistical shape prior for cochlea and inner ear segmentation in micro-ct images[END_REF] and deep learning [START_REF] Lv | Automatic segmentation of temporal bone structures from clinical conventional ct using a cnn approach[END_REF][START_REF] Raabid | Automatic segmentation of inner ear on ct-scan using auto-context convolutional neural network[END_REF][START_REF] Heutink | Multi-scale deep learning framework for cochlea localization, segmentation and analysis on clinical ultra-highresolution ct images[END_REF].

Quantitative comparison of performances is not straightforward due to differences in image modality (CT, µCT or ultra high resolution CT), in metrics (Dice, precision, mean surface error), in subject population (cadaveric vs patient) but also in the target anatomical structures (cochlea vs cochlea labyrinth). In most cases, cochlea segmentation from µCT images are used as ground truth information and a direct comparison between our work with [START_REF] Raabid | Automatic segmentation of inner ear on ct-scan using auto-context convolutional neural network[END_REF] is possible since they used a subset of dataset #3 which is a public database [START_REF] Wimmer | Human bony labyrinth dataset: Co-registered ct and microct images, surface models and anatomical landmarks[END_REF]. We see that our unsupervised approach performs as well as the supervised methods with Dice scores in the range [0.85, 0.91] and outperforms previous unsupervised methods.

Discussion

The proposed approach relies on the definition of a generic shape function S(θ S , x) which can be for instance a statistical shape model, a deformation image template, or an implicit shape equation. In the case of the cochlea, it was defined as a signed distance function of a parametric shape model SDM(S(θ S ), x). This specific choice makes the computation of the shape function and its gradient fairly costly, despite the use of a fully supervised dedicated neural network (DLSDM). There are several ways to optimize its computation time. One could for instance use a supervised appearance model such as a trained neural network which would remove all MI steps in the EM algorithm and would decrease by at least a factor 2 the time of computation. Another way is to use an implicit shape model S(θ S , x) = 0 for instance based on statistical level sets [START_REF] Tsai | A shape-based approach to the segmentation of medical imagery using level sets[END_REF]). The cochlea segmentation example provided in this paper relies on a fully interpretable intensity and shape parameters at the expense of its computational efficiency. Yet, one could train a deep neural regressor for predicting cochlea shape parameters and segmentation by using the segmentations generated by the proposed framework as training set.

For the cochlea segmentation, excellent results were obtained on cadaveric CT images similarly to the supervised methods. Furthermore, to the best of our knowledge, we introduced a first semi-quantitative assessment of cochlea segmentations on clinical CT images acquired on more than 200 patients. However, for a complete study, one would need to assess thoroughly the inter-rater variability of those manual segmentations and ideally combine them with other high resolution image modalities. Finally, an interesting extension of this work would be to segment the scala vestibuli and tympani in addition to the cochlea.

Conclusion

In this paper, we have presented a new probabilistic generative approach for combining shape and intensity models for image segmentation. The resulting segmentation is an interpretable compromise between a fidelity to a parametric shape space (captured in the prior distribution) and an appearance model (captured in the likelihood distribution). The proposed method goes well beyond the concept of shape fitting since it also provides an approximation to the posterior distribution of shape parameters. The use of a logistic shape model allows to control the trade-off between appearance and shape with a single parameter: the reference length. When applied to the recovery of cochlea structures from CT images, we were able to provide accurate segmentations with meaningful shape parameter distributions. Furthermore, we have shown how the approximate shape parameter posterior distribution can be exploited to provide realistic uncertainty maps.

An interesting application of the proposed approach is to perform model selection with Bayes factors, in order to estimate the optimal complexity of a parametric shape model for a given image segmentation task. Future work will also explore the application of this framework to other shape representations than explicit parametric shape models in order to find a reasonable trade-off between computational efficiency and interpretability of shape parameters. For instance, in statistical deformation models (Rueckert et al., 2003a), the computation of shape function gradient ∇ θ S S(θ S , x) is straightforward, but its shape parameters may not be meaningful besides the first modes.

The longitudinal component of the centerline is the sum of an exponentially damped sinusoidal and a linear function:

z(θ c ) =        αe -βθ c cos(θ c + φ) + q 1 θ c if θ c < θ 1 a 2 θ 2 c + a 1 θ c + a 0 if θ c ≥ θ 1 , (B.3)
where β = 0.2 rad -1 , q 1 = 0.225 mm.rad -1 and θ 1 = θ max -π. The polynomial function is used to flatten out the last half turn so that dz(θ)/dθ| θ=θ max = 0 and similarly a 2 , a 1 , a 0 are set to obtain a continuously differentiable curve.

Cross-Sections The cross-sections are modeled by a closed planar shape on which a varying affine transformation is applied along the centerline. The scala tympani and the scala vestibuli are modeled with two half pseudo-cardioids while the cochlear crosssection corresponds to the minimal circumscribed ellipse of the union of the tympanic and vestibular cross-sections. The affine transform of cross-sections is parameterized by a rotation, and a width and height scalings. All cross-sectional parameters are fixed because their variability was found to be small compared to the variability of the centerline.

Shape parameter vector We have chosen a compact description of the cochlea shape to limit as much as possible the correlation between the shape parameters and therefore make them uniquely identifiable. Finally, only 10 free parameters are considered in θ S :

• 6 translation and rotation parameters : t = (tx, ty, tz), r = (rx, ry, rz)

• 2 radial component parameters of the centerline, a and b

• 2 longitudinal component parameters, α and φ Note that there are no free cross-section parameters which implies that θ S can be used to define uniquely the cochlea.

The prior probabilities on the 10 shape parameters were modeled as being an uniform distribution (uninformative prior) such that all regularization terms log p(θ S |α) can be ignored.

Figure 2 :

 2 Figure 1: Expected label posterior probability as function of the normalized signed distance from the reference shape.

Figure 3 :

 3 Figure 3: (a) graphical model for the shape-based generative model; (b) Cochlea segmentation on CT images is shown in solid red with the associated shape model in dashed yellow lines; (c) Evolution of the cochlea shape model during several MS steps shown as 2D contours (from dotted green to solid red) and 3D models.

3. 2 .-Figure 4 :

 24 Figure 4: Example of intensity probability distributions of the foreground ( f g , in red) and the background (b g , in blue) as functions of the Hounsfield unit.

  Figure 5: (Left) Input Ellipse image fitted with a circle shape : initial circle (red), final circle (white) and 0.5 isocontour of posterior label probability for optimal value of l ref (yellow);(Middle) Log likelihood as a function of l ref ; (Right) posterior label probability p(Z n = 1|θ S , θ I ) for optimal value of l ref ;

  one image of dataset #2 and obtained a mean Dice score of 0.81 ±0.1 ( respectively of 0.68 ±0.23 ) for the Gauss-Newton method (resp. the BFGS method).

  : Average surface error for dataset #1

Figure 7 :

 7 Figure7: Average surface error of segmentations generated from dataset #1 resulting from the unsupervised quality control. Red contours correspond to the manual ground truth while yellow ones are segmentation outputs.

Figure 8 :

 8 Figure 8: Distribution plots for shape parameters variance. (b) Average covariance matrix of the 10 shape parameters.

Figure 9 :

 9 Figure 9: Marginal posterior probability p(Z|I, θ I ) (Top) versus posterior probability p(Z|I, θ I , θ S ) (Bottom) computed on patient #1 of dataset #1.

Table 1 :

 1 Computational efficiency proposed methods

	BFGS	VTK SDM DLSDM
	Mean Comput. Time 12h15min	43min	16min

Table 2 :

 2 Performance metrics obtained on dataset #2 and #3.

					Symmetric Hausdorff Distance
	Compared Labels	Dice Score	(voxel size 0.2 mm) Dataset #2 Dataset #3
			Dataset #2	Dataset #3	95% 100% 95% 100%
	CT Manual	SSI SROI	0.74 ± 0.02 0.77 ±0.023 0.53 0.85 ± 0.011 0.91 ±0.015 0.34	1.04 0.82	0.70 0.36	1.91 1.68
		SSI	0.67 ± 0.024 0.76 ±0.068 0.68	1.48	0.67	1.96
	µCT Manual	SROI	0.81 ± 0.04 0.91± 0.019 0.50	1.31	0.36	1.68
		CT Manual 0.70 ± 0.084 0.93± 0.021 0.50	1.34	0.19	0.74

[START_REF] Vercauteren | Diffeomorphic Demons Using ITK's Finite Difference Solver Hierarchy[END_REF]

(as implemented in SimpleITK 1.0.1) to estimate the deformation field. The segmented mask of the template is deformed to match the target image. The average Dice scores are 0.63 for dataset #2 and 0.68 for dataset #3.

Table 3 :

 3 Influence of the hyper parameter: l re f for the segmentation accuracy.

Table 4 :

 4 Dice score for selected segmentation samples from dataset #1 based on the histogram of Fig.7. The ASE are got from automatic quality control algorithm and the DICE score are computed based on manual segmentation.

	Patient ID	213 210 53	164 143
	DICE Score 0.84 0.84 0.84 0.76 0.45
	ASE	0.21 0.30 0.32 0.36 0.50

Table 5 :

 5 Performances of prior work on cochlea segmentation. NL (resp. NT) indicates the number of training (resp. testing) images. Unsup (resp.

	Sup) refers to

Appendix A. Gradient of shape function

In this section, we detail the computation of the shape function gradient ∇ θ S S(θ S , x) when rigid and deformable shape parameters are considered. More precisely, writing the parameters controlling the non-rigid deformation as θ S D , the shape function writes as S(θ S D , Rx n +t). The rotation matrix R is parameterized with rotation vector r, whose norm is the rotation angle and whose direction is the rotation axis. The gradients with respect to the translation and rotation vectors are then given in closed form as :

where ∇ x is the spatial gradient, S x is the 3x3 anti-symmetric matrix associated with vector x. For a deformable parameter θ S D , if the shape function is not given in an analytical form as it is the case for parametric shapes, the shape function gradient can be computed with finite differences based on a parameter increment δθ i S D :

We are interested in the cochlea structure in CT images which is defined as a generalized cylinder, i.e. as cross-sections swept along a centerline.

Centerline The centerline is parameterized in a cylindrical coordinate system by its radial r(θ c ) and longitudinal z(θ c ) components. The range of polar angle θ c is [0, θ max ] where θ max is the maximum polar angle controlling the total number of cochlear turns.

The radial component is defined piecewise with a polynomial function and a logarithmic function of the angular coordinate θ c in the cylindrical coordinate system as:

where θ 0 = 5π/6 and p 0 = 5 mm. Furthermore to obtain a continously differentiable curve, we set :

Appendix C. Initialization of intensity parameters

The 4 * 6 = 24 initial intensity parameters for the mixture of Student's t distributions in datasets #1 and #3 are presented in Table ??.


Brooke Stephenson
email: brooke.stephenson@gipsa-lab.grenoble-inp.fr

Thomas Hueber
email: thomas.hueber@gipsa-lab.grenoble-inp.fr

Laurent Girin
email: laurent.girin@gipsa-lab.grenoble-inp.fr

Laurent Besacier
email: laurent.besacier@univ-grenoble-alpes.fr

Alternate Endings: Improving Prosody for Incremental Neural TTS with Predicted Future Text Input

Keywords: Incremental text-to-speech, prosody, neural language models

Inferring the prosody of a word in text-to-speech synthesis requires information about its surrounding context. In incremental text-to-speech synthesis, where the synthesizer produces an output before it has access to the complete input, the full context is often unknown which can result in a loss of naturalness. In this paper, we investigate whether the use of predicted future text from a transformer language model can attenuate this loss in a neural TTS system. We compare several test conditions of next future word: (a) unknown (zero-word), (b) language model predicted, (c) randomly predicted and (d) ground-truth. We measure the prosodic features (pitch, energy and duration) and find that predicted text provides significant improvements over a zero-word lookahead, but only slight gains over randomword lookahead. We confirm these results with a perceptive test.

Introduction

In incremental text-to-speech synthesis (iTTS), the system starts to output chunks of synthetic audio before the full text input is known [START_REF] Baumann | Evaluating prosodic processing for incremental speech synthesis[END_REF][START_REF] Baumann | Partial representations improve the prosody of incremental speech synthesis[END_REF][START_REF] Pouget | HMM training strategy for incremental speech synthesis[END_REF][START_REF] Pouget | Adaptive latency for part-of-speech tagging in incremental text-to-speech synthesis[END_REF]. The missing input information often hinders the ability to produce a natural sounding speech sequence, mostly because prosodic features that will be determined by the future context (i.e. the remaining words in the sentence ) have not yet been specified. Fortunately, the future input is not completely random; human language is characterized by several lexical and syntactic patterns, which can be statistically learnt and then predicted to a certain extent. Recent advances in language modelling, namely the use of transformer models such as BERT [START_REF] Devlin | BERT: Pretraining of deep bidirectional transformers for language understanding[END_REF] and GPT-2 [START_REF] Radford | Language models are unsupervised multitask learners[END_REF] give us accurate representations of the probability distribution of future words. If this information can be mobilized to fill in the missing data for an iTTS system, it may be possible to retain naturalness while minimizing latency.

Early work in iTTS was conducted in the context of HMMbased models, where linguistic and phonological features extracted from text were used to estimate speech parameters. [START_REF] Baumann | Decision tree usage for incremental parametric speech synthesis[END_REF] studied the effects of missing future features by replacing decision tree split criteria with default values and evaluating the degradation in speech quality: while cepstral and aperiodicity features could be estimated fairly accurately with just a local context, prosodic features (f0 and duration) were found to be more dependent on longer range context. Using a similar default value strategy, [START_REF] Baumann | Partial representations improve the prosody of incremental speech synthesis[END_REF] studied symbolic intonation assignment in the presence/absence of word, phrase and utterance level features. They report that phrase and utterance final words benefit the most from the inclusion of phrase and utterance level features in prosodic assignment determination. [START_REF] Pouget | HMM training strategy for incremental speech synthesis[END_REF] explicitly specified unknown features in the context clustering process and found improvements over a default value strategy. [START_REF] Baumann | Evaluating prosodic processing for incremental speech synthesis[END_REF] tested just-in-time strategies for integrating future chunks in a dialogue system and concluded that incorporating the next phrase once the first word of the current phrase had been processed gave the best latency/quality trade-off. Finally [START_REF] Pouget | Adaptive latency for part-of-speech tagging in incremental text-to-speech synthesis[END_REF] developed an adaptive policy which delayed synthesis when there was high uncertainty regarding POS tags.

Recent research in iTTS has focused on end-to-end neural models. While these models create more natural speech, they are also more difficult to analyze because the relevant features for the task are learnt during training and are subsequently not easily human interpretable. In this new paradigm, the tradeoff between speech quality and synthesis latency has been examined by testing the effects of different degrees of lookahead [START_REF] Stephenson | What the future brings: Investigating the impact of lookahead for incremental neural TTS[END_REF][START_REF] Ma | Incremental text-to-speech synthesis with prefixto-prefix framework[END_REF], reinforcement learning has been used to automatically learn a wait/synthesize policy [START_REF] Mohan | Incremental text-to-speech for neural sequence-to-sequence models using reinforcement learning[END_REF], the effects of synthesis unit have been studied [START_REF] Yanagita | Neural iTTS: Toward synthesizing speech in real-time with end-to-end neural text-tospeech framework[END_REF] and the integration of iTTS into a speech-to-speech translation system [START_REF] Sudoh | Simultaneous speech-to-speech translation system with neural incremental asr, mt, and tts[END_REF] and into a machine speech chain [START_REF] Novitasari | Incremental machine speech chain towards enabling listening while speaking in real-time[END_REF] have been evaluated.

Furthermore, previous studies in conventional (i.e. non incremental) TTS have incorporated language model representations into neural TTS training and found that they could help speed up training time [START_REF] Fang | Towards transfer learning for end-to-end speech synthesis from deep pre-trained language models[END_REF][START_REF] Chung | Semi-supervised training for improving data efficiency in endto-end speech synthesis[END_REF]. In the field of simultaneous translation, where future context is also unknown, [START_REF] Zheng | Speculative beam search for simultaneous translation[END_REF] hallucinated future words to balance the latency/quality trade-off.

In the present work, we propose an iTTS system that incorporates a language model to predict future lookahead. 1 Our approach (described in Figure 1) predicts one word into the future. This limited lookahead was chosen so that the effects of correct and incorrect predictions could be studied. We evaluate our system by contrasting different future word contexts: (a) unknown, (b) language model predicted, (c) randomly predicted (a control group) and (d) ground-truth (see Table 1). Differences are measured at the TTS encoder level and from the generated speech signal through a listening test.

Method

Definitions

For each token in our corpus, we prepare different sequences which are used as input to the TTS model, FastSpeech 2 [START_REF] Ren | Fastspeech 2: Fast and high-quality end-to-end text-to-speech[END_REF].

• x1:n = x1, x2, ..., xn is the sequence of tokens up to n. In the Figure 1: Utilizing language model predictions to improve incremental TTS quality while keeping limited latency.

Table 1: Examples of input sequences with unknown, ground-truth, predicted and random future context. In each sequence, the word in bold is the word which is synthesized from the sequence. proposed iTTS system, the tokenization policy is to split the sentence on space characters, and then synthesis is triggered when a space character is encountered. • k is the lookahead parameter (number of future tokens available when synthesizing token xn).

• s 1:n+k = {x1, x2, ..., xn, xn+1, ..., xn+k } = {x1:n, xn+1:n+k } is
the sequence used for the synthesis of token xn, where for the ground-truth condition (GT) xn+1:n+k = x n+1:n+k , for the prediction condition (Pred) xn+1:n+k is given by the language model, and for the random condition (Rand) xn+1:n+k is random. The random token generation is described in Section 3.1.

• s1:n is the input prompt used to generate language model predictions. • Near the end of the sequence, we replace n+k with min(n+ k, N ) where N is the length of the full utterance.

Models

Language model used for prediction. We use the GPT-2 language model for our study. This is an auto-regressive model trained to predict the next word given a sequence of past words (causal language modeling task), based on a Transformer architecture. The original GPT-2 [START_REF] Radford | Language models are unsupervised multitask learners[END_REF] is large (1.5B parameters) and since our intended use requires fast predictions, we opted to use a smaller version of GPT-2, called "distilled GPT-2" [START_REF] Wolf | Transformers: State-of-the-art natural language processing[END_REF]. 2This model has been trained to produce the same output probability distribution as the original GPT-2 but using fewer layers/parameters. TTS model. For TTS we select a fast and high-quality endto-end model: FastSpeech 2. The implementation we use [START_REF] Hayashi | Espnet-TTS: Unified, reproducible, and integratable open source end-to-end textto-speech toolkit[END_REF], 3trained on the LJ Speech Dataset [START_REF] Ito | The LJ speech dataset[END_REF], takes characters as input and converts them to phonemes. Phoneme embeddings are passed through several self-attention layers before the model makes duration, pitch and energy predictions for each phoneme. These feature predictions and the latent phoneme representations are then passed to the decoder (more self-attention layers) which produces a Mel-spectrogram. 4 The Mel-spectrogram is then input into a Parallel WaveGAN vocoder [START_REF] Yamamoto | Parallel WaveGAN: A fast waveform generation model based on generative adversarial networks with multi-resolution spectrogram[END_REF] (trained on full sentence inputs) for waveform generation. This model is well suited to iTTS because (1) it is fast which is desirable when the objective is to reduce latency (the speed is achieved by predicting all Mel-spectrogram frames in parallel), and (2) it makes explicit duration predictions for each phoneme, which makes it possible to segment words and only synthesize the word(s) of interest.

Incremental synthesis (iTTS)

We implement an incremental synthesis procedure where each token xn is synthesized from the input sequence s 1:n+k . Melspectrogram frames corresponding to individual tokens are identified using the internal duration predictions made by Fast-Speech 2. Successive word-level Mel-spectrograms are input into the Parallel WaveGAN vocoder on a word-by-word basis. Resulting waveforms are concatenated together using a 1ms crossfade to eliminate glitches (synthetic audio samples are available at https://tinyurl.com/ae4nzs).

Experiments

Corpus and predictions

The English corpus we use for analysis consists of 1,000 sentences from LibriTTS [START_REF] Zen | LibriTTS: A corpus derived from LibriSpeech for text-to-speech[END_REF]. Sentence length ranges from 5 to 42 words, with a total of 16,965 tokens and 62,556 phonemes.

For each token xn in the corpus, we sampled five GPT-2 and five random next word predictions (xn+1). The GPT-2 predictions are constrained to the 30 most likely next words (top-30 sampling strategy). The random words were selected from a list of 1,266 of the most common words in English [START_REF] Speer | Luminosoinsight/wordfreq: v2.2[END_REF]. Importantly, we force GPT-2 predictions and random predictions to have comparable lengths in term of characters/phonemes because (1) GPT-2 tends to predict shorter words because they are more frequent, (2) in our previous study [START_REF] Stephenson | What the future brings: Investigating the impact of lookahead for incremental neural TTS[END_REF], we found that longer future words have more influence on the current token's internal representation (in a seq-to-seq model) than shorter ones, (3) otherwise, our results may be biased by the fact that the random condition simply has more future context. To control for word length in the random condition, we (1) took the word-length distribution of GPT-2 predictions, (2) randomly sampled a word-length category from this distribution (e.g. 2-4 characters), (3) limited our most-common list to only words in this category and (4) randomly sampled a word from this list using a uniform distribution.

GPT-2 uses byte pair encoding (BPE) which breaks words down into subword units to better handle out-of-vocabulary tokens. As such, some of its predictions extend the final prompt word rather than predicting a new token (e.g. previous → previously). To avoid such distortions to our input text, we sample until the first character in the predicted text is a space. This also prevents erroneous punctuation marks from being predicted.

Metrics

FastSpeech 2 representations. We aim at evaluating the prosody obtained in the different test conditions: no context (k = 0), ground-truth context (GT), predicted context (Pred), random context (Rand). For this aim, we compare the pitch, duration and energy values produced in those conditions with the values produced in the reference condition (Ref) where the full context (full sentence input) is used. In the present paper, we concentrate on the case k = 1 (one-word lookahead).

As for duration and energy, they are first computed at the phoneme level, using the FastSpeech 2 internal predictions (see Figure 2 for a plot of duration values from an example sentence). A phoneme duration is defined as (the log of) the number of Mel-spectrogram frames of that phoneme. The energy is the squared magnitude of the short-time Fourier transform (STFT), averaged across all frequency bins and across the duration of the phoneme. Then the mean absolute error (MAE) is computed by averaging the absolute value of the difference of duration values obtained in each test condition and in the reference condition across all phonemes of the dataset, and the same for the energy feature. The results are reported in Table 2.

Pitch is evaluated at the sentence level. 5 We first align the Mel-spectrograms obtained in the test and reference conditions with Dynamic Time Warping using the Librosa library [START_REF] Mcfee | Librosa: Audio and music signal analysis in Python[END_REF]. Then we extract the pitch curves from the concatenated audio (see Section 2.3) using Praat/Parselmouth [START_REF] Boersma | Praat: doing phonetics by computer [Computer software[END_REF][START_REF] Jadoul | Introducing Parselmouth: A Python interface to Praat[END_REF] and we compute the MAE in cents between the aligned f0 trajectories:

M AE = 1200/T T t=1 log 2 f T est 0 (t)/f Ref 0 (t) . (1)
Then the sentence-level MAEs are averaged across all sentences of the dataset. The results are reported in Table 3. Perceptive test. Finally, we evaluate the global quality using 40 native English speaking evaluators 6 and a MUSHRA test [START_REF] Itu-R | Recommandation[END_REF]. We selected 20 sentences from our corpus and for each sentence, we presented the listeners with a reference audio clip (generated with the full sentence context) and then asked them to assign a similarity score to five test clips: the hidden reference (identical to the reference and used as the MUSHRA high anchor), k = 0 (used as the low anchor), Ground-Truth k = 1, GPT-2 prediction k = 1 and random prediction k = 1. We then compare the distributions of the similarity scores. The responses from four of the participants were removed because these listeners consistently failed to assign a high similarity score to the high anchor. See Figure 3 for results. 5 We did not evaluate error in the internal FastSpeech 2 pitch predictions because we observed a few extreme prediction values which did not materialize in the resultant audio. 6 Anonymous participants were recruited using Prolific (www.prolific.co). They were compensated at a rate slightly above the UK minimum wage.

Discussion

For all metrics, with regards to the mean, we see a clear ranking in the similarity to the full sentence reference: k = 0 is farthest away, GT k = 1 is the closest and Pred and Rand are in between, the former being slightly closer to full context than the latter. Statistical tests (t-test for pitch, duration and energy measures and Wilcoxon for the listening test) confirmed that Pred and Rand do not belong to the same distribution (p-value < 0.05) and that Pred is better by a small but significant margin. We notice that duration predictions for k = 0 are almost always longer than the other conditions (Figure 2). And as in [START_REF] Baumann | Evaluating prosodic processing for incremental speech synthesis[END_REF], we observe pitch drops for k = 0 words. This is because all words are interpreted as the end of a sentence (as they are the final word in the FastSpeech 2 input, hence sentence final characteristics are predicted by the model). Both the prediction and the random conditions reduce this effect thanks to the additional padding words.

Correct vs. Incorrect Predictions. When we separate the correct from the incorrect GPT-2 next word predictions (see Table 2), we see that the MAE for the incorrect predictions is almost identical to the MAE for the random condition. This suggests that the improved syntactical accuracy gained from the GPT-2 predictions (the POS of the predicted token matches that of the GT next token 43.5% of the time vs. 18.0% for random) does not translate into improved prosodic features.

Since we only see improvement when the exact next word is predicted, it is clear that the minor difference between GPT-2 and random is explainable by the low exact-word prediction rate. We observe that 76% of the GPT-2 sequences have a prediction rate lower than 10%, and 97% have a rate lower than 21% (mean: GPT-2 = 6.8%; random = 0.09%). It is likely that as language models continue to improve [START_REF] Brown | Language models are fewshot learners[END_REF], we will see greater gains in naturalness from the proposed method (improvements in semantic modelling will narrow the range of word choice, resulting in more frequent exact word predictions). However, these potential advances will have a fairly low ceiling if we consider human prediction abilities as the upper limit: [START_REF] Luke | Limits on lexical prediction during reading[END_REF] shows that approximately 5% of context words and 20% of function words are highly predictable by humans. Further prediction gains could be achieved if the language model was fine-tuned on the traits of a specific author [START_REF] Delasalles | Learning Dynamic Author Representations with Temporal Language Models[END_REF]; this would be an advisable step in the use case of assistive technologies for the speech impaired.

Context Sensitivity. Previous studies investigating the impact of lookahead have shown the contrast between k = 0 context and different degrees of ground-truth lookahead. The setup of the present study allows us to investigate where the choice of future context modifies the output the most (i.e. where do prosodic features remain stable irrespective of the future context and where do they vary dependent on the future context). To this purpose, we calculated the range of phoneme duration and pitch feature values predicted by the TTS model in all test conditions except k = 0. More precisely, from the 12 predicted and ground truth conditions (5 × Pred, 5 × Rand, GT k = 1 and full context), we take the max and min values from this set and calculate the difference. This analysis shows that a large portion of phonemes in the corpus alter only slightly when provided with different next word contexts. The pitch range does not exceed 300 cents for approximately 75% of our samples, which falls below the Just Noticeable Difference (JND) threshold for pitch distance found by [START_REF] Hart | Differential sensitivity to pitch distance, particularly in speech[END_REF]. The duration range is limited to a single spectrogram frame (11.75ms) for 40% of phonemes, which, depending on the length of the phoneme, may be imperceptible to the average listener ([33] found a JND of 5%).

We do however see some wide range values in the corpus which explain the large standard deviations in Table 2 and the significant variability of the Pred and Rand scores in the MUSHRA test (Figure 4: the maximum deviation values in a sentence show strong correlation with the mean MUSHRA similarity scores). By examining the corpus, we notice discernible patterns in the locations of large context sensitivity. With respect to pitch, we see large variation when there is a mismatch between predicted and ground truth punctuation at the end of the next word or when there is a reporting verb (e.g. said, exclaimed) rather than the beginning of a new sentence following a punctuation mark. With respect to duration, the largest variance occurs at the beginning of sentences, at punctuation marks and in function words, especially in the coordinating conjunction and.

Conclusion and Perspectives

The results from all metrics show that the language model predicted text does improve prosody when compared to the k = 0 condition. Slight improvements over the random text condition are also observed. We have seen that language model predictions are often incorrect and context mismatches can occasionally cause major distortions compared to the full context prosody. To improve our model, we could a) implement a wait policy that delays synthesis when a context sensitive word is encountered (similar to [START_REF] Pouget | Adaptive latency for part-of-speech tagging in incremental text-to-speech synthesis[END_REF]) or b) retrain the model on both ground truth and predicted input. With this training regime, the iTTS model would find the optimal solution given contextual ambiguity (similar to [START_REF] Pouget | HMM training strategy for incremental speech synthesis[END_REF]); the language model predictions, which frequently differ from the ground truth in terms of relevant contextual features (word choice, next phonemes, number of syllables, POS, and position in the prosodic phrase/utterance) would serve to make the model more robust and perhaps provide more neutral solutions for context sensitive words.

  sentence, k = N -n Do you think that you could manage, Tidy? Unknown (future) k = 0 word s GT 1:n+0 = Do, Do you, Do you think, ... Ground Truth k = 1 word s GT 1:n+1 = Do you, Do you think, Do you think that, ... GPT-2 prediction k = 1 word s Pred 1:n+1 = Do they, Do you agree, Do you think this, ... Random k = 1 word s Rand 1:n+1 = Do dance, Do you until, Do you think art, ...

Figure 2 :

 2 Figure 2: Duration prediction from FastSpeech2 (in number of Mel-spectrogram frames on a log scale) for each phoneme in the sentence "Are there bears up here?" and for the different tested prediction conditions.

Figure 3 :

 3 Figure 3: Violin plots of the distribution of similarity scores between signals generated with full context and signals generated with limited context for the 20 sentences in the MUSHRA test. The middle bars show the mean value.

Figure 4 :

 4 Figure 4: Each point represents a sentence (synthesized under the Pred or Rand condition) from the MUSHRA test. The x-axis shows the scaled and combined (pitch, energy, duration) maximum deviation values (deviation from the full context value) for the phonemes in the sentence. The y-axis shows the mean similarity score for (Pred,Rand) sentences to their full context counterpart, given by the MUSHRA participants. The Pearson correlation coefficient is equal to -0.53.

Table 2 :

 2 MAE (and standard deviation across phonemes) between duration (resp. energy) obtained with full context and with limited context. *unit = number of Mel-spectrogram frames on a log scale; **arbitrary unit: signal is digital, normalized and averaged.

	Input type # phonemes Duration*	Energy**
	k = 0	62,556	0.262 ± 0.297	0.301 ± 0.364
	GT k = 1	62,556	0.077 ± 0.133	0.176 ± 0.241
	Pred k = 1	5×62,556	0.135 ± 0.198	0.247 ± 0.296
	Rand k = 1	5×62,556	0.147 ± 0.208	0.260 ± 0.304
	Correct pred.	38,274	0.086 ± 0.132	0.187 ± 0.239
	Incorrect pred.	274,506	0.142 ± 0.205	0.255 ± 0.301

Table 3 :

 3 MAE between the pitch curves obtained with the full context and with limited context.

	Input type	# sentences Pitch MAE (Cents)
	k = 0	1,000	203.56 ± 45.50
	GT k = 1	1,000	88.57 ± 26.33
	Pred k = 1	5×1,000	120.03 ± 29.34
	Rand k = 1 5×1,000	123.03 ± 30.27

A similar idea was proposed in a contemporary preprint paper[START_REF] Saeki | Incremental text-tospeech synthesis using pseudo lookahead with large pretrained language model[END_REF]. In their work, a context encoder architecture is used.

https://huggingface.co/distilgpt2

https://github.com/espnet/espnet

For implementation details, see https://tinyurl.com/s7p38hcr

Acknowledgements

This work was funded by the Multidisciplinary Institute in Artificial Intelligence MIAI@Grenoble-Alpes (ANR-19-P3IA-0003).