1 introduction 1 1.1 Context and Motivations 1 1.2 Problem definition 5 1.3 Contributions 6 1.4 Thesis Plan 9 2 state of the art in segmentation 11 2.1 Surface-based models 12 2.1.1 Active Contours 12 Active Contours and edge information 13 Active Contours and region information 15 2.1.2 Surface models and shape priors 16 Single reference shape models 16 Deformable templates 20 Superquadrics, hyperquadrics 20 Landmark-based models 21 2.1.3 Implicit deformable models: level sets 24 2.1.4 Level sets and shape priors 27 2.1.5 Graph-based surface models 28 2.2 Region-based methods 30 2.2.1 Continuous methods 30 2.2.2 Graph-based methods 32 2.2.3 Segmentation with Graph Cuts 32 2.2.4 Shape priors and graph-cuts 36 2.2.5 Segmentation with Random Walks 37 2.3 Prior art limitations 37 3 automatic seed placement 39 3.1 Introduction 39 3.2 Related work 40 3.3 Prerequisites 41 3.4 Unlabeled seeds sampling 41 3.5 MRF formulation 42 3.5.1 Energy form 43 3.5.2 Unary potential 44 Geodesic distance potential and graph connectivity 44 Relative orientation potential 46 3.6 Experimental validation 47 3.6.1 MRF Optimization 47 ix x Contents 3.6.2 Experimental Results 47 3.7 Conclusion 48 4 random walks segmentation and prior knowledge 53 4.1 The Random Walks segmentation algorithm 53 4.1.1 Preliminaries 53 4.1.2 Random Walker analogy 54 4.1.3 Probabilistic explanation 54 4.1.4 Transition probabilities 56 4.1.5 Random Walks objective functional 57 4.1.6 Segmentation algorithm 59 4.1.7 Related work 59 4.2 Random Walks with prior knowledge 62 4.2.1 Random Walks segmentation with appearance prior 63 4.2.2 Probabilistic framework 64 4.3 Random Walks segmentation with shape prior 67 4.3.1 Constant prior weight 67 4.3.2 Prior-weighted models 68 4.3.3 Confidence map 69 4.3.4 Experimental validation 70 4.4 PCA model 79 4.4.1 Additional notations 79 4.4.2 Shape space 81 4.4.3 PCA model 82 4.4.4 Experimental validation 83 4.5 Conclusion 84 5 learning parameters for random walks-based segmentation 87 5.1 Related Work 88 5.2 Preliminaries 89 5.2.1 Notations 89 5.2.2 Parameters and Feature Vectors 90 5.2.3 Loss Function 91 5.3 Parameter Estimation using Latent SVM 92 5.3.1 Upper Bound on the Risk 92 5.3.2 Formulation with Latent SVM 93 5.3.3 Annotation Consistent Inference 95 5.3.4 Parameter Update 97 5.4 Experimental validation 98 6 summary and conclusion 105 6.1 Material 105 6.2 Summary of the contributions 105 Contents xi 6.3 Future directions 106 a mri basics 109 a.1 Basics in Nuclear Magnetic Resonance 109 a.2 NMR Imaging 111 b mathematical appendices 115 b.1 Positive definiteness of L U 115 b.2 Random Walks and probability distributions 116 b.3 Random Walks, shape prior and probability distributions

A B S T R A C T Segmentation of magnetic resonance images (MRI) of skeletal striated muscles is of crucial interest when studying myopathies. Diseases understanding, therapeutic followups of patients, etc. rely on discriminating the muscles in MRI anatomical images. However, delineating the muscle contours manually is an extremely long and tedious task, and thus often a bottleneck in clinical research. Typical automatic segmentation methods rely on finding discriminative visual properties between objects of interest, accurate contour detection or clinically interesting anatomical points. Skeletal muscles show none of these features in MRI, making automatic segmentation a challenging problem. In spite of recent advances on segmentation methods, their application in clinical settings is diffcult, and most of the times, manual segmentation and correction is still the only option.

In this thesis, we propose several approaches for segmenting skeletal muscles automatically in MRI, all related to the popular graph-based Random Walker (RW) segmentation algorithm. The strength of the RW method relies on its robustness in the case of incomplete contours and its efficient optimization. Originally, the RW algorithm was developed for interactive segmentation: the user had to pre-segment small regions of the image -called seeds -before running the algorithm which would then complete the segmentation. Our first contribution is a method for automatically generating and labeling all the appropriate seeds, based on a Markov Random Fields formulation integrating prior knowledge of the relative positions, and prior detection of contours between pairs of seeds. A second contribution amounts to incorporating prior knowledge of the shape directly into the RW framework. Such formulation retains the probabilistic interpretation of the RW algorithm and thus allows to compute the segmentation simply by solving a large sparse linear system, like in the original method. In a third contribution, we propose to develop a learning framework to estimate the optimal set of parameters for balancing the contrast term of the RW algorithm and the different existing prior models. The main challenge we face is that the training samples are not fully supervised. Specically, they provide a hard segmentation of the medical images, instead of the optimal probabilistic segmentation, which corresponds to the desired output of the RW algorithm. We overcome this challenge by treating the optimal probabilistic segmentation as a latent variable. This allows us to employ the latent Support Vector Machine (latent SVM) formulation for parameter estimation. v R É S U M É La segmentation d'images anatomiques de muscles striés squelettiques acquises par résonance magnétique nucléaire (IRM) présente un grand intérêt pour l'étude des myopathies. Elle est souvent un préalable nécessaire pour l'étude les mécanismes d'une maladie, ou pour le suivi thérapeutique des patients. Cependant, le détourage manuel des muscles est un travail long et fastidieux, au point de freiner les recherches cliniques qui en dépendent. Il est donc nécessaire d'automatiser cette étape. Les méthodes de segmentation automatique se basent en général sur les différences d'aspect visuel des objets à séparer et sur une détection précise des contours ou de points de repère anatomiques pertinents. L'IRM du muscle ne permettant aucune de ces approches, la segmentation automatique représente un défi de taille pour les chercheurs.

Dans ce rapport de thèse, nous présentons plusieurs méthodes de segmentation d'images de muscles, toutes en rapport avec l'algorithme dit du marcheur aléatoire (MA). L'algorithme du MA, qui utilise une représentation en graphe de l'image, est connu pour être robuste dans les cas où les contours des objets sont manquants ou incomplets et pour son optimisation numérique efficace. Dans sa version initiale, l'utilisateur doit d'abord segmenter de petites portions de chaque région de l'image, appelées graines, avant de lancer l'algorithme pour compléter la segmentation. Notre première contribution au domaine est un algorithme permettant de générer et d'étiqueter automatiquement toutes les graines nécessaires à la segmentation. Cette approche utilise une formulation en champs aléatoires de Markov, intégrant une connaissance à priori de l'anatomie et une détection préalable des contours entre des paires de graines. Une deuxième contribution vise à incorporer directement la connaissance à priori de la forme des muscles à la méthode du MA. Cette approche conserve l'interprétation probabiliste de l'algorithme original, ce qui permet de générer une segmentation simplement en résolvant un numériquement un grand système linéaire creux. Nous proposons comme dernière contribution un cadre d'apprentissage pour l'estimation du jeu de paramètres optimaux régulant l'influence du terme de contraste de l'algorithme du MA ainsi que des différents modèles de connaissance à priori. La principale diffculté est que les données d'apprentissage ne sont pas entièrement supervisées. En effet, l'utilisateur ne peut fournir qu'une segmentation déterministe de l'image, et non une segmentation probabiliste comme en produit l'algorithme du MA. Cela nous amène à faire de la segmentation probabiliste optimale une variable latente, et ainsi à formuler le problème d'estimation sous forme d'une machine à vecteurs de support latents (latent SVM).

Magnetic resonance Imaging

Magnetic Resonance Imaging (MRI) is a medical imaging technique providing 2D or 3D views of the internal structures of the body. MRI makes use of the principle of Nuclear Magnetic Resonance (NMR) which allows to measure some properties of the nucleus of some specific atoms in presence of a powerful magnetic field. MRI is noninvasive and implies non-ionizing radiations (e.g. unlike X-rays) allowing long and repeated acquisitions with no hazard posed for the health of patients.

Neuromuscular Diseases

Neuromuscular diseases may involve disorders in the motor neurons, in cells of the spinal cord, in peripheral nerves or in the muscles themselves. There exist both hereditary -dominant or recessive -and acquired -inflammatory, toxic, etc -disorders. Among the features assessed by clinicians to characterize neuromuscular diseases one finds: time of onset, chronicity, presence or absence of pain, muscle weakness and distribution of the involvement -proximal, distal, facial, symmetric or asymmetric, hypertrophical or atrophical, etc.

Clinical MRI is a non-invasive and powerful tool to characterize skeletal muscle involvement in neuromuscular disorders. This imaging technique has progressively established itself as a reference technique for the initial assessment of pathologies and for the evaluation of therapeutic intervention. In the NMR Laboratory of the Institute of Myology in Paris, several research projects focus on developing biomarkers for skeletal muscle disorders based on MR images: relative image intensity, intensity histogram, signal heterogeneity and more sophisticated methods of texture analysis, geometrical analysis (distance, area, volume, shape), etc. Moreover, some NMR parameters are intrinsic to the different types of tissues and can be quantified: T 1 , T 2 , proton density, T ⋆ 2 , T 1ρ and T 2ρ , D (diffusion constant), etc; We refer the reader to specialized publications -such as [START_REF] Fleckenstein | Muscle imaging in health and disease[END_REF] -for further information on the topic of muscle MR imaging; a brief introduction is also provided in the appendices of this document (see appendix A). In Fig. 1.1.1, we display MRI axial cross-sections of healthy subject and patients, showing the varied visual patterns of some myopathies. By displaying the anatomical structures -bones, vessels, fat and muscle tissues, etc. (see Figs. 1.1.2 and 1.1.1) -anatomical MR imaging allows to visually discriminate between the muscles and other organs, and thus to perform segmentation.

Image Segmentation

"In computer vision, segmentation is the process of partitioning a digital image into multiple segments (sets of pixels). The goal of segmentation is to simplify and/or change the representation of an image into something that is more meaningful and easier to analyze. Image segmentation is typically used to locate objects and boundaries (lines, curves, etc.) in images. More precisely, image segmentation is the process of assigning a label to every pixel in an image such that pixels with the same label share certain visual characteristics" 1 . Fig. 1.1.3 shows an example of image segmentation where a person in the picture is separated from the background.

Segmentation methods are needed to facilitate the identification of a muscle or of a muscle group. This step is particularly crucial when muscles have to be investigated selectively: for instance, when one wants to determine the NMR properties of each muscle group and study the differences between them. It is also necessary to measure the volume of the muscles and monitor their evolution.

Today, segmentation of muscle images is performed manually, by drawing regions of interest on the images. Such task requires enormous amounts of time, and is an extremely tedious and tiring process. Due to the advances in acquisition techniques and (center) same cross-section with all muscles manually segmented; (right) 3D segmentation of the same volume.

computer storage facility, it is common to acquire MRI 3D volumes composed of more than 50 slices. As a result, manually segmenting a volume of a thigh -which comprises more than 15 muscles -represents more than four hours of work. Moreover, during clinical research protocols, patients are often scanned with several sequences and very large amount of images can be generated. As acquiring images becomes easier and faster for large volumes, segmentation of muscles images is becoming more and more a bottleneck in the chain of analysis. Besides, manually processing these data is not only time consuming, it is also highly observer-dependent: different experts produce significantly different segmentations, especially in places where the boundaries of the organs are not clearly visible. Fig. 1.1.4, shows an example of manual segmentation of the muscles in anatomical MRI.

For these reasons, it is highly desirable to possess automatized segmentation tools to process the data more quickly and consistently and thus facilitate their interpretation. This is a crucial objective to sustain the growth of applications of skeletal muscle MRI in the field of neuromuscular disorders.

Whereas many automatized methods have been developed for segmenting images of brains and hearts, very few attempts were made on skeletal muscles. Medical studies on the brain and the heart are more developed due to the existence of widespread diseases associated with theses organs, whereas severe muscle diseases are often genetic and fall in the category of orphan diseases. The specifics of muscle segmentation prevent using directly the methods developed for other organs and require specially adapted methods. Moreover, a large part of the existing publications on muscle seg-mentation focus mainly on the fat-muscle tissue separation rather than on muscle group segmentation itself [START_REF] Salam | Semi-automated segmentation and quantification of adipose tissue in calf and thigh by mri: a preliminary study in patients with monogenic metabolic syndrome[END_REF][START_REF] Barra | Segmentation of fat and muscle from MR images of the thigh by a possibilistic clustering algorithm[END_REF][START_REF] Jp Mattei | Segmentation of fascias, fat and muscle from magnetic resonance images in humans: the dispimag software[END_REF] In this thesis report, we present our work on developing new automatic image segmentation methods adapted to MRI images of skeletal muscles.

problem definition

Any automatized segmentation method relies on (a) finding consistent features of targeted objects -e.g. shape or appearance descriptors -(b) allowing some variability in the target features and (c) defining means of overcoming anticipated difficultiese.g noise in images:

• consistent features of the skeletal muscles in MRI are (cf. Fig. 1.1.4):

fixed relative positions -e.g. muscle a is always top left of muscle b;

interweaved positions, closely entangled -e.g. muscle a is always adjacent to muscle b;

shapes, to some degree -e.g. muscle a can be visually identified based on prior anatomy knowledge;

smooth boundaries, to some degree, especially in the direction parallel to the limbs -i.e. muscles can be visually separated from each other based on visible boundary structures.

muscle/background appearances -i.e. muscles can be visually distinguished from other tissues based on their appearance.

• variability is found in (cf. figure 1.2.1):

shapes, with large inter-subject variations of the muscles size and shapes.

• anticipated difficulties are (cf. figure 1.2.1):

non-discriminative inter-muscle appearance -e.g. one cannot distinguish muscle a from muscle b based on intensity/texture differences;

missing boundaries -e.g. muscle a is sometimes not separated from muscle b by any visible structure;

spurious contours -e.g. muscle a can present internal structures visually similar to boundaries.

Considering these characteristics, our aim is to design semi-automatic or automatic segmentation methods with the following features:

• multi-object; (left) Cross-section of a thigh from a 37 year-old female subject; (right) same height cross-section from a 29 year-old male subject. Notice the differences in shape and proportion of the muscles between the two subjects.

• relying solely on:

detectable boundaries;

prior knowledge of the muscle shapes and topology;

muscle/background intensity differences;

• robust to:

missing boundaries;

missing contours;

shape variations.

contributions

In this thesis report, we present several contributions related to the general Random Walks (RW) for image segmentation algorithm [START_REF] Grady | Random walks for image segmentation[END_REF]. This algorithm has the advantage of relying solely on boundaries while being robust to partial contours, and provides a simple multi-object segmentation framework with no topological issues (that is, no intersection issues).

Originally, the RW algorithm was developed for interactive segmentation: the user had to pre-segment small regions -called seeds -of the image to be segmented be-fore running the algorithm which would then complete the segmentation. Our first contribution is a method for generating and labeling appropriate seeds automatically. This approach exploits the fact that relative positions of the muscles are consistent among individuals, even though the large inter-subject variations of skeletal muscles prevent using a deformable model based on a shape template. Thus, it is reasonable to develop a segmentation method whose prior information is not based on the position of the boundaries and on the shape of the objects, but rather on the approximate relative position of the objects and the existence of boundaries between them. Seeds are compatible with this idea: they are generally placed well inside the target object so that the actual boundaries are determined by the Random Walks algorithm. Thus, our method aims at placing the seeds automatically, based on prior knowledge of the relative positions and on the detection of contours between pairs of seeds. We have the following stages: 1) generate unlabeled seeds within homogeneous regions of the target image; 2) build a graph over the sampled seeds, viewing each seed as a node and adding edges between nearby ones; 3) assign a label to each seed automatically by minimizing an functional defined on the graph nodes and edges; 4) the image is segmented with the RW method using the newly generated seeds. During the labeling stage 3), the node-wise potentials depend on prior-knowledge over label assignments based on a comparison of the target image with a registered reference annotated image. The pair-wise potentials depend on the likelihood of the edge orientation and on the detection of boundaries between the seeds, where both quantities can be obtained by statistical learning over the reference image. This labeling problem is solved using an efficient discrete energy minimization method based on a belief propagation algorithm [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF]. Due to the non-specific (non sub-modular) formulation of the pair-wise potentials, the labeling cost functional can only be minimized approximately; however, we present experiments on real clinical data showing that the approximate solution is better than not using pair-wise potentials. This method was published in the proceedings of the 2012 ISBI conference (Baudin et al., 2012c).

This first method has the inconvenience of being two-stage: if the seed labeling fails, then the segmentation will irremediably fail in the same area. Fortunately, it was shown that it is possible to completely bypass the seed placement stage by incorporating prior-knowledge on the intensity to the RW formulation [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF]. However, intensity is not a reliable feature in the case of muscles separation. Our second contribution consists in proposing several models mixing Random Walks and prior knowledge of the shape. The Random Walks segmentation algorithm is a graphbased method working at the pixel level: each pixel is viewed as a node and edges are defined between each adjacent pixels; the weights of the edges are set according to the strength of the contrast between the connected pixels. The aim is to assign a label (e.g. the index of a muscle) to each pixel. Segmentation is achieved by minimizing a quadratic objective functional based on the combinatorial Laplacian matrix of the image, which amounts to solving a large sparse linear system. In its original version, a unique solution is found only if there exists at least one seed per label. In our approach, a prior term depending on the squared distance of the segmentation to a statistical atlas is added to the RW functional. We point out that such formulation does not break the probabilistic explanation of the RW method, as the added prior term can be seen as adding implicit seed nodes. Noting that each pixel's prior can be weighted independently from the others, we propose several weighting schemes, either based on statistics over the training set or on information from the target image itself. This method was published in the proceedings of the 2012 MICCAI conference (Baudin et al., 2012a). In a second approach, published in the proceedings of the 2012 BMVC conference (Baudin et al., 2012b), we propose to further exploit our training set by building a statistical shape space -a low dimensional space of valid solutions -through Principal Component Analysis. This is motivated by the need to increase the flexibility of our model so that it can adapt to more diverse morphologies. The segmentation is no longer anchored to a unique average atlas like in the previous approach, but is constrained to remain "close" to the shape space, which models many possible valid segmentations. All approaches are validated on real clinical data.

While it is interesting to determine which model gives better results on the muscle data, it is also desirable to combine them together -e.g. using different weighting schemes plus the intensity prior. Moreover, we could also want to use several weighting functions for setting the graph edges in order to capture varied properties of the images simultaneously. However, achieving such generality requires to set balancing weights for a large number of parameters, which would certainly not be optimal if done manually. As a third contribution, we propose to develop a learning framework to determine the optimal set of parameters for balancing the contrast term of the RW algorithm and the different prior models. Our approach is guided by noticing the RW formulation bears important similarities with the more wide-spread pixel-based Markov Random Fields (MRF) [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF] segmentation methods. In recent work, machine learning methods -such as Support Vector Machines (SVM)were used to fit the parameters of the segmentation method to the data type of interest [START_REF] Szummer | Learning CRFs using graph cuts[END_REF]. We propose to adapt the MRF/SVM learning scheme to the RW case. We explain why a direct transposition in the RW case of the previous formulation cannot solve the problem, due to having only a "hard" -i.e. binaryannotation vector for our training data whereas a "soft" -i.e. probabilistic -one is necessary. We overcome this issue by introducing an implicit soft annotation vector, which amounts to solving a latent SVM problem. Due to the large number of variables rendering certain stages of the problem too difficult to solve exactly, we propose approximation schemes in the "Annotation Consistent Inference" (ACI) stage and in the "Loss Augmented Inference" (LAI) stage. Using the same clinical data as previously described, we present preliminary results which demonstrate the potential of this learning approach for the prior-enhanced Random Walks segmentation method.

thesis plan

This report is organized as follows: in chapter 2 we present a detailed state-of-theart of segmentation methods in medical imaging; the next three chapters present our contributions to the field: automatic seed placement in chapter 3, introducing shape priors to the Random Walks framework in chapter 4 and developing a SVM-based learning scheme for improving the latter models in chapter 5; the last part concludes this report.

To address the problem of image segmentation, a wide variety of techniques was presented and one can find a very abundant literature. However, few research works have dealt with muscle image segmentation to the best of our knowledge. In this chapter, we review the existing segmentation methods applied to medical images and their optimization procedures, with a special attention paid on those which were designed to address muscle segmentation.

Almost all segmentation methods have in common is how they rely on a cost functional, or energy functional, E (M, I) which measures how well a proposed segmentation model M explains the target image I. In other words, an energy functional is designed to measure the adequacy of the segmentation model with respect to a combination of criteria, such as alignment of the model surface with image edges, intensity homogeneity within model regions, plausibility of model shape, etc. The better the model M fits the image I, the lower is E (M, I). Depending on the chosen approach, E (•, I) can be convex and have only one global minimum, or non-convex and have many local minima. Determining a model and its associated functional is the most crucial step in developing a new segmentation method. Indeed, the desired segmentation has to correspond to the global minimum, or to an "accessible" local-minimum of E (•, I), otherwise one cannot expect to retrieve such solution relying on E.

When presented with a new image I to segment, the model M is fit to the image via an optimization procedure, which consists on minimizing E (•, I) with respect to the model parameters. Naturally, this optimization step is also crucial in determining the quality of the obtained solution. Even when the global minimum of E (•, I) cor- responds to the desired solution, it sometimes cannot be reached, either because the optimization procedure is mislead to an incorrect result -e.g. a local optimization can be "stuck" in non-desired local minimum -or because it would take too long to reach it -e.g. a NP-hard problem.

In this chapter, we will distinguish between two types of approaches. In a first part, we give an overview of vast family of surface-based models -i.e. models which rely on mathematical representations of the surface of the target objects. Among the most widely used models, Deformable Models, either explicit or implicit, are surface models which are progressively fit to the image in a local optimization process. An optimization process is said local when the solution is found by applying successive local changes from an initial state, a canonical example being the gradient descent. Among them, landmark-based models are surface models based on remarkable, anatomical or statistically consistent image points, which allows introducing statistical shape priors for constraining the surface position. Finally, unlike Deformable Models, which mostly rely on non-optimal local-optimization methods, graph-based surface models can reach a global minimum of the cost functional thanks to efficient discrete optimization algorithms. This advantage is paid by a certain approximation of the solution due to the discretization of the solution space.

In a second part, we present region-based models -i.e. which parametrize every point of the image domain -devoid of intersection issues in the case of multi-object segmentation, but whose topology and shapes are more difficult to constrain. Among these, some use a typical approach based on continuous local optimization, while others rely on a graph formalism: Graph-Cuts and Random Walks. Graph-based models provide convenient and efficient frameworks for multi-object segmentation, for dealing with weak edges, and for accepting simple user interaction in the form of seeds. However there are still few examples of such models incorporating prior knowledge of shapes.

surface-based models

Although they may have very diverse mathematical formulation, or rely on very different types of information, surface-based models all have in common a parametrized model of the surface of the objects of interest. Given a particular state of the model, the image segmentation is determined by the location of the model surface: any point in the image domain is either on one side of the surface, or the other. If the surface is closed, then any point is either inside the volume delimited by this surface -thus assimilated to this object's class, or outside -thus assimilated to another class or to the background class.

The value of the energy functional of a surface model thus depends on the location of the surface with respect to the segmented image. What is measured by the functional can either be related to image boundaries -i.e "is the model's surface matching the edges of the image ?" -or to image regions -i.e. "do the properties of the regions delimited by the surface model match what I know of them a priori ?".

Active Contours

Active Contours is a general term which designates surface-based models whose functional minimization procedure is a suboptimal strategy consisting in evolving an initial surface (or contour in 2D) in the direction of the steepest descent of the functional. In general, the cost functional which is used is not convex and the adopted minimization approaches only reach local minima, highly dependent on the initialization state.

The first and most simple instances of active contours rely only on low-level image features, such as edges, region intensity, texture, as well as low-level shape features, such as contour length, object's surface (object's volume in 3D), contour smoothness, etc. However, a number of factors such as noise, occlusion, complexity of shapes, etc. as well as sensibility to initial conditions often limit the performance of purely low-level methods. For these reasons, higher-level constraints on shapes were later introduced, called shape priors, based on previously seen similar objects.

Let us first consider a boundary object C embedded in the domain space Ω of the image I, such as a planar curve (C ⊂ R 2 ) or a surface (C ⊂ R 3 ) to be deformed from an initial position in order to minimize an energy functional E (C). The energy E (•) is specific to the model, as it is determined by the kind of image information one wants to exploit; in the next sections, we present several formulations for E (•), depending on whether it uses edge information or region information.

The energy functional E can often be minimized using finite elements methods for the variational form, by making C depend on time:

δC δt (s, t) = -∇ C E (C (s, t) , t) . (2.1.1)
This expression is the base evolution equation for most active contours models. Such formulation has the visual effect of making the the surface of the boundary object deform with time, hence the names "active contours" and "deformable models". The literature on active contours is extremely large and, for this reason, we only give a brief and incomplete survey of it. For a deeper review of active contour methods, we indicate to the reader several review articles and books: [START_REF] Mcinerney | Deformable models in medical image analysis: a survey[END_REF]; [START_REF] Montagnat | A review of deformable surfaces: topology, geometry and deformation[END_REF]; [START_REF] Cremers | A Review of Statistical Approaches to Level Set Segmentation: Integrating Color, Texture, Motion and Shape[END_REF]; [START_REF] Paragios | Handbook of mathematical models in computer vision[END_REF]. Active contours can be driven either by the edges of the image or by the properties of the image regions -or by both simultaneously.

Active Contours and edge information

Active Contours for image segmentations based on edge information, also called Snakes, were introduced by [START_REF] Kass | Snakes: Active contour models[END_REF], along with the typical energy functional:

E (C) = ˆ1 0 α 2 |C s | 2 + β 2 |C ss | 2 + E ext (C) ds, (2.1.2)
where C is a 2D parametric curve C : [0, 1] → Ω, C s and C ss denote the first and second derivative of C, α and β are weighting parameters, respectively for tension (C s ) and rigidity (C s ). These first two terms can be viewed as internal energy terms, as they measure properties of the contour itself. E ext is a function based on an edge map of the image I, such that it is minimal on image edges. This term is an external energy term, measuring the adequacy of the contour with the image. A canonical example for E ext is the image squared gradient norm:

E ext (x) = -|∇I (x)| 2 .
(2.1.3) Deriving equation 2.1.2 by following the scheme of 2.1.1 gives:

δC δt = αC ss -βC ssss -∇E ext . (2.1.4)
Such formulation can be extended to 3D straightforwardly [START_REF] Terzopoulos | Constraints on deformable models: Recovering 3D shape and nonrigid motion[END_REF].

Many variations on Snakes were developed, e.g. to decrease the dependance on the initial conditions, the balloon model [START_REF] Laurent | On active contour models and balloons[END_REF] adds a force term to equation 2.1.4:

δC δt = αC ss -βC ssss -∇E ext + νn, (2.1.5)
where ν > 0 sets the strength of the balloon force, and n is the outward normal to the contour. Such expression amounts to maximizing the area of C while keeping the boundaries smooth. Other examples include: adopting a multi-resolution framework using Fourier descriptors for parametrizing the contour [START_REF] Leroy | Multi-resolution algorithms for active contour models[END_REF]; replacing the finite difference numerical optimization approach by a more robust finite element method [START_REF] Cohen | Finite Element Methods for Active Contour Models and Balloons for 2D and 3D Images. Pattern analysis and machine intelligence[END_REF]; using the gradient vector flow in the external energy term in order to extend the capture range of the initial state [START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF].

The mathematical parametrization is also crucial to determine the extent of the solution space; for instance, initial deformable models were based on standard triangular or quadrangular mesh parametrization [START_REF] Terzopoulos | Constraints on deformable models: Recovering 3D shape and nonrigid motion[END_REF][START_REF] Nastar | Fast segmentation, tracking, and analysis of deformable objects[END_REF]; however, this type of parametrization necessitates frequent reparametrization and makes topology changes -that is, changes in the object's Euler characteristic or in the number of objects -difficult. Types of parametrization allowing topology changes have been proposed, such as particle systems [START_REF] Szeliski | Surface modeling with oriented particle systems[END_REF][START_REF] Lombardo | Oriented particles: A tool for shape memory objects modelling[END_REF] -which model the surface with elementary objects whose position is ruled by local forces, simulating the evolution of physical particles in a Newtonian space. The popular Simplex meshes [START_REF] Delingette | General object reconstruction based on simplex meshes[END_REF] -which are topologically dual to triangular meshes -also allow some topology changes and benefit from convenient reparametrization schemes. More recently, Delaunay deformable models have been proposed [START_REF] Pons | Delaunay deformable models: Topology-adaptive meshes based on the restricted delaunay triangulation[END_REF], which model the surface from a Delaunay "tetrahedralization" (that is, the extension of triangulation in 3D) of the entire space based on surface points. This model allows simple topology changes and is free from topology issues, such as undesired intersections of objects.

Active Contours and region information

The models presented in the previous section are considered edge-based, as they are driven by the proximity of edges according to the external energy E ext . Such dependency make these models very sensitive to noise, and inadequate to segment objects which cannot be differentiated by detecting edges between them, e.g. in cases of fuzzy of blurred contours, or of highly textured objects. However, other segmentation approaches are based on image regions rather than edges. The classical approach was proposed by [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]; [START_REF] Blake | Visual reconstruction[END_REF], were the image I is approximated with a piecewise smooth function u, via minimizing the functional:

E (u, C) = 1 2 ˆΩ (I -u) 2 dx + λ 2 1 2 ˆΩ/C |∇u| 2 dx + ν |C| , (2.1.6)
with respect to both u and C. Here, C is a finite set of C 1 -arcs. Function u is designed to be a texture-less approximation of image I, with sharp variations only at the boundaries matched by C and smooth variations inside regions defined by C. When minimizing this functional, the first term ensure the resemblance of u with I, the second term penalizes discontinuities of u except on the boundary C and the third term limitates the length of the boundary to segment only important objects. Minimizing this highly non convex functional is notoriously difficult and necessitates approximate methods and complex optimization procedures [START_REF] Blake | Visual reconstruction[END_REF]. However, a simpler version of the Mumford Shah functional 2.1.6 can be obtained in the limit λ → ∞, yielding the cartoon limit (Mumford and Shah, 1989):

E ({u i } , C) = 1 2 ∑ i ˆRi (I -u i ) 2 dx + ν |C| , (2.1.7)
where {u i } is a partition of I into constant regions {R i } determined by C. This functional was modified by Zhu and [START_REF] Song | Region competition: Unifying snakes, region growing, and Bayes/MDL for multiband image segmentation[END_REF], introducing a Gaussian probabilistic model of the constant gray-level value within each regions: (2.1.8) where P (I|α i ) is the posterior probability of obtaining the intensity values in I inside region R i given distribution parameter α i . A hybrid model, called diffusion snake, was also proposed [START_REF] Cremers | Diffusion Snakes: Introducing Statistical Shape Knowledge into the Mumford-Shah Functional[END_REF]:

E ({α i } , C) = 1 2 ∑ i ˆRi -log P (I|α i ) dx + ν |C| ,
E (u, C) = 1 2 ˆΩ (I -u) 2 dx + λ 2 1 2 ˆΩ/C |∇u| 2 dx + ν ˆ|C s | 2 , (2.1.9)
where the external energy is that of the Mumford Shah functional 2.1.6, while the internal elastic constraint on the boundary is that of the Snake 2.1.2. Such model can then be optimized via standard gradient descent methods.

Surface models and shape priors

The regularization factors in the approaches of the previous sections, rendered necessary due to noise and missing parts, all have the effect of minimizing the contour length [START_REF] Bergtholdt | Variational Segmentation With Shape Priors[END_REF]. These purely geometric constraints can be viewed as shape priors, but they are based only on the intuition that natural shapes have smooth boundaries. Such hypothesis is hardly true in practical cases, and consistently yields excessive regularization, e.g. eroded contours and ignored small details. Shape priors based on previously-seen examples, by imposing the boundary to adopt realistic shapes, allow segmenting complex shapes despite noise and occlusions.

Single reference shape models

The most simple shape prior type in deformable models is enforced when the initial shape of the boundary corresponds to a template of the object to be segmented. The regularization terms are sometimes enough to ensure that the shape remains valid as the model deforms. This reference template boundary can either be manually positioned or automatically positioned onto the target image via rigid registration procedures. This reference shape can be a rough approximation (simple mathematical objects, quickly hand-drawn shapes), or a sample shape obtained from a previously segmented image.

An elastic shape model was introduced by [START_REF] Terzopoulos | Elastically deformable models[END_REF], proposing an internal energy term measuring the changes in the curvature and the length between nearby points -via the first and second fundamental forms -of the deformed body. Spring-mass analogies were also proposed [START_REF] Vasilescu | Adaptive meshes and shells: irregular triangulation, discontinuities, and hierarchical subdivision[END_REF][START_REF] Nastar | Fast segmentation, tracking, and analysis of deformable objects[END_REF], in which the surface model -parametrized as a discrete mesh -has its shape constrained by a set of springs connecting the vertices together. The length of the springs at rest determines the equilibrium position(s) and thus the minimum cost shape configuration(s). The stiffness of the springs determines the elastic properties of the model and thus the range of shapes the model can take. Over the years several models simulating soft tissue deformation have been proposed for segmentation, bridging gaps between computer graphics, biomechanical simulation and computer vision.

muscle segmentation using simplex meshes with medial representations Gilles and Magnenat-Thalmann (2010) presented a skeletal muscle segmentation method based on simplex meshes [START_REF] Delingette | Modélisation, Déformation et Reconnaissance d'objets tridimensionnels a l'aide de maillages simplexes[END_REF]. Considering a 3D surface model, simplex meshes are discrete meshes where each vertex has exactly 3 neighbors. Having a constant connectivity allows to simply parametrize the location of one vertex with respect to its neighbors, and thus parametrize deformation of the shape -translation, rotation, scaling -in a local manner. Indeed, the location of a pixel, denoted as x can be expressed as a linear combination of the locations of the three neighbors plus a local elevation term parallel to the local normal (see figure 2.1.2):

x = ε 1 x 1 + ε 2 x 2 + (1 -ε 1 -ε 2 ) x 3 + hn.
As a result, many local measurements -including curvature and cell surface -can be computed efficiently and global energy terms enforcing local constraints come up naturally.

Here, the authors impose local smoothing via curvature averaging, which does not tend to reduce the surface like 1-order operators typically do. Prior knowledge is imposed by constraining the local scale changes on the elevation parameter with respect to a reference shape. Denoting the surface of the triangle formed by the three neighbors of a pixel as S, given the reference shape parameters ε1 , ε2 , h, S , the new location of the considered pixel is expressed as:

x = ε1 x 1 + ε2 x 2 + (1 -ε1 -ε2 ) x 3 + h S/ S 1/β n, (2.1.10)
where β ∈ [2, +∞[ is a parameter which sets the amount of allowed local deformation: with β = 2 this definition is similitude invariant; with β = +∞ this definition is invariant through rigid transformations only. The model is attached to the target image through either gradient norm maximization in the direction of the gradient at the location of the vertices, or maximization of similarities between the reference and the target images at the vertices location.

A medial representation -similar to the M-reps of [START_REF] Stephen M Pizer | Deformable m-reps for 3d medical image segmentation[END_REF] -is combined with the simplex parametrization to exploit the specific tubular shapes of the skeletal muscles. Medial vertices are added to the model, constrained to remain on the medial axis of the tubular objects. This is achieved by connecting the new vertices to the surface vertices through spring-like forces. This constrains the global structure to resemble its initial reference shape, thus acting as a global shape prior (see figure 2.1.2). This medial axis representation also allows efficient collision handling. The model is fit to the image through an iterative process of successive local evolutions. Such model appear to always yield a valid solution, sometimes at the price of an excessive regularization or lack of adaptability to the specifics of the target image.

muscle segmentation using deformable models and shape matching A shape prior for muscle segmentation in 3D images was derived from a computer animation technique, called shape matching, used to efficiently approximate large soft-tissue elastic deformations [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF]. This method was applied to muscle segmentation with some success. In this approach, discrete meshes are used to parametrize the moving surface. Let x 0 be the vector containing the initial position of the control points of the parametric surface. Clustering is performed on x 0 such that each cluster ζ i contains at least a certain number of vertices (set by the user). During segmentation, the evolution of the active surface is performed according to the following iterative procedure:

1. Shift vertices according to the external force: xt = f ext + x t . The external "force" f ext is computed as the maximal gradient search in the gradient direction.

2. Regularize vertex positions: a) Compute rigid registration for each clusters:

T i = arg min ∑ j∈ζ i T i x 0 j -xt j 2 ,
b) Average target position for each vertex:

x t+1 i = 1 |ζ i | ∑ j∈ζ i T j x 0 j .
Single reference prior models are convenient in that they require only one annotated example of the objects of interest. However, when segmenting a class of objects whose shape varies a lot, such approach becomes too constraining and does not allow the model to adopt valid shapes which are too different from the single reference. (1987). (right) The advantages of generalized gradients (bottom) over standard gradient descent (top) using the same deformable model: a realistic shape is kept all along the deformation process [START_REF] Charpiat | Generalized gradients: Priors on minimization flows[END_REF]. In a different approach, a shape prior was derived from a generalization of the fundamental variational equation 2.1.1 [START_REF] Charpiat | Generalized gradients: Priors on minimization flows[END_REF]. Such generalization is based on the observation that the canonical L 2 inner product is always implied when taking the gradient of the energy. The proposed method relies on replacing this canonical inner product by new inner products, designed to enforce spatially coherent motions with respect to the initial template shape. As a consequence, the trajectory of the boundary is optimized so the energy does not have to include regularization terms. This method can be used both with level sets (see section 2.1.4) and standard Active Contours.

Deformable templates

Deformable templates [START_REF] Yuille | Deformable Templates for Face Recognition[END_REF] are combinations of simple geometric shapes requiring very few parameters -e.g. a center location and a radius for a circle. While their simplicity makes them robust and fast to optimize, these models are not suited to segment complex shapes and large deformations.

Superquadrics, hyperquadrics

Superquadrics [START_REF] Terzopoulos | Dynamic 3D models with local and global deformations: Deformable superquadrics[END_REF][START_REF] Bardinet | Tracking and motion analysis of the left ventricle with deformable superquadrics[END_REF] are geometric shapes defined by parametric equations, generalizing simpler objects such as ellipsoids. Like deformable templates, they require very few parameters, and have similar pros and cons: simple and fast to optimize, but the space of available shapes is too small to model complex shapes. Hyperquadrics [START_REF] Cohen | A hybrid hyperquadric model for 2-d and 3-d data fitting[END_REF][START_REF] Hanson | Hyperquadrics: Smoothly deformable shapes with convex polyhedral bounds[END_REF]) also defined by parametric equations, extend superquadrics by allowing a wider range of shapes. A hyperquadric is defined by the points (x, y, z) satisfying the equation:

H (x, y, z) = N ∑ i |a i x + b i y + c i z + d i | γ i = 1.
(2.1.11)

Unlike superquadrics, hyperquadrics are not restrained to symmetric shapes. Hybrid hyperquadrics [START_REF] Cohen | A hybrid hyperquadric model for 2-d and 3-d data fitting[END_REF] further extend superquadrics by allowing to create some local concavities in hyperquadrics. However, small local details remain inaccessible to these models and many authors have used them as coarse starting points to be refined in a second stage (e.g. using free-form deformations).

Landmark-based models

Statistical shape priors are built from sets of training annotated images, and aim at generating a subspace S of valid solutions capturing most of the shape variability while remaining as small as possible. Once such a subspace is generated, the model is constrained to evolve within it or to remain close to it, thus ensuring a valid shape is always kept. For contour models, such priors generally require to define a set of consistent boundary points, called landmarks (cf. figure 2.1.6).

The Active Shape Model (ASM) [START_REF] Cootes | Active shape models-their training and application[END_REF] uses a Principal Component Analysis (PCA) performed on the position of the landmarks to generate a sub-space S of solution containing the main variation modes (cf. figure 2.1.7). Let {x k } k=1...K be the set of vectors containing the landmark positions for the training data. The average models and its covariance are:

x = 1 K K ∑ k=1 x k , (2.1.12) Σ = 1 K K ∑ k=1 (x k -x) (x k -x) T .
(2.1.13)

Then, PCA on Σ gives:

Σ ≈ U T diag (Λ) U, (2.1.14)
where Λ is the vector of retained eigenvalues and U is the corresponding set of retained eigenvectors stored in columns. Then, any shape lying within the prior subspace of solutions as determined by the set of boundary landmarks writes:

x = x + Up, (2.1.15)
where p is the coordinate vector of x in S. Such formulation assumes that the distribution of realistic shapes can be effectively approximated by a Gaussian distribution. During model evolution, which is similar to a gradient descent, the boundary is regularly projected into subspace S to ensure the model shape always remain realistic with respect to the training examples.

A similar strategy can be adopted, but using the full covariance matrix Σ instead of a subspace via a selection of the main variation modes [START_REF] Cremers | Diffusion Snakes: Introducing Statistical Shape Knowledge into the Mumford-Shah Functional[END_REF]. Shape correctness is enforced by minimizing an energy term in the form of a Mahalanobis distance, in addition to the external energy term:

E shape (x) = (x -x) T Σ -1 ⊥ (x -x) , (2.1.16)
where Σ ⊥ is the regularized covariance matrix, authorizing inversion. Like ASM, such model implies that the distribution of shapes be Gaussian, and penalizes -but does not forbid -more strongly shapes having a low probability. However, for many classes of shapes, the Gaussian hypothesis leading to linear shape statistics is not valid (cf. figure 2.1.8). In order to circumvent this limitation, a typical solution is found in using constraints in a feature space via kernel methods. Let φ : Ω → F be a non-linear mapping from the input space to a higher-dimensional feature space, potentially unknown. Using a Kernel Principal Component Analysis (KPCA), a Kernel Active Shape Model was proposed [START_REF] Romdhani | A Multi-View Nonlinear Active Shape Model Using Kernel PCA[END_REF][START_REF] Twining | Kernel Principal Component Analysis and the Construction of Non-Linear Active Shape Models[END_REF]. As in the linear ASM, the optimization procedure of [START_REF] Romdhani | A Multi-View Nonlinear Active Shape Model Using Kernel PCA[END_REF] consists in regularly projecting the mapped boundary back into the feature space. [START_REF] Twining | Kernel Principal Component Analysis and the Construction of Non-Linear Active Shape Models[END_REF] enforced the shape prior by minimizing an energy based on a "proximity to data" measure in the feature space. Using a formulation similar to that of equation 2.1.16, a non-linear shape prior was proposed (Cremers et al., 2006a) such that the empirical mean and covariance are computed on the mapped training vectors.

muscle segmentation using a hierarchical statistical shape model A hierarchical prior model using Diffusion Wavelets was proposed to segment organs [START_REF] Essafi | Hierarchical 3D diffusion wavelet shape priors[END_REF] -including one calf muscle -in MRI. This model builds on the formulation of the ASM (cf. equation 2.1.15), using a differ-ent basis for the subspace of valid solutions. One of the main drawbacks of ASMs, is that they often require a large number of training data in order to obtain relevant decomposition modes. Indeed, some non-correlated shape features -such as global and local shape configurations -are often modeled by the same deformation modes. Thus, desired shape behaviors are often mixed with unwanted shape behaviors when optimizing the shape parameters for segmenting a new image. The hierarchical approach allows to uncorrelate small and large scale shape behaviors. Moreover, the presented method also uncorrelates long-range shape behaviors, thus ensuring that deformation mode are spatially localized.

We give a brief summary of this method. First, a graph G (V , E ) is built on the set of landmarks: V is the set of nodes and each landmark corresponds to a node in V; E is the set of edges, whose weights are determined through a statistical analysis of the mutual relations between the landmarks in the training set {x k } k=1...K (cf. Shape Maps [START_REF] Langs | Modeling the structure of multivariate manifolds: Shape Maps[END_REF]). As a result, landmarks with independent behaviors will be connected by edges with a small weight, whereas nodes with strongly related behaviors -such as neighboring points -will be connected by large weight edges.

Second, a Diffusion Wavelet decomposition of G is performed. This process involves computing the diffusion operator T of graph G, which is the symmetric normalized Laplace-Beltrami operator, and computing and compressing the dyadic powers of T. The output of this decomposition is a hierarchical orthogonal basis {Γ i } for the graph, whose vectors correspond to different graph scales; considering the vector of landmark positions when decomposed on the new basis:

x = x + Γp,
(2.1.17) global deformations -i.e. global relations between all the nodes -are controlled by some of the coefficients in p, while local interactions -i.e. local interactions between close-by nodes -are controlled by some other coefficients in p. Projecting all the training examples onto this new basis, a PCA is performed at each scale of the decomposition. Finally, during the segmentation process, the landmarks are positioned on the target image in an iterative manner: 1) the position of the landmarks is updated according to a local appearance model; 2) they are projected into the hierarchical subspace defined previously.

Implicit deformable models: level sets

Level set methods, also called Geometric Deformable Models, are a special class of Deformable Models. They were designed to solve certain issues of standard Active Contours (also called Parametric Deformable Models or Snakes): the explicit parametrization of Snakes necessitates frequent re-parametrization in order to preserve accuracy; topo- logical changes in snakes are difficult to achieve; the output depends heavily on the contour's initialization (position and shape); region features (color, texture) are difficult to take into account.

The level set method, proposed by [START_REF] Dervieux | A finite element method for the simulation of a rayleigh-taylor instability[END_REF]; [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF] in the field of fluid dynamics, uses an implicit representation of the evolving boundaries C as the location of the zero level line of an implicit function φ : Ω → R (cf. figure 2.1.9):

C = {x ∈ Ω|φ (x) = 0} .
(2.1.18) Such implicit representation allows natural topological changes of the implied boundary (cf. figure 2.1.9). This representation is then combined with the general energy functional for Active Contour models using parametric curves C : [0, 1] → Ω (Paragios, 2005; Cremers et al., 2006b) (cf. section 2.1.1):

E (C) = ˆ1 0 αE int (C (p)) + βE img (I (C (p))) + γE ext (C (p)) dp, (2.1.19)
where I is the image, E int imposes curve smoothness (regularization term), E img pushes the curve towards boundaries (data term), E ext encodes user knowledge (interaction or prior term), and α, β, γ are balancing parameters. Such cost functional can be minimized within a finite difference approach on the variational form:

∂C ∂t = - ∂E (C) ∂C = F • n, (2.1.20)
where F is a speed function and n is the contour unitary normal vector. Incorporating the level set representation to the previous formula yields:

∂ ∂t φ = -|∇φ| F. (2.1.21)
Such a flow specifies the evolution of the implicit function φ at the location of the boundaries. Among the first applications of the level set formalism, a notable contribution is the Geodesic Active Contours model [START_REF] Caselles | Geodesic active contours. International journal of computer vision[END_REF][START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF], with evolution equation:

∂ ∂t φ = |∇φ| div g (I) ∇φ |∇φ| , (2.1.22)
where g is an edge function, increasing with the strength of the edges in I. Such equation models an inward (respectively outward) flow which shrinks (respectively expands) the initial curve towards the detected edges. For implementation, φ is generally considered as a signed distance function to the evolving curve:

φ (x) =          0 , x ∈ C, D (x, C) , x ∈ C in , -D (x, C) , x ∈ C out , (2.1.23)
where D(x, C) is the Euclidean distance from x to C, C in ⊂ Ω is the region inside C and C out ⊂ Ω is the region outside C.

Having to impose the direction of the flow, as well as relying only on edge detection make such formulation very dependent on the initialization and non-robust to undesired local extrema. To overcome these problems, a global regional term was introduced. For a two-phase segmentation, such model has an energy of the form [START_REF] Chan | Active contours without edges. Image Processing[END_REF][START_REF] Cremers | A Review of Statistical Approaches to Level Set Segmentation: Integrating Color, Texture, Motion and Shape[END_REF]):

E (φ, θ 1 , θ 2 ) = ˆΩ -H (φ) log p (I|θ 1 ) (2.1.24) -(1 -H (φ)) log p (I|θ 2 ) + w |∇Hφ| dx,
where p (I|θ i ) is the posterior probability to observe I given region model parameter θ i , H is the Heaviside step function such that: H (x) = 1, x ∈ C in , and H (x) = 0 otherwise -in practice, H is a smooth differentiable approximation of the Heaviside distribution. Region parameters θ i can either be set, as a prior knowledge of statistical properties of the objects region [START_REF] Paragios | A PDE-based level-set approach for detection and tracking of moving objects[END_REF], or estimated in an alternate optimization method, assuming intensity piece-wise constant objects [START_REF] Chan | Active contours without edges. Image Processing[END_REF]. Many variations on this formulation were developed, for multi-phase seg-mentation [START_REF] Zhao | A variational level set approach to multiphase motion[END_REF][START_REF] Paragios | Geodesic Active Contours and Level Sets for the Detection and Tracking of Moving Objects[END_REF][START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF], for using non-parametric approximation of density functions [START_REF] Kim | Non parametric methods for image segmentation using information theory and curve evolution[END_REF][START_REF] Rousson | A Variational Framework for Active and Adaptative Segmentation of Vector Valued Images[END_REF], instead of Gaussian approximation functions.

Level sets and shape priors

Incorporating shape prior knowledge within the level set framework has been done in multiple ways. Some methods [START_REF] Chen | Using prior shapes in geometric active contours in a variational framework[END_REF][START_REF] Rousson | Shape priors for level set representations[END_REF] use the statistical average curve C M , and associated signed Euclidean distance function φ M to guide the deformation of the model. For example, when added to a level set energy functional, minimizing the term:

E shape (φ, A, s) = ˆΩ δ (φ) (sφ (x) -φ M (A (x))) 2 dx (2.1.25)
imposes the implicit function φ to be similar to φ M , where δ is a smooth approximation of the Dirac distribution, A is a linear or non-linear transform to be estimated, s is the scale factor between φ and φ M to be estimated.

As seen previously with explicit models, single reference models can only recover small deviations from the prior shape. In order to build a model able to recover important shape variations, methods using a PCA on a set of distance functions were developed [START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF][START_REF] Tsai | Model-based curve evolution technique for image segmentation[END_REF], similar in principle to the Active Shape Model (cf. section 2.1.2):

φ λ = φ M + ∑ j λ j U j , (2.1.26)
where λ j is the parameter weighting the eigenvectors U j . Similar to 2.1.25, an energy term enforcing a shape prior can be formulated:

E shape (φ, A, s, λ) = ˆΩ δ (φ) (sφ -φ λ (A)) 2 dx.
(2.1.27)

Such an approach assumes that the distribution of training shapes is Gaussian, which is often not true in practice, and thus necessitates a large number of eigenvectors in order to capture enough details of the shape. Approaches for modeling non-Gaussian shapes distributions were considered, using nonparametric density estimators (Cremers et al., 2006a;[START_REF] Rousson | Efficient kernel density estimation of shape and intensity priors for level set segmentation[END_REF], e.g. yielding prior energy terms of the form [START_REF] Rousson | Efficient kernel density estimation of shape and intensity priors for level set segmentation[END_REF]): .1.28) where λ i = (λ 1 , . . . , λ n ) is the shape vector (e.g. projection into a PCA space) for training shape i, λ is the shape vector of the model, K is a kernel function (e.g. Gaussian kernel K (x) ∝ exp -x 2 ) and σ is a scaling parameter. Such approach ensures the model remains in the vicinity of the training shapes.

E shape (φ, λ) = -log 1 N N ∑ i=1 K λ -λ i σ , ( 2

Graph-based surface models

Unlike most surface-based models, which are usually optimized with continuous local methods, these models are fit to the image via discrete graph-based optimization procedures. The graph frameworks has the advantage of allowing to enforce similarity invariant statistical shape priors while searching for a globally optimal solution.

surface models and pair-wise discrete optimization A landmark-based surface model was proposed by [START_REF] Besbes | Shape priors and discrete mrfs for knowledge-based segmentation[END_REF]; [START_REF] Besbes | Image Segmentation using MRFs and Statistical Shape Modeling[END_REF] for a 2-class segmentation methods with a shape prior based on learning joint probabilities between pairs of control points. Let us call x i the spatial location of the i-th control points, X = {x 1 , . . . , x N } the set of control points. Given an annotated training set of sample images, the joint probability for the relative location of two control points is learned based on the normalized distance between pairs of points:

d ij = x i -x j d , (2.1.29)
where d is the estimated scale of the object. Since prior information relies only on a normalized distance between the control points, it is rotation, translation and scale invariant -although scale-invariance depends on a separate scale estimation. From these distance, a joint probability distribution is estimated: p ij = p d ij given a statis- tical model (e.g. a Gaussian model). A graph G = (V , E ) is then built over the control points. The following functional is defined: (2.1.30) where: V i (•) is the data inference term for each control point. This measures the adequacy of the segmentation with the pixels intensity; V ij (•, •) is the shape prior term, measuring the adequacy of the surface's shape with the prior knowledge:

E (X ) = ∑ v i ∈V V i (x i ) + ∑ e ij ∈E V ij x i , x j ,
V ij x i, x j = -log p ij . (2.1.31)
The most appropriate configuration for the set X can be obtained by minimizing the discrete energy E: each control point is limited to a finite number of displacements. As formulated, this problem is NP-hard and can only be approximately solved, e.g. via an efficient primal-dual algorithm [START_REF] Komodakis | MRF Optimization via Dual Decomposition: Message-Passing Revisited[END_REF]. Due to the limited set of locations for each control point, this optimization method remains non-global, but nonetheless covers a greater range of solutions than variational methods. muscle segmentation using higher-order cliques In the previous method, a shape prior is enforced via learning distribution models for pairs of control points. While such modeling makes a rotation and translation invariant model, it is only approximately scale invariant due to the normalization of distances by an estimation scale of the model. Moreover, the iterative approach for minimizing the functional imposes an initial guess for the control point location, e.g via user interaction.

A surface model using prior distributions over higher-order cliques -as triplets -was proposed by [START_REF] Wang | 3D knowledge-based segmentation using pose-invariant higher-order graphs[END_REF]. This allows to define an exact scale-invariant model and allows global search for control point locations due to the more constraining effect of imposing triplet configurations. Let us denote the spatial distribution for triplet (i, j, k) as p ijk , which is learned for all triplets of points, using a mixture of Gaussians statistical model.

Unlike the previously presented method, control point possible locations are searched in a first stage, within the entire image. Such approach is achieved by capturing the local image appearance around the control points in the training set images. A series of Gabor filters are computed over patches surrounding the landmarks and fed to a Randomized Forest Classifier. Then the trained classifier is run over the target image to determine a set of best matches for each landmark and a score for each match. The k-th match for landmark i represents the assignment label l k i for node v i . Defining x i (l i ) as the location of point i for label l i , the following functional is defined:

E (L) = ∑ v i ∈V V i (x i (l i )) + ∑ e ij ∈E V ijk x i (l i ) , x j l j , x k (l k ) ,
(2.1.32)

where: V i (•) is the data inference term for each control point. This is defined as a decreasing function of the classification score for the corresponding point match;

V ij (•, •, •) is the shape prior term, measuring the adequacy of the surface's shape with the prior knowledge:

V ijk x i (l i ) , x j l j , x k (l k ) = -log p ijk x i (l i ) , x j l j , x k (l k ) . (2.1.33)
This functional contains high-order cliques and can be minimized using a novel Dual-decomposition method [START_REF] Komodakis | Efficient training for pairwise or higher order crfs via dual decomposition[END_REF]. Since the possible location for each landmark are searched for globally, no iterative approach is required and the model does not necessitate an initial guess.

region-based methods

Apart from methods relying on surface representations, there also exists segmentation methods using regions models: each point in the image domain is parametrized by an assignment value, which can be hard -a label -or soft -e.g. a probabilistic label assignment. In general, these methods achieve good accuracy at segmenting the boundaries of the objects, but global shapes and topology are much more difficult to constrain. This is the opposite with surface-based methods which provide natural means for enforcing shape constraints and topology but often fail to segment the boundaries accurately. In the case of multi-object segmentation, surface-based methods necessitate complex computer animation algorithm to deal with interacting objects -e.g. to solve intersection issues. Region-based methods offer a simpler framework in that matter, since one point in the image domain can only belong deterministically to one label. Moreover, fuzzy or probabilistic assignment to different labels is also simpler with such methods.

Continuous methods

A continuous segmentation model based on probabilistic atlases was introduced by [START_REF] Cremers | Shape priors in variational image segmentation: Convexity, Lipschitz continuity and globally optimal solutions[END_REF]. This approach was motivated by the well-known limitations of level sets with shape prior methods: a) the cost functionals are generally non-convex, leading to non-optimal solutions and b) it was noticed that linear combinations of signed distance functions -used in building statistical shape priors -do not yield valid shapes and thus poor shape priors.

This approach is based on a probabilistic assignment model. Let us denote q : Ω → [0, 1] such that q (x) is the probability that x belongs of the target object. Assuming we own a set of annotated training images, denoted T = {q k } k , a PCA is performed on T . This operation involves computing the average probabilistic segmentation:

µ = 1 |T | ∑ k q k , (2.2.1)
and the covariance matrix of the centered samples:

Σ = 1 |T | ∑ k (q k -µ) (q k -µ) T . (2.2.2)
Then a segmentation in the subspace spanned by the first n eigen-modes is:

q γ = µ + Γγ, (2.2.3)
where Γ ∈ R N×n denotes the matrix whose columns are the n retained eigenvectors. With these definitions, the following general cost functional is proposed:

E (q γ ) = E i (q γ ) + αE s (q γ ) . (2.2.4)
where:

E i (q) = ˆf (x) q (x) dx (2.2.5) + ˆg (x) (1 -q (x)) dx + ˆh (x) |∇q (x)| dx,
where f (x) and g (x) are functions binding the model to the image, and h (x) a func- tion constraining the length of the boundary. Such functional is a generalized version of the typical two-phase level set cost functional 2.1.24 [START_REF] Chan | Active contours without edges. Image Processing[END_REF]. However, unlike the level set approaches, this functional is convex in q (whose domain is also convex). The shape prior functional is:

E s (q γ ) = γ T Σ -1 γ, (2.2.6)
which is convex in γ. Such approach has the advantage that the shapes spanned by the PCA are always valid, unlike the level set case. It should be noted that this approach assumes a Gaussian model for the probabilistic assignment variables, which is a rather rough approximation and, as a consequence, the PCA does not generate probability vectors.

muscle segmentation using a continuous region model The previous method was extended to multi-label segmentation by Andrews et al. (2011b), and applied to skeletal muscle segmentation (Andrews et al., 2011a). Before performing the PCA on the training samples, an Isometric Log-Ratio (ILR) transform is applied to the assignment vectors. The reason for using this transform is that multi-label segmentation requires to have probabilities at all times, which the previous method does not achieve. Here, the PCA is performed in the ILR space and its output is projected back into the initial probability space. Denoting η γ = µ + Γγ a segmentation in the subspace of valid solution spanned by the PCA in the ILR space, the following functional is proposed:

E (η γ ) = d (η BG, η γ ) 2 + ˆ(1 -h (x)) |∇η γ | 2 + γ T Σ -1 γ, (2.2.7)
where d (η BG , •) is an intensity prior functional for separating muscle voxels from background voxels, and h (x) is and edge-map of the target image such that the energy is minimal when the boundaries of the model match the edges in the image.

Graph-based methods

In the following we review region-based methods relying on a graph formalism. Digital images are intrinsically discretized: in space into pixels and in time into frames (videos) -amplitude quantization is generally too fine to be considered. For this reason, pixels/voxels decomposition provides a convenient solution domain and are often used as such. However, other sub-domains have been used by further reducing the search domain into sets of pixels with the same assignment value (super-pixels) . These graph-based methods are intrinsically "bottom/up", since they rely on pixelscale features to provide a segmentation of larger scale objects. More importantly, graph-based methods offer means of obtaining globally optimal segmentations, which is not the case for most surface-models, whose results depend heavily on the initial conditions.

notations Let us consider an image I with N pixels. A segmentation on I is formulated as a labeling problem of an weighted graph G = (V , E ), where V is the set of nodes and E ⊂ V × V is the set of edges. The i-th node v i corresponds to the i-th pixel of image I. An edge spanning two vertices v i and v j is denoted e ij , and its weight w ij . The set of edges is composed only of pairs of adjacent pixels -e.g 4-connected neighborhood structure in 2D, and a 6-connected structure in 3D.

Given a set of labels S (e.g. the indices of the muscles), segmenting I consists of assigning a label s ∈ S to each node p ∈ V. In the binary case, S = {0, 1} . We denote x i ∈ S the assignment label of pixel i.

Segmentation with Graph Cuts

The Graph Cuts method was originally developed for binary clustering of images [START_REF] Wu | An optimal graph theoretic approach to data clustering: Theory and its application to image segmentation[END_REF]Shi and Malik, 2000b), and for interactive segmentation [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF] where a user provides annotated pixels, the seeds, as an initialization for the algorithm. In the following, we briefly explain the principles of the graph cut theory for the binary case -i.e. a foreground/background segmentation.

Two terminal nodes, the source S and the sink T are added to the set of nodes V. A cut on graph G is a separation of the nodes into two subsets, one V S containing the source, and the other V T the sink (cf. figure 2.2.1). We call C ⊂ E the subset of edges crossing the separation -i.e. with one vertex in V S and the other in V T . The cost |C| of this cut is the sum of the weight of the edges in C:

|C| = ∑ (i,j)∈C w i,j . (2.2.8)
The Ford-Fulkerson max-flow/min-cut algorithm can be used to find the minimum cut in polynomial time.

Segmenting image I is a (binary) labeling problem which amounts to defining a cut on graph G, where each node in V S are assigned to one label (e.g. the background) and each node in V T are assigned to the other label (e.g. the foreground). An optimal segmentation is one which verifies certain criteria: e.g. intensity homogeneity ("I know this object has the same color everywhere"), intensity prior ("I know this object is dark"). These simple criteria can be expressed at the pixel scale in terms of edge weights:

• homogeneity of intensity between two pixels can be measured using a Gaussian weighting of intensity differences:

∀ (i, j) ∈ E \ {S, T} , w ij = exp -β I i -I j 2 .
(2.2.9)

where β is a contrast parameter;

• prior knowledge over intensity of one pixel can be measured using the negative log-likelihood:

∀v i ∈ V , w iS = -log Pr (I i |T) , w iT = -log Pr (I i |S) , (2.2.10)
where Pr (I i |T) (respectively Pr (I i |S)) is the intensity distribution for the object attached to the sink (respectively the source). Note that if Pr (I i |T) ≈ 0, then w iS ≫ 0, which implies that cutting the edge e iS is very costly.

As a consequence, an optimal segmentation according to pair-wise criteria can be obtained by finding the minimum cut Ĉ, which is a cut of G which has a minimal cost. Indeed, minimizing the cost of the cut implies cutting the edges with the smallest weight on average, which is equivalent of satisfying the selected criteria as well as possible. An optimal segmentation is obtained via minimizing the following functional with respect to the label assignments:

E (x) = ∑ (i,j)∈E w ij x i -x j , , (2.

2.11)

where x = {x S , x T } ∪ {x i } i=1...N and where x i is a binary assignment variable:

x i =    1 if v i ∈ V S , 0 if v i ∈ V T .
(2.2.12)

Under this formulation, the assignment value for terminal nodes is of course already known: x S = 1 and x T = 0. Additional information can be provided in the form of seeds, i.e. used-assigned pixels. Let us call V M the set of marked nodes, then finding the optimal segmentation is achieved by minimizing 2.2.11 with respect to the assignments for un-marked nodes v ∈ V \V M .

In the general case, w ij x i , x j can depend on the labels x i and x j . It was shown that Graph Cuts achieve a globally optimal solution providing that w ij (•, •) is a submodular function [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], i.e. that it verifies:

w ij (0, 0) + w ij (1, 1) ≤ w ij (1, 0) + w ij (1, 0) .
(2.2.13)

When one wishes to segment several objects at once, the multi-label case can be formulated as a multi-way cut problem [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. However, such formulation was proven to be NP-hard and only good approximate solutions can be found [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Komodakis | MRF Optimization via Dual Decomposition: Message-Passing Revisited[END_REF].

A different approach with no required intensity prior [START_REF] Wu | An optimal graph theoretic approach to data clustering: Theory and its application to image segmentation[END_REF]Shi and Malik, 2000b) can be formulated by setting: ∀i, w iS = w iT = constant. Moreover, seeds are no longer required, so the resulting segmentation depends only on the intensity difference between pixels. Such method is called perceptual grouping. However, the output segmentation tends to be biased towards very small segments, because small segments have less edges and induce a lower cut cost [START_REF] Wu | An optimal graph theoretic approach to data clustering: Theory and its application to image segmentation[END_REF]. This problem can be solved with normalized cuts (Shi and Malik, 2000b) where a different functional is minimized, consisting in a normalized cut cost by the cost of all edges within the parts:

(cf. figure 2.2.2),
E nC (x) = ∑ (i,j)∈E w ij x i -x j 1 ∑ (i,j)∈E w ij x i + 1 ∑ (i,j)∈E w ij (1 -x i ) . (2.2.14)
However, minimizing the normalized cut is a NP-hard problem and can only be optimized approximately.

An interactive scheme can be adopted, where the user has to provide seeds to initialize the algorithm [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF] (cf. figure2.2.3). As a result, the segmentation consists in a optimal separation of the seeded regions and is not biased towards small cuts. No intensity prior is required, but one can be estimated simply by computing the histogram of the marked regions, thus enforcing intensity coherence as well as spatial coherence.

Shape priors and graph-cuts

A "star shape prior" was proposed by [START_REF] Veksler | Star shape prior for graph-cut image segmentation[END_REF], allowing to constrain the shape of objects with star shapes. An object O has a star shape if there exists a point c such that any point p in O can be reached in straight line from c without stepping out of O. The shape prior is enforced by setting special weights s ij to edges connecting pixels i and j which both lie on a line passing through c. Assuming pixel j is between c and pixel i, these special weights have the following form:

s ij x i , x j =          0 if x i = x j , ∞ if x i = 1, x j = 0 β if x i = 0, x j = 1, , (2.2.15)
to be added to the original weights w ij . This way, if pixel i belongs to the target object, then pixel j has to belong to the object as well. The center point c is provided by the user.

A "topological" prior was proposed by [START_REF] Delong | Globally optimal segmentation of multi-region objects[END_REF], . Topological constraints are, for instance: containment ("object A contains object B"), exclusion ("object A and object B are not superimposed"), etc. These constraints are enforced in the binary graph-cut framework by solving simultaneous binary problems -one layer per label -and adding edges between the pixels of the different layers. These additional edges enforce the aforementioned topological constraints. Similarly, a topological "location" prior was proposed by [START_REF] Roy Chittajallu | Fuzzy-Cuts: A knowledge-driven graph-based method for medical image segmentation[END_REF], relying on prior knowledge on the location of other objects: e.g. having segmented the lungs in CT-scan of the torso, the position of these organs provides a prior on the location of the heart.

A shape prior based on a single reference segmentation was introduced to the graph-cut framework by [START_REF] Freedman | Interactive Graph Cut Based Segmentation with Shape Priors[END_REF]. Let us denote C the boundary of this template shape, and φ the unsigned distance function whose 0-level set is C: C = {p ∈ Ω|φ (p) = 0}. Then, with p (i) denoting the spatial location of pixel i, the additional binary weight function is used:

s ij x i , x j =    0 if x i = x j , φ p(i)+p(j) 2 otherwise, (2.2.16)
which gives a low cost if the label transition (i.e. the boundary of the object) is located near the boundary of the reference shape.

Segmentation with Random Walks

We refer the reader to chapter 4 for detailed explanations on this approach.

prior art limitations

In this chapter, we have tried to review all existing segmentation methods which could adapt to muscle segmentation, highlighting those which were actually applied to this task. Although there were few attempts at segmenting striated skeletal muscles, the existing methods adopted very different approaches. Among the surface-based models, [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF]; Gilles and Magnenat-Thalmann (2010) adopted the active surface approach: a reference shape is deformed according to internal regularization forces and external image-driven forces in a local optimization process (pp. 17). Shape prior is enforced by a "shape memory" term [START_REF] Gilles | Musculoskeletal MRI segmentation using multi-resolution simplex meshes with medial representations[END_REF] or via a soft tissue simulation algorithm called shape-matching [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF].

Both constraints enforce a strong prior, and the results is always valid in terms of anatomy. However, the local optimization and the possibly unreliable data-terms (gradient driven) sometimes yield over-regularized solutions. The method of [START_REF] Essafi | Hierarchical 3D diffusion wavelet shape priors[END_REF] also uses a surface model, and thus suffers from the same problems. However, due to using a hierarchical statistical prior models, one can expect the surface model to adapt more easily to the specifics of the target image (pp. 23). The method of [START_REF] Wang | 3D knowledge-based segmentation using pose-invariant higher-order graphs[END_REF] has the advantage of being optimized globally, as candidate targets for the surface vertices are searched in all the image. Second, the higher-order clique approach allows to learn a statistical shape prior and to enforce local similitude invariance in a practical way. However, like other landmark-based approaches, it relies on finding consistent detectable anatomical points, which we believe is very difficult with skeletal muscles (pp. 29). Like most surface models, the latter two methods may also suffer from topological issues such as intersection during optimization in a multi-object context, since this case was not addressed. Finally, the method of Andrews et al. (2011a) benefits from its region-based approach: global optimization and no topological issues (pp. 31). Shape prior is efficiently enforced for the entire limb by imposing the shape to evolve in a statistical shape space of valid solutions, built from a training atlas base through principal component analysis. The presented results demonstrate excellent performances. However, one weakness may remain in the image term, which relies on contour detection, an unreliable process in MR images of muscles where contours are often missing or often confused with spurious structures.

In chapter 4, we present some quantitative comparisons of the methods of [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF] and Andrews et al. (2011a) with our own proposed methods.

introduction

The approach presented in this part builds upon the observation that certain regionbased methods can achieve excellent results, providing a substantial manual initialization. Among them, the Random Walks (RW) segmentation algorithm [START_REF] Grady | Random walks for image segmentation[END_REF] has the advantage of addressing the issue of incomplete contours, which is one of the main difficulties in our task. The RW algorithm is a pixel-wise graph-based method which amounts to computing the probabilities of assignment of each pixel to each label (e.g. the index of a muscle). In this approach, a quadratic energy functional measuring the conformity of the segmentation with the image is defined and minimized by solving a large sparse linear system. Since efficient optimization algorithms exist to solve such problem -e.g. the conjugate gradient method -large volumes can be entirely segmented in a few minutes (see chapter 4 for a more in-depth view of the RW method). The interactive phase of the RW method consists in manually annotating some pixels of each label: the seeds (cf. figure 3.1.1). However, while this approach offers a considerable time gain by avoiding to manually annotate all the pixels, manually positioning the seeds remains a time consuming process, especially when large volumes are to be segmented.

In this chapter, we propose a method for generating the seeds automatically. While the large inter-subject variations of skeletal muscles prevent using a deformable model based on a shape template, the relative positions of the muscles are consistent among individuals. Therefore, it is reasonable to develop a segmentation method whose prior information is not based on the position of the boundaries and on the shape of the objects, but rather on the approximate relative position of the objects and the existence of boundaries between them. Seeds are compatible with this idea: it is sufficient that they are placed well inside the target object and the Random Walks algorithm will segment the contours accurately. Through a sampling process, unlabeled seeds are generated in the image domain. The clustering of these seeds with respect to the different muscle classes is achieved through a MRF graph-based approach: each seed is considered as a node and is to be assigned to a label -i.e. a muscle index. The weight of the edges between the seeds is set to reflect the topological information: the presence -or absence -of a boundary between two seeds and the likelihood of the orientation of an edge with respect to the assigned muscle indices. A discrete energy functional is defined to measure the the quality of a global label assignment. Due to being NP-hard, this problem can only be approximately solved; however, we obtained satisfying results using an efficient message passing algorithm [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF]) on a real clinical dataset.

related work

To our knowledge, this is the first method proposing an automatic seed generation based on a graph formulation of the labeling problem. Previous automatic seeding methods may be associated with graph-cuts [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF], seeded region growing [START_REF] L I O G R A P H Y | Seeded region growing. Pattern Analysis and Machine Intelligence[END_REF] or Random Walks. The more straightforward approaches rely on appearance filtering for selecting seeds: only pixels of a certain predetermined intensity or color are kept as seeds [START_REF] Wyatt | Automatic segmentation of the colon for virtual colonoscopy[END_REF][START_REF] Wighton | A fully automatic random walker segmentation for skin lesions in a supervised setting[END_REF]. Similarly, some methods [START_REF] Fan | Automatic image segmentation by integrating color-edge extraction and seeded region growing[END_REF] first detect certain images features (e.g. edges) and generate seeds according to geometrical considerations (e.g. estimated centroids). Integration of knowledge of the anatomy can be based on pre-defined distances and relative positions [START_REF] Wyatt | Automatic segmentation of the colon for virtual colonoscopy[END_REF][START_REF] Maier | Automatic liver segmentation using the random walker algorithm[END_REF], but such approaches require small inter-subject variations or few objects to segment, with simple geometries. Finally, there also exists "pre-segmentation" methods, destined to partition the image into several regions to be latter merged to obtain semantically correct regions [START_REF] Mičušík | Automatic image segmentation by positioning a seed[END_REF]; such methods do not rely on any predefined topology, segmenting an unknown number of objets, with unknown shapes.

prerequisites

We need a binary mask of the muscle tissues, separating them from the background and the fat tissues (see figure 3.3.1). Such a mask can be obtained via intensity thresholding -or in our case, more efficiently, from the "water-map" images we based on for our experiments. This water-map represents the fraction of water protons in the organic tissues; as a result, muscles -made of water for a large part -have high intensity values, whereas fat tissues and bones have intensity values close to zero (cf. appendix A for more details). This mask is used to ensure the seeds are only generated within the muscles. This approach assumes we have a reference atlas rigidly registered with the image to be segmented (see figure 3.3.2). We applied rigid-registration (translation, rotation, scaling) using the free registration software medInria. This software uses the block matching rigid registration algorithm [START_REF] Ourselin | Block matching: A general framework to improve robustness of rigid registration of medical images[END_REF]. To reduce the possibility of registration errors, we used the binary masks of the muscle tissues as input images for the registration software. This way, only the muscle/fat boundaries are taken into account for leading the registration algorithm.

unlabeled seeds sampling

The first step of our algorithm is to generate the seeds that will be used as initialization by the RW algorithm once labeled. We are going to place the seeds well inside the muscles, and not too close to the contours since the task of segmenting the contours accurately is left to the RW method. To achieve this, we first compute an edge map of our volume -e.g. by computing the local variance of the intensity on a small patch centered on each pixel. Then, we sample the seeds on a regular grid, with a high enough density to ensure all the muscles contain several seeds. Finally, we move the seeds in the opposite direction of the gradient of the edge map for a few iterations. More precisely, if the spacing of the grid is k, then the seeds are allowed to move by k/2 from their original position, to prevent any pair of moved seeds to take the same location. This way, seeds are moved away from the close-by edges by a few pixels (see figure 3.4.1).

The seed displacement algorithm can be summed-up by the following algorithm:

1. let us denote the initial location of voxel p as z 0 p ;

2. for i = 1 . . . k 2 do: a) compute normalized displacement vector v i p = -∇I z p / ∇I z p + 0.5 . b) shift seed by v: z i p ← z i-1 p + v i ;
The next step is the label assignment.

mrf formulation

We now describe how to assign a label to the seeds we generated, according to the similarities between the test image and the reference segmentation.

Energy form

We formulate our problem as a labeling problem, that we will solve by minimizing a first-order Markov Random Fields energy functional. Let G = (V , E ) be an undirected graph, where V is the set of nodes and E is the set of edges. Given a set of labels L, we want to assign a label l ∈ L to each node p ∈ V. We denote x p the label assigned to p, and x the collection of all assignments. A first order MRF functional has the form:

E (x) = ∑ p∈V θ p x p + ∑
(p,q)∈V θ p,q x p , x q , (3.5.1)

where the unary potential θ p (•) depends only on the label assigned to node p, and the binary potential θ p,q (•, •) depends on the labels assigned to each nodes of edge (p, q). Unary potentials (or "costs") measure the quality of the label assignment of each node. Similarly, Binary potentials measure the quality of the label assignment of a pair of nodes. The quality measures we are referring to are specific to the desired application. The interest of the MRF formulation is that the assignment of each node depends only on the neighborhood of this node. In the simple case where only pairwise dependencies are considered, efficient optimization algorithm can be used to minimized the previous functional.

In this framework, the nodes in V are the unlabeled seeds and the labels in L are the indices of the muscles in the reference segmentation. In the following, we define the unary potentials so that seeds close to (respectively far from) a muscle in the reference image have a small (respectively large) assignment cost. Binary potentials will reflect the topological properties we seek to exploit: when assigning two different labels (i.e. muscles), unlikely orientations and the absence of boundary between the connected seeds will yield a large binary cost; reciprocally, when assigning the same label to both ends of an edge, it is the presence of a boundary between the connected seeds which will yield a large cost (orientations do not matter is this case). As we will see in a later part, the set of edges will be set to model the spatial proximity of the nodes.

Unary potential

We define a unary cost θ p (•) based on the knowledge that the reference atlas is rigidly registered onto the test image:

θ p (m) = d (p, M m ) , (3.5.2)
where M m is the set of pixels with label m in the reference segmentation and d (•, •) is a distance function between a point and a set:

d (p, R) =    min { p -r , r ∈ R} if p / ∈ R, 0 otherwise. (3.5.3)
The rigid registration ensures that, if we superimpose the test image and the reference segmentation, an unknown muscle in the test image is likely to correspond to a geometrically close muscle in the reference segmentation. Therefore, this term will favor the assignment to a node of a label corresponding to a geometrically close muscle, over the label of a more distant muscle (see figure 3.5.1).

Geodesic distance potential and graph connectivity

The unary term does not take into account the structure of the set of muscles. As we saw previously, the shape of the muscles can vary quite extensively between in-dividuals. On the other hand, the topology of the set of muscles is consistent: muscle a is always "close to" muscle b, but "far from" muscle c, there is most of the time a visible contour between muscle a and b. We designed binary terms to account for this knowledge. The binary potential is the sum of two potentials: θ p,q (•, •) = w geod θ geod p,q (•, •) + w orient θ orient p,q (•, •), weighted by parameters w geod and w orient . The term θ geod p,q (•, •) is derived from the geodesic distance g (p, q) between the seeds p and q. This geodesic distance between two points is understood as the length of the shortest path from one point to another over a curved space. Here, we consider the edge-map of the test image as a curved space: a ground with "hills" -i.e. edges -and "valleys" -i.e. homogeneous regions. The interest of such measure is the following: considering two edges of the same Euclidean length -an edge crossing a boundary has a larger geodesic length than a non-boundary crossing edge. We propose the following definition:

θ geod p,q x p , x q ∝    1/ (1 + exp (γ -g (p, q))) if x p = x q , 1/ (1 + exp (g (p, q) -γ)) otherwise. (3.5.4)
where γ is a free parameter. If x p = x q , i.e. both seeds belong to the same muscle, θ geod p,q

x p , x q will penalize large geodesic distances, as there should be no contour between p and q. On the other hand, θ geod p,q x p , x q will favor large geodesic distances if x p = x q (see figure 3.5.2).

We compute an approximate geodesic distance between the seeds using the Fast Marching algorithm [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] on the test image. The Fast Marching method solves the general static Hamilton-Jacobi equation for a given speed function. Here, the speed is determined by the value of the edge-map (e.g. the gradient magnitude). A propagating discrete front is started at each node and marches in the outward direction at the given spatially-dependent speed. From a technical view, an efficient heap-sort is used to select the best front pixel to update, so that previously evaluated pixels are seen only once. Since there are several fronts evolving simultaneously, the origin point of each front is stored for each visited pixel. Once all pixels have been visited, the arrival time of the front at a given pixel can be seen as a geodesic distance from this pixel to the origin pixel of the front.

We supply the algorithm with the unlabeled sampled seeds. We obtain in return a partition of the image, where each seed in included in a different region, and where the boundary between two regions is equidistant to the seeds in the regions, in the geodesic sense (i.e. a geodesic Voronoi partition). The connectivity of the graph is determined by the Voronoi partition: for each adjacent region, we add an edge between the corresponding nodes. Note that we do not obtain an edge between every pair of seeds, but only between geometrically close seeds. This is not a problem, since θ geod p,q is designed to account for the presence or the absence of a contour between two seeds, and thus is mostly a local term; our unary cost coarsely handles the labeling problem at a larger scale.

Relative orientation potential

The term θ orient p,q (•, •) is intended to ensure that the relative position of a pair of muscles is maintained, e.g. muscle a is located "top right of" muscle b. In a pre-processing stage using the reference atlas, we computed the distribution of the orientations for each pair (m, n) of muscles. More precisely, we estimated the probability density p m,n (•) of random variable U = (P m -P n ) / P m -P n , (3.5.5)

where P m is a randomly chosen point in muscle M m , and

P n = arg min p∈M n P m -p (3.5.6)
(see figure 3.5.3). Thus, p m,n (u) can be seen as the probability that, for any pixel in muscle m, the direction of the closest pixel in muscle n is u. We define the orientation cost as:

θ orient p,q x p , x q ∝    -log P x p ,x q (p -q) if x p = x q , 0 otherwise.
(3.5.7) This cost will be large if the orientation of the edge between two nodes is unlikely given the labels to which the nodes are assigned.

3.6 experimental validation

MRF Optimization

Minimizing the energy 3.5.1 is a NP-hard problem for the binary potential function is arbitrary. However, due to the relatively small number of variables we can expect to reach a good approximate in reasonable time with an efficient message passing algorithm such as the convergent Tree-Reweighed (TRW) Message Passing algorithm [START_REF] Kolmogorov | Convergent tree-reweighted message passing for energy minimization[END_REF]. The TRW algorithm uses a linear programming relaxation technique and operates on the dual of the relaxed problem in order to find a lower bound to the energy of the original problem. Due to its state of the art performance among discrete optimization methods, TRW has been applied to solve many computer vision and medical imaging tasks.

Experimental Results

To evaluate our method, we used a dataset composed of 3D volumes of the right thigh of 15 healthy subjects, covering a wide range of morphologies (8 females, 7 males, ages range: 26 to 60), acquired with a 3T Siemens scanner and using 3pt Dixon sequence (TR=10ms, TE1=2.75 ms TE2=3.95 ms TE3=5.15 ms, rf. flip angle =3°; we used the out-of-phase image), with the an average image size of: 220px × 220px × 64px.

Every volume was manually segmented to obtain the ground truth against which the segmentation results are compared. We focused our evaluation on clinically relevant muscles of the thigh (13 muscles). Sampling the seeds, as described in section 3.4, with grid spacing 7px × 7px × 7px, and 2 iterations to bring the seeds further away from the edges (filter radius = 1px), gave us between 2000 and 2500 seeds per test image. At the MRF stage, detailed in section 3.5, we obtained our best results with the following cross-validation estimated parameters w geod = 0.1, w orient = 0.01, γ = 2 × median p,q∈V g (p, q), and U contains 26 regularly spread unitary orientation vectors. On a 2.8 GHz Intel® processor with 4 GB of RAM, total processing time is under 5 min, which adds up to the RW segmentation time (5 min).

All our results are based on a leave-one-out validation protocol: each volume is used as the reference segmentation for segmenting all the other volumes. First, we evaluated the automatic labeling process by computing the labeling error rate for each muscle (figure 3.6.1) as box-plots1 . For testing the efficiency of the binary potential, we also computed the error rate for an energy only composed of the unary potential. We observe that using the binary potential significantly reduce the error rate. Moreover, the combination of the geodesic potential with the orientation potential is more efficient than using only the geodesic potential, which confirms that both terms exploit complementary types of information.

Then, we computed Dice coefficients to evaluate the segmentation results given by the RW with the labeled seeds (figures 3.6.2 and 3.6.3). The expression of the Dice coefficient is: D = 2 |T ∩ R| / (|T| + |R|), where T and R are the pixel sets of a given muscle in the inferred segmentation and the reference segmentation respectively. The presented method outperforms a simple Atlas-based segmentation method (see chapter 4 for details on this method).

In figure 3.6.4, we show cross-sections of segmentation results. We observe that small muscles tend to be more affected by segmentation errors than large ones, which points out the limitations of our model. In particular, the unary term, based on the superimposition of the registered reference atlas with the target image, is bound to be less effective for small muscles, because the distance term d (., .) is likely to be larger than in the case of large muscles. Besides, many errors are due to the absence of reliable contours between parts of muscles, which, in the case of large errors, indicates that the topology constraints are failing to insure topology or shape correctness.

conclusion

In this chapter we presented a first method relying on the good properties of the RW algorithm when facing weak boundaries. Results on real clinical data show the valid- ity of this approach and the usefulness of the proposed two binary terms. However, an inherent limitation of this method lies in its two-stage formulation: 1) label the seeds;

2) segment the volume from the seeds. Indeed, as shown in figure 3.7.1, if the labeling fails, segmentation will inevitably fail in the same area since seeds act as constraints on the segmentation process (see the Random Walks algorithm detailed presentation in the next chapter). This issue may be revealed by the non-negligible amount of outliers in the Dice plots (Fig. 3.6.2), which means that the method totally failed to segment some muscles (in general, smaller ones, like the biceps femoris or the adductor brevis). In order to avoid such propagation of errors, we propose an integrated method in the next chapter, in which the Random Walks and prior knowledge of shape are combined into one formulation. Although this method is inherently flawed for the aforementioned reasons, some improvements could still be considered. In particular, using non-rigid registration could improve the performance if a better match of the atlas and the segmented image is achieved.

In this chapter, we present different approaches based on the general Random Walks Segmentation algorithm. The state of the art RW algorithm has several benefits over its competitors: (i) unlike discrete methods such as graph-cuts, it provides a probabilistic segmentation (that is, a distribution over all putative labels for each pixel/voxel);

(ii) unlike surface-based methods like Active Contours, it has the ability to handle intersection issues in multi-label cases; (iii) it is robust to partial contours; and (iv) it is computationally efficient, since it only requires solving a sparse linear system of equations. Although it was initially proposed for the interactive setting where a user specifies the seeds of the segmentation, a fully automated extension using the contrast information was proposed by [START_REF] Grady | Random walks for image segmentation[END_REF]. However, contrast information does not permit accurate segmentation in the skeletal muscle case where all muscles share the same appearance. Therefore, we propose to incorporate prior knowledge of the shape of skeletal muscles into the RW framework. In a first published approach (Baudin et al., 2012a), a prior energy functional is added to the RW functional in the form of a proximal regularization of the probabilistic segmentation, whose deviations from an average reference segmentation are penalized. We propose to modulate the strength of the model constraints according to both prior information on the model's local accuracy and to the strength of the contours present in the segmented image. A second published approach (Baudin et al., 2012b), introduces some flexibility to the model by constraining the shape to remain close to a Principal Component Analysis shape space built from training examples. Using the PCA allows us to model complex non-rigid shape variations relying on a few eigen-modes. These methods benefit from the high performance of the Random Walks optimization process, as they require only a simple addition to the original objective function.

4.1 the random walks segmentation algorithm

Preliminaries

From an image I with N pixels, we build a undirected graph G = (V , E ), with V representing the set of nodes, |V | = N, and E the set of edges, where the i-th node v i corresponds to the i-th pixel of image I. We denote the edge connecting the nodes with indices i and j as e ij , and its weight as w ij ≥ 0 . Since the graph is undirected, e ij and e ji denote the same edge and their weight is the same: ω ij = ω ji . The set of edges E is only composed of pairs of adjacent pixels, such that graph G contains only cliques of order 0 and 1. We also denote the neighborhood of pixel i as N i = v j /e ij ∈ E . Given a set of labels S (e.g. the indices of the muscles), segmenting an image I is defined as a graph partitioning procedure, i.e. consisting in assigning a label s ∈ S to each node v ∈ V. We refer to the assignment of node v i to the label s as: l (i) = s.

Random Walker analogy

The Random Walks (or Random Walker) algorithm is a probabilistic method modeling the behavior of a inebriated walker on a non homogeneous ground. At each step, the walker chooses a random direction, where the chosen direction is more likely if it is an easy path -e.g. on a flat ground -than a difficult path -e.g. on a bumpy ground (see figure 4.1.1 for illustration). Suppose the ground is discretized, i.e. the walker can only take a finite number of steps. The ground can be represented by a graph, each node being one possible location. The edges of this graph are the possible paths from one location (node) to another.

The connection with graph partitioning -and thus image segmentation -is established as follows: suppose some of the locations on the ground are marked with a colored flag. The number of colors is finite and there can be several flags of the same color. In graph terms, some nodes in the graph are already assigned to one label, and the number of labels is finite. We refer to the marked nodes as seeds. Then, the Random Walks algorithm computes the probability that the walker reaches one flag before all the others, for each starting location. Finally, each starting location is assigned to the most probably reached-first flag; in graph terms: each node is assigned to the most probable label.

Probabilistic explanation 1

In this framework, the global label assignment is modeled by a random variable. We denote the probability of assignment of voxel i to label s as x s i = Pr (l (i) = s). In this notation, l (i) = arg max s x s i . We denote the transition probability from node i to node j as p ij . In general, we have p ij = p ji .

Assume we possess a set V M of seeds, i.e. pre-labeled nodes for each label (typically, manually marked voxels). We denote the set of unknown nodes as V U , such that

V M ∪ V U = V and V M ∩ V U = ∅.
For convenience in the notations, we denote the sets of indices of the nodes in V U and V M as U and M respectively. Since the label of marked voxels is known for certain, their assignment probability is either 1 or 0

∀i ∈ M, x s i =    1 l (i) = s, 0 l (i) = s. (4.1.1)
We denote to the probability vector for label s as x s , where x s contains the x s i for each voxel i. Then, without loss of generality, we can assume the variables in x s are ordered so that we can write:

x s = x s U x s M , (4.1.2)
where x s U (respectively x s M ) is the vector containing the assignment probabilities for unknown (respectively marked) nodes in V U (respectively V M ). Then, the assignment probability of any unknown node v i ∈ V U can be expressed with respect to the assignment probability of all the nodes (both unknown and marked) in its neighborhood

N i : ∀i ∈ U, x s i = ∑ v j ∈N i p ij x s j . (4.1.3)
The probability of assignment of one pixel depends only on the assignment probabilities of its neighbors. Since all nodes have the same property, the assignment probabilities can all be represented at once in matrix form. We denote the transition probability matrix as Π:

Π =    p ij if e ij ∈ E , 0 otherwise. (4.1.4)
We denote the sub-matrix of Π containing the transition probabilities between two unknown nodes as Π U (respectively: (i)Π M ; (ii)Π UM ; (iii)Π MU ) (respectively: (i) between two marked nodes: (ii) from an unknown node to a marked node; (iii) from a marked node to an unknown node):

Π = Π U Π UM Π MU Π M . (4.1.5)
Thus, we can write:

∀s, x s U = Π U x s U + Π UM x s M , (4.1.6) (I -Π U ) x s U = Π UM x s M , (4.1.7)
with I being the identity matrix of size |U| × |U|.

Transition probabilities

The transition probabilities are to be set accordingly to the desired properties of the segmentation. In general -and such is case of in muscle segmentation -one wishes to segment an image according to the visible boundaries. A boundary in an image is the result of an intensity pattern where many side-by-side pixels of very different intensities form a continuous curve In our probabilistic framework, this leads to have higher probability transitions for pairs of pixels having different intensities and lower probability transitions for pixels having similar intensities (cf. figure 4.1.2). A wellknown choice for transition probability is the Gaussian kernel 2 :

w ij = exp -β I i -I j 2 , (4.1.8)
where w ij is thus defined as the weight for edge e ij . Since p ij is a probability, the transition probability is given by:

p ij = ω ij ∑ i ω ij . ( 4

Random Walks objective functional

Let us define A, the non-normalized transition matrix -also designated as the affinity matrix:

A ij =    ω ij if e ij ∈ E . 0 otherwise, (4.1.10)
and D, the diagonal matrix such that: The Random Walks objective functional is defined as:

D ii = ∑ j ω ij , ( 4
E s RW (x s ) = 1 2 x s⊤ Lx s . (4.1.13)
In the following, we will see that minimizing this functional amounts to solving the probabilistic Random Walks equation 4.1.3, hence determining the assignment probabilities. Since w ij = w ji (it is an undirected graph), the entries of L are:

L i,j =          ∑ k ω kj if i = j, -ω ij if e ij ∈ E , 0 otherwise. (4.1.14)
The Laplacian L has many properties (Von [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], among which the most important ones are:

1. L is symmetric and positive semi-definite; the latter can be seen by reformulating 4.1.13 in a scalar form:

E s RW (x s ) = ∑ e ij ∈E ω ij x s i -x s j 2 , (4.1.15) which gives: ∀a ∈ R N , a ⊤ La ≥ 0; (4.1.16)
2. the smallest eigenvalue of L is 0, and the corresponding eigenvector is the constant one vector 1.

We decompose L into sub-blocks for marked and unknown nodes, hence defining sub-blocks L U , L M and B:

L = L U B B ⊤ L M ,
which allows us to rewrite equation 4.1.13 as:

E RW (x s ) = 1 2 x s⊤ U L U x s U + x s⊤ U Bx s M + 1 2 x s⊤ M L M x s M .
(4.1.17) By differentiating 4.1.17 with respect to the unknown variables, we obtain:

L U x s U = -Bx s M . (4.1.18)
We can show that the solution to this system exists and is unique if their exists at least one path composed of only nonzero-weight edges between each unknown node and a marked node (see appendix B.1).

The system of equation 4.1.18 is equivalent to the probabilistic Random Walks equation 4.1.7. This can be seen by noting that, if D has no zero entry -i.e. each pixel has at least one edge with a non null weight -one can retrieve Π from equation 4.1.9, by normalizing A: where L rw , defined as L rw = I -Π, is called the normalized Laplacian. Decomposing L rw into sub-blocks, we can write:

Π = D -1 A. ( 4
L U = D U (I -Π U ), (4.1.22) B = -D U Π UM , (4.1.23)
where D U is the sub-matrix of D corresponding to unknown indices.

Segmentation algorithm

The Random Walks functional 4.1.13 provides a measure of the probabilistic segmentation x s , and the optimal segmentation is the one which minimizes it, i.e. which has the lower cost according to this measure. We have seen that minimizing the RW functional is equivalent to solving the initial RW system, which we recall here:

∀v i ∈ V U , x s i = ∑ v j ∈N i p ij x s j .
It Once one has inferred the probabilistic segmentation xs for each label s, the "hard" segmentation -i.e. the estimated label l (i) of each node v i -is obtained via:

l (i) = arg max s xs i . (4.1.26)

Related work

The graph Laplacian is a fundamental tool for graph analysis, and is related to many other notions besides Random Walks, as we will see in the following.

spectral clustering In the machine learning community, one finds spectral clustering (Shi and Malik, 2000b;[START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF][START_REF] Meila | A random walks view of spectral segmentation[END_REF], an extremely popular clustering technique. We refer the reader to Von [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] for a detailed overview of many aspects of spectral clustering.

There exists several versions of the spectral clustering technique, involving either the un-normalized Laplacian L, the normalized Laplacian [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] L rw , or the symmetric normalized Laplacian (Shi and Malik, 2000a):

L sym = D -1/2 LD -1/2 . (4.1.27)
Assuming that the data consist of n points x 1 , . . . x n in an arbitrary space, an affinity matrix S -also called similarity matrix -is built from some pairwise similarity function S ij = s x i , x j ≥ 0. A similarity graph G is constructed over the data points and their pair-wise similarities. The first k eigenvectors u 1 , . . . u k of one of the Laplacians of G are computed and stored in a matrix U ∈ R n×k (in [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], the rows of U are normalized to norm 1). Then, denoting the points whose coordinate are the rows of U as y 1 , . . . , y n , clustering with the k-means algorithm [START_REF] Hartigan | Algorithm as 136: A k-means clustering algorithm[END_REF]) is performed on "abstract" points {y i } i .

normalized graph cuts Spectral clustering can be seen as an approximation of some graph-cut problem. Let us denote a partition of the graph with k subsets as A 1 , . . . , A k . With the notation: W (A, B) = ∑ i∈A,j∈B ω ij and Ā as the complement of A, then the min-cut problem consists of minimizing the following quantity:

cut (A 1 , . . . , A k ) = 1 2 k ∑ i=1 W (A i , Āi ) . (4.1.28)
As the solutions to the min-cut problem often tend to create "unbalanced" partitions, selecting excessively small subsets of nodes, some modifications of the objective were introduced; the best-known are the RatioCut [START_REF] Hagen | New spectral methods for ratio cut partitioning and clustering. Computer-aided design of integrated circuits and systems[END_REF] and the normalized cut Ncut (Shi and Malik, 2000b). Denoting the number of node in partition A as|A| and the total weight of the edges in A as vol (A) = ∑ i,j∈A ω ij , the definitions are:

RatioCut (A 1 , . . . , A k ) = 1 2 k ∑ i=1 W (A i , Āi ) |A i | , (4.1.29) Ncut (A 1 , . . . , A k ) = 1 2 k ∑ i=1 W (A i , Āi ) vol (A i ) . (4.1.30)
However, unlike standard min-cut, both problems are NP-hard and can only be solved approximately.

It can be shown (Shi and Malik, 2000a;[START_REF] Hagen | New spectral methods for ratio cut partitioning and clustering. Computer-aided design of integrated circuits and systems[END_REF]) that these quantities can be expressed using the un-normalized graph Laplacian:

RatioCut (A 1 , . . . , A k ) = Tr H ⊤ LH , (4.1.31) with H ij =    1/ A j if v i ∈ A j 0 otherwise
, and the symmetric normalized graph Laplacian:

Ncut (A 1 , . . . , A k ) = Tr H ⊤ L sym H , (4.1.32) with H ij =    1/ vol A j if v i ∈ A j 0 otherwise .
Thus, both problems can be reformulated as trace minimization problems under some constraints. It can be shown that spectral clustering using the un-normalized Laplacian L (respectively the symmetric normalized Laplacian L sym ) solves the relaxation of the RatioCut problem (respectively the Ncut problem (Shi and Malik, 2000a)) when the entries of H are allowed to take arbitrary real values. Please refer to the work of Von [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] for further details and proofs. Although the solutions of the relaxed problems via spectral clustering can be very different from the original problems, spectral clustering remains a popular approximation for its simplicity of use and the quality of its output in practice.

ncut and random walks A relationship exists between Ncut and Random Walks [START_REF] Meila | A random walks view of spectral segmentation[END_REF], which can be seen from the following result:

Ncut (A, Ā) = P ( Ā|A) + P (A| Ā) , (4.1.33)

where P (A|B) is the probability that the random walker jumps from a node in subset B onto a node in subset A. Thus, one can say that minimizing Ncut amounts to finding a partition such that a random walker has the lowest probability of transitioning between parts.

commute time distance

The commute time distance between two vertices v 1 and v 2 is defined as the expected number of jumps of a random walker to travel from v 1 to v 2 and back. Among other interesting properties, the commute distance differs from the shortest path in that it considers all possible paths between the two nodes; as a result, two subsets separated by a partial gap -i.e. whose all but a few inter-connecting edges have low weights -will be considered as "far apart" when using the commute distance, whereas they could be considered as "close-by" with the shortest path distance. The commute time distance can be computed using the generalized inverse of the un-normalized Laplacian L. Recalling the eigen-decomposition of the Laplacian L = UΛU ⊤ , the generalized inverse of L is given by:

L † = UΛ † U ⊤ , (4.1.34) with Λ † ii =    1/Λ ii if Λ ii > 0 0 otherwise .
Then, the commute time distance between nodes v i and v j [START_REF] Saerens | The principal components analysis of a graph, and its relationships to spectral clustering[END_REF]) denoted as c ij , is given by:

c ij = vol (V) L † ii -2L † ij + L † jj . (4.1.35)
As a consequence one can define an embedding of the data points x 1 , . . . , x n , denoted as z 1 , . . . , z n and given by the rows of matrix Λ † 1/2 U, such that the commute time distances between points x i are converted to Euclidean distances between points z i :

c ij = vol (V) z i -z j 2 . (4.1.36)
This embedding is similar to that of the spectral clustering (where the data points are mapped to the rows of U), although it does not yield identical results in clustering.

laplacian eigenmaps The Laplacian eigenmaps algorithm [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] is used to perform dimensionality reduction of low-dimensional manifolds embedded in a high-dimensional space. The embedding of the sample points x 1 , . . . , x n is given by the rows of U which contains the k eigenvectors corresponding to the first k smallest eigenvalues for L sym (same embedding as the spectral clustering of (Shi and Malik, 2000a)). It can be seen [START_REF] Ham | A kernel view of the dimensionality reduction of manifolds[END_REF]) that the Laplacian eigenmap algorithm is equivalent to kernel Principal Component Analysis (PCA) using kernel L † .

random walks with prior knowledge

One can say that the standard formulation of the Random Walks segmentation method relies on prior knowledge in the form of pre-segmented nodes: the seeds. In general seeds are provided by the user, who manually segments a few pixels of each class. While this method is convenient and much more efficient than segmenting the whole image manually, it is still impractical for large datasets. For this reason, one can be interested in completely automatizing the segmentation process by using different, unsupervised types of prior knowledge such as statistical intensity of shape models. Moreover, even in the case of manual seeding, incorporating additional information could improve the segmentation results and thus also decrease the user time.

Random Walks segmentation with appearance prior

Prior knowledge to the Random Walk formulation was introduced by [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF], providing an estimate of the probability distribution of the gray-level intensity for each label. An an example of such technique can be found in figure 4.2.1, which allows using very few seeds to separate organs from the background, using the differences in intensity distributions.

Let us introduce a s (•) as the probability density distribution for label s. Then, denot- ing I i the intensity at pixel i, the Bayes' rule states that the probability of assignment to s for pixel i is:

x s i = a s (I i ) ∑ t a t (I i ) . (4.2.1)
Denoting the vector containing the a s (I i ) for each pixel i as a s = [a s (I 1 ) , . . . , a s (I N )] and A s = diag (a s ), the previous equality can be rewritten as

∑ t A t x s = a s . (4.2.2)
The solution to this system is the minimum of the aspatial functional:

x s⊤ ∑ t =s A t x s + (1 N -x s ) ⊤ A s (1 N -x s ) . (4.2.3)
This prior-appearance functional is combined with the Random Walks functional 4.1.13 with the introduction of a free balancing parameter w app :

E s RWapp (x s ) = x s⊤ Lx s (4.2.4) + w app x s⊤ ∑ t =s A t x s + (1 N -x s ) ⊤ A s (1 N -x s ) .
Intuitively, minimizing this functional will yield a segmentation influenced by both the Random Walk principle -segmentation along boundaries -and by the appearance prior -privileging label assignments that are consistent with intensity distributions.

The optimal segmentation minimizing the functional 4.2.4 is obtained through solving:

L + w app ∑ t A t x s = w app a s . (4.2.5)
Since L is positive semi-definite, system of equations 4.2.5 can be solved even when all entries of x s are unknown (no seeds). However, as noted by [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF], it is still possible -and useful -to use seeds in combination with a prior to obtain more robust segmentations. In the case of muscle segmentation however, the intensity distributions of the labels (the muscles) are extremely similar -except for the background -and thus would make a very weak prior if used alone.

Probabilistic framework 3

In this section, we show that the intensity prior of previous section 4.2.1 is consistent with the probabilistic framework of the Random Walks method, and allows incorporating prior knowledge of shapes as well. As we saw previously in the case of standard Random Walks, prior information is incorporated into the seeds, in the form of presegmented nodes. With the probabilistic notations, this is equivalent to setting the value of the assignment probability for marked nodes to either 0 or 1 -i.e. enforcing the certainty that each one of these pixels belongs to one chosen label. At present, we release these constraints by providing no seed, but add an additional layer of nodes to the graph G with a one-to-one correspondence with the original layer of nodes (cf. figure 4.2.2). The nodes in this second layer V 0 are called the prior nodes, while the nodes in the original layer V are the primary nodes. We set one directed edge from each primary node to its corresponding prior node. We denote the transition probability from primary node v i to prior node v 0 i as λ i , and the assignment probability to label s of prior node v 0 i as x s 0 i ; no edge is set between the prior nodes. Then, the Random Walks node-wise probability (for primary nodes) becomes:

∀v i ∈ V , x s i = (1 -λ i ) ∑ v j ∈N i p ij x s j + λ i x s 0 i , (4.2.6)
where the probability p ij of transition from primary node v i to primary node v j in the original RW equation 4.1.3 is multiplied by the probability to not walk onto the prior node v 0 i : (1λ i ).

Denoting the diagonal transition matrix from primary nodes to prior nodes asΛ = [λ 1 , . . . , λ N ] T and the transition matrix between primary nodes as Π, we reformulate the Random Walks equation 4.2.6 in matrix form

x s = (I -Λ) Πx s + Λx s 0 (4.2.7) [(I -Λ) (I -Π V ) + Λ] x s = Λx s 0 . (4.2.8)
We recall that D is the diagonal matrix such that D ii = ∑ j w ij , and that the com- binatorial Laplacian of the primary graph is given by L = D (I -Π) . Then the new Random Walks equation is equivalent to

[(I -Λ) L + ΛD] x s = ΛDx s 0 , (4.2.9)
which can be solved if the transition probabilities to prior nodes are non-zero, ensuring that (I -Λ) L + ΛD is positive definite.We can rewrite the previous as

(L + ∆) x s = ∆x s 0 , (4.2.10) with the diagonal matrix Ω = (I -Λ) -1 ΛD, (4.2.11)
which excludes setting λ i = 1 without loss of generality, since such setting would make the problem trivial by considering only the prior model and ignoring the Laplacian.

The solution to the last equation 4.2.10 provides the minimum of the following objective:

E s RW+prior (x s ) = x s⊤ Lx s + (x s -x s 0 ) ⊤ Ω (x s -x s 0 ) , (4.2.12) = x s⊤ Lx s + x s -x s 0 2 Ω , (4.2.13)
where . Ω is a weighted norm by Ω. Intuitively, the optimal segmentation minimizing this functional is influenced by both the Random Walk principle -segmentation along boundaries -and by the prior model -privileging label assignments consistent with the prior probabilities. Any solution to this problem is a probability distribution if x s 0 i s is a probability distribution (please refer to appendix B.3):

∑ s x s 0 = 1 ⇒ ∑ s x s = 1. (4.2.14)
Going back to the appearance prior model of section 4.2.1, we set the prior probabilities as: (4.2.15) and Ω as: (4.2.16) retrieving the system for appearance prior 4.2.5.

x s 0 = ∑ t A t -1 a s ,
Ω = w app ∑ t A t ,

random walks segmentation with shape prior

Constant prior weight

In this section we propose a simple shape prior model for the Random Walks based on the general formulation 4.2.12. Suppose that we have an estimate for x s 0 , based on learning a pixel-wise shape over previous segmentations in a training set T (see figure 4.3.1). We set the transition probabilities from the primary nodes to the prior nodes in Λ such that Ω is the constant diagonal matrix:

Ω = w shape I, (4.3.1)
with w shape > 0. Since Ω = (I -Λ) -1 ΛD, this is achieved by setting the transition probability λ i to the prior nodes in the following manner:

∀v i ∈ V , λ i = 1 ∑ j ω ij /w shape + 1 . (4.3.2)
With this choice of Ω, each pixel is biased towards its prior assignment probabilities depending on the relative values of the w ij and w shape : if ω ij ≪ w shape , node v i is influenced more by the prior than by its neighborhood; conversely, if ω ij ≫ w shape , then node v i is more strongly influenced by the label assignment of its neighbors than by the prior probabilities. The corresponding functional is: This model simply penalizes the deviation of vector x s from x s 0 with a uniform weighting scheme.

E s RWconst (x s ) = x s⊤ Lx s + w shape x s -x s 0 2 . (4.3.3)
The prior assignment probabilities can be obtained by computing the empirical mean of the assignment probabilities over the training set T . Since the training T is composed of manually segmented images, the assignment of each pixel is known and the corresponding assignment probability is either 0 or 1. Referring to l k (i) as the assigned label of pixel i in training image k, the binary assignment probability is defined by: (4.3.4) and the corresponding vector is denoted as z s k . Then, x s 0 is computed as:

z s k i =    1 if l k (i) = s, 0 if l k (i) = s,
x s 0 = 1 |T | ∑ k z s k . (4.3.5)

Prior-weighted models

entropy model The previous shape model weights each pixel by the same amount w shape towards the average probability. However, if we own a confidence measure of the quality of this prior model, it could be useful to vary the weights in order to constrain the pixels according to the degree of confidence. For instance, the entropy

e (i) = -∑ s x s 0 i log (x s 0 i ) (4.3.6)
is a measure of how certain the prior model is of the pixel assignment. Indeed, if pixel i was consistently assigned to label s in the training set, then e (i) ≈ 0. Conversely, if pixel i has similar assignment probabilities for all the labels, e (i) ≈ log (|S|) (see figure 4.3.2). Hence, a possible weighting scheme is given by:

Ω ii = w shape 1 - e (i) log |S| . (4.3.7)
gaussian model We observe that the term x sx s 0 2 in the constant shape prior functional can be interpreted as the typical log-likelihood functional for independent Gaussian variables, with mean x s 0 and variance 1:

x s -x s 0 2 2 = -log exp -x s -x s 0 2 . (4.3.8)
Pursuing this Gaussian analogy, we may use the empirical variance as a measure of confidence:

Ω s ii = w shape 1 σ (i) . (4.3.9)
In the previous, we may use the empirical estimate of the variance over the training set:

σ2 (i) = 1 |T | |S| ∑ s ∑ k (z s k i -x s 0 i ) 2 .
(4.3.10) Since x (i, s) is a probability and since the variables in x s are necessarily correlated, such Gaussian modeling can only be a rough approximation but it has the advantage of simplicity. In practice, due to the small number of training examples, it is necessary to regularize the variance term. Indeed, some pixels have a variance that is equal to zero when all training examples agree on this pixel's label. An empirical regularization of the variance can be obtained with the following formula:

Ω s ii = w shape 1 ασ 2 (i) + (1 -α) , (4.3.11) with α ∈ ]0, 1[.

Confidence map

All previous shape models balance the influence of the prior model according to information extracted from the training data. We may also find interesting to reduce the influence of the model according to the strength of the contours in the test image, rather than from the training set: the stronger the contours, the least we should rely on the model (see figure 4.3.2). Assume we possess such a "confidence map" c, with values close to 0 on strong contours, and values close to 1 in homogeneous regions, we can set: Ω s = w shape diag (c) . (4.3.12)

In the following, we propose a simple formula to compute the confidence map, but many other methods would likely be as effective. Denoting the local variance at pixel i computed on a patch with radius r as σ 2 r (i) , we propose the following formula:

c i = exp -k v σ 2 r (i) , (4.3.13)
where k v is a free parameter .

Experimental validation

datasets All our experiments were performed on datasets composed of 3D volumes of the right thigh of healthy subjects, covering a wide range of morphologies (males and females, ages range: 26 to 60), acquired with a 3T Siemens scanner and using 3pt Dixon sequence (TR=10ms, TE1=2.75 ms TE2=3.95 ms TE3=5.15 ms, rf. flip angle =3°; we used the out-of-phase image), with the an average image size of: 220 × 220 × 64px. Four clinically relevant muscles (the quadriceps: vastus lateralis, vastus intermedius, vastus medialis and rectus femoris) were segmented in 30 volumes; in addition, among these 30 volumes, 13 of which had all 13 thigh muscles segmented.

registration For a given test volume, all other volumes served as a training set and were non-rigidly registered onto the test volume. We chose non-rigid registration in order to have the least possible variations of the outer-boundaries of the muscle tissues (the muscle/fat boundaries). The objective is to rely on the RW algorithm only to segment the more difficult inter-muscle boundaries.

As mentionned in the conclusion of the previous chapter, registration of the outerboundaries of the muscle tissues is much easier than registration of the inter-muscle boudaries. Indeed, the latter are often weak or confused with spurious ones, such that trying to register them may lead to larger errors if the registration process is not-rigid enough. For this reason, we generated binary masks of the muscle tissues from the water-map image (see appendix A) and used them to compute the deformation fields. With binary masks, the inter-muscle boundaries are completely ignored by the registration software, and thus only deformed accordingly to the outer-boundary matching. The deformation fields were then applied to the corresponding segmentation atlases in order to generate the average segmentation x 0 .

Registration was achieved using the Drop software [START_REF] Glocker | Dense image registration through MRFs and efficient linear programming[END_REF]. This methods uses a multi-scale discrete modeling approach relying on a Markov random field objective function minimized through efficient linear programming. This method can accept any similarity measure and we chose the default one (sum of absolute differences) to register the binary masks. The software provided a smooth deformation field which we then used to deform the reference atlas.

validation To evaluate the segmentation generated by our methods, we computed Dice coefficients. The expression of the Dice coefficient is:

D = 2 |T ∩ R| / (|T| + |R|),
where T and R are the pixel sets of a given muscle in the inferred segmentation and the reference manual segmentation respectively. Let us denote the binary vector containing the manual segmentation for label s as z s , and the binarized vector containing the inferred segmentation for label s as x s , that is,

x s (i) =    1 if s = arg max t x t (i), 0 otherwise. (4.3.14)
Then, the Dice coefficient for label s of the inferred segmentation x is computed as

D s = 2z s⊤ x s 1 ⊤ (z s + x s ) . (4.3.15)
Thus, if the inferred segmentation is identical to the manual segmentation for a given label, it will have Dice coefficients equal to 1 for all muscles; conversely, if it is entirely different from the manual segmentation, its Dice coefficients will be 0.

Results are presented as boxplots4 , with the Dice coefficients for all the segmented muscles combined into a single boxexcluding the background. Figs. 4.3.3 and 4.3.4show some examples of segmentations obtained with our method (constant weighting scheme, 13-muscle dataset). Large muscles tend to be better segmented that small ones. This is related to the average segmentation x 0 : since the average segmentation and the true segmentation of the target image do not correspond exactly, a larger proportion of the small muscles are incorrectly segmented by x 0 , as compared to the larger muscles. As a result, the assignment probabilities computed by the RW algorithm are generally lower than with large muscles, thus leading to more segmentation errors. This is also true when a muscle has a rarely seen location or shape : large segmentation errors with x 0 may lead to large errors in the resulting segmentation. We conclude that this algorithm only works well for average anatomies. Improvements of this method will necessitate a more flexible model, capable of adapting to a larger range of anatomies. various weighting schemes We tested different weighting schemes based on statistics on the training database: entropy weighting and Gaussian weighting and compared them to the constant weighting scheme. In addition, we also tested the confidence map weighting scheme. Dice coefficient for the different approaches were reported in figure 4.3.5 and show that for all tested parameters w shape , neither weighting scheme improved on the constant weighting; the variance weighting even seems to degrade the performance. The reason for the lack of improvement is possibly due to the limited validity of the statistics on such a small number of training example (29 volumes, the 30th is the one being segmented). However, although no statistically significant differences in terms of Dice coefficients were found, one can observe on the segmented data some visual differences depending on the weighting scheme used. For instance, figure 4.3.6 shows some of the effects of the confidence map scheme on the segmentation results; here, the result is less influenced by the prior shape than with the constant prior, which improves the segmentation quality around the bone; however, the gain in accuracy in one place is lost in another place where the contours are missing, due to the increased freedom of the model. As often, it appears that the compromise between segmenting the contours accurately and keeping a valid shape is hard to find. Since all hyper-parameters were selected by hand -in particular w shape -we expect that an automatized parameter selection would yield more optimal values and maybe find significant differences between the weighting schemes. In the next chapter of this report, we propose such a method for selecting parameters. using the appearance prior As mentioned previously, it is possible to combine several prior terms to benefit from complementary knowledge sources. For instance, let us consider the appearance prior seen in section 4.2.4; this prior term is not usable in itself in the context of muscle separation, since all muscles have the same appearance properties; however, when used in combination with a shape prior term, it could help improving the separation between muscles and fat tissues. In figure 4.3.7, we show identical cross-sections where the left one has been segmented only with the constant shape prior (that is: equation 4.3.3, using w shape = 1.10 -2 ), whereas the right one has been segmented with the same prior term plus an appearance prior (equation 4.2.5, using w app = 1.10 -2 ). One can see the appearance prior has an obvious influence over the muscle/background separation in several places. However, quantitative results -shown in figure 4.3.8 -show that the appearance prior does not improve the segmentation on our dataset. This is due to the fact that, in our test images (see appendix section A.2.2), the average intensity of the fat tissues is close to the average intensity of muscle tissues; moreover, there are notable local variations of the average intensity due to magnetic field inhomogeneities -a notorious artifact in MRI. Although it appears that such an appearance prior is inefficient on our database, the evidence of a discriminative effect shown in figure 4.3.7 leads us to believe it could improve the segmentation results on a different dataset where fat and muscle tissues are more clearly separated (for instance when using T 1 -weighted images).

comparison with previous work We compared the presented method with related methods (see figures 4.3.11 and 4.3.9). For reasons explained in the following, the tests could not be performed on the exact same datasets, but only on overlapping datasets with different number of samples and segmented muscles. Thus, some of these comparisons should be taken cautiously. All tests were performed by the authors of the corresponding methods.

• The method of [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF], briefly presented in section 2.1.2. We recall this method belongs to the deformable model family, using a single reference shape model as initial state and prior term. Dataset: 16 volumes extracted from our 30-muscle dataset, among which one was fully segmented for providing the prior shape; the 15 other volumes were then segmented automatically5 ; we computed the Dice coefficients only on the 4 muscles for which we had the manual segmentations.

• The method of Andrews et al. (2011a), briefly presented in section 2.2.1. This method is a region-based continuous method, with a shape prior derived from PCA performed on training samples projected into a Isometric Log-Ratio space.

Dataset used: our database of 13 volumes with 13 segmented muscles. The images and corresponding segmentations have been cropped and aligned based on the location of the bone, as per the method described in the original paper6 . On this dataset, due to not having a clear difference between the intensity of the muscles and that of the fat tissues, the method did not perform as well as in the published article where the authors present an average global Dice coefficient of 0.92.

• A method called "Registration", which consists of segmenting the volumes using the average segmentation x 0 . This amounts to considering only the performance of the non-rigid registration stage that is applied before constructing the prior model and segmenting the volumes with the RW method. This method can be seen as a multi-atlas segmentation method with majority label voting.

• Our own method for automatic generation of seeds (see chapter 3). Dataset used: our 15-muscle database. Since this method requires only one reference segmentation for providing a shape prior, many different tests could be performed in a cross-validation approach (with 13 volumes, this makes 13 × 12 = 156 tests).

• Our own Random Walks with shape prior method. Dataset used: either our 13muscle or our 4-muscle database. Our first comparison concerns the method of [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF]. We underline the fact that very few volumes of our dataset were segmented using this method, and thus the quantitative comparison with our own method has only limited validity. However some qualitative comments can be made about these results and the cross sections of segmentations obtained with this method (see figure 4.3.10). First, we notice that this method always produces valid shapes, that is, shapes that are realistic and could be actual segmentations. However, due to depending too much on one reference segmentation and also to the local minimization procedure, this method sometimes fails to segment large parts of the volumes even though the contours are well defined, as showed in the figure 4.3.10. These large errors explain the lower figures obtained here.

Second, we compare our method with the method of Andrews et al. (2011a). Here, the obtained Dices on our dataset with the latter method are abnormally low -as compared to the figures produced in the original paper. This can be explained by the fact that the type of contrast of our volumes is ill-adapted to this method, for it relies heavily on a the separation of muscle and fat tissues in a preliminary stage; however, our volumes do not show a clear visual difference between the two types of tissues, as the contrast parameters were rather selected for a good visualization of the contours.

Finally, we compare our method with the automatic seed placement method of chapter 3. As seen in figure 4.3.11, it appears that the autoseed method slightly outperforms7 our current method on the tested dataset (all 13 muscle, 13 samples). It seems that this shape prior using the average segmentation is too constraining and does not [START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF]. This method always produces realistic segmentations, but sometimes fails to detect the contours. allow enough variability to fit the most different shapes. On the other hand, automatic seed placement is less constraining due to having no seeds close to the boundary of the muscles. However, when this latter method fails to label the seed correctly, errors can be extremely large as can be seen from the very low Dice coefficients in figure 4.3.11 per-slice results In practice, certain analyzes of the muscles properties do not require to segment a muscle entirely. In figure 4.3.12 we show that the median Dice coefficient obtained with our method (using constant weighting, as it is the simplest method and provides results that are equivalent to or better than the other weighting schemes) is much greater in the middle slices than in the upper and lower slices. This is due to greater variability of the muscles shapes at their extremities, and to the fact that there are many more missing contours in those places (especially in the upper part of the thigh, where it is extremely difficult to determine the boundaries). These results show that the obtained segmentations are the most accurate in the places most useful for studying the muscles.

pca model

In this section, we propose a segmentation method based on the Random Walks algorithm, in which shape deformation is constrained to remain close to a Principal Component Analysis (PCA) shape space constructed over training examples. Using the PCA allows us to model complex non-rigid shape variations relying on a few eigen-modes. Such formulation does not fit in the probabilistic framework presented in section 4.2.2, and only yields an approximate solution.

Additional notations

Since minimizing 4.1.13 is an independent process for each label, the whole RW process can be equivalently synthesized in one equation. We denote the stacked up labelwise segmentations as x ⊤ = x 1⊤ |x 2⊤ | . . . |x K⊤ , such that x ∈ R KN×1 and the extended Laplacian matrix as:

L =     L • • • 0 . . . L . . . 0 • • • L     ∈ R KN×KN . (4.4.1)
Matrix L is block-diagonal with K identical blocks, whereK is the number of labels. In this notations, the RW objective becomes:

E RW (x) = x T Lx, (4.4.2)

Shape space

It is obvious that the entries of x are not independent from each other, but instead are highly tied to each other, at least for sets of close-by pixels. This implies that there exists an implicit lower-dimensional space in which any true segmentations reside. The principle of a shape space is to design a low-dimensional affine space approximating this implicit space. If we succeed, we expect a valid segmentation to lie "not too far" from the shape space; we also expect the projection of that segmentation onto the shape space to be a good approximation of the segmentation itself. PCA notoriously provides a simple mean to compute such a shape space. Assume we possess a set T of co-registered segmented training volumes. We model a segmentation vector x as a random vector X with normal distribution N (x 0 , Σ), for which we possess a number of samples {z k } k=1...|T | . Let us denote the centered segmentation vectors as x c k = z kx 0 . The expression of the empirical covariance Σ is given by:

Σ ij = 1 |T | ∑ k (z k i -x 0 i ) z k j -x 0 j . (4.4.3)
The previous can be reformulated in matrix form as

Σ = 1 |T | X c X ⊤ c (4.4.4) with X c = x c 1 • • • x c |T | .
Then we compute the eigen-decomposition of Σ:

Σ = U∆U ⊤ , (4.4.5)
where ∆ is a diagonal matrix with the eigenvalues of Σ on the diagonal, and the columns of square matrix U are the eigenvectors of Σ, in the same order as the eigenvalues. Since Σ is symmetric and positive definite, its eigenvectors are orthogonal and its eigenvalues are all positive or zero. Moreover, since the number of training examples is smaller than the size of X: |T | ≪ N, we know we have at most |T | strictly positive eigenvalues, all other eigenvalues being zero. We retain the n < |T | largest eigenvalues and associated eigenvectors, yielding an approximate covariance matrix:

Σ ≈ U n ∆ n U ⊤ n , (4.4.6)
where ∆ n is the n × n diagonal matrix with the n largest eigenvalues of Σ and U n is the non-square matrix whose columns are the n retained eigenvectors.

In practice, due to the large number of variables N, it is not possible to compute Σ. Thus, we use a different approach to perform the PCA. We define the |T | × |T | matrix:

Σ ′ = 1 |T | X ⊤ c X c . (4.4.7)
This matrix being much smaller that Σ, we can perform an eigen-decomposition:

Σ ′ = V∆ ′ V ⊤ (4.4.8)
with obvious notations. Since X c Σ ′ = ΣX c , (4.4.9) the eigenvalues of Σ ′ are all eigenvalues of Σ:

X c Σ ′ V = X c V∆ ′ , (4.4.10) Σ (X c V) = (X c V) ∆ ′ . (4.4.11)
The eigenvectors of Σ are now obtained by computing

U n = X c V n , (4.4.12)
which will be normalized afterward to form an orthonormal basis.

The matrix U n determines a shape space, in which the projection of any segmentation x in the shape space is given by

x = x 0 + U n γ, (4.4.13)
where γ is the coordinate vector of x in the shape space. We expect that x will a good approximation of a valid segmentation x. For any segmentation, we can write:

x = x 0 + U n γ + dx, (4.4.14)
where dx is the deviation of x from the shape space. If the shape space models the space of valid segmentation well, then the norm of dx will be small when x is a valid segmentation.

PCA model

In order to obtain a segmentation which remains close to the shape space, we want to minimize the objective function 4.4.2 with respect to both dx and γ, while keeping dx small. This leads to the following functional: where w dx is a hyper parameter setting a constraint on the norm of dx. We reformulate 4.4.15 as:

E RWpca (dx, γ) = (dx + Uγ + x) T L (dx + Uγ + x) + w dx dx 2 , ( 4
E RWpca (y) = (Ay + x 0 ) T L (Ay + x 0 ) + w dx y T By, (4.4.16)
with

y = dx γ , A = [I KN U] , B = I KN 0 0 0 . (4.4.17)
where I KN is the identity matrix of size KN × KN.

The minimum of 4.4.16 verifies:

A T LA + w dx B y = -A T Lx 0 . (4.4.18)
Obviously, the entries of x being probability distributions, assuming a Gaussian distribution can only lead to a rough approximation. For these reason, we expect the PCA shape space not to be a very good space of realistic solutions. In particular, it is unlikely that any projection into this shape space will yield a probability distributions.

Experimental validation

Our experiments were conducted on the same dataset as with the previous method (see section 4.3.4). This dataset consists of 30 volumes, in which four clinically relevant muscles (the quadriceps: vastus lateralis, vastus intermedius, vastus medialis and rectus femoris) were segmented in 30 volumes. Comparisons are made using Dice coefficients as previously.

In figure 4.4.1, we compare the Dice coefficients obtained with our shape space method to the previously presented simpler method of section 4.3 (using a constant scheme), and to the baseline registration method. It appears that the presented method does not improve over the previous, simpler method. This can be explained by the poor quality of the PCA generative model for probabilistic segmentations.

On figure 4.4.2, we show the effect of the PCA shape prior. We see that the PCA model is deformed from its initial shape (left column) to fit the boundaries of the test image (middle column). Due to the rough approximation of modeling probabilistic atlases with a linear subspace, the boundaries of the muscles are rather fuzzy. However, when allowing a small deviation (called dx) from the model, as does our method, the contours are more accurately defined (right column). Notice there remains segmentation errors as the PCA model is not capable of deforming enough to fit unusual shapes while retaining a realistic topology. In the presented method, the main segmentation errors are due to the muscles with unusual shape -more elongated, smaller, etc. Modeling probabilistic segmentations with PCA model does not allow representing shapes which differ too much from standard shapes.

conclusion

To conclude this chapter, we will say that directly incorporating prior knowledge of shape into the Random Walks framework is valid and compares favorably to state-ofthe-art methods, but also shows some limitations. In the case of the simplest prior (sec. 4.3), the simplicity and the validity of the formulation allow to propose an easyto-use method, quite efficient (less than 30 seconds for 13 muscles in a 220 × 220 × 64 volume) and which does not require a large training dataset (we used between 12 and 29 samples during our experiments). However, errors inevitably appear when the average model x 0 is too different from the segmented image. In particular, smaller muscles with too large inter-subject variations are likely to be incorrectly registered. In a very recent publication, [START_REF] Eslami | Segmentation by retrieval with guided random walks: Application to left ventricle segmentation in mri[END_REF] have proposed a very similar shape prior formulation for the Random Walks framework, called "guided" Random Walks. However, instead of using the empirical mean over the training samples, the author propose to use directly each training segmentation as a reference segmentation x 0 and retain the infered segmentation x that most overlaps with x 0 . This promising approach is aimed at retrieving shapes which are quite different from the average model. For lack of time, we were not able to test this method on our data.

Ideally, we would like to handle large scale deformations by allowing the model to evolve in a low-dimensional shape space of valid segmentations. Small-scale deformations or subject-specific fine details would be segmented like in the simple model through the RW process. The PCA model presented in section 4.4 was an attempt at such flexible shape prior. Indeed, PCA is a well-know way to generate a low- dimensional shape space. The mitigated performance we achieved is likely to be due to the weak modeling of probabilities with a PCA model. Future direction would include finding a different shape space more compatible with probabilities, such as a barycenter model -indeed, a barycenter of a probability distribution remains a probability distribution. In the conclusion of this report, we also propose a different possible direction for dealing with large scale deformations, relying on additional connections between remote voxels.

L E A R N I N G PA R A M E T E R S F O R R A N D O M WA L K S -B A S E D S E G M E N TAT I O N
In the previous chapter, we empirically observed that the accuracy of the RW algorithm relies heavily on the relative weighting between the various contrast and prior terms. For instance, the contrast parameter used to construct the Laplacian matrix was chosen so that it gives the best results on average. In a more general approach, a different edge weighting function could be designed to extract different features of the image, more complex and specific appearance information. Then, instead of selecting the best Laplacian, it would be interesting to see if a particular linear combination of these Laplacians would yield better results than a single one of them. Similarly, we could also combine the different prior models.

At present, the proposed methods rely on a user to hand-tune the parameters or on exhaustive cross-validation. However, both these approaches quickly become infeasible as the number of terms in the RW objective function increase: a human would be unable to hand-tune a large number of parameters, and a thorough cross-validation would require vast and expensive computational resources.

In contrast to the RW literature, the problem of parameter estimation has received considerable attention in the case of discrete models such as conditional random fields [START_REF] Lafferty | Conditional random fields: Probabilistic models for segmenting and labeling sequence data[END_REF]. Recent years have witnessed the emergence of structured-output support vector machine (Structured SVM) as one of the most effective discriminative frameworks for supervised parameter estimation [START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF]. Given a training dataset that consists of pairs of input and their groundtruth output, structured SVM minimizes the empirical risk of the inferred output with respect to the ground-truth output. The risk is determined by a user-specified loss function that measures the difference in quality between two given outputs. The parameters are estimated by minimizing a convex quadratic program, for which many efficient algorithms have been proposed (Joachims et al., 2009;[START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF][START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF]. The structured SVM formulation has been successfully employed for many applications including graph-cuts based segmentation [START_REF] Szummer | Learning CRFs using graph cuts[END_REF], object detection [START_REF] Blaschko | Learning to localize objects with structured output regression[END_REF] and natural language parsing [START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF].

Inspired by the efficacy of structured SVM in discrete models, we would like to discriminatively learn the parameters of the RW formulation. Such an approach would allow us to fully exploit the benefits of RW using numerous contrast and prior terms. To this end, a straightforward application of structured SVM would require a training dataset that consists of pairs of inputs -here, medical acquisitions using MRI scans -as well as their ground-truth outputs -in our case, the optimal probabilistic segmentation. In other words, we require a human to provide us with the output of the RW algorithm for the best set of parameters. This is an unreasonable demand since the knowledge of the optimal probabilistic segmentation is as difficult to acquire as it is to hand-tune the parameters itself. Thus we are faced with a scenario where we cannot directly use structured SVM to estimate the desired parameters.

In order to handle the above difficulty, we propose a novel formulation for discriminative parameter estimation in the RW framework. Specifically, we learn the parameters using a weakly supervised dataset that consists of pairs of medical acquisitions and their hard segmentations. Unlike probabilistic segmentations, hard segmentations can be obtained easily from human annotators. We treat the optimal probabilistic segmentation that is compatible with the hard segmentation as a latent variable. Here, compatibility refers to the fact that the probability of the ground-truth label (as specified by the hard segmentation) should be greater than the probability of all other labels for each pixel/voxel. Clearly, the number of compatible probabilistic, or soft, segmentations for a given hard segmentation (that is, the size of the space for the latent variables) is uncountably infinite. Nonetheless, the resulting representation allows us to learn the parameters using the latent SVM formulation (Yu and Joachims,

2009).

While latent SVM does not result in a convex optimization problem, its local optimum solution can be obtained using the iterative concave-convex procedure (CCCP) [START_REF] Yuille | The concave-convex procedure (cccp)[END_REF]. At each iteration, CCCP performs two steps: (i) estimating a compatible probabilistic segmentation for each training sample using the current set of parameters -commonly referred to as annotation consistent inference (ACI); and (ii) updating the parameters by fixing the compatible probabilistic segmentations to those obtained during ACI. The second step of CCCP involves solving a structured SVM problem, which lends itself to efficient optimization. In order to make the overall algorithm computationally feasible, we propose a novel efficient approach for ACI based on dual decomposition [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF] Komodakis | MRF Optimization via Dual Decomposition: Message-Passing Revisited[END_REF]. We demonstrate the benefit of our learning framework over a baseline structured SVM using a challenging dataset of real 3D MRI volumes1 .

related work

As mentioned earlier, the methods that employ the RW algorithm rely on parameters that have been hand-tuned or exhaustively cross-validated (see chapter 4). In contrast, the parameters of discrete models such as conditional random fields (CRFs) are often estimated using large supervised training datasets. One of the most effective supervised learning formulations is structured SVM [START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF], which generalizes the well-known binary SVM classifier to cases where the de-sired output is a vector of inter-dependent elements. Structured SVM minimizes a regularized upper bound on the empirical risk as defined by a user-specified loss function. The loss function provides a measure of the difference between two different outputs. Using the structured SVM formulation, the problem of parameter estimation reduces to that of solving a convex quadratic problem, for which several efficient algorithms have been proposed in the literature (Joachims et al., 2009;[START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF][START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF]. The strong theoretical foundation and practical feasibility of structured SVM has led to its use in several real-world applications [START_REF] Blaschko | Learning to localize objects with structured output regression[END_REF][START_REF] Szummer | Learning CRFs using graph cuts[END_REF][START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF]. One of the main deficiencies of structured SVM is that it requires a fully supervised training dataset, which consists of pairs of inputs and their desired outputs. Often, it is extremely challenging, or even impossible, to obtain full supervision (for example, in our case where a human annotator cannot specify the optimal probabilistic segmentation). However, it is relatively easy to obtain weakly supervised training samples.

In order to exploit the information present in such datasets, several researchers have independently proposed an extension to structured SVM known as latent SVM [START_REF] Felzenszwalb | A discriminatively trained, multiscale, deformable part model[END_REF][START_REF] Smola | Kernel methods for missing variables[END_REF]Yu and Joachims, 2009). Latent SVM treats any missing information in the human annotation as latent variables, which allows it to upper bound the empirical risk as a difference-of-convex optimization problem. Although its globally optimal solution cannot in general be computed in polynomial time, its local minimum or saddle point solution can be obtained efficiently using the CCCP algorithm (Yu and Joachims, 2009;[START_REF] Yuille | The concave-convex procedure (cccp)[END_REF]. Since its first appearance, latent SVM has been steadily gaining popularity due to its accurate empirical performance [START_REF] Felzenszwalb | A discriminatively trained, multiscale, deformable part model[END_REF][START_REF] Kumar | Self-paced learning for latent variable models[END_REF][START_REF] Kumar | Learning specific-class segmentation from diverse data[END_REF]Yu and Joachims, 2009). Most of the previous instantiations of structured SVM and latent SVM have focused on discrete CRFs, that is, where the desired output is deterministic. In contrast, we propose a novel approach to learn the parameters of the probabilistic RW algorithm, where the latent variables model the unknown optimal soft segmentation that is compatible with a given hard segmentation.

preliminaries

Notations

In this part, the notations will differ slightly from chapter 4, for consistency with the SVM notations. We will refer to a 3D volume as x. We use i to denote the index of a voxel in this volume, and denote the set of all indices as V. In a hard segmentation, each voxel is assigned a label s ∈ S (for example, the index of a muscle). We will use z to represent the human annotation (that is, the class labels of the voxels in x) in binary form:

z (i, s) =    1 if voxel i is of class s, 0 otherwise.
(5.2.1)

In other words, the binary form z of the annotation specifies delta distribution over the putative labels for each voxel. Our training dataset is a collection of training images x and hard segmentations z: D = {(x k , z k )} k . Note that we use subscript k to denote the input index within a dataset, and parenthetical i to denote a voxel within a particular input.

We recall that the RW algorithm provides a probabilistic -or soft -segmentation of an input x, which we denote by y, that is,

y (i, s) = Pr [voxel i is of class s] , ∀i ∈ V , s ∈ S.
(5.2.2)

When using one contrast term and one prior model, the RW algorithm amounts to minimizing the following convex quadratic objective functional:

E (y, x) = y ⊤ L (x) y + w prior y -y 0 2 Ω 0 (x) , (5.2.3) = y ⊤ L (x) y + E prior (y, x) . (5.2.4)
Here, y 0 is a reference prior probabilistic segmentation dependent on appearance [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF] or shape (Baudin et al., 2012a), and Ω 0 (x) is a diagonal matrix that specifies a voxel-wise weighting scheme for x. It can be shown that the energy functional 5.2.4 is convex and can be minimized over all possible probabilistic segmentations y in polynomial time. In fact, [START_REF] Grady | Random walks for image segmentation[END_REF] showed that the RW algorithm amounts to solving a sparse linear system of equations, which lends itself to efficient optimization. We refer the reader to section 4.2.2 and to the work of [START_REF] Grady | Random walks for image segmentation[END_REF] for further details.

Parameters and Feature Vectors

In the above description of the RW algorithm, we restricted ourselves to a single Laplacian and a single prior. However, our goal is to enable the use of numerous Laplacians and priors. To this end, let {L α } α denote a known family of Laplacian matrices and E β (•) β denote a known family of prior energy functionals. In section 5.4, we will specify the family of Laplacians and priors used in our experiments. We denote the general form of a linear combination of Laplacians and prior terms as:

L (x; w) = ∑ α w α L α (x) , w α ≥ 0, ∀α, (5.2.5) E prior (•, x; w) = ∑ β w β E β (•, x) , w β ≥ 0, ∀β.
(5.2.6)

Each term E β (•, x) is of the form:

E β (y, x) = y -y β 2 Ω β (x) , (5.2.7)
where y β is the β-th reference segmentation and Ω β (x) is the corresponding voxel- wised weighting matrix (which are both known). We denote the set of all parameters as w = w α , w β . Clearly the RW energy 5.2.4 is linear in w. and can therefore be formulated as:

E (y, x; w) = y ⊤ L (x; w) y + E prior , (5.2.8)

= w ⊤ ψ (x, y) , (5.2.9)

where ψ (x, y) is known as the joint feature vector of x and y. Note that by restrict- ing the parameters to be non-negative (that is, w ≥ 0), we ensure that the energy functional E (•, x; w) remains convex.

Loss Function

As mentioned earlier, we would like to estimate the parameters w by minimizing the empirical risk over the training samples. The risk is specified using a loss function that measures the difference between two segmentations. In this work, we define the loss function as the number of incorrectly labeled voxels. Formally, let bin (y) denote the underlying hard segmentation of the soft segmentation y, that is bin (y) (i, s) = δ (s = arg max s∈S y (i, s)), where δ is the Kronecker function. The loss function is defined as

∆ (z, y) = 1 - 1 |V | bin (y) ⊤ z, (5.2.10)
where V is the set of all voxels and |•| denotes the cardinality of a set. If both segmen- tations are equal, then ∆ (z, y) = 0; if they are disjoint, ∆ (z, y) = 1 . Note for later that this function in not concave in y, and that infinitely many soft segmentations y are compatible with the hard segmentation z, that is, there is an infinite number of solutions y to the equation ∆ (z, y) = 0.

Using the upper bound on the empirical risk, the latent SVM optimization problem is specified as follows:

min w, ξ k ≥0 λ w 2 + λ ′ w -w 0 2 + 1 N ∑ k ξ k , (5.3.9) s.t. min ∆(x k , ŷ)=0 w ⊤ ψ (x k , ŷ) ≤ w ⊤ ψ (x k , ȳ) -∆ (z k , ȳ) + ξ k , ∀ ȳ, ∀k,
where the slack variable ξ k represents the upper bound of the risk for the k-th training sample. Note that we have added two regularization terms for the parameters w. The first term w 2 , weighted by hyper-parameter λ, ensures that we do not over-fit to the training samples. The second ww 0 2 , weighted by hyper-parameter λ ′ , ensures that we do not obtain a solution that is very far away from our initial estimate w 0 . The reason for including this term is that our upper bound to the empirical risk may not be sufficiently tight. Thus, if we do not encourage our solution to lie close to the initial estimate, it may drift towards an inaccurate set of parameters. In section 5.4, we show the empirical effect of the hyper-parameters λ and λ ′ on the accuracy of the parameters.

A local minimum or saddle point solution for problem 5.3.9 can be found using the iterative CCCP method, which is outlined in Algorithm 5.1. It consists of two main steps at each iteration: (i) step 3, which involves estimating a compatible soft segmentation for each training sample -known as annotation consistent inference (ACI); and (ii) step 4, which involves updating the parameters by solving problem 5.3.10. In the following subsections, we provide efficient algorithms for both the steps.

Algorithm 5.1 The CCCP method for parameter estimation using latent SVM. Input: Dataset D, λ, λ ′ , w 0 ,ε 1: Set t = 0. Initialize w t = w 0 .

2: repeat 3:

Compute y ⋆ k = arg min ŷk , ∆(z k , ŷk )=0 w ⊤ t ψ (x k , ŷk ) , ∀k.

4:

Update the parameters by solving the following problem

w t+1 = arg min w, ξ k ≥0 λ w 2 + λ ′ w -w 0 2 + 1 N ∑ k ξ k , (5.3.10) s.t. w ⊤ ψ (x k , y ⋆ k ) ≤ w ⊤ ψ (x k , ȳ) -∆ (z k , ȳ) + ξ k , ∀ ȳ, ∀k, 5: t = t + 1 6:
until The objective function of problem 5.3.9 does not decrease below tolerance ε.

where ρ t m are the dual variables satisfying ∑ m ρ t m = 0. The dual variables are initialized as ρ 0 m = 0, ∀m, and updated at iteration t as follows:

ρ t+1 m ← η t y ⋆ m -1 M M ∑ n=1 y ⋆ n ,
where η t is the learning rate at iteration t. Under fairly general conditions, this iterative strategy converges to the globally optimal solution of the original problem, that is,

y ⋆ = y ⋆ m , ∀m.
We refer the reader to the works of Bertsekas (1999); Komodakis et al. (2007) for details.

In order to specify our slave problems, we divide the set of voxels V into subsets V m , m = 1, . . . , M, such that each pair of neighboring voxels (i, j) ∈ N appear together in exactly one subset V m . Given such a division of voxels, our slave problems correspond to the following: min (5.3.17) where L m (x; w) is the Laplacian matrix whose non-zero entries correspond to the voxels in V m , which leads to

y m ∈C(V m ) y ⊤ m L m (x; w) y m + E prior m (y m , x; w) + ρ t m y m ,
L = ∑ m L m . (5.3.18)
The prior prior energy functions E prior m modifies the original prior E prior by weighting each voxel i ∈ V m by the reciprocal number of subsets that contain i. This weighting scheme allows us to decompose the prior term exactly, that is

E prior (y, x; w) = ∑ m E prior m (y, x; w) , (5.3.19)
which is the consequence of the following trivial equality: .3.20) where R is the diagonal matrix such that R ii = 1/ |{V n , s.t. i ∈ V n }|. In other words, the prior term for each voxel i ∈ V is multiplied by the number of times voxel i appears in a slave problem V n . Thus, the slave problems defined above provide a valid reparametrization of the original problem 5.3.11, that is g (y) = ∑ m g m (y). Decomposing the set of constraints C (V m ) for each subset of variables is straightforward. By using small subsets V m we can optimize each slave problem in every iteration independently using a standard quadratic programming solver. The sub-problems -which are quadratic programscan be made sufficiently small to be solved in parallel with a dedicated solver. In our experiments, we used the Mosek solver. To the best of our knowledge, this is the first application of dual decomposition to solve a probabilistic segmentation problem under linear constraints.

y -y β 2 Ω β = ∑ m y -y β 2 RΩ β , ( 5

Parameter Update

We now describe how the parameters can be efficiently updated by solving problem 5.3.10 for a fixed set of soft segmentations y ⋆ k . Note that while problem 5.3.10 is convex, the number of constraints is of the order of the number of possible soft segmentations of an input. In other words, it consists of an infinite number of constraints. Nonetheless, the problem can be solved efficiently using the cutting plane method outlined in Algorithm 5.2. The method starts by specifying no constraints for any of the training samples. At each iteration, it finds the most violated constraint for each sample (step 3), and updates the parameters (step 5) until the increase in the objective function is less than the tolerance.

Algorithm 5.2 The cutting plane method for updating the parameters.

Input: Dataset D, y ⋆ k , λ, λ ′ , w t , ε 1: Set W k = ∅ for all k. Initialize w = w t . 2: repeat 3:
Compute ȳk = arg min y w ⊤ ψ (x k , y) -∆ (z k , y) , ∀k.

4:

Update W k = W k y k .

5:

Update the parameters by solving the following problem:

min w, ξ k ≥0 λ||w|| 2 + λ ′ ||w -w 0 || 2 + 1 N ∑ k ξ k , (5.3.21) s.t. w ⊤ [ψ (x k , y ⋆ k ) -ψ (x k , y)] ≤ ξ k -∆ (z k , y) , ∀y ∈ W k , ∀k .
6: until The objective function of problem 5.3.21 does not increase above tolerance ε.

The main ingredient of the cutting plane method is the computation of the most violated constraint, that is,

y = arg min y w ⊤ ψ (x, y) -∆ (z, y) (5.3.22)
for a given sample (x, y). While in most of the previous instantiations of structured SVM and latent SVM, the above problem can be solved optimality, in our case it presents a computational challenge. The reason for this is that the loss function, defined in equation 5.2.10, is not concave in y, which makes the above problem non-convex. One obvious solution to this challenge would be to replace the loss function by its concave approximation, that is,

∆ ′ (z, y) = 1 - 1 |V | y ⊤ z.
(5.3.23)

However, the resulting convex approximation of problem 5.3.10 would not lend itself to the RW algorithm, as it does not guaranty that the output is a probability distribution as required by the RW probabilistic framework. In other words, we would need to resort to a dual decomposition strategy to obtain a constraint at each iteration of the cutting plane method. Such an approach would render the overall CCCP algorithm computationally infeasible. Instead, in this work, we obtain the approximate most violated constraint as the predicted segmentation, that is, y = arg min y w ⊤ ψ (x, y) .

(5.3.24)

The above problem only requires the minimization of the energy (without any loss terms), which can be performed efficiently using the RW algorithm. As our results indicate, coupled with the regularization of the parameters, this simple cutting plane strategy provides an accurate set of parameters.

experimental validation

We demonstrate the efficacy of our approach on a challenging dataset. Specifically, we test the hypothesis that it is important to model the unknown optimal soft segmentations of the training samples using latent variables in order to learn an accurate set of parameters. In what follows, we describe the experimental setup and our results in detail.

Dataset

We use the 30 volume dataset described in section 5.4. This datasets consists of volumes of the thigh region of dimensions 224 × 224 × 100. The various segments correspond to 4 different muscle groups together with the background class. We randomly split the dataset into 80% for training and 20% for testing. In order to reduce the training time for both our method and the baselines, we divide each volume into 100/2 volumes of dimension 224 × 224 × 2. During dual decomposition, we further divide the volumes into dimensions 10 × 10 × 2 in order to specify the slave problems.

Laplacians and Prior Terms

We used 4 different Laplacians based on the standard weighting function given in section 4.1, 3 of them corresponding to β ∈ {10, 50, 100} -to be used with volumes which have been normalized with respect to their standard deviation. The last one to the alternate expression given in footnote 2 with β = 100. Furthermore, we use two shape priors and one appearance prior: the first one as the shape prior with constant weighting as described in section 4.3, the second one using the entropy weighting also described in section 4.3, and the third one being the contrast based prior of [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF]. This results in a total of 7 parameters to be estimated.

Methods

The main hypothesis of our work is that it is important to represent the unknown optimal soft segmentation using latent variables. In order to test this hypothesis, we compare our method with a baseline structured SVM that replaces the latent variables with the given hard segmentations. In other words, our baseline estimates the parameters by solving the following optimization problem:

min w, ξ k ≥0 λ||w|| 2 + λ ′ w -w 0 2 + 1 N N ∑ k=1 ξ k , (5.4.1) s.t. w ⊤ [ψ (x, z k ) -ψ (x, y)] ≤ ξ k -∆ (z, y) , ∀ ȳ, ∀k,
which is the same as 5.3.8 with an added proximal regularization term. The above problem can be solved using Algorithm 5.2, where the imputed soft segmentations y ⋆ k are replaced by the hard segmentations z k .

During our experiments, we found that replacing the hard segmentation with a pseudo soft segmentation based on the distance transform systematically decreased the loss of the output. Thus, the method referred to as "Baseline" uses a structured SVM with distance-transform "softened" segmentations. The transformed segmentation, denoted as zk , is computed in the following manner:

zk (i, s) = exp -γdist k (i; s) 2 ∑ t exp -γdist k (i; t) 2 , (5.4.2)
where dist k (i, s) is the Euclidean distance from pixel i to the closest pixel with label s in the hard segmentation z k ; if pixel i already belongs to class s, then dist k (i; s) = 0. This term is normalized by a sum over all labels such that zk is a probability distribution. Parameter γ > 0 determines the amount of spatial smoothing of the distribution, that is, the smaller γ is, the more progressive are the transitions between labels. Fig. 5.4.1 illustrates the softening process of the segmentation.

Results

Figure 5.4.2 shows the test loss for three different methods: (i) the initial hand-tuned parameters w 0 ; (ii) the baseline structured SVM with distance transforms; and (iii) our proposed approach using latent SVM. As can be seen from Fig. 5.4.2, latent SVM provides significantly better results than the baselines -even when using the distance transform (which is the baseline method displayed in Fig. 5.4.2).

For the 4 x 5 hyper-parameter settings that we report (that is, four different values of λ and 5 different values of λ ′ ), latent SVM is significantly better than SVM in 15 cases, and significantly worse in only 2 cases. Note that latent SVM provides the best results for very small values of λ ′ , which indicates that the upper bound on the empirical risk in tight. As expected, for sufficiently large values of λ ′ , all the methods provide similar results. For the best settings of the corresponding hyper-parameters, the percentage of incorrectly labeled voxels as follows: (i) for w 0 , 13.5%; (ii) for structured SVM, 10.0%; and (iii) for latent SVM, 9.2%. Fig. 5.4.3 shows some example segmentations for the various methods.

Conclusion

Results show that our approach allows to estimate good parameters for the Random Walks segmentation method using shape priors. At the moment, since tests were performed on a small dataset (6 test volumes, 24 training volumes, 4 muscles), we can only claim it to be a proof of concept. Our results provide empirical evidence for our main hypothesis, namely the importance of treating the underlying soft segmentations as latent variables. This can be seen not only in the overall improvement achieved λ -that is, the SVM regularization parameter. Curves are plotted with respect of λ ′ -that is, the proximal regularization parameter, which penalizes deviations from the hand-tuned parameters w 0 . "Baseline" corresponds to the standard struct-SVM method using distance-transformed pseudo-soft segmentations zk as annotation vectors.

by our approach, but also by the relative robustness of latent SVM to the choice of hyperparameters. This allows us to envision more extensive estimation processes. By carefully choosing more Laplacians and shape priors and learning their relative weighting, our framework may be able to provide bigger improvements. In particular, we are interested in learning separate parameters for different image regions. Indeed, some regions have edges which are consistently visible, while others do not. Thus it would be useful to change the balance between the model and the contrast term based on this knowledge, a task that would be impossible to fulfill manually. Moreover the latent variable formulation opens the door for more semi-supervised learning: authorizing the use of partially or coarsely segmented ones would allow us to increase the size of the learning set by orders of magnitude. An improvement of the method could be tested in replacing the implied L 2 norm of the regularization termsw 2 and ww 0 2 -by a L 1 norm. Such change could promote sparsity in the parameter values, such that only the most relevent terms are selected, rather that searching for a combination of all terms. However, our experiments revealed an important drawback: computing time. Indeed SVM learning is an iterative process and in each iteration, all the training samples have to be segmented (during the LAI process). Each segmentation task is taking several dozen of seconds for the full image on a powerful desktop computer. Even more problematic is the ACI stage, where we had to solve a large quadratic program using dual decomposition. Even though it allowed us to estimate the latent soft segmentation, this process takes several dozens of minutes for a single volume. All our experiments were conducted on the computer cluster of the École Centrale Paris. Even there, a single experiment; for a given set of parameters λ, λ ′ , etc., took more than 2 hours for estimating seven parameters. However, we point out that the code written in Python and not optimized for efficient computation. Besides, all presented algorithms (cutting plance, CCCP, ACI) are highly parallelizable.

A second limitation is the apparent weakness of the improvements in the segmentation results (as measured by loss decrease). The limited improvements may be due to the small number of parameters we estimated and by the fact that there were little room for improvement with the segmentation method considered. As always, parameter learning can only improve the results of a method within the limits of the method.

Here, the single reference shape model may not be able to provide results which are much better than what was obtained by hand-tuning the parameters. To conclude, we say that this learning method appears to be effectively learning and could be re-used and extended to more interesting cases, but we also believe that the segmentation model is even more crucial and should be improved as well.

material

Most of the presented work was coded using the open source and highly portable Python programming language, in combination with the scientific libraries Numpy and Scipy. The author would like to contribute to publicize to the scientific community the great versatility and power of this free software. All development and computations were performed on standard desktop machines, except for the learning process presented in chapter 5 which had to be performed on the computer cluster of École Centrale Paris.

All MRI volumes used for testing and training the segmentation methods were acquired on the 3 Tesla Siemens scanner of the Institute of Myology.

Moreover, all manual segmentations where performed using the free and open source ITK Snap software [START_REF] Yushkevich | User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability[END_REF].

summary of the contributions

In this thesis report, we have presented several contributions towards solving the difficult problem of muscle MRI segmentation. All presented methods require a require a registration stage -either rigid or non-rigid -of a training set of atlases onto the target image.

Our first contribution takes advantage of the good performance of interactive Random Walks segmentation algorithm [START_REF] Grady | Random walks for image segmentation[END_REF] by automatically computing the seeds, which are normally drawn manually. This approach consists in generating unlabeled seeds over the target image, and labeling them automatically by detecting the contours between the muscles and based on prior knowledge of the relative position of the muscles extracted from one reference atlas. Labeling is achieved through a graph-based discrete energy approximate minimization procedure.

In a different approach, we proposed to integrate prior-knowledge shape directly to the Random Walks formulation. In a first contribution, similar to an intensity prior formulation [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF], we added a term to the Random Walks objective that prevents the segmentation to deviate too much from a known average segmentation. We show this shape prior term does not break the probabilistic framework of the RW formulation with an easy-to-minimize objective functional -which amounts to solving a large sparse linear system. We propose several weighting schemes for the shape term based on reasonable intuitions. Obtained results showed that our method is favorably comparable to the current state-of-the-art. In a second approach, we try to expand the range of accessible shapes in order to improve the segmentation performance on morphologies far from the average. We build a low-dimensional shape space through Principal Component Analysis on the training set. Then we propose an objective in which the segmentation is constrained to remain in the vicinity of the shape space while it minimizes the RW objective. However, the results indicate that this approach suffers from the poor modelization of probabilistic segmentations using a PCA as well as the small number of samples in the training set.

In a last part, we proposed a novel discriminative learning framework to estimate the parameters for the probabilistic Random Walks segmentation algorithm. Indeed, the information in the training dataset had not been exploited to finely tune the proposed models up to this point. In order to deal with the lack of full supervision, we represented the optimal soft segmentation that is compatible with the hard segmentation of each training sample as a latent variable. This allowed us to formulate the problem of parameter estimation using latent SVM, which upper bounds the empirical risk of prediction with a difference of convex optimization program. Using a challenging clinical dataset of MRI volumes, we demonstrated the efficacy of our approach over the baseline method that replaces the latent variables with the given hard segmentations. Our formulation allows for a large number of parameters to be estimated using a training dataset.

future directions

Although it appears that the automatic seed method presented in chapter 3 slightly outperforms the Random Walks-based model of chapter 4, we believe that the latter method has more potential than the former for several reasons: (i) the automatic seeding method requires one more computing stage than the RW-based method (registration, seeding, segmentation for the former, registration and segmentation for the latter) and is more likely to propagate errors from one stage to the next. Integrating the shape prior within the segmentation method as does our RW-based method in a very simple way is likely to be more robust -which the smaller number of outliers for this method seems to confirm (in Fig. 4.3.11). (ii) Many possible extensions of the RW-based method can be explored, for which we propose two: a. Incorporate long-range connections between pixels instead of using only adjacent edges -this would help in situations where the prior shape model is too different from the target image. This could be achieved in the following manner: generate random edges between voxels, targeting an optimal edge length based on a small world perspective [START_REF] Kleinberg | The small-world phenomenon: an algorithm perspective[END_REF]. On the training set, learn typical distributions of the voxel intensity along the edges, in order to build a contour detector. The geodesic distance or the commute-time distance can be used as features to the classifier. Run the classifier on the tested image to gener-ate the weights for the long-range edges: voxels with no contours between them would have a transition probability close to one. The added edges do not change the formulation of the Random Walks algorithm. In addition, if one wishes to increase the probability of label change when a contour is found between the connected voxels, it can also be done by introducing non-zero inter-label weights in the global Laplacian (that is, the Laplacian which gather all RW equations, as seen in section 4.4). It can be shown that this formulation retains a probabilistic meaning, and therefore also resolves in a large sparse linear system.

b. Instead of a PCA (see sec. 4.4), build a shape space from a barycenter of training segmentations -a formulation which may better fit within the probabilistic framework since a barycenter of probability distributions is still a probability distribution.

(iii) Apart form shape priors, different types of priors can be easily integrated to the RW-based method, in particular appearance-based priors involving intensity [START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF].

Regarding parameter estimation, as mentioned in chapter 5, future work could involve exploring extensions where the parameters specify spatially varying terms in the RW objective function. This would allow us to encourage the presence of edges in between different segments, while at the same time not penalize the absence of edge in the background. The latent SVM framework can also be used to estimate parameters from partially segmented data. Such an approach would allow us to scale the size of the training dataset by orders of magnitude.

Although their is plenty of room for improvement in automatic techniques, one should not loose sight that segmentation methods are to be used by researchers and clinicians in practical settings. In particular, considering no fully automatic method will provide perfect results in a near future, one should still consider user interaction as an acceptable mean. Thus, it is important that segmentation methods remain simple to use, and that learning methods remain practical, the latter of which the presented work did not fully achieve. More importantly, it is crucial to develop new interactive methods for correcting the imperfect output of automatic methods. Otherwise, they may remain useless if it takes as much time to correct an automatic segmentation than to segment the image manually from the start. The aim of this section is to give a brief insight into the origin of the contrast in MRI. It is based on the basics of Magnetic Resonance Imaging by [START_REF] Fleckenstein | Muscle imaging in health and disease[END_REF]. Nuclear Magnetic Resonance is a quantum phenomenon affecting certain atomic nuclei having a non-zero nuclear spin -→ S . The spin is an inherent property of protons and neutrons, and of nuclei with an odd number of these particles. Some of the nuclei which have a non-zero spin can be found in the human body, such as: 1 H, 13 C, 23 Na, 31 P. Among them, Hydrogen is the most prevalent and therefore is the most adequate for MRI.

The nucleus of a 1 H hydrogen atom consists of a single proton, which has a positive electric charge. The spin of the proton generates an electrical current loop and hence a magnetic field. The proton has a magnetic moment parallel to the direction of its spin:

-→ µ = γ -→ S , where γ is the gyromagnetic ratio of the nucleus of interest (γ H = 267.513 × 10 6 rad.s -1 .T -1 ). The gyromagnetic ratio γ is specific for each nuclear species.

Without a magnetic field, the orientation of the magnetic moment is random. When placed in a magnetic field -→ B 0 , -→ µ tends to align with the direction of B 0 and, due to a resonance phenomenon, begins a precessional movement of the spinning axis around the direction of B 0 . The angular velocity ω 0 -called Larmor angular velocity -of the precession movement is proportional to the amplitude of B 0 and the gyromagnetic ratio γ: ω 0 = γB 0 . (A.1.1) When a sample containing many protons is placed in a magnetic field -→ B 0 , the individual protons tend to align with the direction of the field which results in a net magnetization -→ M of the sample. However, due to thermal movements which prevent a complete alignment, the resulting magnetization vector is very weak. A NMR experiment aims at measuring -→ M. Let us consider the rotating frame of reference (x, y, z) and assume -→

B 0 = B 0 -→ e z .
We can decompose the net magnetization vector as:

-→ M = M x -→ e x + M y -→ e y + M z -→ e z .

The measured signal S is proportional to the transverse component: S ∝ M 2 x + M 2 y . During an experiment in NMR, -→ M is tilted away from its equilibrium position, using short pulses of radio-frequency (RF) radiation -→ B 1 , perpendicular to -→ B 0 , centered at the middle of the NMR spectrum. Applying -→ B 1 on -→ M produces a magnetic moment which causes -→ M to precess about -→ B 1 as well. This additional moment is: -→ m = -→ M × -→ B 1 . This results in a magnetization vector -→ M which deviates from the z axis by a flip angle α. When the pulse ceases, the out-of-equilibrium magnetization vector returns to equilibrium state after a certain delay. This phenomenon is exponential and is characterized by two relaxation constants (cf. figure A.1.1):

• the longitudinal magnetic -or spin-latticerelaxation time T 1 is relaxation constant of the longitudinal component M z of the net magnetization vector returning to its initial state: M z (t) = M 1e -t/T 1 ; (A.1.2)

• the transverse -or spin-spinrelaxation time T 2 is relaxation constant of the transverse component M xy of the net magnetization vector when returning to zero. It is the average time for the spins of precessing nuclei to dephase, due to mutual interactions: M xy (t) = M sin (α)e -t/T 2 . (A.1.3)

The relaxation constants T 1 and T 2 depend the amplitude of the magnetic field -→ B 0 and on the physicochemical environment: the biological composition of the tissue of interest, as well as controlled parameters such as the temperature and the pressure. However in practice, the local field seen by the nuclei is not spatially homogeneous, due to -→ B 0 being non-homogeneous itself, but also to the chemical environment created by the biological tissue. This makes the resonant frequency of the spins vary continuously across space. The phase shift between the spins increasing with time, this results in an different exponential decrease in transverse magnetization, modeled by different relaxation time T ⋆ 2 < T 2 (cf. figure A.1.1). The values of T 1 , T 2 and T ⋆ 2 change signifi-cantly depending on the tissue in which the protons are located, and are the source of the contrast in MRI.

The visibility of protons in MRI depend mainly on their relaxation time. Therefore, proton density reflects water and fat concentration in tissues. For example, a healthy muscle contains about 80% water. By choosing a sequence of RF pulses and adequate parameters and -→ B 0 gradients, one can have the different biological structures -e.g. muscle tissues and fat tissues -to take different intensities, allowing to separate them visually. For instance, a T 2 -weighted sequence provides images in which the intensity is a increasing function of the T 2 value of the tissue in this location. There exist many sequences for acquiring T 1 -weighted, T 2 -weighted or T ⋆ 2 -weighted images, in which the biological tissues are shown with different intensity distributions.

a.2 nmr imaging

In a volume, it is possible to make localized measurements of the amplitude of the components of -→ M, using magnetic field gradients and strategies involving certain Fourier space encoding, resulting in an multidimensional image (2D or 3D).

For the task of segmenting images of muscles, it is convenient to use. a choice of sequence and sequence parameters which reveals as much as possible the boundary of the muscles. The muscles are surrounded by a thin bio-membrane of connective tissue -the fasciae (singular: fascia) -and generally separated by a layer of fat tissues. Fasciae are made of collagen fibers, with the property of having an ultra-short T 2 which, in most typical modalities yields very little signal and show in black in images: if the transverse component of -→ M is fully relaxed, the output signal is zero whatsoever. Muscle and fat tissues have different and well know relaxation constants, and some methods have been developed to give them different intensities in images, for easy visual discrimination. For instance, some T 1 -weighted images can be used for separating muscle and fat tissues, even thought the fasciae are difficult to see (cf. figure 1.1.1).

Some modalities reveal the muscle/fat transition more strongly than others. For example, gradient-echo sequences produce T * 2 -weighted images using the fact that the spin of fat protons (i.e. protons in fat-tissues) precess in a different rate than muscle (water) protons. In regions containing both kinds of protons, this causes a phase difference which results in a decreased amplitude of the transverse component at the acquisition time. In these images, muscle boundaries show a deep contrast (cf. figure 1.1.1). Among T ⋆ 2 -weighted modalities, the Dixon method combines several acquisitions in order to provide separate images for fat-composed tissues and a water-composed tissues (muscle tissues and vessels). This method provides convenient material for segmenting muscles (cf. figure A.2.2):

• a map of muscle tissues only (vessels are negligible volume). We may use this water-map for easy separation of the muscles from fat tissues; • an opposed-phase image, acquired are the instant when fat protons and water protons have opposed-phase spins, which enhances the boundaries. We will use this image in segmentation techniques relying on having contrasted boundaries for the objects of interest.

One should note that healthy muscles tissues cannot be discriminated from their appearance, since all muscles are made of the same biological tissue, which results in identical intensities and textures. The differences of intensities observable in the presented muscles are due to spatial inhomogeneities of the magnetic fields.

∀ x ∈ R |U | , x ⊤ L U x > 0.
Proof. We recall the definition of L U :

L U (i , j ) =          d i if i = j ,
w i , j if e i j ∈ E , 0 otherwise, where d i = ∑ j ∈U ∪ M w i j . Consider the following: Let us consider the index k of any nonzero entry of x: ∀ k ∈ {i ∈ U / x i = 0}. There exists a path P k l ⊂ E of edges with a nonzero weight between node v k and a marked node v l ∈ V M . Then (a)∀ e i j ∈ P k l ∩ ( V U × V U ) equation B.1.6 gives: x i = x j ; (b) ∀ e i j ∈ P k l ∩ ( V U × V M ) equation B.1.7 gives: x i = 0. Both results lead to ∀ k ∈ U , x k = 0, which contradicts the first hypothesis B.1.5.

∀ x ∈ R |U | , x ⊤ L U x = ∑ i ∈U d i x 2 i -2 ∑ i , j ∈U

b.2 random walks and probability distributions

Proposition 2. The Random Walks system 4.1.18 can only yield a solution {x s } s that is a probability distribution: Proof. We use the fact that L U is positive definite and thus invertible, and that 1 = [1, . . . , 1] ⊤ is a solution to Lx = 0 and thus

(L U 1 U + B1 M ) = 0.
Consider: In order to prove the non-negativity of x, consider the RW objective and its gradient at location 0 = [0, . . . , 0] ⊤ :

L U x s U = -Bx s M , (B.
E RW (x s U ) = Since all entries of the gradient of E RW (•) at location 0 are non-positive:

(Bx s M ) [i] = -∑ j∈M w ij x s j ≤ 0, (B.2.9) we know that the direction of the solution is towards the positive numbers:

∀i ∈ U, ∀s, x s i ≥ 0.

b.3 random walks, shape prior and probability distributions Proposition 3. Assume the entries of {x s 0 } s are probabilities. Then, the Random Walks system incorporating a prior knowledge term 4.2.10 can only yield a solution that is a probability: Thus the entries of ∇E RW+prior (0) are non-positive and we know the direction of the solution is towards the positive numbers:

∀i, ∀s, x s i ≥ 0.

Figure 1

 1 Figure 1.1.1 -Various axial cross-sections: (left) healthy subject (Dixon imaging); (middle) patient with centronuclear myopathy (Dixon imaging) -mutation of the MTM1 gene coding for the myotubularin; (right) patient with centronuclear myopathy (T 1 -weighted imaging) -unknown mutation.

Figure 1 .

 1 Figure 1.1.2 -(left) Several axial cross-sections of a 3D MRI volume of the lower limbs (T ⋆ 2

Figure 1

 1 Figure 1.1.4 -(left) Axial cross-section of an MRI volume of a right thigh (Dixon imaging);(center) same cross-section with all muscles manually segmented; (right) 3D segmentation of the same volume.

Figure 1

 1 Figure 1.2.1 -Missing and spurious contours in muscle image. (left) Cross-section of athigh from a 37 year-old female subject; (right) same height cross-section from a 29 year-old male subject. Notice the differences in shape and proportion of the muscles between the two subjects.

Figure 2 .

 2 Figure 2.1.1 -Segmentation with snakes: (left) model evolution: the boundary model converges towards the edge of the pear (Kass et al., 1988); (middle) tracking a person's mouth with a snake in a video (ibid); (right) snake-segmented cell in an EM photo-micrograph (McInerney and Terzopoulos, 1996).

Figure 2

 2 Figure 2.1.2 -(left) Local parametrization of simplex meshes (Delingette, 1999); (right) medial representation for constraining the model shape (Gilles and Magnenat-Thalmann, 2010).

Figure 2

 2 Figure 2.1.3 -Muscle segmentation using deformable models with a shape matching prior (Gilles and Pai, 2008): from the first iteration (left) to the last (right).

Figure 2 .

 2 Figure 2.1.4 -Single reference shape priors. (left) An elastic cube model Terzopoulos et al.(1987). (right) The advantages of generalized gradients (bottom) over standard gradient descent (top) using the same deformable model: a realistic shape is kept all along the deformation process[START_REF] Charpiat | Generalized gradients: Priors on minimization flows[END_REF].

Figure 2 .

 2 Figure 2.1.5 -Deformable templates: (top) a template eye made of a few simple mathematical shapes (Yuille, 1990); (bottom) segmentation of a person's mouth, using deformable mouth template (ibid).

Figure 2

 2 Figure 2.1.6 -Examples of landmarks: (left) 32 manually placed landmarks on the surface of a resistor Cootes et al. (2001); (middle) manually placed landmarks on a person's face (Romdhani et al., 1999); (right) manually and interpolated landmarks on the surface of lungs in an X-ray image (Ginneken et al., 2006).

Figure 2

 2 Figure 2.1.7 -Active Shape Model: (left): scatter plot of the landmarks of the training examples (Cootes et al., 1995); (middle) synthesized models using: x = x + Up for different values of p [0] (ibid); (right) segmentation of a heart ventricle in an echo-cardiogram using ASM (ibid).

Figure 2

 2 Figure 2.1.8 -(left) Invalid Gaussian hypothesis: variations of the first mode of a linear ASM of a hand; the most probable shape (center) is the least realistic (Cremers et al., 2003). (right) A kernel-ASM trained on both right and left hand images allows to segment a right hand in a new image (ibid).

Figure 2 .

 2 Figure 2.1.9 -Level set segmentation: (left) illustration of the level set method (from Wikipedia); (right) topological changed during segmentation of synthetic images (Paragios, 2005).

Figure 2

 2 Figure 2.2.1 -Graph-cut: (left) a directed graph G; (right) a cut on G (Boykov and Kolmogorov, 2004).

Figure 2

 2 Figure 2.2.2 -Automatic image partition using normalized cuts (Shi and Malik, 2000b).

Figure 2

 2 Figure 2.2.3 -Two examples of graph-cut segmentation (red curve) using user-provided seeds (green and blue marks) (Sinop and Grady, 2007).

Figure 3

 3 Figure 3.1.1 -(left) Axial cross-section of a right thigh with manually drawn seeds; (right) segmentation.

Figure 3

 3 Figure 3.3.1 -Binary mask (right) generated from water map image (left). This mask is used for registering the reference image onto the test image.

Figure 3

 3 Figure 3.3.2 -(left) Axial cross-section of a right thigh; (right) manual segmentation. The reference segmentation will serve as a reference atlas when segmenting a new image.

Figure 3

 3 Figure 3.4.1 -Generating the unlabeled seeds: (left) points are sampled regularly over the test image with a spacing of k pixels; (center) edge map; (right) unlabeled seeds, after removal of the points lying outside of the mask of the muscles, and after displacement of k/2 pixels in the opposite direction of the gradient of the edge map.

Figure 3

 3 Figure 3.5.1 -The reference atlas is registered onto the test volume, providing priorknowledge on label assignment likelihood for each node.

Figure 3

 3 Figure 3.5.2 -Visual representation of the geodesic distance between seeds: the color-map goes from blue to red as the geodesic distance increases. Notice that edges crossing a visible boundary have a large geodesic distance whereas edges crossing a homogeneous area have a small geodesic distance.

Figure 3 .

 3 Figure 3.5.3 -(left) the probability density of the orientations between muscles M m and M n is estimated by sampling points within M m and determining the closest point in M n . (right) the cost of an edge between nodes p and q depends on the probability of the orientation of vector -→ pq.

Figure 3

 3 Figure 3.6.1 -Labeling error rate (in percent) as box-plots. "Nearest" is when using only the unary potential (w geod = w orient = 0). The label of the seeds is only determined by the geometrically nearest muscle in the reference atlas. "Autoseed (geo. only)" is the error rate when using only the geodesic distance term (w orient = 0, w geod = 0.1). "Autoseeds" refers to the proposed method. Using both binary potentials significantly reduces the error rates.

Figure 3 .

 3 Figure 3.6.2 -Dice coefficients for all muscles as box-plots "Multi atlas registration" refers to a simple atlas registration-based segmentation method (see next chapter for more details). The expression of the Dice coefficient is: D = 2 |T ∩ R| / (|T| + |R|), where T and R are the pixel sets of a given muscle in the inferred segmentation and the reference segmentation respectively.

Figure 3 .

 3 Figure 3.6.3 -Detailed Dice results for all muscles separately.

Figure 3 .

 3 Figure 3.6.4 -Examples of cross-sections of 3d volumes of right thighs automatically segmented with the proposed method.

Figure 3

 3 Figure 3.7.1 -Incorrectly labeled seeds (diamond markers) prevent the RW method to segment the corresponding regions correctly (white regions)

Figure 4

 4 Figure 4.1.1 -(left) Random walker on a discretized ground with target seeds (flags). The probability of the walker of stepping onto an adjacent square/pixel can be different for each pixel, as figured by the red arrows of various sizes. (Right) Segmented ground/image: the label assigned to one pixel is the label of the most probably first reached seed during a random walk starting from that pixel.

  .1.9) 2 we also implemented the following alternate formulation w ij = 1 β|I i -I j |+ε .

Figure 4

 4 Figure 4.1.2 -Image contrast is chosen as the source for the transition probabilities (figured by red arrows). The green curve, following the sequence of lowest transition probabilities, is more likely to be selected as a boundary during the segmentation.

  .1.11) and L, the un-normalized combinatorial Laplacian matrix L = D -A.(4.1.12)

Figure 4 .

 4 Figure 4.2.1 -Appearance model for RW in heart segmentation: (left) original image; (center) seeds as black lines; (right) segmented image. The appearance model reduces the amount of required user interaction[START_REF] Grady | Multilabel Random Walker Image Segmentation Using Prior Models[END_REF]

Figure 4

 4 Figure 4.2.2 -Prior information is introduced to the Random Walks framework by adding a second layer of nodes (the prior nodes); each of these prior nodes has a known assignment probability to each label (here, each color on the second layer represents the most probable prior label). The transition probabilities (figured as vertical arrows) from the first layer to the second layer can vary along the pixels.

Figure 4

 4 Figure 4.3.1 -Label assignment probability of prior model x s 0 , for two different labels s : (left) background; (right) vastus lateralis.

Figure 4 .

 4 Figure 4.3.2 -Weighting schemes, displaying the diagonal of Ω as an image: (left) entropy weighting; (right) confidence map.

Figure 4

 4 Figure 4.3.3 -Examples of segmentations obtained with the RW algorithm with shape priors (constant weighting scheme).

Figure 4

 4 Figure 4.3.4 -Two 3D views of a segmentation (right thigh) obtained with our method.

Figure 4

 4 Figure 4.3.5 -Comparing performances for various weighting schemes; the boxplots are generated from Dice coefficients for all labels (4 muscle database). "Cmap" designates the confidence map weighting scheme.

Figure 4 .

 4 Figure 4.3.6 -Comparison of segmentations using the constant weighting scheme (left)and the confidence map weighting scheme (right).

Figure 4 .

 4 Figure 4.3.7 -Local segmentation differences due to adding an appearance prior term

Figure 4 .

 4 Figure 4.3.8 -Comparing performances with and without the appearance prior; the boxplots are generated from Dice coefficients for all labels (13 muscle database).

Figure 4 .

 4 Figure 4.3.9 -Comparison of our method (bottom box) with the "registration" method (top box) and the method from (Gilles and Pai, 2008) (middle box). Dataset: subset of the 4-muscle dataset with only 15 test volumes.

Figure 4

 4 Figure 4.3.10 -Cross-sections of segmentations obtained with the method from[START_REF] Gilles | Fast musculoskeletal registration based on shape matching[END_REF]. This method always produces realistic segmentations, but sometimes fails to detect the contours.

Figure 4 .

 4 Figure 4.3.11 -Comparison of our method (bottom box) with the "registration" method (top box), the method from (Andrews et al., 2011a) (second box) and our automatic seeding method (third box). Dataset: 13-muscle dataset with 13 samples.

Figure 4 .

 4 Figure 4.3.12 -Per-slice plot of the Dice coefficients for the segmentation method with constant weighting. The darker bars represent the average of the lower quartile, the intermediary bars the median and the lighter bars the average of the upper quartile.

Figure 4

 4 Figure 4.4.1 -Comparison of our RW with shape space method (bottom box) with the simpler RW with prior method with constant weighting (middle box, see section 4.3) and the registration-based baseline method (top box, see section 4.3.4).

Figure 4

 4 Figure 4.4.2 -Segmentation examples. Cross-sections of segmented volumes with: (left column) the registration method x = x 0 ; (middle column) the shape-space component of the inferred segmentation x = x 0 + Uγ; (right column) the inferred segmentation x = x 0 + Uγ + dx.

Figure 5

 5 Figure 5.4.1 -Probability map of one muscle using the distance transform to "soften" the probabilistic segmentation spatially.

Figure 5

 5 Figure 5.4.2 -Test loss for three different methods. Each figure corresponds to a differentλ -that is, the SVM regularization parameter. Curves are plotted with respect of λ ′ -that is, the proximal regularization parameter, which penalizes deviations from the hand-tuned parameters w 0 . "Baseline" corresponds to the standard struct-SVM method using distance-transformed pseudo-soft segmentations zk as annotation vectors.

Figure 5

 5 Figure 5.4.3 -(left column) segmentation using w 0 (hand-tuned parameters); (right column) segmentation using learning w using latent structured SVM. Segmentations with learned parameters are more accurate.

Figure A. 1

 1 Figure A.1.1 -Longitudinal (left) and transverse (right) relaxations. The transverse component decreases with a time constant T 2 in a homogeneous field -→ B 0 , and T ⋆ 2 in an inhomogeneous field.

Figure A. 2

 2 Figure A.2.1 -Axial-cross sections of thigh: (left) T 1 -weighted; (right)T ⋆2 -weighted with fat-suppression (a special method was applied during the acquisition of this volume to ignore signal from fat tissues).

Figure

  Figure A.2.2 -Dixon method, axial-cross sections of thigh: (left) water-map; (rig) opposedphase image.

BM

  AT H E M AT I C A L A P P E N D I C E S b.1 positive definiteness of L U Proposition 1. Assume there exists at least one path composed of nonzero-weight edges between each unknown node and one seed. Then L U is positive definite:

Using

  Reductio ad absurdum, suppose:∃ x ∈ R |U | / { 0 } s.t. x ⊤ L U x = 0. (B.1.5)For this hypothesis to stand, we must have:∀i , j ∈ U , w i j x i -

  We use the fact that (L + ∆) is positive definite -and thus the solution to (L + ∆) x = b is unique -and that 1 = [1, . . . , 1] T is a solution to Lx = 0. Consider:(L + ∆) x s = ∆x s 0Considering the RW objective and its gradient at location 0 = [0, . . . , 0] T :E RW+prior (x s ) = x sT Lx s + x s -

Wikipedia, on 01.01.2013

1.2 problem definition

Box-plot presentation: the boxes contain the middle 50% of the data and the median value, and the extremities of the lines indicate the min and max values, excluding the outliers (for more details, see the documentation of Matplotlib).

Thanks to Danny Goodman for his help in this section.

Thanks to Danny Goodman for his help in this section.

Box-plot presentation: the boxes contain the middle

50% of the data and the median value, and the extremities of the lines indicate the min and max values, excluding the outliers (for more details, see the documentation of Matplotlib).

I would like to thank B. Gilles for performing these experiments

I would like to thank S. Andrews for performing these experiments

unlike what we published in(Baudin et al., 2012a). This is due to using a slightly different method here, where the labels are not weighted differenly, for remaining true to the probabilistic framework presented in section 4.2.2.

This part was developed with the help of Danny Goodman and Pawan Kumar.

This section was developed with the help of Puneet Kumar and Pawan Kumar

This idea was suggested by DannyGoodman.

parameter estimation using latent svm

Upper Bound on the Risk

Given a dataset D = {(x k , z k ) , k = 1, . . . , N}, which consists of inputs x k and their hard segmentation z k , we would like to estimate parameters w such that the resulting inferred segmentations are accurate. Here, the accuracy is measured using the loss function ∆ (•, •). Formally, let ỹk (w) denote the soft segmentation obtained by minimizing the energy functional E (•, x k ; w) for the k-th training sample, that is, ỹk (w) = arg min y w ⊤ ψ (x k , y) .

(5.3.1)

We would like to learn the parameters w such that the empirical risk is minimized over all samples in the dataset. In other words, we would like to estimate the parameters

The above objective function is highly non-convex in w, which makes it prone to bad local minimum solutions. To alleviate this deficiency, the latent SVM formulation upper bounds the risk for a sample (x, z) as follows:

The first inequality follows from the fact that the prediction ỹk (w) has the minimum possible energy (see equation 5.3.1). Thus, its energy has to be less than or equal to the energy of any compatible segmentation ŷ with ∆ (z k , ŷ) = 0. The second inequality is true since it replaces the loss augmented energy of the prediction ỹk (w) with the minimum loss augmented energy.

relation to the crf case In some cases -for instance, when learning conditional random fields [START_REF] Szummer | Learning CRFs using graph cuts[END_REF][START_REF] Komodakis | Efficient training for pairwise or higher order crfs via dual decomposition[END_REF] -the loss function is concave, so that the minimum loss is only reached when the prediction ŷ is equal to the annotation z k : min

As a result, the upper bound on the empirical risk becomes

This upper bound is convex in w, since the minimum of a collection of affine functions is concave (which becomes convex with thesign). As a result, one obtains a learning algorithm by minimizing the upper bound for each sample in the training set:

where the slack variable ξ k represents the upper bound of the risk for the k-th training sample. The regularizing term ||w|| 2 , weighted by hyper-parameter λ, ensures that we do not over-fit to the training samples. The previous algorithm is a convex quadratic program and corresponds to the formulation of Structural SVM with rescaled margin [START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF], for which may efficient algorithms have been proposed (Joachims et al., 2009;[START_REF] Shalev-Shwartz | Pegasos: Primal estimated sub-gradient solver for svm[END_REF][START_REF] Taskar | Max-margin markov networks[END_REF][START_REF] Tsochantaridis | Support vector machine learning for interdependent and structured output spaces[END_REF]. However, as mentioned earlier, the loss function ∆ (•, •) is not concave in w, which in turns leads to a non-convex upper bound on the empirical risk. Equation 5.3.6 is not verified, and solving the baseline Struct-SVM problem 5.3.8 is bound to fail in our case.

Formulation with Latent SVM

Let us recall the upper bound of the empirical risk derived above:

While this upper bound is not convex, it is a difference of two convex functions in (Yu and Joachims, 2009). This observation allows us to obtain a local minimum or saddle point solution using the CCCP algorithm (Yu and Joachims, 2009;[START_REF] Yuille | The concave-convex procedure (cccp)[END_REF], which iteratively improves the parameters starting with an initial estimate w 0 (the detailed optimization algorithm will be described shortly).

Annotation Consistent Inference 2

Given an input x and its hard segmentation z, ACI requires us to find the soft segmentation y ⋆ with the minimum energy, under the constraint that it should be compatible with z (see step 3 of Algorithm 5.1). We denote the ground truth label of a voxel i by s i , that is, s i = arg max s z(i, s), and the set of all voxel indices by V. Using our notation, ACI can be formally specified as

y ⊤ L (x; w) y + E prior (y, x; w) .

(5.3.11)

Here, C(V ) is the set of all compatible probabilistic segmentations, that is,

y (i, s) ≥ 0, ∀s ∈ S, (5.3.12) ∑ s∈S y (i, s) = 1, ∀i ∈ V , (5.3.13) y (i, s i ) ≥ y (i, s) , ∀i ∈ V , ∀s ∈ S.

(5.3.14)

Constraints 5.3.12 and 5.3.13 ensure that y ⋆ is a valid probabilistic segmentation. The last set of constraints 5.3.14 ensure that y ⋆ is compatible with the hard ground truth z. Note that in the absence of constraint 5.3.14, the other constraints are not necessary (see section 4.2.2) and the above problem can be solved efficiently using the RW algorithm. However, since the ACI problem requires the additional set of compatibility constraints, we need to develop a novel efficient algorithm to solve the above convex optimization problem. To this end, we exploit the powerful dual decomposition framework [START_REF] Bertsekas | Nonlinear programming[END_REF][START_REF] Komodakis | MRF Optimization via Dual Decomposition: Message-Passing Revisited[END_REF]. (5.3.15) that is, any problem which can be exactly decomposed into a sum of (preferably simpler) sub-problems with functionals g m (y). (5.3.16)