Consistent Atlas Estimation on BME Template Model: Applications to 3D Biomedical Images. 1

2

3

4

S. Allassonni`ere , E. Kuhn , J.T. Ratnatather , A. Trouv´e 1

CMAP, Ecole Polytechnique, Palaiseau, 2MIA, INRA, Jouy-en-Josas, 3 CIS, JHU, Baltimore, 4CMLA, ENS Cachan http://www.cmap.polytechnique.fr/∼allassonniere

[email protected]

Abstract: This paper aims at validating a methodology which has been proposed in [1, 2] to estimate a Bayesian Mixed Effect (BME) Atlas, i.e. the coupled templates and geometrical metrics for estimated clusters, in a statistically consistent way given a sample of images. We recall the generative statistical model applied to the observations which enables a simultaneous estimation of the clusters, the templates and geometrical variabilities (related to the metric) in the population. Following the idea of [1, 2, 3], we work in a Bayesian framework, use a Maximum A Posteriori estimator and approach its value using a stochastic variant of the Expectation Maximisation (EM) algorithm. This method is validated with a data set of 3D biomedical images of dendrite spines from a mouse model of Parkinson’s disease. We show the performance of this estimation on the estimated template and on the geometrical variability as well as on the clustering.

1 Introduction

6 The Bayesian Model

Estimated template with 30 images from only 3 types:

In the field of Computational Anatomy, one aims at segmenting images, detecting pathologies and analysing the normal versus abnormal variability of segmented organs. The most widely used techniques are based on the comparisons between subjects and a prototype image (usually called template in the literature). Such a prototype is an image whose biological properties (pathologic or control, segmentation of the tissues and structures, areas of functional activations, etc) are known and which - in a sense to be defined - characterises the population which is studied. This template contains common features of the population which would not be highlighted by multiple inter-subject comparisons. Our special interest is the construction of a statistically consistent atlas, called Bayesian Mixed Effect (BME) Atlas, as the estimation of the templates and their global geometric variabilities in estimates cluster for a given population in a statistically consistent way.

Left: Continuous template estimated from the algorithm. Right: Thresholded estimated template to better visualise the shape.

b One component geometry 7 MAP Estimation and Theoretical results

2 Advantages of this modelling • These models do not only explain data but are also able to randomly generate new data: can exhibit new features. • Population is heterogeneous: automatic decomposition into several sub-groups. • The deformation is applied to the observations instead of the template: neither interpolation nor non-iid noise. • The deformations are considered as hidden variables: not as nuisance parameters which have to be optimised.

3 The Observation Model We consider a population of n gray level images which we aim to automaticaly cluster in a small number T of groups called components later. We assume that each observation yi belongs to an unknown component τi. We work within the small deformation framework [4] so that conditional on the unobserved image membership t, there exists an unobserved deformation field z : R3 → R3 of a continuously defined template It : R3 → R and a Gaussian centred white noise ǫ of variance σt2 such that where Λ y(s) = It(xs − z(xs)) + ǫ(s) = zIt(s) + ǫ(s) , is a discrete grid of pixels and the pixels location is denoted by (xs)s∈Λ.

4 The Template and Deformation Model The templates It and the deformation z belong to Vp and Vg , two RKHS with respective kernels Kp and Kg : Given (pk )1≤k≤kp and (gk )1≤k≤kg ∃αt ∈ Rkp and (β (1), β (2)) ∈ Rkg × Rkg such that: It(x) = Kpαt(x), = z(x) = (Kg β)(x) =

kp P

k=1 kg P

k=1

Kp(x, pk )αt(k) ,

Kg (x, gk )(β (1)(k), β (2)(k)).

a One component template

Simulated (synthetic) shapes obtained using the previous estimated template and then thresholded.

ˆ ρ) ˆ = arg max(θ,ρ) q((θ, ρ)|y1n ). Estimator: (θ, Properties: • Existence of the Maximum A Posteriori estimator given any n−sample • Consistency of the MAP estimator as n goes to infinity

• Convergence of the algorithms (for either one or more components)

8 Estimation with the MCMC-SAEM algorithm ˆ ρ) ˆ = arg max(θ,ρ) q((θ, ρ)|y1n ). (maximum likeliTo compute (θ, hood in the setting of unobserved variable) −→ EM algorithm. BUT : • E step untractable −→ SAEM algorithm.

• Direct simulation of the missing variable w.r.t. the posterior: not possible −→ MCMC-SAEM algorithm.

The obtained shapes are similar to these of the training set (up to a smoothing due to the choice of some hyper-parameters in the model, for example kp, kg ). The algorithm managed to capture the large geometrical variability of the population.

c Experiments: Two component model Estimated templates (after thresholding) using 30 images from 3 types.

8b MCMC-SAEM algorithm (one component) This yields to the three following steps for iteration l: • Draw the missing data using a transition probability of a convergent Markov Chain having the posterior distribution as stationary distribution: Simulation step: β l+1 ∼ Πθl (β l , .).

• Approximate the complete log-likelihood using the previous simulations: Stochastic approximation step: Ql+1(θ) = Ql (θ) + ∆l [log q(y, β l+1|θ) − Ql (θ)] where (∆l ) is a non increasing sequence with limit 0 of positive step-size.

The estimated templates show the large heterogeneity of the population which has been captures. However, the small sample size does not lead to a very sharp estimation of the three cluster.

d Experiments: Two component model Estimated templates (after thresholding) using 50 images from the 6 types.

• Parameter update: Maximisation step: θl+1 = arg max Ql+1(θ). Remark: In the multi-component algorithm, the parameters are (θ, ρ) and the missing variables are (β, τ ). Due to a numerical (trapping state) problem, it is required to use several Markov Chains in parallel to simulate the missing variables (Simulation step).

5 Parameters and Hierarchical Model Model parameters: θ = (θt = (αt, σt2, Γtg ))1≤t≤T where T = number of components, weight of the different mixtures: ρ = (ρt)1≤t≤T . For each observation yi we consider the pair of unobserved variables ξi = (βi, τi). The hierarchical Bayesian structure of our model is :  2 , Γt ) T (ν ⊗ ν )  ρ ∼ ν , θ = (α , σ ∼ ⊗ ρ g t t g 1≤t≤T  t=1 p   T  P    τ1n ∼ ⊗ni=1 ρtδt | ρ t=1

n ∼ ⊗n N (0, Γτi )| τ n  β g  1 1 i=1       y n ∼ ⊗n N (z I , σ 2 Id ) | β n, τ n βi αi τi Λ 1 1 1 i=1 where the density functions are given by a Bayesian model.

Introduction of weakly informative priors:  !ag  1  −1 0  νg (dΓg ) ∝ exp(−hΓg , Γg i/2) p dΓg , ag > 2kg + 1    |Γg |       ! !ap   2  σ0 1  2 × νp(dσ , dα) ∝ exp − 2 √ 2 2σ σ      t(Γ )−1(α − µ ) dσ 2dα  exp (α − µ )  p p p    aρ    T  Y     ρ(τ ) . νρ(ρ) =    τ =1

9 Experiments: Data set Training set: Murine dendrite spines images: 50 binary images of microscopic structures, tiny protuberances found on many types of neurons termed dendrite spines. The images are from control mice and a population of mice which have been genetically modified to mimic human neurological pathologies like Parkinson’s disease. The acquisition process consisted of electron microscopy after injection of Lucifer yellow and subsequent photo-oxidation. The shapes were then manually segmented on the tomographic reconstruction of the neurons. The images are labelled by experts as belonging to six different categories (called types): double, filopodia, long mushroom, mushroom, stubby and thin.

The large geometrical variability is higher with 6 types and leads to two very different shapes. However, one would need more images of all types and estimate 6 components to expect to get the six clusters. With only 50 images, more than two components would not be relevant for a statistical estimation of high dimensional parameters.

10 References

References [1] St´ephanie Allassonni`ere, Yali Amit, and Alain Trouv´e. Toward a coherent statistical framework for dense deformable template estimation. Journal of the Royal Statistical Society, 69:3–29, 2007. [2] St´ephanie Allassonni`ere, Estelle Kuhn, and Alain Trouv´e. Bayesian deformable models bulding via stochastic approximation algorithm: A convergence study. to appear in Bernoulli Journal. [3] St´ephanie Allassonni`ere and Estelle Kuhn. Stochastic algorithm for bayesian mixture effect template estimation. to appear in ESAIM Probab.Stat. [4] Y. Amit, U. Grenander, and M. Piccioni. Structural image restoration through deformable templates. Journal of the American Statistical Association, 86:376–387, 1989.