Image Segmentation using MRFs and Statistical Shape Modeling

Keywords: Shape Modeling, Segmentation, Markov Random Fields, Discrete Optimization xiv List

In this thesis, we introduce a new statistical shape model and use it for knowledge-based image segmentation. The model is represented by a Markov Random Field (MRF). The vertices of the graph correspond to landmarks lying on the shape boundary, whereas the edges of the graph encode the dependencies between the landmarks. The MRF structure is determined from a training set of shapes using manifold learning and unsupervised clustering techniques. The inter-point constraints are enforced using the learned probability distribution function of the normalized chord lengths.

This model is used as a basis for knowledge-based segmentation. We adopt two approaches to incorporate the data support: one is based on landmark correspondences and the other one uses image region information. In the rst case, correspondences between the model and the image are obtained through detectors and the optimal conguration is achieved through combination of detector responses and prior knowledge. The second approach consists of minimizing an energy that discriminates the object from the background while accounting for the shape prior. A Voronoi decomposition is used to express this objective function in a distributed manner using the landmarks of the model. Both algorithms are optimized using state-of-the art ecient optimization methods.

We validate our approach on various 2D and 3D datasets of images, for computer vision applications as well as medical image analysis.

La structure du champ de Markov est déterminée à partir d'un ensemble de formes, en utilisant des techniques d'apprentissage de variétés et de groupement non-supervisé. Les contraintes entre les points sont assurées par l'estimation des fonctions de densité de probabilité des longueurs de cordes normalisées.

Dans une deuxième étape, nous construisons un algorithme de segmentation qui intègre le modèle statistique de forme, et qui le relie à l'image grâce à un terme région, à travers l'utilisation de diagrammes de Voronoi. Dans cette approche, un contour de forme déformable évolue vers l'objet à segmenter.

Nous formulons aussi un algorithme de segmentation basé sur des détecteurs de points d'intérêt, où le terme de régularisation est lié à l'a priori de forme.

Dans ce cas, on cherche à faire correspondre le modèle aux meilleurs points candidats extraits de l'image par le détecteur. L'optimisation pour les deux algorithmes est faite en utilisant des méthodes récentes et ecaces.

Nous validons notre approche à travers plusieurs jeux de données en 2D et en 3D, pour des applications de vision par ordinateur ainsi que l'analyse d'images médicales. I am indebted to my oce mates Charlotte, Martin, Régis, Noura, Olivier J., Daniel, Fabrice, Pacale and Bo for creating a friendly atmosphere and an inspiring environment. I was also very lucky to be surrounded by other outstanding people including Salma, Aris, Olivier T., Loïc, Panos, Katerina, Maxime, Mickaël, Chaohui, Iasonas, Dimitris, etc. We shared very nice moments and I enjoyed traveling with them to attend conferences, after some long nights of work spent at the lab to meet the submission deadlines. We also had the opportunity to go sightseeing and to visit wonderful places together, and the 2010 Easter in Greece will remain among my best memories. Thank you for your team spirit and for helping me out in the dicult times. It was a great pleasure to work with you. Special thanks go to Fabrice for reading my thesis manuscript and for his very helpful comments and suggestions. I would also like to thank

Sylvie and Annie for their help and kindness, for so nicely taking care of us in the lab, and for making things so easier.

Thanks to my wonderful friends Zied, Mathieu, Benoit, Mourad, Aymen, Imen, Atef, Tarek, Asma, Leila, Hamza, Majed, Dorra, Sami and Heyfa, my life in Paris was lled with joy and happiness. I am grateful to them for their support, understanding and company throughout the years. My thanks also go to my dear friends back home in Tunisia, Jihen, Soane and Jihen, who always wished me all the best for my studies. To Souene my longtime great friend, thank you for your kindness and your faithfulness.

I would like to express my extreme gratitude to my dear friend and colleague Radhouène. We made this long journey together, from the very beginning to the very end, and he was always there for me. His rigor as a scientist and his uprightness as a man always inspired me respect and admiration. On top of the numerous discussions where his bright intuition always proved to be true, I gratefully acknowledge him for being my unocial thesis reader. His precise and constructive input and constant moral support during the time pressure of thesis dissertation was invaluable. Coming back to Paris from abroad to attend my defense truly warmed my heart. Thank you Radhouène for being this pillar that never shivered throughout the years.

This thesis would have remained a dream had it not been for my beloved girlfriend Memia. She shared with me this wonderful experience and backed me at every moment, with patience and passion. Thank you Memia for everything.

Last but not least, I would like to pay a tribute to my family: my bothers Slim and Raouf, my father Abdelwaheb and my mother Salsabil. Thanks to their unconditional love, their eternal care and and their endless support, I coped with the most dicult and tough situations. Without the high ideals of Knowledge and Science that my parents instilled in me since my little childhood, without the enormous sacrices they made for my well-being and my education, this work would not have been possible.

To you, Dad and Mom, is dedicated this thesis.

Chapter 1

Context and Motivation

Image segmentation is a fundamental task in computer vision and medical image analysis. It can be dened as the extraction of the boundaries of an object of interest from the image. Alternatively, it can be viewed as a labeling problem, where image pixels are assigned either to the object or to the background. Segmentation is of a great importance because it provides qualitative and quantitative information about the image. For instance, semi-automatic photo-editing is a popular application that benets from segmentation techniques. Behavior analysis is another domain where segmentation represents a data pre-processing step. Note that the accuracy of the extracted boundaries would highly inuence the quality of the resulting model. Image/videobased automatic translation of hand signs is a concrete example of gesture recognition that could benet from ecient segmentation as well. Such an application can hardly be achieved without precisely outlining the hand and detecting the relative nger poses. In medical image analysis, clinical routines have integrated segmentation algorithms, to assist diagnosis. These tools allow physicians to gain precious time, alleviate the need for tedious manual processing, while providing satisfactory accuracy. For instance, assessing the left ventricular function of the heart is done through the estimation of the left ventricle ejection fraction. The latter is a quantitative indicator that can be computed from segmented magnetic resonance cardiac images.

A common methodology to address computer vision tasks denes a mathematical model that represents the problem, and then extracts the model parameters from image observations using optimization methods. Low-level segmentation approaches exclusively rely on neighboring image features and grouping methods. These techniques have limited success when compared to high-level methods that integrate shape priors in their modeling. In elds where a prior knowledge is available (like medical imaging), such methods carry on great potentials since the domain knowledge can be used to intro-duce constraints and improve the reliability and accuracy of the segmentation result. In this challenging task, one rst has to determine a model representing these constraints and then an inference process which aims to combine the visual support with the prior knowledge.

Dierent classes of shape models have been used for segmentation purposes. These models are linear or non-linear, global or local, represent the shape in an explicit or implicit manner, and may benet from statistical learning. Although global and linear statistical shape modeling approaches are very popular in the computer vision community, and have been applied

to a wide variety of segmentation problems, we believe that non-linear local shape priors can be very powerful. They can be further enhanced if they are combined in a structured manner, which can be derived from shape populations using manifold learning techniques. Such local object parts interactions are also very suitable for the use of Markov Random Fields (MRF) optimization techniques.

In this work, we made the choice to use discrete optimization methods. Despite more modeling constraints with respect to the continuous optimization setting, discrete methods have very appealing properties in general. Variational methods are known to converge to local minima in the general case, without giving any quality evaluation of the obtained solution. They require the computation of the gradient of the energy functional, which might not be dierentiable in some cases, and also depend on a careful choice of an approximate numerical scheme, which very often leads to slow algorithms. On the contrary, the recent developments in the discrete optimization of MRFs provide tools that guarantee convergence in certain cases, provide suboptimality bounds beforehand, and per-instance tighter bounds in practice. Although they converge in general to a local solution as well, these methods could yield a better result as they permit to "jump" over local minima. They do not require the dierentiation of the objective function, exhibit numerical robustness and are computationally ecient. If such a method gives an unsatisfactory result, then it is probably due to an erroneous modeling rather than to numerical instabilities.

We introduce in this thesis a new statistical shape model, in order to use it for segmentation tasks. This goal adds to the natural requirement of proper shape variation representation the necessity of a tractable inference procedure.

Our statistical shape model (SSM) relies on three building blocks that encapsulate geometric, manifold learning and graph theoretic aspects. Essentially, we adopt a landmark-based representation of shapes, that is embedded in a MRF, with an appropriate structure, namely a k-fan structure. Probability density distributions of the normalized inter-point distances are associated to the MRF to form the SSM, leading to a representation that is invariant with respect to translation, rotation and scale changes. In our development, we ensure the coherence of the chosen MRF structure with respect to the geometric aspect, namely the normalized distance representation. We also propose a framework that learns a particular instance of this structure for a particular class of objects, given a representative population of shapes. Hence, the nal structure is directly related to the information encoded by the data. Furthermore, we show that the k-fan has desirable graph properties that permit exact inference only in some limited cases, due to the curse of dimensionality.

In practice, we resort to state-of-the art ecient discrete MRF optimization techniques.

The second step of our work attempts to prove the relevance of our model for object segmentation. Our SSM enables the use of general image cues, and is suitable for segmentation algorithms in 2D and in 3D. We propose a method that relies on image region intensities to extract the object boundaries, and we show its robustness with respect to noise and occlusions. We also suggest another framework that benets from multi-scale image features and trained classiers in combination with our SSM. Experiments are performed in 2D and 3D with medical and natural images.

The remainder of this thesis is organized as follows. In chapter 2, we rst review the shape modeling techniques related to segmentation applications.

Then, we develop our statistical shape model by emphasizing our threefold contribution. We discuss the Euclidean distance geometry issues related to our k-fan structure, by making the link with a more general class of graphs called k-trees. We explain then how we use the shape maps technique to embed the shape population in a Euclidean space, where we apply an unsupervised linear programming-based clustering method to derive the particular graph structure that represents the data. The obtained graph being chordal, exact inference is possible using dynamic programming, if the search space is small enough. After demonstrating the adequacy of the shape modeling technique on some practical examples, we proceed to the development of the segmentation algorithm in chapter 3. We rst give an overview of the discrete MRF optimization methods. Then, we describe the segmentation techniques in computer vision, and especially the graph-based methods, discussing the relation of our method to the previous work. We particularly notice that incorporating statistical shape priors in MRFs is a dicult task. Then, we develop our aforementioned segmentation algorithms. Furthermore, extensions of our segmentation approach to other application settings are presented in appendix A.

More precisely, we propose a method to incorporate dynamic shape priors in the MRF formulation in order to perform object tracking. In a dierent application, we use graph matching techniques to extract discriminatively-matched pairs of landmarks towards image registration. We conclude this work by

Introduction

Bottom-up computer vision methods rely on comparing and grouping pixel features in local image neighborhoods with little prior knowledge. Top-down approaches is an alternative that gained a signicant popularity in computer vision. These techniques allow the combination low-level cues from the image in a principled manner, based on a priori constraints. The learning of such prior information is central to a number of computer vision applications, such as recognition, segmentation and tracking of objects. In particular, shape models proved to be very useful in the context where one can determine a priori knowledge of the object of interest. Such models incorporate information about the geometry of the object, and can also describe its visual image features. This method applies to a wide range of applications in general computer vision (for instance faces, eyes, hands, human bodies...) and medical image processing as well (modeling of the heart, liver, corpus callosum...).

Building a shape model involves in general two components:

• A mathematical shape representation is dened.

• Constraints are added to this representation such that an admissible set of shapes is constructed. These constrains can be for instance learned as a statistical model that represents the variations of the shape.

The literature of shape modeling is very rich, and dierent representations were used in the context of computer vision. In particular, continuous eorts have been made to incorporate shape priors in segmentation algorithms, since the early 1990's. In [START_REF] Heimann | Statistical Shape Models for 3D Medical Image Segmentation: A Review[END_REF], a review of statistical shape models (SSM) for three dimensional medical image segmentation is provided.

Another review of statistical shape priors for level set segmentation is available in [Cremers et al. 2007]. We now review the major shape modeling techniques previously proposed in the literature.

A Review of Shape Models

Point-Distribution Models -Linear Models

One of the most widely used algorithms in the area of shape modeling is the active shape model (ASM), also called point distribution model (PDM).

Point-based or landmark-based representations have been studied, developed and successfully applied in vision applications as well as medical image analysis following the pioneering work of [Cootes et al. 1995]. In this approach, the shape is represented as a vector of concatenated point coordinates in a space of dimension n×d, where n is the number of points, and d is the dimension of the points space. These points are in general boundary points, and the shape is easily recovered by connecting them, to form a meshed surface or a 2D curve.

Then, a (linear) statistical model is learned from a training set of t shape instances, where correspondences between the landmarks are known. Very often, the correspondences are obtained through the manual labeling of the training shapes. However, such data are not always available (especially for 3D medical applications). Moreover, the manual precision of the correspondences is questionable in general, and can vary from one person to another.

For these reasons, dierent methods for automatic correspondence nding were developed. For instance, [Davies et al. 2002] propose a population-based formulation to the problem by minimizing an objective function that features the minimum description length (MDL) of the statistical model.

The learning phase of ASM is based on principal component analysis (PCA). First, an alignment of the training shapes is performed, to lter out translation, rotation and scale change eects. This can be done using various rigid alignment methods like Procrustes analysis. Next, a mean shape s and a covariance matrix C s are computed. An eigendecomposition of C s gives the modes of variations of the shape with respect to the mean shape. These modes are represented by the eigenvectors e i of the covariance matrix. By projecting the residual of an instance shape s from the mean on the obtained basis of eigenvectors, it can be expressed as a linear combination of the eigenmodes.

More explicitly, we have s = s + n×d i=1 b i e i . In the last step, the dimensionality of the problem is reduced and only the eigenvectors corresponding to the largest m eigenvalues are kept, where in general m is chosen such that a large part of the trace of C s is retained. This yields the compressed expression of a shape s = s + m i=1 b i e i with a small residual error. The admissible set of shapes is generally dened by constraining the coecients b i in the range [-3 √ λ i , 3 √ λ i ], where λ i is the eigenvalue corresponding to the eigenvector e i .

Point distribution models are hence linear models that make the Gaussianity assumption of the data, and rotate the coordinate axes to the directions of maximum variance. They are global shape models of local variations (with respect to the mean).

Because of their global eects, eigenmode variations of the PDM result in a change of all the landmarks. Each single mode controls the variations of the whole shape. In some applications, shape model sparsity is a desired property. For instance, PDM can handle the variations of a healthy organ, since the global shape assumption is reasonable in this case. However, for pathologies that aect the organ (and hence its shape) locally, the success of a PDM might be limited. This observation motivated the work of [START_REF] Stegmann | [END_REF]] who looked for variation modes that aect the shape locally. This is achieved through the Orthomax method that rotates the eigenmodes such that a sparsity criterion is optimized. This leads to a similar linear orthogonal decomposition of the shape as in PDM, yet with local inuence properties. Independent component analysis (ICA) [START_REF] Hyvärinen | [END_REF] is a dierent approach that tackles sparse modeling. It supposes that the measured signal is a linear mix of non-Gaussian source signals that are statistically mutually independent (without orthogonality constraints). It was for example used by [Üzümcü et al. 2003] to segment the left ventricle (LV) from magnetic resonance (MR) images in an active appearance model (AAM) [Cootes et al. 2001] framework. A variation on the PCA eigenmodes rotation is also introduced by [START_REF] Alcantara | [END_REF]. Following the ideas presented in [START_REF] Chennubhotla | Chakra Chennubhotla and Allan Jepson. Sparse PCA Extracting Multi-scale Structure from Data[END_REF], they dene an energy function that is a trade-o between a variance term and a locality term. The optimization is performed by a PCA-based initialization, followed by the rotation of pairs of eigenmodes in their plane, such that orthonormality is preserved. These methods however suer from the lack of natural order in the variation modes, in contrast with PCA, where the eigenmodes order is specied by variance.

Hence, the authors used dierent methods and criteria in order to cope with this limitation.

Another problem arises in practice when dealing with PDM. The number of available shapes in the training set is often insucient, especially in the 3D case, where a manual labeling of medical data can be very tedious while establishing correspondences can also be very challenging. The maximum number of degrees of freedom that can be described by the training set is t -1, and so the number of possible eigenmodes. In this case, the learned PDM model would result in incorrect tting during the search step. [START_REF] Cootes | [END_REF] attempt to solve this problem by directly acting on the covariance matrix, and adding (synthetic) variance and covariance values for neighboring curve points, to correlate their displacements.

[ Wang & Staib 1998] use the same idea to build a exible deformable model.

They keep the mean shape, learned from the training data, and discard the covariance matrix, replacing it by the so-called smoothness covariance matrix.

The latter contains only variance values, and covariance term for neighboring landmarks on the curve. This model does not depend on statistical shape variations, and allows the deformation of the shape under smooth boundary constraints.

Linear shape models have also been used in the variational framework. A common criticism that is formulated to PDM regards their limited ability to model shapes with relative rotations, or bending variations, which makes the assumption of linearity (or locality on the Riemannian manifold of shapes) not valid any more. This limitation was indeed described by [Cootes et al. 1995], and motivated the development of non-linear models.

Non-Linear Models

Dierent approaches were considered to model non-linearity of shape variations. Among others, [Heap & Hogg 1997] build a hierarchical shape model by clustering the data into (overlapping) clusters, the number of clusters being specied by the user. They dene afterwards a PDM for each cluster, with bounds on the shape parameters, such that the dened shape region is valid. [START_REF] Cootes | [END_REF]] also use a supervised method to deal with non-linearities. After projecting the training data on the PCA eigenmodes, they estimate a probability density function (PDF) of the shape as a mixture model of Gaussians, where the number of Gaussians is specied beforehand.

The estimation is done using the expectation maximization (EM) algorithm.

In an alternative approach, a non-linear mapping of the data to a feature space is introduced, and PCA is performed in the feature space, supposed to be linear. This is the principle of kernel principal component analysis (KPCA) [START_REF] Schölkopf | [END_REF]] that was applied by [START_REF] Romdhani | [END_REF] to the shape modeling problem in the case of multiple 2D views of a 3D object,

and by [START_REF] Twining | [END_REF] for the modeling of large-amplitude bending objects (nematode worms). The central idea of this approach lies in Mercer's theorem, where the non-linear mapping does not have to be explicitly dened, but can be specied in terms of inner products, expressed using a Mercer kernel. The feature space is called a Reproducing Kernel Hilbert Space (RKHS).

A related method is [Cremers et al. 2002a] where the authors extend their previous Gaussian prior Mumford-Shah formulation [START_REF] Cremers | [END_REF] to the non-linear case. They estimate the distribution of the mapped training data by a Gaussian PDF in the RKHS, and discuss the close relation of their method to the Parzen windows estimator and KPCA. Gaussian kernels were often chosen in the above methods.

[ [START_REF] Charpiat | [END_REF]] develop a method to compute statistics of a training set of contours and to use them as a non-linear shape prior in a variational framework. Considering a dierentiable approximation of the Hausdor distance between shapes, they dene the mean shape as the one that minimizes the sum of distances to all the training samples up to rigid motions, and propose a gradient descent framework to compute it. Then, statistics on the deformation elds that drive this mean shape to each sample are learned using PCA. They are plugged into a variational framework for image segmentation, and dene priors that are invariant to rigid motions.

Implicit Representations

Implicit representation is another typical approach to model shapes, where the contour of the shape is embedded in a high dimensional space and considered as the zero level of a representing surface. The signed distance transform or signed distance map (SDM) of the contour is commonly considered as the embedding function in this case [Osher & Sethian 1988]. [START_REF] Leventon | [END_REF] augment this representation with a statistical model to capture shape variations.

After a rigid alignment of the training contours, they also use PCA, and apply it on the distance transform image collection obtained from the training contours. After dimensionality reduction, eigen distance maps represent the variations of the shapes. This same representation is also used by [START_REF] Tsai | [END_REF], in an ecient level-set segmentation algorithm.

[ Chan & Zhu 2005] dene a distance between two SDM using the Heaviside distribution.

Similarly to the PDM case, criticism was addressed to the above approaches regarding the linearity assumption. Hence, the underlying approximation of the manifold of signed distance maps as a linear vector space is questionable. Linear combinations of signed distance maps do not necessarily give a signed distance map, which can compromise the use of PCA in this case. Moreover, as pointed out in [Cremers et al. 2007], the benets of PCA dimensionality reduction can have a downside. The principal eigen distance maps encode most of the variations of the training set on the space of embedding functions, which does not necessarily correspond to most of the variations of the embedded contours in the learning population. This may result in selecting a large number of eigenmodes to build a satisfying shape model. These considerations motivated the use of the following alternative statistical models.

Implicit representations were also used in [START_REF] Rousson | [END_REF]] and [START_REF] Rousson | [END_REF]]. The prior information was however learned pixelwise (grid-wise) as a Gaussian probability density function. [START_REF] Taron | Uncertainty-Driven Non-parametric Knowledge-Based Segmentation: The Corpus Callosum Case[END_REF] adopt a similar approach and represent shapes as signed distance functions.

However, they deal with the correspondence problem by registering the training shape implicit representations using a topology-preserving free-from deformation algorithm. Then, the registration uncertainty is estimated by comput-ing a covariance matrix at the surface zero-level. The following estimation of the probability density function (PDF) of the shape, using variable bandwidth kernels, takes into account the registration uncertainties. [Pohl et al. 2006] and [START_REF] Pohl | [END_REF] dealt dierently with the problem as they tackled the representation part. They introduced the logarithm of the odds ratio (Lo-gOdds) representation to the area of computer vision. The LogOdds function maps a probability p ∈ P (P being the probability space) to a real value through the logarithm of the ratio between p and its complement 1 -p (the inverse map is the logistic function). The LogOdds maps form a vector space (L n = R n ), and the authors deduce a vector space structure on P n . This allows probabilistic interpretation for addition and scalar multiplication in the LogOdds space L n . A shape silhouette (binary map) can be hence transformed to a LogOdds map, where the zero-crossing corresponds to the object boundary. This framework provides an embedding of the signed distance maps manifold into the LogOdds vector space. Hence, SDMs can be considered as LogOdds maps (while the reverse is in general not true) and interpreted in terms of space conditioned probabilities. Determining the statistics of SDMs, as LogOdds maps, using PCA is hence possible, and the authors incorporate it in a Bayesian classication algorithm.

[ [START_REF] Etyngier | [END_REF]] addressed the problem dierently by using manifold learning techniques. They rst compute the distance between two training shapes as the Sobolev W 1,2 -norm of the dierence between their SDM. These pairwise distances are plugged in the diusion maps algorithm [Coifman & Lafon 2006] to build an approximation to the shape manifold through an isometric mapping of the data. A Delaunay triangulation is then computed in the reduced manifold space, and used to nd the nearest neighbors of a new shape instance, in order to estimate its projection on the manifold. This projection is introduced in a shape prior term in the variational framework, to attract the segmentation result to the dened manifold. [ [START_REF] Metaxas | [END_REF]] combined superquadrics with splines to couple global and local shape parameters, and applied them to the tracking of human limbs. The global parameters enable to t approximately the target shape in the image, and the local degrees of freedom allow the segmentation to be rened.

Parametric Geometric Representations

[ [START_REF] Baba | Multi-Resolution Stochastic 3D Shape Models for Image Segmentation[END_REF] represent the deformable superquadric model in an orthonormal wavelet basis to inherit multi-resolution capability, and apply it to the reconstruction of human brain structure from MRI images. Then [Bardinet et al. 1995] and [START_REF] Bardinet | [END_REF]] also used a deformable superquadric model to segment the left ventricle (LV). The model was combined with free-from deformations (FFD), which allows one to rene the initial guess of the object and to improve the reconstruction.

Medial Models

The medial axis model is a shape representation introduced by [Blum 1973].

It is dened as the set of centers of all internal bi-tangent spheres to the object boundary that are not inscribed in any other such sphere. The medial transform of an object describes it by its medial axis as well as the corresponding radii of the bi-tangent spheres. This approach is used by [START_REF] Pizer | [END_REF], who represent the interior of an object by sheets of medial atoms from which the boundary can be synthesized. This model is called mrep, and besides this discrete 2D version, a three-dimensional extension is proposed in [Pizer et al. 2003]. A continuous formulation of m-reps is presented in [START_REF] Yushkevich | [END_REF]]. M-reps are very suitable for representing shape variations such as thickness change, bending, and widening. However, their non-Euclidean parameters prohibit the use of linear statistical approaches.

[ Fletcher et al. 2004] circumvent this problem by developing a generalization of PCA to the Riemannian manifold setting. In particular, they dene the intrinsic mean, compute the principal components by maximizing a variance criterion, dene the geodesic sub-manifolds (counterparts of linear subspaces)

and devise a way to project on them. They apply the so-called principal geodesic analysis (PGA) to a population of hippocampi in a schizophrenia study.

Frequency-Domain Representations

An alternative approach that considers the Fourier transformation of a contour and uses a probability distribution on the Fourier coecient was adopted in [START_REF] Staib | [END_REF] (and extended by [Chakraborty et al. 1996]). A Gaussian probability distribution was learned for each parameter. This representation was also used by [Székely et al. 1996] who however modeled the variations of the shape using principal component analysis (PCA) in the Fourier coecient space. This work is extended in [Kelemen et al. 1999], where the parametric representation of shapes uses series of spherical harmonics (SPHARM)

and PCA is performed in the shape parameters space as well. The applications dealt with the segmentation in 2D and 3D of brain and cardiac structures. [START_REF] Davatzikos | [END_REF]] represent shape using a wavelet transform of the landmark coordinates. This frequency analysis permits to obtain coarselevel shape information, as well as ne-level details, thanks to the wavelets coecients in each band. A statistical hierarchical model is then built, by applying PCA to the coecients in each band separately. An ASM-based search strategy was applied to the segmentation of the corpus callosum in MR images, as well as the segmentation of hand images. [Nain et al. 2007] and [Yu et al. 2007] extended this approach to 3D via spherical wavelets. Their method is however limited to 3D shapes that have a spherical topology (genus-0 manifolds). [Essa et al. 2009] overcome this limitation by utilizing diusion wavelets [Coifman & Maggioni 2006] as a hierarchical shape representation.

They apply the orthomax algorithm to the obtained coecients which leads to a sparse shape parametrization.

Shape Descriptors-based Models

Other approaches associate to landmarks individual shape descriptors that can have geometric-invariant properties. For instance, [START_REF] Coughlan | [END_REF] dene a deformable template of the hand that is invariant to translation and rotation. The landmarks are associated with the normal vector to the contour at their location. When matching the template to an image, two geometric Gaussian constraints (with mean values from the template) are imposed to the model, and combined with the extracted image features (corners, histograms of gradient orientations,etc...). The rst one concerns the angle between the normals associated to two successive points on the contour. The second one is related to the angle between two successive points on the contour and the local normal vector. The problem is dened as a maximum a posteriori (MAP) optimization in a Bayesian framework. The structure of the prior is a Markov chain which allows the use of dynamic programming to obtain the global solution without a need for initialization. A similar approach is considered in [Coughlan & Ferreira 2002], where the structure of the prior is not a tree, but a complete graph. The minimization of the Markov Random Field (MRF) dened energy is done using loopy belief propagation.

Attribute vectors were also used to describe shape geometry by [START_REF] Shen | [END_REF]] (as a 2D version followed by a 3D version in [Shen et al. 2001]). For each landmark, they used normalized areas of local triangles (respectively normalized volumes of local tetrahedrons in 3D)

formed by connecting each control point to its local neighboring points. The snake internal energy reects the dierence between the model attributes and those of the evolving contour.

[ [START_REF] Belongie | [END_REF] design a dierent descriptor for each landmark, called "shape context", by attaching to each point a log-polar coordinates histogram of all the relative positions with respect to the other shape landmarks. A distance is dened based on this descriptor to estimate the similarity between shapes. It is possible to obtain a similarity-invariant version of the descriptor using a scale normalization, and rotation-invariance is ensured through a local Frenet frame instead of the absolute frame for coordinates calculations.

The work of [Manay et al. 2006] can be viewed as a continuous generalization of "shape context". The authors dene integral invariants as shape descriptors and provide a distance for shape matching. [Hong et al. 2006] use this descriptor for object segmentation.

From Single Templates to Statistical Pictorial Structures

Incorporating priors in the form of a (deformable) shape template is an idea that retained the attention of many researchers. It goes back to the work of [Fischler & Elschlager 1973] who used a mass-spring template as a represen-tation of faces (by modeling the inter-parts dependencies). [Yuille et al. 1989] matched "hand-crafted" models of face parts (eye, mouth) to the image to extract objects. Two-dimensional templates were used as a shape representation in [Tagare 1997]. These templates were closed contours, similar to the object to be segmented. They were associated with orthogonal pre-computed curves This framework allows for learning the model parameters statistically in [Felzenszwalb & Huttenlocher 2005] using a maximum likelihood method. Hence, the algorithm chooses the parameters that best represent the appearance of each part of the model in the training set, as well as the optimal tree structure that best explains the data.

Our Proposed Method

In this work, we endeavor to propose a shape model that overcomes the shortcomings of the existing approaches. First, global models often suer from their lack of exibility. For instance, while aiming to prevent invalid shapes, PDM will fail to capture new shape variations if the training set is too small, and the object projection on the learned subspace can be too constraining.

Adding local degrees of freedom is subsequently a desired property. Moreover, we would like to be able to model non-linear variations, such that the representation is usable for a large class of shapes. A generic formulation that widens the range of applications to an arbitrary object dimension (2D or 3D) is advantageous, while the ability to properly represent the data is of equal importance. Another relevant aspect to be addressed is the tractability of the inference related to the shape model. Eventually, it is widely accepted that shape is the remaining information of the data after ltering location, orientation and scale eect.

Subsequently, we propose a statistical shape model (SSM) that is similarity-invariant. We adopt a landmark-based representation, where the global shape information is expressed as a combination of local (sparse) interactions. Our representation is associated with a particular model structure that properly captures the geometry of the object and exhibits useful mathematical properties. The model parameters are deduced from the data, through manifold learning techniques. Our derivation permits to dene a shape prior as an energy term that is decomposable into pairwise terms. When exact inference is not possible (curse of dimensionality), the latter property opens up new horizons for ecient optimization using state-of-the art MRF methods.

More precisely, the chord lengths of the object, dened by pairs of landmarks, and normalized by the object scale, are used as a description of the shape. These normalized inter-point distances bequeath similarity-invariance to our representation. A statistical model of the shape is then dened by estimating a probability density function for each normalized chord length (NCL) from a training set of images. Akin to the PDM setting, we require correspondences between the control points. However, there is no need to register the shape since the invariance property is intrinsic to our representation. Hence, the model encodes the -possibly nonlinear -variations of the shape, regardless of its location, orientation and scale.

In the following step, we observe that considering the whole set of chords is redundant. A smaller set of such edges, with the distance information, is sucient to capture the shape geometry. However, a particular structure of the retained edges is needed to ensure the geometric modeling coherence. It is hence interesting to see that the k-tree structure (and also the k-fan structure -denitions are given later, see Fig. 2.3 for examples) is equivalent to the fully connected structure in the case of perfect distance matching. Nonetheless, for a given set of control points, dierent instances of such a structure can be dened. On top of these numerous possibilities there is also the variability that a change in the parameter k generates. Since dierent k-trees (whether they have the same k or not) would lead to dierent models, with various degrees of data delity, we aimed to learn the k-fan structure as well as the number of fans k from the data, such that it is best represented by the model.

To this end, we use a manifold learning technique, called shape maps [START_REF] Langs | Modeling the Structure of Multivariate Manifolds: Shape Maps[END_REF], that combines the strengths of diusion maps [Coifman & Lafon 2006] and model description length. This technique embeds the training shape population into a metric space, where the Euclidean distances between the embedded landmarks account for their interdependency in the data. Then, we use an unsupervised clustering technique (based on linear programming (LP)) [Komodakis et al. 2008a] to group the similarly-behaving control points in an automatically-dened number k of clusters. Then, the k-tree construction stems from the obtained clusters.

From a graph-theoretic point of view, k-tree graphs belong to the larger class of chordal graphs, that inherit the property of the perfect elimination scheme. This paves the way to the use of dynamic programming approaches for exact inference, as well as for sampling from the posterior probability distribution. Moreover, since our formulation is based on normalized interpoints distances, it exhibits pairwise energy terms, that are adequate for the use of ecient MRF optimization techniques.

Related Work

Our approach relates to previous art through the use of inter-points-based representations, k-trees and graph structure learning.

Chord-based Representations

The use of chord-length representations of shapes has been to some extent explored in the literature. Distributions of random chord lengths were used for shape matching and recognition in the work of [START_REF] Smith | [END_REF]] and [You & Jain 1984]. They compute for each shape all inter-boundary points distances and normalize them by the length of the longest chord. They estimate then a probability distribution using all the computed normalized chord lengths (NCL). Then they use a Kolmogorov-Smirnov statistical test to compare two shapes. Despite their interesting results, [Mallows & Clark 1970] had showed earlier that this shape descriptor is not information preserving, as the global distributions of chord lengths can be equivalent for two dierent (non-identical) shapes. These early methods conrm the richness of this mathematical representation. Nevertheless, note that our statistical model A representation based on the distances between pairs of landmarks was also considered in the context of PDM [Cootes et al. 1992] [ Taylor & Cooper 1990][Cooper et al. 1991]. A shape is hence represented as a vector of all chord lengths, and PCA is applied to obtain a shape model as a mean shape and eigenmodes. The authors explain in [Cootes et al. 1995] the ability of this rigid-motion-invariant representation to capture nonlinear shape variations (as bending shapes), opposite to conventional PDM. Before further commenting this approach, we need to introduce the following concepts.

Denition 2.1. Euclidean Distance Matrices (EDM).

A matrix D = [d 2 ij ] ∈ R N ×N is an EDM if there exist N points (x 1 , . . . , x N ) ∈ R n , for some n ≥ 1, such that ∀ 1 ≤ i, j ≤ N, d ij = x i -x j 2 . The set of points (x 1 , . . . , x N ) is called a realization of D in R n .
EDM are invariant with respect to translation, rotation and reexion [START_REF] Dattorro | Convex Optimization & Euclidean Distance Geometry[END_REF]]. Hence, two sets of points (x 1 , . . . , x N ) and (p 1 , . . . , p N ) have the same EDM if and only if they are related by an isometry. This observation justies the use of EDM (as being a "strong constraint") to represent shapes. According to [START_REF] Dattorro | Convex Optimization & Euclidean Distance Geometry[END_REF]], the set of EDM is not a linear space, but is a convex cone. Moreover, note the following classical result of Euclidean geometry. Theorem 2.1 ([Schoenberg 1935][Young & Householder 1938]). A matrix

D = [d 2 ij ] ∈ R N ×N admits a realization (x 1 , . . . , x N ) in R n but not in R n-1 if and only if the (N -1) × (N -1) symmetric matrix A = [a ij ] 1≤i,j≤N -1 dened by: a ij = 1 2 (d 2 1(i+1) + d 2 1(j+1) -d 2 (i+1)(j+1) )
is positive semidenite and has rank n.

The proof of this theorem can be found in [So 2007], and is reported in Appendix B. We can see that the set of EDMs that admit a realization in R n is a sub-manifold of a convex cone. For this reason, applying PCA on the set of EDMs may be questionable. With this linearity assumption, the computed chords may be non-Euclidean in n-D, as pointed out in [Cootes et al. 1992].

Furthermore, the produced matrices may violate the EDM assumption. The authors try to cope with this problem by approximating the produced chord lengths using a steepest gradient descent optimization. This local search may generate plausible shapes for simple cases, but is doomed to fail for more complex objects.

Our approach is also related to another line of work by [Caetano et al. 2006] that deals with graph matching.

Based on the

Euclidean geometry theory, they dene a matching criterion between two graphs with respect to the induced distortion, using the inter-graph nodes distances. While their approach does not deal with statistical shape modeling, they show an interesting result that is inherent to the graph structure. Before proceeding, let us give a few denitions. We denote a complete graph with n vertices by K n . A k-tree is dened recursively as follows:

Denition 2.2. k-trees, base k-cliques.

K k is a k-tree. A k-tree T n+1 on n + 1 vertices can be generated from a k-tree

T n on n vertices, by adding a new vertex to T n and connecting it with k edges to some k-clique (a clique with k vertices) of T n . This latter clique is called then a base k-clique. [ Caetano et al. 2006] prove the global rigidity of k-trees in R k-1 , and subsequently show that the k-tree model is equivalent to the fully connected (complete graph) model for the case of exact graph matching in R k-1 (1) .

Note that the use of a shape model based on this structure to segment an image can be thought of as matching a "statistical graph" to some graph ex-tracted from the image. Hence, the latter property is of great interest to our approach. [Caetano et al. 2006] explain however that in case of position jitter (inexact matching case), the k-tree model is only an approximation to the complete model, while being extremely robust. They also discuss the problem of the k-tree selection as dierent k-trees can result in dierent matching accuracies in the inexact matching case. In our method, we do not restrict ourselves to k-trees in R k-1 , but consider the n-trees structure (n ≥ k). We learn n as well as the structure of the tree from the data in a (statistical) manifold learning framework. The work of [Caetano et al. 2006 the boundary of the shape). This representation has the following nice property: the dual graph of the triangulated polygon (which is a 2-tree) is a tree which bequeaths to the shape vertices the perfect elimination order property (see Fig. 2.3(c)). This property permits the denition of segmentation of objects in images as a minimization problem (using edge information as data term and constraining the deformation of each triangle from the template to the image to be as close to a similarity transform as possibleas smoothness / prior term) that can be solved exactly and globally using dynamic programming, by quantizing the possible locations of the vertices (on a grid for instance). However, complexity is cubic to the number of grid points and the model can not be extended in a straightforward manner to 3D (for instance for medical imaging applications). In a second step, the shape model is extended to a statistical representation, by learning for each triangle a mean and a covariance using a local Procrustes analysis (for each triangle). This work is similar in some aspects to the pictorial structure formulation [Felzenszwalb & Huttenlocher 2005] as both structures (in both works) are chordal graphs, which permits to solve the MAP problem in polynomial time where the complexity is related to the size of the maximum clique in the graph. Thanks to this structure, [Felzenszwalb 2005] dene a Bayesian formulation and estimate the parameters and the tree structure from the data in a maximum likelihood (ML) fashion. They are also able to sample shapes from the posterior probability distribution thanks to this same property. The suggested model is similarity-invariant, deals naturally with occlusions and noise and can handle multiple segmentations. In the same line of development, [Crandall et al. 2005] extend the pictorial structure [Felzenszwalb & Huttenlocher 2005] by dening another class of chordal graphs, namely k-fans. These graphs are exactly k-trees with a single base (also called reference part), and we use in our approach the same graph structures (see Fig. 2.3(d), Fig. 2.3(f )). They also extend their Bayesian framework to this case, and learn the statistical parameters (Gaussian distributions on the relative landmarks positions) as well as the reference part (and hence the k-fan structure) from the data. They apply their approach in 2D object recognition. However, in practice, they did not consider similarity invariance (while it can be dened for k ≥ 2 with respect to the reference part, as they precisely explain), and restricted application to k ∈ {0, 1, 2}. This choice is probably guided by the wish of controlling the polynomial complexity of the MAP problem. In our case, the geometric aspects dictate to use k ≥ n + 1 in R n (global graph rigidity issue as explained in the previous section), k is learned from the training set of shapes and ecient approximate optimization techniques [Komodakis et al. 2008b] [Kolmogorov 2006] are used to cope with the curse of dimensionality.

A related segmentation approach in medical image analysis is presented by [Seghers et al. 2007a]. They dene a 2D shape prior using landmarks that form a closed curve. This chain structure enables the use of dynamic programming. Three dierent prior statistics were tested. The rst shape prior estimates a Gaussian distribution of the relative positions of successive landmarks on the curve, and is hence translation-invariant. The second model estimates a Gaussian distribution of the relative positions of successive landmarks on the curve in a local frame. The third model estimates a conditional Gaussian distribution of the relative positions of successive landmarks on the curve with respect to the previous pair of landmarks. The last two models are also rotation-invariant. While the search using the rst prior is quadratic in the number of grid points, it is cubic for the other two priors. This framework is extended to 3D in [Seghers et al. 2007b] using the rst proposed prior with a mesh of the shape. A heuristic search method is used for the optimization, and is called iterative dynamic programming. More details of this segmentation approach will be discussed in the next chapters. Here again, the authors make an explicit assumption about the graph structure, and this issue is not addressed.

Other Related Approaches

We now comment a few additional related papers. [Gu et al. 2007] assume a Gaussian Markov Random Field (GMRF) prior over the landmarks. This prior is fully represented by the mean landmark positions as well as the precision matrix (the inverse of the covariance matrix) over the landmarks. The zero entries in the precision matrix fully determine the graph structure, and account for the conditional independence between the corresponding landmarks.

Unfortunately, the precision matrix is often ill-dened in practice (where the number of training examples can be smaller than the problem dimension).

They circumvent this structure learning problem by using a Lasso regression technique [Meinshausen & Bühlmann 2006]. This leads to an interesting sparse representation. However, the Gaussian assumption is a shortcoming of this model in case of non-linear shape variations.

[ A 2-tree that is not a 2-fan -the reference cliques are emphasized. (f ) A 4-fan.

manner using an incomplete graph, where the cluster centers form a complete graph, and where each cluster forms a complete graph as well. Our approach diers from their in two aspects. We use shape maps [START_REF] Langs | Modeling the Structure of Multivariate Manifolds: Shape Maps[END_REF] to model the joint behavior of two control points, and this measure takes into account the interactions between this pair of landmarks and the remaining ones in the training set. Hence the shape maps technique is able to bring more information about the co-dependencies than the Bhattacharyya divergence. Then, a major dierence lies in the rationale behind the nal structure of the graph. We obtain a k-fan graph, and this choice takes into consideration geometric constraints that will be explained in section 2.2.2. The graph proposed by [Glocker et al. 2009] is not guaranteed to have the same global rigidity property, which is important for our shape modeling. On the contrary, the FFD registration requirements are dierent from these shape considerations. Their model permits to deal with the ill-posedness of the registration problem while allowing discontinuities for the deformations of neighboring control points belonging to dierent clusters.

A Normalized Chord-Length Statistical

Shape Model

We describe here our statistical shape model (SSM) and detail the inherent geometric, manifold learning and graph-theoretic aspects.

Second-Order Approximation and Chord-Length

Knowledge-based segmentation methods are based on the denition of a model which is then combined with image support towards object extraction. Traditional approaches consist of representing the shape using a number of landmarks and learning their behavior using a training set.

In this work, we consider that a shape is dened by its boundary (a curve in 2D, or a surface in 3D) as well as a set of control points. Knowing the positions of these controls points, the shape boundary can be (approximately) retrieved by the use of an interpolation method. The control points can in general belong to the boundary, or not. In practice however, it is natural that such landmarks are chosen to lie on the contour of the object, which exhibits very often distinctive features. Therefore, we will assume in the following such a hypothesis, meaning that the control points are boundary points.

We represent a shape as a Markov Random Field (MRF). Let G = (V, E)

be a graph where the nodes V = {x 1 , . . . , x n } represent the control points positions, and where the set of edges E represents the conditional dependencies between the landmarks. We assume for the time-being that this graph is complete, that is we take into account all the pairwise dependencies between the landmarks. Hence, E = {(i, j) | 1 ≤ i < j ≤ n}. Pursuing the reasoning in [START_REF] Cremers | [END_REF]] and [START_REF] Seghers | [END_REF], we express the shape probability in terms of co-occurrence probabilities. Using the Bayes rule, we can for example write:

p (x 1 , . . . , x n ) = p (x 1 ) p (x 2 . . . , x n |x 1 ) = p (x 1 ) p (x 2 |x 1 ) p (x 3 . . . , x n |x 1 , x 2 ) = . . . = p (x 1 ) p (x 2 |x 1 ) p (x 3 |x 1 , x 2 ) • • • p (x n |x 1 , . . . , x n-1 ) . (2.1)
First, seeking a representation that does not depend on an absolute position, the probability of co-occurrence of a pair of control points can be then written as the probability of their dierence vector, or:

p (x i , x j ) = p (x i -x j ) . (2.2)
This hypothesis also leads to a uniform probability distribution for the absolute points positions or, a constant p (x i ), for all 1 ≤ i ≤ n. It also means that the conditional probabilities are equal for any pair of landmarks, [START_REF] Seghers | [END_REF]] proceed similarly to induce independence with respect to absolute location.

or p (x i |x j ) = p (x j |x i ), ∀ 1 ≤ i, j ≤ n.
We will consider the second order approximation of the shape probability, assuming that the positions of two landmarks is independent from any third landmark, or:

p (x i , x j | x k ) ≈ p (x i , x j ) , ∀k = i, j . (2.3)
Then, applying this property to equation (2.1) by choosing particular point pairs leads to:

p (x 1 , . . . , x n ) ∝ p (x 2 |x 1 ) p (x 3 |x 2 ) • • • p (x n |x n-1 ) .
(2.4) Equation ( 2.4) is one particular choice of decomposing the joint probability. This choice corresponds to specifying a spanning tree of the graph G.

Considering all such possibilities is equivalent to performing the previous simplication using all possible spanning trees of G. Let us multiply all these resulting equations, we obtain thanks to the simplication (2.3):

p (x 1 , . . . , x n ) τ (G) ∝ i<j p (x i , x j ) ν , (2.5)
where τ (G) is the number of all possible spanning trees of G, and ν is the number of appearances of each pair in all these possible equations. This number is obviously the same for each pair in a complete graph. Let us simplify the above equation using double counting. A spanning tree is a sequence of n -1 elements (or edges). Then the total number of elements used to form all the spanning trees is (n -1) × τ (G). There are n(n-1)

2 edges in a complete graph, and each edge appears ν times in all the spanning trees, which makes a total of n(n-1) 2

× ν elements. Then, ν = τ (G) × 2 n
. This gives a simple form to the product (2.5), or:

p (x 1 , . . . , x n ) ∝ i<j p (x i , x j ) 2 n
.

(2.6)

We would like to enforce invariance properties to our shape model, with respect to translations, rotations, and global scale changes. Such similarityinvariance characteristics enable to dismiss object changes due to the pose variations and permit to learn the intrinsic variability of the object. Discarding the transformation parameters in the inference step from the unknowns is also an advantage that makes an invariant shape model more tractable and easier to use.

Invariance with respect to rotation is equivalent to ignoring any orientation direction, that to consider the Euclidean norm of the dierence vector, or:

p (x i , x j ) = p ( x i -x j ) .
(2.7)

Eventually, let us consider d = 2 n(n-1) i<j x i -x j the mean Euclidean distance between the control points. This mean chord length can be considered as an estimate of the shape scale. Hence, dividing the chord lengths by d gives a representation that is invariant to scale changes as well, or:

p (x i , x j ) = p x i -x j d .
(2.8)

Strictly speaking, d depends on all the control points. Therefore, the ambiguity carried by equation ( 2.8) will be explained later in the algorithmic sections.

By combining equations (2.6) and (2.8) we dene the graph probability as:

p (x 1 , . . . , x n ) ∝ i<j p x i -x j d 2 n , (2.9)
Upon taking the negative log of the above equation, the graph energy is dened, up to an additive constant ν , as:

E shape (x 1 , . . . , x n ) ≈ - 2 n i<j log p x i -x j d + ν .
(2.10)

Let us consider now a set S = {s 1 , . . . , s m } of m instances of the object, where each example is represented using n control points, i.e. s u = {x u 1 , . . . , x u n }, ∀u ∈ {1, . . . , m}. Hence, ∀i ∈ {1, . . . , n}, the set X i = {x 1 i , ..., x m i } represents instances of the i th control point of the shape. In practice, this training set is obtained by manually labeling the landmarks for each instance of the shape, or by deducing the landmarks from the registration between a labeled shape and a set of non-labeled shapes. Note that an alignment of the shape is not needed before the learning phase since our representation is intrinsically invariant to similarity transforms.

Then, given a statistical model, we learn from the training set the probability density distributions of the relative positions of the control points

p ij ≡ p (x i , x j ). These n(n-1)
2 densities enable us to describe the information contained in the training set. In practice, we used a Gaussian kernel density estimation to compute the probability distributions.

One can question the use of the NCL representation for shape modeling.

Recall the observation made in section 2.1.3 with respect to the ability of EDM to model shape in an isometric-invariant manner [START_REF] Dattorro | Convex Optimization & Euclidean Distance Geometry[END_REF]]. The scaling factor d that we add simply extends the result to the similarity-invariant case.

However, this representation suers from redundancy and could be computationally expensive and hard to optimize during inference. We describe in the following the learning technique we use to remove the redundancy from our model, to obtain a sparser representation. But before proceeding, we discuss some geometric aspects that will enlighten the redundancy encapsulated by the complete graph.

Graph Rigidity and the Exact Matching Problem

This section is built on the work of [Caetano et al. 2006]. We rst recall the result that they provide, and then extend it to a more general case that is of particular interest in our modeling.

The Point Pattern Matching Problem

Let P s = (p 1 , . . . , p s ) a set of s "source" points in R n , and Q t = (q 1 , . . . , q t ) a set of t "target" points in R n where in general s ≤ t. The point pattern matching problem consists in nding a mapping m : P s → Q t that minimizes a distortion measure between P s and m(P s ). This can be formulated mathematically as minimizing:

D Ks (m) = s i=1 s j=1 ρ(d s ij , d t m(i)m(j) ) , (2.11)
where ρ is a dissimilarity measure, and

d s ij = p i -p j , d t ij = q i -q j are
Euclidean distances. In general, due to the noisy data, it is not possible to achieve a perfect matching, i.e. an exact matching. An exact matching is a mapping m such that d s ij = d t m(i)m(j) for all 1 ≤ i, j ≤ s. Hence, in the perfect matching case problem, ρ is chosen to be a binary function that checks the equality of its arguments. The pattern point matching can be written as a weighted complete graph matching, where the edge weights correspond to the Euclidean distance between the embedded vertices in R n . In the general case, if only a subset E s of the pairs of source points in P s are used to evaluate the matching distortion, this corresponds to the graph matching problem formulated on the graph G s = (P s , E s ), or minimizing:

D Gs (m) = (i,j)∈Es ρ(d s ij , d t m(i)m(j) ) .
(2.12)

In the remainder of this section, we will use indierently the terms graph or graph embedding.

Global Graph Rigidity [Caetano et al. 2006] dene the global rigidity of a graph in an intuitive manner.

Denition 2.5. A graph G = (V, E) in R n is said to be globally rigid if the lengths of the edges in E determine uniquely the lengths of the edges of the complement graph (the "missing" edges).

Then, they show that a k-tree with all base k-cliques in general position in R k-1 (i.e. do not lie in a k -2-dimensional vector subspace) is globally rigid in R k-1 . They use this property to prove that in the perfect matching case in

R k-1 , the k-tree model is equivalent to the complete model. Theorem 2.2. In the exact matching case in R k-1 , if G s = (P s , E s ) is a k-tree,
then a mapping m which minimizes D Gs (m) also minimizes D Ks (m).

We now extend this property to a larger case of graphs that will appear to be useful later in the manifold learning part. To this end, we invoke the following result, the proof of which is reported in Appendix C: Lemma 2.1 ([Cai & Maray 1993]). Every k -tree T with k ≥ k ≥ 1 possesses a spanning k-tree T with the property (P) that every k -clique of T contains a k-clique of T .

Then, it follows that the property proven by [Caetano et al. 2006] for the k-trees in R k-1 is generalizable. We suppose that all base k-cliques are in general position in R k-1 . Then: Proposition 2.1. In the exact matching case in

R k-1 , if G s = (P s , E s ) is a k - tree, with k ≥ k, then a mapping m which minimizes D Gs (m) also minimizes D Ks (m).
Proof. The mapping m minimizes D Gs (m). Then for all (i, j)

∈ E s , d s ij = d t m(i)m(j) . G s being a k -tree, with k ≥ k, then (lemma 2.1) there exists a k-tree G k s = P s , E k s , such that E k s ⊂ E s . Subsequently, for all (i, j) ∈ E k s , d s ij = d t m(i)m(j) , which is equivalent to saying that the mapping m minimizes D G k s (m). G k s being a k-tree in R k-1 , then (theorem 2.
2) the mapping m also minimizes D Ks (m), which completes the proof.

To summarize, a k-tree with all subsets of size n + 1 of the base k-cliques in general position in R n , with k > n, is equivalent to the complete graph model in the exact matching case problem. This result shows that such structures encode sucient information and dependencies to correctly model the shape. It also suggests that the complete graph model is redundant, and hence can be made sparser. Moreover, as [Caetano et al. 2006] show in their experiments, dierent k-tree models lead to varying accuracy in the non exact graph matching case. For this reason, a method of selecting a particular ktree with optimal properties is advantageous in their case. In the next step of our modeling we tackle the redundancy problem, and we build a k-tree graph structure where the parameter k is learned from the data, and where the particular structure of the graph is related to the training shape population.

2.2.3

Removing Redundancy from the Model

The task of eliminating the redundancy in the model, while preserving its ability to represent the data, is related to the minimum description length principle on one hand, and can be thought as a spectral clustering problem on the other hand. We aim to obtain as compact a model as possible assuming that the high dimensional data space can be approximated by a lower dimensional embedded manifold, which reduces the dimension of the problem signicantly. In other words, we would like to probe the population of training shapes for statistical dependencies between the control points variations.

Shape maps [START_REF] Langs | Modeling the Structure of Multivariate Manifolds: Shape Maps[END_REF] handle precisely these two aspects, and are learned from the data in a way closely related to the diusion maps [Coifman & Lafon 2006], but using the compactness of models that describe subsets of the entire data instead of the spatial distances or similarities between individual points. Therefore, we compute the shape map of the control points, using the training set {X 1 , . . . , X n }, and then we cluster the control points according to their mutual shape map distances. We start rst by reviewing the theoretical foundations of diusion maps.

Diusion Maps

Diusion maps is a technique based on diusion processes that aims to dene a metric, namely the diusion distance, that measures meaningfully the connectivity of data sets, with respect to their geometry. An ecient representation of the geometric structures complexity is given in [Coifman & Lafon 2006],

through the eigenfunction of Markov matrices. Formally, let X be the (dis-crete) data set, and let κ : X × X → R be a symmetric and positivity preserving kernel. Then a reversible Markov chain on X is constructed, with η(x, y) = κ(x,y) δ(x)

being the anisotropic transition kernel, and δ(x) = y κ(x, y).

The associated Markov matrix P (with entries η(x, y)) encode geometric information about the data and permits its analysis in a multi-scale fashion. The kernel η t (x, y) (which corresponds to the entries of the t th power of P ) represents the transition probability from the node x to the node y in t time steps.

The diusion process reects random walks on the Markov chain according to these transition probabilities. After the diusion process, the posterior distribution u → η t (x, u) characterizes the connectivity of the node x with respect to the other nodes. A weighted L 2 distance on these distributions denes the diusion distance between two nodes x and y at time t, or:

D t (x, y) 2 = η t (x, .) -η t (y, .) 2 L 2 = u (η t (x, u) -η t (y, u)) 2 π(u) , (2.13)
where

π(u) = δ(u) v δ(v)
is the probability of u in the unique stationary distribution. The larger the number of high probability paths of length t connecting the nodes x and y, the lower will be D t (x, y). Then, the authors relate the diusion distance to the eigendecomposition of the operator P . Let

λ 1 ≥ • • • ≥ λ i ≥ .
. . be the eigenvalues of P , and let {Ψ 1 , . . . , Ψ i , . . . } be its eigenfunctions. Then, a tractable expression of the diusion distance is given by:

D t (x, y) 2 = i λ 2t i (Ψ i (x) -Ψ i (y)) 2 ≈ N i=1 λ 2t i (Ψ i (x) -Ψ i (y)) 2 , (2.14)
and a good approximation is obtained by considering the top N eigenfunctions.

Shape Maps

Shape maps [START_REF] Langs | Modeling the Structure of Multivariate Manifolds: Shape Maps[END_REF]] bridge the gap between behavior analysis and manifold learning. It uses the diusion map technique with an appropriate kernel κ. Let us rst recall the notations of section 2.2.1.

The data are represented by a graph, where the vertices/control points are V = {x 1 , . . . , x n }. The pool of training data is {X 1 , . . . , X n }, where X i = {x 1 i , ..., x m i } represents instances of the i th landmark of the shape.

Shape maps use the description length (DL) principle to assess the mutual coherence of the landmarks behavior. The authors compute the compactness of small sub-models of landmarks subsets, that capture the non-rigid local variations, in order to dene a measure of similarity between landmarks, in terms of statistical dependence. The description length consists of the cost L of transmitting a model M and the data D encoded with the model:

L(D, M) = L(M) + L(D|M).
For each pair of landmarks, (x i , x j ), and for a dened kernel size s,

µ s (i, j) = min M (L(M)|(i, j) ⊂ M & |M| = s) (2.15)
is the minimum description length of models encompassing the two landmarks and s -2 other landmarks. The considered sub-models were Gaussian, and were determined using PCA on the corresponding s landmarks after alignment of the data. If x i and x j belong to a cluster of coherently behaving landmarks, then this will be detected by a low value of µ s (i, j). In practice, we considered s = 3 and evaluated the DL of all sub-models of this size.

Next, a Markov chain that reects the model structure is created, with a node corresponding to every landmark x i , through the denition of the kernel κ(i, j) = e -µs(i,j) ε . The diusion maps technique provides an embedding Ψ t : V → R N of the landmarks in an N dimensional Euclidean space:

x i → Ψ t (x i ) =     λ t 1 Ψ 1 (i) . . . λ t N Ψ N (i)     .
(2.16)

Thereby, the diusion distance in equation (2.14) can be expressed as:

D t (x i , x j ) = Ψ t (x i ) -Ψ t (x j ) .
(2.17 In terms of the joint model, low distance values suggest high compactness for a representation that encompasses landmarks x i and x j . We considered in practice a time step t = 1. Having obtained these behavior-explaining distances, we would like to cluster the control points accordingly. A recent clustering algorithm [Komodakis et al. 2008a] was used for this nal task, and is described in the section 2.2.4. The obtained clusters reect the interdependencies between the control points, and refer to the parts of the object that have highly-correlated relative displacements.

Unsupervised Clustering Using Linear Programming

Clustering refers to the process of organizing a set of objects into groups, where the members of each group are as similar to each other as possible. More formally, a common denition for clustering is the following one: suppose we are given a set of objects V = {v 1 , . . . , v n } endowed with a distance function d that can measure the dissimilarity between any two objects v i , v j ∈ V. In such a case, the goal of clustering is to choose K objects, say, {c 1 , c 2 , . . . , c K } from V (these will be referred to as the clusters centers), so that the obtained sum of distances between each object and its closest center is minimized, or:

min c 1 ,c 2 ,...,c K v i ∈V min c k d(v i , c k ) .
(2.18)

A common drawback of many popular clustering techniques (such as the Kmeans algorithm) is that they need to be given the number K of clusters beforehand. However, this is problematic as this number is very often not known in advance. To address this issue, we will let this number vary as well and change the objective function of clustering so as to assign a penalty (denoted by d(v i , v i )) whenever an object v i is chosen as a cluster center, or:

min K,c 1 ,c 2 ,...,c K v i ∈V min c k d(v i , c k ) + c k d(c k , c k ) . (2.19)
Another very bad symptom of many clustering techniques is that they are particularly sensitive to initialization. For instance, the K-means algorithm (which is one of the most commonly used clustering techniques) is doomed to fail if its initial cluster centers happen not to be near the actual cluster centers. To address this very important issue, we have used a novel clustering method. The main idea behind our method is to rst formulate the clustering as a linear integer program as follows:

min n i=1 n j=1 d(v i , v j )x ij (2.20) s.t. n j=1 x ij = 1, ∀i (2.21) x ij ≤ x jj , ∀i = j (2.22)
x ij ∈ {0, 1}, ∀i, j

(2.23)

In the above formulation, the binary variable x ij (with i = j) indicates whether object v i has been assigned to cluster center v j or not, while the binary variable x jj indicates whether object v j has been chosen as a cluster center or not. It is then very easy to prove that the above linear integer program is actually equivalent to minimizing the objective function (2.19) of clustering. To this end, it suces to observe that (2.21) simply expresses the fact that each object v i can be assigned to exactly one cluster center v j , while (2.22) simply expresses the fact that if any object v i has been assigned to an object v j , then v j must be chosen as cluster center. To obtain an approximately optimal solution to the above integer program, we will then rely on rst solving its linear programming relaxation and then rounding the relaxed solution in an appropriate manner. More details about the formulation of the problem and its optimization are given in [Komodakis et al. 2008a]. In the validation section of [Komodakis et al. 2008a], it is added that a constant penalty cost, roughly set to the median of the distances d (v i , v j ) is used in the experiments.

We also considered the same penalty value for our tests.

In our case, the set of objects correspond to the control points {x 1 , . . . , x n }, and the distance d corresponds to the aforementioned diusion distance D t in section 2.2.3. As an example of application, in the previously introduced synthetic case (Fig. 2.4(c)), this method is able to nd the actual number of clusters used to generate the data, and clusters it correctly.

The Shape Model

The unsupervised clustering technique described in section 2.2.4 is applied using the aforementioned shape maps distances (section 2.2.3). The process output is a number k of clusters (in practice, k ≥ 3 in 2D, and k ≥ 4 in 3D), and a set of k cluster centers, chosen among the control points. Then, the nal step to build our shape model consists in dening the structure of the graph, in a principled manner.

A cluster of landmarks is dened as a set of control points of highly correlated behavior. Thus, such a group of points carries a high amount of redundant information in terms of shape variability. This observation drives us to factor out the redundancy using the obtained clusters, through the center landmarks. The cluster centers being the best representatives of their respective groups, we make the mild assumption that for a given control point, that is not a cluster center, the joint landmark-cluster information can be encoded as landmark-cluster center information. Hence, the landmarks will be considered to be conditionally independent, knowing the cluster centers. The global coherence of the model is kept thanks to the reference part, that is the cluster centers, that form a clique of size k. The obtained graph is a k-tree with a single base k-clique. It is called a k-fan, following [Crandall et al. 2005].

The novelty here consists in the method that denes automatically from the training data the number of clusters and their centers. To each one of these pairs (x i , x j ) we associate a probability density distribution p ij learned from the training set as previously stated in section 2.2.1. Our k-fan graph, being chordal, benets from the simplicial elimination ordering (SEO) property, and permits theoretically the development of a Bayesian framework for segmentation for instance. Dynamic Programming-based approaches ensure nding a global optimum. However, the maximum clique size of a k-tree being k + 1, these polynomial time algorithms are not usable in practice for huge memory requirements and extreme slowness. For these reasons, we resort to approximate optimization techniques for inference on the k-fan graph.

We will denote the k-fan graph by G k = (V, E k ). Without loss of generality, we can assume that the base k-clique corresponds to the k rst nodes in V.

Let R k be the subset of edges in E k that correspond to the base k-clique. The k-fan model assumes that a node x i is conditionally independent from a node

x j (with k + 1 ≤ i < j ≤ n) knowing the base k-clique. To derive the joint probability of our shape model, we use the same reasoning presented in section 2.2.1, by taking into account this conditional independence assumption.

Similarly to equation (2.1), using the Bayes rule, we have:

p (x 1 , . . . , x n ) =p (x 1 ) p (x 2 |x 1 ) • • • p (x k |x 1 , . . . , x k-1 ) p (x k+1 |x 1 , . . . , x k ) • • • p (x n |x 1 , . . . , x n-1
) .

(2.24)

The aforementioned independence assumption allows the following simplication of equation ( 2.24):

p (x 1 , . . . , x n ) =p (x 1 ) p (x 2 |x 1 ) • • • p (x k |x 1 , . . . , x k-1 ) p (x k+1 |x 1 , . . . , x k ) p (x k+2 |x 1 , . . . , x k ) • • • p (x n |x 1 , . . . , x k ) .
(2.25)

Then, by applying the second order assumption in (2.3), we can for example write:

p (x 1 , . . . , x n ) ∝ subset of R k p (x 2 |x 1 ) p (x 3 |x 2 ) • • • p (x k |x k-1 ) p (x k+1 |x 1 ) p (x k+2 |x 1 ) • • • p (x n |x 1 ) subset of E k R k . (2.26)
Equation (2.26) corresponds to specifying a spanning tree of the base k-clique, and choosing for each node that does not belong to the base k-clique a neighbor x i , with 1 ≤ i ≤ k. Upon multiplying all the equations that result from this construction, we obtain the following expression (which is a modied version of the equation (2.5)):

p (x 1 , . . . , x n ) ξ ∝ (i,j)∈R k p (x i , x j ) ν (i,j)∈E k R k p (x i , x j ) ν .
(2.27) Note that the pairs (i, j) with k + 1 ≤ i < j ≤ n do not appear in this expression due to their conditional independence assumption. A node in the base k-clique has a degree n -1 (i.e. n -1 neighbors), whereas a node that does not belong to the base k-clique has a degree k. Therefore, ν and ν are dierent in general. Let ζ be the number of spanning trees of the base k-clique.

Then we clearly have that ξ = ζk n-k , since there are n-k non-reference nodes having a degree equal to k. Let η be the number of appearances of an edge in R k in the total number of spanning trees of the base k-clique. Then, as showed in section 2.2.1, η = 2 k ζ. Now, it is easy to count the number of edge appearances ν and ν in the product (2.27). On one hand, we can see that Algorithm 1 Build the NCL Shape Model 1: Apply shape maps to the training data.

2: Cluster the control points x i using the LP-based unsupervised clustering algorithm, through their pairwise diusion distances D t (x i , x j ). The output of this step is k clusters, and k prototype centers. 3: Build the k-fan graph G k = (V, E k ) that represents the data by linking each landmark x i to all the k cluster centers x j . 4: For every pair of control points (x i , x j ) such that (i, j) ∈ E k , learn a Gaussian kernel density estimation of the probability density function

p (x i , x j ) = p x i -x j d . ν = ηk n-k , which leads to ν = 2 k ξ. On the other hand, ν = ζk n-k-1 which gives ν = 1 k
. Then we can write equation ( 2.27) more explicitly:

p (x 1 , . . . , x n ) ∝ (i,j)∈R k p (x i , x j ) 2 k (i,j)∈E k R k p (x i , x j ) 1 k . (2.28)
To simplify the notations, we introduce the weighting γ ij dened by:

γ ij = 2 k if (i, j) ∈ R k 1 k if (i, j) ∈ E k R k (2.29)
The shape energy in equation (2.10) can now be approximated as:

E shape (x 1 , . . . , x n ) ≈ - (i,j)∈E k γ ij log p x i -x j d .
(2.30)

Algorithm 1 summarizes the construction of our shape model. Recall that the thread leading to our shape model is threefold:

• Geometry: the k-fan MRF has an optimal geometric property in the exact matching case, that is related to its global rigidity.

• Manifold learning: the k-fan MRF particular structure, and the number k of fans are deduced from the population of training shapes.

• Graph theory -optimization: although with very limited usefulness, due to the curse of dimensionality, the k-fan MRF benets from the advantages of the larger class of chordal graphs. The projection of the control points on the rst 3 shape map dimensions. The density of the landmarks is estimated by a Gaussian kernel and is color-coded.

(c) Control points clustered in 9 clusters: centers are represented by squares.

Experimental Validation

Drawing Samples from the Learned Distribution A rst application that is usually done after building a statistical shape model consists in generating synthetic shapes. This can be done by estimating the probability distribution over a grid, and then drawing sample from this distribution. We formulate the probability estimation similarly to [Felzenszwalb & Huttenlocher 2005]. The probability of our shape model is given by equation (2.28). We will rst compute the marginal probability of the reference part (the base k-clique) p (x 1 , . . . , x k ). Then we will sample a location for the landmarks belonging to the base k-clique. Next, we will sample from p (x i |x 1 , . . . , x k ) a location for the non reference landmarks. We start by marginalizing (2.28) over the non reference part: The exponential number of operations described by the above equation can be reduced to polynomial, by rearranging the terms following a SEO, or:

p (x 1 , . . . , x k ) ∝ k+1 • • • n p (x 1 , . . . , x n ) . (2.31) (a) (b) (c) (d)
p (x 1 , . . . , x k ) ∝ k+1 k i=1 p (x k+1 , x i ) γ i,k+1 • • • n k i=1 p (x n , x i ) γ in   (i,j)∈R k p (x i , x j ) γ ij   .
(2.32)

The previous expression is fairly easy to obtain, as the depth of the dual tree corresponding to the k-fan (or its junction tree) is equal to 2. The conditional probability of a landmark knowing the reference points is:

p (x i |x 1 , . . . , x k ) ∝ k j=1 p (x i , x j ) γ ij . (2.33)
Sampling on a grid of size m runs in O(nm k+1 ) polynomial time, and requires O(nm k ) memory space. This is why in practice it is unfeasible to perform such a sampling. We resort to another method to assess the ability of our method to capture the shape information, by applying a priori-constrained perturbations to a randomly located set of random points. Let (p 1 , . . . , p n ) be a set of random points in R N . We propose to verify the information encoded by our shape prior through constrained displacements of this set of points. Formally, according to the shape model refers to the displacements that minimize:

let D i = {d i1 , . . . , d im }, 1 ≤ i ≤ n be
-

(i,j)∈E k γ ij log p x i + d ii -x j -d jj d .
(2.34)

This quadratic assignment problem min i ,j

(i,j)∈E k V ij (l ii , l jj
) is NP-hard in the general case, and can be solved approximately using [Komodakis et al. 2008b] or [Kolmogorov 2006]. An implementation of these methods is available online 3,4 .

The minimization problem (2.34) can be solved iteratively. The solution of the iteration i initializes the iteration i + 1. After a few iterations, the shape does not evolve any more.

In practice, applying shape prior constraints to an initial set of random control points leads to a "mean" instance of the learned object. The use of such prior in the segmentation framework is explored in the following. We used all these landmarks to train our shape model, and we obtained the shape maps embedding in Fig. The last step of the validation consisted of applying the shape prior to a set a randomly located points. We used Gaussian distributions to models the pairwise probability functions. As expected, the shape prior guides the landmarks towards a corpus callosum shape. Some iterations of this contour evolution are presented in Fig. 2.8.

Right Lung and Left Lung dataset We applied our method to the shape modeling of the right and left lungs from 2D radiographic (X-ray) images. We used a publicly-available dataset Gold standard segmentation masks are hence available as well as corresponding landmark positions lying on the contour. We used the 44 available landmarks to learn the NCL shape model. In terms of computational complexity, we executed the "perturbation moves" algorithm 50 successive times using a 3 GHz dual-core processor. We report the running time, using respectively the complete graph and the k-fan graph, for the hand example and the left ventricle example in the Table 2.1.

(a) (b) (c) (d)

Introduction

The use of prior models is very popular in image segmentation. The most common approach consists of associating the model with the observations and then seeking for the optimal model parameters, that is a compromise between the observed image and the prior constraints. In this section, we This chapter is structured as follows: we review in section 3.1.1 the optimization techniques related to pairwise Markov random elds (MRFs). In section 3.1.2, we detail dierent segmentation techniques, and describe specifically the graph-based methods. We will particularly discuss the approaches that combine MRFs with shape priors, showing that it is is a tedious procedure in general, due to the fact that graph-based methods assume weak connectivity while shape priors often require global knowledge of the shape parameters.

In an attempt to incorporate a statistical shape prior model within the MRF framework, we propose in section 3.2 a segmentation method that combines our learned prior knowledge with a region-based data term, through the use of Voronoi diagrams. The resulting algorithm features an object contour (or surface) that evolves towards the solution, in the spirit of deformable contours and surfaces. Experiments of this method are shown in a 2D computer vision application, and in a 3D medical image analysis setting. In order to illustrate the exibility of our framework, we detail in section 3.3 a slightly dierent approach that shares ideas with graph matching and point pattern matching. In this algorithm, the data term is related to the responses of classiers based on distinctive image feature vectors. In both cases, we show that our statistical shape model guarantees regular and robust segmentation.

Discrete MRFs and Optimization

MRFs were introduced to the computer vision domain by [START_REF] Geman | [END_REF].

Since their seminal work, MRFs have been successfully applied to the modeling of a wide range of vision applications, such as image restoration, image completion, segmentation, stereo vision, tracking and motion estimation, registration, etc [Boykov & Funka-Lea 2006].

These applications often seek the estimation of the model variables with respect to some observations. This is often formulated as a maximum a posteriori (MAP) estimation problem, and under mild assumptions, it leads to minimizing the following MRF energy:

E(l) = p∈V V p (l p ) + (p,q)∈E V pq (l p , l q ) , (3.1)
where G = (V, E) is the considered MRF, V is the set of graph variables, E is a local neighborhood system associated with this graph that can be spatially varying, and l p ∈ L is a possible state of the variable p ∈ V. This problem is also interpreted as a labeling problem in the discrete case, where the unary potentials V p (l p ) indicate/measure the quality of the labeling with respect to the observations for a given node, and the pairwise potentials V pq (l p , l q ) aim to establish coherence between the labeling of connected nodes. Minimizing such function is a challenge in the general case, as this problem is proven to be NPhard. However, numerous optimization techniques were developed to tackle this problem from dierent points of view. For instance, some methods guarantee an exact optimal solution of the problem (3.1) for a very restricted set of cases. Others provide nearly optimal solutions, with predened condence, if some constraints are satised. We present in the following an overview of these methods.

Early Methods

Iterated conditional modes (ICM), simulated annealing (SA), highest condence rst (HCF) and mean eld annealing (MFA) are among the rst algorithms considered to solve the aforementioned MRF minimization problem.

ICM [Besag 1986] is a greedy approach that iteratively updates the labels of the MRF nodes. At each iteration, a node is allowed to change its label to the one that provides the local largest decrease of the objective function assuming the remaining of the labels constant from the previous iteration. These updates are repeated for all nodes until convergence to a local minimum, i.e.

until no label update gives any decrease of the energy. This method is sensitive to initialization. HCF [START_REF] Chou & Brown | [END_REF]] is a deterministic algorithm that attempts to improve the ICM by updating the nodes in a certain order.

It introduces an additional label to describe an "uncommitted" state. At the beginning of the algorithm, all nodes are "uncommitted". The algorithm denes the so-called stability for each node, and at each iteration, only the node with the smallest stability is allowed to commit, or update its label. Once a node is committed to a label, the potentials of the connected nodes are also updated and decisions for them can be taken with more condence. SA [START_REF] Kirkpatrick | [END_REF][Kirkpatrick et al. ][erný 1985] is also an iterative algorithm, where a temperature parameter T plays an important role (in analogy with the slow cooling/annealing phenomenon in metallurgy). Iteratively, the SA heuristic choses whether to stay in the current state or to move to a neighboring one in a probabilistic manner. If moving to a new state reduces the energy, then the move is accepted. Otherwise, the move is accepted with a certain probability that depends on the temperature T , and decreases over time, i.e. as the temperature tends to zero. Intuitively, for large values of T , the system is allowed to move with certain freedom, such that it is able to escape local minima. This freedom decreases when T is close to zero, and the probabilistic moves become almost deterministic. Asymptotically, SA nds the global optimum, which in practice would require an innite time. Note also that SA implementations converge to a local minimum. Moreover, [Greig et al. 1989] gave binary labeling examples where the solution provided by SA is far from the optimum. MFA [Bilbro et al. 1989] is an algorithm that is similar in spirit to SA. It proceeds however by estimating the mean eld of each node at decreasing temperatures, rather than with a sampling scheme and probabilistic moves. The mean eld of the MRF converges to the MAP when T tends to zero. This popular minimization algorithm resorts in practice to approximations as it needs to estimate the partition function, which is computationally intractable. This is another instance of algorithms with asymptotic optimal properties, which converge however to local minima when used in real-world applications. We turn now to present another class of methods that belongs to combinatorial optimization techniques and exhibits eciency properties as well as the ability to achieve a global optimum under certain conditions.

Graph-Cut -Binary Case [Greig et al. 1989] were the rst to use a s/t graph cut algorithm in the context of MAP estimation of a MRF in computer vision. They successfully applied it to the restoration of binary images and showed that it outperforms SA. A decade later, the rediscovery of graph cut techniques in the computer vision community started a new research trend. This was due to the seminal work of [START_REF] Boykov | [END_REF] where an ecient implementation of the maxow/min-cut problem was introduced in the context of vision. These ideas were popularized in [Boykov & Jolly 2001] where the use of the graph cut technique for interactive n-dimensional image segmentation was introduced.

In this subsection, we start rst by describing the basics of graph cut algorithms. In order to simplify the presentation, let us rst consider a simple case where the label set consists of two labels.

The binary labeling problem is represented by a directed graph where each node represents a variable, and where edges connect neighboring nodes.

These nodes are called internal nodes, and the corresponding edges are called internal links or n-links. In addition, the graph contains two terminal nodes called the source s and the sink t. These two special nodes correspond to the two classes of the binary labeling. The graph contains also terminal edges, or t-links, that connect each internal node with both terminal nodes. Each edge is associated with a cost, or a capacity. The cost of a n-link denes the penalty of assigning the corresponding two nodes to dierent classes, whereas the cost of a t-link represents the penalty of assigning the node to the corresponding terminal. For the case of the MAP estimate of a binary MRF, dened by (3.1), the n-link costs correspond to V pq (l p , l q ) and the t-links to V p (l p ), where L = {s, t}.

An s/t cut on such a graph partitions the nodes in two disjoint classes S and T such that S contains the source terminal s and T contains the sink terminal t. The cost of such a cut is equal to the sum of the costs of directed edges (p, q) such that p ∈ S and q ∈ T . Hence, a MAP binary labeling is equivalent to a cut with a minimum cost, or a min-cut. The theorem of [Ford & Fulkerson 1962] proves that the min-cut problem can be equivalently solved by nding the maximum ow from the source s to the sink t, through the capacitated edges. It also states that this max-ow saturates a set of edges that correspond to a min-cut. The max-ow and min-cut values are equal. Dierent low (polynomial) complexity algorithms for solving the mincut/max-ow problem exist, and can be classied in two major classes: "pushrelabel" methods [START_REF] Goldberg | [END_REF] and "augmenting paths" methods [Ford & Fulkerson 1962]. [START_REF] Boykov | [END_REF] give a description of these methods, and compare their performance on a set of computer vision applications. They also present a modied "augmenting paths" algorithm that produces the best results in their experiments. As far as the quality of the obtained solution is concerned, [Kolmogorov & Zabih 2004] discuss the set of energy functions that can be minimized using graph cuts. Hence, a global optimum is only guaranteed when the energy function (3.1) is submodular,

meaning that V pq (0, 0) + V pq (1, 1) ≤ V pq (0, 1) + V pq (1, 0), ∀(p, q) ∈ E.
For non-submodular functions, the problem remains NP-hard in general. However, [START_REF] Kolmogorov | [END_REF]] review a method called quadratic pseudo-Boolean optimization (QPBO) [START_REF] Hammer | [END_REF]] that provides partial optimal solutions in the arbitrary case. Essentially, the algorithm builds an appropriate s/t graph, where each variable p ∈ V is represented by two nodes, z and z. Ideally, the label of z should be the negation of the labeling of z, i.e. l z = 1 -l z . The new energy corresponding to the constructed graph is submodular, and is hence minimized using min-cut/max-ow. A consistent labeling of z and z (i.e. z ∈ S and z ∈ T , or the converse) provides a labeling for p. The variables corresponding to an inconsistent labeling of z and z remain unlabeled. This algorithm provides hence a partial solution, but with the guarantee that the labeled variables belong to the optimal solution.

One can put the importance of graph cut methods into perspective as they are limited to the binary case, and come with restrictions on the energy function. However, as we shall see in the following, the binary case is closely related to the multi-label setting.

Graph-Cut -Multi-label Case

Several works in the late 90's showed the high relevance of s/t graph cut with respect to the multi-label problem. [Roy & Cox 1998] were the rst to use an s/t graph cut for the general labeling case. Building on their work, and using appropriate graph constructions, [Boykov et al. 1998] proposed a method to exactly minimize the problem (3.1) when the pairwise potentials are linear functions of the labels. [Ishikawa 2003] shows subsequently that convex pairwise interactions can be handled as well. [Boykov et al. 2001] further generalize the result to the case where V pq is a metric, showing that a good approximation of the optimal solution can be found within guaranteed bounds.

The algorithms that they proposed, α-expansion and α -β-swap are ecient and perform well in practice, even in the non-metric case [Rother et al. 2005] (yet losing optimality guarantees). We give now some details about these algorithms.

The proposed expansion move and swap move algorithms in [Boykov et al. 2001] try to iteratively improve the current labeling. In an α-expansion move, any node is only allowed to change its label to α. In an α -β swap, only nodes with labels α (respectively β) are allowed to change their label to β (respectively α), and β (respectively α) only. Both algorithms nd in each iteration the move that minimizes (locally) the energy.

The number of these moves being exponential, the authors translate this minimization problem into a min-cut / max-ow problem in an appropriately constructed graph. Hence, during each of these iterations, an s/t graph cut is performed. A complete set of iterations over all the labels L is called a cycle. The algorithms are guaranteed to converge in a nite number of cycles to a local minimum. For α-expansion, this local minimum is within a known factor of the global optimum.

Loopy Belief Propagation

Belief propagation (BP) [START_REF] Pearl | Judea Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]] is an iterative algorithm that propagates local messages between the nodes of the MRF. We present here an intuitive description of its mechanism. At each iteration, a node sends messages to its neighbors, and receives messages from them. The algorithm keeps repeating these iterations until a xed point is reached, i.e. all the messages stabilize.

In each iteration, a node p sends a message {m pq (l q )} lq∈L to the node q.

This message expresses how likely the node q should be assigned a label l q according to the node p. This message considers two aspects. First, using the interaction term V pq (l p , l q ), the node p examines the likelihood of a label l q for the node q if it gets the label l p . Then, it assesses the likelihood of getting a label l p , based on the evidence V p (l p ) and the messages m rp (l p ) r =q it receives from its neighbors (other than q) about the likelihood of l p .

After exchanging such messages suciently, the messages converge. Then a belief b p (l p ) of assigning a label l p to a node p can be computed, as a combination of the evidence V p (l p ), and the messages m rp (l p ) that p received from its neighbors. One can distinguish between two kinds of BP algorithms:

the sum-product algorithm that is used to estimate marginals, and the maxproduct algorithm that gives MAP estimates. In both cases, BP is exact on acyclic graphs, or trees, and only needs to pass messages in a forward way (from the leaves to the root) and then in a backward way to converge.

However, when the graph contains loops, there are no such guarantees, and the algorithm may not converge [Yedidia et al. 2003] [Kschischang et al. 2001]. In practice, loopy belief propagation (LBP) has successfully been applied even with cyclic graphs.

Linear Programming Relaxation

The linear programming (LP) relaxation of the MAP-MRF problem (3.1) is a prominent approach to nd approximate solutions. The latter can be cast to an equivalent integer program:

min p∈V l∈L V p (l)x p (l) + (p,q)∈E l,l ∈L V pq (l, l )x pq (l, l ) , (3.2)
such that:

l x p (l) = 1 ∀p ∈ V (3.3) l x pq (l, l ) = x q (l ) ∀l ∈ L, ∀(p, q) ∈ E (3.4) l x pq (l, l ) = x p (l) ∀l ∈ L, ∀(p, q) ∈ E (3.5)
x p (.), x pq (., .) ∈ {0, 1} .

(3.6)

The binary variables x p (.) and x pq (., .) that are introduced in the above functions act as indicator functions, and specify which label is assigned to a node or a pair of nodes. Hence, x p (l) = 1 is equivalent to assigning label l to node p, and x pq (l, l ) = 1 is equivalent to assigning label l to node p and label l to node q. The constraint in equation (3.3) guarantees that each node is assigned one and only one label, and constraints (3.4) and (3.5) ensure consistency between the unary and pairwise binary variables.

The LP relaxation consists in relaxing the integer constraints (3.6) to the constraints x p (.) ≥ 0 and x pq (., .) ≥ 0. This LP relaxation was rst introduced by [START_REF] Schlesinger | [END_REF]], as pointed out in [Werner 2007]. The dual program that corresponds to the (primal) linear program can be written as:

max p∈V y p , (3.7)
such that:

y p ≤ min l∈L h p (l) ∀p ∈ V (3.8) y pq (l) + y qp (l ) ≤ V pq (l, l ) ∀(l, l ) ∈ L, (p, q) ∈ E , (3.9)
with h p (.) = V p (.)+ q (p,q)∈E y pq (.). The quantity h p (.) is called height variable.

Optimizing the obtained linear program using conventional techniques such as the simplex algorithm or interior-point methods is often not practical in the context of computer vision problems because of their very large scale. For this reason, the above-mentioned methods would be very slow. This motivated the development of approximate but ecient LP optimization techniques. These methods either use the primal-dual paradigm, or proceed in the dual domain to maximize a lower bound of the energy. We give in the following a description of the most prominent approaches.

Primal-Dual approaches and Fast-PD [START_REF] Komodakis | Nikos Komodakis and Georgios Tziritas. Approximate Labeling via Graph Cuts Based on Linear Programming[END_REF] apply the primal-dual principle to solve approximately the linear program (i.e.

the relaxed version of (3.6)), while providing suboptimality bounds. Their algorithm requires that the pairwise costs are positive (which is less restrictive than the metricity requirement of α-expansion for instance). By weak duality, the optimal cost of the primal-dual problem lies between the costs of any pair of primal and dual feasible variables. Therefore, the primal-dual principle states that controlling their cost gap by a factor f provides an fapproximation to the optimal integral solution. Moreover, feasible primal dual solutions that satisfy the relaxed complementary slackness condition provide an f -approximation to the optimal integer solution. Subsequently, the primaldual schema can be applied to nd a solution that is an f -approximation to the optimum. This schema is devised to generate pairs of integral primal and dual solutions until the last pair forms a primal and a dual feasible solutions that satisfy the relaxed primal complementary slackness conditions. [START_REF] Komodakis | Nikos Komodakis and Georgios Tziritas. Approximate Labeling via Graph Cuts Based on Linear Programming[END_REF] propose an algorithm that follows the primal-dual schema. They maintain during optimization an active primal labeling x p (.), and dual variables y pq (., .), such that the complementary slackness conditions are satised. To this end, they iteratively update the primal and dual variables to meet in particular the height complementary condition h p (x p ) = min l∈L h p (l). The update of the dual variables takes place in a selective manner. In each c-iteration, only c-label variables y pq (c) are changed, to satisfy the constraints as much as possible, which in turn would produce new h p (c) height values. The optimal update of the c-variables is found by solving a min-cut max-ow problem in an appropriately constructed capacitated graph. At this point, there might still be nodes violating the height complementary condition with respect to the label c. This is taken care of by updating accordingly the primal variables corresponding to those nodes.

After each update, the dual cost is guaranteed to decrease. The algorithm converges to a pair of primal-dual solutions that satisfy the complementary slackness conditions (when no update is done after |L| c-iterations), and hence provide an f -approximation to the solution. The value of f can be estimated from the pairwise costs, which are (only) assumed to be positive, as f = max pq max l =l Vpq(l,l ) min l =l Vpq(l,l ) . Besides the theoretical guarantees on the approxima-tion quality, per-instance suboptimality bounds are provided for each application. These bounds proved to be in practice much tighter than the theoretical ones. It also turns out that this Primal-Dual method is a generalization of α-expansion.

[ Komodakis et al. 2008b] further extend their method to obtain more efciency and provide a very fast primal dual algorithm, called FastPD. They essentially build upon their work in [START_REF] Komodakis | Nikos Komodakis and Georgios Tziritas. Approximate Labeling via Graph Cuts Based on Linear Programming[END_REF], by additionally making the link between the primal-dual gap and the number of path augmentations during the min-cut/max-ow optimization. Indeed, the speed of the algorithm depends heavily on the eciency of solving these max-ow problems, which in turn is related to the number of augmenting paths per max ow. They show that the primal-dual gap is an approximate upper bound on the number of such paths. A modication of their previous algorithm ensures that the primal-dual gap always decreases (and is never destroyed, even temporarily). This rst adjustment leads to a clear decrease in the number of augmentation paths and hence to a substantial speed up of the algorithm, while maintaining it optimality properties. Moreover, they design a max-ow algorithm that is adapted to their graph (that has few s-links) and provide an incremental graph construction method which yields an additional acceleration. In practice, FastPD performed 3 to 9 times faster than α-expansion.

Message-Passing Approaches Tree-reweighted (TRW) methods are another class of approaches that aim to approximately solve the MAP-MRF problem. While relying on the passing of messages similarly to BP, these methods are closely related to LP. Indeed, they attempt to solve the LP relaxation of the MRF, by maximizing the dual of the LP. [START_REF] Wainwright | [END_REF] introduced the tree-reweighted message passing algorithm (TRW-MP). They

show that by combining concave subproblems, they can provide a concave lower bound on the MRF energy. These subproblems can be chosen as trees.

They show that the message passing algorithm leads to a tight bound if and only if the collection of trees shares a common optimum. They also prove that the maximization of this lower bound (that is independent from the choice of the trees) is equivalent to optimizing the dual of the LP relaxation of the MAP MRF problem. However, TRW-MP does not have any guarantee of convergence, or increase of the lower bound. [Kolmogorov 2006] proposes a sequential version of the tree-reweighted algorithm, TRW-S, that has convergence guarantees. However, the author shows that the algorithm may converge to a point that is not the optimum of the LP lower bound. The optimal properties are only shown for the binary submodular case [Kolmogorov & Wainwright 2005].

[ Komodakis et al. 2007] introduced a new message-passing method based on dual decomposition (DD) [Bertsekas 1999], in order to address the abovementioned limitations. More explicitly, they aimed to solve the dual LP problem, i.e. maximize the lower bound. The dual decomposition technique consists of decomposing the Lagrangian dual of the intractable original MAP MRF problem, into less complex subproblems, thanks to a Lagrangian relaxation. [Komodakis et al. 2007] consider tree subproblems that can be easily solved. Then, the optimization (using the projected subgradient method)

amounts to an iterative "message passing" between the "master" problem and the "slave" problems. Iteratively, each easy subproblem is solved according to the updated MRF potentials, and the result is sent to the master problem that updates the potentials according to the new results. Moreover, the authors show that the Lagrangian relaxation problem is equivalent to the LP relaxation of the MAP MRF problem. Thus, they establish the link between their method, the TRW methods and LP. Unlike TRW approaches, this DD-based algorithm is guaranteed to maximize the lower bound of the LP relaxation. [Werner 2007] to impose necessary optimality conditions. When using this algorithm, the lower bound on the MRF energy does not decrease.

Tightness of the Relaxation

These iterations are repeated until there are no more inconsistent cycles. This can be thought of as a series of relaxations of gradually increasing tightness.

In practice, this algorithm leads to better optima. It especially performed well in dicult MRF problems, with respect to state-of-the art techniques.

[ [START_REF] Szeliski | [END_REF]] provide a comparison between dierent MRF energy minimization methods, using benchmark problems. The algorithms that are compared are ICM, max-product LBP, a LBP implementation derived from TRW-S, α-β-swap, α-expansion and TRW-S.

Image Segmentation Techniques

We briey sketch the most prominent image segmentation techniques in the literature. We note that in addition to the methods that we will detail in the following, there is another class of segmentation techniques that is mainly applied in medical image analysis, namely atlas-based segmentation. Often described as segmentation by registration, their principle consists of registering a reference image, which is called atlas, to a query image. A ground truth segmentation is associated with the reference image, and a segmentation of the query image is obtained by a warping operation, based on the resulting deformation eld. A more sophisticated approach that goes beyond the use of simple single-image atlas considers the multi-atlas segmentation problem.

Typically, every atlas is tted to the query image, and a nal segmentation is obtained in a multi-classier labeling manner by voting for instance. To capture the shape changes due to anatomical and morphological variations, more advanced atlas construction techniques were developed, and involved incorporating probabilistic information. Research about atlas-based segmentation has mainly focused on registration algorithms and atlas construction. A review of these dierent methods is beyond the scope of this work. However, the reader can be referred to a detailed description in [START_REF] Rohlng | [END_REF].

We give now an overview of deformable models and level sets for segmentation, active shape and active appearance models, and eventually graph-based segmentation approaches. We end this last part by discussing the endeavors that researchers made to introduce shape priors in the MRF framework, and commenting the approaches that are particularly related to our work. curve regularity and imposes smoothness constraints to the evolving contour, whereas the external energy stands for the image-driven forces that attract the curve to the image location with high gradient values. This 2D edge-based approach was then generalized to the 3D case in [START_REF] Terzopoulos | Constraints on deformable models: recovering 3D shape and nongrid motion[END_REF]].

This method is sensitive to initialization, and can get trapped in local min-ima. [START_REF] Cohen | [END_REF]] introduces an additional term in the energy that represents "balloon" forces. This inating (or deating) term aims to avoid the local minima by helping the curve to pass them. Another criticism that is commonly formulated to active contours is related to their inability to allow for topology changes. [START_REF] Delingette | General Object Reconstruction Based on Simplex Meshes[END_REF]] attempt to overcome this limitation by introducing the deformable simplex meshes. In such meshes, each vertex has exactly the same number of neighbors. This non-parametric surface representation is able to model objects of arbitrary topologies. The author proposes an algorithm that adapts the topology of an initial surface to the object to be extracted from the image. [Montagnat et al. 1999] use a 2-simplex mesh to segment 3D cylindrical echocardiographic images. This work is extended in [Montagnat et al. 2003] to handle 3D ultrasound image sequences. By using reparametrization tricks, [McInerney & Terzopoulos 1999] are also able to propose a topology-adaptable version of snakes. Alternatively, level sets [Osher & Sethian 1988] are an elegant framework that naturally handles topology changes.

The work of [START_REF] Malladi | [END_REF],

[ Kichenassamy et al. 1995] and [Caselles et al. 1997] benet from these aspects, as they introduce the geodesic active contours, which can be viewed as geometric (parameter-free) counterparts of classical snakes.

Opposite to the previously-mentioned methods, other variational approaches rely on image region characteristics. [Mumford & Shah 1989] dene a functional which allows to segment the image into homogeneous parts. They approximate the image by a piece-wise smooth function. The discontinuities of the approximating function are only allowed on a curve with its length being minimized within the functional. [Chan & Vese 2001] incorporate this approach in the level set framework, and design the so-called "active contours without edges". [START_REF] Zhu | [END_REF]] use statistical techniques of region growing by considering an MDL criterion. Their approach has also links to the snakes/balloon formulation. [START_REF] Yezzi | [END_REF] implement this method in the level sets framework.

[ [START_REF] Paragios | Coupled Geodesic Active Regions for Image Segmentation: A Level Set Approach[END_REF] introduce the geodesic active regions Bayesian formulation that combines image region statistics with edge-based terms within the level set framework. The use of both region and edge information was considered earlier in [Chakraborty et al. 1996] without beneting however from the power of the level-set framework. They apply the Green's theorem to the region term, such that is can be expressed on the boundary of the evolving contour. Their approach features a shape model that uses the probability distributions of the Fourier coecients representing the contour. However, their method depends on the initial conditions as it is optimized using a conjugategradient algorithm. The same idea of applying the divergence theorem was used by [START_REF] Unal | [END_REF] where an active polygonal deformable model was dened based on image region information.

Continuous eorts have been made since then to improve these methods by developing knowledge-based variational formulations. Representing shape using a signed distance maps (SDMs) was a natural way of achieving this goal within the level set framework (section 2.1.1). Linear shape models (section 2.1.1) as well as their non-linear kernel-based counterparts (section 2.1.1)

were also used in a variational framework.

Active Shape and Active Appearance Models

Active shape models (ASM) have become increasingly popular since their introduction by [START_REF] Cootes | [END_REF]]. This technique learns the shape of a class of objects using PCA (section 2.1.1). A shape s is hence approximated is assessed by the Mahalanobis distance. Proles of derivatives were also used in practice, and these dierent proles were often normalized. ASM is an iterative algorithm that uses a local search method, due to the large size of the search space in general. At an iteration t, the evolving shape s t in the image is expressed using its parameters b t and a similarity transformation T , or s t = T t (s + Eb t ). Then, for each landmark, the locally optimal candidate position along the normal to the contour is identied using the local appear-ance model (intensity prole). The output of this step is a vector of optimal displacements ds t . The transformation T t and then the shape parameters b t are updated accordingly (while applying limit constraints on b t ). Hence, a new valid instance of the shape s t+1 is obtained. These steps are repeated till convergence. The convergence depends on the initial conditions. Since ASM is a local search method, the initialization needs to be close enough to the solution. The robustness and the speed of the algorithm are improved through the use of coarse-to-ne search. For each level of the image corresponding pyramid, a dierent appearance model is built in the training stage.

Alternative ASM appearance models The exibility of the ASM framework allows the use of dierent features, which paves the way to substantial improvements. For instance, Gabor wavelets were used in [Jiao et al. 2003],

and the resulting feature distributions were modeled using Gaussian mixture models. [START_REF] Langs | [END_REF]] describe the object appearance using steerable features [START_REF] Freeman | [END_REF].

[van Ginneken et al. 2002] suggest to use locally orderless image features [Koenderink & Van Doorn 1999]. In this method, a given image is fed to a bank of multiscale Gaussian derivative lters. Then, rst statistical moments are extracted from local histograms in the ltered images (the histograms extents are related to the corresponding scale).

These moments represent the considered features. [START_REF] Van Ginneken | [END_REF] used in practice all derivatives up to the second-order, ve scales for Gaussian ltering, and the two rst moments, and applied ASM to the segmentation of chest radiographs. An optimal set of features per landmark is then extracted using a forward-backward feature selection scheme. reformulate the segmentation task as a maximum likelihood (ML) problem, that is optimized using particle ltering [START_REF] Isard | [END_REF]]. Iteratively, shape samples are drawn from the PDM. Then, particle weights are computed using a likelihood term (a kNN classier is used for pixel probability estimation akin to [START_REF] De Bruijne | [END_REF]). According to these weights, a random sampling of new particles is performed, followed by a small random perturbation of duplicated shapes. This process results in the multiplication of likely particles, whereas the bad shapes disappear gradually. Asymptotically, the algorithm converges to the maximum likelihood estimate. However, the optimal shape is approximated in practice by the strongest local mode of the particle distribution (that can be identied using the mean shift algorithm [START_REF] Comaniciu | [END_REF]) after a sucient number of iterations. The use of a large number of initial samples bequeaths to the algorithm additional robustness against local minima and more independence with respect to the initial conditions, with the downside of increasing the computational cost. The authors extend their approach in [de Bruijne & Nielsen 2005] to the segmentation of multiple shape instances in the same image, by dening interaction constraints between neighboring shapes, and adjusting the data term. [START_REF] Tu | Face Localization via Hierarchical CONDENSATION with Fisher Boosting Feature Selection[END_REF]] and [Qu et al. 2008] are other examples of successful combination of PDM with particle ltering approaches.

An alternative research direction deals with handling outlier candidates during the ASM search. In fact, the standard ASM search algorithm considers a least-square criterion with respect to the residual between the model and the data. The minimization of such criteria leads to optimal results only when the Gaussian noise assumption holds, which is not always true in practice. [START_REF] Rogers | [END_REF] tackle this problem by considering two options: M-estimators and random sampling techniques. M-estimators aim to alleviate the eect of outliers. This leads to a weighting strategy of the landmarks that depends on the size of their residual. The random sampling approach selects repeatedly small random subsets of landmarks, and uses them to estimate the shape parameters. Random sample consensus (RANSAC) [Fischler & Bolles 1981] is an algorithm that follows this rationale. It evaluates the obtained parameter estimation by the amount of its consistency with respect to the data, and returns as nal estimation the parameters with the largest consensus set. [Abi [START_REF] Nahed | [END_REF]] handle the outliers issue in a dierent way. They extract from the image to be segmented a set of candidate points for each landmark. Then, they use a robust point matching algorithm [Chui & Rangarajan 2000] to nd the best correspondences between a legal shape instance and the pool of candidates, while allowing outlier rejection.

Incorporating MRF regularization in the search step is another method to tackle the independent estimation of the landmarks optimal displacements in the standard ASM. [Behiels et al. 1999] introduce a regularization term that constrains the relative displacements of consecutive landmarks on the object curve. The optimal landmarks are the one determined by minimizing a cost function that combines the aforementioned pairwise regularization terms, and individual landmark terms. A minimum cost path search algorithm (using dynamic programming, thanks to the tree structure of the open 2D contour) solves the problem. The method exhibits improvements compared to ASM, as far as outliers handling is concerned. More recently, [START_REF] Tresadern | [END_REF] use a hybrid shape model that combines PDM and MRF. The unary MRF potentials correspond to a data term that describes the goodness of the candidate features, and the pairwise potentials account for the relationship between the relative displacements of pairs of landmarks. This regularization term is given by a Mahalanobis distance, where the mean and covariance parameters are learned from the data. In both methods, the solutions are then regularized using the global model. Alternative developments of ASM are related to the shape modeling part and question the Gaussianity assumption (section 2.1.1)

or seek a sparser representation (section 2.1.1).

Active Appearance Models AAM [Cootes et al. 2001] can be thought of as ASM with region-based image features. Besides modeling the shape variations, AAM encode the texture of the object and its variations. Hence, similarly to PDM that can synthesize new shapes, AAM are able to generate new images of the learned object. However, this generative model combines shape and appearance description into one linear system, such that the parameter vector controls both shape and appearance variations.

The rst step consists in warping the training shape images to a common mean frame, using a triangulation algorithm. The intensity values are then stacked into texture vectors, which are used to perform PCA after intensity normalization. This leads to expressing the shape appearance as g = g + E g b g , where g is the mean normalized texture vector, E g is the eigenbasis and b g is the vector of texture parameters. The shape variations are expressed akin to the ASM setting as s = s + E s b s . Then, the shape and texture (independent) parameters are combined in a common parameter vector b = W s b s b g , where W s is a diagonal matrix that weights each shape parameter, accounting for the units dierences between the shape and texture models.

Such a vector is generated for each training example, and then PCA is applied on the obtained pool of vectors. This yields the appearance model b = E c c, where the appearance parameter vector c controls both shape and texture variations: The key approximation of AAM is that the matrix R is assumed to be xed.

s = s + E s W -1 s E cs c g = g + E g E cg c ,
This matrix is learned from the training set using numeric dierentiation, and hence provides a priori knowledge on the model parameters adjustment during the search. This approximation leads to a tractable algorithm that solves an initially-tough high-dimensional optimization problem.

Dierent modications to the standard AAM exist in the literature, we briey report a few of them. [Matthews & Baker 2004] introduce a new computationally-ecient analytical gradient descent algorithm to optimize the AAM search and demonstrate its applicability for face tracking.

[ [START_REF] Andreopoulos | [END_REF] extend this latter approach, initially presented in 2D, to the 3D case and apply it to the segmentation of cardiac MRI images.

[ Donner et al. 2006] use the canonical correlation analysis algorithm (CCA) in order to capture the correlation between the parameter update and the image residual. This alternative method to compute the matrix R yields more precise parameter prediction than the standard numeric dierentiation method and allows speeding up the search by a factor of four. [Gross et al. 2006] build on their previous work in [Matthews & Baker 2004] to develop a robust AAM with respect to occlusions.

Graph-based Methods and Related Work [Boykov & Jolly 2001] introduce an interactive method of n-dimensional image segmentation based on graph cut. The user provides hard segmentation constraints by specifying some pixels as belonging to the object or the background. Then, the segmentation problem is cast into a min-cut/max-ow problem, where the internal nodes of the graph represent the image pixels, and the source and sink node refer to the object and background classes.

The internal graph links represent the neighboring system in the image, and their weights dene a discontinuity penalty that is related to the pixel intensity dierences (a boundary-related term). The terminal weights represent the penalty of assigning a pixel to one class or the other. These penalties represent the likelihood of the pixel intensity given the intensity distribution of the labeling class. The hard constraints are enforced by pairs of (0, K) weights, where K is a large penalty. The segmentation of the image is given by the min-cut.

This work is further reviewed in [Boykov & Funka-Lea 2006]. The authors survey the developments and variations of their method and discuss the historical steps that led to a ubiquity of the graph cut approach. Although the multilabel case is NP-hard, an extension using multi-way graph cut can be done in practice using α-expansion for instance. Dierent improvements were added to the standard graph cut formulation [Boykov & Jolly 2001]. For instance, [Blake et al. 2004] propose to learn the regional parameters from color images using a Gaussian Mixture Markov Random Field (GMMRF) generative model.

This idea is combined in the "Grabcut" algorithm [START_REF] Rother | [END_REF]] with an iterative minimization, where from an initial loose initialization provided by the user, the parameters of the model are updated and then the segmentation is estimated using graph cut, until convergence. This algorithm alleviates the user interactions while providing good performance. [Li et al. 2004b] use watershed pre-segmentation to increase the eciency of the algorithm. Such pre-segmentation provides super-pixels that are used to build the graph.

In a dierent line of work, [Li et al. 2006] propose a method to segment multiple 3D "terrain-like" surfaces. Their algorithm can combine edge-based costs and region based weights. They dene several geometric constraints to enforce surface smoothness and inter-dependencies between the dierent objects. They use an image unfolding operation to handle tubular structures. The optimization is performed by an s/t min-cut in an appropriately constructed graph to take into account both image costs and geometric constraints. They applied their method to 3D medical image segmentation of the diaphragm as well as the inner and outer walls of vascular structures.

[ Zabih & Kolmogorov 2004] propose a pixelwise segmentation method based on clustering. Their framework is general and applies clustering to image pixels according to their feature vectors, while imposing spatial coherence.

Their algorithm is applicable to parametric as well as non-parametric clustering methods, as long as the cluster quality function meets a linearity criterion with respect to the image pixels. The algorithm is twofold, similarly to the expectation maximization algorithm (EM). Iteratively, the clusters parameters are xed, and the best corresponding pixel labeling is found by solving a graph cut using α-expansion. Then, given the pixel labels, new cluster parameters are estimated.

[ Shi & Malik 2000] point out that the graph cut segmentation method is biased towards small region cuts that minimize the length of the partition boundary. They tackle this problem using a graph spectral approach, by dening a new grouping criterion. This introduced dissociation measure is called normalized cut. They obtain an approximate solution to the problem by dening a relaxed formulation, expressed as a generalized eigenvalue system. [Grady 2006] introduced a new graph-theoretical algorithm for multilabel interactive image segmentation, known as the random walker algorithm.

The method requires regions seeds dened by the user. The algorithm can be intuitively thought of as a random walker that starts from each pixel and arrives to each seed with a given probability. Then the largest probability values determine a labeling (a segmentation) of the image. The random walk is biased by the image structure, which is encoded by the weights of the edges connecting the pixels. [Grady 2006] provides a tractable means of solving the random walker problem by formulating it as a Dirichlet problem. By the use of combinatorial operators, the problem boils down to solving a system of linear equations. However, the size of this system depends on the size of the corresponding image, which can be very large in practice, in the 3D medical imaging setting for instance. Such an observation prohibits the use of direct linear system solvers, and motivates the use of iterative solvers. The random walks algorithm provides a global solution to the problem, and unlike graph cut, does not suer from the metricity bias or the "shrinking bias".

However, for the seeded image segmentation problem, this algorithm is biased towards the spatial locations of the seeds. [Grady 2005] extends the previous work to alleviate the need for user initialization by using prior intensity prob-ability densities. This modication allows the segmentation of disconnected objects belonging to the same class. A further improvement is presented in [Singaraju et al. 2008] where directed edges are used to construct the graph, resulting in directional asymmetric weights. This new formulation allows the improvement of the seeded segmentation results, by balancing the bias towards seeds locations.

Although the random walker formulation seems to be dierent from the graph cut formalism, it turns out that these two approaches are closely related.

[ [START_REF] Sinop | [END_REF]] present a unifying framework for seeded segmentation algorithms, the p-norm of the spatial gradient of a potential function is minimized, subject to Dirichlet boundary conditions, specied by the seeds. They

show that using the 1 norm is equivalent to graph cut, whereas using the 2 norm leads to the random walker algorithm. [START_REF] Singaraju | [END_REF] build on this work and interpret the potential function as an MRF with continuous valued variables. The MAP estimate corresponds to the minimization of the the dened general cost function in [START_REF] Sinop | [END_REF]. They tackle this problem by using interactive least square techniques. Moreover, by choosing an appropriate value of p (between 1 and 2), a trade-o between the graph cut behavior and the random walker behavior can be achieved. They show in the experiments that with p = 1.5, drawbacks of both methods are minimized.

As we have already noted, one drawback of graph cut methods for image segmentation is related to their bias towards shorter boundaries, which is known as the "shrinking bias". In particular, common techniques miss thin elongated structures during segmentation. This observation motivated the work of [START_REF] Vicente | [END_REF], who tackled this problem by adding constraints to the energy minimization problem. These constraints are intuitively translated into an interactive segmentation framework. Starting from an imperfect s/t segmentation, where these structures are missing, the user would click on pixels that must be connected to the object (through thin components).

The ensuing minimization problem turns out to be NP-hard. However, the authors develop a practical optimization algorithm called DijkstraGC, that combines aspects from the Dijkstra algorithm and graph cuts. The proposed approach performs well in practice, and the resulting interactive segmentation framework can handle thin image structures. An alternative (and slower) minimization strategy based on dual decomposition (DD) was considered as it provides a lower bound on the optimal energy in general, and the optimal solution in a special case. This property was used to verify that DijkstraGC also achieved an optimal performance in these cases.

An alternative way to improve the graph cut technique involves incorporating shape priors. [Funka-Lea et al. 2006] apply the binary graph-cut technique to the segmentation of the entire heart in computed tomography (CT)

images. This entire heart segmentation serves visualization purposes as it allows a 3D reconstruction that exhibits the coronary arteries. An automatic seed detection sets the hard unary potential constraints. A Potts interaction model is used, and an additional pairwise constraint is introduced. This socalled "blob" energy term encourages graph cuts that produce a shape which edges are orthogonal to the rays coming from the seed region. A dilation of the graph-cut output image mask gives the nal result.

[ Veksler 2008] adds a new pairwise regularization term to the standard s/t graph cut image segmentation formulation. This term corresponds to a generic shape prior that is related to a large class of shapes, namely the starshaped objects. The prior term assumes knowing a point inside the object that is a center of the star shape. This seed point is often provided by the user in interactive image segmentation tasks. Subsequently, she shows that the dened prior can be used as a "ballooning" force that biases the segmentation towards longer partition boundaries. This can be achieved by an appropriate parameter choice and allows the author to counter the "shrinking bias".

[ Das et al. 2009] adopt a similar approach aiming to incorporate a shape prior in the graph cut segmentation. They dene a class of shapes called "compact" shapes. These shapes satisfy some regularity properties on a 4neighborhood grid with respect to an interior point, which is assumed in practice to be close to the center of the object and specied by the user. Then, the standard boundary-related pairwise potential is modied such that shapes that violate the "compactness" condition are highly penalized. Moreover, a "bias" parameter is introduced similarly to [Veksler 2008]. By varying its values, one can encourage the cut towards longer object boundaries, and hence compensate for the "shrinking bias".

The work of [START_REF] Zeng | [END_REF]] is in the line of continuous eort to incor-porate global priors into graph cut segmentation. Their approach modies the max-ow algorithm to enforce topological consistency between an initialization and the nal segmentation. Their work is based on the concepts of topology of digital images and simple points [Bertrand 1994]. They prove that the dened topology cut problem is NP-hard. However, they propose an s/t graph cut is performed such that the regional term is related to the computed mean intensity values, and the shape prior term is enforced by the binary elliptic mask, using the t-weights. Then, a new ellipse is tted to the result of the graph cut, until convergence. Similar ideas are also used in the work of [Lin et al. 2005].

An alternative attempt to enforce prior knowledge in the graph cut interactive segmentation setting is presented in [START_REF] Freedman | [END_REF]. The authors use a xed template shape prior, represented by an unsigned distance map Φ. Then, for every neighboring pixels p and q in the image, an additional n-link weight is dened as Φ( p+q 2 ). This penalty encourages the cut to take place on the zero level of the distance map Φ, i.e. the learned contour. To deal with the pose estimation problem, [START_REF] Freedman | [END_REF] rely on the user input to compute a rigid transformation that maps the distance map to the image, using Procrustes Analysis. However, the user input is very often insucient, which impairs the quality of the rigid alignment. Then, scale variations are handled in a brute force manner. The use of this template prior shows improvements in the result of the graph cut-based segmentation. [Zhu-Jacquot & Zabih 2007] and [START_REF] Zabih | [END_REF]] introduce a method that combines statistical or parametric shape priors with graph cut in an EM framework, and apply it to the segmentation of the LV in 2D MR image.

In [Zhu-Jacquot & Zabih 2007], they apply PCA to the SDMs of training shapes to learn a statistical shape prior, whereas in [START_REF] Zabih | [END_REF], they use a geometric parametrized prior (concentric circles). In both cases, the shape prior is represented by a set of parameters, and the image intensity is modeled using a mixture of Gaussians. In the expectation step, the shape parameters as well as the photometric model are updated using a gradient descent, and in the maximization step, the pixels of the image are labeled according to the shape prior and the intensity model constraints using graph cut (or α-expansion if there are more than 2 classes).

The process is iterated until convergence.

Another approach that incorporates a statistical shape prior into an iterative graph cut segmentation process is proposed in [Malcolm et al. 2007a].

They use KPCA to learn a statistical shape model from binary images representing the training set. The algorithm repeats the following operations.

First, object and background histograms are computed from the current pixel labeling. Then, a pre-image of the projection in the feature space of the current segmentation is formed, using an approximate algebraic expression.

This pre-image is used to introduce a Bayesian shape prior in a unary regional image term. Then, a graph cut is performed and a new pixel labeling is produced. The algorithm iterates these operations, starting from an initial segmentation provided by the user, until convergence.

[ Ali et al. 2007] propose a method to segment 2D Dynamic Contrast Enhanced Magnetic Resonance Imaging (DCE-MRI) slices. Their approach combines graph cut and shape priors. They build from the aligned training slices a "shape image" or a template that is composed of three regions: the kidney, the background, and an area of variability, where the probability of belonging to one class decays exponentially as the distance to the considered class increases. Then, this probabilistic template is used during segmentation, by incorporating a shape prior term for each pixel, as a terminal weight. To this end, a preliminary alignment of the query image with respect to the training images is necessary, and is performed using a rigid transformation.

[ [START_REF] Vu | [END_REF] propose a segmentation algorithm based on a 2D shape prior within the binary graph cut framework. They use a template shape prior (a binary mask) and enforce the corresponding constraint using a shape distance derived from [Chan & Zhu 2005]. Essentially, in the discrete setting, this distance evaluates the region of disagreement between two binary shape images. This representation allows arbitrary topologies. Hence, the distance to the template of some labeling can be decomposed over the image pixels.

Thus, [START_REF] Vu | [END_REF] dene the shape prior term as terminal weights.

They also resort to a normalization technique to obtain invariance properties with respect to rigid motion in their model. Their framework allows the segmentation of multiple overlapping instances of the object, through the socalled multi-phase graph cut. This algorithm solves iteratively a graph cut for each object instance, while accounting for the overlap between objects through an appropriate choice of the data term (which is also encoded as terminal weights). Each iteration is exactly solved by the min-cut/max-ow algorithm, due to the submodularity of the regularizing pairwise terms. However, the multi-phase iterative graph cut algorithm is not guaranteed to converge to the global optimum. An extension of the approach allows the use of multiple shape templates to encode the prior knowledge. Successful segmentation results

under noise and occlusion are achieved.

Another template-based approach is introduced in [Schoenemann & Cremers 2007]. This method benets from global optimization based on graph theory. Their method seeks a 2D curve such that it is close to strong gradient locations and similar to a template curve.

They also seek to enforce invariance with respect to similarity transformations. They dene an energy that is the sum of three terms. The rst one integrates a decreasing function of the norm of the gradient along the curve.

The second term is translation-invariant, and measures how well the curve is aligned with the template. In the perfect case, the obtained mapping should be a dieomorphism. Then, a third term is introduced as a penalization of the mapping. The problem is then cast as nding the best curve and the best mapping that minimize the dened energy functional.

All these terms are actually normalized by the length of the sought curve to eliminate bias toward shorter boundaries. This ratio of integral along the curve paves the way to the use of the minimum ratio cycles (MRC) minimization technique [Jermyn & Ishikawa 2001].

Then, inspired by the work in [Ishikawa & Jermyn 2001], scale invariance is enforced through the use of a product graph. In this appropriately built graph, the authors solve for a MRC optimally and eciently, by parallelizing computations. Their model is not however rotation-invariant. This is addressed by matching to dierent rotated instances of the template. Using this s/t min-cut formulation, they show improvements over similar variational approaches. More specically, they dene geometric and photometric constraints using distribution matching, rather than pixelwise costs. They dene the optimal segmentation (or binary labeling) in an unseen image frame as the one that maximizes the Bhattacharyya coecient of the color and shape object distributions with respect to the one of the manually segmented frame.

The image intensity distribution is dened using a kernel density estimate.

The shape of the object is described by the collection of normalized distances separating all its interior points and a centroid point. A distance distribution is also dened using a kernel density estimation. In addition to these distribution matching terms, a smoothness constraint that is related to the length of the partition boundary is considered. The obtained minimization criterion is then relaxed using a rst-order approximation of the Bhattacharyya coecient. The resulting energy is dened as the sum of unary and pairwise potentials, that are suitable to s/t graph cut minimization. Moreover, a global optimum is guaranteed as the pairwise penalties are submodular.

Discussion & Related Work The aforementioned graph-based segmentation methods that made use of a (possibly statistical) shape prior essentially model the image pixels using a MRF. Hence, the MRF variables account for pixel labels. This representation narrows the possibil-ities of incorporating a shape model while being computationally expensive.

With the exception of very few techniques ([Ben Ayed et al. 2009], [Schoenemann & Cremers 2007]), the overwhelming majority of the approaches dene the shape prior at the pixel level, through the t-links weights.

To this end, an implicit shape representation through a distance map [START_REF] Freedman | [END_REF], or the use of a binary shape image [START_REF] Vu | [END_REF] is required. A single template representation is insufcient to capture large shape variations, and the use of numerous templates to enforce variability might be limited compared to statistical approaches, and comes with a large computational burden. This template representation applies also to the statistical approaches, where a probabilistic shape image [Ali et al. 2007], or a pre-image of the projection of the segmentation on learned subspace [Malcolm et al. 2007a] were used. We also notice that in many cases, the coupling between the statistical shape parameters and the MRF variable is complex, and is therefore accounted for through the use of iterative optimization schemes, leading to suboptimal solutions. We note also that the need of aligning the template to the image is an ill-posed problem because the pose of the object to be segmented is a priori unknown. In addition, the feasibility of extending these methods to the 3D case is not always guaranteed. Nevertheless, the exibility of the graph cut framework allows these dierent methods to combine several image cues, from regional features to boundary evidence. A dierent modeling approach uses the MRF to represent parts of the object, or boundary landmarks. These methods are more closely related to our work.

[ Zhang et al. 2004] introduce a technique to segment human bodies using a graphical model. They represent the articulated shape using a triangulated silhouette, with a specied vertex elimination order. They augment this representation by a parametrization that species articulation joint angles. The graphical model encodes the dependencies among the nine considered angles, whereas each silhouette landmark is considered to be dependent on its parent triangle. The root triangle corresponds to the face part. After learning the parameters from a large database of training data, the obtained translationinvariant model is combined with several visual cues (edge gradient map, background subtraction, skin color, etc.) to perform human body segmenta-tion. The inference is realized using particle ltering. Although the presented results are appealing, this method suers from its application-specic design.

Hence, it is not generalizable in a straightforward manner to other applications, or to higher dimensions. [Felzenszwalb 2005] adopts a dierent approach, that is related to the spirit of the pictorial structures work [Felzenszwalb & Huttenlocher 2005]. A triangulated polygons representation of 2D shapes is proposed, its dual structure is a tree (the model is in fact a 2-tree, which is a chordal graph). Such a template is used to segment an unseen image, using boundary information. Beneting from the perfect elimination scheme property of their graph, they are able to provide a global solution to the segmentation problem. The shape prior is enforced in each triangle clique through a constraint on the local induced transformation, which is required to be as close to a similarity as possible. The resulting segmentation algorithm is similarity-invariant and handles naturally occlusions.

[ Seghers et al. 2007a] propose a segmentation method that shares ideas with the approach in [Felzenszwalb 2005]. A 2D contour is represented by a set of landmarks structured in a chain.

For each landmark, a classier is trained based on a locally orderless image feature vector, following [START_REF] Van Ginneken | [END_REF]. Segmentation is then cast into the problem of selecting the best candidate points found by the classier, while satisfying the statistical shape prior constraints. The latter are dened by the Gaussian distributions of the relative positions of two successive landmarks on the chain.

A higher order prior that conditions this relative position with respect to the one of the previous landmark pair on the chain was also used. Inference was done using a dynamic programming approach, namely the minimal cost path (MCP) between the "root" and the "leaf" of the chain. Successful segmentation results of the lungs in standard posterior-anterior chest radiographs were shown. A more complex model, like for instance a closed contour, or a connected graph, makes the use of dynamic programming impossible. Moreover, this exact inference is not extendable as it is to the 3D case (without making approximations), to handle minimal cost surfaces ( [Grady 2010] has recently surveyed the work of [Sullivan 1990] and proposes a method that tackles the problem). Therefore, [Seghers et al. 2007b] resort to approximate optimiza-tion techniques to extend their work to the 3D case, and showed results of liver segmentation in CT images. A heuristic approach called iterative dynamic programming is used, where the problem is solved on randomly extracted paths from the meshed surface using DP, and the solutions are combined using a voting scheme. In [START_REF] Seghers | [END_REF]] mean eld annealing (MFA) was also considered for optimization.

Another classier-based discrete detection method is proposed in [START_REF] Donner | [END_REF]. Their approach is based on symmetry interest points and descriptors, that are extracted from gradient vector ow (GVF) elds [Xu & Prince 1998]. Given a binary edge map, the gradient vector ow is nearly equal to its gradient when the latter is large, and varies smoothly in homogeneous regions. GVF minima locations are of a great interest as they correspond to local symmetry points. [START_REF] Donner | [END_REF]] use this property to dene their MRF model. The MRF nodes correspond to landmarks at local symmetry locations, and the edges encode inter-landmark dependencies. Local descriptors using local GVF around the interest points are used to dene the unary potentials, as well as the pairwise shape prior. The latter uses relative edge / descriptor orientation as well as the length of the edges. Therefore, the shape prior is invariant with respect to rigid motion. During inference, candidate points are extracted from the unseen image using its GVF. The model is matched to the image as a MAP-MRF problem that is solved using the MAX-SUM diusion algorithm introduced by [START_REF] Schlesinger | [END_REF]] and reviewed in [Werner 2007]. They apply their method to the segmentation of 2D hand X-rays and 2D spine MRIs [Donner et al. 2007b]. To make the training phase easier and lessen the manual annotation needs, they propose in [START_REF] Donner | [END_REF]] a technique to learn the model in a weakly-supervised manner. Note however that their approach does not handle the object boundaries. Thus, it can be considered as techniques to be used more for detection than for segmentation.

[ Kumar et al. 2010] present a knowledge-based segmentation approach that bridges the gap between graph-cut methods and pictorial structures. On one hand, the image is modeled by a MRF, and a binary labeling l that species the object and the background is sought. On the other hand, they use a layered pictorial structure (LPS) model to represent articulated objects, where the parts are connected in a complete graph (for a better spatial consistency).

The layer term refers to a depth number associated with each part, in order to allow overlap between the dierent components. A non-articulated object was modeled as a single part. Each part was specied by a set of exemplars, forming the so-called set of exemplars model (SOE), as well as a pose parameters vector. The authors provide a Bayesian formulation to the problem using a graphical model representation, that adds essentially a shape variable S (referring to the LPS) on top of the MRF representing the image. In addition to the standard s/t graph cut regional and boundary terms, the dened energy features a shape prior term, that is computed using a SDM of the prior shape instance. Inference in a new image I is performed by maximizing the posterior p(l|I), instead of the joint probability p(l, S|I). This is done by marginalizing over the shape prior. This marginalization is performed approximately and eciently through a sampling strategy. Indeed, the authors suggest to use the generalized expectation maximization algorithm (EM). In the E-step, BP is applied to estimate the probability of the shape knowing the current image labeling, in order to perform sampling. During the M-step, maximizing the posterior turns out to be eciently approximated by solving only one s/t graph cut, where the shape prior cost corresponds to a linear combination of the costs induced by the sampled shape instances in the E-step, weighted by their respective likelihoods. The so-termed OBJCUT algorithm converges to a local minimum.

The aforementioned line of research considers top-down modeling using MRFs or graphical models. While [Felzenszwalb 2005] succeeds in exactly optimizing the energy functional he denes, the method does not generalize to higher dimensions. More complex modeling that would overcome this limitation as in [Kumar et al. 2010] makes the use of approximations inevitable.

Moreover, the complexity of the method would highly increase in the 3D case since it is related to the number of image pixels, and the search space for the pose estimation will be substantially bigger. The approaches proposed in [Seghers et al. 2007a] and [START_REF] Donner | [END_REF]] rely exclusively on edge-based cues or local features. In some applications, the use of regional cues can be of a great interest.

In the following sections, we introduce a nearly fully-automatic segmenta-tion framework that permits to plug-in various image cues, as in the standard graph cut setting [Boykov & Jolly 2001]. The use of regional image information is introduced in section 3.2, and classier-based segmentation is explored in section 3.3, in a relatively similar manner to the work of [Seghers et al. 2007a] and [START_REF] Donner | [END_REF]. Our formulation does not make any assumption about the dimensionality of the image. It is hence extendable to the 3D case, at a relatively low cost. The complexity of our method is related to the number of landmarks, and as we do not explicitly model the pose parameters, the increase in the dimensionality of the search space does not aect the tractability of the problem. Moreover, for the special case of regional image features, we benet from the Green's theorem to substantially speed up computations. As far as optimization is concerned, we resort to approximate inference and our method leads to a local minimum in general. However, the experiments show that in practice, satisfactory results are obtained, and that the segmentation algorithm is robust to noise and occlusions thanks to the shape prior. In a particular case where the dimension of the problem is small, we show that exact inference is tractable using dynamic programming.

We also mention that our work is related to the registration algorithm of [START_REF] Glocker | [END_REF] (initially introduced in [Glocker et al. 2007]) and to the one of [Tang & Chung 2007]. Both methods are based on MRFs where the labeling corresponds to quantized displacements of control points. This idea will be used in the following.

Region-based Segmentation

The main challenges of knowledge-based segmentation are: (i) appropriate modeling of shape variations, (ii) successful inference between the image and the shape model. Let us consider the simplest possible scenario that aims to detect an object of particular interest from the background in an image I. Then, the prior-based segmentation task corresponds to the estimation of the shape model parameters from the observation. Formally, ones seeks to maximize the posterior probability p(G k |I), where G k = (V, E k ) is the k-fan statistical shape model introduced in section 2.2. By a simple manipulation using the Bayes rule, the segmentation is expressed as the MAP estimation of the MRF parameters, or:

max p (G k ) p (I|G k ) = max p (x 1 , . . . , x n ) p (I|x 1 , . . . , x n ) .
(3.12)

By taking the negative log of the latter equation, we formulate the segmentation problem as an energy minimization problem. We seek optimal landmark positions (x * 1 , . . . , x * n ) such that:

(x * 1 , . . . , x * n ) = argmin (x 1 ,...,xn)∈Ω n E (x 1 , . . . , x n , I) = argmin (x 1 ,...,xn)∈Ω n E shape (x 1 , . . . , x n ) + E image (x 1 , . . . , x n , I) . (3.13)
The energy E (x 1 , . . . , x n , I) corresponding to the landmark positions (x 1 , . . . , x n )1 in the image domain Ω is the sum of a data-related term E image (x 1 , . . . , x n , I) expressing the local image cost conditioned on the control points, and a prior term E shape (x 1 , . . . , x n ) expressing the cost of deforming the landmark pairs in E k with respect to the prior learned distributions.

We explain in this section how we dene these two energy terms, and we detail in particular the relationships between the control points and the image domain. We also develop an optimization procedure that enables to solve (3.13) using an ecient discrete optimization algorithm.

Regional Statistics and Image Segmentation

By applying the data-related cost, we seek the optimal separation of the object from the background in terms of visual properties. Let p O and p B be the prior densities for these two hypotheses. We suppose that the control points form a closed boundary, or that there is a tractable way of inferring such a boundary from the control points positions, using an interpolation scheme for instance. Then, they partition the image domain Ω into an object domain O and a background domain B. Note that we consider this simplifying hypothesis for the sake of a simple presentation and lighter notations. In the general case, the image domain Ω is partitioned into several O r objects and a background domain B. The following derivation holds in this case too.

Upon considering the negative log of the posterior probabilities, we express the cost E image (x 1 , . . . , x n , I) using the regional statistics [START_REF] Zhu | [END_REF] as follows:

E image (x 1 , . . . , x n , I) = O -log (p O (I y )) dy + B -log (p B (I y )) dy = O -log p O (I y ) p B (I y ) dy + O -log (p B (I y )) dy + B -log (p B (I y )) dy = O -log p O (I y ) p B (I y ) dy + Ω -log (p B (I y )) dy constant . (3.14)
In order to evaluate this component and associate it with the proposed shape representation, we decompose the image domain Ω according to the control points (x 1 , . . . , x n ) by considering their Voronoi diagram (a survey on Voronoi diagrams is available in [Aurenhammer 1991]). The Voronoi tessellation partitions Ω into cells Ω i that are associated with the control points x i . Every cell Ω i consists of all image pixels closer to x i than to any other landmark. Fig. By intersecting these Voronoi cells with the "object" domain, we obtain the partition O = ∪ n i=1 O i that relates each pixel of the object part of the image to one control point, with O i = O ∩ Ω i . Note that thanks to the rearrangement of equation (3.14), we only need to consider the "object" domain pixels, and this is equivalent to using all the image information. Thanks to the Voronoi tessellation, one can decompose the global image term (3.14) into sub-terms which are dened at the partition cells as follows, up to an additive constant:

E image (x 1 , . . . , x n , I) = n i=1 V i (x i , I) (3.15) with V i (x i , I) = O i -log p O (I y ) p B (I y ) dy (3.16)
being the image-related cost associated with the control point position x i .

Clearly, the notation in equation (3.16) is abusive, because the data term V i (x i , I) depends on the positions of all the control points, through the use of the Voronoi cells. The independence hypothesis will be considered under some mild assumptions that will be detailed later. The above dened terms can be calculated very eciently in practice by combining rasterization techniques and fast integral computing over polygons [START_REF] Unal | [END_REF]]. This observation stems from Green's theorem, that permits to compute the integral in (3.16) using the boundaries ∂O i of O i , or:

V i (x i , I) = ∂O i I ω dω , (3.17)
where I is an integral of the image -log p O (I) p B (I) . Green's theorem had also been applied in a similar context in earlier segmentation algorithms, as in [Chakraborty et al. 1996] for example. We should note that this term is equivalent to using the entire image domain to determine the image support and can be replaced either using more complex descriptors (this will be explored in section 3.3), or through edge-driven support. Note that equation (3.16) generalizes to the case of multiple object domains O r in a straightforward manner:

V i (x i , I) = r O r i -log p O r (I y ) p B (I y ) dy , (3.18)
where O r i = O r ∩ O i , and p O r is the prior distribution of the color values in O r . In practice, we used simple Gaussian models in some applications and mixture of Gaussians models in others to model the object and the background intensity distributions. Such a component will perform well if the data support is strong but will fail to deal with noise, clutter, missing parts, etc. The use of prior knowledge on the expected geometry of the shape could address the above mentioned limitations.

Prior Knowledge and Image Segmentation

In the context of our approach, we have dened the shape model as a k-fan graph. Furthermore, we were able to determine an approximate density of this model using a small number of joint probabilities. In order to impose the prior, we minimize the cost E shape (x 1 , . . . , x n ) that we decompose over all the pairs in E k , as derived in equation (2.30):

E shape (x 1 , . . . , x n ) = (i,j)∈E k V ij (x i , x j ) , (3.19) with V ij (x i , x j ) = -γ ij log p x i -x j d . (3.20)
Note that equation (3.20) carries a notation abuse as well. Indeed, the shape scale d = 2 n(n-1) i<j x i -x j depends on the positions of all the control points. This issue will be accounted for in the following development. This model allows for the encoding of global dependencies as local combinations of individual pairwise densities.

One can now combine the data term with the prior term towards knowledge-based segmentation, by combining (3.13), (3.15) and (3.19):

E (x 1 , . . . , x n , I) = α n i=1 V i (x i , I) + (i,j)∈E k V ij (x i , x j ) . (3.21)
The factor α that we introduced in the above equation weights the relative importance of the data term with respect to the prior term. Note however that the relative weightings of the base k-clique landmark pairs with respect to the "non-reference" pairs are dened by γ ij (in equation (2.29)).

Energy Minimization

The optimization of the cost function (3.21) in the continuous domain is rather problematic. One can expect that it is not convex and therefore a gradientbased optimization is prone to failure. In order to optimize such an energy functional, we consider recent results from discrete optimization, and adopt a quantization approach similar to those presented in [Glocker et al. 2007] and in [Tang & Chung 2007].

We make two assumptions that are most often veried in practice. First, the initial positions of the control points are within the image domain, and that is why we can assume an upper bound on the maximum displacements that would lead to the solution. Second, we consider that the precision required about the solution is specied, which enables to choose a quantization step.

Let us suppose that we are given an initialization of the shape (x 0 1 , . . . , x 0 n ).

Then, a segmentation of the image refers to an optimal displacement of the control points, that minimizes the dened energy. Then, we can approximate the continuous deformations of our shape model towards the solution by a nite set of displacements vectors D = {d 1 , . . . , d z }. Let L = {1, . . . , z} be the set of labels associated to the quantization D of the displacements. Then, displacing the control point x 0 i by the vector d l i is equivalent to assigning the label l i ∈ L to x 0 i . The minimization of the energy in (3.21) can be hence written as a labeling problem, or:

(l * 1 , . . . , l * n ) = argmin l i ∈L E x 0 1 , . . . , x 0 n , I, (l 1 , . . . , l n ) , (3.22)
where the data unary terms and the prior pairwise terms in (3.21) are adjusted as follows:

V i (x i , I, l i ) = O i +d l i -log p O (I y ) p B (I y ) dy (3.23) V ij (x i , x j , l i , l j ) = -γ ij log p d ij (l i , l j ) d , (3.24)
where d ij (l i , l j ) = x i + d l i -x j -d l j . We illustrate in Intuitively, such an approximation of a higher order term generates small inaccuracies if the displacement range is relatively small. On the other hand, equation (3.24) uses an approximation of the higher order quantity d. Indeed, when estimating the costs of pairs of labels, the shape scale is supposed to be xed at a value estimated from the initial attributes (x 0 1 , . . . , x 0 n ). This assumption is also acceptable in the case of small displacements. As we shall see next, starting from an initial shape position, we iteratively solve sequentially several labeling problems based on the formulation in (3.22), where the solution of one step initializes the next iteration. Although the shape scale is xed within each iteration, it is updated from one step to the other, and the trade-o between the regional image information and the shape prior allows for scale changes. The resulting approach hence deals with translations, rotations and scaling without explicitly modeling them. Let us now detail the iterative search scheme that we propose.

The cardinality of the label set is quite important since on one hand it denes the accuracy of the search, while on the other hand increases the complexity of the algorithm. In order to address these issues, we rst consider an approach that is incremental in terms of displacements. To this end, we cope with the accuracy issue, that is closely related to the quantization of D, by using a pyramidal coarse-to-ne approach. Each level of the pyramid corresponds to a quantization step that is rened in the following level. To speed up the convergence in each level of the pyramid, and be consistent with the local displacements assumptions that we mentioned earlier, we also adopt an incremental approach in terms of the label set, where in each iteration t we look for the set of labels that will improve the current solution. This label set refers to a local search neighborhood. More explicitly, we repeatedly minimize:

E t (l 1 (t), . . . , l n (t)) = α n i=1 V i (x t-1 i , l i (t)) + (i,j)∈E k V ij (x t-1 i , x t-1 j , l i (t), l j (t)) , (3.25)
with

x t i = x 0 i + t τ =1 d l i (τ ) , (3.26)
until convergence, with l i (τ ) being the optimal label associated with the i th control point at iteration τ .

Recovering the optimal solution of this objective function is known to be an NP-hard problem and the complexity is dependent mostly on the pairwise potentials function, that is not submodular in our case. Hence, we consider an approximate solution to the labeling problem. We presented in section 3.1.1 several methods that permit to tackle the minimization problem, such as α-expansion [Boykov et al. 2001], FastPD, TRW-S [Kolmogorov 2006] and loopy belief propagation (LBP). We used in practice the computationally ecient Fast-PD algorithm [Komodakis et al. 2008b] Compute the energy terms: the unary potentials and the pairwise potentials.

12:

Run FastPD to solve (3.25).

13:

Update the control points according the the optimal labeling, t ← t + 1. 14:

until Convergence: the labeling stabilizes. 15: end for Implementation details The algorithm we propose uses a local search method. Hence in principle, the initialization is required to be quite close to the solution. However, we derived a method to deal with an initialization that is relatively far from the object of interest. We consider a subset of a few control points, that will be used to roughly locate the object by estimating similarity transformations (typically 3 points in 2D, 4 points in 3D).These control points are connected in a complete clique. We apply the same procedure described in algorithm 2 using these landmarks, without restricting the search to a small neighborhood. However, to prevent the distortion of the shape, we project the transformation resulting from each inner loop on 2 An implementation is available online: http://www.csd.uoc.gr/~komod/FastPD/ the space of similarity transformations, in order to update the landmarks positions. Then, the shape boundary is interpolated using thin plate splines (TPS). After this initialization step, the model is close to the solution, and we apply algorithm 2 using our shape model. An example of this process is shown in Fig . We also note that our method does not impose any explicit constraints with respect to the local transformations that the shape model undergoes during the iterations. While this leads to a simple formulation without any need to quantize rotations for instance, the resulting formulation does not prevent control points displacements that would cause boundary self-intersection or surface ipping. However, thanks to the geometric properties of our shape model, such an undesired behavior should be unlikely, which is veried in practice. Nevertheless, due to our construction, our algorithm is not suitable to handle self-occlusions of an articulated object, as permitted in the layered pictorial structure [Kumar et al. 2010] for instance.

Another observation concerns the energy minimization with respect to the iteration scheme. Due to the approximation made at the level of the data term, the estimated energy within the FastPD algorithm will dier from that computed once an optimal labeling is found. Hence, we carefully controlled these two quantities in order to ensure that the energy decreases monotonically. Any violation of this condition leads to the termination of the algorithm.

Experimental Validation

After preliminary testing on synthetic data, we validated our method by segmenting hands in 2D images, and segmenting the left ventricle in 3D computed tomography (CT) images.

Synthetic data As a rst step to assess the performance of the proposed region-based segmentation method, we tested our algorithm on synthetic images that are corrupted by Gaussian and salt and pepper noise, and where the object is occluded. In these preliminary experiments, we considered the simple case of a white rectangle on a black background. We modeled this object using 8 landmarks, and considered a complete graph in this case. Gaussian probability distributions were used to model the pairwise distances, as well as the object and background intensities. Moreover, the scale of the object was supposed to be known in these experiments. Fig. 3.4 shows the result of a simple case. Although rotations are not modeled within our search algorithm, the shape model is able to recover the rotated object through local landmark displacements. A more dicult example in shown in Fig. 3.5 where the object is partially occluded. This is also handled by our method naturally, without a dedicated modeling. Fig. 3.6 simulates the tracking of a moving object that is occluded during its displacement. No dynamic prior has been used.

Hand segmentation In a subsequent experiment to validate the performance of our method, we considered the application of modeling the hand using a 2D 40-example dataset of annotated left hands, showing dierent relative nger positions, hand sizes, and textures [Stegmann & Gomez 2002]. On each hand contour, 56 landmarks were used to describe the structure. We per- formed clustering in the diusion maps space as described in section 2.2.4, using shape maps [START_REF] Langs | Modeling the Structure of Multivariate Manifolds: Shape Maps[END_REF]]. The clustering provided 11 clusters shown in Fig. 2.5(c). The constructed model was used as a shape constraint as shown in Fig. 2.6(a)-2.6(d), and we hereafter apply it in a segmentation setting. We considered a multi-scale implementation of the approach using a gradually increasing number of control points to accelerate convergence, and handle non-local initialization, as explained in the previous section. shown in Fig. 3.8(a). We also compared quantitatively our method to AAM segmentation Fig. 3.7 based on the dice overlap coecient with respect to the ground truth data. We can see in these boxplots that our algorithm performs better quantitatively with examples where the forearm is hidden by a sleeve. In the case of bare forearms, the data term drives the model to "oversegment" the hand in comparison with the ground truth, which explains our results. These "oversegmentations" are visually correct (especially the ngers are correctly segmented) as we can see in Fig. 3.8(a). The three examples where our method did not succeed are particularly dicult to handle for our algorithm because they exhibit occlusion between the ngers, which can cause folding in the evolving contour. The rst two segmentation examples shown in Fig. 3.8(b) demonstrate that our method can sometimes handle such difculties, but this is not guaranteed. Towards checking the robustness of the method with respect to occlusions, we removed some hands parts for several examples, and despite the important missing structure, the results were quite satisfactory as shown in Fig. 3.8(b). The weight of the prior in these cases was increased, and enforced the correct segmentation, as the data term was less reliable. To further assess the robustness of our method, we added severe Gaussian noise to the database images. The segmentations obtained in Fig.

3.8(a) are completely or almost recovered, thanks to the prior knowledge, as shown in Fig. 3.8(c). Eventually, we used our segmentation method in a real setting, on hand video frames, with a cluttered background and partial occlusion cases. Fig. 3.8(d) gives some examples of the obtained segmentations.

We could reproduce the result we obtained on the noisy images using AAM segmentation [Cootes et al. 2001], but this algorithm could not reproduce our results for the occlusion cases. [Gross et al. 2006] developed an AAM algorithm that addresses the occlusion problem. Note that this is a dedicated development of the standard AAM, whereas our method naturally performs well for the segmentation of a partially occluded object.

Left Ventricle segmentation In this experiment, we did not use our k-fan graph model. Instead, we considered the complete graph to enforce shape prior. In order to evaluate the performance of our method in the three- While no expert manual segmentation was available for this dataset, we could compare to the results given by the Random Walker algorithm described in [Grady 2006] 3 . First, to build the model, we selected randomly 11 CT images as a training set. We placed manually on the surface of the left ventricle of one example from the training set 90 control points, that we will call P 90 .

The remaining 10 examples of the training set were then registered to the labeled example using the method described in [START_REF] Taron | Registration with Uncertainties and Statistical Modeling of Shapes with Variable Metric Kernels[END_REF] and correspondences between the control points instances were consequently deduced.

After smoothing the segmentation mask obtained by [Grady 2006] on the labeled example, we generated a meshed surface S of the myocardium and the blood pool. We learned next the probability density distributions p ij of the normalized chord lengths (section 2. As a rst application, we tested the consistency of the learned prior by only minimizing the shape term of the designed energy, and canceling the data term, starting from randomly perturbed positions of the control points. The prior constraints made the initial shape converge in all cases to the "mean" shape. Next, the segmentation was performed following algorithm 2 by following two steps: (i) the surface S is deformed from its initial position by using the control points P 4 , and by applying a similarity transformation to the mesh after each iteration. This is repeated till convergence. (ii) The control points P 90 are introduced and the mesh is transformed using a Thin Plate Spline (TPS) [Bookstein 1989] deformation. For the segmentation experiments, the During implementation, we noticed that computing the Voronoi diagram in 3D, and performing the intersections of the resulting cells with the myocardial and blood pool surfaces was time consuming. Hence, in order to increase the computational eciency, we introduced an additional approximation. In addition to the control points and the meshed surface, we computed for each landmark x i the voronoi cells O myo i and O bp i oine. Then after each inner- iteration of the algorithm, these cells were updated similarly to the mesh using TPS. Thereby, complex geometric computations were no longer required.

We compared our results quantitatively for the whole dataset to [Grady 2006] and compiled them in Tab. 3.1.

Detector-based Segmentation

We presented in section 3.2 a region-based shape-driven segmentation algorithm and demonstrated its performance in a variety of settings. However, regional image cues are not always suitable in practice. An approach that uses such image features makes the implicit assumption that image regions are smooth or have homogeneous intensity distributions, at least locally. In applications where objects have a particularly textured appearance, this hypothesis does not hold, and region-based statistics would lead to a poor segmentation performance. For instance, Fig 3.11 shows posterior-anterior chest radiographs (PA X-Rays). This is an example of images that can be hardly segmented using region-based costs. Moreover, in addition to shape variations, object appearance can dramatically dier from one case to another. This is also the case in computer vision and in particular in video-based surveillance where segmenting and tracking people can be very dicult as human body appearance is almost arbitrary, depending on the clothes that are worn. As an alternative to image intensity distributions, one can consider feature vectors extracted from the image. Such high dimensional descriptors benet from the neighboring data around a particular image location, and can exhibit strong discriminative power. In our approach, we seek optimal landmark positioning.

Therefore, a discriminative image support that isolates highly probable control point locations in the image is of a great interest. Indeed, combined with a shape prior, feature-based detectors can yield very powerful segmentation tools. [Seghers et al. 2007a] pursue this line of thought as they use a chainstructured shape prior to nd the best landmark positions among candidates extracted from the image, based on locally orderless image (LOI) features.

In the same spirit, [START_REF] Donner | [END_REF] focus on symmetry interest points detected using GVF-based cues, and enforce a shape prior through a MRF formulation. Similarly, we propose to build image classiers that discriminate the control points from the background. These classiers provide landmark candidates, with associated condence values. The segmentation task boils down to nding the trade-o between shape prior penalties, and candidate selection costs, in an integrated manner. In practice, we adopt a similar approach to the LOI for feature extraction. Then, we apply the Adaboost algorithm to the computed local descriptors in order to train the landmark classiers. We observed in our experiments that in some dicult cases, the classiers may fail to provide good candidates to the segmentation algorithm, which leads to a suboptimal result. We propose to account for the missing correspondence problem by incorporating an additional label in the segmentation formulation to account for outliers. Although this idea has already been mentioned in [Donner et al. 2007a] and [START_REF] Seghers | [END_REF], the authors do not provide a clear method to post-process the outlier landmarks, and the impact of the procedure is not shown in the experiments. Our approach is detailed in the next section, and the results that we obtain suggest that accounting for the outliers in the segmentation process can lead to improved results qualitatively and quantitatively in the case of the complete graph.

The remainder of this section is organized as follows: rst, we describe the features that we use to build the control points descriptors and train the corresponding classiers. Then, we explain how we embed this information in our segmentation framework. In the next subsection, we propose a method to account for outliers during segmentation, and "repair" their positions in a post-processing step. The last part is dedicated to the segmentation results we obtained on PA chest X-ray images. We show results for both cases: the case without handling outliers, and the case where outliers are accounted for.

Landmark-based Image support

We will explain in the following how we learn a classier for each control point of the object, using features that are computed from the image. These trained classiers are then used as detectors in our algorithm, to segment unseen images. Our method consists of extracting from the image a set of candidates for each landmark, and then choosing the best ones in terms of detection response on one hand, and global delity to the learned shape prior on the other hand.

Boosted Landmark Detectors

We present our method by describing the features and the classiers that we use in practice. However, the framework is exible, and dierent types of detectors can be plugged in, using for instance the gradient vector ow (GVF) descriptors that were considered in [Donner et al. 2007a], or the locally orderless image (LOI) features that were adopted in [Seghers et al. 2007a]. The recent boundary detection framework introduced in [Kokkinos 2010] can also represent an alternative choice to dene the image features. The author considered the distribution of image gradient orientation in the neighborhood of the landmarks at multiple scales through Scale-Invariant Feature Transform (SIFT) descriptors [Lowe 2004], and then performed discriminative dimensionality reduction using the Sliced Average Variance Estimation algorithm (SAVE) [Cook & Lee 1999].

In our work, we rst process the image I using a lter bank of derivatives for dierent Gaussian smoothing values σ (σ ∈ {0.5, 1, 2, 4, 8}). The output of this process is a feature image F as showed in Fig. 3.12. Such a feature extraction was previously considered by [START_REF] Van Ginneken | [END_REF] and [Seghers et al. 2007a]. Then, they obtained locally orderless image (LOI)

descriptors [Koenderink & Van Doorn 1999]

j ) 1≤i≤N , with the class label y (i) j ∈ {-1, 1} referring to the background/object, we train the classiers f (i) T in order to discriminate the control points from their neighbor- hoods. Before explaining the relationship between these classiers and the data-related energy term, we describe briey the Boosting framework.

Boosting is a classication method that is based on machine learning. It uses labeled data to learn a decision function. Let our training data be a set of N pairs (ϕ i , y i ) ∈ Φ × {-1, 1}, 1 ≤ i ≤ N , where Φ is the feature vector space and y i are the binary class labels. Boosting aims to nd a highly accurate decision rule f T as a linear combination of many "weak" or basic hypotheses h t , or:

f T (ϕ) = T i=1 α t h t (ϕ) .
L (f T ) = N i=1 λ (f T (ϕ i ) , y i ) . (3.28)
There is an extensive literature on Boosting, and several reviews on the subject exist. Among others we refer to [Schapire 2002] and [Meir & Rätsch 2003].

Dierent Boosting algorithms had been proposed. Their dierences lie in the denition of the cost function, the selection of the weak classiers and the computation of the step size. We restrict our description to the algorithm that we used in the experiments.

Adaboost A well-known Boosting algorithm is Adaboost which was introduced by [Freund & Schapire 1996]. This machine learning algorithm was then generalized by [Schapire & Singer 1998]. In this case an exponential loss function is used, or λ (f T (ϕ i ) , y i ) = exp (-y i f T (ϕ i )). Essentially, at each Boosting step a weak classier that best discriminates "positives" and "nega- Training: Learn weak classier h t : X → R using distribution D t .

5:

Compute the weight: α t ∈ R.

6:

Update:

D t+1 (i) = D t (i) exp (-α t y i h t (ϕ i )) Z t ,
where Z t is a normalization factor.

7: end for 8: Output: The strong classier:

H (ϕ) = sign (f T (ϕ)) , f T (ϕ) = T t=1 α t h t (ϕ) .
Let us now summarize the learning phase. For each landmark x i , "posi-

tive" examples and "negative" examples are extracted from the training images to form a learning set (ϕ

(i) j , y (i)
j ) 1≤j≤N . The latter is used to train an Adaboost classier f (i)

T (ϕ) = T t=1 α (i) t h (i) t (

Data Energy Term

We adopt the same Bayesian approach presented in section 3.2 to formulate the segmentation task. The equivalent MAP estimation problem is dened in equation (3.13) as minimizing:

E (x 1 , . . . , x n , I) = E shape (x 1 , . . . , x n ) + E image (x 1 , . . . , x n , I) , (3.29)
with respect to (x 1 , . . . , x n ) ∈ Ω n , Ω being the image domain. Recall that the prior term E shape refers to our k-fan shape model and is expressed in terms of pairwise landmark interactions, similarly to equation(2.30):

E shape (x 1 , . . . , x n ) = - (i,j)∈E k γ ij log p x i -x j d .
(3.30)

The energy term E image is based on the image cues, and is dened thanks to the learned classiers. Using a logistic function, we turn the strong classier response into a probability measure as p i f (i)

T (ϕ) = 1 1+exp -f (i) T (ϕ)
. Then, by assuming independence of the feature vectors at dierent image locations, we can write:

E image (x 1 , . . . , x n , I) = n i=1 -log p i f (i) T (ϕ (x i )) = n i=1 -log   1 1 + exp -f (i) T (ϕ (x i ))   , (3.31)
where ϕ (x i ) is the feature vector at the landmark x i location. Let us consider now that each position in the image domain Ω corresponds to a label l ∈ L, L being the labeling set. We will note x i (l) placing the landmark x i at the location labeled by l. Then, minimizing the energy (3.29) is equivalent to solving the labeling problem: 3.32) where the unary and pairwise potentials are dened as: For eciency, we used an approximate optimization methods to solve (3.29), namely the sequential tree-reweighted message passing algorithm (TRW-S) [Kolmogorov 2006]. In TRW-S, a lower bound on the optimal energy is maximized. We observed in our experiments the lower bound value was very often equal to the graph energy at the end of the optimization, meaning that a global optimum was found. The results we obtained are presented in the next section. Considering a relatively small number of candidate points allows solving the problem eciently. The results can be however let down if there are no good landmarks in the candidate pool. We address this limitation in the following.

min l i ∈L 1≤i≤n E (l 1 . . . , l n ) = α n i=1 V i (l i ) + (i,j)∈E k V ij (l i , l j ) , (
V i (l i ) = -log p i f (i) T ϕ x i (l i ) V ij (l i , l j ) = -γ ij log p x i (l i ) -x j (l j ) d . ( 3

Segmentation Outliers

The above formulation can be in some dicult cases insucient to obtain satisfactory segmentations because of the limitation of the detected candidates.

The proposed algorithm relies heavily on the quality of the detectors and the provided costs, and hence a poor generalization performance of the classier may aect the whole segmentation process. To account for the missing correspondences problem, we consider the idea of an additional articial candidate per landmark, as suggested in [Donner et al. 2007a]. The m + 1 th candidate will refer to nding an outlier. The critical aspect is then to dene the corresponding unary and pairwise costs. [Donner et al. 2007a] dene these costs as being proportional to the mean of the one computed for the real candidates.

On one hand, this cost choice introduces an additional parameter (of proportionality) that needs to be adjusted. On the other hand, using the current test image to compute these costs may be unsuitable as these computations will be dependent on the number of candidates. A dierent number of candidates may result in dierent costs, and hence dierent outliers detection. Moreover, the generalization performance of the classier will have a direct impact on the outlier cost in this case. An alternative way to dene these costs relies on their statistical estimation from the training set, where the energy terms of the available solutions (segmentations) can be easily evaluated. Hence, we compute from the training set the mean and standard deviation of each unary and pairwise cost, and set the outliers penalties to the mean plus two standard deviations.

The post-processing of the detected outlier landmarks is another important aspect when handling such dicult cases. Hence, one has to dene a strategy to "repair" the outlier landmark positions, and obtain an allowable segmentation. It is not clear in [Donner et al. 2007a] how this step is accounted for.

[ [START_REF] Seghers | [END_REF] proposed to relax the positions of landmarks that have high shape energy terms to lower shape energy locations. However, the used method is not clearly explained. We provide now a way to solve for the missing outlier positions.

Augmenting the candidate set by an outlier label provides a partial segmentation of the object. In practice, after minimizing (3.32)

{d * kk } i∈O = argmin d kk ∈D k E outliers (d ij ) , (3.34)
where: MISCP algorithm [Seghers et al. 2007a] 93.9(±3.1) 1.49(±0.63) the graph, and the pairwise normalized distance distributions were learned using Gaussian kernel density estimation. The boosted classiers f (i) T were also learned using the 124 images of the training set according to the described scheme in section 3.3.1. The Fig. 3.13(c) shows also that we used for each node regions of interests (ROIs) to speed up the candidate search step. The limits of these ROIs were learned from the training data, as done in [Seghers et al. 2007a]. In order to evaluate the performance of our method, we compared the segmentation results that we obtained using a complete graph against using the k-fan learned graph. In both cases, the maximum number of candidates m was set to 20. As in [van Ginneken et al. 2006] and [Seghers et al. 2007a], we used the overlap coecient Γ = T P T P +F P +F N (where TP stands for true positive, FP for false positive, and FN for false negative) and the mean curve distance error (MDE) to assess the results quantitatively.

E outliers (d kk ) = (i,j)∈O×O (i,j)∈E -log p ij (x i + d ii , x j + d jj ) + (i,j)∈O×I (i,j)∈E -log p ij (x i + d ii , x j ) . ( 3
They are summarized in Table 3.2.

We also report the outlier label, and evaluated its impact on the segmentation results. Table 3.2 and the boxplots in Fig. 3.17 show that the quantitative performance of the segmentation algorithm is improved when using the complete graph and handling outliers.

We show in Fig. 3.15 and Fig. 3.16 examples of images where outliers were detected, and where the segmentation using the complete graph was improved. The examples in Fig. 3.15 also suggest that the performance of the k-fan graph can be improved by accounting for the outliers. However, as we can see in Fig. 3.16, the method performs better in the case of the complete graph. This observation is conrmed by the quantitative results in Table 3.2 and the boxplots in Fig. 3.17. They show that the computed performance indicators are slightly impaired when using outliers with the k-fan graph. This dierence in the performance of the k-fan graph and the compete graph may be explained as follows: selecting a landmark as being an outlier can be thought as eliminating it from the graph, and destroying all its connections. Suppressing these geometric constraints is clearly more critical for the k-fan graph than for the complete graph, especially if landmarks of the base k-cliques are labeled as outliers. Thereby, we alleviate the need of an extensive pool of training data, which is needed in higher dimension estimation.

Then, observing the computational burden induced by a complete graph structure, and the redundant information that is carried by the complete set of edges, we tackled the problem of graph structure learning. Seeking eciency and model compactness, we approximated the complete graph by a k-fan graph which structure is determined from a training set of shapes using manifold learning and unsupervised clustering techniques. We showed that this approximation is tight in the case of the exact matching problem, by using concept of graph rigidity. We qualitatively veried this strong geometric property in experiments, where we recovered 2D and 3D shapes from random points thanks to the shape prior constraints.

The obtained k-fan graph belongs to the more general class of chordal graphs, and hence inherits the perfect elimination order property. Theoretically, this property allows the use of dynamic programming to exactly solve optimization problem dened on the graph in polynomial time. Unfortunately, the dimensionality of the applications that we considered made this approach intractable in practice.

In the next step, we proposed dierent methods to combine the shape prior with various image cues in a maximum a posteriori MRF optimization framework towards image segmentation. The obtained energy minimizations are cast into labeling problems that are eciently solved using recent approximate optimization techniques. In the experiments, we noticed that we were able to recover the global optima in some cases.

More specically, we integrated statistical regional image information in the objective energy as unary potentials using a Voronoi decomposition of the image domain. We proposed an iterative coarse-to-ne procedure to solve the problem, in a deformable model fashion. We also considered the use of landmark-based cues through boosted classiers in a Bayesian framework. To account for poor generalization performance of the classiers, we proposed a method that detects outliers during inference, and repairs their positions according to the shape prior. We also proposed a tracking algorithm that uses edge-based image information and combines the static shape model with dynamic priors. In addition, we applied graph matching techniques in the image registration setting to select optimal mutually-salient corresponding landmarks. Although we used these dierent image cues seperately in the applications that we presented, our framework is exible and allows the integration of cues of dierent nature at the same time.

Despite the encouraging results that we obtained, we think that our approach has some limitations. Our shape model handles naturally partial object occlusions. However, it does not explicitly account for self-occlusions. Hence, our region-based segmentation algorithm may be unsuitable to process objects with overlapping parts. Moreover, the scale-invariance of our model is enforced through normalization, which supposes knowing the object scale.

We resort in practice to approximations in order to estimate the unknown scale values. Computing the regional image information as unary potential is another approximation that we made.

Future Research Directions In our future research work, we will endeavor to address the above mentioned downsides. The use of pairwise terms to model the regional support will reduce the shortcomings introduced from the Voronoi approximation. In particular using the Stokes' theorem, regional integrals could be expressed as inner products on the local 2D segments forming the 2D contour. Then, an explicit and direct estimation of the regional properties would be feasible. Furthermore, the combination of image-based (edges or regional support) and landmark based segmentation will be also feasible. The idea will be to have a two layer graph with interconnected variables, one acting on the landmarks and the second acting on the image. The interconnections will guarantee that the solution to both problems will be unique.

The case of 3D will be more complex since it will require higher order MRFs. The Green's theorem will decompose the regional image term on the boundary elements (triangles in 3D) in an exact manner. However, such a formalism will inherit higher order MRFs. This is also a natural extension of our approach since it will introduce the notion of complete shape invariance to the process. The use of third order potentials inherits invariance with respect to scale changes, translation and rotation, while being able to cope with arbitrary image support. The main challenge in this context will be to minimize the energy but either the use of methods that convert higher order MRFs into pairwise MRFs or the ones based on dual decomposition and master-slave approaches will be considered.

As far as optimization is concerned, we have already mentioned that beneting from the structure of our model through the use of dynamic programming is cursed by the dimensionality of the considered problems. Developing optimization algorithms that take into account the particular structure of the k-fan graph would represent an interesting research direction.

The case of spatio-temporal shape modeling is also of great interest in particular for medical image analysis when studying organs with periodic motion.

Inferring the temporal connectivity of the graph can be extremely challenging and benecial as well since the notion of k-fan model is not applicable.

Re-formulating the inference process as an unsupervised clustering problem with terms accounting for spatial and others for temporal connectivity might be the right direction to go.

Last but not least in all above approaches the selection of the model was based on purely geometric properties and characteristics. Then, during inference, the model was associated with image cues. A promising research direction would bridge the gap between these two aspects. An interesting modeling scheme would not only learn the structure of the graph from the data, but it would also benet from training image information to dene the control points as salient image points. Combining shape and appearance learning is a possible extension of our work. Introduction

Deformable image registration is one of the most challenging problems in medical imaging. We distinguish two classes of prior work: iconic methods and geometric methods. Iconic (or voxel-based) methods optimize a similarity criterion (for instance mutual information) that is dened on the images in order to nd an optimal transformation with respect to this criterion. These methods generally perform well for intra-modal registration. However, when images belong to dierent modalities, such approaches can lead to failure.

In the geometric methods, landmarks are extracted from the images and the transformation that maps the source landmarks to the target landmarks is estimated.

Establishing landmark correspondences is of a great importance in many landmark-based non-rigid registration algorithm. This task is often performed in two successive steps, namely, landmark detection (e.g., [START_REF] Kadir | Timor Kadir and Michael Brady. Saliency, Scale and Image Description[END_REF]) and landmark matching (e.g., [Joshi & Miller 2000]).

One can formulate two criticisms to these separate two-step methods.

First, landmarks are often separately detected in two images. Although having individual saliency property, they might be not discriminative for matching. As pointed out in [Ou & Davatzikos 2009], salient points in one image are not necessarily present, or uniquely present in the other image.

Therefore, an alternative approach that detects pairs of landmarks with a preferably unique matching score across images would yield better registration performance.

Second, individually optimal correspondence for each landmark is not necessarily globally optimal for the deformation eld. For instance, such an approach might result in displacement vectors at two nearby landmarks that point to completely opposite directions, leading to a dense deformation eld with self-intersections. This undesired behavior can be avoided by seeking a trade-o between the global smoothness of the resulting transformation and the landmarks similarities when establishing correspondences.

In this section, we tackle the above-mentioned drawbacks. First, instead of extracting individual interest points, we simultaneously detect and match that, the template image space Ω 2 is not partitioned, because at this stage, no transformation is conducted, therefore no corresponding regions can be assumed.

As we observed earlier, the top ranked pair is a locally optimal pair, that might be not globally optimal with respect to a smooth deformation eld be- In the rst case, one can cite for instance methods like the Kanade-Lucas-Tomasi tracker [Tomasi & Kanade 1991], correlation-based methods [Shi & Tomasi 1994], the mean-shift algorithm [START_REF] Comaniciu | [END_REF] and its numerous variants [Collins 2003]. These methods use a similarity criterion to compare the previous appearance of the object with possible candidates in the new images towards recovering the most probable position. These methods are a good compromise in terms of computational complexity but may fail with severe deformations. The tracking task has been formulated in dierent optimization frameworks. For instance, [START_REF] Paragios | [END_REF] proposed a solution to the problem in the level set framework. [Cremers 2008] adopted a variational approach as well incorporating dynamical shape priors. Other approaches dynamic shape prior, and then, we formulate the edge-based energy term. In the last subsection, we show the tracking results that we obtained.

A.2.2 Weak Edges -Static and Dynamic Shape Priors

Let I 1 , . . . , I N be a sequence of image frames where an object is sought. Let x i (t) be the position of the i th control point in the image frame I t . We extend the MAP estimation problem dened in section 3.2 to the temporal domain by considering an energy term that accounts for the interactions between the control points in the current frame, and those in the previous ones. Then, we rewrite the equation (3.13) as minimizing:

E (x 1 (t), . . . , x n (t), I(t)) = βE dyn (τ, x 1 (t), . . . , x n (t))

+ E shape (x 1 (t), . . . , x n (t))

+ αE image (x 1 (t), . . . , x n (t), I(t)) ,

(A.7)
where E dyn (τ, x 1 (t), . . . , x n (t)) is a dynamic shape prior of order τ , as explained in the following, and α and β are weighting parameters that modulate the contribution of each energy term. Recall that the static shape prior energy E shape is related to our k-fan shape model and expressed as pairwise interaction between the landmarks (2.30). We consider such pairwise interactions in the temporal domain as well.

Dynamic Prior

Let us also consider the inter-point normalized distance between two landmarks in dierent frames separated by an index τ . We suppose that the global motion is accounted for through translation of the gravity center of the object:

d ij (t, τ ) =
x i (t) -x(t) -x j (t -τ ) -x(t -τ ) d(t, τ ) , x k (t) -x(t) -x l (t -τ ) -x(t -τ ) , and x(t) (resp. x(t -τ )) being the gravity center of the object at the frame t (resp. at the frame t -τ ). Given a training set of image frames where the object of interest is segmented, and correspondences between the control where we use the notation x n+1 (t) = x 1 (t). Although this low-level feature is sensitive to noise, and would often fail, its association with the static and dynamic prior overcomes its usual limitations. To summarize, the objective energy (A.7) is the sum of unary terms that account for the dynamic prior, and two types of pairwise terms: an edge-based one and a static prior-related one.

To minimize the objective energy, we consider the same approach described in section 3.2.3. At each image frame I t , we initialize the shape model using the result of the previous frame. We consider a quantization of the search space, and we solve iteratively several labeling problems, where each label refers to a possible control point displacement. The obtained MRF minimization is an NP-hard problem. We consider approximate and computationallyecient solutions given by the FastPD algorithm [Komodakis et al. 2008b].

We also use a coarse-to-ne search to speed up computations.

  Segmentation d'Images avec des Champs de Markov et Modélisation Statistique de Formes Résumé : Nous présentons dans cette thèse un nouveau modèle statistique de forme et l'utilisons pour la segmentation d'images avec a priori. Ce modèle est représenté par un champ de Markov. Les n÷uds du graphe correspondent aux points de contrôle situés sur le contour de la forme géométrique, et les arêtes du graphe représentent les dépendances entre les points de contrôle.

  Mots clés : Modélisation de formes, Segmentation, Champs de Markov, Optimisation discrète Paragios for his guidance and support during the four years of my PhD. He gave me the opportunity to explore a challenging research topic, and constantly oered to me his help and encouragements. He always had an open door for discussions whenever I needed support or advice. This PhD was a great experience, and it would not have been possible without his insightful input and continual trust. I was honored to have Prof. Horst Bischof, Dr. Leo Grady, Dr. Grégoire Malandain, Prof. Mads Nielsen, Dr. Patrick Pérez, Prof. Christian Saguez and Prof. Ramin Zabih as members of my thesis committee. I would like to thank them for traveling from far away to attend my defense, despite their busy schedule, and for providing me with timely feedback. Prof. Horst Bischof and Dr. Patrick Pérez reviewed this manuscript, and their pertinent comments contributed to improve its quality. I owe my deep gratitude to all of them. I also would like to thank Prof. James Duncan for having welcomed me at the Biomedical Engineering laboratory at Yale University, during two months in the summer of 2007. Under the supervision of Prof. Nikos Paragios, the work that was conducted during this stay was very important for the remainder of my PhD studies. This thesis would not also have been possible without the collaboration and help of my colleagues. I particularly acknowledge Dr. Nikos Komodakis and Dr. Georg Langs for their discussions and for kindly sharing their code. My thanks also go to Yangming Ou and Yun Zeng for their contribution to this work. I am grateful as well to the members of the Applied Mathematics laboratory of Ecole Centrale Paris, and especially to the members of the computer vision and medical imaging team.

[

  [START_REF] Cremers | [END_REF] introduced shape priors into the Mumford Shah functional. They used a spline parametrization of the contours and learned a multivariate Gaussian distribution over the control point positions. They introduced in[Cremers et al. 2002b] invariance to translation rotation and scale changes in the variational framework. To overcome the non-invertibility of the covariance matrix in the case of a small training set, the latter was regularized by replacing the zero eigenvalues by the smallest eigenvalue. Hence, a shape variation that is not represented by the training set has a smaller probability than any shape variation encoded by the training set.

  that constrained the deformations of the contour during the search. An average contour model have been used to represent shapes in [Chen et al. 2002], after rigid alignment of the training examples. In the case of a class of shapes that exhibits large variations, the training examples were classied in dierent clusters and an average shape model was considered for every cluster. More recently, [Schoenemann & Cremers 2007] developed a graph-theoretical algorithm for segmentation and alignment of shapes that produces a global minimum to the optimization problem. The algorithm uses an elastic shape prior that is invariant to scale and translation changes, and is able to match a single object template to the image based on edge information. The algorithm can benet from parallelization. It is also possible to induce rotation-invariance by matching successively the image to the template rotated by dierent angles. 1973]). The best match of the model to an image is the one that gives good condence for each part while satisfying inter-parts constraints. They require the object model to have a tree structure, which allows the use of dynamic programming to minimize the objective function in linear time (rather than quadratic time, thanks to a restriction on the pairwise cost function) and obtain a global optimum. A person model and a car model are designed and applied for detection. The paper introduces also a Bayesian formulation of the problem as a maximum a posteriori (MAP) estimation.

(

  individual statistics of NCL over a training population) is dierent from what they used as we will show later.

Fig. 2 .Figure 2 . 1 :

 221 Fig. 2.1 illustrates the previous denition, and Fig. 2.3 shows, among other chordal graphs, examples of k-trees. As it will be detailed later,

  (a) The steps of building a 3-tree are shown. In each step, the current base 3-clique is depicted in red. (b) The obtained graph has 2 base 3-cliques. A simplicial elimination order is given.

Figure 2 . 2 :

 22 Figure 2.2: (a) A chordless graph with 8 vertices. (b) A chordal graph with 8 vertices, and a simplicial elimination order (SEO).

Figure 2 . 3 :

 23 Figure 2.3: Several chordal graph structures with 8 vertices, and a SEO for each graph. (a) A tree. (b) The complete graph. (c) A triangulated polygon and its dual graph. (d) A 2-fan graph -the reference vertices are shown. (e)

Figure 2

 2 Figure 2.4: A synthetic example illustrating the use of shape maps and unsupervised clustering. (a) The generated synthetic data: 4 boxes with varying aspects rotate around each other independently. (b) The shape maps embedding of the data is shown using the rst 3 dimensions. We can see that 4 clusters appear, which corresponds to the generated data. (c) The unsupervised LP-based clustering (section 2.2.4) is performed using the obtained diusion distance. It nds automatically the number of clusters and correctly identies them, which permits to separate the boxes in the data.

Figure 2

 2 Figure 2.5: Hand model construction. (a) The considered training set. (b)

Figure 2 .

 2 Figure 2.6: Hand shape prior applied to random points.(a)(d) Deformation of a random point cloud according to the shape prior term.

Figure 2

 2 Figure 2.7: Corpus callosum model construction. (a) The considered training set. (b) The projection of the control points on the rst 3 shape map dimensions. The density of the landmarks is estimated by a Gaussian kernel and is color-coded. (c) Control points clustered in 13 clusters: centers are represented by squares.

Figure 2 . 8 :

 28 Figure 2.8: Corpus callosum shape prior applied to random points.(a)(d) Deformation of a random point cloud according to the shape prior term.

  2.7(b). Next, we performed the unsupervised LP-based clustering method that provided 13 clusters as shown inFig 2.7(c).

Figure 2 . 9 :

 29 Figure 2.9: Joint lungs model. (a) The considered training set. (b) The projection of the control points on the rst 3 shape map dimensions. The density of the landmarks is estimated by a Gaussian kernel and is color-coded. (c) Lungs landmarks are clustered in 12 clusters: centers are represented by squares.

Figure 2 . 10 :

 210 Figure 2.10: Joint lungs shape prior applied to random points.(a)(d) Deformation of a random point cloud according to lungs the shape prior term.

Figure 2 . 13 :

 213 Figure 2.12: A separate model of the left lung. (a) Shape maps embedding (left) and unsupervised clustering (right). We obtain 9 clusters. (b) Applying the learned shape model to random points.

Figure 2 . 14 :

 214 Figure 2.14: Left ventricle shape prior applied to random points.(a)(d) The mesh is updated through thin plate spline (TPS) interpolation of the control point (in green) displacements. These displacements are guided by the shape prior.

  will exploit the normalized chord length shape model towards model-based segmentation. The underlying graph structure is determined from the data, and the shape model parameters are learned from the training shape population as well. Given such a model, one has rst to introduce data support and then ecient optimization techniques for the inference of the shape parameters, given an unseen image (an image that does not belong to the training dataset). Our motivation is to propose a exible framework that allows the use of dierent image cues, and permits ecient optimization thanks to the recent developments in the eld.

  by its parameters vector b which is the orthogonal projection of the residual from the mean on the eigenbasis E, or b = E T (s -s). An appearance model per landmark is associated to the shape model, and is used to enforce image-driven prior knowledge during search. This appearance model is related to the boundary of the object: intensity proles, centered in the point location, and orthogonal to the contour (and hence capturing internal and external object features) are extracted from the training data. Mean proles and covariance matrices are then computed, and the quality of a new prole

  Moreover, other non-linear methods apart from Gaussian mixtures were used to model the appearance features.[START_REF] De Bruijne | [END_REF]] propose a method based on k-nearest neighbors (kNN) to evaluate the probability of a given prole. During the training phase, proles are sampled inside and outside the object in the vicinity of the landmarks positions. Hence the pool of proles contains "true" proles and "false" proles, corresponding respectively to correct and wrong landmark positions. Then during the search phase, the kNN of a candidate prole assesses its quality by voting, and a probability value can be determined as the fraction of "true" votes. A similar approach was adopted by[START_REF] Van Ginneken | [END_REF] to evaluate the goodness of candidate positions during search. The dierence from[START_REF] De Bruijne | [END_REF] lies in the use of a weighted kNN voting (the weights are inversely proportional to the distance to the landmark position), which establishes for a given point the probability of belonging to the inside or to the outside of the object. Then, the optimal position is dened by minimizing a cost function that combines such individual probabilities along proles orthogonal to the contour and centered at the tested positions. The same rationale is applicable with dierent classiers. Successful methods are presented in[Li et al. 2004a] and[START_REF] Li | Shape Parameter Optimization for AdaBoosted Active Shape Model[END_REF] and benet from the strengths of the Adaboost algorithm[START_REF] Freund | [END_REF]].Alternative ASM search schemesOther eorts to develop the ASM algorithm mainly considered the search part. [de Bruijne & Nielsen 2004] use a point distribution model (PDM) in a stochastic inference framework. They

  algorithm is fairly dierent from its ASM counterpart.The AAM search corresponds to an image interpretation problem. It aims to minimize the dierence between the given image and a version synthesized by the appearance model. An instance of the model in the image is dened by a similarity transformation, an appearance parameter vector c and a texture transformation parameter vector. All these parameters are described in a unique parameter vector p. Let g be the model-generated texture, and let g img be the part of the image texture inside the shape region. Then, minimizing the residual r(p) = gg img amounts after a rst order Taylor expansion to expressing the parameters update as a function of the residual, or: δp = -R.r(p) .

  an algorithm that converges to a local minimum of the energy function. In the experiments, they successfully apply their method and show examples of topology-preserving segmentations. Although the last three methods account for shape characteristics, their prior remains generic, and they do not allow specifying a particular class of shapes. This is exactly what other researchers worked towards. [Slabaugh & Unal 2005] use novel terminal weights as a formulation of an elliptic shape prior. This geometric primitive adds regularity constraints to the segmentation. They propose an iterative framework, where a binary elliptic mask is repeatedly formed knowing the current ellipse position in the image, mean image intensity inside and outside the ellipse are computed,

  More recently,[Ben Ayed et al. 2009] introduce a novel graph-cut-based algorithm to segment 2D left ventricle (LV) images in magnetic resonance (MR) sequences. In their method, the user manually segments the rst frame, by drawing a single contour that partitions the image into two disjoint sets (the object and the background), and the algorithm propagates the segmentation to the remaining ones. They argue that this user intervention alleviates the need for intensive shape training of the geometric and photometric features.

3. 1

 1 gives an example of a Voronoi diagram. Note that the Voronoi tessellation is the dual of the Delaunay triangulation. Hence we can write Ω = ∪ n i=1 Ω i .

Figure 3 . 1 :

 31 Figure 3.1: Voronoi decomposition of the model domain.

Fig 3 . 2

 32 the way we compute the unary potentials(Fig 3.2(a)) and the pairwise potentials(Fig 3.2(b)). On one hand, equation(3.23) expresses the approximation that we make with respect to the data term computations. The Voronoi cell associated with the control point x i displaced by d l i is estimated by translating the cell O i by the same displacement vector, which is noted as O i + d l i .

Figure 3 . 2 :

 32 Figure 3.2: Computing the energy terms: grey squares refer to the image grid, and numbered disks refer to the control points possible labels. The label 1 refers to the current landmark position, i.e. a null displacement. (a) Unary data term: the cost of the label 4 of the landmark x i is evaluated by translating the corresponding Voronoi cell O i (depicted in red) by d 4 and then computing the integral in (3.23). The translated cell O i +d 4 is shown in green. (b) Pairwise shape term: computing the cost V ij (4, 2) is done by applying the displacements d 4 and d 2 respectively to x i and x j and then estimating the resulting distance d ij (4, 2).

Figure 3 . 3 :

 33 Figure 3.3: Segmentation algorithm inner loop: pairs of images showing the evolving contour (left image of the pair) and the Voronoi diagram (right image of the pair). An "object" Voronoi cell is emphasized in yellow. Two cases: a simple case (a) & (b) and a case with missing parts (c) & (d). (a) & (c) Initialization: only 6 landmarks are used. The resulting transformation is projected on the similarity transformations space. (b) &(d) Convergence using 56 points.

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Easy synthetic case. (a) Initialization. (b)(c) Intermediate inneriterations. (d) Convergence.

Fig 3. 3 showsFigure 3 . 6 :

 336 Figure 3.6: Synthetic tracking with occlusion. (a) Initialization. (b)(j) Results on subsequent frames.

Figure 3 .

 3 Figure 3.7: Boxplots of dice overlap coecients comparing our method to AAM.

  dimensional medical image analysis setting, we used a dataset of 28 computed tomography (CT) cardiac volumes, having an approximate mean size of 512 × 512 × 250 voxels, where the voxel size is about 0.35 × 0.35 × 0.35mm 3 . (a) Database Examples: Successful segmentations (b) Finger Collusion -Missing Part Examples. Two rst images: dicult examples because of ngers collusions. Three last images: segmentation of hands with missing parts. (c) Severe Noise Added: The prior knowledge highly contributes in correctly segmenting very noisy images. (d) Video Frames -Partial Occlusions. Real video frames: cluttered background and occlusions. (e) AAM results: succeeds with the learning examples but fails with occlusions. Initialization on the left -result on the right.

Figure 3 . 8 :

 38 Figure 3.8: Model-based segmentation of the hand. Initialization is shown in white, segmentation in red, and the nal control points positions in blue.

Figure 3 . 9 :

 39 Figure 3.9: Segmentation results: 3 unseen examples (do not belong to the training set). Initialization in yellow, shape after ane transform in blue, nal segmentation after TPS deformation in red. Random Walker result in green.

Figure 3 . 10 :

 310 Figure 3.10: Our model: a deformable shape associated with control points. (a) In red: control points P 4 (in the apical and basal parts) used to dene ane transforms of the shape. In green: the set of control points P 90 that denes the TPS deformation. (b) The apical control point with the associated Voronoi cell, intersected with the blood pool and the myocardium.

Figure 3 . 11 :

 311 Figure 3.11: Example of posterior-anterior chest X-Rays with manual expert segmentations.

  are usually easy to train, but have quite poor individual performance. They are "boosted" by combining them, which yields a "strong" classier f T . Boosting works iteratively by trying to improve the function f T at each step. During each iteration, the best weak classier h t is added with

Figure 3 . 12 :

 312 Figure 3.12: Features computation using a lter bank.

  tives" is learned from weighted training data, and a corresponding step size is determined. Then, poorly classied training examples are identied. Subsequently, their weights are increased in order to emphasize their importance in the next step. A typical choice for the weak learner family is decision stumps. Algorithm 3 summarizes these dierent steps. Algorithm 3 Adaboost 1: Input: N training examples {(ϕ 1 , y 1 ) , . . . , (ϕ N , y N )} where ϕ i ∈ Φ and y i ∈ {1, -1}.

  ϕ). The obtained classier responses are incorporated into a prior-based segmentation algorithm as explained in the following.

. 33 )

 33 However, searching in the whole image domain Ω increases the complexity of the problem. Therefore, we restrict the search in practice to the image locations where the classiers provide their highest m responses (for each landmark), m being a parameter that is xed by the user. Hence, the problem boils down to nding the best conguration among the top m control points candidates.

(Figure 3 . 13 :

 313 Figure 3.13: Modeling the right lung using the training set. (a) The obtained shape maps embedding. (b) The unsupervised linear programming-based clustering lead to 7 clusters. (c) A region of interest is used for each landmark to speed up the segmentation algorithm.

Figure 3 . 14 :

 314 Figure 3.14: Segmentation results of the right lung, using the complete graph (top row) and the k-fan graph (bottom row).

Figure 3 . 15 :

 315 Figure 3.15: Eect of the outlier label on the segmentation results. Outlier landmarks are in red. The crosses correspond to the top candidate position that is used as initialization for the perturbation algorithm. (a) Complete Graph (b) Complete Graph + Outliers (c) k-fan Graph (d) k-fan Graph + Outliers.

Fig. 3 .Figure 3 . 16 :

 3316 Fig. 3.14. The quantitative and qualitative results back our claims regarding the learned graph structure: we obtain a sparse model, and we maintain a performance that is equivalent to the complete graph. Moreover, as described earlier, we learned from the training set unary and pairwise costs for the

Figure 3 . 17 :

 317 Figure 3.17: Quantitative evaluation of the segmentation of the right lung: boxplot representation of the distributions of (a) the overlap coecient Γ and (b) the mean curve distance error (MDE). The box lines represent the lower quartile, the median, and upper the quartile values. The whiskers extend to a range of 1.5 times the interquartile range. The values beyond this range are drawn with a "+" sign.

Figure A. 2 :

 2 Figure A.2: Examples of dierent voxel pairs having dierent mutual-saliency values.

  applied to nding mutually-salient landmark pairs in brain and cardiac images. Results in the following sub-sections aim to demonstrate the advantage of unifying landmark extraction and landmark matching in our framework. Results for Detecting Mutually-Salient Landmark Pairs Fig. A.5shows a landmark pair detected across subjects based on the mutual-saliency measure. Alternative methods that consider only image intensity to describe voxels and that use edge/corner detectors to extract landmarks separately in the images would probably fail to detect this pair of voxels: these voxels are simply not edges or corners or surface boundaries. We succeed in detecting this pair for two reasons. On the one hand, Gabor attributes are powerful descriptors that provide a distinctive characterization of voxels, as shown in similarity maps (c,d). On the other hand, the mutual-saliency measure eectively quanties the matching uniqueness, as shown in similarity maps (e,f ).Overall, this example shows the advantage of detecting landmarks in pairs other than one-by-one and separately from the two images.

Figure A. 4 :

 4 Figure A.4: Dense deformation elds generated by (a) M1 no MRF regularization and (b) M2 with MRF regularization.

Figure A. 5 :

 5 Figure A.5: An example landmark pair (denoted by red and blue crosses) detected based on mutual-saliency measure. (a) Source image I 1 (b) Target image I 2 . Similarity maps are generated (c) between the red cross point and all voxels in Ω 1 (d) between the blue cross point and all voxels in Ω 2 (e)between the red point and all voxels in Ω 2 and (f ) between the blue cross point and all voxels in Ω 1 .

  Dynamical systems are a promising alternative to encode motion dynamics. One can cite numerous examples using such model: Kalman ltering[Kalman 1960], condensation[START_REF] Isard | [END_REF]], or multiple hypotheses testing and particle lter tracking[Arulampalam et al. 2002]. These methods benet from an explicit model that represents the motion of the object, and can thereby perform better than blob-based methods.One should also mention methods that aim to track articulate models.In such a context, objects are represented with parts, and then constraints between the relative positions of these parts are introduced [Felzenszwalb & Huttenlocher 2005] [Black & Jepson 1998] [Urtasun et al. 2006] [Wang et al. 2008]. These methods can be very ecient but assume an explicit hierarchical representation of the model, and impose constraints on their dynamics. Therefore, their use in the context of tracking arbitrary highly deformable objects is not straightforward.

  Figure A.6:Frames extracted from a tracking sequence using our method with the static prior. Note that the tracking quality is better when using the dynamic prior for the frame 66 for example (see Fig. A.7(b)): here one leg is missed.

Figure

  Figure A.7: Frames extracted from a tracking sequence using our method with the dynamic prior. Note that the tracking quality is better than in the static prior case for the frame 66 for example (see Fig. A.6(b)).

Figure A. 9 :

 9 Figure A.9: Frames extracted from a tracking sequence using our method with the dynamic prior.

  , . . . , v n of its vertices such that v i is simplicial in the subgraph induced by {v i , . . . , v n }.

	] discusses ge-
	ometric aspects as well as graph-theoretic optimization aspects. The k-trees
	being chordal graphs, they permit to solve general NP-hard problems in poly-
	nomial time (precisely, they use the junction tree algorithm to solve a MAP
	problem).
	Chordal Graphs
	Let us rst give a few denitions:
	Denition 2.3. Chordal Graph
	A graph is chordal if it does not contain any chordless cycle of length at least
	4. A chord is an edge that connects two non successive vertices in a cycle.
	Denition 2.4. Simplicial Elimination Ordering (SEO)
	A vertex of a graph is called simplicial if its neighbors form a clique (a clique
	on n vertices is K n , i.e. a complete subgraph with n vertices). A simplicial
	elimination ordering of a graph (SEO), or a perfect elimination scheme, is an
	ordering v 1 An important graph theory result is that a graph is chordal if and only if
	it has a perfect elimination scheme. An example of a chordless graph and a
	chordal graph, with a possible SEO, are given in Fig. 2.2. Other example of
	chordal graphs are shown in Fig. 2.3.
	[Felzenszwalb 2005] introduces a representation of 2D shapes (with
	piecewise-smooth boundaries and no holes) using triangulated polygons (by
	the means of constrained Delaunay triangulation on the vertices that form

Table 2

 2

		Hand example LV example
	Complete Graph runtime (s)	34(±4)	51(±5)
	k-fan Graph runtime (s)	11(±1)	26(±2)
	Number of Labels	361	105
	# edges Complete Graph	1540	4005
	# edges k-fan Graph	459	1079

.1: Running time comparison between the complete graph and the k-fan graph.

  The above-mentioned class of methods try to optimize the dual of the LP relaxation. Interestingly, this dual energy is a lower bound on the MAP MRF energy.

	show that the QP relaxation and the SOCP relaxation are equivalent in terms
	of tightness. Furthermore, they show that they do not compete favorably with
	the LP relaxation, meaning that the latter leads to better approximations of
	the MAP-MRF problem.
		[Komodakis & Paragios 2008] present a work that deals with the relax-
	ation tightness. They act on the lower bound of the MRF energy by dening
	tighter LP-relaxations, through the use of a new set of pairwise potentials that
	satisfy certain constraints, called virtual potentials. They give a characteri-
	zation of dual LP energies that are not tight, as the existence of the so-called
	inconsistent node cycles with respect to the labeling. They propose an iter-
	ative algorithm to repair these inconsistent cycles, by successively adjusting
	the virtual potential to appropriate values, and applying the augmenting DAG
	algorithm
		Hence, one can hope that
	by solving the LP, a good approximation to the original integer program
	can be obtained. This observation raises the question of the tightness of
	the relaxation. Indeed, a loose relaxation would lead to a poor approxi-
	mation. Researchers have tried other types of relaxations, besides the LP
	one.	For instance, a quadratic programming (QP) relaxation is used by

Table 3 .

 3 1: Comparison of our method with the Random Walker Algorithm

	Correctly Segmented Voxels Incorrectly Segmented Voxels Segmentation Time
	(True Positives)	(False Positives)	(3 GHz, 3 GB)
	85.12%±7.3	15.3%±10.2	90s±10
	axial data slices were downsized to 128 × 128 pixels.

  by computing the rst statistical moments of F in local neighborhoods.[Seghers et al. 2007a] modeled the feature vector probability function using a multivariate Gaussian distribution for each blurring level σ. Our derivation is slightly dierent as we build for each control point of the object a "boosted" classier fSchapire 1996]. To this end, at each control point location x i , we extract a patch of size 17 × 17 from the feature image F, amounting to a K-

	(i) T , using Adaboost
	[Freund & (i) j , y

dimensional descriptor ϕ i ( with K =8959). Hence, given a set of N p training images, where the positions of the landmarks are known, we obtain for each control point x i a set of "positive" examples of size N p . We also extract N n feature vectors from the background to form a set of "negative" examples. Then, using the training set of N = N p + N n examples (ϕ

  , most of the landmarks positions are determined, and only a few of them are set to be outliers. Then, the next step consists to locate these outliers according to the shape prior, and to the set of xed (already solved) landmarks. Let O be the set of outliers, and let I be the set of xed landmarks. Then iteratively, we search the optimal set of displacements d *

i , i ∈ O of the outlier landmarks such that a shape energy is minimized. More formally, we consider the approach presented in section 2.2.6, and more precisely the equation

(2.34)

. Let D i = {d i1 , . . . , d ir }, i ∈ O be a quantization of the possible displacements d ii . Hence, we iteratively look for:

Table 3 .

 3 2: Overlap coecient and Mean Distance Error between the manual gold standard and the obtained right lung segmentation. Mean value over the testing set of 123 images and standard deviation are reported.

		Γ (%)	MDE (pixels)
	Complete Graph	94.77(±2.37) 1.39(±0.64)
	k-fan Graph	94.74(±2.31) 1.40(±0.62)
	Complete Graph + Outliers	94.95(±2.48) 1.34(±0.65)
	k-fan Graph + Outliers	93.87(±3.74) 1.62(±0.87)

Table A .

 A 1: Radiologist's evaluation on results generated by two methods (M1 and M2 for without and with MRF optimization).(or ) means "obviously better (or worse)", > (or <) means "slightly better (or worse)", and ≈ means "almost equivalent". .s. M2). Four datasets are used: Dataset1 (size 256 × 256 × 171) and Dataset2 (size 192 × 236 × 171) for two dierent sets of intra-modality brain MR images across-subjects; Dataset3 (size 192 × 236 × 171) for a pair of multi-modality brain MR images across-subjects; and Dataset4 (size 150 × 150 × 49) for a pair of intra-modality cardiac MR images across-subjects. For each dataset, 25 random pairs out of all the nally selected pairs (typically hundreds) are evaluated in terms of uniqueness and accuracy of correspondences. Evaluation results in Table A.1 show the advantage of incorporating MRF regularization (M2) in all datasets. tracking is one of the fundamental problems of motion analysis in computer vision. Despite enormous work in the eld, the case of highly deformable objects is still an open problem. Prior work on tracking involves blob-based appearance methods, dynamical systems and part-based models. We briey mention some examples of these lines of work.

	A.2	Dynamic Shape Prior and Tracking
	A.2.1	Introduction
			M1	M1 >	M1 ≈	M1 <	M1
			M2	M2	M2	M2	M2
		Dataset1			68%	12%	20%
		Dataset2		4%	80%	16%
		Dataset3	4%	4%	68%	20%	4%
		Dataset4		4%	60%	24%	12%

vObject

Chordal graphs will be dened subsequently in denition

2.3 (section 2.1.3), and the notion of rigidity will be clearly introduced in denition 2.5 (section 2.2.2)

http://www.csd.uoc.gr/~komod/FastPD/

http://vision.middlebury.edu/MRF/code/

[START_REF] Ravikumar | [END_REF], and a second order cone programming (SOCP) relaxation was proposed by[Muramatsu & Suzuki 2003] for the binary case and extended to the general case by[Kumar et al. 2006].[START_REF] Kumar | [END_REF]

This is a notation where no dierence is made between the random variables and their values.

This dataset is courtesy of Siemens Corporate Research

Acknowledgements

I would like to express my sincere gratitude to my advisor Prof. Nikos

• Chapter A.2 is based on a joint work with Yun Zeng.

A.1 Graph Matching for Registration

A.1.1 a number of mutually-salient candidate landmark pairs, i.e., pairs that are uniquely corresponding to each other across images. This pairs nding step relies on the mutual-saliency measure introduced in [Ou & Davatzikos 2009],

which assesses the matching uniqueness of a pair of voxels. In order to alleviate the second limitation and seeking global optimality of the correspondences, we formulate the problem in the Markov-random-eld (MRF) framework. We combine similarity measures with geometric constraints in an energy functional. Its minimization leads to consistent landmark displacements, thereby avoiding potential self-intersections in the resultant deformation eld. Radiologists evaluated the results that we obtained, suggesting the advantage of the two components in our framework.

The remainder of this section is organized as follows: we present our framework in section A.1.2 and provide our results in section A.1.3.

A.1.2 Detecting Landmark Pairs

In this section, we rst give a brief denition of mutual-saliency measure [Ou & Davatzikos 2009]. Then, we describe the two components of our framework: the detection of mutually-salient landmark pairs, and the MRF formulation to nd globally optimal correspondences.

Denition of Mutual-Saliency Measure

Given two images I 1 : Ω 1 → R and I 2 : Ω 2 → R in the 3D image domains

they are similar to each other and only to each others in their neighborhoods.

In the perfect case that is shown in Fig. A.1, the similarity map between u and all voxels in the neighborhood of v should exhibit a delta-shaped function centered at v. We denote the mutual-saliency between u and v as µ(u, v).

We dene this quantity as the ratio of the mean similarity in the central part of the neighborhood of v (denoted as CN (v)), by the mean similarity in the peripheral neighborhood of v (denoted as P N (v)), or: where ρ(u, w) is a similarity measure between the voxels u and w. It is dened on the d-dimensional attribute vectors A i (•) characterizing each voxel (i = 1, 2

for image I i ), or:

Hence, the smaller their attribute vector Euclidean distance, the higher their similarity will be. The central neighborhood CN (•) and peripheral neighborhood P N (•) are dened as concentric rings around the voxel of interest, and their radii are determined according to the scales of the descriptors A i [Ou & Davatzikos 2009]. More specically, we construct the attribute vector A i (•) by incorporating the multi-scale and multi-orientation Gabor attributes. We opted for Gabor attributes to characterize each voxel because dierent application, including registration (e.g., [START_REF] Liu | [END_REF] [ Ou & Davatzikos 2009]), proved their success in practice. Moreover, their ability to capture multi-scale and multi-orientation cues provides them with a great potential to describe each voxel more distinctively than the intensity attribute or other texture features [START_REF] Kadir | Timor Kadir and Michael Brady. Saliency, Scale and Image Description[END_REF]. tween the images. Therefore, we keep the K top-ranked pairs and seek the globally optimal ones using a Markov-random-eld (MRF)-based optimization. This step is explained in the next sub-section.

Finding Globally Optimal Landmark Pairs by MRF Formulation

Let us denote all the J × K landmark pairs detected in the previous step as a set P = {(p k j , q k j )|j = 1, 2, . . . , J; k = 1, 2, . . . , K}. One should note that K

is not necessarily the same number for dierent regions. In practice, we keep K = 10 top-ranked candidate pairs in each region.

We aim in this section to select one pair (out of K pairs) from each region, i.e., P = {(p l j j , q l j j )|j = 1, 2, . . . ; l j ∈ {1, 2, . . . , K}}, such that the obtained set of pairs is globally optimal, in that they altogether maintain high mutualsaliency as well as the smoothness of the resultant deformation eld.

We can formulate this task as a Markov-random-eld (MRF)-based labeling problem. We seek a set of optimal labels l = (l 1 , . . . , l J ), where l j ∈ {1, 2, . . . , K} is the label (or index) of the globally optimal pair out of K candidate pairs in the j th region. Therefore, we construct a graph G = (V, E),

where each node in V represents a region and each edge in E expresses the constraint of a local geometric compatibility between the landmark pairs.

In this formulation, our goal can be achieved by minimizing the following labeling energy on the graph, l = arg min l E(l 1 , . . . , l J ),

where the energy E(l 1 , . . . , l J ) consists of two terms, E(l 1 , . . . , l J ) = E data (l 1 , . . . , l J ) + α E reg (l 1 , . . . , l J ).

(A.4)

Those two terms express the two criteria for selecting the globally optimal landmark pairs. The rst term, E data , is related to the quality of the landmark pair, as dened in the pair extraction step, or:

where σ is a scaling factor, estimated as the standard deviation of µ(., .)×ρ(., .)

values of all the candidate pairs.

The second term, E reg is a regularization term that favors the smoothness of the resultant deformation eld. The regularization imposes constraints on both the spatial positions and the displacement directions on the selected landmark pairs. Specically, spatial position constraint encourages those selected pairs to scatter in the image space other than being close to each other. Displacement direction constraint encourages displacement vectors on adjacent pairs to be consistently oriented, therefore avoiding self-intersection of the resultant deformation eld. Those two constraints are similar to the distortion terms in [START_REF] Leordeanu | A Spectral Technique for Correspondence Problems Using Pairwise Constraints[END_REF] and [START_REF] Torresani | [END_REF]]. Mathematically, the regularization term is expressed as:

One should note that such a framework is not invariant to scale changes.

However, since we have assumed that a rigid registration step was performed beforehand, this is not a main concern.

The dened labeling problem is NP-hard in general. However, recent approximate solvers perform well in practice. The total energy in Eqn. (A.3) is minimized using the sequential tree-reweighted message passing (TRW-S)

algorithm [Kolmogorov 2006].

are based on Markov Random Fields (MRFs).

[ Kohli et al. 2008] introduced an algorithm that performs human pose estimation and segmentation.

Their method estimates in each frame the pose parameters of the stickman shape model, which is subsequently used in a binary graph cut segmentation as a distance map prior. An alternative approach is presented in [Malcolm et al. 2007b], where the solutions in the previous frames are used to predict the poses of the tracked objects. Then, binary shape maps that are centered at the predicted positions are used to enforce shape prior in a multi-label graph cut. Another method of multiple object tracking was recently proposed in [Wang et al. 2009]. The authors dene a generative image model that naturally handles occlusions in the MRF framework. Minimizing the objective energy allows them to simultaneously segment and track objects, and determine their depths. In another recent work, [START_REF] Bugeau | Detection and segmentation of moving objects in complex scenes[END_REF] address the problem of tracking in dynamic scenes with moving camera. The authors proposed an algorithm that is threefold. They determine rst a set of points which motion is dierent from the camera's (the camera motion is assumed to be the predominant one in the scene). A descriptor vector that combines multiple cues such as motion and photometric features is then associated with each extracted point. In the second step, these points are clustered based on their descriptors using a variable-bandwidth mean-shift algorithm.

Then, the joint segmentation and tracking are achieved by performing a graph cut on the pixel level, where the cluster information is used as a constraint.

In this section, we present an extension of our segmentation framework to the temporal domain. We describe a tracking algorithm that is based on image boundary information, and that benets from static and dynamic shape priors. The static shape prior refers to the normalized chord length shape model introduced in section 2.2. Hence, the spatial consistency is encoded through probability densities on the relative positions of pairs of control points. Similarly, we encode temporal priors through relative deformations of dierent control points in time once motion has been implicitly accounted for.

Tracking is then reformulated as an energy minimization problem, where the positions of the control points are deformed in the current frame in order to solve a MAP-MRF problem.

The remainder of this section is organized as follows: rst, we dene the points are guaranteed, one can learn a probability distribution p d ij (t, τ ) for each of the "temporal" distances d ij (t, τ ), and for dierent dynamic orders τ . We consider such interactions not only for the the same landmark (simple case d ii (t, τ )), but for the relative dynamics of the dierent landmarks.

Then, given the position of the object at the previous frames (t -τ, . . . , t -1), one can dene dynamic constraints that account for the correlation between the landmark positions at the current frame and their motion history. More formally, the dynamic shape energy term can be written as:

. Hence the dynamic shape prior is expressed as a sum of unary potentials. Now that we have dened the static and dynamic prior model, the next task consists in dening relationships between the graph and the image sequence (I 1 , . . . , I N ), which enables us to infer the optimal positions of the control points.

Weak Edge Support

We dene the data-term in an image I t of the sequence using its edges. Independently from the previous frames, the optimal positions of the control points should superimpose the boundary of the object to its edges in the image. We suppose that the boundary of the object can be simply approximated by connecting the successive control points in 2D (the 3D counterpart is a mesh of the control points, which would lead to high order MRF terms).

Let M t be the distance map to the edges of the considered image (we computed in practice the chamfer distance to the canny edges). Then, since the object contour can be approximated by the landmark segments, its distance to the edges can be computed using the same approximation. Supposing that the control points are indexed according to their order on the boundary, this

A.2.3 Experimental Validation

To test the applicability of our model and algorithm, we considered a commonly-used example in tracking: walking people. The problem of tracking walking people provides a deformable object with interesting dynamics. We used in our experiments video sequences from the Georgia Institute of Technology database (http://www.cc.gatech.edu/cpl/projects/ hid/Description.html).

We rst selected a total of 455 frames of walking persons from the database with dierent gaits. We labeled manually these frames placing the landmarks at the corresponding positions. Then we used these labeled images to learn a static prior and a dynamic prior. Next, for the testing, we applied our trained priors to three video sequences. The model points were initialized close to the walking target. For the qualitative evaluation, we compare the results that we obtain using the static prior, and then by using the static and the dynamic prior. From our experiment we observe that the results with dynamic priors are more robust with respect to the noisy edges, and have a better quality than those from the static prior tests. Proof of Theorem 2.1. Suppose that there exists N points (x 1 , . . . , x N ) in R n such that d ij = x i -x j 2 for 1 ≤ i, j ≤ N , and that span(x 1 , . . . , x N ) = R n . Without loss of generality, we may assume that x 1 = 0. Upon writing x i -x j = (x i -x 1 ) + (x 1 -x j ), we see that:

Subsequently, for any y ∈ R N -1 , we have:

meaning that A is positive semidenite. Moreover, if we let V to be the n × (N -1) matrix whose i th column is the vector x i+1 (where 1 ≤ i ≤ N -1), then we see that A is the Gram matrix of V , or A = V T V . Since span(x 1 , . . . , x N ) = R n , we have rank(A) = n.

Conversely, suppose that A is positive semidenite and has rank n ≤ N -1.

Then we can write A = U ΛU T , where U is a (N -1) × (N -1) orthogonal matrix and Λ = diag(λ 1 , . . . , λ n ), where λ

are the eigenvalues of A. Now, let (e 1 , . . . , e N -1) be the standard basis vectors of R N -1 , and set

, where P n is the orthogonal projection of R N onto the rst n coordinates. Then for 1 ≤ i, j ≤ N -1 we have:

Thus, (0, x 2 , . . . , x N ) denes the desired realization, which completes the proof.

Appendix C

Spanning k-tree of a k -tree Proof of Lemma 2.1. Without loss of generality, we may assume that k > k.

Let T be a k-tree with t vertices v 1 , . . . , v t , and let v be the last vertex in a

given SEO of T . The proof of the lemma is given by induction on t.

If t = k , then T is a clique. Let T be the k-tree that is built as follows (according to denition 2.2): start from the clique {v 1 , . . . , v k }, and for k+1 ≤ i ≤ t, connect v i to {v i-k , . . . , v i-1 }. Then T solves the problem and (P) holds clearly.

If t > k , we apply the induction hypothesis to T -v. Then there exists a spanning k-tree T of T -v, such that every k -clique of T -v contains a k-clique of T . In particular, the neighborhood of v in T -v is a k -clique Q , and so it contains a k-clique Q of T . Let T be the graph obtained by adding to T the vertex v and all the edges between v and Q. Clearly T is a spanning k-tree of T . Moreover, T satises (P). Indeed, we already know that every k -clique of T that does not contain v contains a k-clique of T , and thus of T . Furthermore, every k -clique C of T that contains v excludes exactly one vertex u in the neighborhood of v in T (because v is simplicial and hence has exactly k neighbors in T ). Thus C = Q + v -u.If u is not a vertex of Q then Q ⊂ (Q -u) ⊂ C , and (P) is proved. Else, the remaining k -1 vertices of Q together with v form a k-clique of T contained in C .