S Durrleman

X Pennec

A Trouvé

N Ayache

J Braga

Comparison of the endocast growth of chimpanzees and bonobos via temporal regression and spatiotemporal registration

In this paper, we aim at characterizing and quantifying the differences between the growth of bonobos (Pan paniscus) and chimpanzees (Pan troglodytes). We use a collection of endocasts of wild-shot animals of both species. Each sample has been associated with a dental age, as a common temporal marker. To compare the endocasts, we used the current-based metric which allows us to quantify the shape differences without the need to find homologous landmarks on the surfaces. First, we perform a temporal shape regression, which estimates a typical growth scenario of the endocast for the bonobos and the chimpanzees. Then, a spatiotemporal registration scheme is used to quantify the differences between these two growth scenarios. The variations are decomposed into one morphological deformation and one time warp. The morphological deformation accounts for the anatomical differences independently of the age. The time warp accounts for the change of the dynamics of growth. It shows that the growth speed of the bonobos at juvenility is more than twice less than the one of the chimpanzees. This estimation gives more insights into the developmental delay observed in the bonobos growth.

Introduction

This paper aims at characterizing the differences of growth patterns between the two closest human relatives: the bonobo (Pan paniscus) and the chimpanzee (Pan troglodytes), two species of the genus Pan. According to phylogenetic studies [START_REF] Won | Divergence population genetics of chimpanzees[END_REF], these two species share a common ancestor, which used to live at least two millions years ago. By comparison, the common ancestor of humans, chimpanzees and bonobos used to live at least ten millions years ago.

Bonobos were discovered to science in 1929. Since then, both morphological and behavioral studies tend to show that the bonobo is a "juvenilized" version of the chimpanzee [START_REF] Schwarz | Das Vorkommen des Schimpanzen auf then linked Kongo-Ufer[END_REF]. Indeed, body growth is considered to be retarded in bonobos compared to chimpanzees, a feature corresponding to their delay in motor development during the first years of postnatal life [START_REF] Kuroda | Developmental retardation and behavioural characteristics of pygmy chimpanzees[END_REF]. The adult bonobo skull shows a decreased facial prognathism and teeth with a reduced sexual dimorphism [START_REF] Shea | Heterochrony in human evolution: The case for neoteny reconsidered[END_REF]. Bonobos are also characterized by a longer dependency of the child on the mother [START_REF] De Waal | Bonobo sex and society[END_REF]. These studies lead to put forth the "bonobos hypothesis", which assumes the bonobo to be a dwarfed version of the chimpanzee. To investigate to which extent this hypothesis may be true, one must carefully analyze the difference between the growths of both species, called ontogenesis. Unfortunately, the lack of reliable data, in particular for the bonobos, has made difficult so far the study of the ontogenesis of both species and their comparison.

In this study, we focus on the endocasts of the specimens, which is a mould of the endocranium. The surface of this mould provides a replica of the inner surface of the skull, and therefore has often played an important role for the study of the evolution of the brain in fossil mammals [START_REF] Radinsky | Brains of early carnivores[END_REF][START_REF] Macrini | Cranial endocasts from a growth series of Monodelphis domestica (Didelphidae, Marsupiala): A study of individual and ontogenetic variation[END_REF]. We use here a cohort of about 60 endocasts of wild-shot animals of each species, which is among the largest collection available. It has been observed in [START_REF] Kinzey | The dentition of the pygmy chimpanzee[END_REF] that the sequences of teeth emergence in bonobos and chimpanzees are essentially identical. This gives a way to estimate the "dental age" of each skull, as a common growth marker.

In this paper, we detect and compare geometrical growth patterns in these data using the methodology recently introduced in [START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]. First, a shape regression method enables to estimate a typical growth scenario for each species. Second, a spatiotemporal registration scheme is used to compare the two growth scenarios both in terms of morphological changes and in terms of different growth speeds. This yields a quantitative estimation of the possible morphological delay of the bonobos with respect to the chimpanzees growth, which is the main feature we want to measure in regards to the bonobos hypothesis.

Last but not least, the implementation of this methodology uses the metric on currents to quantify shape dissimilarity [START_REF] Vaillant | Surface matching via currents[END_REF][START_REF] Glaunès | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] (although other metric could be also compatible with the proposed methodology). This allows us to compare endocasts without relying on the choice of landmarks points on the surface. This differs from usual methods used in physical anthropology [START_REF] Lieberman | A geometric morphometric analysis of heterochrony in the cranium of chimpanzees and bonobos[END_REF][START_REF] Neubauer | The pattern of endocranial ontogenetic shape changes in humans[END_REF], which uses the so-called geometric morphometric analysis as initiated in [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF]. Here the shape comparison is expert-independent. The subsequent analysis is not biased by the choice of the landmarks and the inherent inter-expert variability. The geometry of the endocast is considered as a continuous surface and includes the direction of the normals. Moreover, the analysis of shape differences can be used afterwards to drive the search of significant anatomical landmarks like the ones which highlight the most striking differences between species.

The paper is organized as follows. In Section 2, we present the methodology used for the analysis of the endocasts data. In Section 3, we show that this method enables to estimate a typical growth scenario for the two species and to measure both their morphological differences and their different growth speeds.

Methodology

As explained in [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF], the proposed methodology is subject-specific. This means that every subject, here every species, is supposed to evolve in the same atemporal morphological space (a given reference frame provided by its own coordinate system). The growth of each species is described by a time-varying function χ t (x) (for any time t and point in space x) in its own coordinate system (See Fig. 2). The growth scenarios of two species may differ in two ways. First, different species have different morphologic features, independently of the age of the individuals. Second, the dynamics of the growth scenarios may be different for different species: one anatomical development may occur at different speed, thus leading to developmental delay or advance. The first difference is described by a morphological function φ(x) (spatial deformation) which maps the reference frame (i.e. morphological space) of one species into the one of the second species. The change of dynamics is described by a time-warp ψ(t) which is a 1D-diffeomorphism (smooth deformation which prevents time-reversal) which maps the age of one species t to the corresponding age of the second species ψ(t). This function changes the speed the growth of the species in a non-linear way. Eventually, the differences between two growth scenarios is described by the spatiotemporal deformation (φ, ψ) which best aligns the scenarios.

Two subjects have been scanned twice (one subject per rows, one time-point per column).

Fig. 1. In the subject-specific approach, one considers that one subject is "circle" and the other is "square": the difference between both subjects is a single function φ, which maps circles to squares, independently on time. The evolution of the first subject is described by a function χ which maps a small circle to a big circle. The evolution of the second subject is then described by the function χ φ which maps a small square to a big square. Once registered into the "morphological space" of the squares, the circle evolution may still differ from the square one by a change of the speed of the evolution. This variations will be captured by the time-warp ψ.

Temporal shape regression

The growth scenario of one species is described by the growth function χ t , a 3D-deformation of the underlying space. Given a baseline shape S 0 , this function induces a continuous shape evolution S(t) = χ t (S 0 ) for t in the time interval of interest. Given a discrete set of observed shapes (S i ) associated to time-points t i and a baseline S 0 , one aims at estimating the continuous growth function χ t in a least square sense (i.e. by minimizing the sum of squared differences between each shape S i and the growth scenario at time-point t i : S(t i ) = χ ti (S 0 )):

J(χ) = ti d(χ ti (S 0 ), S i ) 2 + γ χ Reg(χ) (2.1)
where d is a similarity measure between shapes, Reg(χ) a regularity term and γ χ a trade-off between regularity and fidelity to data. Here we use the distance on currents as the metric d, which allows us to measure the discrepancy between shapes like curves or surfaces without the need for introducing landmarks on shapes. This metric depends on a scalar parameter λ W , which is the typical spatial scale under which geometrical variations are considered as noise (see [START_REF] Vaillant | Surface matching via currents[END_REF][START_REF] Durrleman | Inferring brain variability from diffeomorphic deformations of currents: an integrative approach[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] for more details). The regularity term is given as the total kinetic energy of the deformation [START_REF] Glaunès | Transport par difféomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]. In the framework of the RKHS [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF], the regularity of the deformation is controlled by a scalar parameter λ χ , which is the typical scale under which points move in a correlated manner. This selects a group of deformations whose spatial frequencies fall within a given band. The more concentrated the "spectrum" of the deformation at the band limit, the higher the regularity term Reg(χ). Once the scale λ χ is fixed, the trade-off γ χ controls the goodness of fit achieved by the deformation within the group which minimizes the criterion. For γ χ → ∞, the resulting deformation will be identity map (no deformation) independently of the scale λ W . For γ χ → 0, data will be fit as closely as possible without penalizing irregular deformations within the group. The baseline S 0 is usually chosen as the smallest shape associated to the earliest age.

As an illustrative example, we applied this temporal shape regression method on a set of 2D profiles of hominid skulls (segmented from images obtained at www.bordalierinstitute.com). Results are shown in Fig. 2.

Spatiotemporal registration

Once the temporal shape regression has been performed separately for each species, we end up with two growth scenarios: B(t) and C(t) (described by two growth functions χ b t and χ c t ). The morphological deformation φ maps the chimpanzee morphological space to the bonobos one: is deforms the shapes C(t) into φ(C(t)) for any age t. The time-warp ψ(t) change the dynamics of the chimpanzee growth by mapping the age t to the age ψ(t). The spatiotemporal deformation of the chimpanzee growth is given as: C(t) = φ(C(ψ(t))). One aims at estimating the spatiotemporal deformation (φ, ψ) which best aligns the chimpanzees growth to the bonobos one. For this purpose, we sample the evolution of B as a collection of B(τ i ) for linearly spaced time-points τ i . A Maximum A Posteriori estimation leads to the minimization of the criterion [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]:

J(φ, ψ) = τi d(φ(C(ψ(τ i ))), B(τ i )) 2 + γ φ Reg(φ) + γ ψ Reg(ψ) (2.2)
where d is the distance between shapes given by the metric on currents like for the regression function, Reg(φ) and Reg(ψ) the regularity terms of the morphological deformation and the time-warp, γ φ and γ ψ two scalar parameters which balance the weight of the fidelity-to-data and the regularity terms. The regularity terms are the total kinetic energy of the diffeomorphism in the 3D domain for φ and 1D domain for ψ. In the same framework as the regression function, they depend on two scalar parameters: λ φ which determines the spatial scale under which two points move consistently according to φ and λ ψ which determines the scale under which the time-points are moved in a correlated manner along the time axis. In [START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF][START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF], we propose a systematic way to construct 1D-diffeomorphisms which prevents time-reversal. Note that ψ(t) is greater than t if the deformed scenario is delayed with respect to the original one, it is greater than t if the deformed scenario is in advance. The slope of the graph of ψ at age t (ψ ′ (t))

indicates whether the deformation is accelerating (ψ ′ (t) ≥ 1) or slowing down (ψ ′ (t) ≤ 1) the growth scenario.

As an illustrative example, we artificially divide the set of profiles of hominid skulls into two groups, which play the role two different subjects (see Fig. 3). The result of the spatio-temporal registration between these two subjects is illustrated in Fig. 4a. The time-warp is a 1D-function which can be plotted as a graph as shown in Fig. 4b. The slope greater than 1 between the target data indicates an acceleration of the target with respect to the source. Its value provides an an estimated measure of the acceleration. 3, second row: the source evolution is deformed by the estimated "morphological change" function (the same 3D deformation is applied to every frame of the sequence); third row: the source evolution is then accelerated according to the estimated time-warp (the frames are shifted along the time axis, as shown by the dashed black lines, while the geometry of the shapes remain unchanged). The morphological deformation shows that the skulls in the target evolution are rounder and the jaw less prominent than in the source evolution. The time warp shows an acceleration of the target evolution w.r.t. the source evolution. b-Graph of the time-warp (dashed line indicates the no-dynamical change axis, i.e. x = y axis). The slope of the curve measures the acceleration of the evolution: here the target evolution occurs 1.66 times faster than the source evolution.

Experimental results

Material

As part of the collaborative project ARC 3D-Morphine, we used a set of endocasts of 59 chimpanzees and 60 bonobos. The original skulls are from the collection of "Musée de l'Afrique centrale" in Tervuren, Belgium (curator: E. Gilissen). They represent wild-shot individuals with approximately equal numbers of male and female. They have been scanned using a Siemens Somatom Esprit Spiral CT, with slice thickness between 0.33 and 0.50mm. The segmentation of the endocasts using itkSNAP [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF] leads to surface meshes, as shown in Fig. 5. These surfaces have been rigidly co-registered using gmmreg [START_REF] Jian | A robust algorithm for point set registration using mixture of Gaussians[END_REF].

The analysis of the permanent teeth development of the skulls provides an estimate of the age of the death of the samples, which we call here "dental age". It has been observed in [START_REF] Kinzey | The dentition of the pygmy chimpanzee[END_REF] that the sequences teeth emergence in bonobos and chimpanzees are essentially identical. Each skull is therefore associated to one the 6 dental ages defined in [START_REF] Shea | Heterochrony in human evolution: The case for neoteny reconsidered[END_REF]: infant, child, young juvenile, old juvenile, subadult and adult. To refine the classification, some skulls have been associated the intermediate class 'child/young juvenile' by the experts. Age distribution is shown in Fig 5.

Typical growth scenario estimation

Without loss of generality, we assume that each of the dental age lasts the same amount of time, namely 5 time steps. Therefore the time interval of interest has been divided into 30 time-steps and the endocasts are associated to the time-point t i = 5, 10, 15, 20, 25, 30 according to their dental age. The age child/young juvenile has been associated to the time-point 13. We choose the smallest endocast within the child class as the baseline S 0 and associate it to the time point t = 2. Then, we perform a temporal shape regression of the endocasts (independently for each species) with respect to the dental age by minimizing the cost function (2.1). We set the typical spatial interaction between currents λ W = 20mm, the spatial scale of deformation consistency λ χ = 50mm and the trade-off between fidelity-to-data and regularity γ χ = 10 -4 mm 2 (unit of time). The diameter of the endocasts are typically between 60 and 70mm.

Results are shown in Fig. 6. The scenarios reveal that the endocast growth is not a simple scaling but involves non-linear and anisotropic effects. The most salient effect, besides the increase of volume, is an elongation along the posterior/anterior axis and a slight contraction along the superior/inferior axis. As a consequence, the geometry of the endocast, which is almost spherical at birth, becomes more and more ellipsoidal. These observations hold for both species, although it seems that the chimpanzees endocasts have a strongest anisotropy and that this anisotropy increases faster in time. The subsequent spatiotemporal registration will measure these differences more precisely.

By contrast, the two growth scenarios seem to differ a lot at infancy and childhood. This difference is due to the little amount of data in infancy. The only two infant chimpanzees have a larger endocasts compared to both the infant bonobos and the children chimpanzees. To have a more relevant estimation of the growth in infancy, we expect to scan more infant chimpanzees skulls in the future. Note that in the next section we will not take the infancy data into account and will consider the growth scenarios starting at childhood.

We can deduce from the growth scenario an estimation of the evolution of the endocranial volume across ages, as shown in Fig. 6 (Note that we have not performed a regression of the volume but of the shapes). Besides the evolution in infancy, one intriguing feature is the apparent decrease of bonobos endocranial volume at sub-adulthood. This feature is also present in the endocranial volume distribution in the original endocasts (mean and standard deviation are shown in Fig. 6): the mean of the volume at sub-adulthood is smaller than the one of old juveniles. However, the Mann-Whitney U test gives a p-value of 0.47 when comparing the volume distribution of old juveniles and sub-adults: the median of the two distributions are not proved to be statistically different. The test run for every pair of consecutive distributions shows a significant increase of volume only between infancy and childhood for the bonobos (p-value: 9 10 -3 ) and between old-juvenility and sub-adulthood for the chimpanzees (p-value: 0.02).

Spatiotemporal registration between the growth scenarios

We perform a spatiotemporal registration between the two growth scenarios starting at childhood for the reasons explained in the previous section. We consider the chimpanzee growth scenario as the reference scenario. The bonobo scenario is sampled every 2 time-steps. These samples play the role of the target shapes. We set the scale of currents to λ W = 20mm as for the regression estimation. We run the registration for different set of parameters and pick the ones which enable to achieve the smallest discrepancy term in (2.2). This gives the consistency scale of the morphological deformation: λ φ = 10mm, the consistency scale of the time-warp λ ψ = 1 unit of time, morphological trade-off γ φ = 25 10 -9 mm 2 and the temporal trade-off γ ψ = 4 10 -9 mm 4 /(unit of time).

The morphological deformation changes the shape of each frame of the chimpanzee growth as shown in Fig. 7a (this is the equivalent figure to the first and second row in Fig. 4a). It shows that, independently of the age, the bonobos endocast are rounder than the chimpanzees ones. The movie of this deformation clearly shows a twist at the anterior and posterior part of the endocast. This effect could be quantify, for instance by measuring the ratio between the elongation in the superior-inferior and the posterior-anterior direction. This shows how our analysis can be used to drive the search of relevant anatomical features, in contrast to methods like [START_REF] Bookstein | Morphometric tools for landmark data: geometry and biology[END_REF] which select the information a priori.

The graph of the estimated time-warp is shown in Fig. 7b. Its shows the ages of the bonobos match the chimpanzees ones. An important speed reduction of the growth of the bonobos with respect to the chimpanzees occurs between old-juvenility and sub-adulthood. The almost constant slope of the curve during this period of time indicates that the bonobos growth speed is 0.42 times that of the chimpanzees. The graph shows also that the bonobos seem to be slightly in advance with respect to the chimpanzees at childhood and that the delay of the bonobos growth at sub-adulthood seems to be reduced at adulthood.

Discussion and conclusion

The experiments on this set of endocasts showed that the proposed methodology enables to detect, measure and compare growth patterns in a set of timeindexed shapes. The shape regression averaged the inter-individual variability and estimated a typical growth scenario for each species. The comparison of the growth scenarios using the spatiotemporal registration gave more insights into the "bonobos hypothesis". The measure of the morphological differences based on non-linear 3D deformation allowed us to detect anisotropic variations, which cannot be seen by the analysis of some simple geometrical features like the volume for instance. The non-linearity of the time-warp allowed us to detect not only the expected developmental delay of the bonobos but also more subtle effects like the developmental advance of bonobos at childhood or the delay reduction at adulthood. If we assume an ancestor relationship between both species, this measure suggests that the analysis of heterochrony cannot be reduced to a black-and-white decision between paedomorphosis and peramorphosis [START_REF] Lieberman | A geometric morphometric analysis of heterochrony in the cranium of chimpanzees and bonobos[END_REF]. This analysis is based on a generic methodology which reduces the dependence on expert choices like the locations of landmarks. The anatomical features of interest are not selected a priori but appears as the result of the analysis. Future work should focus on the extraction and quantification of the features which seems the most discriminative, like the "roundness" of the endocast for instance. The method depends on exactly 7 parameters, which have a physical or statistical interpretation [START_REF] Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF]. The parameters of the spatiotemporal registration have been set so that it reduces the dissimilarity between the deformed source and the target. Note that due to the modeling assumptions, there may not be any set of parameters which achieves a perfect data fit (i.e. zero residuals). We proposed here to measure the inter-species growth differences using the estimated growth scenarios. This has some drawbacks: one does not account for the fact that one has a different number of samples at each age and this approach is not symmetric as we chose the chimpanzee as the reference species. To overcome these issues, we will use the unbiased atlas construction scheme of [START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF] in the future. One also needs to investigate the variability of the results with respect to the parameters and the samples. For instance, bootstrap procedure can be used to estimate a confidence interval around the growth scenarios. and S. Prima for their active role in the project 3D-Morphine and G. Gerig for his insightful comments. We thank the anonymous reviewers for their relevant suggestions. This work has been partly funded by the INRIA ARC 3D-Morphine, the European IP project Health-e-child (IST-2004-027749) and Microsoft Research.

b-Time Warp

Graph of time-warp showing that the bonobos development is in advance with respect to the chimpanzees one at childhood and then that it drastically slows down during juvenility (almost linearly by a factor 0.42 between old-juvenility and sub-adulthood). This delay seems to decrease at adulthood.

Fig. 2 .

 2 Fig.2. Example of shape regression on profiles of hominid skulls (in red). We choose the australopithecus profile as the baseline S0. The temporal regression computes a continuous shape evolution S(t) (in blue) such that the deformed shape matches the observations at the corresponding time-points. If several shapes are associated to the same time-point, they are averaged during the estimation of the regression.

Fig. 3 .Fig. 4 .

 34 Fig.3. We aims at comparing the evolution {Homo habilis-erectus-neandertalensis} (red shapes) to the evolution {Homo erectus-sapiens sapiens} (green shapes). The regression scheme estimates the growth scenario of the source (sequence of blue shapes). The spatiotemporal deformation of the source evolution which best matches the target shapes is shown in Fig.4.

Fig. 5 .

 5 Fig. 5. Samples of the original endocasts of both species (top) and their distribution according to their dental age (bottom)

  The morphological deformation is applied to the chimpanzees endocasts at 3 ages. Left column: the endocasts from the chimpanzees growth scenario. Right: the deformed endocasts. Middle: an intermediate deformation. This shows that, on average, the endocast of a chimpanzee is more elongated and less round than the one of a bonobo (more clearly visible on the movie of the deformation). Note that the deformed endocasts do not match the ones of the bonobo growth at the same age, but at the age given by the correspondence graph shown in b-The dashed black line represents the non dynamical changes axis (x = y). When the time-warp (magenta curve) is above (resp. below) the dashed line, bonobo growth is in advance (resp. delayed) w.r.t. the chimpanzee growth. Slope greater than 1 (resp. smaller than 1) indicates that the bonobo growth accelerates (resp. slows down) w.r.t. to the chimpanzee growth.

Fig. 7 .

 7 Fig. 7. Spatiotemporal registration between the chimpanzees and the bonobos growth scenario.

Acknowledgments We thank B. Combès for data pre-processing, Wim Van Neer as the former curator of the 'Musée de l'Afrique Centrale' at Tervuren, G. Subsol

Temporal regression of endocasts of bonobos


Christof Seiler

Xavier Pennec

Mauricio Reyes

PARAMETRIC REGRESSION OF 3D MEDICAL IMAGES THROUGH THE EXPLORATION OF NON-PARAMETRIC REGRESSION MODELS

Keywords: Non-Parametric Regression, Parametric Regression, Log-Euclidean Framework, Diffeomorphic Deformations, Femur

Currently there is an increase usage of CT-based bone diagnosis because low-radiation and cost-effective 2D imaging modalities do not provide the necessary 3D information for bone diagnosis. The fundamental objective of our work is to build a model connecting 2D X-ray information to 3D CT information through regression. As a first step we propose an univariate non-parametric regression on individual predictor variables to explore the non-linearity of the data. To later combine these univariate models we then replace them with parametric models. We examine two predictors, shaft length and caput collum diaphysis angle on a database of 182 CT images of femurs. We show that for each predictor it is possible to describe 99% of the variance through a simple up to second order parametric model. These findings will allow us to extend to the multivariate case in the future.

INTRODUCTION

Up to now, 3D bone anatomy has been generated from X-ray images using computational tools. However these tools have focused only on 3D bone shape reconstruction while little attention has been given to 3D reconstruction of bone mineral density, which is important in analysis of bone fragility, orthopedic surgery, orthopedic implant design, etc. We hypothesize that full 3D bone anatomy (i.e. bone shape and bone mineral density) for X-ray based diagnosis can be achieved through development of computational and statistical tools, making use of vast amount of femur CT images.

The fundamental objective is to perform multivariate regression on the anatomy of bones. As predictors, the regression model uses patient-specific metadata (e.g. age, weight, body mass index, etc.), and image features extracted from patient radiographs. The variables to predict are patient-specific 3D CT images of bones. As a first step, in this paper we propose univariate parametric regression models based on the exploration of non-parametric regression results for femur bones on a selected subset of clinically meaningful morphological parameters: Femoral shaft length depicted in Fig. 1 (a,b) and caput collum diaphysis (CCD) angle in Fig. 1 (c). Currently in medical image regression, parametric (e.g. [1, 2, 3]) and non-parametric (e.g. [START_REF] Davis | Population Shape Regression from Random Design Data[END_REF]) approaches have been employed.

• On one hand, parametric-based approaches can be utilized for multivariate regression without encountering the curse of dimensionality problem, but its use preimposes a statistical structure on the data.

• On the other hand, non-parametric approaches do not impose a certain structure, but they suffer from the curse of dimensionality when trying to optimize for hyper-parameters in high-dimensional spaces. However, provided low-dimension subspaces can be found, efficient optimization could be performed. Nevertheless, this is still an open research question [START_REF] Gerber | On the Manifold Structure of the Space of Brain Images[END_REF][START_REF] Steinke | Non-Parametric Regression Between Riemannian Manifolds[END_REF].

Therefore, we selectively combine these two approaches to solve the problem in a low-dimensional space, without preimposing a statistical data structure. For this goal we propose a two-step approach. First, we explore the data structure through univariate non-parametric regression. This step enforces no assumptions on the data structure while avoiding the curse of dimensionality problem. Second, we parametrize the explored structures.

In [START_REF] Davis | Population Shape Regression from Random Design Data[END_REF] the authors showed kernel regression formulated with Fréchet weighted means to take into account the non-Euclidean nature of diffeomorphisms endowed with a right invariant metric (LDDMM) and applied it to images of the brain. In contrast, in this work we formulate kernel regression in an Euclidean way in a Log-Euclidean framework. This simplifies and speeds-up the process significantly while still taking into account a large part of the non-Euclidean nature of the manifold-valued data. Furthermore, this simplification allows for other computations that are out of reach to LDDMM, such as determining the optimal kernel bandwidth through cross-validation. The Log-Euclidean framework uses stationary velocity fields to parametrize a diffeomorphic deformation, whereas in [START_REF] Davis | Population Shape Regression from Random Design Data[END_REF] non-stationary velocity fields are used. Although the theory shows that not all diffeomorphic deformation can be reached with stationary velocity fields; there is no indication so far that this affects the anatomical shape analysis in any way.

In the following, we develop the methodology and show results obtained on femur CT images.

METHODS

The methods description will be divided in five parts: Logdomain registration, Log-Euclidean statistics, non-parametric kernel regression, cross-validation and parametrization of principle component (PC) scores.

Log-Domain Registration

To setup correspondences between anatomical images, a set of images are registered to a reference. We use the novel symmetric diffeomorphic registration approach described in [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF]. What is new in this registration framework is the efficient optimization in the log-domain. As a consequence, the results of the registration are so-called stationary velocity fields. These velocity fields can be looked at as generators for diffeomorphic deformations through the group exponential map that can be very efficiently computed using the scaling and squaring method [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF].

Log-Euclidean Statistics

Applying the Log-Euclidean framework [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF] on these fields allows us to compute statistics, e.g. averages, and still preserve diffeomorphism. In the Log-Euclidean framework, velocity fields are regular elements in a vector space; this allows us to use simple Euclidean statistics instead of more complex non-linear techniques, which we needed when working in the LDDMM space of diffeomorphic transformations. To map resulting velocity fields into diffeomorphic transformations the exponential is calculated. To go from diffeomorphic transformation back to velocity fields, a logarithmic mapping is performed. However, in many cases the intrinsic parametrization of the transformation by its log in the log-domain registration allows to avoid this numerically unstable step. For a detailed survey of the methodology we refer to [START_REF] Pennec | Statistical Computing on Manifolds: From Riemannian Geometry to Computational Anatomy[END_REF].

Non-Parametric Kernel Regression

We use a kernel regression method to compute the deformation of the template that best predicts the images based on prediction variables x, in our case shaft length and CCD angle. Our kernel regression function is

mσ (x) = exp N i=1 K σ (x -x i )v i N i=1 K σ (x -x i ) , (1)
where N , x i , K σ , exp and v i are the total number of images, shaft length or CCD angle for image i, a Gaussian kernel function with σ bandwidth, the mapping from velocity fields to diffeomorphic deformations and the i th velocity field, respectively.

Cross-Validation

The quality of kernel regression methods strongly depends on the selection of bandwidth parameters. To select a bandwidth parameter we apply cross-validation with penalty functions. The penalty and corresponding weighting functions penalize very small bandwidth values. Bandwidth values equal to zero are not interesting because they are just a nearest neighbor interpolation of the data. In our case we solve the following minimization problem:

σ = argmin σ∈R N i=1 || log( mσ (x i )) -v i || 2 Ξ(W σ,i (x i )), (2)
where N , log, σ, Ξ are total number of images, mapping from diffeomorphic deformations to velocity fields, bandwidth and penalty function, respectively, and

W σ,i (x i ) = K 1 (0)/ N j=1 K 1 (σ -1 (x i -x j ))
is the weighting function. For details we refer to [START_REF] Härdle | Applied Nonparametric Regression[END_REF]. By solving this optimization problem we obtain a kernel bandwidth greater than zero which minimizes the regression function's prediction error for all images N .

Parametrization of Univariate Kernel Regression

To parametrize the non-linear regression function that we have established via kernel regression, the regressed velocity fields are reduced in dimension using principle component analysis (PCA). The data points are then projected onto a low-dimensional space covering 99% of the variance and evaluated for possible parametrization. In all cases that we have observed so far it is possible to fit a polynomial to each PC. The procedure can be summarized as follows:

1. Regress velocity fields 2. PCA on regressed velocity fields 3. Plot scores for each PC

Fit parametric function to each score plot

The new parameterized regression model is then:

m(x) = exp µ + M i=1 p i (x)z i , ( 3
)
where µ is the mean regressed velocity field, p i (•) is the ith polynomial function representing ith PC scores, z are PCs and M the number of PCs describing 99% of the variance.

Regression is done on 50 predictor values drawn from a normal distribution we obtained on the original data. This way we avoid conducting PCA on the original data and still get a reasonable representation of the original data structure.

RESULTS

Considering generalized cross-validation (GCV), Ξ(u) =

(1 -u) -2 , as the penalty function in (2) we obtain 2.1 for shaft length and 1.6 for CCD angle. With these bandwidth values, prediction is performed on a dataset of 182 left femur CT images.

To validate the results, 50 velocity fields were predicted with prediction values drawn from the normal distribution of the original data. Then the exponential was taken of these fields and the shaft length and CCD angle were measured in the exponentiated fields.

Shaft Length Regression

To find the best parametric model, different polynomials were evaluated. Our evaluation showed that the best fit is linear for all PCs. In Fig. 2 the score plots are depicted and in Tab. 1 four different parametrization and the root mean square (RMS) prediction error are listed. Parametric model 4 performs similar to the non-parametric model.

In Fig. 3

Caput Collum Diaphysis Angle Regression

Similar to the shaft length low order polynomial give better prediction results. The best prediction results are reached using a linear parametrization for the first PC and quadratic for the second (Tab. 2). This parametric model even outperforms the non-parametric model. See Figs. 4 and5 for score and prediction error plot.

DISCUSSION AND CONCLUSIONS

Only three PCs for shaft length and two for CCD angle are needed to describe 99% of the variance in the velocity fields. This fact enables us to evaluate each predictor individually. Surprisingly, in both experiments low order polynomial parameterizations provided the best results. One possible reason for this could be the behavior of polynomials at the boundary data points. In both cases the parametric slightly outperforms the non-parametric model. It seems that the non-parametric model overfits on the training data. Therefore, we gain in three ways from parametric models: Simplicity of the model, computational efficiency and reduction of overfitting.

We plan to evaluate other parametric functions to avoid possible problems at boundary data points. The shown validation can only test the consistency of the method, therefore we will further validate our method by comparing the prediction to manual segmentation results. In future work, simplicity and efficiency, will allow us to use univariate exploration methods to build multivariate regression models with two and more predictors.

[1] Yuhui M. Yang, Daniel Rueckert, and Anthony M. J.

Bull

Fig. 1 :

 1 Fig. 1: (a) Starting point of shaft length morphological parameter at the greater trochanter. (b) End point between condyles. (c) CCD angle.

Fig. 2 :Fig. 3 :

 23 Fig.2: Score values for the first three PCs.

Fig. 4 :

 4 Fig. 4: Score values for the first two PCs.

Fig. 5 :

 5 Fig. 5: Comparison of parametric and non-parametric regression results for CCD angle.

Table 1 :

 1 Shaft length parametrization comparison. PCn columns show the order of the polynomial fit used for the corresponding PC.

	the prediction values are compared in a more
	qualitative manner.
		PC1 PC2 PC3 RMS [mm]
	Parametric 1	1	2	5	1.37
	Parametric 2	1	2	3	2.26
	Parametric 3	2	2	5	0.97
	Parametric 4	1	1	1	0.83
	Non-parametric -	-	-	0.84

Table 2 :

 2 , "Predicting the Shapes of Bones at a Joint: Application to the Shoulder," Computer Methods in Biomechanics and Biomedical Engineering, vol. 11, no. 1, pp. 19-30, February 2008. CCD parametrization comparison. PCn columns show the order of the polynomial fit used for the corresponding PC.

		PC1 PC2 RMS [degrees]
	Parametric 1	2	3	1.14
	Parametric 2	1	2	0.88
	Parametric 3	1	1	0.98
	Non-parametric -	-	1.00


François Chung
email: francois.chung@inria.fr

Jérôme Schmid
email: jerome.schmid@miralab.ch

Nadia Magnenat-Thalmann

Hervé Delingette

Franc ¸ois Chung

Franc ¸ois Chung

Miralab

Comparison of statistical models performance in case of segmentation using a small amount of training datasets

Keywords: model based segmentation, statistical models, principal component analysis, clustering

Introduction

Building meaningful priors for model-based image segmentation purposes is an important topic in medical imaging. To account for a large variety in both shape and appearance, numerous datasets are usually required in a training stage. Atlas-based methods are designed in the hypothesis that a significant number of training datasets may yield a reasonable mean model (i.e. in the sense of a population mean) as well as meaningful modes of variation. Those methods usually define a reasonable estimation of the prior using the principal modes of variation (i.e. variations may be controlled with few parameters). They have been successfully used in coarse-to-fine approaches where the algorithm starts with a rough approximation of the prior (i.e. with few degrees of freedom) and evolves adding more variation until reaching a steady state [START_REF] Heimann | A shape-guided deformable model with evolutionary algorithm initialization for 3d soft tissue segmentation[END_REF][START_REF] Schmid | Mri bone segmentation using deformable models and shape priors[END_REF].

In the literature, shape variations are often described using Principal Component Analysis (PCA) [START_REF] Schmid | Mri bone segmentation using deformable models and shape priors[END_REF][START_REF] Seim | Automatic segmentation of the pelvic bones from ct data based on a statistical shape model[END_REF][START_REF] Heimann | A shape-guided deformable model with evolutionary algorithm initialization for 3d soft tissue segmentation[END_REF][START_REF] Behiels | Evaluation of image features and search strategies for segmentation of bone structures in radiographs using active shape models[END_REF][START_REF] Cootes | Data driven refinement of active shape model search[END_REF][START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF] of Point Distribution Models (PDM) [START_REF] Cootes | Statistical models of appearance for computer vision[END_REF]. From a sufficient number of shapes with known point correspondence, a mean shape is calculated. During the segmentation, a mesh initialized as the mean shape is deformed until it fits a desired target shape. PCA-based methods allow the mesh to be deformed only in the range of estimated variations. In this way, deformations occur in a predictable way since they originate from the principal variations observed from the datasets used to build the statistical model. Regarding the appearance, methods such as Active Shape Models (ASM) [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF] and Active Appearance Models (AAM) [START_REF] Cootes | Active appearance models[END_REF] were proposed to account for the main intensity variation around and within structures of interest respectively.

However, the number of datasets required to account for both the shape and appearance of a structure may be an issue. First, medical imaging acquisitions require time and resources, and this may limit the number of datasets available for training. Second, the large number of acquisition protocols and hardware characteristics (especially in case of versatile modalities such as MRI) tend to produce images with a large variety of intensity distribution for the same structure of interest. Additionally, noise and artifacts (e.g. patient movement and partial volume effect) are likely to corrupt the image intensity. Those factors strongly affect the construction of priors, as they bring meaningless intensity information into the appearance priors. Also, the manual segmentation of images by an expert, which is required for the extraction of shapes, is a tedious task and is a limiting factor for the availability of training datasets. Finally, the large natural variability of shape and appearance cannot be well represented by a Gaussian distribution assumed by PCA and thus capturing all the shape variations is still very challenging. As a result, the number of required datasets often seems insufficient to fully capture variations in both shape and appearance, especially in case of 3D modeling [START_REF] Heimann | Statistical shape models for 3d medical image segmentation: A review[END_REF]. Various works such as FEM vibrational modes [START_REF] Cootes | Data driven refinement of active shape model search[END_REF] have been proposed to artificially produce additional modes of variation, but it seems that this approach is mainly restricted to cope with intra-subject variability.

In this paper, we propose to study various shape and appearance priors in the context of lower limb bones segmentation using MR images and few training datasets. Two training sets are tested: one with only three datasets featuring a rather homogenous intensity distribution and the other with three more datasets featuring MRI artifacts. For shape modeling, both PCA and shape memory strategies are tested. PCA-based methods are known to need several datasets to be meaningful, while shape memory method requires in practice only one dataset. For appearance modeling, methods based on intensity profiles are tested, namely mean intensity profiles, multivariate Gaussian distributions of profiles and multimodal profiles from EM clustering. Our objective is to find the most efficient strategy, i.e. the strategy that is robust against the low number of datasets and their intensity inhomogeneities. This strategy would have the advantage to be more easily integrated in a clinical environment where the need of quick results is vital, regardless of the number of datasets.

MRI data

Acquisitions

For this study, six acquisitions were performed on six different subjects (4 females and 2 males, aged between 25 and 35). Protocols used for each acquisition are detailed in Table 1. The acquisitions took place in two different locations. Three subjects were scanned at the St Mary's Hospital, London, UK on a GE Medical Systems 1.5T MRI device (subjects #4, #5 and #6) and the three other subjects at the University Hospital of Geneva, Switzerland on a Philips Medical Systems 1.5T MRI device (subjects #1, #2 and #3). Both acquisitions feature 256x256 matrices but with a different slice thickness (i.e. 2 mm at London and 5 mm at Geneva). An institutional medical-ethical committee approved this study and subjects gave their written consent. Some MR images present strong artifacts like subject motion, surface coil, slice-to-slice interference and bias field (see Fig. 1). Subject motion is an artifact created by the displacement of a structure, which arises when the subject slightly moves between two consecutive acquisitions. Surface coil is characterized by a very strong signal due to the close proximity of the subject with the surface coil. Sliceto-slice interference artifact is due to the cross-excitation of adjacent slices with contrast loss in reconstructed images. Finally, bias field is a very common artifact in MR images, which may be induced by a number of factors such as poor radio frequency coil uniformity, static field inhomogeneity and radio frequency penetration.

Fusion of MRI slabs

An MRI acquisition of the whole lower limb cannot be performed in one scan due to a limited Field of View (FoV) of the machine. Consecutive scans, known as slabs, are necessary to cover its entire length. The registration of those slabs is then required to generate a complete MR image of the lower limbs. The registration is computed thanks to a sufficient overlap, which is used to compute the transformation matrix between slabs. In our case, we use a rigid registration based on the manual placement of landmarks. However, the registration is not enough to create a satisfactory MR image of the lower limbs. Indeed, intensity distribution may vary between registered slabs (see top of Fig. 2), which means the intensity histogram for the same structure may slightly differ between slabs. This is because the intensity range of the MR image may change between slabs (i.e. the minimum and maximum of intensity value is different). Also, the presence of strong artifacts at the image boundaries is likely to affect the intensity distribution (in addition to the common bias field). These artifacts are unpredictable and may differ in their number and importance between slabs.

Our solution consists in putting artifact intensity values into the background, so that when normalizing intensity between slabs, artifacts do not bias the correction. The normalization is performed using an intensity scaling factor. To this end, we calculate for each slab the histogram of the main structures (i.e. bones, muscles and fat) and compute the scaling factor so that histograms are similar. Finally, a bias correction is applied on the generic MR image of the lower limbs to correct intensity inhomogeneities [START_REF] Styner | Parametric estimate of intensity inhomogeneities applied to mri[END_REF], i.e. to correct the bias field and to smooth the intensity scaling between slabs (see bottom of Fig. 2).

Creation of appearance and shape priors

In this paper, we propose to exploit priors built from both the appearance and shape of structures in model-based image segmentation. A set of P training shapes {S 1 , . . . , S P } with corresponding images {I 1 , . . . , I P } is necessary to model the priors.

Each shape S is modeled as a 2-simplex mesh [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF] defined by N points x i with normals n i :

S = {(x 1 , n 1 ), . . . , (x N , n N )}
The shapes are produced by a supervised segmentation approach in which the point correspondence is established (i.e. landmarks on all P training shapes are located at corresponding positions). This point correspondence, which is the first necessary step when building shape models with PDM, is ensured through a registration between training shapes. Then, the P training shapes are aligned in a common coordinate frame. The most popular method to solve this problem is the General Procrustes Analysis (GPA) [START_REF] Goodall | Procrustes methods in the statistical analysis of shape[END_REF][START_REF] Gower | Generalized procrustes analysis[END_REF], which aligns the set of P training shapes to their unknown mean by minimizing the mean squared distance between two shapes in an iterative procedure. After alignment, dimensionality of the training set is reduced to find a small set of modes that best describes the observed variation. This is accomplished using PCA [START_REF] Jolliffe | Principal Component Analysis[END_REF]. For more information on the issues of building training sets, we refer the interested readers to the recent review of Heimann et al [START_REF] Heimann | Statistical shape models for 3d medical image segmentation: A review[END_REF].

Shape priors

For the shape prior construction, we propose to use both PCA-based priors (PCA) and shape memory (SMEM).

PCA

After the alignement of the P training shapes into a common coordinate frame with GPA, a Statistical Shape Model (SSM) is built [START_REF] Cootes | Active appearance models[END_REF].

An arbitrary shape S is approximated from the computed statistics by:

S ≈ T (S + Φ b)
where vector S is the mean shape, Φ is a matrix of M (M ≤ P ) principal modes (with respective variances λ i ), b is a vector of shape parameters and T denotes the alignment transform.

To estimate the unknown parameters b and T , an iterative procedure is used [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF]. To ensure the SSM specificity, two kinds of constraints are considered: hard or soft.

Hard constraint is defined as:

-3 λ m ≤ b m ≤ 3 λ m , ∀m ∈ [1, M ] Soft constraint scales b m to lie inside a hyperellipsoid: b 2 m /λ m ≤ C, ∀m ∈ [1, M ]
where C is computed from the χ 2 M distribution [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF].

When a shape is replaced by its closest shape counterpart, we qualify this replacement as PCA regularization.

SMEM

A single shape can also express a basic prior by using some of its geometrical properties. In case of a 2-simplex mesh, a local description of each vertex with respect to its three neighbors is computed [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF]. Only three independent simplex parameters are necessary, and those are similarity transform invariant. By using these prerecorded parameters, shape can be locally recovered. Although this prior is based on a single representative shape, the similarity invariance property confers more flexibility than encoding the shape as 3D points. This simple prior based on the simplex parameters is often denoted as shape memory and is similar to the notion of strain energy in mechanics.

Appearance priors

To take appearance into account, intensity profiles p i are built by sampling the image intensity at each point x i along the normal direction [START_REF] Cootes | Active appearance models[END_REF][START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF]. From these profiles, various approaches to build a prior have been presented in the literature. In this paper, we consider mean intensity profiles (PROF), multivariate Gaussian distributions of intensity profiles (MGD) and multimodal profiles (MPAM), which are built from an EM clustering of intensity profiles (see Fig. 3).

PROF

Mean intensity profiles constitute the simplest appearance prior. At each corresponding point through all datasets, a mean intensity profile is computed as:

µ i = P j=1 p i /P
Though faster to compute due to its simplicity, this prior does not make any assumption about the variance. As a result, mean intensity profiles are rather poor priors. However, they have been successfully exploited in previous works [START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF][START_REF] Schmid | Mri bone segmentation using deformable models and shape priors[END_REF] when combined with robust similarity measures such as the Normalized Cross Correlation (NCC) [START_REF] Holden | Voxel similarity measures for 3d serial mr brain image registration[END_REF].

MGD

Cootes et al. [START_REF] Cootes | Active appearance models[END_REF] proposed to model the prior of the (normalized) intensity profile at each corresponding point by a multivariate Gaussian of mean µ i and covariance matrix Σ i . This model has been successfully exploited in various works [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF][START_REF] Heimann | Statistical shape models for 3d medical image segmentation: A review[END_REF]. As µ i and Σ i are provided by the prior, the Mahalanobis distance is usually used as similarity measure.

However, such a distance needs the computation of the inverse of Σ i . This computation is known to be problematic, since the inversion of the covariance matrix may lead to singular matrices. This is especially the case with intensity profiles that feature a high dimensional space (i.e. profiles are large to efficiently describe appearance), but without enough profiles to represent those dimensions (so called curse of dimensionality). In our case, this effect is even more important as we study statistical models performance when using few training datasets.

In the literature, several techniques have been proposed to regularize the covariance matrix, and thus to avoid singularities [START_REF] Tadjudin | Covariance estimation with limited training samples[END_REF]. Cootes proposed an alternative approach [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF] consisting in first performing a PCA on Σ i , and then computing the Mahalanobis distance directly from the PCA statistics using the principal modes and associated eigenvalues. With such approach, no matrix inversion needs to be performed. In addition, a dimensionality reduction can be applied by selecting only the first modes, which leads to an approximated Mahalanobis distance.

MPAM

The clustering of intensity profiles is formulated in the context of a probability density estimation using Gaussian Mixture Models (GMM) [START_REF] Chung | Multimodal prior appearance models based on regional clustering of intensity profiles[END_REF]. Intensity profile classes are estimated using the Expectation-Maximization algorithm (EM) for each mesh but not for each point (i.e. without the need for any registration). EM is initialized with the Fuzzy C-Means algorithm (FCM), which is itself initialized with random cluster centers.

A large profile length (to efficiently describe appearance) combined with a coarse sampling of meshes (to speed up computation) is likely to lead to singular matrices when inverting EM covariance matrices during the E-step. To solve this problem, three distinct methods to regularize the covariance matrix are proposed: Spectral, Diagonal and Constant Regularization. To determine the number of classes that best represents the data, a novel non parametric model order selection criterion called Overlap Separation Index (OSI) inspired by cluster validity indices [START_REF] Kim | On cluster validity index for estimation of the optimal number of fuzzy clusters[END_REF] is used.

Spatially, the EM classification does not take into account the neighborhood information of intensity profiles. This leads to a non spatially smooth distribution of intensity profiles into the different clusters, which may impair the fusion of appearance regions. To take the connectivity between profiles into account, the Neighborhood EM algorithm (NEM) [START_REF] Ambroise | Clustering of spatial data by the em algorithm[END_REF] is used. Fig. 3 outlines an example of those different steps used for intensity profile clustering.

After the classification, each mesh M p is associated with K p clusters. MPAM are created by the fusion of these K p clusters that may vary among meshes. The K p clusters are compared between meshes and possibly merged. In order to merge similar clusters, the Jaccard index is used as a sim-ilarity measure. A threshold J between 0 and 1 is used to decide whether two clusters are equivalent.

The last step consists in providing a geometric embedding for the new clusters. For a set of meshes without point correspondence, each mesh M p needs to be registered to a reference mesh M , on which each posterior probability is resampled using a closest point approach [START_REF] Chung | Multimodal prior appearance models based on regional clustering of intensity profiles[END_REF]. However, since we use a set of P meshes M p with point correspondence, this registration is not necessary.

At the end, several clusters, or intensity profile priors, may be assigned at each vertex of reference mesh M , which leads to a multimodal approach. From an image and a mesh (a), a segmentation is performed (b). Intensity profiles are sampled (c), then classified (d) into classes (here, K={2,3,4}). Using an appropriate criterion, best classification is chosen (here, for K=4). Finally, the classification is spatially smoothed (e).

Segmentation based on priors

The proposed segmentation is based on dynamic deformable models [START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF][START_REF] Schmid | Mri bone segmentation using deformable models and shape priors[END_REF]. In such segmentation, a deformable template evolves until reaching an equilibrium. Internal forces regulate its evolution while external forces drive it towards anatomical boundaries. This section discusses the effect of abovementioned priors into this deformable model framework.

Evolution

Dynamic deformable models behave like a particle system, in which each particle corresponds to a lumped-mass vertex subject to internal and external forces. The dynamics of the system follow the Newtonian laws of motion, which express particle state (i.e. position and velocity) with respect to forces. The resulting time-discretized differential equation system is solved by an implicit Euler scheme.

A multiresolution scheme is used to produce various levels of detail (LOD) of the shapes [START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF] (see Fig. 5). The LOD are then exploited in a coarse-to-fine fashion, improving the robustness and accuracy of the segmentation evolution. In case of simultaneous segmentation of more than one structure, efficient collision response and detection are applied to prevent interpenetrations [START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF]. Alternatively, a postprocessing correction method can be used [START_REF] Schmid | Musculoskeletal simulation model generation from MRI datasets and motion capture data[END_REF]. Forces are expressed based on the image information, the current model configuration and the pre-computed priors.

Forces at point x i are expressed as the force of a Hookean spring:

f i = α * (x i -x i )
where xi denotes the target point and α is a weighting coefficient specific to each type of force. We will see in the following how the target point is computed given the different forces. This procedure is hereupon referred to as sourceto-target approach. Hard or soft constraints are applied to discard illegal configurations. Then, Ŝ = { x1 , . . . , xN } eventually becomes the target shape and the source-to-target approach is applied:

f pca i = α pca ( xi -x i )

SMEM

In case of a single reference shape used as a prior, both precomputed and current simplex parameters are used to derive new target point positions [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF][START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF]. A force f smem i is then produced at each point x i .

External forces based on appearance priors

At each iteration of the evolution, a number W of intensity profiles {p i 1 , . . . , p i W } are sampled along the normal n i at point x i . Among them, a target profile pi is chosen, whose corresponding position is xi . Usually, this target profile is chosen so that it maximizes a similarity criterion or minimizes a distance.

PROF

In case of mean intensity profiles, a target profile pi is selected if it maximizes the Normalized Cross Correlation N CC with the mean intensity profile µ i :

pi = argmax N CC(p i j , µ i ) where j ∈ [1, W ] 4.3.2 MGD
When using a multivariate Gaussian distribution of intensity profiles, the information from the covariance matrix Σ i is also taken into account. In this case, a target profile pi is selected if it minimizes the Mahalanobis distance d M derived from the computed distribution [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF]:

pi = argmin d M (p i j ) = (p i -µ i ) T Σ -1 i (p i -µ i ) where j ∈ [1, W ].
As previously said in section 3.2.2, a regularization of the covariance matrix Σ i is necessary to avoid any singularities [START_REF] Schäfer | A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics[END_REF]. However, the Mahalanobis distance d M may be computed without inverting the covariance matrix, i.e. by using an alternative approach proposed by Cootes [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF].

Given a PCA performed on P profiles p i at vertex x i , expressed by the principal matrix Φ i , the m eigenvalues λ il and the mean intensity profile µ i , Mahalanobis distance d M is defined as:

d M (p i ) =
pi = µ i + Φ i b i 4.3.3 MPAM
In case of multimodal profiles, the comparison with multiple intensity profile priors is considered. Indeed, since the prior is multimodal, more than one intensity profile prior may be assigned to x i . In practice, only intensity profile priors whose posterior probability is higher than a threshold (e.g. 10 -3 ) are considered. This can be done because EM classification leads to sharp posterior probabilities whose values are either very high or very low (i.e. quite below 10 -3 ).

Similarly to MGD, the Mahalanobis distance is used as similarity measure. But this time, a profile p i is compared with the K i intensity profile priors associated to x i (i.e. remember that each point x i is associated with a series of clusters, whose number may not be the same for every point).

Each point x i is associated with K i clusters of center µ k i and covariance matrix Σ k i .

The target profile pi of a given point x i is chosen as one of the W profiles sampled along the normal that has the smallest Mahalanobis distance d M with one of its K i associated clusters:

pi = argmin d M (p i j,k ) where j ∈ [1, W ] and k ∈ [1, K i ] .
Unlike MGD, only the diagonal terms {σ k i,1 , . . . , σ k i,L } of the L × L covariance matrix Σ k i are taken into account:

d M (p i k ) = (p i -µ i ) T diag(1/σ k i,1 , . . . , 1/σ k i,L )(p i -µ i )
where diag(1/σ k i,1 , . . . , 1/σ k i,L ) depicts a diagonal matrix.

Using only the diagonal terms has the advantage to speed up computation, which may be considerable as the comparison is now multimodal. This makes also sense since the covariance matrix regularization performed during MPAM creation strongly reduces the influence of non-diagonal elements.

Experimental setup

Two training sets, D3 and D6 with three and six lower limb MRI datasets respectively (see Fig. 4 for more details), are used with Leave-One-Out (LOO) cross validation. A multiresolution scheme is used to produce four levels of detail for each structure. For femur bone, the four resolutions consist in N = 514, N = 2056, N = 8224 and N = 32896 vertices respectively. For hip bone, they consist in N = 814, N = 3256, N = 13024 and N = 52096 vertices respectively (see Fig. 5).

Three appearance models are compared: mean intensity profiles (PROF) with Normalized Cross Correlation measure, multivariate Gaussian distributions of intensity profiles (MGD) and multimodal profiles (MPAM), both with Mahalanobis distance. With MGD, PCA takes 95% of the total intensity into account to compute the approximated Mahalanobis distance. To regularize the covariance matrix during EM iterations, MPAM is created using Constant Regularization method coupled with parameter h = 1.0 [START_REF] Chung | Multimodal prior appearance models based on regional clustering of intensity profiles[END_REF]. Due to the limited number of training datasets, we prefer not to merge any mode (i.e. Jaccard index J = 1.0).

Regarding intensity profiles, thirty-one intensities are sampled every 0.5 mm from 12.5 mm inside to 2.5 mm outside mesh surface at each vertex for PROF and MGD, since those values were successfully used for the segmentation of bones

#1 #4 #2 #5 #3 # 6
Fig. 4 The six datasets that are used for the experiments with their image and corresponding reference meshes (i.e. hip bones in blue and femur bones in green). The three datasets on the left constitute the D3 training set while all datasets together constitute the D6 training set. Note how rather linear is the intensity distribution between datasets of D3. Conversely, the three datasets on the right feature a quite heterogeneous intensity distribution. Among them, dataset #4 and dataset #6 are strongly subject to MRI artifacts.

in MR images [START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF][START_REF] Schmid | Mri bone segmentation using deformable models and shape priors[END_REF]. For MPAM, appropriate values need to be defined since the appearance prior construction method is different. Experiments showed that eleven intensities sampled every mm from 5 mm inside to 5 mm outside mesh surface give reasonable results in terms of clustering and segmentation. Those values are thus used in the framework of this work. To have a fair comparison between methods, same initialization and internal forces (i.e. PCA-based prior or shape memory) are used for every appearance model. The initialization is based on the manual placement of landmarks coupled with a shape interpolation approach [START_REF] Gilles | Musculoskeletal mri segmentation using multi-resolution simplex meshes with medial representations[END_REF]. Segmentation accuracy is assessed with DICE measure (DSC) [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF] based on reference manual segmentations. In the results, femur and hip bones are both considered as one structure, though they are both represented by two instances (i.e. left and right). To simplify the statistical analysis, the DICE measure on each structure is actually a mean on its both instances. Also, only datasets #1, #2 and #3 are segmented, since they feature a more homogeneous intensity distribution (see Fig. 4).

Results

First results consist in averaging the DICE measure on all structures and on all segmented datasets. Using D3 (see Table 2), PROF gives the best results, with both SMEM (DSC = 91.98) and PCA (DSC = 89.51), followed closely by MPAM with PCA (DSC = 89.40). Results are quite similar using D6 (see Table 3), PROF performs the best with both SMEM (DSC = 90.91) and PCA (DSC = 89.14), followed by MPAM with PCA (DSC = 88.08). MGD coupled with PCA gives by far the worst results, when using both D3 and D6. Regardless of the appearance prior used, SMEM is more accurate than PCA except when the latter is coupled with MPAM. When comparing global results on D3 and D6, we notice that D3 always gives better results.

In a second step, we study in more details the results on the two separate structures (see some delineations on Fig. Reference mesh is depicted in red, initialization in green, SMEM shape prior in blue and PCA shape prior in orange. and Fig. 7 for femur and hip bones respectively). Using D3 (see Table 4), PROF gives the best results with SMEM for femur bone (DSC = 92.51). This combination performs also the best for hip bone (DSC = 91.45). Using D6 (see Table 5), MPAM coupled with PCA gives the most accurate segmentation for femur bone (DSC = 92.02). However, for hip bone, PROF with SMEM gives once again the best results (DSC = 90.09). A clear tendency shows that results for femur bone are clearly better than for hip bone, regardless of the training set, shape and appearance prior used.

When comparing the influence of shape priors on femur bone segmentation with D6, we notice very slight differences between appearance priors, except for MPAM. SMEM performs the best for hip bone segmentation though (with both D3 and D6), but once again except for MPAM. When comparing the influence of the training set, D3 gives better results on hip bone, regardless of the shape and appearance prior used. Conversely, D6 gives better results on femur bone when PCA is used with all appearance priors.

Discussion

Generating statistics from a segmentation framework using a small amount of training datasets is a quite challenging work and a difficult task. Indeed, statistics are usually relevant when numerous data are at disposal. However, results on both training sets D3 and D6 show some trends that we discuss in this section. In general, using only three datasets with D3 gives better results than six datasets with D6. This may appear as counter-intuitive: more datasets should lead to better results. However, the three additional datasets from D6 are mostly corrupted by noise and artifacts. They also feature a quite heterogeneous intensity distribution. In this case, adding more information, which is corrupted, does not improve the appearance prior but rather weakens it. As a result, external forces have more difficulties to find the right boundaries.

In general, using shape memory with SMEM gives better results than PCA-based priors. This is because SMEM is a local approach. Given a shape that does not exhibit large inter-subject variation (i.e. like bones compared to soft organs) and that is sufficiently initialized, SMEM is a quite robust approach that can deal with few training datasets. However, PCA is known to give better results when using more datasets. In fact, a PCA-based evolution tends to be less sensitive to initialization when a sufficient number of datasets is available to provide a richer information about the shape [START_REF] Schmid | Extreme leg motion analysis of professional ballet dancers via mri segmentation of multiple leg postures[END_REF]. This is explained by the fact that PCA af-fects all the shape vertices in a global manner compared to the shape memory force, which confers a better robustness against local artifacts.

As depicted in both Table 4 and Table 5, MGD coupled with PCA gives very bad results when segmenting hip bone (DSC = 67.07 with D3 and 48.77 with D6). This would suggest that this combination is the worst: PCA as a shape prior (in case of few training datasets) and MGD as an appearance prior. External forces generated by MGD are thus less efficient than those generated by PROF and MPAM. This would explain the huge difference of DSC between D3 and D6, knowing that D6 contains datasets with images corrupted by noise and artifacts. As a result, those external forces are more sensitive to internal forces. But since internal forces based on PCA are also weak due to the few number of datasets, the segmentation is doomed to give bad results.

Regarding MPAM, the best results are for femur bone segmentation when combined with PCA and when using D6. Though in theory MPAM should need less datasets than MGD (i.e. a PCA-based appearance method), MPAM seems to feature a certain sensitivity to noise. Indeed, results for femur bone are clearly better than for hip bone (when using both D3 and D6). This is because intensity distribution is more heterogeneous at hip level, as hip bones are located close to image top border. This would suggest that the similarity measure (i.e. the Mahalanobis distance) should be optimized to cope with non linear intensity. Moreover, the intensity profile computation is different for MPAM than for PROF and MGD. Additional experiments should be performed to determine an optimal intensity profile length for the segmentation of bones in MR images.

In general, PROF gives the best results (i.e. compared to MGD and MPAM). We believe there are two explanations for that. First, reference meshes are produced by a supervised approach that uses a semi-automatic segmentation controlled by manually defined constraints [START_REF] Schmid | Extreme leg motion analysis of professional ballet dancers via mri segmentation of multiple leg postures[END_REF]. This semiautomatic segmentation exploits the same deformable model evolution coupled with the NCC similarity measure reported in this paper. This creates a bias, which is likely to give an advantage to PROF-based appearance prior (i.e. PROF also uses NCC during the segmentation). Second, NCC similarity measure is known to be more robust to intensity change (i.e. NCC is affine invariant). Using NCC would thus help PROF to be more efficient in presence of intensity change (i.e. especially with the three additional datasets corrupted by noise and artifacts that are used with D6).

Conclusion

When using a small amount of training datasets for a segmentation, results tend to show that local and simple methods perform the best. As a shape prior, shape memory (SMEM) gives very good results. As an appearance prior, mean intensity profiles (PROF) gives the best results. We believe these good results are partly due to the Normalized Cross Correlation (NCC) similarity measure, which is more robust to intensity change in MR images.

When increasing the number of training datasets, results tend to show that statistical methods feature the best results: PCA-based shape prior (PCA) and multimodal profiles (MPAM) as an appearance prior. Both methods capture more and more information while the number of training datasets is increased and we believe that better segmentation results would be obtained when increasing this number of datasets. However, this hypothesis only holds when the training data is of enough quality to produce meaningful and efficient priors.

Fig. 1

 1 Fig. 1 Some MRI artifacts: (a) surface coil artifact on three distinctive spots, (b) slice-to-slice interference artifact at hip level and (c) subject motion artifact at the overlap between two consecutive slabs where the femur bone is clearly shifted.

Fig. 2

 2 Fig. 2 MRI saggital view of the lower limbs before (top) and after (bottom) intensity correction. Note how the intensity has been corrected at foot level.

Fig. 3

 3 Fig. 3 Construction of the MPAM appearance prior of a femur bone.From an image and a mesh (a), a segmentation is performed (b). Intensity profiles are sampled (c), then classified (d) into classes (here, K={2,3,4}). Using an appropriate criterion, best classification is chosen (here, for K=4). Finally, the classification is spatially smoothed (e).

4. 2

 2 Internal forces based on shape priors 4.2.1 PCA As depicted in section 3.1.1, the shape priors are expressed by a SSM built on a PCA. At each iteration, a closest shape Ŝ is found by projecting the current deformable shape S into the PCA space. An iterative procedure computes the adequate transformation T and appropriate constrained shape parameters b = b 1 , . . . , b M [6].

  where b i = (b i1 , . . . , b im ) T is the model parameter vector of the best fit pi of p i given the PCA model:

Fig. 5

 5 Fig. 5 Detail of the four increasing resolutions (from left to right) for both femur (top) and hip (bottom) bones. For femur bone, the four resolutions are N = 514, N = 2056, N = 8224 and N = 32896 vertices respectively. For hip bone, the four resolutions are N = 814, N = 3256, N = 13024 and N = 52096 vertices respectively.

Fig. 6

 6 Fig. 6 Segmentation of left femur bone from dataset #1 using D3 training set. Results of PROF (left) and MPAM (right) appearance priors.Reference mesh is depicted in red, initialization in green, SMEM shape prior in blue and PCA shape prior in orange.

Fig. 7

 7 Fig. 7 Segmentation of left hip bone from dataset #1 using D3 training set. Results of MGD (left) and PROF (right) appearance priors. Reference mesh is depicted in red, initialization in green, SMEM shape prior in blue and PCA shape prior in orange.

Table 1

 1 MRI protocols used to scan the 6 subjects.

	Subject	TR/TE(ms)	FOV(cm)/Matrix	Resolution(mm)
	#1	4.15/1.69	35/256x256	1.37x1.37x5
	#2	4.15/1.69	35/256x256	1.37x1.37x5
	#3	4.15/1.69	35/256x256	1.37x1.37x5
	#4	5.06/2.23	43/256x256	0.84x0.84x2
	#5	4.34/1.56	40/256x256	0.78x0.78x2
	#6	5.09/2.22	43/256x256	0.84x0.84x2

Table 2

 2 Mean DICE measure (i.e. on all structures and on all segmented datasets) when combining appearance and shape priors on D3 training set.

		MGD PROF MPAM
	SMEM 88.48	91.98	87.39
	PCA	78.82	89.51	89.40

  [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF]

Table 3

 3 Mean DICE measure (i.e. on all structures and on all segmented datasets) when combining appearance and shape priors on D6 training set.

		MGD PROF MPAM
	SMEM 87.64	90.91	87.19
	PCA	70.37	89.14	88.08

Table 4

 4 Mean ± standard deviation of DICE measure for each structure and on all segmented datasets when combining appearance and shape priors on D3 training set. ± 0.63 92.51 ± 1.36 90.51 ± 0.51 PCA 90.57 ± 1.36 90.52 ± 0.91 90.71 ± 0.90 ± 9.84 91.45 ± 2.46 84.27 ± 6.58 PCA 67.07 ± 30.6 88.50 ± 2.07 88.10 ± 3.22

	Femur	MGD	PROF	MPAM
	SMEM 92.28 Hip MGD	PROF	MPAM
	SMEM 84.67

Table 5

 5 Mean ± standard deviation of DICE measure for each structure and on all segmented datasets when combining appearance and shape priors on D6 training set. ± 0.89 91.74 ± 1.72 90.46 ± 0.37 PCA 91.97 ± 1.42 91.62 ± 1.71 92.02 ± 0.96 ± 5.24 90.09 ± 2.22 83.92 ± 5.65 PCA 48.77 ± 33.1 86.66 ± 2.29 84.13 ± 7.20

	Femur	MGD	PROF	MPAM
	SMEM 91.92 Hip MGD	PROF	MPAM
	SMEM 83.36

Acknowledgment This work is supported by the 3D Anatomical Human project (MRTN-CT-2006-035763) funded by the European Union. The authors would like to thank the six subjects that took part in the MRI acquisitions and both St Mary's Hospital and University Hospital of Geneva. Special thanks to Mitchelle Chen and Prof. Andrew Todd-Pokropek from University College London for their help and collaboration during the MRI acquisitions in London.