Stanley Durrleman

Xavier Pennec

Alain Trouvé

Guido Gerig

Nicholas Ayache

Nicholas Ayache Spatiotemporal

Spatiotemporal Atlas Estimation for Developmental Delay Detection in Longitudinal Datasets

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Spatiotemporal Atlas Estimation for Developmental Delay Detection in Longitudinal Datasets

S. Durrleman 1,2 , X. Pennec 1 , A. Trouvé 2 , G. Gerig 3 , and N. Ayache 1 1 INRIA -Asclepios Team-Project, Sophia Antipolis, France 2 Centre de Mathématiques et Leurs Applications (CMLA), ENS-Cachan, France 3 Scientific Computing and Imaging (SCI) Institute, University of Utah, USA

Abstract. We propose a new methodology to analyze the anatomical variability of a set of longitudinal data (population scanned at several ages). This method accounts not only for the usual 3D anatomical variability (geometry of structures), but also for possible changes in the dynamics of evolution of the structures. It does not require that subjects are scanned the same number of times or at the same ages. First a regression model infers a continuous evolution of shapes from a set of observations of the same subject. Second, spatiotemporal registrations deform jointly [START_REF] Qiu | Parallel transport in diffeomorphisms distinguishes the time-dependent pattern of hippocampal surface deformation due to healthy aging and the dementia of the alzheimer's type[END_REF] the geometry of the evolving structure via 3D deformations and (2) the dynamics of evolution via time change functions. Third, we infer from a population a prototype scenario of evolution and its 4D variability. Our method is used to analyze the morphological evolution of 2D profiles of hominids skulls and to analyze brain growth from amygdala of autistics, developmental delay and control children.

Methodology for Statistics on Longitudinal Data

Many frameworks has been already proposed in medical imaging to analyze the anatomical variability of 3D structures like images, curves or surfaces. Less attention has been paid to the variability of longitudinal data (several subjects scanned several times). In [START_REF] Qiu | Parallel transport in diffeomorphisms distinguishes the time-dependent pattern of hippocampal surface deformation due to healthy aging and the dementia of the alzheimer's type[END_REF], the evolution between two shapes is modeled by a geodesic deformation, which cannot be used for more than two data per subjects. In [START_REF] Gorczowski | Statistical shape analysis of multi-object complexes[END_REF], shape growth is measured via the evolution of extracted features like volumes, shape or pose parameters. In [START_REF] Khan | Representation of time-varying shapes in the large deformation diffeomorphic framework[END_REF][START_REF] Davis | Population shape regression from random design data[END_REF], a temporal regression is proposed globally for a population, but this does not allow inter-subject comparisons. In cardiac motion analysis [START_REF] Chandrashekara | Construction of a statistical model for cardiac motion analysis using nonrigid image registration[END_REF][START_REF] Peyrat | Registration of 4D Time-Series of Cardiac Images with Multichannel Diffeomorphic Demons[END_REF], spatiotemporal registration relies on 3D-registrations between images of the same moment of the cardiac cycle and between two consecutive time-points. These works rely on time-point correspondence and do not call the labels of the time-points into question. By contrast, in longitudinal studies, subjects are scanned at ages which do not necessarily correspond. Moreover, evolutions may be delayed or advanced within a population, a key feature that we precisely aim at detecting. In [START_REF] Declerck | Definition of a 4D continuous planispheric transformation for the tracking and the analysis of LV motion[END_REF][START_REF] Perperidis | Spatio-temporal free-form registration of cardiac MRI sequences[END_REF], deformation of cardiac motion are proposed both in space and time but they require a fine temporal sampling of the motion, whereas only few acquisitions per subjects are available in most longitudinal studies. In this paper, we propose to use a regression model to estimate a continuous evolution from data sparsely distributed in time and spatiotemporal deformations which register jointly both the 3D geometry and the scenario of evolution. Geometrical data are modeled as currents to avoid assuming point correspondence between structures. Large deformations are used which gives a rigorous framework for statistics on deformations and atlas construction [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF][START_REF] Durrleman | Inferring brain variability from diffeomorphic deformations of currents: an integrative approach[END_REF][START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF]. From longitudinal data, we estimate consistently the most likely scenario of evolution and its spatiotemporal variability within the population.

In this paper, we call longitudinal data a set of geometrical data (curves or surfaces, called here shapes), acquired from different subjects scanned at several time-points. We assume that the successive data of a given subject are temporal samples of a continuous evolution. We propose therefore a regression model which computes a continuous evolution which matches the data of the subject at the corresponding time points (Fig. 1). This continuous evolution allows us to compare two subjects at a given age, even if one subject has not been scanned at this age. We can also analyze how the shape varies near this age to detect possible developmental delays. We define then the spatiotemporal deformation of a continuous evolution, which consists of two deformations: (1) a morphological deformation (of the 3D space) which changes the geometry of every frame of the evolution independently of the time point and (2) a time change function (deformation of the time interval) which changes the dynamics of the evolution without changing the geometry of shapes. To avoid time-reversal, the time change function must be smooth and order preserving: it is a diffeomorphism of the time interval of interest. A 4D registration between two subjects looks for the most regular spatiotemporal deformation, such that the deformation of the continuous evolution inferred from the first subject maps the successive target data (Fig. 2). Eventually, we use this 4D registration framework to estimate a spatiotemporal atlas from a population, based on an 4D extension of the statistical model of [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF]. We look for a template and a continuous evolution of this template (called mean scenario of evolution), so that data of each subject are temporal samples of a spatiotemporal deformation of the mean scenario. A Maximum A Posteriori estimation enables to estimate consistently the template, the mean scenario and the spatiotemporal deformations of this mean scenario to each subject.

We present the regression model in Sec. 2, 4D registration in Sec. 3 and atlas construction in Sec. 4. We focus here on the methodology and refer the reader to [START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF] for more computational details. In Sec. 5, we apply our method on 2D profiles of hominids skulls and amygdala of children scanned at 2 and 4 years.

Regression Model for Shape Evolution

We want to fit a continuous shape evolution to a set of shapes (S i ) of the same subject acquired at different time points (t i ). Without loss of generality, we can assume that t min = 0 and t max = T . This evolving shape is equal to the baseline M 0 at time t = 0, which may be the earliest shape of this subject or a template as in Sec. 4. The evolution has the form: M t = χ t (M 0 ) where t varies continuously in the time interval [0, T ]. For each t, χ t is a diffeomorphism of the 3D space, such that χ 0 = Id (which leads to χ 0 (M 0 ) = M 0 ). The regression (M t ) must match the observation S i at the time-points t i , while a rigidity constraint controls the regularity of the regression. This is achieved by minimizing:

J(χ) = ti d(χ ti (M 0 ), S i ) 2 + γ χ Reg(χ) ( 1
)
where d is a similarity measure between shapes, Reg(χ) a regularity term and γ χ a trade-off between regularity and fidelity to data. Among other possible choices, we use here the large deformations of [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF], and model curves or surfaces as currents [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]. Therefore, d is the distance between currents and χ is the solution of the partial differential equation : ∂χ t (x)/∂t = v χ t (χ t (x)) with initial condition χ 0 = Id. The regularity term Reg(χ) is given by T 0 v χ t 2 dt. In the setting of [START_REF] Vaillant | Surface matching via currents[END_REF], J depends on time-dependent momenta α k (t) at each samples of M t 's, which are used as variables for the gradient descent [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF][START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF].

As a result, for all t, χ t depends on all the constraints in the past and future. This differs from pairwise registrations between consecutive time-points. The function χ is piecewise geodesic and can be extended at all times by assuming v χ t = 0 (and hence χ t constant) outside [0, T ]. This is useful to compare this evolution with another subject which may have data outside the time interval.

Spatiotemporal Pairwise Registration

We assume now that we have successive shapes for the source subject (S(t i )) i and for the target (T (t j )) j . As in Sec 2, we perform a regression on the source shapes which leads to a continuous evolution S(t) t∈[0,T ] . Our goal is to find a diffeomorphism of the 3D space φ and a diffeomorphism of the time-interval ψ which deform the source evolution S(t) into S (t) = φ(S(ψ(t))) such that S (t j ) match T (t j ). Thanks to the regression function, no correspondence is needed between the time points t i and t j . Formally, we minimize:

J(φ, ψ) = tj d(φ(S(ψ(t j ))), T tj ) 2 + γ φ Reg(φ) + γ ψ Reg(ψ) (2)
The spatial (resp. temporal) deformation φ (resp. ψ) is solution at parameter u = 1 of the flow equation

∂ u φ u (x) = v φ u (φ u (x)) (resp. ∂ u ψ u (t) = v ψ u (ψ(t))
). The norm of the speed vector fields v φ u and v ψ u integrated for u ∈ [0, 1] defines the regularity terms Reg(φ) and Reg(ψ) respectively. As in Sec 2, the geometrical (resp. temporal) deformation is parametrized by momenta α (resp. β) at the points of S(t j ) (resp. at the t j ), which are used as variables for the gradient descent. ∇ α J(α, β) is computed as for the registration of the collection of φ(S(ψ(t j ))) to the (T tj ) [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]. To compute ∇ β J, we need to derive the matching term with respect to the time-points t j = ψ(t j ) [START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF]. We use here a centered Euler scheme: ∇ t φ(S(t)) -T 2 ∼ φ(S(t)) -T, φ(S(t + δt)) -φ(S(t -δt)) /δt, where δt is the discrete time step. We used here the fact that the shapes are embedded with a vector space (the space of currents) provided with an inner product. Note that we minimize J with respect to the geometrical and the temporal parameters jointly. We do not performed alternated minimization.

Spatiotemporal Atlas Construction

We assume now that we have a set of N subjects (S i ), provided each with temporal observations (S i (t i j )) j . We are looking for a template M 0 and a mean scenario of evolution of this template M (t) = χ t (M 0 ), such that the observations correspond to particular moments of a spatiotemporal deformation of the mean scenario. This means that φ i (M (ψ i (t i j ))) match S i (t i j ) for each subject i and time t i j . Maximum A Posteriori estimation in the same setting as in [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF], leads to the minimization of

J(ψ i , φ i , χ, M 0 ) = N i=1    t i j d(φ i (χ ψ i (t i j ) M 0 ), S i (t i j )) 2 + γ φ Reg(φ i ) + γ ψ Reg(ψ i ) + γ χ Reg(χ)   
We perform a 3 step alternated minimization. If the template M 0 and the regression χ are fixed, the minimum is achieved for N registrations of the mean scenario χ t (M 0 ) to each subject's set of data S i (t i j ), as in Sec. 3. If we fix the N spatiotemporal deformations (φ i , ψ i ) and the regression χ, we need to minimize i,j d(Φ i,j (M 0 ), S i (t i j )) 2 , where the Φ i,j = φ i • χ ψ i (t i j ) are 3Ddiffeomorphisms. This is exactly the estimation of an unbiased template in the setting of [START_REF] Durrleman | Statistical models of sets of curves and surfaces based on currents[END_REF], when the deformations are given by Φ i,j . When the template M 0 and the N spatiotemporal deformations (φ i , ψ i ) are fixed, we need to minimize i,j d(φ i (χ ψ i (t i j ) M 0 ), S i (t i j )) 2 + γ χ Reg(χ). This is not the regression problem stated in Sec. 2 because of the deformation φ i in the matching term. To turn it into regression, we approximate the matching term d(φ i (χ

ψ i (t i j ) M 0 ), S i (t i j )) by d(χ ψ i (t i j ) (M 0 ), (φ i ) -1 (S i (t i j )
)) (subject's shapes are matched back to the mean anatomy). This approximation is valid only for diffeomorphisms φ i whose Jacobian is close to the identity. To initialize the minimization, we set M 0 as the mean current of the earliest data and set the diffeomorphisms χ, φ i , ψ i to identity.

Numerical Experiments

Experiments on 2D curves relates to Sec.2 and 3, those on surfaces to Sec. 4.

Evolution of 2D Curves

We have five 2D-profiles of hominids skulls which consist of six lines each (source: www.bordalierinstitute.com). Our regression framework infers a continuous evolution from the Australopithecus to the Homo sapiens sapiens which matches the intermediate stages of evolution in Fig. 1.

Then, we register the evolution {Homo habilis-Homo erectus-Homo neander-talensis} to the evolution {Homo erectus-Homo sapiens sapiens} in Fig 2 . The geometrical deformation shows that during the later evolution the jaw was less prominent and the skull larger and rounder than during the earlier evolution. The time change function shows that the later evolution occurs at a speed 1.66 times faster than the earlier evolution. This value is compatible with the growth speed of the skull during these periods (See Fig. 3): between Homo erectus and sapiens the skull volume growths at (1500 -900)/0.7 = 860cm 3 /10 6 years, whereas between Homo habilis and neandertalensis, it growths at (1500 -600)/1.7 = 530cm 3 /10 6 years, namely 1.62 times faster.

Evolution of 3D surfaces

We use here meshes of amygdala of the right hemisphere from 4 autistics, 4 developmental delay and 4 control children scanned twice [START_REF] Hazlett | Magnetic resonance imaging and head circumference study of brain size in autism[END_REF]. Age distribution is shown in Fig 5-a. From these data registered rigidly, we infer a template, a mean scenario of evolution of this template and the spatiotemporal evolution of this mean scenario to each subject. In the setting of [START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF][START_REF] Durrleman | Inferring brain variability from diffeomorphic deformations of currents: an integrative approach[END_REF], the diffeomorphisms are controlled by the standard deviation of Gaussian kernel set to 15 mm for χ, φ and 1 year for ψ; the typical scale on currents is set to 3 mm. Trade-off γ χ , γ φ were set to 10 -3 and γ ψ to 10 -6 . An amygdala is typically 10 mm large. The discrete time step is set to 0.2 years.

By inspection of the companion movie of the mean scenario, one distinguishes 4 phases during growth (See also Fig. 4). Preliminary tests do not show correlations between the morphological deformations and the pathology. From the time change functions shown in Fig. 5, we cannot conclude that a subject with pathology is systematically delayed or advanced compared to controls, even at a given age. However, the curves show that the growth speed seems to follow the same pattern, mainly an acceleration between age 2.5 and 3.5 for the autistics and between age 4 and age 5 for controls. The developmental delay also have such pattern but it occurs at a very variable age. These results suggest that the discriminative information between classes might not be inferred from the anatomical variability at a given age, but rather from variations of the growth process. These results, however, must be strengthen using larger database. The more time-points per subjects, the more constrained the mean scenario estimation. The more subjects, the more robust the statistics. Arrows measures the differences between age t+0.2 and age t in the space of currents as in [START_REF] Durrleman | Sparse approximation of currents for statistics on curves and surfaces[END_REF]. From age 2 to 2.8, the evolution is mainly a torque at the posterior part; then the structure becomes thicker, mostly at the superior part between age 2.8 and 4 and at the inferior between age 4 and 6; from age 6 the evolution is a mainly a torque at the anterior part.

Discussion and Conclusion

In this paper, we present a generic framework to analyze variability of longitudinal data. A regression model fits a continuous evolution to successive data of one subject. 4D registrations decompose the difference between two sets of longitudinal data into a geometrical deformation and a change of the dynamics of evolution. The more acquisitions per subjects, the more constrained this decomposition. However, no constraint is imposed in terms of number and correspondence of measurement points across subjects. These pairwise registrations are used for group-wise statistics: ones estimates consistently a template, the mean evolution of this template and the spatiotemporal variability of this evolution in the population. Then, statistical measures can be derived, like the first mode of temporal deformation in Fig. 5. Further experiments have still to be performed to give more quantitative measures of variability. However, these first results suggest that pathologies might be characterized more by a particular scenario of evolution than by the anatomy at a given age. Our methodology can be used therefore to drive the search of new anatomical knowledge and to give characterization of pathologies in terms of organ growth scenario. This may be applied to the study of degenerative diseases or cardiac motion disorders.

Fig. 1 .

 1 Fig. 1. Skull profile of five hominids (in red). The regression model estimates a continuous evolution (in blue) of the Australopithecus, which closely matches the data.

Fig. 2 .

 2 Fig. 2. Registration of the evolution {Homo habilis-erectus-neandertalensis} (in red) to the evolution {Homo erectus-sapiens sapiens} (in green), shifted to start at the same time. Top row: Regression of the source data (red) gives the continuous evolution in blue. Middle row: The geometrical part φ is applied to each blue frame. This shows morphological changes: the skull is larger, rounder and the jaw less prominent. Bottom row: The time change function ψ is applied to the evolution of the second row. The blue shapes are moved along the time axis (as shown by dashed black lines), but they are not deformed. Black arrows that a better alignment is achieved when one accounts both for morphological changes and a change of the evolution speed.

Fig. 4 .

 4 Fig.[START_REF] Davis | Population shape regression from random design data[END_REF]. Mean Scenario of the right Amygdala (right lateral part). Arrows measures the differences between age t+0.2 and age t in the space of currents as in[START_REF] Durrleman | Sparse approximation of currents for statistics on curves and surfaces[END_REF]. From age 2 to 2.8, the evolution is mainly a torque at the posterior part; then the structure becomes thicker, mostly at the superior part between age 2.8 and 4 and at the inferior between age 4 and 6; from age 6 the evolution is a mainly a torque at the anterior part.

Fig. 5 .

 5 Fig. 5. Temporal deformation of the mean scenario Left: distribution of original (top) and registered (bottom) ages. Middle: time change functions for the 12 subjects. Right: First mode of variation at ±σ of the time change functions for each class. Autistics and controls show the same evolution pattern, but shifted in time.

Acknowledgments This work was partially supported the European IP project Health-e-child (IST-2004-027749) and Microsoft Research.


Florence Billet
email: florence.billet@inria.fr

Maxime Sermesant

Hervé Delingette

Nicholas Ayache

Cardiac Motion Recovery and Boundary Conditions Estimation by Coupling an Electromechanical Model and Cine-MRI Data

We present a method for cardiac motion recovery using the adjustment of an electromechanical model of the heart to cine MRI. This approach is based on a proactive model which consists in a constrained minimisation of an energy coupling the model and the data. The presented method relies on specific image features in order to constrain the motion of the endocardia and epicardium and impose boundary conditions at the base. Thus, image intensity and gradient information are used to constrain the motion of the myocardium surfaces while a 3D block matching technique leads to the motion estimation of base vertices. Finally, we show that the implicit time integration of those forces and personalised boundary conditions lead to a better cardiac motion recovery from cine-MR images.

Introduction

Cardiovascular diseases are the number one cause of death in the industrialised world and is projected to remain so. That's why the modeling of the heart's electromechanical activity is an active research area [START_REF] Hunter | Computational Electromechanics of the Heart[END_REF][START_REF] Sachse | Electro-mechanical modeling of the myocardium: Coupling and feedback mechanisms[END_REF][START_REF] Frangi | Three-Dimensional Modeling for Functional Analysis of Cardiac Images: A Review[END_REF]. Moreover, the coupling of these models and clinical data allows us to personalise them and then to simulate a number of pathologies or the effect of therapeutic actions, and to analyse the cardiac function of the considered patient.

There exist several approaches in the literature to create such patient-specific models of the heart from clinical images. The concept of deformable models is often used to adjust a geometrical model of the heart on time series of medical images [START_REF] Papademetris | 3D Cardiac Deformation from Ultrasound Images[END_REF][START_REF] Montagnat | 4D Deformable Models with temporal constraints: application to 4D cardiac image segmentation[END_REF][START_REF] Qian | A Learning Framework for the Automatic and Accurate Segmentation of Cardiac Tagged MRI Images[END_REF]. In such approaches, the surface or the volumetric mesh is fitted to the apparent boundaries of the heart in the image by minimising the sum of an image energy and a regularising or internal energy. Those models may produce proper segmentations and are rather efficient, however they rely on little a priori knowledge of the heart motion which limits their ability to recover tangential motion due to the aperture problem.

Other approaches strongly constrain the estimation of motion by coupling an electromechanical heart model using data assimilation methods. In such a framework, dynamical systems representing electromechanical models of the heart are adjusted to time series of images by adding a filter that constrain the state variables (i.e the position and the velocity) and the model parameters to match the available observations. P.C Shi et al. introduced such techniques by integrating cardiac models and Kalman filters [START_REF] Shi | Stochastic finite element framework for simultaneous estimation of cardiac kinematic functions and material parameters[END_REF][START_REF] Wong | Physiome Model Based State-Space Framework for Cardiac Kinematics Recovery[END_REF]. However, the use of extended Kalman filtering is often limited since they implies the use of full covariance matrices whose size are equal to the square of the number of state variables augmented by the number of parameters to estimate. To avoid the curse of dimensionality, some authors [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] have devised non optimal filters which do not involve any matrix inversion (unlike Kalman filtering). This approach leads to much faster computations and has been the inspiration of the work presented in this paper.

We showed in [START_REF] Billet | Cardiac Motion Recovery by Coupling an Electromechanical Model and Cine-MRI Data: First Steps[END_REF] the equivalence under some hypothesis between the method proposed in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF] and the use of a proactive deformable model introduced in [START_REF] Sermesant | An Electromechanical Model of the Heart for Image Analysis and Simulation[END_REF].

In this paper, we use the proactive deformable model which extends the concept of deformable models by adding an active electromechanical model of the heart motion as the regularising energy term. In this paper, we improve from our previous work with two additional contributions. First, we provide an implicit formulation of image forces that can better predict the image forces between two time steps. Second, we estimate the motion of the base of the ventricles, where the ventricles are connected to the atria and outflow tracts. This represents a significant advance since most previous work on electromechanical modelling of the heart consider that ventricle bases are fixed or linked to fixed points through springs. To recover the base motion from image, we use a block matching technique at the base points that uses intensity correlation to track points along the cardiac cycle.

Electromechanical Model

As we want the complexity of the model to match the relatively sparse measures available from imaging data, we consider in this paper a fairly reduced electromechanical model. For instance, our model relies on linear elasticity for myocardium mechanics although a nonlinear material undergoing large strain would be better suited. However, we found the current reduced model to be a good trade-off in terms of computational efficiency and physiological realism. Furthermore, this coarse level of modeling makes the estimation of the mechanical state and parameters tractable and allows us to test the behaviour of the model on several heart beats. Our cardiac model couples two biophysical components : electrophysiology and biomechanics. Due to its efficiency, we use a multi-front anisotropic Eikonal approach solved with a Fast Marching algorithm on a volumetric tetrahedral mesh [START_REF] Sermesant | An anisotropic multi-front fast marching method for real-time simulation of cardiac electrophysiology[END_REF].

The depolarisation time t d of the electrical wave for a given vertex of the volumetric mesh is computed by solving the anisotropic Eikonal equation v 2 (∇t T d D∇t d ) = 1, where v is the local conduction velocity parameter and D is the tensor defining the conduction anisotropy. In the fibre coordinates, D = diag(1, ρ, ρ), where ρ is the conduction anisotropy ratio between longitudinal and transverse directions.

The biomechanical model, derived from a multi-scale modelling of the myocardium detailed in [START_REF] Bestel | A Biomechanical Model of Muscle Contraction[END_REF], is composed of an active contractile element controlled by the electrophysiology and a parallel one which represents the properties of the tissue and is anisotropic linear visco-elastic.

The active contractile element is controlled by the depolarisation time t d computed with the Eikonal equation and creates a stress tensor σ C f ⊗ f where f is the 3D fibre orientation, ⊗ the dyadic product and σ C the strength of the contraction defined by the equation 1:

σ C (t) = σ 0 1 -e k AT P (t d -t) during depolarisation t d ≤ t < t r σ C (t r )e k RS (tr-t) during repolarisation t r ≤ t < t d + HP (1)
where t r = t d + AP D is the repolarisation time with AP D the active potential duration, σ 0 the maximum of the contraction, k AT P and k RS the contraction and relaxation rates respectively and HP the heart period. This tensor results in a 3D contraction force f C for each vertex of each element. Finally, the simplified dynamic law is:

M Ÿ + C Ẏ + KY = F P V + F C + F B
where M, C are the mass, Rayleigh damping matrices, K the stiffness matrix, F C , F P V , F B the force vectors respectively created by the active contractile elements, corresponding to the pressure forces in the ventricles and corresponding to other boundary conditions and Y , Ẏ = dY /dt, Ÿ = d 2 Y /dt 2 respectively the position, velocity and acceleration.

The four cardiac cycles (filling, isovolumetric contraction, ejection, isovolumetric relaxation) were simulated as presented in [START_REF] Sermesant | An Electromechanical Model of the Heart for Image Analysis and Simulation[END_REF]. We compute the arterial pressures using a 3-element Windkessel model described in [START_REF] Stergiopulos | Total arterial inertance as the fourth element of the windkessel model[END_REF].

The pressure forces are computed differently in the cardiac phases. In the filling phase, they are set to the values of the pressures in the atrias. In the ejection phase, they are set to the arterial pressures computed with the 3-element Windkessel model. In the isovolumetric phases, the ventricular pressures are computed in order to counterbalance the contraction forces and to keep the ventricle volumes constant.

Cine MR Image Processing

Cine-MRI data consists in a sequence of 20 to 30 high resolution MR images of the heart that describes a whole cardiac cycle. We build a myocardium mesh that fits the mid-diastolic MRI image of the sequence when the filling of the ventricles is almost finished, just before the atrial contraction. As described in [START_REF] Billet | Cardiac Motion Recovery by Coupling an Electromechanical Model and Cine-MRI Data: First Steps[END_REF], we select this mid-diastole image based on blood pool volumes curves. Then, an interactive tool is used to delineate the epicardium and left and right ventricles endocardia in order to create three binary masks of the epicardium and the endocardia. We combined these masks to obtain the binary mask of the myocardium and create the mesh with an isosurface extraction and tetrahedral mesh generation.

Furthermore, synthetic fibre orientations are computed by linearly interpolating the elevation angle between the fibre and the short axis plane, from 80 o on the endocardium to -80 o on the epicardium.

Proactive Deformable Model

Image Forces

In deformable model approaches, the optimal shape corresponds to the minimization of the sum of an image energy and a regularizing energy. The use of an electromechanical model as a regularising term allows us to introduce some a priori knowledge about the active cardiac motion and thus to recover some movement (like the tangential motion) which cannot be obtained from classical geometrical tracking methods. Thus, we add to the previous dynamic equation

M Ÿ + C Ẏ + KY = F P V + F C + F B an image force F img that derives from an image energy.
Of course, the image forces have no physiological meaning, and then the estimated motion is not physiological. But by coupling the model and the data and by estimating the model parameters (which is the next step of this work), the motion generated by the model should converge to the one observed in the images. Thus, the intensity of image forces should decrease along the estimation and the estimated motion should be more and more physiological.

The image energy is a metric that compares the simulated motion and the observed one. Thus, for each vertex Y i of the epicardium or endocardium surfaces, if we define as Z i the closest point on the apparent image boundary along the normal direction [START_REF] Billet | Cardiac Motion Recovery by Coupling an Electromechanical Model and Cine-MRI Data: First Steps[END_REF][START_REF] Sermesant | An Electromechanical Model of the Heart for Image Analysis and Simulation[END_REF]). In this equation γ i is the image gain that depends on the confidence we have in the corresponding image point Z i . In fact, with high gains, the estimated state will rely more on image data information than in the electromechanical model. Conversely, with gains equal to 0, the estimation does not take into account the data and is equivalent to the electromechanical model alone. Thus, the choice of the gains depends on the relative confidence in the model and the data. Furthermore, we showed in [START_REF] Billet | Cardiac Motion Recovery by Coupling an Electromechanical Model and Cine-MRI Data: First Steps[END_REF] how to set these gains in order to ensure convergence. For each vertex Y i of the epicardium or endocardium surfaces, the definition of Z i is a function of the image intensity and image gradient. More precisely, as shown in Fig. 1 (a) and (b), we first extract the intensity profile at Y i along N i and then determine the closest voxel Z i whose image gradient norm is greater than a threshold and whose intensity is within a given range value. Those threshold and range value are set on the first image and kept constant during the cardiac cycle. In Fig. 1 (a) and (b), the way to find the point Z i is shown in 2D for more simplicity, but it is performed using the all 3D image.

N i at Y i , then the image force at Y i is F img = 2γ i (Z i -Y i ) (see

Implicit Image Forces

The image forces then may be temporally integrated either with an explicit or implicit scheme. In [START_REF] Billet | Cardiac Motion Recovery by Coupling an Electromechanical Model and Cine-MRI Data: First Steps[END_REF], we used an explicit scheme for time integration. Thus, at each time step t, we solved the approximated equilibrium equation KY t+dt = F + F t img to compute the position Y t+dt of the mesh at the following time step, where K is the generalised stiffness matrix, F the external physical forces and F t img the image force as 2γ i (Z t i -Y t i ) if Y i belongs to the heart surfaces and 0 otherwise. As we know the position of the mesh at time t, we could find the point Z t i as explained above. This explicit scheme is computationally efficient since the generalised stiffness matrix does not change other time. However with such scheme, the proactive deformable model tends to lag behind the image observation which introduces a bias especially for large time steps.

Therefore, in this paper, we introduce an implicit scheme solving KY t+dt = F + F t+dt img , which requires predicting the next image force F t+dt img . To this end, we use a Taylor expansion of the image point Z t+dt i :

Z t+dt i = Z t i + ∂Z i ∂Y i t (Y t+dt i -Y t i ) = Z t i + B(Y t+dt i -Y t i ) (2)
The B matrix is a 3×3 matrix that indicates how point Z i moves as a function of Y i . We provide here a first order approximation of this matrix by discarding the effect of surface curvature. Thus, we consider that Y i lies on a plane with normal N i while Z i lies on a plane oriented by normal N img i . With those hypotheses, upon a displacement dY i , Z i is determined as the intersection of the line passing through Y i + dY i and directed by N i through the plane P img i . Obviously, only components of dY i orthogonal to N i entails a change in the point Z i . Similarly, only components orthogonal to N img i would entail a change in point Z i . Thus, skipping a detailed proof, one can show that the matrix B can be written as :

B = 1 cos 2 β (I 3×3 -N img i ⊗ N img i )(I 3×3 -N i ⊗ N i ) (3)
where β is the angle between N i and N img i , I 3×3 is the identity matrix and a ⊗ b = ab T is the dyadic product between the vectors a and b.

Finally, we need to solve the following equation: ( K + K img )Y t+dt = F + Fimg with Fimg a vector of size 3 × N whose 3D component f img,i is equal to 2γ i (Z t i -Y t i ) + 2γ i (I 3×3 -B)Y t i if Y i belongs to the myocardium surfaces, and 0 otherwise and with K img a block diagonal matrix whose 3 × 3 component (K img ) i is equal to 2γ i (I 3×3 -B) if Y i belongs to the myocardium surfaces, and the null matrix otherwise.

Boundary Constraints Estimation

Due to the modeling of the electromechanical behaviour of the heart, we have a lot of a priori knowledge on the motion of the ventricles. But for the base of the ventricles, near the atria, we have little a priori knowledge of the motion, and boundary constraints are set simply by attaching base vertices to fixed springs. Thus, we introduce here a way to estimate these boundary constraints by using a more complex external energy.

We cannot rely on intensity profiles to find Z i as described in the previous section. Indeed this would only estimate displacement of base vertices along its normal direction which is ill-defined and physically incorrect. Instead, we use a 3D search method based on 3D block matching algorithm, as in [START_REF] Ourselin | Block Matching: A General Framework to Improve Robustness of Rigid Registration of Medical Images[END_REF]. We hope therefore to better recover the real motion of the base, and not only its apparent one. Finally, we can use the similarity measure between the image blocks as a way to control the confidence in the point Z i .

Therefore, we attach to each base vertex of the initial mesh a surrounding block of voxels from the reference image, which is the first image of the simulation (see section 3). Then, we search at each time step the most similar block around the current position in the following image of the cardiac cycle (the target image) (see Fig. 1 (c) and(d)). This is based on the assumption that the intensity around a given vertex of the base of the mesh should be similar in the two images, the difference coming from noise and from the movement of the vertex between the two images, which is not very important at the base of the heart, near the valve. The matching voxel selected is the center of the block in the target image maximising the chosen similarity criterion. Then, in order to take into account that the motion at the base is mainly a vertical one, we project the vector connecting the current point and the marching voxel on the normal of the mesh at this current point. The external force is then proportional to the norm of this projection and directed from the normal of the mesh at the considered point. The search range and the search step can be changed in a coarse to fine approach, as in [START_REF] Ourselin | Block Matching: A General Framework to Improve Robustness of Rigid Registration of Medical Images[END_REF]. In this paper, we used a similarity based on the sum of square differences, but several other similarity measures can be used depending of the noise and the used image modalities. Moreover, these estimated boundary constraints could be set as inputs in further simulations with the personalised model. 6 Results

Validation on Synthetic Data

We generated synthetic cine-MR images using the electromechanical model with standard values. We considered 29 instants of the second simulated cycle and we generated the corresponding segmented 3D images, using tetrahedra rasterisation. We use for the different state estimation a cardiac model with same parameters as the model used to generate the synthetic data except for the reference position and the damping value. We choose to estimate the state using a low damping because the effect of the implicit forces is more visible with low damped systems. Indeed in a highly damped system all applied forces are smoothed to a large extent which entails a delay in the reaction of the cardiac model to new image information.

We then compare the state vector estimation (position and velocity) with no image forces and with an explicit and implicit time integration of image forces. The gain γ of the image forces is set to the same value for the two estimations. For each method, we computed the root mean squared error (RMSE) between the estimated positions and the positions of the reference simulation (see Fig. 2 (a)). The RMSE with implicit forces is lower than the one with explicit forces (15% lower), and much lower than the one computed without image forces. Indeed, the average of all RMSE without image forces, with explicit and implicit image forces are respectively 14.9mm, 4.5mm and 3.6mm. Furthermore, explicit integration with low damping leads to large oscillation in the right ventricle (see Fig. 2 (b)) due to undamped high frequencies in image forces. Finally, we can see that the mesh obtained with implicit forces better fits the segmented contour than the mesh obtained with explicit forces.

Results on Clinical Data

We used cine-MRI data of a patient with a left bundle branch block (LBBB) to test the effect of the block matching forces. The sequence consists in 30 images that describe a total heart cycle. The resolution of this dataset is 1.56 × 1.56 × 10mm 3 . In order to increase the robustness of the algorithms on clinical routine cardiac MRI, images are pre-processed. First, images are trilinearly resampled to get isometric voxel spacing (and then the resolution becomes 1.56 × 1.56 × 1.56mm 3 ). Next, the contrast is enhanced by removing the extrema quantiles of the grey-level histogram. Finally, the histogram is linearly rescaled between 0 and 255 to have consistent intensities between images. We manually adjusted the electrical parameters using the ECG and, in this case, with some clinically measured propagation times of the patient. If we have no measures of this propagation times, we initialize the electrical model with the volumes curves extracted from the images, as shown in [START_REF] Billet | Cardiac Motion Recovery by Coupling an Electromechanical Model and Cine-MRI Data: First Steps[END_REF]. The passive mechanical parameters used are taken from the literature [START_REF] Shi | Stochastic finite element framework for simultaneous estimation of cardiac kinematic functions and material parameters[END_REF]. In order to calibrate the active mechanical parameters, we use the volume curve to compute the ejection fractions, wich are closely related to these parameters. For our data, σ 0 was set to 0.09 MPa/mm 2 for the left ventricle, 0.07 MPa for the right ventricle and 0.06 for the scars. The cardiac cycle duration is 1.035s.

We did three different estimations, all using the implicit forces described above for the epicardium and the endocardium surfaces. The first estimation uses no image forces applied to the base and no boundary constraints. The second one uses no image forces applied to the base, but applies boundary constraints involving fixed springs. The last one applies block matching forces on the base and image forces in the endocardium and epicardium surfaces.

Fig. 3 shows some short and long axis views with the contours of the three estimated positions superimposed with the corresponding MR image.

With no boundary constraints and no image forces there is no attach points in the model of the heart which entails that the generalised stiffness matrix may become singular. Also we have no a priori knowledge of the base motion. Thus, large errors in the estimation may occur, as shown in Fig. 3 (d, e andf) (the motion of left part of the left ventricle and of the top of the right ventricle is not well recovered).

With boundary constraints on the base, the model may be too much constrained, specially during the contraction as shown in Fig. 3 (c andd). For instance, the motions of the septum (see the short axis views of Fig. 3 (c andd)) and at the top of the lateral part of the right ventricle (see the long axis views of Fig. 3 (c andd)) are not well recovered.

Finally, the block matching forces achieve a good compromise between these two approaches. The motion of the base vertices are determined by block matching instead of being fixed with springs. We can see that the right ventricle motion has been well recovered, even during the contraction of the heart, as shown in Fig. 3 (b, c andd). However, we can see that the motion of the top of the septum is not well recovered, due to a not so good intitial segmentation of the top of the septum. This segmentation error makes the initialisation of the block matching algorithm not good, and then, the founded appariements try to keep this error. We can see here the need of a good segmentation for the use of the block matching approach.

Conclusion

We presented a method which relies on specific image features to recover the motion of the myocardium and estimate boundary constraints at the base. The use of implicit time integration of image forces allows us to ensure a better accuracy and smoother motion tracking for low damped systems. Furthermore, the use of a 3D block matching algorithm for the estimation of the boundary constraints at the base appears to be a good trade off between cardiac motion recovery with no boundary constraints and with classical boundary constraints consisting of fixed vertices or fixed springs and set manually. Finally, these promising results must be extended in order to perform parameter estimation of the model which would lead to a patient-specific electromechanical model of the heart.

Fig. 1 .

 1 Fig. 1. (a, b) Surface image forces and corresponding Zi point in the image. (c,d) Block matching algorithm (c) Reference image: initialisation in the reference image of the 3D block. (d) Target image: most similar block search around the current position.

Fig. 2 .

 2 Fig. 2. (a) RMSE between the positions of the reference simulation and the positions estimated with implicit forces (solid) and explicit forces (dashed). (b) Image and superimposed meshes estimated with implicit (blue) and explicit (yellow) image forces.

Fig. 3 .

 3 Fig. 3. MR image short axis (top) and long axis (bottom) slices at different times of the heart cycle, with the contrours of the estimated model position without boundary constraints (magenta), with boundary constraints (cyan) and with a block matching algorithm (yellow)

Acknowledgements: INRIA Research Action CardioSense3D, European project euHeart (FP7/2007-2013 no 224495), Microsoft Research.


Barbara André

Tom Vercauteren

Aymeric Perchant

Anna M Buchner

Michael B Wallace

Nicholas Ayache

Introducing space and time in local feature-based endomicroscopic image retrieval

Keywords: B., k=1 B., k=2 B Correct Outlier

Interpreting endomicroscopic images is still a significant challenge, especially since one single still image may not always contain enough information to make a robust diagnosis. To aid the physicians, we investigated some local feature-based retrieval methods that provide, given a query image, similar annotated images from a database of endomicroscopic images combined with high-level diagnosis represented as textual information. Local feature-based methods may be limited by the small field of view (FOV) of endomicroscopy and the fact that they do not take into account the spatial relationship between the local features, and the time relationship between successive images of the video sequences. To extract discriminative information over the entire image field, our proposed method collects local features in a dense manner instead of using a standard salient region detector. After the retrieval process, we introduce a verification step driven by the textual information in the database and in which spatial relationship between the local features is used. A spatial criterion is built from the co-occurence matrix of local features and used to remove outliers by thresholding on this criterion. To overcome the small-FOV problem and take advantage of the video sequence, we propose to combine image retrieval and mosaicing. Mosaicing essentially projects the temporal dimension onto a large field of view image. In this framework, videos, represented by mosaics, and single images can be retrieved with the same tools. With a leave-n-out cross-validation, our results show that taking into account the spatial relationship between local features and the temporal information of endomicroscopic videos by image mosaicing improves the retrieval accuracy.

B.

Introduction

With the recent technology of probe-based confocal laser endomicroscpy (pCLE) [START_REF] Wallace | Probe-based confocal laser endomicroscopy[END_REF], endoscopists are able to image tissues at microscopic level with a miniprobe, and in real time during ongoing procedure. However, as the acquired pCLE images are relatively new for them, the physicians are still in the process of defining a taxonomy of the pathologies in the images, for instance to differentiate benign tissues and neoplastic, i.e. pathological, tissues of colonic polyps, see Fig. 1 and 2 for an illustrative example of such images. To face this clinical challenge, a valuable aid to the physician in establishing a diagnosis would be to provide endomicroscopic images that have a similar appearance to the image of interest and that have been previously diagnosed by expert physicians. Knowing that pathological tissue is characterized by some irregularities in the cellular and vascular architecture, we aim at retrieving texture information coupled with shape information by using local operators on pCLE images. To serve that purpose, we decided to investigate a modern method for content-based image retrieval (CBIR), the bag-of-visual words (BVW) method [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF]. BVW has been successfully used in many applications of computer vision. For example, on a well-defined non medical application, by using this method on a large variety of images of natural or artificial textures, the authors of [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF] obtained excellent recognition results that are close to 98%.

The standard BVW method detects salient regions in the images and extracts information only on these specific regions. However in pCLE images, the discriminative information is distributed over the entire image field. Contrary to classical methods that apply sparse detectors, we use a dense detector to collect densely the local features in the images. This overcomes the information sparseness problem. Moreover, pCLE images contain characteristic pattern at several scales, in particular the microscopic scale of individual cells and the mesoscopic scale of groups of cells. For this reason, we perform a bi-scale description of the collected image regions. Another problem is that the spatial relationship between the local features is lost in the standard BVW representation of an image, whereas the spatial organization of cells is highly discriminative in pCLE images. So we looked at measuring a statistical representation of this spatial geometry. This was achieved by exploiting the co-occurence matrix of the visual words labeling the local features in the image. After the retrieval process, we introduce the measured spatial criterion in a verification step that allows to remove outliers from the retrieved pCLE images, which are given by the most similar to queried images. Taking into account the spatial relationship between local features is the main contribution of our study, it can be used as a generic tool for many applications of CBIR. Besides, we noticed that the FOV of single still pCLE images may not be large enough for the physicians to see a characteristic global pattern and make a robust diagnosis. As this limitation cannot be solved by the standard methods, we decided to take into account the time information of pCLE video sequences by considering them as objects of interest instead of still images. More precisely, we use image mosaicing [START_REF] Vercauteren | Robust mosaicing with correction of motion distortions and tissue deformation for in vivo fibered microscopy[END_REF] to project the temporal dimension of video sequences onto a large FOV image, cf. some resulting mosaics in Fig. 4 and5. With a leave-n-out cross-validation, classification experiments on the pCLE database serve the validation of the methodology: our method outperforms other methods taken as references, by improving the classification accuracy and by providing more relevant training images among the first retrieved images.

The bag-of-visual words method

As one of the most popular method for image retrieval, the BVW [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF] method aims at extracting a local image description that is both efficient to use and invariant with respect to viewpoint changes, e.g., translations, rotations and scaling, and illumination changes, e.g., affine transformation of intensity. Its methodology consists in first finding and describing local features, then in quantizing them into clusters named visual words, and in representing the image by the histogram of these visual words. The BVW retrieval process can thus be decomposed into four steps: detection, description, clustering and similarity measuring, possibly followed by a classification step for image categorization. The detection step extracts salient regions in the image, i.e. regions containing some local discriminative information. In particular, corners and blobs in the image can be detected by the sparse Harris-Hessian (H-H) operator around keypoints with high responses of intensity derivatives. Other sparse detectors like the Intensity-Based Regions (IBR) and the Maximally Stable Extremal Regions (MSER) are also specialized for the extraction of blob features in the images. We refer the interested reader [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF] for a survey of these detectors.

Then, each local region can be typically described by the Scale Invariant Feature Transform (SIFT) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] descriptor. We refer the reader [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF] for a survey of this and other powerful descriptors. At the description step, the SIFT descriptor computes, for each salient region, a description vector which is its gradient histogram at the optimal scale provided by the detector, the gradient orientations being normalized with respect to the principal orientation of the salient region. As a result, the image is represented in a high dimensional space by a set of SIFT description vectors that are invariant by translation, rotation and scale.

To reduce the dimension of the description space, the clustering step, for example based on a standard K-Means, builds K clusters, i.e. K visual words, from the union of the description vector sets gathered from all the N images of the training database. Since each description vector counts for one visual word, an image is represented by a signature of size K which is its histogram of visual words, normalized by the number of its salient regions.

Given these image signatures, it is possible to define a distance between two images as the χ 2 distance [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF] between their signature and to retrieve the closest training images as the most similar to the image of interest. The relevance of the similarity results can be quantified by a further classification step, for instance based on a standard nearest neighbors procedure that weights the votes of the k-nearest neighbors by the inverse of their χ 2 distance to the signature of the queried image, so that the closest images are the most determinant. Besides, performing image classification is a way to validate a new retrieval method by comparing it with other methods. observed that a large number of salient regions sparsely extracted by these standard detectors do not persist between two highly correlated successive images taken from the same video. This is consistent with the rather poor results of the standard sparse detector shown in Section 4. To overcome the persistence problem and take into account all the information in the images, we use a dense detector contrarily to the standard method. Furthermore, a dense detector is relevant because local information appears to be densely distributed over the entire field of the pCLE image. The dense detector is made of overlapping disks localized on a dense regular grid, such that each disk covers a possible image pattern at microscopic level.

We also noticed that the endoscopists establish their diagnosis on pCLE images from the regularity of the cellular architecture in the colonic tissue [START_REF] Buchner | High-resolution confocal endomicroscopy probe system for in vivo diagnosis of colorectal neoplasia[END_REF], where goblet cells and crypts are both round-shaped characteristic patterns, but where a crypt has larger size than its surrounding goblet cells so it must not be recognized as the same object. In order to be sensitive to scale changes, our method looked at describing local disk regions at various scales that are not automatically computed, for example by choosing a microscopic scale for individual cell patterns and a mesoscopic scale for larger groups of cells. This leads us to represent an image by several sets of description vectors that are scale-dependent, resulting in several signatures for the image that are then concatenated into one larger signature. This previous observation also suggests that the spatial organization of the goblet cells must be included in the retrieval process because it is substantial to differentiate benign tissues from neoplastic tissues. The authors of [START_REF] Jegou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] previously proposed adding a geometrical verification to take spatial information into account, however their method is based on the assumption that they want to retrieve images of the exact same scene, which is not the case for our application. In like manner, our idea is to introduce a geometrical verification process after the retrieval process, but based on the assumption that the spatial relationships between the local features are only statistically the same in the images with similar appearance. To introduce spatial information, we took advantage of the dense property to define the adjacency between two visual words as the 8-adjacency between the two disk regions that are labeled by them on the detection grid. Thus, we are able to store in a co-occurence matrix M of size K × K the probability for each pair of visual words of being adjacent to each other. In order to best differentiate the images of the benign class from the images of the pathological class, we looked at the most discriminative linear combination W of some elements m of M . This is achieved by a linear discriminant analysis (LDA) which uses the textual diagnostic information in the database. The LDA weights are given by W = Σ -1 (µ 1 -µ 2 ), where Σ is the covariance matrix of the elements m of M in all training images and µ i is the mean of the elements m of M in the training image belonging to the class i. From these weights W , we computed the spatial criterion α = W m for each retrieved image and compared it with the α value of the image of interest. By thresholding the α value during a verification process, outliers are rejected and the first retrieved training images are more relevant.

Expert physicians pointed out that some characteristic global patterns are too partially visible on single still pCLE images to make a robust diagnosis: two still images may have a very similar appearance but be attached to contradictory diagnoses [START_REF] Becker | High resolution miniprobe-based confocal microscopy in combination with video-mosaicing[END_REF]. Fig. 3 shows some bad retrieval results caused by the small-FOV problem. To address this problem, the time dimension of pCLE videos needs to be exploited, by including in the retrieval process the temporal relationship between successive images from the same video sequence. The study reported by [START_REF] Sivic | Efficient visual search for objects in videos[END_REF] proposes a method for video retrieval using the spatial layout of the image regions, but this method has been designed for object matching, which is not our objective. Since successive frames from pCLE videos are only related by viewpoint changes, our approach uses the image mosaicing of [START_REF] Vercauteren | Robust mosaicing with correction of motion distortions and tissue deformation for in vivo fibered microscopy[END_REF] to project the temporal dimension of a video sequence onto one image with a larger FOV and of higher resolution. Thus, mosaics can be queried and retrieved in the same way as still images.

Experiments and discussion

At the Mayo Clinic in Jacksonville, the Cellvizio ® system, MKT, Paris, was used to image colonic polyps during surveillance colonoscopies in 54 patients. On each acquired video sequence, expert physicians established a pCLE diagnosis [START_REF] Buchner | High-resolution confocal endomicroscopy probe system for in vivo diagnosis of colorectal neoplasia[END_REF] that differentiates pathological sequences from benign ones. The video sequences contain from 5 to over a thousand frames and each frame is an image of diameter 500 pixels corresponding to a FOV of 240 µm. To build our pCLE database, we considered a subset of these sequences by discarding those whose quality was insufficient to perform a reliable diagnosis. In each of the remaining 52 video sequences, we selected groups of successive frames according to the length of the sequence. The resulting database is composed of N = 1036 still pCLE images and N ′ = 66 pCLE mosaics, half of the data coming from benign sequences and half from pathological ones.

In pCLE images, the disk regions containing information at mesoscopic scale have a radius value ρ 1 = 40, while the radius value of those containing information at microscopic scale is ρ 2 = 15. For the dense detector, we then chose δ = 20 pixels of grid spacing in order to get a reasonable overlap between adjacent regions. Among the values from 30 to 1500 found in the literature for the number K of visual words provided by the K-Means clustering, the value K = 100 yielded satisfying classification results. To prevent overfitting, as the size of our pCLE database is still rather small, especially concerning the number of mosaics, the number of LDA weights in the computation of the spatial criterion α had to be restricted. For the elements m of the co-occurence matrix M , we only considered the K diagonal elements of the matrix M build from the visual words of large radius 40, observing that the overlapping regions of radius 40 have a sufficient spatial correlation, better than those of radius 15. The good values of the threshold θ α were chosen by analysing the distribution of α across the benign and pathological images: 2 when retrieving still images from a queried still image and 0.5 when retrieving still images from a queried mosaic.

The classification results of our method are presented in Fig. 6 on the left and compared with the following methods taken as references: the standard sparse scale invariant SIFT method, the statistical approach of Haralick features [START_REF] Haralick | Statistical and structural approaches to texture[END_REF][START_REF] Srivastava | Computer-aided identification of ovarian cancer in confocal microendoscope images[END_REF] and the texture retrieval method of Textons [START_REF] Leung | Representing and recognizing the visual appearance of materials using three-dimensional textons[END_REF]. To ensure a non-biased classification, our validation scheme retrieves k nearest images in the training set with training images not belonging to the video sequence of the image being queried, i.e. a leave-n-out cross-validation where n is the number of frames in the video of the queried image. According to the accuracy, sensitivity and specificity rates yielded by each method on the still images of the pCLE database, our retrieval method including spatial information is the most efficient, with an accuracy rate of 78.2% for k = 22 neighbors, which is 11.5 points better than the standard SIFT method. The gain of accuracy can be decomposed in 10.2 points for the choice of a dense detector and a bi-scale SIFT description, and 1.3 points for the verification process on the spatial criterion. It is also worth mentioning that with the spatial verification, fewer nearest neighbors are necessary to classify the query at a given accuracy. For k = 4 neighbors, some illustrative examples of the image retrieval results are shown in Fig. 1, 2 and 3.

Moreover, when including both spatial and temporal information by querying mosaics, our classification results are much better, see Fig. 6 on the right. Since mosaics contain more information than single images, their content-based neighborhood is more representative of their pathological neighborhood, so they can be better classified by a smaller number k ′ of nearest neighbors. Indeed, if we retrieve still images for queried mosaics, the classification accuracy is 83.3% for k ′ = 10 neighbors, which demonstrates the robustness of our retrieval method applied on heterogeneous data with different resolution. For the retrieval of still images from queried mosaics, the poor specificity can be explained by the fact that a mosaic annotated as neoplastic may contain some benign patterns which induce the retrieval of single benign images and classify it as benign. However the expert physicians diagnose a pCLE video sequence as neoplastic as soon as it contains neoplastic patterns, even when some benign tissue is imaged. Besides, if we retrieve mosaics for queried mosaics, the classification accuracy is 89%. Thus, even though we only have a small number of mosaics, including time dimension in mosaics provides us proof of concept results for endomiscroscopic video retrieval. The relevance of mosaic retrieval can also be qualitatively appreciated in Fig. 4 and5. For the retrieval of mosaics from queried mosaics, including the spatial information does not improve the classification results because of the overfitting phenomenon: indeed, the number of LDA weights, 100, is bigger than the total number of mosaics, N ′ = 66.

Conclusion

Using visual similarity between a given image and medically interpreted images allowed us to provide the physicians with semantic similarity, and thus could potentially support their diagnostic decision. Although our experiments are focused on a relatively small training dataset, the classification results constitute a validation of our generic methodology. By taking into account the spatiotemporal relationship between the local feature descriptors, the first retrieved endomicroscopic images are much more relevant. For future work, a larger training database would not only improve the classification results if all the characteristics of the image classes are better represented, but also enable the exploitation of the whole co-occurence matrix of visual words at several scales. Besides, the learning step of the retrieval process could leverage the textual information of the database and incorporate the spatial information of multi-scale co-occurence matrices into descriptors. The co-occurence matrix could also be better analyzed by other tools more generic than linear discriminant analysis. For example, a more complete spatial geometry between local features could be learned by estimating the parameters of Markov Random Fields [START_REF] Descombes | Estimation of Markov random field prior parameters using Markov chain Monte Carlo maximum likelihood[END_REF]. We also plan to perform a leave-one-patient-out cross-validation to ensure fully unbiased retrieval and classification results. As for introducing the temporal information, a more robust approach would not only consider the fused image of a mosaic but the 2D + t volume of the registered frames composing the mosaic to work on more accurate visual words and better combine spatial and temporal information.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Typical image retrieval results provided by our method from two benign queries. B. indicates Benign and N. Neoplastic. From left to right on each row: the queried image, and its k-NNs on the top layer, and their respective colored visual words on the bottom layer. An outlier is indicated by Outlier if it has been rejected by the spatial verification process, and by Error otherwise. FOV of the images: 240 µm.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Typical mosaic retrieval results provided by our method from one benign query. B. indicates Benign and N. Neoplastic. From left to right on each row: the queried mosaic, and its k-NNs on the top layer, and their respective colored visual words on the bottom layer. FOV of the mosaics: from 260 µm to 1300 µm.

Including spatial and temporal information

When we applied the standard BVW method on pCLE images [START_REF] André | Endomicroscopic image retrieval and classification using invariant visual features[END_REF], we obtained rather poor classification results and the presence of many retrieval outliers. To improve the accuracy of endomicroscopic image retrieval, we decided to include both spatial and temporal information contained in the pCLE images. While testing on pCLE images the numerous sparse detectors listed by [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF], we first


B André

T Vercauteren

A Perchant

A M Buchner

M B Wallace

N Ayache

ENDOMICROSCOPIC IMAGE RETRIEVAL AND CLASSIFICATION USING INVARIANT VISUAL FEATURES

Keywords: Endomicroscopy, content-based image retrieval, Bag of Visual Words (BVW) method, k-nearest neighbors classification

This paper investigates the use of modern content based image retrieval methods to classify endomicroscopic images into two categories: neoplastic (pathological) and benign. We describe first the method that maps an image into a visual feature signature which is a numerical vector invariant with respect to some particular classes of geometric and intensity transformations. Then we explain how these signatures are used to retrieve from a database the k closest images to a new image. The classification is finally achieved through a procedure of votes weighted by a proximity criterion (weighted k-nearest neighbors). Compared with several previously published alternatives whose maximal accuracy rate is almost 67% on the database, our approach yields an accuracy of 80% and offers promising perspectives.

INTRODUCTION AND AIMS

Probe-based confocal laser endomicroscopy (pCLE) is a new technology which enables dynamic microscopic imaging of tissues in vivo with a miniprobe during ongoing endoscopy. Our study investigates the application of image retrieval and classification methods to pCLE images of colonic polyps, which would aid the physician in differentiating benign tissues and neoplastic (pathological) tissues. Due to the very specific nature of endomicroscopic images (see Fig. 1 and Fig. 2), where membranes and nuclei are not always clearly visible, typical criteria used by classical computer aided diagnosis methods (e.g. nucleocytoplasmic ratio in histological slices) cannot be applied. Moreover, the taxonomy of pathologies in endomicroscopic images is still under active construction by the physicians, who are discovering their complex underlying semantic. To face this difficulty, we investigated some methods successfully applied in the field of computer vision, where important progress has been recently achieved using local operators for invariant image description and classification. Indeed, the authors of [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF] reported excellent recognition results on textured images with a Bag of Visual Words [START_REF] Sivic | Efficient visual search for objects in videos[END_REF] (BVW) method: the method reaches classification results close to 98% on a large variety of images of natural or artificial textures at various scales. This motivated us to adapt their approach to our pCLE images which also tend to contain discriminative texture information (coupled with shape information) at various scales. Thus, our objective is to find the images which are the most similar to a given image by exploring a content based image retrieval (CBIR) approach, and to quantify the relevance of the similarity results by performing a supervised binary classification of the database.

Section 3 explains the contributions of our methodology: the use of a dense detector for salient regions to describe the information over the entire image field, the concatenation of signatures of the same image that is described at different scales (each scale corresponding to a physical group of image patterns), and the supervised selection of the most discriminative visual words to improve the classification results. The performance comparisons between several methods applied to the training set of pCLE images is presented in Section 4; in particular we show that, with a leave-n-out cross-validation, our method outperforms the method of [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF] mentioned above.

MATERIALS

For our study, colonic polyps were imaged at the Mayo Clinic in Jacksonville using the Cellvizio R system (MKT, Paris), with a prior administration of fluorescein, during surveillance colonoscopies in 54 patients. On each acquired video sequence, the expert physicians established a pCLE diagnosis [START_REF] Buchner | High-resolution confocal endomicroscopy probe system for in vivo diagnosis of colorectal neoplasia[END_REF] differentiating pathological sequences from benign ones, according to the presence or not of neoplastic tissue which is characterized by some irregularities in the cellular and vascular architectures. The video sequences contain from 5 to over a thousand frames (images with a circular shape of diameter 500 pixels and with a field-of-view of 240 µm). We considered a subset of these sequences by discarding those whose quality was insufficient to perform a reliable diagnosis, or whose pCLE diagnosis was not the same as the "gold standard" (the "gold standard" for the polyp was established by a pathologist, after histological review). In each of the 52 video sequences that were retained, we selected groups of successive frames according to the length of the sequence.

The resulting database is composed of N = 1036 endomicroscopic images, half of the data coming from benign sequences and half from pathological ones.

METHODS

The primary goal of local methods for object recognition is to ensure extraction and description of features invariant w.r.t. viewpoint changes (e.g., translations, rotations and scaling) and illumination changes (e.g., affine transformation of intensity). One of the most popular methods for image retrieval using invariant features is the BVW [START_REF] Zhang | Local features and kernels for classification of texture and object categories: a comprehensive study[END_REF][START_REF] Sivic | Efficient visual search for objects in videos[END_REF] method, that couples for example the sparse Harris-Hessian (H-H) detector of salient regions with the Scale Invariant Feature Transform (SIFT) [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] descriptor of these regions. Given the success of this method, we decided to develop our methodology by adapting the BVW to pCLE images. Our approach is composed of three steps: the detection step, the description step and the classification step.

1) The detection step consists in selecting salient regions in an image, i.e. regions containing some local discriminative information. It is worth noticing that the physicians establish their diagnosis from the regularity of the cellular architecture in the colonic tissue, where goblet cells and crypts are roundshaped or tubular-shaped. For this reason, we looked at extracting blob features in the images by some sparse detectors, like the H-H, the Intensity-Based Regions (IBR) and the Maximally Stable Extremal Regions (MSER) detectors (see for example [START_REF] Mikolajczyk | A comparison of affine region detectors[END_REF] for a survey of these methods). In particular, the H-H operator detects corners and blobs in the image around key-points with high responses of intensity derivatives for at least two distinct gradient directions. In our pCLE images, we observed that a large number of salient regions extracted by these sparse detectors do not persist between two highly correlated successive frames, although some of them were localized in the center of crypts or in the neighborhood of other colonic patterns. The weak robustness of sparse detectors applied to pCLE images explains the "poor" classification results presented in Section 4; it could be detailed in another study. To take into account local information over the entire image field, the idea is to use a dense detector made of overlapping disks of constant radius distributed on a dense regular grid, such that each disk covers a possible image pattern at microscopic level.

2) For the description step, the standard SIFT descriptor computes a 128-bin gradient histogram as a description vector for each salient region centered on the key-points at optimal scales, the gradient orientations being normalized with respect to the main orientation of the salient region. As a result, each image is represented by a set of description vectors in a high dimensional space. To reduce the dimension of the description space, a standard K-Means clustering builds K clusters named "visual words" from the union of the vector sets representing all the N images of the database. Thus, an image is represented by a signature of size K which is its normal- ized histogram of visual words, each description vector counting for one visual word. The standard BVW method is at the same time translation, rotation and scale invariant. However, the level of required invariance depends on the problem at hand: we desire invariance by translation and rotation, but we do not want invariance by scaling because in colonic tissue, round-shaped crypts have larger size than round-shaped goblet cells and thus must not be recognized as the same object. For the feature description to be sensitive to scale changes, we decided to perform multiple SIFT descriptions, each of them being associated to a different value of disk radius: for instance, two different radius values lead to represent an image by two sets of description vectors, hence by two signatures that are then concatenated into one larger signature.

3) The classification step is a standard nearest neighbors procedure: given these image signatures, it is now possible to define a distance between two images as the χ 2 distance between their signature. Besides, we weighted the votes of the k-nearest neighbors by the inverse of their χ 2 distance to the tested image signature, so that the closest images are the most determinant. We also improved the classification performance by selecting in a supervised way the most discriminative visual words, i.e. those minimizing the intra-class distances while maximizing the inter-class distances. For each class C of images, we considered the distribution p(w|C) of the number of occurrences of a visual word w in the images belonging to the class C. The discriminative power f (w) of the visual word w is chosen by using the Fisher criterion which can be expressed as the Mahalanobis distance between the two distributions p(w|Benign) and p(w|Pathological):

f (w) = (µ 1 -µ 2 ) 2 /(0.5 (σ 2 1 + σ 2 2 ))
, where µ i and σ i are respectively the mean and the variance of the distribution of w in the images belonging to class i. Furthermore, by reducing the number of visual words, the size of image signatures is decreased, so the image retrieval and classification processes run faster.

RESULTS

For the dense detector, a large disk radius of ρ 1 = 40 pixels is relevant to cover groups of cells with a disk; a smaller disk of radius ρ 2 = 15 pixels allows to cover at least one cell in the images (see Fig. 1). Given these radii values, we chose δ = 20 pixels of grid spacing in order to get a rea- To compare our approach with other methods, global as well as local, the validation scheme consists in the following non-biased classification: the k nearest images in the training set are retrieved, with training images not belonging to the video sequence of the image being tested (i.e. leave-nout cross-validation, where n is the number of frames in the video of the tested image); then the votes of the closest image signatures are weighted by their distance to the tested signature. For performance comparison, the following methods are taken as references: the sparse scale invariant SIFT method [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF], the statistical approach of Haralick features [START_REF] Haralick | Statistical and structural approaches to texture[END_REF][START_REF] Srivastava | Computer-aided identification of ovarian cancer in confocal microendoscope images[END_REF] and the texture retrieval method of Textons [START_REF] Leung | Representing and recognizing the visual appearance of materials using three-dimensional textons[END_REF]. The Haralick method computes global statistics from the co-occurrence matrix, such as contrast, correlation or variance, so as to represent an image by a vector of statistical features; this method is worth being compared with, because of its global scope. The last reference approach, the Textons method, defines for each image pixel p a "texton", as the response of a patch centered on p to a texture filter which is composed of orientation and spatial-frequency selective linear filters. While only tex- ture information is extracted by this method, the fact that its extraction procedure is dense makes it interesting for method comparison.

After the classification process has been applied, accuracy, sensitivity and specificity rates are computed for each method, as shown in Fig. 4. According to these values, our method is the most efficient on the training database, with an accuracy rate of 80.12%, which is 13.42 points better than the standard SIFT method: the gain of accuracy can be decomposed in 9.85 points for the choice of a dense detector, 1.45 points for the bi-scale SIFT description and 2.12 points for the supervised selection of the "best" visual words. These results point out that, for our training data, the standard SIFT method, the Textons method and the Haralick method are clearly less efficient than our approach.

CONCLUSION

To our knowledge this is the first attempt to classify endomicroscopic images by adapting a recent and powerful local image retrieval method, the Bag of Visual Words method, introducted for recognition problem in computer vision. Although our study has been focused on one relatively small training set of pCLE images of colonic polyps, the classification results show that our objectives have been successfully reached by the methodology, which consists first in densely collecting feature information, secondly in describing this information at two different scales, one corresponding to microscopic features (cells) and the other to mesoscopic features (groups of cells), and in selecting the most discriminative visual words.

Besides, our generic framework could be reasonably extended to many other image retrieval applications.

As for future work, efforts will be made to strengthen the validation of our method, for example by using ROC curves to better estimate the parameters. Concerning the classification procedure, probabilistic models with hidden variables [START_REF] Hörster | Comparing local feature descriptors in pLSA-based image models[END_REF] could be more efficient than the voting algorithm. Moreover, the classification results would be improved if a larger training database was considered, where all the characteristics of the image classes are better represented. Improvements could also be made for the treatment of outliers in similar images by including the temporal dimension of video sequences, which would address at the same time the problem of noise, motion distortions or partially visible macroscopic features (e.g. crypts) in a given frame. A more complete study focused on content based image retrieval could take into account, in the image description, spatial relationships between the salient regions to close the gap between highly local image features and more global image features.

Fig. 1 :

 1 Fig. 1: Disk regions (of radii 15 and 40 pixels) associated to the most discriminative visual words (represented by colors). Left: Benign pCLE image. Right: Neoplastic (pathological) pCLE image. Field of view of the images: 240 µm.

Fig. 2 :

 2 Fig. 2: Four rows of similar pCLE images provided by our method. From left to right on each row: the tested image, and its first, second, third and fourth most similar images. F indicates the frame number of the image and S the number of the video sequence. Field of view of the images: 240 µm.sonable overlap between adjacent regions. The number K of visual words provided by the K-Means clustering was selected among values from 30 to 1500 given by the literature: the value K = 100 yielded satisfying classification results. The K ′ most discriminant visual words are selected by applying on their discriminative power a hand-picked threshold, θ = 0.25, which gives good classification results when testing the whole training set without cross-validation. This threshold θ is applied inside the cross-validation loop to select a certain number of discriminant visual words, and the mean value of K ′ for all cross-validations is 40. Finally, the value chosen for the number k of nearest neighbors is the global value maximizing the classification accuracy on the graph of accuracy rates presented in Fig.3(where the best result is reached for k = 42). Some illustrative examples of the resulting image similarities are shown in Fig.2.

Fig. 3 :

 3 Fig. 3: Accuracy rates with leave-n-out cross-validation.

Fig. 4 :

 4 Fig. 4: Classification results (accuracy, sensitivity, specificity) with leave-n-out cross-validation for k = 42 nearest neighbors.