Keywords: Segmentation, Level-Sets, Finite Element Method, Local Statistics, Anatomical Head Model Segmentation, Bayesian formulation, Local statistics, Level-set, Mumford-Shah, GPU

Titre : Segmentation d'Image par M éthode des Ensembles de Niveaux et Application à la Cr éation de Mod èles Anatomiques de la T ête

R ésum é : L'apparition des techniques d'imagerie par résonance magnétique (IRM) à la fin du XXe siècle a révolutionné le monde de la médecine moderne, en permettant de visualiser avec précision l'intérieur de structures anatomiques de manière non invasive. Cette technique d'imagerie a fortement contribué à l'étude du cerveau humain, en permettant de discerner avec précision les différentes structures anatomique de la tête, notament le cortex cérébral.

Le discernement de structures anatomiques de la tête porte le nom segmentation, et consiste à "extraire" des régions dans les IRMs. Plusieurs méthodes de segmentation existent, et cette thèse portent sur les méthodes à base d'évolution d'hypersurfaces: une hyper-surface (surface en 3D) est progressivement déformée pour finalement épouser les frontières de la region à segmenter.

Un modèle de tête correspond au partitionnement de la tête en plusieurs structures anatomiques préalablement segmentées. Un modèle de tête classique comprends en général 5 structures anatomiques (peau, crane , liquide céphalo rachidien, matière grise, matière blanche), imbriquées les unes dans les autres à la manière de "poupée russes". Néanmoins de par la complexité de leurs formes, segmenter ces structures manuellement s'avère pénible et extrêmement difficile.

Cette thèse se consacre à la mise en place de nouveaux modèles de segmentation robustes aux altérations d'IRMs, et à l'application de ces modèles pour la création automatique de modèles anatomiques de la tête. Apres avoir brièvement survolé un état de l'art des differentes méthodes de segmentation d'image, deux contributions à la segmentation par évolution d'hypersurface sont proposées. La première constitue une nouvelle représentation et un nouveau schémas numérique pour la méthode des ensembles de niveaux, en utilisant des éléments finis quadrilateraux. Cette representation vise à améliorer la qualité et la robustesse du modèle. La seconde contribution constitue un nouveau modèle de segmentation basé sur des statistiques locales, robuste aux altérations présentes dans les IRMs. Ce nouveau modèle vise à unifier plusieurs modèles "état de l'art" en segmentation d'image. Enfin, un cadre pour la création automatique de modèle de tête est proposé, utilisant princpalement le precédent modèle de segmentation par statistiques locales.

Mots cl és : segmentation, ensembles de niveaux, m éthode des él éments finis, statistiques locales, mod èle anatomique de la t ête Title: Image Segmentation with Level Set Method and Application to Anatomical Head Model Creation

Abstract: Magnetic Resonance Images (MRI) have been introduced at the end of the XXth century and have revolutionized the world of modern medicine, allowing to view with precision the inside of anatomical structures in a non-invasive way. This imaging technique has greatly contributed to the study and comprehension of the human brain, allowing to discern with precision the different anatomical structures composing the head, especially the cerebral cortex.

Discernment between these anatomical structures is called segmentation, and consist in "extracting" structures of interest from MRIs. Several models exists to perform image segmentation, and this thesis focus on those based on hypersurface evolutions: an hypersurface (surface in 3D) is incrementally adjusted to finally fit the border of the region of interest.

A head model corresponds to the partitioning of the head into several segmented anatomical structures. A classic head model generally includes 5 anatomical structures (skin, skull, cerebrospinal fluid, grey matter, white matter), nested inside each other in the manner of "Russian nested dolls". Nevertheless because of the complexity of their shapes, manual segmentation of these structures is tedious and extremely difficult. This thesis is dedicated to the creation of new segmentation models robust to MRI alterations, and to the application of these models in the purpose of automatic creation of anatomical head models. After briefly reviewing some classical models in image segmentation, two contributions to segmentation based on hypersurface evolution are proposed. The first one corresponds to a new representation and a new numerical scheme for the level-sets method, based on quadrilateral finite elements. This representation aims at improving the accuracy and robustness of the model. The second contribution corresponds to a new segmentation model based on local statistics, and robust to standard MRI alterations. This model aims at unifying several 'state-of-the-art' models in image segmentation. Finally, a framework for automatic creation of anatomical head models is proposed, mainly using the previous local-statistic based segmentation model. CHAPTER

INTRODUCTION (VERSION FRANC ¸AISE) CONTEXTE

Les neurosciences sont un domaine fascinant. Elles ont gagné en popularité au cours des dernières décennies, principalement gr âce aux récentes avancées dans le domaine de l'imagerie médicale et fonctionnelle. Les premiers travaux de recherche consistaient à analyser les formes du cerveau pour en déceler certaines caractéristiques, mais n'étaient d'aucune aide pour comprendre l'activité neuronale. Avec l'apparition d'équipements permettant des études fonctionnelles du cerveau, plusieurs méthodes ont émergées pour étudier en temps réel l'activité neuronale (TEP-scan, IRMf, Optical Imaging, EEG, MEG . . . ). Chacune de ces méthodes est basé sur une hypothèse sous-jacente de cette activité:

• la Tomographie pas Émission de Positon (TEP-scan), et l'Imagerie par Résonance Magnetique fonctionnelle (IRMf) se basent sur le fait que les changements de circulation sanguine et d'oxygénation du sang dans le cerveau (connu sous le nom d'hémodynamique) sont étroitement liés à l'activité neuronale. Le TEP-scan mesure le flux sanguin à l'intérieur du cerveau, par l'intermédiaire de traceurs radioactifs préalablement injecté dans le sang. L'IRMf se base sur le même principe que les IRMs classiques, mais se différencie en étant plus sensible à l'oxygène. Cette méthode est donc capable de détecter les grandes quantités d'oxygène libérées par le sang au niveau des neurones actifs.

• l'ElectroEncephaloGraphie (EEG) et la MagnetoEncephaloGraphie (MEG) se basent sur le fait que l'activité neuronale produit un courant électrique, en raison de la nature chimique des connexions synaptiques. L'EEG mesure les potentiels électrique au niveau du scalp, tandis que la MEG mesure le champ magnétique induit par cette activité électrique [START_REF] Wolters | Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: A simulation and visualization study using high-resolution finite element modeling[END_REF].

Ces méthodes ont chacune leurs propres avantages et leurs propres inconvénients. Elles sont caractérisées par leurs résolutions spatiales et résolutions temporelles, et visent à minimiser ces 2 grandeurs. Les méthodes fournissant la meilleure résolution temporelle sont celles basées sur la mesure de l'activité électrique neuronale (MEG/EEG).

BESOIN DE MODELES ANATOMIQUES

Cependant, la résolution spatiale des méthodes d'EEG/MEG dépend de la précision de reconstruction du problème inverse, qui consiste à estimer l'emplacement des sources électriques à l'intérieur du cerveau à partir de la série de mesures sur le scalp. Cette reconstruction est réalisée gr âce aux équations de Maxwell, qui décrivent des propriétés reliant les champs électriques et magnétiques à leurs sources électriques respectives (on cherche donc à résoudre ces équations " à l'envers"). Toutefois, il est dit dans ces équations que les champs électriques et magnétiques sont très sensibles à la nature des zones qu'ils traversent (surtout pour l'EEG).

La tête humaine est très complexe et comprend divers types de structures anatomiques conduisant différemment l'électricité. Par conséquent, afin de réaliser correctement la reconstruction inverse, l'emplacement des différentes structures anatomiques de la tête doit être connu à l'avance. Estimer l'emplacement de ces structures consiste à créer un modèle de tête, qui partitionne la tête suivant la conductivité des tissus. En général, ces modèles comprennent 5 structures anatomiques, qui sont imbriquées les unes et dans les autres à la manière de "poupées russes" (Fig. 1.1). De l'extérieur vers l'intérieur, ces modèles sont constitués de : la peau, le cr âne, le liquide céphalo-rachidien (LCR), la matière grise (MG), et la matière blanche (MB) [START_REF] Blinkov | The Human Brain In Figures and Tables[END_REF].

La structure la plus importante dans ces modèles est sans nulle doute la matière grise, qui constitue le cortex cerebral. Anatomiquement, le cortex est une "couche" de quelque millimètres d'épaisseur, fortement convoluée et de topologie sphérique. C'est dans cette couche que reposent tous les neurones du cerveau, la matiere blanche ("intérieur" du cerveau) ne contenant que les fibres nerveuses reliant les neuronnes entre elles. La detection de l'activité neuronale se fait a l'intérieur du cortex, ce qui implique que la qualité de son modèle doit etre primordiale.

Plusieurs méthodes existent pour effectuer ces types de modèles, et appartiennent au domaine de la segmentation d'image.

SEGMENTATION D'IMAGE

La segmentation d'image consiste à partitionner une image en "zones significative". Ces zones correspondent en général à des objets à extraire et peuvent être considérées comme des régions séparées par des frontières distinctes ou des régions caractérisées par des propriétés globales. Dans notre cas, ces zones correspondent aux différentes structures anatomiques présentes dans la tête, différenciées entre elles par leur intensité dans les IRMs anatomiques. Toutefois, il est fréquent que les IRMs anatomiques souffrent d'altérations conduisant a des problèmes de bruits et d'inhomogénéités, rendant la t âche de segmentation non triviale.

Plusieurs méthodes existent pour effectuer des segmentations de régions, et nous nous sommes principalement concentrés sur celles dérivées d'évolutions d'hypersurfaces: une hypersurface (courbe en 2D ou surface en 3D) est initialisée dans une image, et est progressivement déformée pour finalement épouser les frontières de la region à segmenter. Cette méthode permet d'intégrer des connaissances a-priori sur la forme de la structure à segmenter, en particulier des contraintes topologiques requises pour la segmentation du cortex.

ORGANISATION ET CONTRIBUTION DE CETTE TH ÈSE

Ce manuscrit est organisé en 4 chapitres, qui reflètent les différentes contributions de cette thèse, et 2 annexes présentant des détails supplémentaires:

Le chapitre 3 couvre un état de l'art des méthodes de segmentation d'images. Après un bref aperc ¸u de certains modèles discrets, ce chapitre se concentre sur les méthodes à base de propagation de front pour la segmentation, notamment ceux utilisant la méthode des ensembles de niveaux.

La méthode des ensembles de niveaux consiste à représenter le front comme étant le niveau zéro d'une fonction de dimension supérieure, en général choisie comme fonction distance au front. Cette méthode permet de gérer facilement les changements de formes du front, notamment ceux entraînant un changement de topologie. L'évolution du front est ensuite exprimée à travers une équation aux dérivées par-tielles, découlant de la minimization d'une fonctionelle définie:

• sur l'hypersurface elle-même dans le cas des modèles basés contours. Ces modèles consistent à faire évoluer librement le front dans les régions homogènes, et à stopper son évolution au niveau des zones de contours dans l'image, correspondant en général à des zones de fort gradient.

• sur les régions délimitées par l'hypersurface dans le cas des modèles basés régions. Ces modèles consistent à faire évoluer le front de fac ¸on à séparer significativement les 2 régions qu'il délimite. Dans le cas d'une image scalaire, cette évolution revient à séparer les régions claires des regions sombres.

Le chapitre 4 présente un nouveau schéma numérique pour la méthode des ensembles de niveaux à l'aide d'éléments finis quadrilatéraux. Similairement à l'approche de "la bande étroite", le calcul de la fonction distance est restreint aux éléments qui sont proches du niveau zéro, et peut donc être représentée par un ensemble connexe d'éléments finis. Pour garder une trace du niveau zéro au cours de son évolution, certains éléments sont insérés dans la bande, et d'autres retirés.

Cette implémentation a pour but d'améliorer la robustesse du modèle, en résolvant simultanément l'équation d'évolution et l'équation Eikonale gr âce à un schéma en temps implicite. Toutefois, cette implémentation alourdit la charge de calcul, et augmente de fac ¸on non négligeable la complexité du modèle. Quelques détails concernant cette implémentation peuvent être trouvés en Annexe A.

Le chapitre 5 présente un modèle de segmentation rapide basé sur des informations locales. Ce modèle combine des statistiques locales à l'intérieur et à l'extérieur du front pour le faire évoluer. Comme il est basé sur des décisions locales, ce modèle est plus robuste aux variations des régions d'intérêt (en intensité, contraste, bruit, . . . ). Ces statistiques locales sont exprimées à travers des convolutions Gaussiennes, ce qui permet une implémentation efficace gr âce à des filtres récursifs. Ce modèle généralise plusieurs modèles état-de-l'art, gr âce à un nouveau paramètre de voisinage dont le rôle est étudié. Ce modèle introduit toutefois de nouvelles difficultés qui sont inhérentes au fait de baser une propriété globale (la segmentation), sur de pure décisions locales. La seconde partie de ce chapitre explore certaines de ces difficultés et propose quelques corrections. Enfin, une implémentation efficace est proposé, en utilisant la puissance de calcul de dispositifs modernes hautement parallèles tels que les GPU (Graphic Processing Units).

Le chapitre 6 porte sur la création de modèles anatomiques de la tête à partir d'IRMs bi-modales T1 et T2. Après avoir brièvement présenté les différentes structures anatomiques de la tête, un cadre est proposé pour créer automatiquement un modèle de tête en segmentant successivement les différentes structures anatomiques:

• La peau est segmentée en 2 phases: une segmentation grossière avec contrainte de regularité importante, puis plus fine une avec moins de régularité pour rentrer dans des détails d'importance.

• Le liquide c éphalo rachidien (LCR) est ensuite segmenté, et constitue une enveloppe englobante du cerveau ("skull-stripping").

• La mati ère blanche est ensuite segmentée, en partant de la segmention du LCR et en faisant rétracter le front sur lui même.

• Puis la mati ère grise est segmentée, en partant de la segmentation de la matiere blanche et en faisant "gonfler" le front. Une contrainte topologique est proposée pour maintenir une topologie sphérique.

• Finalement le crane est finalement segmenté en partant d'une couche créée à partir de la segmentation du LCR, puis en faisant évoluer le front dans la zone comprise entre la peau et le LCR.

L'Annexe A présente quelques détails sur l'implémentation de la méthode des ensembles de niveaux gr âce aux éléments finis, proposée dans le chapitre 4.

L'Annexe B présente la dérivée de la fonctionnelle introduite dans le chapitre 5. Cette dérivée est calculée gr âce au gradient de forme [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients[END_REF], et inclue quelques termes supplémentaires qui augmentent la robustesse du modèle. CHAPTER

INTRODUCTION CONTEXT

Neuroscience has always been a challenging and fascinating domain. It has gained in popularity in the last decades mainly due to recent advances in medical imaging and functional imaging. Early researches were based on anatomical and morphological studies of the brain, and were of no help for understanding neural activity. With the apparition of functional imaging devices, several methods have been proposed to study in real-time neural activity (PET-scan, fMRI, optical imaging, EEG, MEG, . . . ). Each of these methods is based on an underlying assumption of this activity:

• Positron Emission Tomography (PET-scan) and functional Magnetic Resonance Imaging (fMRI) are based on the fact that changes in blood flow and blood oxygenation in the brain (collectively known as hemodynamics) are closely linked to neural activity. PET-scan measures the blood flow inside the brain, by injecting radioactive tracers into the blood stream and detecting them using nuclear imaging. fMRI follows the same principle as MR images, but is more sensitive to oxygen and is thus able to detect large amounts of oxygen released by blood inside active neurons.

• ElectroEncephaloGraphy (EEG) MagnetoEncephaloGraphy (MEG) are based on the fact that active neurons produce some electrical activity, due to the chemical nature of synaptic transmissions. The EEG method measures the electric potential created by the induced electrical field, while the MEG measures the magnetic field created by this electrical activity [START_REF] Wolters | Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: A simulation and visualization study using high-resolution finite element modeling[END_REF].

Each methods has its own advantages and drawbacks, and are characterized by their spatial resolutions and temporal resolutions. The methods providing the best temporal resolution are without a doubt those based on electrical activity measurement (MEG/EEG).

NEED FOR ANATOMICAL MODEL OF THE HEAD

However, spatial resolution in the MEG/EEG methods strongly depends on the accuracy of the inverse problem reconstruction, which consists in computing locations of electrical sources inside the brain from the set of measurements on the scalp. This inverse problem is performed using the Maxwell's equations, that describe the properties of the electric and magnetic fields and relate them to their electrical sources. However, it is stated in these equations that both fields are highly sensitive to the conductivity of areas they cross (especially for EEG).

The human head is very complex and includes various types of anatomical structure that conduct differently electricity. Therefore, in order to accurately perform the inverse problem, one has to know precisely the different anatomical structures present in the head. Differentiating between these structures consist in creating a model of the head, which partition the head according to the local conductivity of these structures. In general, these models include 5 structures, that are nested inside each others in the manner of "Russian nested dolls" (Fig. 2.1). They include, from outer to inner structures : the skin, the skull, the CerebroSpinal Fluid (CSF), the Grey Matter (GM), and the White Matter (WM) [START_REF] Blinkov | The Human Brain In Figures and Tables[END_REF].

The most important structure in these models is without a doubt the grey matter, that constitutes the cerebral cortex. Anatomically, the cortex is a highly convoluted layer of grey matter, with spherical topology. It includes all the neurons of the brain, the white matter ("inside" of the brain) including only the nervous fiber linking these neurons to each others. Neural activity is detected inside the cortex, and accurary of its model is therefore crucial.

Several methods exist to perform these kind of models, that have been widely studied in the image segmentation domain.

IMAGE SEGMENTATION

Image segmentation consists in partitioning an image into "meaningful areas". These areas correspond in general to objects to be extracted from images, and can be seen as regions separated by distinct boundaries or regions characterized by some global properties. In our case, these areas correspond to different anatomical structures present in the head, differentiated according to their intensities. However, several alterations coming from the imaging device can perturb the MRI, leading in most of the case to noise and intensity drifts. These alterations make the segmentation non trivial, and models robust to these alterations are required Several models exist to perform image segmentation, but we focused on those based on hypersurface evolutions: an initial curve or surface is initialized in the image, and is incrementally adjusted to finally fit the border of the region of interest. This method allows to incorporate some knowledge about the shape of the anatomical structure, such as topological constraints in order to segment the cortex.

ORGANIZATION AND CONTRIBUTIONS OF THIS THESIS

This manuscript is organized in four chapters, which reflect the different contributions of this thesis, and two appendices presenting some additional details: Chapter 3 covers state-of-the-art methods for image segmentation. After briefly reviewing some discrete models such as Markov random field and graph cuts, this chapter focuses on front propagation methods for image segmentation, especially on those using the level-set method. These models are then classified in two separate categories : edge-based segmentation models, and region-based segmentation models.

Chapter 4 presents a novel implementation of the level-set method using the finite element method with quadrilateral elements. Similarly to the narrow band approach, the computation of the distance function is restricted to the elements that are close to the zero level-set, and thus can be represented as a sparse set of connected elements. To keep track of the zero level-set during its evolution, some elements are inserted to the set, and others removed from the set depending on the motion of the curve.

Chapter 5 presents a fast and robust segmentation model based on local informations. It combines local statistics of the regions along the contour (inside and outside) to drive the segmentation procedure. Since it is based on local decisions, this model is more robust to local variations of the regions of interest (contrast, noise, blur, . . . ). In order to avoid the computational burden of a direct estimation, these local statistics are expressed as the result of Gaussian convolutions. This makes an efficient implementation via recursive filters possible, and gives a complexity of the same order as methods based on global statistics. This model generalizes several state-of-the-art segmentation models, using a neighborhood parameter whose role is investigated. This model also introduces some new difficulties which are inherent to the fact of basing a global property (the segmentation) on pure local decisions. The second part of this chapter explores some of those difficulties and propose some possible corrections. Finally, an efficient implementation is proposed, using the computational power of highly parallel devices such as modern GPUs (Graphic Processing Units).

Chapter 6 presents a method for the creation of anatomical head models from bimodal T1-and T2-weighted MRIs. After briefly reviewing the anatomical structures present in the head, a framework is proposed to segment successively the different anatomical structures included in the model by using the local-statistics model introduced in the previous chapter.

• The scalp is first segmented by starting from a bounding box of the MRI.

• The CerebroSpinal Fluid (CSF) is then segmented by initializing the front as a sphere inside the brain.

• The White Matter (WM) is then extracted by starting from the previously segmented CSF.

• The Grey Matter (GM) is then extracted by starting from the previously segmented WM.

• The skull is finally segmented by initializing the front as a layer between the CSF and the scalp.

Appendix A presents some details about the finite-element implementation of the level-set method proposed in chapter 4.

Appendix B presents the derivative of the local-statistics functional introduced in chapter 5. This derivative is performed using the shape gradient method [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients[END_REF], and include some additional terms that increase the robustness of the model.

OVERVIEW

This chapter covers state-of-the-art methods for image segmentation. After explaining the purpose of segmentation, we briefly review some famous discrete models for image segmentation such as Markov random field and graph cuts. We then focus on front propagation methods for image segmentation, and particularly on those using the level-set representation. Finally, we give a detailed insight on segmentation models using the level-set method, by classifying them in two separate categories : edge-based segmentation models, and region-based segmentation models.

Keywords: Graph-cuts, Snakes, Level-set, Mumford-Shah, Chan-Vese

INTRODUCTION

Image segmentation is a domain that has been widely studied in the last twenty years, and still remains a challenging task [START_REF] Morel | Segmentation of images by variational methods: A constructive approach[END_REF][START_REF] Morel | Variational Methods in Image Segmentation[END_REF]. It consists in partitioning an image into "meaningful areas". These areas correspond in general to objects to be extracted from images, and can be seen as regions separated by distinct boundaries or regions characterized by some global properties.

Early approaches in image segmentation were based on simple heuristics, and consisted in several different steps [START_REF] Perkins | Area segmentation of images using edge points[END_REF]. Recent approaches are now mainly based on optimization methods, and consist in minimizing a given energy to obtain the desired partition. These energies differ depending on the type of segmentation desired, and are in general formulated as a weighted sum of different terms. Most of the times, these energy include two terms : a data attachment term, standing for how the image will be partitioned depending on its properties (intensity, color, texture ...), and a partition term, defining some prior knowledge about the shape of the desired partitions (smoothness, topology, global form...). These optimization-based approaches can be classified in two distinct categories : spatially discrete and spatially continuous representations.

Spatially discrete optimization methods rely usually on expressing images as graphs, where the pixels are seen as nodes and adjacent pixels are linked by weighted edges (depending on a given connectivity). These weights denote in general the likelihood of the two linked pixels to belong to the same region. The segmentation is then performed by minimizing a discrete energy defined on the edges of the graph, by using some well-known algorithms from the graph theory domain. For example, one of the most famous discrete optimization method is the graph cuts method, first proposed by Greig et al. [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF] and recently popularized by Boykov et al. [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]: the segmentation is seen as a cut in the graph which has minimal cost.

Spatially continuous optimization methods rely on evolving closed curves or closed surfaces in images, to finally fit the border of the regions to extract.

Curve evolution for image segmentation was first proposed in the seminal work of Kaas et al. [START_REF] Kass | Snakes: Active contour models[END_REF], often referred as Snakes : the evolution is performed by minimizing a continuous energy defined along the contour. This technique was a breakthrough in early development of segmentation models, and have been rapidly extended to overcome some of its important limitations. These evolving contours were initially represented as explicit models, which means they were parameterized by a set of points along the contour, and therefore seen as polygons. Although this representation is simple and efficient, it also brings some important drawbacks, especially concerning the re-parameterization process.

An important step has been reached by using implicit models to represent these contours. This representation, also called level-set representation, was first introduced by Dervieux and Thomasset [START_REF] Dervieux | A finite element method for the simulation of Rayleigh-Taylor instability[END_REF] and later re-discovered and largely popularized by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF]. It consists in expressing the contour as the zero level of a higher dimensional function, most of the time chosen as the distance function to the contour. While being computationally more expensive, this representation overcomes some important limitations of explicit models, mostly due to the absence parameterization of the contour.

Several energies have been proposed to make these contours evolve, and can be globally classified into two categories : edge-based models, and region-based models.

Edge-based models consist in evolving the contour in homogeneous areas, and locally stopping it when it reaches high image gradients [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF][START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF]. Among these models, the geodesic active contour model [START_REF] Caselles | Geodesic active contours[END_REF] is certainly the most elegant one: the segmentation problem is reformulated as finding cyclic geodesics (paths of minimal length) in a Riemanian space whose metric is inversely proportional to the norm of the image gradient.

Region-based models consist in evolving the contour in order to separate significantly the statistics of the regions it delimits [START_REF] Chan | Active contours without edges[END_REF][START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF][START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF][START_REF] Paragios | Geodesic active regions: a new paradigm to deal with frame partition problems in computer vision[END_REF][START_REF] Cremers | A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape[END_REF]. For example, separating bright areas from dark ones, or extracting areas matching a certain histogram [START_REF] Kim | Nonparametric methods for image segmentation using information theory and curve evolution[END_REF]. These models take into consideration the whole region and are thus more robust to noise and to initialization than edge-based models.

In the following, we give a detailed insight of all the preceding notions. After briefly reviewing discrete approaches in image segmentation, we will present in detail techniques based on front evolution for image segmentation. We will particularly focus on the level-set method, and how it can be efficiently implemented. Finally, we will review several geometric flows using the level-set method,categorized into two different sets: regions-based and edge-based models.

DISCRETE MODELS: MARKOV RANDOM FIELDS AND

GRAPH CUTS

Although discrete optimization models for image segmentation are not really the scope of this thesis, they have become extremely popular in the last few years, mostly due to the popularization of the graph cuts optimization method [START_REF] Boykov | An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision[END_REF]. Despite some inconvenients, these discrete models have proved to be very powerful in particular cases, and are indeed a true alternative to curve evolution based segmentation methods.

Markov Random Field

Markov Random Field (MRF) based image restoration was first introduced by Geman and Geman [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF], and deals elegantly with site interactions. It is based on viewing pixels as sites S = {s 1 , • • •, s N }, and specifying a neighborhood system N = {N s ∈ S } between these sites.

• A random variable X i is associated to each site s i , and the set of all these random variable is a Markov Random Field, defined by the Markov property:

P (X i = x i |X j = x j , i = j) = P (X i = x i |X j = x j , j ∈ N i ) ,
with x i and x j realizations of the random variables X i and X j . This means that the solution at a site i only depends on its neighborhood N i .

• A clique C is defined as a set of mutual neighbor sites, and

C = {C 1 , • • •, C M }
is the set of all possible cliques in the field. A potential function V (C, x) is then defined, associating to each cliques C ∈ C a potential according to the realization x = {x 1 , • • •, x N } of the MRF. This potential specifies interactions between sites inside a clique, and defines the problem to solve.

• The global energy U (x) is expressed as the sum of the potentials of all cliques:

U (x) = C∈C V (C, x) -→ min ,
and the Gibbs distribution of the MRF is defined as:

P (X = x) = Z -1 exp (-U (x)) -→ max ,
where Z -1 is a normalization term.

In summary, besides mathematical definitions, the main idea is to formulate a given problem in term of local interaction between adjacent pixels. If the cliques are pairwise, the image can then be seen as a weighted graph: pixels are seen as nodes, and adjacent pixels are linked by an edge whose weight depends on the local interaction between them. These interactions define the problem, and the energy to minimize is defined as the sum of these weights.

Graph Cuts

The "graph cuts" method is a combinatorial optimization method whose purpose is to minimize a discrete energy defined on a graph. It has been introduced by Greig et al. [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF] for the bi-partitioning segmentation case, and has been somehow forgotten for 10 years, to finally become popular in 1999 [START_REF] Boykov | An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision[END_REF]. Since then, it has been successfully applied to several computer vision problems such as restoration, segmentation, . . . The main idea is to consider the image as a graph G = V, E , where the pixels are seen as a set of nodes V, and adjacent pixels are linked by a set of weighted edges E. Two special nodes are also present : a source node and a sink node, that are in general linked to all pixels of the image.

A cut C ⊂ E is a set of edges removed from the graph, such that there is no more paths connecting the source node to the sink node. A minimal cut is a cut whose cost (sum of all its edges) is minimal. It is then proved that a discrete energy defined on a graph (ie. sum of all the edges of the graph) can be minimized by finding a minimal cut, which is performed by maximizing a flow from the source to the sink (duality min-cut/max-flow) using some well known algorithms from the graph theory domain [START_REF] Ford | Flows in Networks[END_REF][START_REF] Goldberg | A new approach to the maximumflow problem[END_REF].

One of the main advantages of this method is that under certain conditions, a solution corresponding to the global minimum of the energy can be found, which outperforms usual methods finding only local minimas (this may however be an advantage in segmentation). Unfortunately, only a few kind of energies can be minimized this way [START_REF] Boykov | An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision[END_REF][START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], which restricts the range of application of this method.

For example in the case of segmentation, the two subsets of pixels separated by the cut correspond to the two regions for the bi-partitioning segmentation. Various types of segmentation models have been proposed using this approach [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Juan | Active graph cuts[END_REF][START_REF] Boykov | Interactive graph cuts for optimal boundary and region segmentation of objects in n-d images[END_REF][START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF], and this is indeed a true alternative to active contour based segmentation models.

Random Walker

The random-walker method [START_REF] Grady | Random walks for image segmentation[END_REF] is another quite popular segmentation method based on graph representation. It consists in labelling manually a set of pixels, and automatically determining the probability that a random walker starting at each unlabeled pixel will first reach one of the prelabeled pixels. Then, each unlabelled pixel is assigned the label for which the greatest probability is calculated. These probabilities are not computed by simulating a finite set of random walks, but are instead computed by solving a dirichlet problem. This require solving a linear system with graph Laplacian matrix, that can be performed using an efficient GPU implementation allowing interactivity and real time segmentation [START_REF] Grady | Random walks for interactive organ segmentation in two and three dimensions: Implementation and validation[END_REF]. The random walk processes are performed in term of transition probabilities to move from one pixel to another: for a given unlabelled pixel, they consider (infinitely many) random walkers starting from this pixel and all moving according to these transition probabilities.

CONTINUOUS MODELS: SNAKES AND ACTIVE CON-

TOURS

In opposition to discrete optimization models, continuous optimization models consider images as continuous functions, and rely in evolving curves or surfaces 1in images, to finally fit the border of the regions to extract. Evolution of these fronts is done by minimizing a continuous energy defined:

• along the front for boundary based models.

• inside the regions the front delimit for regions based models.

This method was first introduced by Kaas, Witkins and Terzopoulos with the wellknown snake model [START_REF] Kass | Snakes: Active contour models[END_REF], which belong to the boundary-based category.

The Snake Model

Let I : Ω ∈ R 2 → R + be a scalar image. The Snake model is based on evolving a parameterized 2D contour Γ : [0, 1] → Ω by minimizing the following functional:

E(Γ) = 1 0 α|Γ ′ (p)| 2 dp + 1 0 β|Γ ′′ (p)| 2 dp internal terms -λ 1 0 |∇I(Γ(p))| dp external term , (3.1)
where α, β and λ are parameters balancing between the different terms.

This functional is composed of different terms, each of them contributing to the evolution of the contour. As most of the segmentation models, these terms are classified in two categories :

• The first two terms are regularization terms corresponding to internals properties of the contour, and contribute to the desired shape of the segmentation.

• The last term is an external property standing for a data-attachment term, and depends on image characteristics. This energy is intended to be minimum when the contour is smooth enough (internal terms) and overlaps image boundaries corresponding to high image gradient (external term). This model brought several advantages over existing ones, but was quite limited in several aspects and numerous extensions have been proposed.

One of the main limitations is that this model was only restricted to 2D parameterized contours (curve). This type of contour, referred as explicit contour, suffers from several inconvenients that are discussed in the next section.

FRONT REPRESENTATION: EXPLICIT MODELS VS

IMPLICIT MODELS

Before considering the functional defining the segmentation model, one has to choose a type of front representation. Several representations exist in the literature, and each of them has been used in various domains (fluid dynamics, computational geometry, etc. . . ). We can gobally classify these representations in two separate categories: explicit representations and implicit representations.

Explicit Representation

These evolving fronts were initially represented explicitly, which means they were parameterized by a set of values along the front (Fig. 3.1). This representation allowed simple and efficient implementation, but suffered from several inconvenients:

• Topology changes were not allowed due to the parameterized nature of the front.

Several heuristics have been proposed to overcome this issue [START_REF] Mcinerney | Topology adaptive deformable surfaces for medical image volume segmentation[END_REF][START_REF] Mcinerney | T-snakes: Topology adaptive snakes[END_REF][START_REF] Montagnat | A review of deformable surfaces: Topology, geometry and deformation[END_REF].

• A re-parameterization was needed during the contour evolution, that was likely to perturb position of the contour. To solve this issue, different types of basis functions have been proposed in the literature [START_REF] Zienkiewicz | The Finite element method[END_REF][START_REF] William | Numerical Recipes in C[END_REF]5,[START_REF] Ruppert Curwen | Dynamic contours : Real-time active splines[END_REF].

• Generalization to higher dimensions (3D and more) was tedious due to the parameterization.

However, explicit models are often favored when efficiency is preferred to robustness and theoretical correctness [START_REF] Montagnat | A review of deformable surfaces: Topology, geometry and deformation[END_REF], or when open-contours are needed. Recent models based on Delaunay triangulation have also been proposed [START_REF] Pons | Delaunay deformable models: Topology-adaptive meshes based on the restricted delaunay triangulation[END_REF] and are quite promising as they benefit from the robustness of the Delaunay triangulation.

Implicit Representation: The Level-Set Method

In order to overcome the limitations of explicit models, a new way to represent evolving fronts has been introduced with the level-sets method. This method was introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] in 1988, although the same idea was previously used 10 years earlier by Dervieux and Thomasset [START_REF] Dervieux | A finite element method for the simulation of Rayleigh-Taylor instability[END_REF] for fluid dynamics purpose.

This method consists in expressing the interface between two regions as the zero level of a higher dimensional function Φ, defined on the whole image domain. This implicit function is required to be a Lipschitz function, and most of the time is chosen as the signed distance function to the front, such as:

Φ(x) =          -D(x, Γ) x ∈ Γ in 0 x ∈ Γ +D(x, Γ) x ∈ Γ out
where Γ represents the interface (front) separating two regions Γ in and Γ out , and D(x, Γ) the Euclidean distance from the point x to the contour Γ (Fig. 3.2).

Through this representation, geometric properties of the front can easily be expressed as differential terms defined on the level-set function:

• The unit normal vector to the front is -→ n = ∇Φ |∇Φ| .

• The mean curvature of the front is κ = ∇ ∇Φ |∇Φ| .

These properties allow to express front evolution as a time-dependent partial differential equation on the level-set function: consider a front Γ evolving under a speed F in its normal direction -→ n with no tangential velocity, so that ∂Γ ∂t = F -→ n .

As the front Γ corresponds to the zero level-set of Φ, we have, Φ(Γ(t), t) = 0 ∀t , by applying the chain rule, we get

∂ ∂t Φ(Γ(t), t) = ∂Φ ∂t + ∇Φ ∂Γ ∂t = 0 .
In level-sets notation, the evolution equation can be rewritten as,

∂Γ ∂t = F - → n ⇐⇒ ∂Γ ∂t = F ∇Φ |∇Φ| .
Finally, putting all the terms together, we get the evolution equation: In summary, the evolving front is embedded in a higher dimensional function, and its evolution is done through the evolution of this function using a time-dependent partial differential equation.

∂Φ ∂t + F |∇Φ| = 0 . (3.2)
Since it is not necessary anymore to worry about the parameterization, the levelsets framework has naturally become the standard way of expressing contour & surface evolution techniques [START_REF] Osher | Geometric Level Set Methods in Imaging, Vision and Graphics[END_REF]. It brings a lot of advantages compared to explicit models:

• Generalization to higher dimensions is easy and natural.

• Topological changes are automatically handled, since the front is the zero level-set of the level-set function (Fig. 3.2).

• Theoretical correctness based on viscosity solution have been widely studied, providing a good and robust mathematical basis [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF][START_REF] Sethian | Theory, algorithms, and applications of level set methods for propagating interfaces[END_REF].

• Stable and efficient numerical schemes have been proposed to solve the PDE.

• Properties of the curve such as normal vectors, or curvature can be easily and accurately expressed through differential terms.

One of the main inconvenient is that this framework is computationally more expensive than explicit models. Indeed, the evolution of the front is done through the evolution of the level-set function, which is defined on the whole image domain. This requires to compute the solution in the whole image domain, with a computational cost that not depend on the actual size of the front. Several methods have However been proposed in the literature in order to improve its efficiency.

Narrow Band Implementation

Since the front corresponds the zero level-set of the distance function, computing the solution of the PDE (Eq. (3.2)) far from the zero level-set does not contribute to the evolution of the front. In order to decrease the complexity of the method, an alternative is to compute the solution of the PDE only at the neighborhod of the zero level-set, in a "band" thick enough to keep track of the front; this is known as the Narrow Band Approach. This approach decreases dramatically the complexity, and allows computation times close to those obtained using explicit models.

Numerical Schemes

Several numerical scheme have been proposed to solve Eq. (3.2), and the most common choice is to compute spatial derivatives using a finite differences ENO scheme, and time derivative with an explicit scheme. This choice is driven by the fact that finite differences scheme are quite efficient, and combined with the narrow band approach speeds up considerably the level-set method. However finite difference schemes are nonetheless limited in accuracy, and some problems of stability may arise if CFL conditions2 are not satisfied. The chapter 4 investigates a new way to implement the level-set method, using the finite elements method.

Reinitialization

One of the limitations of the level-set method is that the implicit function has to remain a distance function during its evolution, which means it must satisfy a viscosity solution of the Eikonal equation. In general, this is not the case, and the slope of the function may become steep or flat, depending on the speed chosen. The easier way to deal with this limitation is to periodically "reinitialize" the function to a distance function again, without moving its zero level-set (which is not very different from the reinitialization of explicits models). Several methods have been proposed to do so:

• a reinitialization PDE [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] have been proposed to "re-slope" the implicit function, without moving its zero level set:

∂Φ ∂t = sign (Φ) (1 -|∇Φ|) ,
the complexity of this reinitialization depends on the number of iterations to converge, and is in general relatively expensive.

• a re-distancing algorithm, the fast marching method [START_REF] Sethian | Theory, algorithms, and applications of level set methods for propagating interfaces[END_REF], has been proposed to recompute a new distance function from the previous one. It is based on initializing seed points on the zero level-set, and create a new distance map using some incremental updates of an evolving front. This algorithm is quite famous in the computer vision community, as it is an easy and efficient way to compute a distance map in Euclidian or Riemanian spaces, and for example to track geodesics of interest. The complexity of this algorithm is in O(n log n), where n corresponds to the number of points of the grid. However, some algorithms have been proposed to reach a linear complexity in O(n), by using specialized data structures [START_REF] Yatziv | O(n) implementation of the fast marching algorithm[END_REF]. Let us note that the same idea has been used in the fast sweeping method [START_REF] Kao | Fast sweeping methods for a class of static hamilton-jacobi equations[END_REF], which allows a parallelizable scheme.

• The work of Gomes and Faugeras [START_REF] Gomes | Reconciling distance functions and level sets[END_REF] aims at preserving the distance function during its evolution. It is based on computing the solution of the PDE at the points lying on the the zero level-set, and propagating this speed in the whole function thanks to an advection PDE. This method is however expensive, and a re-distancing method is often favored.

In practice, a mix of the first two method is commonly used: a reinitialization PDE is applied on the points very close to the zero level-set, and a re-distancing with a fast marching is performed to recreate a complete narrow band.

In this section, we briefly reviewed methods to represent evolving fronts by first introducing explicit models, and then implicit models through the level-set method. In the next sections, we will focus on segmentation models that will drive the front, and characterized by a functional to minimize. These models can be classified in two distinct categories : Boundary-based models and Region-based models.

BOUNDARY BASED SEGMENTATION MOD- ELS

The main idea behind boundary-based segmentation models is:

1. choose a "boundary model" to characterize boundaries in images.

2. move freely the front in image, and locally stop it when it reaches a boundary.

In general, boundaries are chosen as high image gradient, but several models has been proposed in the literature.

Boundary-based models can be seen as simple "closed edge detectors", and are quite effective on simple objects whose boundaries are well defined. However, these models are "myopic" since they only provide a local criterion to stop the evolution, and do not look at the whole image structure.

The first model of this category was introduced by the snake model, where boundaries were seen as high image gradient. Since then, several extensions have been proposed and the most famous one is without a doubt the geodesic active contour.

Geodesic Active contour

The geodesic active contour model was introduced in [START_REF] Caselles | Geodesic active contours[END_REF] and [START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF][START_REF] Kichenassamy | Conformal curvature flows: from phase transitions to active vision[END_REF], and can be seen as an extension of the classical snake model. Rather than expressing the segmentation functional as a sum of regularization terms and data attachment terms, the authors introduced a functional composed of one single term, as follows:

E(Γ) = L 0 g (|∇I (Γ(s))|) ds = 1 0 g (|∇I (Γ(p))|) ∂Γ ∂p (p) dp , (3.3)
where ds is the arc-length element, and L the Euclidian length of the front Γ.

This functional represent the length of the active contour in a Riemanian space, and minimizing this functional implies finding a contour of minimal length, corresponding to a Geodesic3 . The metric of this space is induced by image characteristics, and has to be strong in homogeneous areas and weak on image boundaries. The function g(|∇I|) define the metric, and is a monotonically decreasing function defined on R + → R + , such that g(0) = 1 and lim x→+∞ g(x) = 0.

A very common choice is given by:

g(|∇I|) = 1 1 + |∇I| r , r = 1 or 2 . (3.4)
Derivative of this functional is done through Euler-Lagrange equations, and minimization of this functional is performed thanks to a steepest gradient descent. This gives us the following contour evolution :

∂Γ ∂t = gκN -(∇g . N ) N , (3.5)
where N denotes the inward normal vector of the front, κ the curvature of the front, and g = g(|∇I|) the metric of the Riemanian space.

The functional (3.3) does not depends on a parameterization any more, and can thus be applied with implicit models, which was the principal limitation of the initial snake model. Moreover, this allows to easily extend the model to higher dimensions. The contour evolution can thus be rewritten in level-set notations:

∂Φ ∂t = g |∇Φ| div ∇Φ |∇Φ| + ∇g , ∇Φ , (3.6)
where Φ represent the level-set function. This evolution is composed of two terms having a different incidence in the evolution of the contour:

• the first term corresponds to a mean curvature motion of the front, that locally stops when g → 0, ie. when the contour reaches high image gradient. It can be seen as a mix between a regularization term (mean curvature motion) and a data attachment term (stop on image gradient). In order to speed up the convergence, a balloon pressure force [START_REF] Cohen | On active contour models and balloons[END_REF] can be incorporated into this term, such that it becomes g |∇Φ| div ∇Φ |∇Φ| + ν , where ν is whether positive (contour shrinks) or negative (contour inflates).

• the second term is a less classical term, and plays the role of an advection term that speeds up the evolution of the front when it is close to an image boundary. Moreover, it is also a refinement term that pushes back the contour on the boundary if it has slightly missed it. Imagine it is a force that make the contour "fall" into the "saddle" created by the metric of the Riemanian space (Fig. 3.3).

Note that previous boundary models were proposed in [START_REF] Caselles | A geometric model for active contours[END_REF][START_REF] David | Computing minimal surfaces via level set curvature flow[END_REF][START_REF] Malladi | A topology independent shape modeling scheme[END_REF], but were not based on an optimization framework. The authors proposed an evolution equation for implicit functions that was similar to the first term of Eq. (3.6). The geodesic active contour model is a more interesting model since it is based on an optimization framework, elegantly incorporate a Riemanian space in the formulation, and its derivation introduces a new term that reinforces the robustness of the model.

Gradient Vector Flow

Other boundary-based models have been introduced, but most of them are relatively similar to the geodesic active contour model (GAC), providing only some slight extensions. However, the Gradient Vector Flow [START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF] model introduces a new approach, and overcomes some limitations of the GAC model. It is based on successive steps, defined as follows:

• First, an edge map f (x, y) is computed similarly to the GAC model, by applying a function g() on the image gradient, except that it has to be large near image boundaries, and low in homogeneous areas (Fig. 3.4, top row).

• A vector field V(x, y) = (u(x, y), v(x, y)) is then obtained from the edge map f (x, y), by minimizing the following functional:

E(V) = µ |∇v| 2 + |∇u| 2 + |∇f | 2 |V -∇f | dx dy .
This vector field has to be smooth when there is no image boundary, and similar to ∇f close to image boundary. Its purpose is to extend the gradient of the edge map in the whole image (Fig. 3.4, top row).

• Finally, an implicit contour is initialized in the image and is driven by the vector flow V(x, y) using the following evolution equation:

∂Φ ∂t = -∇V.∇Φ .
This equation is an advection equation similar to the second term in the derivative of the GAC model. However, the front will be attracted from a further distance due to the smoothness of the field V(x, y) (Fig. 3.4, bottom row).

This model is interesting as it is more robust to initialization than the GAC model. However, it is computationally more expensive since the vector fields computation requires a lot of iterations in order to correctly "extend" the edge-map gradient to the whole image.

In this section, we briefly reviewed boundary-based segmentation models, by first introducing the Geodesic Active Contour Model, and then the Gradient Vector Flow model. These models are quite effective in simples cases, but become quickly limited as soon as image boundaries are not well defined, when it need to be looked at a more global scope. In order to overcome these difficulties, region-based segmentation models have been introduced, and are the scope of the next section.

REGION BASED SEGMENTATION MODELS

Although boundary-based segmentation models are quite powerful to segment simple objects, they suffer from two important limitations:

• These models are based on very local information, and are thus very sensitive to local minima. This is especially a problem when images are noisy: the active contour remains stuck in areas where noise is too important.

• Initialization of the front has to be finely tuned in order to correctly segment the desired object. In most of the cases, the initial contour has to be either completely inside the object, or outside around the object. This is a major drawback especially in the case of segmenting complex 3D structures in medical images.

To overcome these difficulties, region-based models have been introduced. The main idea of these models is to take the whole image into account, by expressing the functionals as integrals over the different regions delimited by the front. In the bi-partitioning case of two regions Ω in and Ω out , these functionals will look like:

E(Γ) = Ω in f in (x) dx + Ω out f out (x) dx , (3.7)
where Ω in and Ω out are two regions partitioning the image domain Ω, such that

Ω in ∪ Ω out = Ω and Ω in ∩ Ω out = ∅.
Γ is the front separating the two regions, and f in and f out are two "penalty functions" defining a pixel-wise energy for both regions Ω in and Ω out .

If the model has to be applied with implicit representations, a common choice is to reformulate the energy in terms of level-set notations, and to derive it in order to obtain a level-set evolution, using a gradient descent:

E(Γ) -→ E(Φ) =⇒ ∂Φ ∂t = - ∂E(Φ) ∂Φ ,
where Φ is the implicit representation of the front Γ. This reformulation can be done very easily in the bi-partitioning case, by assigning a "positive" and a "negative" region, using a regularized heavyside on the level-set function to differentiate between the regions:

E(Φ) = Ω H α (Φ(x))f in (x) + (1 -H α (Φ(x))) f out (x) dx .
The function H α is a regularized heavyside function defined by:

H α (x) =          0 x < -α 1 2 1 + x α + 1 π sin πx α |x| < α 1 x > α
where α is the amount of regularization.

The functional (3.7) can be extended in the case of several regions Ω 1 , • • •, Ω n (multi-phase segmentation), by expressing it as a sum of integral across each region with a corresponding penalty function, such as:

Ω = Ω 1 Ω 2 • • • Ω n , Ω i Ω j = ∅, ∀i = j E = Ω 1 f 1 (x) dx + Ω 2 f 2 (x) dx + • • • Ω n f n (x) dx
However the level-set formulation is more challenging, since it requires several levelset functions to define these regions. A direct approach is to associate a level-set function to each region, by using only their "positive areas". Unfortunately, this approach clearly increase the computational complexity, and some heuristics may be needed to forbid overlapping between the region, and suppress the vacuum between them [START_REF] Zhao | A variational level set approach to multiphase motion[END_REF].

Another approach introduced by Chan and Vese [START_REF] Chan | Active contours without edges[END_REF] is to use the overlapping between level-set to define separate regions. This approach is quite interesting, as it decreases the number of level-set to log 2 (n) in order to model n regions. For example, it would require 2 level-set function to define 4 regions.

In the following, we will focus on the case of the simplified bi-partitioning segmentation, with two regions Ω in and Ω out .

In summary, region-based models consist in choosing adequate penalty functions, that will specify the pixels composing the regions. Several regionbased models have been proposed, but we can again classify them in two categories: those expressed directly through a meaningful functional; and those derived from a Bayesian formulation. However, these two different approaches can lead in some cases to the same functional, under some specific conditions.

Functional-based models

The first category of region-based models are those whose functional is directly expressed in terms of image properties.

Chan-Vese Model

One of the most famous region-based model is the Chan&Vese model (Active contour without edges) [START_REF] Chan | Active contours without edges[END_REF], that gained in popularity mainly due to its simplicity and efficiency. The main idea is to see the segmentation as a piecewise constant partition of the image. A scalar value is associated to each region, and the penalty functions are expressed as the distance between pixels to their corresponding region-wise scalar value. For the bi-partitioning case, this functional is defined by:

E(Γ, µ in , µ out ) = Ω in |µ in -I(x)| 2 dx + Ω out |µ out -I(x)| 2 dx + ν |Γ| ,
and its corresponding level-set formulation:

E(Φ, µ in , µ out ) = Ω |µ in -I(x)| 2 H Φ + |µ out -I(x)| 2 (1 -H Φ ) + ν |∇H Φ | dx ,
where µ in and µ out are the two scalars corresponding to regions Ω in and Ω out .

H Φ = H α (Φ(x))
the Heaviside of the level-set function Φ, and ν the amount of regularization of the front 4 .

In order to minimize this functional, derivatives are performed:

• keeping Φ fixed, derivative of this functional with respect to µ in and µ out is straightforward, and gives us the values µ in (Φ) and µ out (Φ):

             ∂E(Φ, µ in , µ out ) ∂µ in = 0 =⇒ µ in (Φ) = Ω I(x)H Φ dx Ω H Φ dx ∂E(Φ, µ in , µ out ) ∂µ out = 0 =⇒ µ out (Φ) = Ω I(x) (1 -H Φ ) dx Ω (1 -H Φ ) dx .
• keeping µ in and µ out fixed, the derivative of this functional with respect to Φ is done using Euler-Lagrange equations, and gives us the evolution equation:

∂Φ ∂t = - ∂E ∂Φ = δ Φ |µ in -I(x)| 2 + |µ out -I(x)| 2 + νdiv ∇Φ |∇Φ| ,
where δ Φ = δ α (Φ(x)) a regularized Dirac of the level-set function Φ.

The evolution of the front is then performed by alternating between the computation of µ in and µ out which correspond to average means of the two regions, and the level-set evolution. The front will likely evolve in order to separate the bright areas from the dark ones, and will not be influenced by the sharpness of image boundaries. This model is particularly suited to segment objects having a different "color" than the background, but will become limited when both share the same average intensity.

Mumford-Shah functional

The Mumford-Shah (MS) functional [START_REF] Mumford | Boundary detection by minimizing functionals[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] introduced in 1989 was one of the first model proposed to segment regions. Its initial purpose was to perform both segmentation and restoration, by computing a piecewise smooth partition (µ, Γ) of images thanks to the minimization of the following functional:

E(Γ, µ) = Ω |µ(x) -I(x)| 2 dx mean square data term + α Ω-Γ |∇µ(x)| 2 dx smooth regions + ν |Γ| smooth boundaries
, where µ(x) is a piecewise smooth function approximating the image I(x), and Γ a set of discontinuities in µ(x). This functional is similar to the Chan&Vese functional when the middle term is neglected, in the case of piecewise constant function.

Minimization of this functional is quite difficult, and several approaches have been proposed [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF][START_REF] Chambolle | Discrete approximations of the Mumford-Shah functional in dimension two[END_REF]. Two methods using the level-set framework have been proposed independently in [START_REF] Chan | A level set algorithm for minimizing the Mumford-Shah functional in image processing[END_REF] and [START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF]. For the bi-partitioning case, the MS functional is reformulated in level-set notation:

E(Φ, µ in , µ out ) = Ω |µ in (x) -I(x)| 2 H Φ + |µ out (x) -I(x)| 2 (1 -H Φ ) dx + α Ω |∇µ in (x)| H Φ + |∇µ out (x)| (1 -H Φ ) dx + ν Ω |∇H Φ | dx .
with the same notations as before, and µ in and µ out two functions defined on Ω. Similarly to the previous approach, derivatives of the functional are performed:

• keeping Φ fixed, derivative of the functional with respect to µ in (x) and µ out (x) provide two PDEs, that can be solved using the following gradient descent:

           ∂µ in ∂t = α∆µ in + (I -µ in ) in Ω in ; ∂µ in ∂ - → n = 0 on Γ ∂µ out ∂t = α∆µ out + (I -µ out ) in Ω out ; ∂µ out ∂ - → n = 0 on Γ .
• keeping µ in and µ out fixed, the derivative of this functional with respect to Φ

gives the following evolution equation:

∂Φ ∂t = δ α -|µ in (x) -I(x)| 2 + |µ out (x) -I(x)| 2 -α |∇µ in | 2 + α |∇µ out | 2 + νdiv ∇Φ |∇Φ| ,
with the same notations as before.

This approach can be seen as an extension to the previous Chan&Vese model, providing segmentation of piecewise smooth regions instead of piecewise constant regions. However this model is not as robust as the Chan&Vese model concerning front initialization. Indeed, the initial front has to be rather close to the object to segment in order to evolve toward it (the same problem is evoqued in section 5.5.1).

Another inconvenient is that this approach is computationally expensive, since several PDEs need to be solved at each iteration of the level-set evolution. These PDEs are solved using a gradient descent scheme, and therefore the whole minimization process can be seen as gradient descent (PDEs solution) included inside another gradient descent (level-set evolution). The number of iterations for each descent may be rather large, and this can lead to a very expensive scheme.

In Chapter 5, a new segmentation model that minimize the Mumford-Shah functional is proposed, with a computational complexity largely improved using recursive filters to perform the smoothing process.

Bayesian-based models

The second category of region-based models are those derived from a Bayesian formulation. The main idea is no more to define regions as constant or smooth functions, but to define them in term of probability density function. These functions correspond to the probability of a given pixel to belong to the region, and their log-likelihood play the role of the "penalty function" previously used. The maximum a posteriori (MAP) segmentation is then obtained by minimizing functionals similar to Eq. (3.7). You can refer to section 5.2.2 for a more detailed insight about Bayesian formulations.

For example, in the bi-partitioning case, we will obtain functionals of the form:

E(Ω in , Ω out ) = Ω in -log p in (I(x)) dx + Ω out -log p out (I(x)) dx + ν |Γ| ,
where p in (I(x)) and p out (I(x)) are the probability functions for the pixel I(x) to belong respectively to regions Ω in and Ω out . The functional is then reformulated in terms of level-set notations:

E(Φ) = Ω -log (p in (I(x))) H Φ -log (p out (I(x))) (1 -H Φ ) + ν |∇H Φ | dx ,
and evolution equation is obtained by deriving the functional with respect to Φ:

∂Φ ∂t = - ∂E ∂Φ = δ Φ -log p in (I(x)) p out (I(x)) + νdiv ∇Φ |∇Φ| , (3.8)
The main idea is to consider the front as a set of "trial pixels", and to compute for each of them the probability to belong to regions Ω in or Ω out . Each pixel is then assigned to the region having the highest probability, and the front will evolve according to these successive classifications.

Bayesian models differ depending on the probability density function chosen to define regions. Two kind of function can be used: parametric density functions and non-parametric density functions.

Parametric density functions

Parametric density functions rely in expressing distributions with a finite number of parameters. A common, and simple choice is given by the Gaussian distribution [START_REF] Rousson | A variational framework for active and adaptative segmentation of vector valued images[END_REF][START_REF] Rousson | Active unsupervised texture segmentation on a diffusion based space[END_REF]:

p(x|µ, σ 2 ) = 1 √ 2πσ 2 e -(x-µ) 2 2σ 2
, defined by 2 parameters : the mean µ, and the variance σ 2 .

The functional is then expressed according to the parameters of the two regions:

E(Φ) -→ E(Φ, µ in , µ out , σ 2 in , σ 2 out ) ,
and derivative of the functional,

• with respect to µ α and σ 2 α , α ∈ {in, out} gives us the following formulas:

             ∂E(Φ, µ in , µ out , σ 2 in , σ 2 out ) ∂µ α = 0 =⇒ µ α = 1 |Ω α | Ω α I(x) dx ∂E(Φ, µ in , µ out , σ 2 in , σ 2 out ) ∂σ 2 α = 0 =⇒ σ 2 α = 1 |Ω α | Ω α (I(x) -µ α ) 2 dx.
• with respect to the Φ, gives us the following evolution equation:

∂Φ ∂t = δ Φ (I(x) -µ in ) 2 2σ 2 in - (I(x) -µ out ) 2 2σ 2 out + log σ 2 in σ 2 out + νdiv ∇Φ |∇Φ|
The front will then evolve by alternating between computation of region mean and region variance, and level-set evolution. Computation of the mean and variance can moreover be optimized, by incrementally adding or removing contributions of the pixels moving between the two regions, from one iteration to another. This model is simple and efficient, and generalizes the Chan&Vese model [START_REF] Chan | Active contours without edges[END_REF], which can be identified to the case of constant variances σ 2 α = 1 2 . However, using this density implies that regions follow a Gaussian distribution, which may not be true in all cases. Moreover, some numerical problems may arise in the case of almost homogeneous regions, due to division by zero (σ 2 α → 0). In Fig. 3.5, the evolution of the front is performed using the Bayesian Gaussian model.The regions-based segmentation model allows for a more flexible initialization, and a common choice is to initialize the region Ω in as the inside of several circle aligned on a regular grid, and the region Ω out as the background located between the circles. The segmentation finally converges by separating a bright region from a dark one.

Note that if the two initial regions share the same distribution, the front will not move. A slight difference in the distributions of the two regions is required, in order to separate them significantly (consequently, the Chan&vese model as well).

Non-parametric density functions

A more general approach is to consider non-parametric density functions [START_REF] Zhu | Region competition: unifying snakes, region growing, and Bayes/MDL for multiband image segmentation[END_REF][START_REF] Paragios | Geodesic active regions: a new paradigm to deal with frame partition problems in computer vision[END_REF][START_REF] Kim | Nonparametric methods for image segmentation using information theory and curve evolution[END_REF]. The main idea is, given an observed dataset (the regions), to estimate an underlying probability density function that does not belong to a family of parametric function anymore. A common choice is to use a kernel density estimation, known as Parzen windows method. This method estimates the density functions of the regions by "summing" all the pixels inside regions, with a kernel function K(x) applied on each of them. This "summation process" results in convolving region histograms with a kernel function K(x). In the case of generic kernels K(x), the probability density function of a region Ω α is defined by:

p α (x) = 1 |Ω α | y∈Ω α K(x -I(y)) = 1 |Ω α | k∈B H α (k)K(x -k) = 1 |Ω α | (H α * K) (x) ,
where H α (k) is the histogram of the image I(x) inside region Ω α , composed of a finite set "bins" B. The main idea is thus to define probability density functions of regions through their corresponding histograms.

However, a problem concerning the choice of the kernel function remains open. A common choice is to choose a Gaussian kernel, and therefore to express the probability density function as a simple Gaussian convolution of the region histograms. However the variance of these Gaussian kernel ("kernel size") has to be chosen a priori, and its incidence on the segmentation is important. For example, using very large Gaussian kernels somehow generalizes the parametric Gaussian model, as the probability density function become similar to Gaussian functions 5 . Several models have been proposed to automatically choose a correct size of kernel, but they suffer from an important increase in computational complexity.

Differentiating the functional using this kind of probability density function is not straightforward, and in general additional terms appear in the evolution equation [START_REF] Kim | Nonparametric methods for image segmentation using information theory and curve evolution[END_REF][START_REF] Cremers | A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape[END_REF]. In general, these additional terms can be neglected and the evolution is done using Eq. (3.8). Evolution of the front is then performed by alternating between histogram computation, and front evolution.

This model has proved to be very powerfull in a majority of cases, especially on natural images. Discrimination between objects and background through their corresponding histograms appears to be a rather good choice, and this approach has been extensively used in the last few years.

In the geodesic active region model [START_REF] Paragios | Geodesic active regions: a new paradigm to deal with frame partition problems in computer vision[END_REF], Paragios and Deriche proposed to express the segmentation functional as the sum of a region term (non-parametric based model), and a boundary term (geodesic active contour model). This model combines the advantages of both region and boundary models by adding the contributions of two different evolutions in the final evolution equation. But, this model requires some parameters that may be difficult to adjust.

In Fig. 3.6, evolution of the front using a non-parametric model is performed on a simple synthetic image. Probability density functions of the inside and outside regions are computed using their corresponding histograms (Fig. 3.6, bottom row), and the contour evolves until it separates significantly both regions (Fig. 3.6, top row). Note that the parametric Gaussian model would not have worked in this case (consequently, the Chan&Vese model as well), as both regions share the same global mean.

In Fig. 3.7, segmentation of the inside of the skull (brain + cerebrospinal fluid) is performed on a T2 MRI, using a non-parametric model. The non-parametric model is well suited to segment this kind of anatomical structure, as the histogram of the brain is composed of mainly two bright intensities (white for the CSF, and grey for the brain), while the skull is rather dark.

CONCLUSION

In this chapter, we presented an overview of different techniques for image segmentation. A diagram of all these notions can be seen in Fig. 3.8. We first briefly reviewed discrete models for image segmentation, that gained in popularity in the last few years, due to the popularization of the graph-cuts optimization method. These methods exhibit several advantages (globals optimas, fast schemes, numerical sta- Then, we reviewed continuous models for image segmentation which consist in evolving contours in images. Two kinds of contour representations can be used: explicit models that depend on a parameterization; and implicit models introduced by the level-set method, that overcome some important drawbacks but increase the computational complexity.

Evolution of these contours is done through the minimization of a given functional, defining segmentation models. Two different types of models have been considered: boundary-based models, that provide only a very local criteria to stop the evolution of the contour; and region-based models, that provide a more global criteria in order to separate significantly the two regions delimited by the contour.

Region-based models can furthermore be classified in two separate categories: those whose functional is directly expressed in term of image characteristics, and those derived from a Bayesian formulation. However, several links exist between these two categories (Fig. 3

.8, blue curves).

Each model brings its own set of advantages and inconvenients, and no model can pretend to be robust to all alterations present in images (noise, intensity drifts, etc). In chapter 5, a new model based on local statictics is proposed, that aims at unifying both region-based and boundary-based models.

• Clustering

• Graph cuts Optimization • MRF • Discrete • Graphs Discrete Models Optimization Bayesian based Functional based Parametric Models Non-Parametric Models • Gaussian model • Parzen windows model • Level-Sets • Delaunay-based • Fourier descriptors • B-Splines Explicit Models • 3D deformable
Region based models defined on regions

• Gradient Vector Flow

• Geometric/Geodesic active contour

• Snake model

Boundary based

Probability density functions

Continuous Models

Implicit Models

• Finite difference • Reinitialization ... • Finite elements Numerical Scheme "Penalty functions" Continuous • Mumford-Shah • Chan-Vese

Geodesic Active Regions

• Active contours

OVERVIEW

In this chapter, we propose an novel implementation of the level-set method using the finite element method with quadrilateral elements. After briefly explaining the basis of the finite element method and some notions about the type of elements and the discretization, we show how to solve the level-set evolution PDE using finite elements. Similarly to the narrow band approach (Section 3.4.2), the computation of the distance function can be restricted to the elements that are close to the zero level-set, and thus can be represented as a sparse set of connected elements. To keep track of the zero level-set during its evolution, some elements are inserted to the set, and others removed from the set depending on the motion of the curve. Finally some results are shown in 2D and 3D. Keywords: Level-set, Numerical scheme, Finite elements method, Mesh, P1 elements, Q1 elements, Galerkin

INTRODUCTION

The level-set method (Section 3.4.2) is a simple and efficient method used to perform the evolution of a closed contour through a time-dependent partial differential equation (PDE). It thus transform the problem of evolving a curve into the one of solving a specific PDE of the form:

∂Φ ∂t + F |∇Φ| = 0 , (4.1)
where Φ is the level-set function and F the speed of the front in its normal direction. Several numerical scheme have been proposed to solve this PDE [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Barth | Numerical schemes for the hamiltonjacobi and level set equations on triangulated domains[END_REF], and the most common choice is to use the finite difference method (FDM) with an explicit scheme (Section 3.4.2). While being relatively fast and simple to implement, the FDM is also one of the least accurate and least stable numerical scheme. It is often required to take relatively small time step in the evolution of the level-set, to ensure the stability of this method.

Solving the evolution PDE alone does not guarantee that the level-set function remains a distance function during its evolution, which is important to accurately compute geometric properties of the contour such as the curvature or the normal. As a consequence, another PDE has to be alternately solved in order to "reinitialize" the level-set into a distance function again. This equation, known as Eikonal equation, ensures that the slope of the distance function remains 1, and is defined as follows:

|∇Φ(x)| = 1 ∀x ∈ Ω, (4.2)
where Ω the domain of definition of the level-set function Φ.

Solving this equation is equivalent to re-computing a distance map to the active contour, and this can be approximated very efficiently with iterative schemes such as the fast marching method or the fast sweeping method (Section 3.4.2). In summary, usual numerical schemes for the the level-set method suffer from two limitation :

• they are based on the FDM with explicit scheme which often produce instabilities if the time step is too big.

• they have to reinitialize the level-set function each k iterations by performing a reinitialization PDE or an iterative redistancing using the fast marching method.

Several methods have been proposed to overcome these limitations.

In [START_REF] Gomes | Implicit representations of evolving manifolds in computer vision[END_REF], the authors proposed a method to evolve the level-set function without the need of a reinitialization process. It is based on computing the speed of the levelset function only at the points neighboring the zero level-set, and propagating this speed in the entire function by following the normal direction of the front using an advection PDE. Unfortunately, this framework is still based on the FDM, and thus the stability and the convergence of the level-set may not be guaranteed.

In [START_REF] Li | Level set evolution without re-initialization: A new variational formulation[END_REF], the authors formulated the energy of the evolution with an additional term which penalize the deviations of the level-set function from a distance function. The derivation of this energy lead to a new term in the evolution equation which is supposed to preserve the distance function without the need of a reinitialization process. Unfortunately this formulation bring another competitive term in the derivative, whose weight has to be accurately chosen. Moreover, being based on an explicit scheme, this method does not allow to choose large time steps in order to preserve the distance function, thus increasing the total number of iterations.

In [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF] and [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF], the authors proposed a stable and efficient numerical scheme for the level-set method based on the finite elements method (FEM) instead of the FDM. They solved both level-set evolution equation and Eikonal equation simultaneously, using an implicit scheme which allow to choose large time step in the level-set evolution without having to reinitialize it. These new method bring several advantages while introducing some new drawbacks. These are investigated in this chapter.

In the following, we explain some basis of the finite elements method, and how it can be used to implement the level-set method. Following [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF] and [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF], we propose a novel framework for implementing the level-set method using the FEM with quadrilateral elements (Q 1 ) instead of the triangular elements (P 1 ) used in these works.

FINITE ELEMENT METHOD

Several methods exist to find approximate solutions of PDEs. The finite difference method (FDM) is one of the most common choice, and consist in expressing the solution at discrete points on a regular grid and approximating differential operators as finite differences between neighbor points. It has the advantage of being extremely simple to implement and computationally fast, but it is known to be the least stable and the least accurate method. Others methods such as finite volumes method or finite elements method are a little more complicated to implement, but give significantly better results. f

The finite elements method [START_REF] Zienkiewicz | The Finite element method[END_REF] is a powerful numerical scheme for approximating the solutions of PDEs. This method is based on partitioning the domain of the solution into several elements of arbitrary forms (polygons in 2D, polyhedra in 3D), thus leading to the creation of a mesh composed by vertice and elements. The solution of the PDE is then computed at the vertice positions, and interpolated inside each elements, providing a continuous piecewise polynomial approximation of the solution. Because finite elements methods can be adapted to problems of great complexity and unusual geometry, these are extremely powerful tools in the solution of important problems in heat transfer, fluid mechanics, and mechanical systems.

Let's see step by step how to use this method in order to solve a PDE of the form:

       A(u(x)) = 0 ∀x ∈ Ω u(x) = f (x) Dirichlet or ∂u ∂ -→ n (x) = f (x) Neumann ∀x ∈ ∂Ω
where u is the solution to approximate, A denotes any linear operator of the function u, Ω the domain of definition of u, ∂Ω the boundary of Ω, and f a smooth function defined on ∂Ω. A border condition is necessary in order to correctly define the function u, and various type can be chosen (Dirichlet, Neumann, . . . ).

Variational Formulation

The first step of the FEM is to rephrase the original problem into its weak, or variational form. This reformulation is based on the following functional analysis theorem :

f (x) = 0 ∀x ∈ Ω ⇔ Ω f (x)g(x) dx = 0 ∀g ∈ H 1 0 (Ω) .
If a function f equal zero in a domain Ω, then its scalar product with any differentiable function g defined on Ω equals zero 1 . The second part is called the weak formulation of the first part.

In our case, we get

A(u(x)) = 0 ∀x ∈ Ω ⇔ Ω A(u(x)) w(x) dx = 0 ∀w(x) .
The purpose of this notation is to transform the initial PDE into a system of equation, by choosing a finite number of independent test functions w (each of them creating a new equation). Then, the main idea is to discretize the solution u by using a finite number of values, and finding these values by solving a system of n equations (induced by the different test functions w) with n unknowns (the unknown representing the solution u).

In general, PDEs with high differential orders are difficult to solve, as they impose some constraints on the polynomial order of the discretization basis of the solution 2 . Fortunately, the weak formulation allows to decrease the differential order of the equations, by applying some integral transformations based on integration by parts. For example, the Green theorem consists in transforming an integral of a second order differential function by two integrals of first order differential functions, and is defined as follows:

Ω -(∆a b) (x)dx = Ω (∇ a∇ b) (x)dx - ∂Ω ∇ a(x), n(x) b(x)dx (4.3)
where a and b are any functions defined on a domain Ω, ∂Ω referring to the border of Ω, and n(x) the normal vector to Ω at point x oriented toward the exterior (Fig. 4.1).

Moreover, this transformations allows to include the border condition of the PDE in the weak formulation, by replacing the correct values in the boundary integral. In summary, little to no computation is usually required for this step, and all the transformations are done by hand on paper.

Discretization of the solution

The second step is the discretization of the solution. It raises the following question: how to approximate a continuous function u defined on a specific domain Ω with a finite set of discrete values ?

The most common way is to partition the domain Ω into a mesh (Fig. 4.2). Only the values at the vertice position of the mesh are sufficient to approximate the solution. These values are then interpolated inside the elements of the mesh, and the solution becomes a piecewise polynomial approximation of u.

Ω Ω

Choice of the type of elements

In general, the partition of the domain is done by using an unique type of element, but it's not necessary. Any type of element can be chosen, but this choice is generally limited between two types: elements of type P n and elements of type Q n .

• Elements of type P n (Fig. 4.3) are triangles in 2D, tetrahedra in 3D, and more generally n-simplex in nD. These elements bring several advantages concerning the mesh generation, as usual triangulation methods (Delaunay) can be used to perform the mesh. Therefore, these elements are the most adapted if the problem is defined on a geometrically complex domain.

Elements P 0

Elements P 1 Elements P 2 • Elements of type Q n (Fig. 4.4) are squares in 2D, cubes in 3D, and more generally hypercubes in nD. Depending on the problem, these elements can facilitate the mesh creation. For example, if the mesh must have its vertice positioned on a regular grid, the mesh creation using quadrilateral elements becomes trivial.

Elements Q 0 Elements Q 1 Elements Q 2 Figure 4.4: Elements of type Q n .
The 'n' in the P n and Q n notations denotes the polynomial order in each dimensions of the approximation inside the elements. It is linked to the number of nodal points (discretization points of the function) related to the element. These nodal points are in general located on the element boundary, but it's not necessary (Fig. 4.3, 4.4). In the following, we will essentially focus on P 1 and Q 1 elements, where each nodal points correspond to a vertex of the mesh.

Basis functions and interpolation

To interpolate the discrete values of u at the nodal points inside the elements of the mesh, basis functions are introduced. These functions, often called "hat functions" due to their shape, are functions associated to each nodal points of the mesh. The i th basis function ϕ i associated to the i th vertice V i is defined as follows:

• ϕ i : Ω → R continuous on Ω • ϕ i (V j ) = δ ij (1 if i = j, 0 if i = j)
In the case of P 1 and Q 1 elements, these functions are respectively piecewise affine and piecewise bilinear (Fig. 4.5) on the elements of the mesh. Notice that the gradient of these function in the Q 1 case is not constant per element, unlike in the P 1 case (These functions are called "pagodas", due to their particular shape). The purpose of these functions is to define a basis of the space where the piecewise polynomial approximations are defined.

Therefore, the piecewise polynomial function û approximation of the continuous function u on Ω can be defined as follows:

u(x) ≈ û(x) = n i=1 u i ϕ i (x) ∀x ∈ Ω ,
where i are the vertice of the mesh, u i the values of the solution at vertex i, and ϕ i (x) the basis function associated to the i th vertex.

In summary, this step consists in choosing a type of element (ie. shape of element with polynomial order), partitioning the domain into a mesh using this kind of element, and computing the basis functions associated to each vertex of the mesh.

Discretization of the weak formulation

The third step consists in discretizing the weak formulation. This commonly means that if u is discretized with n values, finding a set of n independent test functions for the weak formulation in order to transform the PDE into a system of n equation with n unknown. Several methods exists (points-subdomain collocation [START_REF] Zienkiewicz | The Finite element method[END_REF], etc...), but the most common discretization scheme is the Galerkin discretization: the basis functions associated to each vertex are chosen as the test functions in the weak formulation. Finally, if A is a linear operator, we obtain the linear system:

A(u(x)) = 0 =⇒                Ω A(u 1 ϕ 1 (x) + . . . + u n ϕ n (x)) ϕ 1 (x) dx = 0 Ω A(u 1 ϕ 1 (x) + . . . + u n ϕ n (x)) ϕ 2 (x) dx = 0 . . . Ω A(u 1 ϕ 1 (x) + . . . + u n ϕ n (x)) ϕ n (x) dx = 0
Basically, each of these equations denotes the local deviation of the solution at the neighborhood of a vertex in the mesh. This leads to a system of n equations with n unknowns, and solving this system3 will give values u 1 , u 2 , . . . u n of the solutions at specific locations. Interpolating these values inside the elements will finally lead to the piecewise polynomial approximation.

FINITE ELEMENT IMPLEMENTATION OF THE LEVEL SET METHOD

Let's see how to implement the level-set method using the finite element method.

Formulation of the Problem

Following [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF] and [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF], we use the finite element method in order to approximate the solution of two equations :

• the level-set evolution equation, for any normal speed

F ∂Φ ∂t + F (Φ) |∇Φ| = 0 , (4.4)
• the Eikonal equation

(∇Φ) 2 = 1 . (4.5)
By solving simultaneously these two equations, the level-set function will remain a distance function during its evolution, and no auxiliary reinitialization scheme will be needed. Solving the Eikonal equation is more like a constraint that must be satisfied, thus there is no need of border condition for this PDE. Notice the difference between the reinitialization PDE of the level set (Section 3.4.2), that is formulated as a time dependent PDE and this formulation : in the framework of Section 3.4.2, one has to alternate between evolution and reinitialization; in our case, both evolution and reinitialization are done simultaneously.

Time Discretization : implicit scheme

The time derivative appearing in equation 4.4 has to be discretized in a special way. We choose to discretize it using an implicit scheme: this means that the time derivative is approximated using first order finite difference, thus introducing a time step; and the spatial derivatives appearing in the PDEs are calculated at the end of the time step. Compared to explicit schemes, where spatial derivatives are calculated at the beginning of the time step, implicit schemes allow to choose much larger time step ∆t, and are numerically more stable. However, implicit schemes involve the coupling between unknowns, and therefore require solving a linear system.

Actually, due to our formulation, this implicit scheme is required in order to satisfy the Eikonal equation at the end of the time step. Using an explicit scheme would not have allowed us to include the distance constraint.

So by including time as a parameter, and replacing the continuous function Φ by its discrete approximation u, Eq. 4.4 becomes:

u(x, t + ∆t) = u(x, t) -∆t F (u(x, t + ∆t)) |∇u(x, t + ∆t)| , (4.6)
where the normal speed F is calculated at time t + ∆t.

Thanks to the constraint imposed by Eq. 4.5, the gradient appearing in the second term of Eq. 4.6 has to be 1, and thus can be removed. Moreover, the "incremental function" that has to be added to the function u at time t in order to get u at time t + ∆t is actually the unknown of this equation. By expressing this unknown as a variable v, we get:

u(x, t + ∆t) = u(x, t) + v(x, t) , with v(x, t) = -∆t F (u(x, t) + v(x, t)) .
By expressing spatial derivatives at the end of the time step, Eq. 4.5 becomes:

(∇ (u(x, t) + v(x, t))) 2 = 1 . (4.7)
Finally, the system of Eq. 4.4 and Eq. 4.5 becomes:

{

• the level-set evolution equation consists in finding the "step function"

v v + ∆tF (u + v) = 0 , (4.8)
• the Eikonal equation calculated at the end of the time step

(∇u + ∇v) 2 -1 = 0 . (4.9)
Solving this system will give us the "step" v, that makes the level-set evolve while preserving it as a distance function. Note that the normal speed F is expressed at the end of the time step.

Variational Formulation

The next step is to write the weak formulation of equations 4.8 and 4.9 and apply some integral transformations to reduce the differential order.

The weak formulation of Eq. 4.9 does not need an integral transformations, since it is a first order differential equation. 4 . Therefore, we get:

Ω (∇u + ∇v) 2 w - Ω w = 0 ∀w . (4.10)
The weak formulation of Eq. 4.8 depends on the normal speed F applied of the level-set. We choose to use the Gaussian region based model from section 3.6. By expressing F at time t + ∆t, we get:

F (u + v) = ǫ∇. ∇ (u + v) |∇ (u + v)| + (I -µ in (u)) 2 σ 2 in (u) - (I -µ out (u)) 2 σ 2 out (u) . ( 4

.11)

The first term corresponds to the regularization term with a weight ǫ, while the second and third terms stand for data attachment terms. We choose to express the data attachment terms at the beginning of the time step in order speed up the gradient descent. Consequently, our model is a semi-implicit model, and we can note these data attachment terms as a constant α precomputed at the beginning of each timesteps. Moreover, due to the distance constraint imposed by Eq. 4.7, the regularization term can be simplified as a Laplacian.

Therefore, the normal speed (Eq. 4.11) becomes:

F (u + v) = ǫ ∆(u + v) + α , (4.12)
and the weak formulation of Eq. 4.8 is then defined as

Ω v w -∆t ǫ Ω ∆(u + v) w + Ω α w = 0 ∀w . (4.13)
This weak formulation is problematic since it includes a second order differential term (Laplacian), which prevents us from using first order polynomial elements such as P 1 or Q 1 elements. In order to decrease the differential order of this equation, we apply the Green theorem (Eq. 4.1). Thus, the weak formulation becomes:

Ω v w -∆t -ǫ Ω ∇(u + v) ∇w + ǫ ∂Ω ∇(u + v) , n w + Ω α w = 0 . (4.14)
where ∂Ω represents the boundary of the domain Ω, and n is the normal to Ω oriented toward the exterior.

Discretization of the Solution & Algorithmics

The discretization of the solution involves choosing a domain for computing the solution, and partitioning this domain using an appropriate type of element.

P 1 elements
In [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF] and [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF], the authors restricted the domain of computation to a thin band located at the neighborhood of the zero level-set, similarly to the narrow band approach (Section 3.4.2). To keep track of the zero level-set during its evolution, they introduced some dynamics to insert elements or remove elements from the band. They used regular P 1 elements (triangles in [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF], and tetraedra in [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF]) to partition this band, thus approximating the level-set as a piecewise affine function (Fig. However, by using such P 1 elements, the authors introduce some problems:

• By using anisotropic meshes, some bias may appear in the evolution of the level-set. Indeed, as the Galerkin discretization measures the local deviation of the solution at the neighborhood of a vertex, the accuracy of the solution is directly linked to the "shape" of this neighborhood. If this neighborhood is not isotropic, the final solution may be affected by this anisotropy, and we may see some slight perturbations in the solution. For example, the evolution of a circle under a mean curvature leads (in the intermediate steps of the evolution) to an ellipsoid oriented by the "mesh aspect" (Fig. 4.8). • While the triangulation of a regular grid is trivial in 2D as each square produce 2 triangles, it becomes more complicated in 3D, where each cubes produce 6 tetraedra (Fig. 4.9). This difficulty is even more problematic in our case, as the mesh has to evolve incrementally by adding and removing elements at each iterations.

Figure 4.9: 3D Triangulation of a cube: creation of 6 tetraedra per cube.

Q 1 elements

To avoid these problems, we choose to discretize the solution using Q 1 quadrilateral elements. This choice seems more natural in both 2D and 3D, and while it requires more calculations (Cf. Appendix A), it overcomes some important drawbacks of the previous methods [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF] and [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF]:

• because quadrilaterals tile the image space, there is no need to introduce arbitrary simplicial subdivisions of the image elements (pixels or voxels), thus leading to a simpler method.

• As the vertice of the mesh are in direct correspondence with image pixels, the standard multilinear interpolation scheme can be directly used to provide the continuous interpolant of the discrete solution.

• Using these elements correct the bias induced by anisotropic triangular meshes. We use the same dynamics as in [START_REF] Weber | Sparse finite elements for geodesic contours with level-sets[END_REF] and [START_REF] Weber | Sparse finite element level-sets for anisotropic boundary detection in 3d images[END_REF], by computing the solution only in a band defined at the neighborhood of the zero level-set, and evolving this band by adding and removing elements to keep track of the zero level-set.

However, adding an element in the band involves extrapolating some values in order to have a good initialization for the next iteration. This extrapolation is performed using a fast marching (Section 3.4.2), but in a simpler way as only a few elements are added at each iterations.

Thanks to this discretization, we can express the discrete approximation u(x, t) of the level-set function Φ, as well as the step function v(x, t) by:

u(x, t) = n i=1 u i (t) ϕ i (x) , (4.15)
and

v(x, t) = n i=1 v i (t) ϕ i (x) , (4.16)
where u i (t) and v i (t) are the values located at the i th vertex at time t of the level-set function and the step function ; and ϕ i (x) the basis functions associated to the i th vertex and defined on its incident Q 1 elements (Fig. 4.6).

Discretization of the weak formulation

The next step is to discretize both of the weak formulations Eq. 4.10 and Eq. 4.14, by choosing a set of n orthogonal tests functions w(x). We choose the Galerkin discretization, which means that we take the basis functions ϕ i (x) as test functions. This leads to the following two systems of equations.

Eikonal equation

The Galerkin discretization of Eq. 4.10 gives us the system:

Ω (∇u + ∇v) 2 ϕ i - Ω ϕ i = 0 i ∈ {1, 2, . . . , n} (4.17)
By expressing u and v through their linear combination of basis functions, we get:

∇u = ∇
(u a + v a ) ∇ϕ a , ∇ϕ b (u b + v b ) , therefore Ω (∇u + ∇v) 2 ϕ i = Ω n a=1 n b=1 (u a + v a ) ∇ϕ a , ∇ϕ b (u b + v b ) ϕ i = n a=1 n b=1 (u a + v a ) Ω ∇ϕ a , ∇ϕ b ϕ i (u b + v b ) = ( u + v) T Q i ( u + v) ,
and finally, Eq. 4.17 can be reformulated as:

( u + v) T Q i ( u + v) - Ω ϕ i = 0 i ∈ {1, 2, . . . , n} , (4.18)
with

u =     u 1 . . . u n     , v =     v 1 . . . v n     and Q i n×n matrix defined by Q i ab = Ω ∇ϕ a , ∇ϕ b ϕ i

Evolution equation

Similarly, the Galerkin discretization of Eq. 4.14 gives us the system:

∀i ∈ {1, 2, . . . , n} Ω v ϕ i -∆t -ǫ Ω ∇(u + v) ∇ϕ i + ǫ ∂Ω ∇(u + v) , n ϕ i + Ω α ϕ i = 0 ,
and this system can be reformulated as:

P i v -∆t ǫ -R i + B i ( u + v) + P i α = 0 i ∈ {1, 2, . . . , n} , (4.19)
with

u =     u 1 . . . u n     , v =     v 1 . . . v n     and                P i vector of size n defined by P i a = Ω ϕ a ϕ i R i vector of size n defined by R i a = Ω ∇ϕ a , ∇ϕ i B i vector of size n defined by B i a = ∂Ω ∇ϕ a , n ϕ i
Thanks to the Galerkin discretization, we can compute the step function v by solving the two systems of equations (4.18) and (4.19).

Solving the system

In order to solve both of the systems 4.18 and 4.19, we use the conjugate gradient method [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. This method is the most popular iterative method for solving large systems of linear equations of the form

Ax -b = 0 , (4.20)
where x is the unknown, b is a known vector, and A is a known, square, symmetric positive-definite matrix. In our case, we first have to reformulate our systems, by expressing them in a similar form as Eq. 4.20. The unknown is the "step function" v that needs to be computed in order to update the level-set function. Therefore, we have to solve 2 systems of the form:

   A 1 v -b 1 = 0 A 2 v -b 2 = 0 (4.21)
where A 1 and A 2 are respectively the matrices of systems 4.18 and 4.19 and b 1 , b 2 the corresponding known vectors. The system 4.18 is not a linear system, as quadratic terms of the unknown v appears in the equations. However, since the deviation is almost affine in the unknown v, the quadratic term v T Q i v can be removed, and the system can indeed be linearized. Therefore, the system becomes:

2 u T Q i A i 1 v + u T Q i u - Ω ϕ i -b i 1 = 0 ∀i ∈ {1, 2, . . . , n} , (4.22)
where A i 1 is a vector corresponding to the i th row of the global matrix A 1 and b i 1 a scalar corresponding to the i th component of the known vector b 1 . By developing the system 4.19, we find

P i -∆tǫ -R i + B i A i 2 v -∆tǫ -R i + B i u -∆tP i α b i 2 = 0 , (4.23)
where A i 2 is a vector corresponding to the i th row of the global matrix A 2 and b i 2 a scalar corresponding to the i th component of the known vector b 2 .

Unfortunately, both of the matrix A 1 and A 2 are not symmetric, and therefore the conjugate gradient method can not be applied directly. Therefore, instead of solving the linear system 4.21, we solve the following system:

   A T 1 (A 1 v -b 1 ) = 0 A T 2 (A 2 v -b 2 ) = 0 (4.24)
and finally, by combining these two system:

A T 1 A 1 + A T 2 A 2 v -A T 1 b 1 + A T 2 b 2 = 0 (4.25)
Matrices A 1 and A 2 have to be recomputed at each iteration, because they both depend on the level-set function u, constantly updated. Finding the elements composing matrix A 1 and A 2 is purely combinatorial, and actually a simple visit on the elements of the mesh is necessary to compute these values.

In summary, the system 4.25 needs to be solved at each iterations of the level-set evolution, in order to compute the step v that is needed update the level-set function. This step will then make the level-set evolve, while preserving it as a distance function. The computation of this step is performed by using the efficient conjugate gradient method, and therefore is computed in O(n).

The algorithm below explains the different steps of our implementation. It can be seen as a gradient descent performed inside another gradient descent: the inside one is meant to compute the step v needed to evolve the level-set function, while the global one is the evolution of the level set function, which is done by minimizing a given functional.

RESULTS AND COMPARISONS

Actually, as the data attachment term in the evolution equation are computed at the beginning of the time step, our finite element implementation only manages to compute accurately the curvature of the level-set, with a distance preservation constraint. Therefore the following results present mean curvature evolutions with distance preservation using Q 1 elements. In Fig. 4.13 can be seen the mean curvature motion of a cube using a finitedifference implementation with no reinitializations (top row), and our Q 1 finiteelement implementation (bottom row). The finite-difference implementation does not provide intrinsic sub-voxelic representation of the distance function, but some methods exist to extract sub-voxelic iso-contours, such as the marching cube algorithm. By using finite-differences with no reinitializations, the different level of the level-set functions can be seen to get away from each others, meaning that the distance functions is not preserved during its evolution. Moreover, finite-differences are known to be less accurate than finite-elements, and we can indeed see the isocontours extracted being less regular than the finite element-implementation. Our Q 1 finite-element implementation performs the evolution very well, by providing a good sub-voxelic representation of the solution at very low resolution, and indeed preserves the distance function during its evolution. In Fig. 4.14 can be seen the good preservation of the distance function during its evolution, with several level-sets staying at equal distance from each others. Without the distance constraint, the level-set of the functions would get away from each other at the corner of the cube, due to higher curvature in these areas. Moreover we can see that, at very low resolution (only a few pixels), we manage to represent very well the zero level-set of the distance function.

CONCLUSION & DISCUSSION

In this chapter, we presented a novel and stable implementation of the level-set method using Q 1 finite elements. This implementation has the advantage of being more accurate than finite-differences based implementation, and provides a good subvoxelic representation of the zero level-set without the need of auxiliary iso-surface extractions methods. This implementation is optimized by restricting the domain of computation to a thin band located at the neighborhood of the zero level-set. This band is constantly updated to keep track of the front, and this allows to omit useless calculations at positions far from the zero level-set.

Thanks to the Green theorem, the differential order of the evolution equation has been decreased: the computation of the second order differential term has been transformed into first order differential terms. This decreases the error in the approximation by preventing numerical computations of second order derivatives, with all the problems brought by the discretization (ENO schemes, etc...). Moreover, this transformation allows to decrease the thickness of the narrow band. Indeed, finite difference implementations requires a computational stencil relatively large for computing sec- ond order derivatives (stencil of 5 pixels wide), and thus influence the thickness of the narrow band. In our case, no computational stencil is needed, and all the derivatives are computed by integration inside elements. This allows to choose a much thinner band, with a minimum effective thickness of 1 element. By using an implicit scheme instead of an explicit one, this implementation allows to choose larger time steps in the evolution, while guaranteeing the stability. Moreover, the use of an implicit scheme allows to incorporate a distance constraint in the evolution, in order to preserve the distance function of the level-set at the end of the time step. This allows to skip the auxiliary reinitialization process in the level-set evolution which may perturb the zero level-set, such as reinitialization PDE or fast marching (Section 3.4.2).

By using Q 1 elements instead of P 1 , some bias introduced by the anisotropy of the P 1 elements is corrected, and the sub-voxelic representation of the solution is improved. Moreover, this kind of element facilitate the partition of the space, and the evolving mesh embedding the zero level-set can be updated more easily, with adding new cubes in the mesh by extrapolating some values.

However, this implementation also bring some drawbacks. By using an implicit scheme, a large system of equation need to be solved at each iterations of the levelset evolution. While being done by using an efficient conjugate gradient descent, this is of course a bottleneck in the implementation compared to usual finite difference based implementation. However, while the conjugate gradient method ensures to find the solution of the system in n step, one possible solution would be to stop the descent before convergence, depending on some properties of the quadratic form (energy) of the system. But this introduces some new problems, such as detecting if the approximation is good enough to stop the descent.

Moreover, using an implicit scheme is a bit incoherent with having a thin narrow band ... The implicit scheme allows to choose large time steps in the evolution, while having a thin narrow band prevents from using large time steps in order to keep track of the zero level-set. This implicit scheme is only used to include the distance constraint, and the good stability of this scheme with large time step can unfortunately not be used with thin bands. Of course, one can indeed increase the thickness of the band and choose large time steps, but in this case the conjugate gradient descent computed at each iteration would be very slow, as much more equations (one per vertex) would appear in the system.

While the distance constraint allows us to skip reinitialization process of the level set, a fast marching is indeed needed at each iterations to extrapolate the new values in the narrow band. This fast marching is less expensive than usual fast marching performed in common finite difference implementations, but requires the same initialization and data structures which are not negligible.

At last, the sub-voxel representation may not be very important, as the grid of image pixel is relatively dense, and the common marching cube method performs iso-surface extraction quite well.

This approach brings some good advantages, but is also quite difficult to implement and requires a lot of calculus in order to accurately compute a Laplacian with distance preservation. It indeed compute accurately the regularization term of the level-set, but common finite-difference methods do not actually suffer from the stability of this computation. Image segmentation is a difficult problem, and the data attachment term in the evolution of the level-set plays a much more important role than accurately computing the regularization term.

OVERVIEW

In this chapter, we propose a fast and robust segmentation model based on local informations. It combines local statistics of the regions along the contour (inside and outside) to drive the segmentation procedure. Since it is based on local decisions, this model is more robust to local variations of the regions of interest (contrast, noise, blur, . . . ). In order to avoid the computational burden of a direct estimation, we express these local statistics as the result of Gaussian convolutions. This makes an efficient implementation via recursive filters possible, and gives a complexity similar to methods based on global statistics. The model generalizes several state-of-the-art segmentation models, thanks to a neighborhood parameter whose role is investigated in this chapter. This model nonetheless introduces some new difficulties which are inherent to the fact of basing a global property (the segmentation) on pure local decisions. The second part of this chapter explores some of those difficulties and propose some possible corrections. Finally, an efficient implementation is proposed, using the computational power of highly parallel devices such as modern GPUs.

INTRODUCTION

Several segmentation models based on curve evolution have been proposed, each of them bringing its own advantages and its own drawbacks (Cf. Chapter 3). Activecontours based segmentation models can be classified in two distinct categories: edgebased models and region-based ones.

Edge-based models consist in evolving a contour in homogeneous areas, and locally stopping it when it reaches high image gradients [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Kichenassamy | Gradient flows and geometric active contour models[END_REF][START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF]. These models have several advantages like local coherence and robustness to region inhomogeneities, but they also have important drawbacks that make them inefficient on noisy images [START_REF] Cremers | A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape[END_REF], or when the contour initialization is not completely inside or outside the region to segment. Moreover, they are only able to segment regions with rather sharp edges, and this can result in "leaks" when the region edges are too diffuse.

Region-based models consist in evolving a contour in order to separate significantly the regions it delimits [START_REF] Chan | Active contours without edges[END_REF][START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF][START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF][START_REF] Paragios | Geodesic active regions: a new paradigm to deal with frame partition problems in computer vision[END_REF][START_REF] Cremers | A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape[END_REF]. For example, separating bright areas from dark ones, or extracting areas matching a certain histogram [START_REF] Kim | Nonparametric methods for image segmentation using information theory and curve evolution[END_REF]. These models take account of the whole region and are thus more robust to noise and to initialization than edge-based models. However, they are inefficient in handling images with regions which statistics are spatially varying across the domain. and contour between white matter and grey matter may be locally diffuse.

One difficulty comes from the fact that segmentation models often relies on global criteria1 , whereas images are most of the time subject to local variations (in intensity, contrast, noise . . . ). In opposition, it is admitted that the human visual system is more sensitive to local variations of the information [START_REF] Wohrer | Model and large-scale simulator of a biological retina with contrast gain control[END_REF], and the Adelson illusion is one of the many example showing this property (Fig. 5.2). This explains why it is often difficult to create segmentation models that are able to extract regions matching our perception.

Medical images are a typical example of image that suffer from alterations that makes usual segmentation models most of the time inefficient. For example, MRIs coming from high field MR devices suffer from biases and intensity drifts, and therefore prevent the direct use of classical region-based models. In addition, MRIs can be very noisy depending on the acquisition parameters, and the partial volume effect coming from the voxelic discretization often produce locally diffuse contours between tissues. These two alteration make classical edge-based models inefficient, as the contour may be stuck in regions with high level of noise, while some leaks may appear in the segmentation where contours are too diffuse (Fig. 5.1).

To overcome these problems, new segmentation models have recently emerged to segment regions by computing their statistics locally, in the neighborhood of the active contour [START_REF] Piovano | Efficient segmentation of piecewise smooth images[END_REF][START_REF] An | gamma-convergence approximation to piecewise smooth medical image segmentation[END_REF][START_REF] Brox | On the statistical interpretation of the piecewise smooth Mumford-Shah functional[END_REF][START_REF] Brox | On local region models and a statistical interpretation of the piecewise smooth mumford-shah functional[END_REF][START_REF] Benoit Mory | Variational Segmentation using Fuzzy Region Competition and Local Non-Parametric Probability Density Functions[END_REF][START_REF] Lankton | Hybrid geodesic region-based curve evolutions for image segmentation[END_REF][START_REF] Li | Minimization of region-scalable fitting energy for image segmentation[END_REF][START_REF] Li | A variational level set approach to segmentation and bias correction of medical images with intensity inhomogeneity[END_REF][START_REF] Taron | Border detection on short axis echocardiographic views using a region based ellipse-driven framework[END_REF]. All these models are somehow similar, and have been developed in parallel. These models bring several advantages like robustness to noise and inhomogeneities, but introduce some new drawbacks that are investigated hereafter.

In this chapter, we propose such a model a segmentation model based on local statistics. It consists in computing regions statistics differently for each point of the contours, inside a Gaussian window of fixed size at the neighborhood of the active contour (Fig. 5.3). This model introduces a new parameter which corresponds to the size of the Gaussian neighborhood, and its effect on the segmentation is very interesting: when a large neighborhood is chosen, the model is equivalent to standard region-based models, whereas a small neighborhood results in a model that behaves like standard edge-based ones. Moreover, choosing an intermediate neighborhood generalizes the well known Mumford-Shah functional, as it provides segmentation for piecewise smooth regions. Thus, this model can be seen as an elegant unification of edge-based and region-based models. Experimental results show that this model is well suited for medical images, especially on those suffering from high intensity drifts.

But basing a segmentation on local decisions brings a whole new set of difficulties: locally the information might be insufficient to make the contour evolve, and it might need to be looked at different scales through a scale space approach; local decisions in neighboring regions might be contradictory, as the notions of inside and outside might be different depending on the initialization and position of the active contour. The second part of this chapter explores these problems and describes two rules that help in defining a meaningful segmentation procedure based on local contrast information.

Finally, the third part of this chapter describes an efficient implementation for computing these local statistics, based on Gaussian recursive filters, and highly parallel devices such as GPUs (Graphic Processing Units).

REGION SEGMENTATION USING LOCAL STATIS- TICS

Piecewise Smooth Regions and Local Means

The proposed model can be seen as a local extension to the Chan & Vese model. Instead of minimizing the distance between regions and their corresponding mean value, we propose minimizing the distance between regions and their local means, different for each point and computed inside a Gaussian window (Fig. 5.3). These local means are expressed through Gaussian convolutions, and correspond to smoothed versions of segmented regions.

Assume that the image is covered by a set of regions Ω i , i = 1...n. Denoting by µ i (x) the local mean of the region Ω i at point x, this quantity is defined by:

µ i (x) = Ω i g ρ (x -y)I(y)dy Ω i g ρ (x -y)dy , (5.1)
where I is a given image and g ρ a Gaussian kernel of standard deviation ρ.

The denominator of this expression is a normalization factor which is important for the points that are close to the region boundaries, when the Gaussian kernel does not overlap completely with the regions (Fig. 5.4).

Based on the Chan&Vese model, this model consists in minimizing the distance between the regions and their corresponding smooth approximation, as follows:

E(Ω 1 , • • • , Ω N ) = i Ω i ((I(x) -µ i (x)) 2 dx + ν|C| , (5.2)
where µ i (x) is the local mean of the region Ω i at point x, and |C| the total length of the interface between regions, which stand for a regularization term. The derivation of this energy can be found in Appendix B, and gives us the following gradient descent2 , for the bi-partitioning case in level-set notation:

                       µ in (x) = Ω in g ρ (x -y)I(y)dy Ω in g ρ (x -y)dy µ out (x) = Ω out g ρ (x -y)I(y)dy Ω out g ρ (x -y)dy ∂Φ ∂t (x) = (I(x) -µ in (x)) 2 -(I(x) -µ out (x)) 2 + ν∇. ∇Φ |∇Φ| . (5.3)
where µ in (x) and µ out (x) denote respectively the local means of the inside and outside regions delineated by the zero level of the level-set function Φ.

The minimization of this energy is then performed by alternating contour evolution and local mean computation. Note that the contour is only driven by local means computed at its neighborhood, and those computed inside the region have no incidence in the evolution. This raises some problems which are discussed further.

In Fig. 5.5, comparisons between Chan & Vese model (top), and our model (bottom) is done through the segmentation of the liver from a 2D anatomical MRI. Evolution of the contour can be seen, as well as their corresponding piecewise constant segmentation for the Chan & Vese model, and piecewise smooth segmentation for the local where |K| denote the total length of the interface, and ν the regularization parameter.

Several simplifications need to be done for decomposing the term P (I|K). Assuming that there is no correlations between regions Ω i , as well as independence between pixels inside each region, we can write,

P (I|K) = P (I|Ω 1 ) P (I|Ω 2 ) • • •P (I|Ω n ) = x∈Ω 1 P (I(x)|Ω 1 ) x∈Ω 2 P (I(x)|Ω 2 ) • • • x∈Ω n P (I(x)|Ω n )
where P (I(x)|Ω i ) corresponds to the probability of the pixel I(x) to belong to the region Ω i .

In order to maximize the probability Eq. (5.5), we apply the negative logarithm and minimize the obtained energy using usual gradient descent. Therefore, we get

P (K|I) -→ max ⇐⇒ -log (P (K|I)) -→ min (5.7) with -log (P (K|I)) = E = i Ω i -log(P (I(x)|Ω i ))dx + ν|K| . (5.8)
The choice of the distribution P (I(x)|Ω i ) depends on the type of regions desired in the segmentation. Several distributions have been proposed 3 , and a classical choice is to model regions using a global Gaussian distribution so that:

P (I(x)|Ω i ) = 1 √ 2Πσ i exp -(I(x) -µ i ) 2 2σ 2 i (5.9)
where µ i and σ 2 i are respectively the mean and the variance of the region Ω i . This distribution generalizes the Chan&Vese model, as the energy obtained is similar when

σ 2 i = 1 2 .
In order to model regions through local statistics, we choose to use local Gaussian distributions. It takes the position into account, and assumes that each pixel in the region has a different probability, induced by its own neighborhood. This distribution is defined as:

P (I(x), x|Ω i ) = 1 √ 2Πσ i (x) exp -(I(x) -µ i (x)) 2 2σ i (x) 2 (5.10)
where µ i (x) and σ i (x) 2 are respectively the local mean and local variance of the region Ω i at position x.

These local statistics can be computed in a Gaussian windows around each pixel, and thus expressed as Gaussian convolutions:

             µ i (x) = Ω i g ρ (x -y)I(y)
dy

Ω i g ρ (x -y)dy σ i (x) 2 = Ω i g ρ (x -y) (I(y) -µ i (x)) 2 dy Ω i g ρ (x -y)dy
.

(5.11)

In summary, the probability for each pixel to belong to a given region is modeled through a Gaussian distribution, locally computed in a spatially Gaussian window.

Putting all the term together, the energy to minimize becomes,

E = i Ω i (I(x) -µ i (x)) 2 2σ i (x) 2 + 1 2 log(2πσ i (x) 2 )dx + ν|K| , (5.12)
and its derivatives can be found in the same way as previously. This leads to the following gradient descent for the bi-partitioning case,

∂Φ ∂t (x) = (I(x) -µ in (x)) 2 2σ in (x) 2 - (I(x) -µ out (x)) 2 2σ out (x) 2 + log σ in (x) 2 σ out (x) 2 + νκ , (5.13)
where κ is the curvature of the front, and µ in (x), µ out (x) and σ in (x) 2 , σ out (x) 2 are respectively the local means and local variance of the inside and outside regions.

This model is similar to the previous piecewise smooth model, with the addition of the local variance. It is useful to model the local variation of image intensity, such as local noise, or proximity to an image edge. However, this term might be problematic if the size of the Gaussian windows is too small. Indeed in homogeneous areas, the local variance has to be close to zero, and therefore dividing by this term in the evolution may produce some numerical instabilities.

In Fig. 5.6 can be seen a template for computing these local statistics at each point of the contour. The local means are first computed by blurring each region separately, and then normalized by a blurred version of regions characteristic functions. Local variances are then computed in the same way using the previously computed local means. See section 5.6 for more details about the implementation.

VARYING THE NEIGHBORHOOD

Varying the size of the Gaussian neighborhood impacts significantly the behavior of the curve evolution. However, several generalizations of other segmentation models can be made. 3 cases are investigated :

• small neighborhood: the variance of the Gaussian window is close to zero.

• intermediate neighborhood: any Gaussian window of medium size.

• large neighborhood: the variance of the Gaussian window tends to infinity.

Small Neighborhood

Since the size of the Gaussian window is close to zero, the local variance has to be zero, and this is thus problematic. Therefore, we consider in the following the piecewise smooth model, without the local variance. We have, for the bi-partitioning case of two regions Ω in and Ω out ,

∂Φ ∂t = (I(x) -µ in (x)) 2 -(I(x) -µ out (x)) 2 + ν∇ ∇Φ |∇Φ| , (5.14)
with

µ i (x) = Ω i g σ (x -y)I(y)
dy

Ω i g σ (x -y) dy and g σ (x) = 1 √ 2Πσ exp -x 2 2σ 2 .
(5.15)

When σ → 0, we have:

lim σ→0 g σ (x) = δ(x) , (5.16)
where δ(x) denote the Dirac function.

Therefore, for a point x lying on the boundary of the open subspace Ω i ,

lim σ→0 µ i (x) = Ω i δ(x -y)I(y)
dy

Ω i δ(x -y) dy = lim y→x I(y) , y ∈ Ω i = lim α→0 I(x + α - → n ) ,
where -→ n denote the normal to Ω i , oriented toward the interior.

Therefore, as the regions Ω in and Ω out are complementary, we get:

lim σ→0 ∂Φ ∂t (x) = lim σ→0 (I(x) -µ in (x)) 2 -(I(x) -µ out (x)) 2 + ν∇ ∇Φ |∇Φ| = lim α→0 I(x) -I(x -α - → n ) 2 -I(x) -I(x + α - → n ) 2 + ν∇ ∇Φ |∇Φ| = (∆-→ n I) 2 -(∇-→ n I) 2 + ν∇ ∇Φ |∇Φ|
where -→ n is the normal to the region Ω in , and ∆-→ n and ∇-→ n denotes respectively the backward and forward difference derivative in the normal direction -→ n .

In continuous areas, both of the backward and forward derivatives are similar, therefore the level-set function Φ will only evolve under a mean curvature motion. If we are infinitesimally close to a discontinuity in the normal direction, one of the two derivatives has to be very large, and the other one can then be neglected. This leads to

∂Φ ∂t (x) =              ν∇ ∇Φ |∇Φ| x ∈ continuous areas fmax (∇-→ n I) 2 , -(∆-→ n I) 2 + ν∇ ∇Φ |∇Φ| x ∈ discontinuity
where 'fmax' denotes the maximum of the absolute value.

Finally, by taking only one derivative in presence of discontinuity, the directional derivative can be rewritten in term of scalar product:

∂Φ ∂t =              ν∇ ∇Φ |∇Φ| x ∈ continuous areas - → n , ∇I 2 + ν∇ ∇Φ |∇Φ| x ∈ discontinuity
where ∇I is the image gradient.

If we assume that image edges correspond to discontinuities, the behavior of the model when σ → 0 is equivalent to standard edge-based models.

Indeed, the term appearing in presence of a discontinuity is strongly similar to the advection term coming from the geodesic active contour, or the gradient vector flow close to image edges: the curve is "attracted" by the areas where the image gradient is very strong, and then remains stuck in the "discontinuity", oscillating around the image edge.

Therefore, with very small Gaussian neighborhoods, the piecewise smooth model is somehow equivalent to an edge-based model: there are no data attachment term contribution in homogeneous areas 4 , and the contour remains stuck at image edge positions.

Intermediate Neighborhood

Choosing an intermediate neighborhood allows the segmentation of piecewise smooth regions. However, one famous model has been defined for segmenting piecewise smooth region : the Mumford-Shah model [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. This model aims at finding a piecewise smooth approximation of a given image by minimizing the given functional:

E(µ, C) = Ω (I -µ) 2 dx + λ Ω\C |∇µ| 2 dx + ν|C| (5.17)
where µ is the piecewise smooth function, C the interface between smooth regions, I the original image, and λ and ν are constant weighting parameters. The minimization of this energy gives the piecewise smooth approximation u, as well as the interface Γ between regions.

As depicted in chapter 3, the minimization of this functional can be done using an active contour framework thanks to the level-set method. However, this minimization is computationally expensive, since a Poisson PDE needs to be solved for each region at each iteration of the contour evolution.

In [START_REF] Brox | On the statistical interpretation of the piecewise smooth Mumford-Shah functional[END_REF], the authors highlighted a link between the piecewise smooth model (sect. 5.2.1), and the Mumford Shah functional. The main idea is to express the piecewise smooth regions through the minimization of an energy instead of the result of a convolution. Based on the work of [START_REF] Nielsen | Regularization and scale space[END_REF][START_REF] Nielsen | Regularization, scale-space and edge detection filters[END_REF][START_REF] Yuille | A computational theory for the perception of coherent visual motion[END_REF], they found the following link:

µ i (x) = Ω i g λ (x -y)I(y)
dy

Ω i g λ (x -y) dy with g λ (x) = exp (-λx 2 ) E(µ i ) = Ω i   (I(x) -µ i (x)) 2 + ∞ k=1 λ k k! j 1 +j 2 =k d k µ i ∂x j 1 ∂y j 2 (x)   dx -→ min .
By "inserting" this energy into the piecewise smooth energy, it can be rewritten as:

E(K) = i Ω i ((I(x) -µ i (x)) 2 dx + ν|C| E(µ, K) = i Ω i   (I(x) -µ i (x)) 2 + ∞ k=1 λ k k! j 1 +j 2 =k d k µ i ∂x j 1 ∂y j 2 (x)   dx + ν|C| where K = {Ω 1 , • • •, Ω N } and µ = {µ 1 , • • •, µ N }.
Finally, neglecting the penalizers of order k > 1 which have a lower incidence, yields

E(µ, K) = i Ω i (I(x) -µ i (x)) 2 + λ |∇µ i | 2 dx + ν|C|
which states exactly the Mumford Shah functional.

Therefore, with intermediate Gaussian neighborhood, the piecewise smooth model is somehow equivalent to Mumford Shah functional: both energies are similar under some simplifications, and the piecewise model greatly improves the minimization process by filtering each region with an efficient Gaussian filter scheme (see section 5.6) instead of solving a Poisson PDE for each region and at each iteration.

Large Neighborhood

Using a large neighborhood makes the model behave like usual global regions models. There are no more locality appearing in region statistics, and each point of the contour shares a unique global distribution for each region of the segmentation. Using a Gaussian distribution generalizes the Chan&Vese model, and any global distributions can be used to represent regions. However, using large neighborhoods makes the model behave like global models, and it thus looses its robustness to inhomogeneities, and its good locality properties.

COHERENCE IMPROVEMENT

Unfortunately, the local models presented in the previous sections introduce new difficulties. One of these difficulties comes from a problem of coherence in the segmentation: as there is no more one single distribution for a given region, the notion of regions is not well defined. Distant points of the contour evolve rather independently, and there is no more coherence in the segmentation of regions. In the following, we focus on the Bayesian model using the local Gaussian distribution (with the local variance).

Coherence Problem

For each point of the contour, a local distribution for the inside region and a local distribution for the outside region is estimated. The segmentation model aims at separating these two distributions, and thus distant parts of the contour evolve independently. Later in the evolution, these incoherences lead to the creation of "singularities" corresponding to unwanted local stable minima. When a singularity is formed, no topological change is possible, and the segmentation cannot be turned into something coherent anymore. This is illustrated in Fig. 5.7: the red circle segments the interior of the dark object, while the green circles segment its exterior according to their initialization. No topological change are possible between these two growing regions, and a singularity is indeed created. Other circles collapse under the effect of regularization, because of the absence of contribution due to the data attachment term in homogeneous areas. Note that this problem also occur with simpler initializations, as soon as the contour crosses the edges of the object to segment (Fig. 5.8).

Coherence Solution

A lot of works have been done to avoid the kind of local minima depicted in Fig. 5.7. In [START_REF] Juan | Stochastic motion and the level set method in computer vision: Stochastics active contours[END_REF], the authors added some noise in the evolution to perturb the active contour, and thus to prevent it from being stuck in local minimas. Unfortunately, this stochastic evolution is not well adapted to our problem, as it would create as much incoherence in the evolution.

In [START_REF] Charpiat | Generalized gradients: Priors on minimization flows[END_REF], [START_REF] Sundaramoorthi | Sobolev active contours[END_REF], the authors smoothed the evolution term of the active contour along the contour, in order to add some correlation in the evolution of neighbor points. This allows the contour to evolve in a more coherent way, avoiding local minima. Unfortunately, the local minima we tackle cannot be avoided in this way, as the evolution can be seen as coherent until a singularity is created. However, if we make the assumption that the objects to segment are characterized by boundaries where image intensity changes "in the same way" from the inside to the outside 5 , we can use this constraint to maintain the coherence in the segmentation (Fig. 5.8.c). This constraint can be defined by having, for each point x of the contour, the local interior mean µ in (x) either inferior or superior to the local exterior mean µ out (x). Thus, if the contour is locally in contradiction to this constraint, the evolution direction needs to be changed. This can be explained by looking at the energy of the model. Actually, the final segmentation corresponds to a local minimum of the energy, which strongly depends on the initialization of the contour. If the initial contour is close to a local minimum, then the gradient descent will minimize the energy by making the contour converge toward this minimum. However, if this local minimum is not the one desired, one should instead try to go away from it, in order to find another one. In our case, these local minima are characterized by the local means of the inside and outside regions for all points of the contour. Fig. 5.9 shows the energy restricted to a point x of the contour Γ. According to its initialization, the contour must evolve in order to minimize the energy by separating the two local regions with µ in (x) > µ out (x) (red path). However, if we impose the constraint of having µ in (x) < µ out (x), then the contour has to evolve in the opposite direction, thus increasing the energy. When the constraint of having µ in (x) < µ out (x) is indeed satisfied, the gradient descent start to minimize the energy by making the contour converge toward the desired local minimum (green path). This behavior can be seen as a local simulated annealing, as the energy is locally increased in order to avoid an undesired local minimum. Therefore, with the incorporation of this behavior, the evolution equation becomes for the bi-partitioning case:

µ in (x) > µ out (x) µ in (x) < µ out (x) E(Ω in , Ω out )| x x ∈ Γ Initialization
             ∂Φ ∂t (x) = S(x) (I(x) -µ in (x)) 2 2σ in (x) 2 - (I(x) -µ out (x)) 2 2σ out (x) 2 + log σ in (x) 2 σ out (x) 2 + νκ S(x) = ± sign (µ in (x) -µ out (x))
where '±' depends on the type of the global constraint. The energy is therefore increased or decreased, depending on the sign of S(x). Note that this simulated annealing can be used in a similar approach with multiple regions, by increasing the energy if the local features of the contour separating two different regions are not the ones desired.

In Fig. 5.10, the convergence of the active contour can be seen with and without this coherence constraint, on a synthetic image with high intensity drifts. Without the constraint (Fig. 5.10, top row), the segmentation obtained is "locally correct", as it locally separates the regions of the interior and exterior of the contour. However, it has no sense if we look at the global scope, and this kind of segmentation is of course not acceptable. With the constraint of having the inside region locally darker than the outside region (Fig. 5.10, middle row), the contour finally manages to segment the whole grey circle, even under strong inhomogeneities. By looking at the energy of the evolution with coherence constraint (Fig. 5.10,bottom), the simulated annealing slightly increases the energy at the beginning, and that's why the evolution of the top-right circle is somehow not natural. However, the energy starts decreasing when the evolution become more natural, and finally reaches a minimum of lower energy than without the coherence constraint.

Absence of evolution in locally homogeneous areas is natural since the local statistics of the inside and outside regions are similar. This is another problem introduced by the use of local statistics, and a solution to this problem is proposed in the next section.

MULTISCALE IMPROVEMENT

Homogeneous Regions Problem

It is easy to understand that there is no data term contribution in homogeneous areas, as the local statistics in the region inside and outside are the same (see Eq. (5.13)).

At convergence, the small circles in homogeneous areas collapse because of the regularization term, which leads to some kind of phantom contour, resulting in "segmentation layers" that are of course not desired (Fig. 5.11.bottom). These layers come from the fact that the only contribution of the data attachment term is close to image edges, where the difference of local region statistics make the contour evolve.

However, it can be observed that the "thickness" of these segmentation layers depends on the size of the Gaussian windows of the model. With a small Gaussian windows, the layers will be thin, whereas a large windows will result in a wider layer. Unfortunately, increasing the size of the Gaussian windows makes the model behave like global models, and it thus looses its good locality properties.

In Fig. 5.11, the left blue circles collapse because they are too far from an image edge. This edge proximity is linked to the neighborhood parameter of the model, and a larger neighborhood would have indeed made these circles evolve. Note that the right circles are collapsing under the effect of the data attachment term, as the "middle layer" influence their local statistics.

Multiscale Approach

As mentioned previously, the thickness of the "segmentation layers" is correlated with the neighborhood parameter of the model, but using large neighborhoods give results equivalent to standard global models (section 5.3.3).

Thus, we would like to use small neighborhoods close to an image edge, and large neighborhoods in homogeneous areas, depending on the proximity of the closest image edge. The question is how to detect this edge proximity.

In [START_REF] An | gamma-convergence approximation to piecewise smooth medical image segmentation[END_REF], the authors introduced a scale selection framework in a local statistics model by using two given scales of local statistics, and combining these scales differently depending on the proximity to an image edge. To detect this proximity, they used a simple edge detector similarly to edges based models. This approach is interesting, however using an edge detector may not give accurate scales for locally diffuse contours.

In [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF], the authors tackled similar problems with the purpose of extracting edges and ridges from images. Their method is based on detecting features at various scales of the same image, and keeping the k most salient features depending on a saliency measure. It allows them to extract sharp edges as well as smooth ones by increasing the locality of the estimations of these features.

Our problem is somehow similar. If we assume that the neighborhood parameter is actually a scale used to compute local statistics, the problem becomes the one of finding the most salient scale for each point of the contour to make the con-tour locally evolve, even in locally homogeneous areas. Actually, the local model of the previous sections is not restricted to have a unique size of Gaussian window for all points, and we can indeed use local window of different size depending on local features. As we want to favor sharp edges over smooth ones, and make the contour evolve in locally homogeneous areas, the scale is chosen as the smallest one inducing an evolution speed superior to a given threshold :

ρ(x) = inf ρ ∂Φ ∂t ρ (x) > ǫ x ∈ Γ , (5.18)
where ∂Φ ∂t ρ

is the evolution equation computed with a Gaussian windows of variance ρ, x a point of the contour Γ, and ǫ a constant corresponding to a threshold for the scale selection.

By choosing this scale selection scheme, the contour will evolve even in homogeneous areas. Indeed, at each iteration and for each point, the optimal scale is chosen by slowly increasing the neighborhood from the minimum scale to the maximum, until it reaches some edge of an object, thus giving information about how the contour should locally evolve.

Unfortunately, we removed a parameter (the windows size) by introducing another parameter (the scale selection threshold). However the speed of the level-set is somehow normalized by the local variance Eq. (5.13), and therefore the range of the speed induced by the model is independent of the nature of the images to segment. Experimentations have shown that by fixing this threshold to '1', the behavior of the model is correct, and all the points of the contour evolve even in homogeneous areas. Figure 5.12: Top: initialization and evolution of the contour; Bottom: Scales chosen at each point of the contour, with ǫ = 1 and using 30 scales. Figure 5.12 shows the evolution of the active contour (top row), as well as the different scales chosen at each point (bottom row). This image is interesting as it shows three fundamental problem in medical image : large intensity drifts (top corner of the triangle), noisy regions (right corner of the triangle) and smooth contours (left corner of the triangle). Where standard edge-based or region-based model would fail to provide good results, the multiscale local model behave quite well, choosing accurately the scales depending on the proximity and the "smoothness" of the triangle edges. The robustness of the model to noisy regions comes from the fact that the speed (Eq. (5.13)) is normalized by the local variance, so the noisier the region is, the smaller the speed is, and the higher the scale chosen will be. This allows to choose larger scales if the contour is in a noisy regions, thus improving the accuracy of computation of local statistics.

EFFICIENT IMPLEMENTATION

Parallelized Finite-Differences based Level-Sets

These local-statistics models are implemented using a finite-difference based levelset method with a narrow band approach. This choice is driven by the fact that finite-difference based numerical schemes are very easy to implement with acceptable accuracy, and are highly parallelizable. Indeed, we parallelized this method using the OpenMp library 6 . This library allows us to efficiently distribute the speed computation for each point of the narrow band on the different cores of modern multicore CPUs. The reinitialization of the level-set is performed at each iteration using a fast-marching method, optimized to have a complexity in O(n) (section 3.4.2).

Recursive Gaussian filter

All the local statistics used in the model are Gaussian convolutions, and can be computed very efficiently using recursive Gaussian filters.

Several recursive filters have been proposed to perform Gaussian convolutions [START_REF] Tan | Performance of three recursive algorithms for fast space-variant gaussian filtering[END_REF]. We chose the Deriche recursive filter [START_REF] Deriche | Recursively implementing the gaussian and its derivatives[END_REF] for its efficiency and highly parallelizable implementation. This filter is based on approximating the causal part (i.e. x 0) of the Gaussian function by a function of the form:

h + σ (x) =k a 0 cos w 0 σ x + a 1 sin w 0 σ x exp - b 0 σ x + k c 0 cos w 1 σ x + c 1 sin w 1 σ x exp - b 1 σ x
where k, a 0 , a 1 , b 0 , b 1 , c 0 , c 1 , w 0 , w 1 are constants depending on the Gaussian parameter. The same approximation is done for the anti-causal part (i.e. x < 0). The z-transform of these functions is known to be an exact ratio of polynoms, which leads to the following notation for the causal part with a third order approximation:

H + σ (z) = n + 0 + n + 1 z -1 + n + 2 z -2 1 + d + 1 z -1 + d + 2 z -2 + d + 3 z -3 , (5.19)
where the coefficients n + i and d + i are functions of the parameter k, a 0 , a 1 , b 0 , b 1 , c 0 , c 1 , w 0 , w 1 . The same transformation is done for the anti-causal part.

The convolution of an input f (x) by a kernel h(x) in the spatial domain is written as a multiplication of their z-transform in the z-domain [START_REF] Deriche | Recursively implementing the gaussian and its derivatives[END_REF]. Therefore, using the previous z-transform of the causal and anti-causal parts of the Gaussian approximation leads to the following recursive equations:

               g + (x) = n + 0 f (x) + n + 1 f (x -1) + n + 2 f (x -2) -d + 1 g + (x -1) -d + 2 g + (x -2) -d + 3 g + (x -3) g -(x) = n - 0 f (x) + n - 1 f (x + 1) + n - 2 f (x + 2) -d - 1 g -(x + 1) -d - 2 g -(x + 2) -d - 3 g -(x + 3) . (5.20)
where g + (x) corresponds to the result of the left-to-right recursion, and g -(x) to the result of the right-to-left recursion. The complete result of the filtering g(x) is the sum of the two recursions:

g(x) = g + (x) + g -(x) .
Due to the separability of the Gaussian function, an image can be filtered by a 2D Gaussian kernel in 2 pass, by first filtering all the rows separately, then filtering all the column separately (Fig. 5.13). This lead to a computational complexity in O(n), where n is the number of pixels.

In addition, this scheme is highly parallelizable since each of the rows and columns of the image can be filtered independently, with for each of them, the backward and forward recursion being done simultaneously. On highly parallel devices, this can lead to significant gains in speed.

Let's see how many convolutions are needed to compute the local statistics for the bi-partitioning case. By noting (. * .) | Ω the convolution operator on a domain Ω we have,

• The local mean is actually a ratio of convolutions:

µ i (x) = Ω i g ρ (x -y)I(y)
dy

Ω i g ρ (x -y)dy = (g ρ * I) | Ω (x) (g ρ * 1) | Ω (x)
where 2 convolutions are needed: the numerator corresponds to a filtered image of the region Ω, and the denominator to a blurred characteristic function of Ω.

• The local variance can be computed using the previous local mean:

σ i (x) = Ω i g ρ (x -y) (I(y) -µ i (x)) 2 dy Ω i g ρ (x -y)dy = Ω i g ρ (x -y)I(y) 2 dy Ω i g ρ (x -y)dy -2µ i (x) Ω i g ρ (x -y)I(y)
dy

Ω i g ρ (x -y)dy + µ i (x) 2 = Ω i g ρ (x -y)I(y) 2 dy Ω i g ρ (x -y)dy -µ i (x) 2 = g ρ * I 2 | Ω (x) (g ρ * 1) | Ω (x) -µ i (x) 2
Only one additional convolution is needed which corresponds to a filtered squared image of the region Ω.

Therefore, the total number of convolutions needed for computing the local statistics of a given region is 3. For the bi-partitioning case, the total number of convolutions needed is 6. However, an extra convolution can be saved by expressing the blurred characteristic function of the complementary region (g ρ * 1) | Ω by 1 -(g ρ * 1) | Ω , but special border conditions are needed in this case.

Scale Selection using GPU

With the multiscale approach, various scales are incorporated in our model, so each region has to be blurred several times using different Gaussian kernels to obtain a scale space of local statistics. This allows to retrieve, for each point of the contour, the optimal scale using the criteria of Eq. (5.18).

For efficiency purposes, as the local statistics need to be recomputed at each iteration, the scale space is limited to only a few scales. In the experiments, the Our model was tested on several synthetic and real images in 2D and 3D.

In Fig. 5.14, our model is compared to the global Gaussian model (Chan & Vese model), on a synthetic image created as a simple cortex model with large inhomogeneities. The kind of initialization used here does not facilitate the global model as the contour crosses the edges of the structure to segment. Moreover, the distribution of intensities of the internal and external regions are the same. Due to its sensitivity to inhomogeneities, the global model (Fig. 5.14 top row) results in separating globally the bright areas from the dark ones, with no preference to sharp edges. Moreover, global non parametric models such as Parzen windows [START_REF] Kim | Nonparametric methods for image segmentation using information theory and curve evolution[END_REF] give similar results.

By contrast, our model (Fig. 5.14 middle row) behaves quite well. It extracts accurately the details of the region's shape, it is more robust to inhomogeneities, and improves notably the behavior of the model to initialization. Note that the scales ( Fig. 5.14 bottom row) are chosen accurately, depending of the proximity to the closest image edge. Due to the proximity of the initialization to the solution, only a few scales are needed. As can be seen, the segmentation without coherence constraint produce some singularities that correspond to unwanted local minima. These singularities are visible on the piecewise smooth image as well, and slightly increase the energy of the segmentation.

As pointed out in section 5.3.2, the local model is similar to the Mumford Shah functional, that's why the result obtained on this image is strongly similar to the results obtained by minimizing the Mumford Shah functional. This gives us a partition of the image into piecewise smooth regions.

However, minimizing the Mumford Shah functional using our model gives us local minimas of the functional, due to active contour framework. This can be an advantage as it allows to give an initialization of the evolving contour, and to focus on a region in particular.

A previous segmentation model has been proposed in [START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF] to minimize the Mumford Shah functional using an active contour framework. However, this model is computationally expensive as it requires to solve a Poisson PDE inside each region, and at each iteration. Our model greatly improve the complexity, by filtering each region with a complexity in O(n).

In Fig. 5. [START_REF] Caselles | Geodesic active contours[END_REF], we applied the model on structures of interest in retina. These structures are quite difficult to segment, as they suffer from locally diffuse contour, and inhomogeneities coming from the imaging device. Moreover, large blood vessels crosses these structures, with similar image intensity. As can be seen the local model behave quite well by locally separating the inside region from the outside.

We incorporated in the model a rather simple shape prior, penalizing the points of the contour that are far from the centroid of the region, and thus preventing the contour from being attracted by the blood vessels. This shape prior is formulated in a Bayesian form, by maximizing the probability of a point to be close to the centroid of the region it belongs. This term alone result in the creation of circles, and coupled with the data attachment term prevent the deviation of the region from an regular circle.

The model has been implemented in 3D, in order to create a complete 3D head model from anatomical MRI, and is the basis of the next chapter.

CONCLUSION

In this chapter, we proposed a segmentation model based on local statistics. This model is positioned between classical boundary models (very local criteria), and classical region models (very global criteria). Since it is based on local decisions, this model is more robust to local variations of the regions of interest (in contrast, noise, blur . . . ). In order to avoid the computational burden of a direct estimation, these local statistics are expressed as the result of Gaussian convolutions, with a fixed Gaussian kernel.

The incidence of the size of the Gaussian windows has been analyzed, and we pointed out that our model generalizes several state-of-the-art segmentation models:

• With a small Gaussian window, the model is somehow equivalent to standard edge based models.

• an intermediate Gaussian window is similar to minimizing the Mumford-Shah functional.

• using a large Gaussian window is equivalent to usual region based models.

The model nonetheless introduce some new difficulties which are inherent to the fact of basing a global property (the segmentation) on pure local decisions. We emphasized two problems in the model, and proposed solutions to both:

• a problem of coherence in the evolution, that have been overcome by incorporating a constraint using a local simulated annealing. Although this constraint restricts the nature of objects to segment, it is quite effective on objects that satisfy it. Actually this constraint is not so restrictive, as the majority of structures of interest in medical images satisfy it.

• an absence of evolution in locally homogeneous areas, that have been overcome by incorporating a scale-space approach in the computation of the speed for all the points of the contour. The scale selection is a computational burden, but an efficient implementation using highly parallel GPU is used.

At last, we proposed an efficient implementation using recursive Gaussian filters to compute the statistics, and efficient parallel library (OpenMP and Cuda), to fully use the potential of modern multicore CPUs and GPUs. Future work includes using other local distributions of intensities, such as nonparametric distributions. However, incorporating the local variance in the model produced some numerical instabilities if the Gaussian window is too small. Therefore using locally complex distribution may not be meaningful, especially if the area of computation of local statistics is too small.

Another extension is to investigate the advantages of using more than two regions to segment images. However, our local formulations does not include one single distribution per region, so two regions may be sufficient to segment any image composed of several different structures ( Fig. 5.15). It has been shown that 4 regions are sufficient to segment any kind of images due to the "four color theorem" [START_REF] Chan | A level set algorithm for minimizing the Mumford-Shah functional in image processing[END_REF].

We proposed a simple scale selection which gives more importance to sharp edges, but other scale selections can be used depending of the type of segmentation desired. An optimal scale coming from the derivative of the energy with respect to the scale could be formulated. However, it would certainly require to look at the neighborhood of each point of the contour, hence our scale selection method is not so different.

Finally, an extension of this model to 3D vectorial images is proposed in the next chapter, using registered T1 and T2 MRIs.

OVERVIEW

This chapter presents a method for the creation of anatomical head model from bimodal T1 and T2 MRIs. After briefly reviewing the anatomical structures present in the head, a framework is proposed to segment successively the different anatomical structures included in the model by using the local-statistics model introduced in the previous chapter. The scalp, the CerebroSpinal Fluid, the White Matter, the Grey Matter and the skull are segmented consecutively, with a topological constraint is proposed for the segmentation of the Grey Matter.

Keywords: Segmentation, Head model, Topological constraint, Bimodal MRIs

INTRODUCTION

Creation of complete and detailed models of the head is extremely useful in the study of brain functionality and anatomy. It consists in automatically delimiting a set of anatomical structures from medical images of the head. The main purpose of these models is to be used for simulating electromagnetic propagation in the head, such as in the MEG/EEG inverse problem, where one tries to retrieve locations of neural activity by measuring the electro-magnetic field induced by the electrical activity of the cortex [START_REF] Faugeras | Variational, geometric, and statistical methods for modeling brain anatomy and function[END_REF][START_REF] Wolters | Influence of tissue conductivity anisotropy on EEG/MEG field and return current computation in a realistic head model: A simulation and visualization study using high-resolution finite element modeling[END_REF]; or for computing the effects of electromagnetic fields induced by cellular phones on the brain. Moreover, the mapping of measurements from positron emission tomography (PET-scan) or functional magnetic resonance imaging (fMRI) on highly detailed cortical models helps a lot in the understanding of neural activity organization. At last, the morphometric study of sulcal and gyral patterns can lead to clinical diagnosis of pathological diseases.

In general, these models include 5 structures nested inside each others like "Russian nested dolls". They include, from outer to inner structures: the skin, the skull, the CerebroSpinal Fluid (CSF), the Grey Matter (GM), and the White Matter (WM) [START_REF] Blinkov | The Human Brain In Figures and Tables[END_REF].

In this work, the segmentation of these anatomical structures is performed on couples of co-registered T1-and T2-weighted MRIs. These two images provide rather complementary informations concerning the localization of anatomical structures, but most of the methods create head models from single T1 MRIs, and neglects the T2 MRIs. This chapter reflects the benefits from using both T1 and T2 images, as it allows the segmentation of finer structures, especially concerning the segmentation of the skull. All the couples of T1 and T2 MRIs used in this chapter comes from the Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany.

The segmentation of the scalp is relatively easy to find as there is a clear separa-tion in the T2 MRI between the fat present in the skin and the air. However, there are some areas where the fat layer is too thin, which result in no visible separations between the air and the skull. Fortunately, these areas are very localized (nasals cavities and acoustic canals), and increasing the regularity of the deformable surface prevents potential leaks.

The segmentation of the skull from a single MRI (T1 or T2) is difficult and challenging as such images only provide incomplete informations. Most of the time, the skull is approximated by inflating a convex envelop of the cortex, which may be acceptable for the MEG inverse problem but not for the EEG one.

There exists some approaches based on mathematical morphology that manage to extract the skull from single MR data [START_REF] Richard | Segmentation of skull and scalp in 3-d human mri using mathematical morphology[END_REF][START_REF] Rifai | Segmentation of the skull using deformable model and taking partial volume effect into account[END_REF]. Unfortunately, due to the limitation of the skull definition in single MRIs, these approaches are quite approximative and only manage to extract correctly the upper part of the skull.

Several approaches performs the skull segmentation by using CT data instead of MR data. While CT data provide a better definition of skull, it does not give any informations about interfaces between soft tissues. However, an accurate and complete head model can be constructed by extracting the skull from CT data, and other interfaces from MR data [START_REF] Studholme | Automated 3d registration of truncated mr and ct images of the head[END_REF]. Unfortunately, having two different kind of acquisitions for the same subject is often inconvenient and uncommon, and expose the subject to nocive X-ray radiations. On the contrary, acquiring both T1 and T2 MRIs can be done in a single MRI session, and provide a rather good information about the skull localization.

Extraction of the inner and outer cortical surfaces from T1 MRI data is a difficult task due to its highly convoluted nature [START_REF] Davatzikos | Using a deformable model to obtain a shape representation of the cortex[END_REF][START_REF] Wells | Adaptative segmentation of MRI data[END_REF][START_REF] Macdonald | Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI[END_REF]. ATLAS based approach do not typically perform well the segmentation of the cortex due to high intersubject versatility. In general, most of the methods are based on successive steps to segment the cortex [START_REF] Faugeras | Variational, geometric, and statistical methods for modeling brain anatomy and function[END_REF][START_REF] Han | CRUISE: Cortical reconstruction using implicit surface evolution[END_REF][START_REF] Shattuck | Brainsuite: An automated cortical surface identification tool[END_REF][START_REF] Mangin | Robust brain segmentation using histogram scale-space analysis and mathematical morphology[END_REF][START_REF] Cointepas | Brain-VISA: Software platform for visualization and analysis of multi-modality brain data[END_REF]]:

• It is common to first remove non-brain tissue from the MR volume in order to facilitate the segmentation. This approach is called skull-stripping and consists in removing the skull and all other extracranial tissue from the brain. Several methods exist to perform this, and most of them are based on morphological operations or atlas-based registrations [START_REF] Souplet | A comparative study of skull stripping methods in relapsing-remitting multiple sclerosis: Emergence of a new automatic segmentation algorithm[END_REF][START_REF] Shattuck | Brainsuite: An automated cortical surface identification tool[END_REF][START_REF] Han | CRUISE: Cortical reconstruction using implicit surface evolution[END_REF][START_REF] Ségonne | A hybrid approach to the skull stripping problem in MRI[END_REF].

• The second step consists in image enhancements and low level segmentations.

Given the skull-stripped brain, this task consist in removing inhomogeneities present in MRIs, and classifying the voxels into three different classes: WM, GM and CSF. This labellization is in general performed using Markov Random Fields methods (MRF) or some clustering algorithms (EM, fuzzy C-means, . . . ) [START_REF] Shattuck | Brainsuite: An automated cortical surface identification tool[END_REF][START_REF] Wells | Adaptative segmentation of MRI data[END_REF][START_REF] Held | Markov random field segmentation of brain MR images[END_REF][START_REF] Kapur | Segmentation of brain tissue from mr images[END_REF][START_REF] Zhang | Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm[END_REF].

• The last step consists in removing topological inconsistencies from the previous labellization of the brain. In general, this step is performed by evolving deformable surfaces [START_REF] Davatzikos | Using a deformable model to obtain a shape representation of the cortex[END_REF][START_REF] Malladi | A real-time algorithm for medical shape recovery[END_REF][START_REF] Baillard | Segmentation of 3d brain structures using level sets and dense registration[END_REF] with topological constraints [START_REF] Han | Topology correction in brain cortex segmentation using a multiscale, graph-based algorithm[END_REF] inside the brain toward the interfaces separating the GM and the CSF. This approach allows to get a sub-voxellic representation of the brain, with the correct toplogy.

An interesting approach has been proposed in [START_REF] Zeng | Segmentation and measurement of the cortex from 3D MR images using coupled surfaces propagation[END_REF] and [START_REF] Macdonald | Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI[END_REF], and later in [START_REF] Goldenberg | Cortex segmentation: a fast variational geometric approach[END_REF] 1 and [START_REF] Gomes | Implicit representations of evolving manifolds in computer vision[END_REF]: the authors used a system of coupled surfaces initialized as two concentric spheres inside the white matter, to extract inner and outer cortical surfaces simultaneously. They based their methods on the assumption that the cortex is a layer of almost constant thickness, and used a traditional boundary based evolution with incorporation of a force that prevents the coupled surfaces to move away, or being too close to each other. This approach has the advantage of being robust to locally diffuse contours due to the partial volume effect. However, this assumption is not valid for all cases as the cortex thickness may vary a lot from one subject to another.

In addition to the complexity of these structures, some problems arise due to the limitations of medical imaging equipment:

• Inhomogeneities are quite common in MRIs due to the non uniformity of the MR field. It results in low frequency alterations, or intensity drifts, that smoothly perturb the intensity of the anatomical structures.

• Most of the MRIs suffer from noise coming from the imaging equipment. This noise follow a Rician distribution [START_REF] Sijbers | Signal and Noise Estimation from Magnetic Resonance Images[END_REF], but can be simplified as a Gaussian distribution. This noise can be quite strong depending on the acquisition parameters.

• The partial volume effect occurs when a single voxel contains a mixture of multiple tissue. It results in locally diffuse boundary between distinct anatomical structures.

Because of these alterations, most of the methods cannot use directly active contours to extract the cortex from MRI and some preprocessing is needed. As discussed in chapter 5, these alterations make classical edge-based and region-based model for active contours inefficients.

In this chapter, we propose a framework for creating a complete head model by using the segmentation method proposed in chapter 5, using bimodal registered T1and T2-weighted MRIs. Our goal is to segment anatomical structures directly from the initial MRIs, without any intermediate step that would require additional parameters.

ANATOMICAL REPRESENTATION OF THE HEAD MODEL

Before presenting our method, let's see a brief overview of the anatomical structures present in the head. All the following anatomical pictures comes from the Grant's Atlas of Anatomy, 12 th edition. The head can be simplified as a set of consecutively nested anatomical structures (Fig. 6.1). From outer to inner structure, it includes: Figure 6.1: The head can be seen as consecutively nested anatomical structures.

• The scalp, which is formed of several thin layer besides the skin (Fig. 6.1). It covers the whole head, apart from the nasal cavities and acoustic canals.

• The skull, constituted of several bones that englobe and protect the brain (Fig. 6.11). The skull covers entirely the brain and includes one big hole corresponding to the spinal cord, and some tiny holes (called foramens) in which nervous fibers and veins pass through.

• The CSF, which constitute the area between the brain and the skull. This area is mainly composed of liquid, but includes various type of anatomical tissues (dura, cerebral falx, vascularization, . . . ). The CSF also goes deep inside the brain, in areas that correspond to the brain ventricles.

• The GM, which correspond to the cerebral cortex (Fig. 6.3). The cortex is a highly convoluted thin layer of grey matter, of 1mm to 5mm thickness, with an average of nearly 3 mm [START_REF] Zilles | Cortex[END_REF][START_REF] Blinkov | The Human Brain In Figures and Tables[END_REF]. It is composed of two separate hemispheres, has a spherical topology, and is located beetween the WM on the inside, and the CSF on the outside [START_REF] Blinkov | The Human Brain In Figures and Tables[END_REF].

• The WM, which is the last structure corresponding to the "inside of the brain" (Fig. 6.3). It is mainly composed of nervous fibers linking the neurons of the cortex to each others. It includes the corpus callosum that constitute the only link between the two hemispheres.

CREATION OF HEAD MODEL FROM BI-MODAL ANATOMICAL MRIS

Our method consists in successively segmenting the five anatomical structures present in the model, using bimodal registered T1-weighted and T2-weighted MRIs. This bimodal approach is driven by the fact that a combination of T1-weighted or T2-weighted MRI brings significantly more information about tissues, and especially about the skull interfaces. Indeed, these two images are rather complementary, and while in some places the T2-weighted MRI provide a better localization of inner skull interface, T1-weighted MRI provide a globally better localization of the outer skull interface. Moreover, it does not need much additional efforts to obtain several MRIs of the same subject, with different weightings. Some centers are even starting to systematically acquire both these images.

MRIs Acquisition & Registration

All the couples of T1 and T2 MRIs used in this chapter have been acquired from a SIEMENS Trio 3T imaging device, from the Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany.

The T1-weighted MRIs were acquired using a 3D MP-RAGE [START_REF] James | Three-dimensional magnetizationprepared rapid gradient-echo imaging (3d mp rage)[END_REF] (magnetization-prepared rapid gradient echo) sequence with selective water excita- The T2-weighted MRIs were acquired using a 3D Tourbo Spin Echo sequence, using 8 channel array head coil. The registration between the T1-and T2-weighted MRIs is done as follows: first, the T1-anatomy is aligned with the stereotactical coordinate system proposed by Talairach and Tournoux [START_REF] Talairach | Co-planar Stereotaxic Atlas of the Human Brain[END_REF]. The T2 image is then co-registered with the T1-anatomy using rigid-body transformations [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF], implemented in FSL 2 . Despite some artifacts, the obtained results are globally satisfactory and the two registered images can indeed be used as a single multimodal image.

A couple of registered T1-and T2-weighted MRIs can be seen in Fig. 6.4 through their coronal, frontal and sagittal views. One of the main difference between those two images is that the cerebrospinal fluid appears very dark in the T1 MRI, while being very bright in the T2. It is thus easier to extract the inner interface of the skull and the outer interface of the cortex from the T2-weighted MRI, as there is a clear separation between the skull, the CSF and the brain. On the other hand, the marrow response in the T2 image is also very strong, and this leads to high contrast between the marrow and the bone. As a consequence, the evolving surface will be stuck in these areas, and the skull won't be segmented correctly. However the T1 MRI provides a marrow response close to the bone response, and there are in some parts a better definition of the interface between the skin and the skull. Moreover, the contrast between the GM and the WM tissues is much more pronounced than in the T2 MRI. Therefore the T1 MRI is a better candidate to segment correctly the outer interface of the skull and the inner interface of the cortex. The T1 and T2 MRIs are thus complementary, and a combination of both bring significantly more information about localization of anatomical structures.

Scalp Segmentation

The segmentation of the scalp is done by initializing the deformable surface as a cube around the head, and making it evolve using a simple boundary-based model [START_REF] Caselles | Geodesic active contours[END_REF] on the T2-weighted MRI. Indeed, the T2 MRI provide a better definition of the fat present in the skin than the T1 (Fig. 6.4), and thus allows the extraction of a more detailed scalp surface. The cube will slowly shrink on itself, and finally converge toward the skin (Fig. 6.5, left). Note that the regularity constraint of the model is rather large, in order to avoid the surface from being stuck in noisy areas of the background. The scalp will act as a "mask " in order to constraint the skull segmentation inside it. Indeed, there are some areas where it is impossible to find a clear separation between the skull and the air, like in the acoustic canals or in the nasal cavities. Consequently, the scalp has to include deep details in order to accurately prevent the skull from leaking in the air parts. Once the surface has converged toward the skin, some additional iterations3 are made while reducing the regularity constraint and using the local-based model of chapter 5 with a small fixed neighborhood of 2. The surface will then enter in deep parts of the scalp, like the acoustics canals and the nasal cavities (Fig. 6.5, right).

CSF Segmentation

The CSF is segmented by initializing a sphere inside the brain, and evolving it using the local-based model of chapter 5 with scale selection (2 scales: 4 and 16).

To perform the CSF segmentation, the front has to evolve from the WM through the cortex, to finally fit the interface separating the CSF from the skull. However, the cortex is quite visible in both T1 and T2 MRIs, and without any "tweaking" the surface will likely converge toward the WM/GM or the GM/CSF interfaces. Thus, an a-priori parameter need to be fixed in order to define the desired behavior of the front.

The segmentation model of chapter 5 aims at separating the local distributions of the regions inside and outside the front. Since the model is locally-based, there is not one single distribution for each region, hence the a-priori parameter cannot simply be the type of regions desired (the main advantage of the local approach is its robustness to inhomogeneities). In the bimodal case, these distributions are local bivariate normal distributions. In order to make the front evolve through the cortex, the local distributions of the GM, WM and CSF have to be as similar as possible to be considered as a single anatomical structure. A simple way to ensure this similarity is For example, Fig. 6.6-left shows the histogram of the skull, CSF, GM and WM in a bias corrected bimodal MRI. In order to extract the interface separating the skull from the CSF, the front need to evolve from the WM to the GM to the CSF, without stopping its evolution. By tuning the "orientation" of the distributions (linked to correlation), the GM, WM and CSF can be considered as one single anatomical structure different from the skull. In our case, we computed this correlation from a set of previously segmented head models (good enough to be considered as "ground truth") by estimating the correlations of the region including the CSF, GM and WM. These correlations do not vary a lot from one image to another (an average of -0.7 ± 0.1), and can indeed be used as an a priori parameter to segment the CSF. The front will finally evolve from the WM through the cortex, and will stop on the interface separating the skull from the CSF (Fig. 6.6, right)

The CSF segmentation plays the role of a "skull-stripping " that separates the brain from other extracranial structures. There are various method to perform skullstripping in the literature, and some of them are actually rather complex [START_REF] Souplet | A comparative study of skull stripping methods in relapsing-remitting multiple sclerosis: Emergence of a new automatic segmentation algorithm[END_REF]. Our method manages to segment the CSF (corresponding to the inside of the skull) with a lot of details, especially the veins in the dura and the sagittal sinus (the vein that goes along the separation of the two cortex hemisphere).

WM Segmentation

The WM is then extracted by initializing the front using the previously segmented CSF, and evolving it thanks to the local-based model of chapter 5 with a fixed neighborhood of 8. The WM is a structure composed of very sharp gyri (ridges of the brain), and very smooth sulci (fissures of the brain). In order to accurately segment the sharp gyri of the WM, a common way is to initialize the deformable surface outside the brain (eg. skull-stripping), and makes it deflate in order to fit the ridges that might be difficult to reach from the inside (Fig. 6.7).

The interface between the GM and the WM is particularly visible in the T1 MRI, but the front has to evolve from the CSF through the cortex once again, and thus might be stuck on the interface separating the CSF from the GM. As for the CSF segmentation, we cannot simply fix the type of tissue desired due to the locality of the segmentation model. The WM segmentation is performed in two steps:

• In order to emphasize the WM from the GM and CSF, the segmentation is first performed on a scalar image I 1 corresponding to the combination of the T1 and T2 MRIs (Fig. 6.4, right), such that:

I 1 (x) = max(T 1(x) -T 2(x), 0) ∀x ∈ Ω (6.1)
where I is the combined image, T 1 and T 2 are the two MRIs, and Ω the image domain. By doing this simple difference, the intensity of the CSF and the GM will be close to zero, while the white matter will appear rather bright (Fig. 6.7, left). The front will thus evolve through the cortex, and get closer to the interface separating the WM from the GM.

• The segmentation is then performed on the T1 MRI alone, with a small fixed neighborhood of 2 in order to slightly refine the front position.

GM Segmentation

The GM is then segmented by initializing the front using the previously segmented WM, and evolving it thanks to the local-based model of chapter 5 with a fixed neighborhood of 4.

In opposition to the WM, the GM is a structure composed of very smooth gyri, and very sharp sulci. In order to accurately segment the sharp sulci of the GM, a common way is to initialize the deformable surface inside the brain, and make it inflate in order to fit the fissures that might otherwise be difficult to reach from the outside.

The interface between the GM and the CSF is particularly visible in the T2 MRI, and the front has to evolve from the WM through the cortex once again. Similarly to the WM segmentation, the GM segmentation is performed in two steps:

• In order to emphasize the CSF from the GM and the WM, the segmentation is first performed on a scalar image I 2 corresponding to the combination of the T2 and T1 MRIs (Fig. 6.4, right), such that:

I 2 (x) = max(T 2(x) -T 1(x), 0) ∀x ∈ Ω (6.2)
where I is the combined image, T 1 and T 2 are the two MRIs, and Ω the image domain. Similarly to the WM segmentation, this combination emphasize the CSF from the other tissues, and the front will evolve toward the interface separating the CSF from the GM.

• The segmentation is then performed on the T2 MRI alone, with a small fixed neighborhood of 2 in order to slightly refine the front position.

Topological Constraint

However, one problem concerning the shape of the cortex remains. In some areas, the sulci of the GM are so thin that the CSF lying inside is no more visible in the MRIs, due to the partial volume effect. Therefore, these sulci are not segmented, which impacts significantly the global shape of the brain (Fig 6 .9). This is a major drawback concerning functional studies of the brain, since it creates some "proximity" between cortical areas that might be completely wrong. Shape priors based on training sets of brains [START_REF] Rousson | Shape priors for level set representations[END_REF] are not suitable to segment these sulci, due to the high intersubject variability of the cortex.

To overcome this limitation, classical methods incorporate some topological constraints in order to preserve the spherical topology of the cortex. These constraints aims at preventing any fusion of the front with itself, thus preserving the sulci and gyri. Note that this goes against one of the two strengths of the level-set method (the other one being the absence of parameterization).

Several methods exist to preserve topology. Most of them are quite simple, based on very local criteria that monitor the points in the level-set that could change the topology of the front (non-simple points) [START_REF] Han | A topology preserving level set method for geometric deformable models[END_REF][START_REF] Ségonne | A novel level set framework for the segmentation of medical images under topology control[END_REF]. However, these approach just freeze the front at some specific areas, creating undesired "stair effects". Recently, new methods have emerged to preserve topology, based on variational formulation of higher order active contours [START_REF] Sundaramoorthi | More-than-topology-preserving flows for active contours and polygons[END_REF][START_REF] Rocha | 3d topology preserving flows for viewpoint-based cortical unfolding[END_REF][START_REF] Flynn | Self-repelling snakes for topology-preserving segmentation models[END_REF]. The main idea is, for each point of the front, to integrate over the whole contour in order to look for the points that might get too close. These approaches allow the contour to change its geometry progressively, and no more freezing it "at the last moment". They corrects the "stair case" effect, and provide results of very good accuracy. However, these methods are computationally expensive, since an integral over the contour is required for each point of the contour.

We propose a novel and efficient method for preserving the topology in level-sets. It is based on detecting skeletons of the regions inside and outside the front and "repelling" the front if it gets too close to these skeletons. In level-sets function, skeletons correspond to "shocks" or singularities of the distance function, that can be detected very easily by looking at the gradient of the level-set (however, more accurate methods exists to detect them [START_REF] Shi | Hamilton-jacobi skeleton on cortical surfaces[END_REF]). If a shock is detected, a force is created that locally repels the front, and thus preserves its topology.

The evolution is based on the following equations:

S ǫ (x) = sign(Φ(x)) f ǫ (∇Φ(x)) , ∂Φ ∂t (x) = αS ǫ (x) |∇Φ(x)| , with f ǫ (x) =    1 x < ǫ 0 x > ǫ
where S ǫ correspond to a "skeleton map" that is 1 on the positive skeleton, -1 on the negative skeleton and 0 elsewhere, and α, ǫ are two non-negative scalars. The skeleton corresponds to areas where the gradient norm of Φ is inferior to a certain threshold ǫ. This skeleton map just plays the role of a "conditional balloon force".

This topological constraint is rather simple, and corrects the undesired "stair effect" of the previous methods. It worked quite well for ǫ = 0.5, and the weight α has to be chosen rather large, in order to prevail from other evolution terms. The implementation of this constraint is straightforward and efficient, since it does not need to compute the skeletons in the whole image: for a given point x in the narrow band, the only computation needed to detect if it lies in the skeleton is |∇Φ(x)| (with Neumann boundary conditions). If so, a repelling force appear, preventing the front to evolve towards the skeleton direction.

However, this method suffers from some inconvenients:

• The main inconvenient is that some restrictions on the evolution speed are required: the absolute value of the speed can't be greater than one (one pixel per time step) in order to efficiently preserve topology. A way to overcome this limitation is to incorporate some CFL conditions on the different evolutions terms.

• The two self repelling parts of the fronts are separated by a thin space of one pixel wide corresponding to the detected skeleton (Fig. 6.10). This may be an advantage as well if some methods require this type of constraint (in order to mesh some particular domains, etc). A way to overcome this limitation is to detect skeletons and evolve the front in an image of finer scale.

As can be seen in Fig 6 .9, this constraint works quite well for the segmentation of the GM, as the sulci are clearly visible in the topology-constrained segmentation (However, these sulci are "thicker" than in reality due to the 1-pixel vacuum evoked). In Fig. 6.10 can be seen a 3D segmentation of the cortex under this topology constraint. The exterior skeleton (yellow structure) prevents the fusion of the front in sulci areas, and thus preserve the topology of the front. This skeleton can be seen as a "frontier" that cannot be crossed. Moreover, it corresponds to cortical folds, that are of great interest for morphometric studies of the brain.

As can be seen in Fig. 6.10-Bottom, this constraint is quite easy to implement, as a simple gradient norm is required for the pixels inside the narrow band.

Skull Segmentation

The skull is segmented by initializing the front as a layer between the skull-stripping and the scalp (by inflating the CSF), and evolving it using the local-based model of chapter 5 with a small fixed neighborhood of 2.

The skull is a structure rather difficult to segment since anatomical MRIs (T1 and T2) only provides incomplete informations about its position. However, as stated in Section 6.3.1, the T1-and T2-weighted MRIs are complementary concerning the localization of the skull, and a combination of both images bring more information.

The first approach we used was to combine both T1 and T2 MRI into one scalar image I 3 , such that:

I 3 (x) = T 1(x) . T 2(x) ∀x ∈ Ω
This combination indeed emphasized the skull from other tissues (Fig. 6.4), and its segmentation provided good initial results. Unfortunately, this worked on relatively few subjects, and was not good enough to be included in an automatic framework for head model creation.

Then, we followed the same approach used for the CSF segmentation, by performing the segmentation on the bimodal T1 T2 image, and fixing the correlation to emphasize the skull from the other tissues. This combination provided better results (Fig. 6.11), and have been positively reproduced on different subject (Fig. 6.16). However, it is still impossible to find a clear separation between the skull and the air in regions where the skin is too thin, like the acoustic canals and nasal cavities. Therefore, we used the previously segmented scalp as a mask to locally stop the evolution of the deformable surface. In Fig. 6.12, the scalp and the skull are shown as well as their projection in the T2 MRI. It can be seen that the skull does not leak in the acoustics canals or nasal cavities, thanks to the highly detailed scalp segmentation. Figure 6.12: The skull does not leak in the background thanks to the scalp.

There are unfortunately some segmented structures that do not belong to the skull, and one of the main drawback is that some veins are also segmented. Indeed they appears very dark in the T2, and a bit brighter in the T1 (Fig. 6.4) and the combination of both image make their intensities close to the one of the skull. As a consequence, the vein close to the inner skull interface (superior sagittal sinus) as well as some carotids are falsely segmented. Note that, for example in Fig. 6.11,right,second column, the zygomatic arch (Fig. 6.11) is a structure that is difficult to segment in single MRI, as it appears to be part of the darker structure on its top in the T1 MRI, and on its bottom in the T2. The combination of both help extracting it, without leaking in other neighboring tissues.

VALIDATION

We validated our model according to the GM and WM segmentations. We did not manage to validate the skull segmentations since our results include several other anatomical tissues (bone marrow, veins, air in bucal cavity) that strongly perturb the validation process.

In order to validate results of medical segmentations, a ground truth corresponding to a labellization is required. Most of the ground truth are done by experts, that manually labellize voxels in MRIs to particular anatomical structures. An interesting and accurate validation method can be performed using the visible human database 4 . It consist in several acquisitions of a complete human (CT, T1 MRI, T2 MRI, etc), as well as true colored photos of a sliced version of the human. These photos are so detailed that differentiation between anatomical structures becomes trivial, and help defining with precision a ground truth.

In our case, we validated our method using the BrainWeb simulations 5 . It consists in simulating MRI acquisitions, given specific parameters and percentage of synthetic alteration (noise and bias). These simulations starts from a brain phantom, which plays the role of ground truth for the MRIs simulations. The validation of our model can be found in Fig. 6.13. We tested our model on 3 simulated images : one free from alterations (0% noise & 0% bias) ; one slightly altered (5% noise & 20% bias); and one strongly altered (9% noise & 40% bias). We found an average of almost 0,90 in Dice coefficient, which is quite good given the strong alterations that occur in the last image. Possible errors comes from the topology constraint, that require a vacuum voxel between each segmented sulci.

Images

Results

Ground Truth

.14: Dice coefficient: a ratio between the TP, and the FP and FN is computed. Our framework have been tested on 8 different subjects, and for each of them, the 5 anatomical structures segmented (Fig. 6.16). All these segmentations have been done automatically using the same set of parameters. The only interaction corresponds to the initialization of the front inside the brain (for the CSF segmentation): a sphere can be initialized almost anywhere in the brain, and our model will be robust enough to make it evolve correctly. Moreover, since both T1 and T2 MRI are aligned with the Talairach stereotactical coordinate system, a general initialization can be found that will be valid for all the registered subjects. The results are particularly good, especially for the skull that is quite uncommon in standard methods.

RESULTS

All these segmentation benefit from the GPU implementation of the local-statistic based model from chapter 5, allowing very fast computations (less than 1 hour to segment the 5 structures).

In Fig. 6.15 can be seen slices of the 5 segmentations in both T1 and T2 MRIs. The 5 structures are indeed nested inside each other, and the GM satisfy its spherical topology. In Fig. 6.17 can be seen the 5 interfaces simultaneously, on 2 different subjects. Finally, Fig. 6.18 shows zoom parts of the segmentation on the T1 and T2 MRIs. As can be seen, the topology constraint works quite well, as it prevents self intersections of the front.

CONCLUSION

This work has shown that rather good head models can be obtained using combinations of T1 and T2 MRIs. These images provide complementary informations and can conveniently be acquired in a single MR session.

The segmentation method is based on local statistics of images, which allows for thin structures extractions and naturally copes with the intensity drifts that arise in MR images.

However, instead of a being a global multi-partitioning segmentation model (like EM classifiers for example), this model remains a local bi-partitioning segmentation model, and is not quite suited to segment several regions simultaneously. Hence, our framework is based on consecutive steps to make the front evolve correctly. Our segmentation model works very well in the majority of cases, but is much more powerfull in "refining process", in order to adjust locally the contour from a close initialization. Therefore, we think our framework would benefit from using first a more global criteria, and then refining the results with our model on bimodal T1-and T2-weighted MRIs. However, avoiding classical preprocessing steps (skull-stripping, classification) is a strong advantage of our approach. In our framework, 5 structures are segmented consecutively: the skin, the Cere-broSpinal Fluid (CSF), the White Matter (WM), the Grey Matter (GM), and the skull (Fig. 6. [START_REF] Chan | Active contours without edges[END_REF]).

• The skin has been successfully segmented from the T2 MRI, and provides highly detailed results especially for thin structures such as the nose and the hears. Segmentation of these structures from the T1 MRIs is rather approximative since they confound themselves with the background.

• The CSF has been successfully segmented from a combination of T1 and T2

MRIs, and provides highly accurate results of the inside of the skull. Very nice details are extracted, especially vascularizations of the brain. This segmentation plays the role of a skull-stripping.

• The WM has been segmented from combined ''T1-T2" first, then on the T1 alone. Maybe one problem would be that the sharpest gyri are not completely segmented due to the regularity constraint (Fig. 6.18). Otherwise, the results are pretty good.

• The GM has been segmented from combined ''T2-T1" first, then on the T2 alone.

A topological constraint based on skeleton extraction has been proposed, and is efficient and easy to implement. Moreover, it allows the simultaneous extraction of sulcal folds, which are quite important for morphological studies of the brain.

• Finally, the skull has been segmented from a combination of T1 and T2 MRIs.

The skull is the most difficult structure to extract and maybe the more interesting for electromagnetic simulations (EEG, MEG, etc). Few models manage to extract it correctly from MRIs, due to rather incomplete informations. Our method produce very encouraging results, since the skulls are almost perfectly segmented. However, the proposed method still mixes the skull with other compartments in some places (eg air, veins and ligaments). In many cases (air and ligaments around the mouth, veins in the bottom part of the head), these defaults are not very important for the head modelling as they occur far from the area of interest and have thus little influence. It might still be interesting to remove them for the simplicity of the model or for other applications. However, the superior sagittal sinus (the vein that goes along the separation of the two cortex hemisphere) is currently segmented as being part of the skull6 .

Our framework has been successfully applied to 8 different real acquisitions and shows quite nice details. The timings are quite reasonable thanks to an efficient GPU implementation.

Future works include incorporation of a more global multi-region segmentation for the brain, followed by a refinement process using our local segmentation model. This approach would take in consideration the prior knowledge that the brain is composed of only three separate tissues having a different MR response: WM, GM, CSF; which is somehow neglected in our approach.

An interesting addition would be to incorporate a more elaborate shape prior term, that would characterize the particular shape of the cortex: it would ensure its layered nature (but with a stronger flexibility than classical coupled surface approaches [START_REF] Zeng | Segmentation and measurement of the cortex from 3D MR images using coupled surfaces propagation[END_REF][START_REF] Macdonald | Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI[END_REF][START_REF] Goldenberg | Cortex segmentation: a fast variational geometric approach[END_REF]), as well as its spherical topology and the fact that cortical folds are espaced with a rather constant distance (Fig. 6.10). The cortical folds would be detected, and evolved in order to take into account this spacing constraint. Their evolutions would influence evolution of the GM segmentation, and thus would produce more accurate results. We think this can be done with a higher order active contour model [START_REF] Rochery | Higher order active contours[END_REF][START_REF] Horvath | A higher-order active contour model of a 'gas of circles' and its application to tree crown extraction[END_REF], that is quite powerful to elaborate more complex shape constraints.

Finally, the segmentation of the skull provides quite encouraging results, but some artefacts still remain. Correcting these artefacts could be done using some post-processing method, eg an evolution starting with the segmentations obtained with the proposed method but based on other T1 and T2 image combinations. Some preliminary tests show that this seems feasible.

Aknowledgements:

We would like to thanks Alfred Anwander from Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany for providing all the couples of registered T1 and T2 MRI. head model creation. This consists in segmenting five structures of interest from anatomical MRIs: the White Matter (WM), the Grey Matter (GM), the CerebroSpinal Fluid (CSF), the skull and the scalp.

Creation of head models from MRIs has been widely studied in the last decade, and most methods are based on separate consecutive steps. Moreover, standard models cannot segment easily the skull. We propose a semi-automatic framework that segments each of the anatomical structures successively from bi-modal registered T1-and T2-weighted MRIs. This framework is robust to classical alterations in MRIs thanks to the previous segmentation model, and manage to extract the 5 structures of interest with a very good accuracy.

In summary, the main contributions of this thesis are:

• Implementing the level set method using quadrilateral finite elements.

• Proposing a novel and efficient local-statistic based segmentation model robust to MRI alterations. Pointing out two drawbacks of the local model, and proposing two solutions.

• Creating anatomical head model using the previous local-based model. Segmenting the skull from bi-modal MRIs. Proposing a simple and efficient topological constraint for the cortex segmentation. Some perspective would be:

• Investigating another, more efficient, way to implement level-sets using finite elements (for example, avoiding to include a gradient descent inside another gradient descent).

• Segmenting other kind of images with the local-statistic based model, such as textures data, or high dimensional data.

• Incorporating some global criteria for the segmentation of the brain, and adding some prior shape knowledge for the cortex segmentation, using higher order active contours.

Neurosciences have dramatically improved in the last few years, mainly due to the apparition of more and more detailed images, and new non-invasive functional techniques allowing to distinct with precision the various neural activities in the cortex. Moreover, a very recent imaging technique called Diffusion Tensor Imaging has been introduced. It allows to discern with precision the different nervous fibers composing the white matter. Coupling these techniques opens a wide range of new research concerning the understanding of the brain, since it allows to visualize the different actived cortical zones, and to "backtrack" this activity by following their nervous fibers through the white matter. Head model creation plays a very important role in these research, as it corresponds to the "first step": localizing the differents structures of interest, especially the cerebral cortex. Accuracy of the cortex is crucial to determine precisely in which cortical zone belong a specified neural activity.

Although the brain is without a doubt the most complex anatomical structure of the human body, these neuroscience techniques are very encouraging, and are certainly at the premises of a complete understanding of the brain. CHAPTER 8

CONCLUSION (VERSION FRANC ¸AISE)

Dans cette thèse, nous avons proposé quelques contributions à la segmentation d'images, dans le but de créer des modèles anatomiques de la tête à partir d'Images à Résonnace Magnétique (IRMs).

Nous avons d'abord brièvement présenté un aperc ¸u des différentes méthodes en segmentation d'images. Bien que les modèles discrets (Graph-cut, etc), sont aujourd'hui largement utilisés en segmentation, les modèles basés sur des évolutions d'hypersurface sont encore le meilleur choix pour incorporer certaines connaissances a-priori sur la forme de la segmentation souhaitée. Dans notre cas, la préservation sphérique de la topologie du cortex est cruciale, et ne peut être mise en place facilement en utilisant des modèle discrets.

Nous avons ensuite insisté sur le fait que les modèles basés sur des évolutions d'hypersurfaces sont classés en deux catégories (basés contours et basés régions), et aucune de ces catégories n'est robuste à l'ensemble des altérations présentes dans les IRMs (bruit, inhomogénéités, contours diffus).

Nous avons ensuite proposé une nouvelle implémentation pour la méthode des ensemble de niveaux, en utilisant des éléments finis quadrilatéraux. Cette représentation vise à améliorer la robustesse et la précision du modèle tout en offrant une bonne représentation sous-voxellique du niveau zéro. Toutefois, sa mise en oeuvre est très difficile, et augmente largement la complexité du modèle. Nous n'avons pas utilisé cette implémentation dans le reste de cette thèse en raison de son importante complexité. À la place, nous avons opté pour un modele classique à base de differences finies, avec quelques routines de parallélization.

Nous nous sommes ensuite focalisés sur des modeles de segmentations d'image, et avons proposé un nouveau modèle basé sur des statistiques locales, robuste à toute les altérations dont souffrent les IRMs. Sa formulation unifie les modèles basés contours et les modèles basés régions, et son implémentation est particulièrement rapide à base de filtres récursifs et d'implémentation GPU. Nous avons ensuite mis en évidence deux inconvénients majeurs, et avons proposés deux solutions. Ce modèle a fournis de bons résultats, en accomplissant simultanément la segmentation et la restauration d'images (les résultats produit sont des images lisses par morceaux).

Enfin, nous avons appliqué le précédent modèle de segmentation dans le but de créer des modèles anatomiques de la tête à partir d'IRMs. Cela consiste à segmenter cinq structures anatomiques d'intérêt dans les IRMs: la matière blanche (MB), la matière grise (MG), le liquide céphalo-rachidien (LCR), le cr âne et le cuir chevelu.

Les créations automatiques de modèles anatomiques de tête à partir d'IRMs ont largement été étudiés au cours des dernières années, et la plupart sont fondées sur différentes étapes consécutives. En outre, les modèles classiques ne peuvent facilement détecter le cr âne à cause d'un manque de donnée concernant sa localisation. Nous avons proposé un cadre semi-automatique qui segmente chaque structure anatomique de manière successive à partir d'IRMs bi-modales T1 et T2 recalées. Ce cadre est robuste aux altérations classiques des IRMs, et permet d'extraire avec une très bonne précision les 5 structures d'intérêt.

En résumé, les principales contributions de cette thèse sont:

• Implémentation de la méthode des ensemble de niveaux à l'aide d'éléments finis quadrilatéraux.

• Introduction d'un nouveaux modèle de segmentation basé sur des statistiques locales. Mise en évidence de 2 inconvénients majeurs dans le modèle, et proposition de 2 solutions.

• Création de modèle anatomique de la tête en utilisant le modele local. Segmentation du crane à partir d'images bi-modales. Introduction d'une nouvelle contrainte topologique simple et rapide pour la segmentation du cortex.

Cette thèse ouvre quelque perspectives, qui sont:

• Implémenter de manière plus rapide et efficace la méthode des ensemble de niveaux à l'aide d'élements finis (par exemple, éviter d'inclure une déscente de gradient à l'interieur d'une autre déscente de gradient).

• Segmenter d'autre types d'images à l'aide du modèle de segmentation à base de statistiques locales, comme des images texturées, ou des images à haute dimension.

• Inclure des critères un peu plus globaux pour la segmentation du cerveau.

Utiliser des contraintes de formes plus élaborées pour la segmentation du cortex, en utilisant des contours actifs d'ordre superieurs.

Les neurosciences se sont amélioré de fac ¸on spectaculaire au cours des dernières années, principalement en raison de l'apparition d'images anatomiques de plus en plus détaillées, et de nouvelles techniques non-invasives permettant de discerner avec précision les différentes activités neuronales du cortex. En outre, une technique d'imagerie appelée Imagerie par Tenseur de Diffusion a été très récemment introduite, et permet de discerner avec précision les différentes fibres nerveuses qui composent la matière blanche. Le couplage de ces technique ouvre un large éventail de recherches concernant la compréhension du cerveau, en permettant de visualiser les différentes zones corticales en activité, et de "suivre" cette activité à travers les fibres nerveuses de la matière blanche. Les modeles anatomiques de tête jouent un rôle crucial dans ces recherches, car il permettent de localiser les différentes structures anatomiques d'intérêt, en particulier le cortex cérébral. La précision du cortex est primordiale pour déterminer les zones corticale auxquelles appartiennent les activité neuronales détectées.

Bien que le cerveau soit sans doute la plus complexe des structures anatomiques du corps humain, ces techniques de neurosciences sont très encourageantes et sont certainement les prémices vers une compréhension globale du cerveau.

For example, if the element e is a 2D triangle defined by vertices V 1 = (1, 0), V 2 = (1, 1), V 3 = (0, 1), we will have:

A.2 Q1 BASIS FUNCTIONS

We used the method of appendix A.1 to compute the basis functions of 2D (square) and 3D (cube) regular Q 1 elements of size 1. Note that the gradient of these basis functions inside the elements is not constant, in opposition to P 1 elements.

The formulation of these basis function will help computing the terms of matrice A 1 and A 2 depicted in Section 4.3.6.

A.2.1 2D basis functions

The 2D basis functions ϕ associated to the vertex V = (a, b) is defined of the 4 elements adjacent to V . The restriction of ϕ inside these 4 elements is shown in

A.2.2 3D basis functions

The 3D basis functions ϕ associated to the vertex V = (a, b, c) is defined of the 8 elements adjacent to V . The restriction of ϕ inside these 8 elements is shown in

A.3 CALCULUS OF DISTANCE CONSTRAINT SYS- TEM

The distance constraint equation:

(∇u + ∇v) 2 -1 = 0 (A.1)
gives us the following system to solve, after some approximations: The next step is to use the basis functions of appendix A.2 to compute these specific integral. These computations are different for the 2D case and 3D case.

2 u T Q i A i 1 v + u T Q i u - Ω ϕ i -b i 1 = 0 i ∈ {1, 2

A.3.1 2D case

In the case of a 2D mesh with only one regular Q 1 element (therefore with 4 vertices), the matrices Q i can be computed using integrals previous 2D Q 1 basis functions. The matrix Q 0 associated to the vertex 0 is defined as follows:

Q 0 ab = Ω ∇ϕ a , ∇ϕ b ϕ 0 0 3 2 1 Q 0 = 1 12       3 -1 -1 -1 -1 2 0 -1 -1 0 1 0 -1 -1 0 2      
By developping the term 2 u T Q 0 v using the matrix above, it appears that we can reformulate it in terms of difference between all possible pairs of vertex, which is directly linked to edges of the element. This formulation is a lot easier to implement, as the system can then be retrieved by looking at all possibles edges of a given mesh. In the following, we will use the notation u ij = (u iu j ) and v ij = (v iv j ).

For the vertex 0, this lead to the following notation: that can be computed by looking at the vertical, horizontal and diagonal edges.

This help us find the term A 0 1 v of system A.2. The first term of b 0 1 is found in the same way and the second term ( Ω ϕ i ) corresponds to the number of elements adjacent to vertex 0. These notations are valid for vertex 0, but calculus confirms that this can be generalized for all the other vertices of the mesh, by looking at the same "adjacency" of edges. Moreover, calculus shows that if vertices 0 is adjacent to more elements, each element simply add its contribution in the same way to the equation induced by vertex 0 in the system A.2.

A.3.2 3D case

In the case of a 3D mesh with only one regular Q 1 element (therefore with 8 vertices), the matrices Q i can be computed using integrals previous 3D Q 1 basis functions. The matrix Q 0 associated to the vertex 0 is defined as follows:

Q 0 ab = Ω ∇ϕ a ,
-1 -5 1 -3 -1 1 7 -1 1 -5 -3 1 -1 1 -1 7 1 -3 -1 -1 1 -1 1 1 3                
By developping the term 2 u T Q 0 v using the matrix above, the sames reformulations as the 2D case can be made, by expressing it in terms of difference between all possible pairs of vertex. In the following, we will use the same notation u ij = (u iu j ) and v ij = (v iv j ).

However, the 3D case is more complicated since several type of "edges" contribute differently to the reformulation of the system of equation. In a single element, there is a total of 28 "edges": 12 vertical & horizontal edges ; 12 diagonal edges inside the 6 faces ; 4 "internal edges".

For a given vertex, there are 4 types of contributions to its associated equation, depending on a specific "adjacency" of these edges: the 9 horizontal & vertical edges not adjacent to it; the 6 diagonal edges adjacent to it ; the 6 diagonal edges not adjacent to it ; and the other edges (hotizontal, vertical, and internal edges). For the vertex 0, This lead to the following notation:

2 u T Q 0 v = 1 144 -Q 0 a + 5Q 0 b + Q 0 c + 3Q 0 d ,

with

• Horizontal & vertical edges not adjacent to vertex 0 This help us find the term A 0 1 v of system A.2. The first term of b 0 1 is found in the same way and the second term ( Ω ϕ i ) corresponds to the number of elements adjacent to vertex 0.

As for the 2D case, these notations are valid for vertex 0, but calculus confirms that this can be generalized for all the other vertices of the mesh, by looking at the same "adjacency".

Therefore, all these formulae help us computing the system A.2 very easily, by looking at all elements in the mesh, and for each of them adding its contribution to its 8 vertices they include, by looking at all their possible "edges".

A.5 CONJUGATE GRADIENT METHOD

The conjugate gradient method is an iterative method to approximate the solution of a system of equations by minimizing an energy induced by the system [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. This method is effective for systems of the form:

Ax = b
where x is a unknown vector, b is a known vector and A is a square, symmetric, positive definite matrix. If the matrix A is not symmetric, one can still solve the following system:

A T Ax = A T b
The basic idea of gradient descent is to take an initial solution as a "starting point", and to consecutively improve this solution by updating its values in order to decrease the energy induced by the system. This update is performed differently depending of the method:

• in the steepest gradient descent, this update is performed by going in the opposite direction of the gradient of the energy, with constant steps. By following this direction, the solution will converge toward the minimum of the energy, and thus will solve the initial system.

• the steepest gradient descent with optimal step proceeds in the same way, but the difference is that each step are optimally chosen to converge more rapidly to the solution (the further the solution is from the minimum of the energy, the larger the step will be) (Fig A .3). Figure A.3: Minima of the energy obtained with the steepest gradient descent with optimal steps (Left) and with the conjugate gradient method (Right).

• the conjugate gradient method proceeds in a similar way to the previous one, but the difference is that the direction is no more the gradient of the energy, but a more complicated direction depending on all the directions already taken in the descent. It has the advantage of converging to the optimal solution in at most n step, where n is the number of equations in the system (Fig A .3).

The conjugate gradient method is the most efficient algorithm, as it guarantee that the solution will be reached in at most n steps.

The algorithm below explains how to implement the conjugate gradient method, in order to solve a system of the form Ax = b. with

f s (x, Ω, V) = f G 1 (x, Ω, G 1 , G 2 ) G ′ 1 (x, Ω).V + f G 2 (x, Ω, G 1 , G 2 ) G ′ 2 (x, Ω).V ,
where f G 1 and f G 2 denote the partial derivative of (B.1) with respect to G 1 and G 2 .

They can be expressed as: Putting all the term together yields Ω g ρ (xz)dz dx minimization. It has been shown that without these terms in the derivative, the energy sometimes slightly increase, but nonetheless converge toward the same local minima. Therefore, for the sake of simplicity and efficiency, these terms can be removed from the derivative.

           f G 1 (x, Ω, G 1 , G 2 ) = - 2 G 2 (x, Ω) I(x) - G 1 (x, Ω) G 2 (x, Ω) f G 2 (x, Ω, G 1 , G 2 ) = 2G 1 (x, Ω) G 2 (x, Ω) 2 I(x) - G 1 (x, Ω)
f s (x, Ω, V) = f G 1 G ′ 1 .V + f G 2 G ′ 2 .V
At last, the bi-partitioning energy can be minimized using the following gradientdescent of coupled equation, in level-set notation:

                      

Figure 1 . 1 :

 11 Figure 1.1: Création d'un modèle de la tête à partir d'IRMs.

Figure 2 . 1 :

 21 Figure 2.1: Creation of anatomical head model from MRIs.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Discrete Models: Markov Random Fields and Graph Cuts . . . 3.2.1 Markov Random Field . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Graph Cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Random Walker . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Continuous Models: Snakes and Active Contours . . . . . . . . 3.3.1 The Snake Model . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Front Representation: Explicit models vs Implicit models . . 3.4.1 Explicit Representation . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Implicit Representation: The Level-Set Method . . . . . . . . 3.5 Boundary Based Segmentation Models . . . . . . . . . . . . . . . 3.5.1 Geodesic Active contour . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Gradient Vector Flow . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Region based segmentation models . . . . . . . . . . . . . . . . . 3.6.1 Functional-based models . . . . . . . . . . . . . . . . . . . . . . 3.6.2 Bayesian-based models . . . . . . . . . . . . . . . . . . . . . . . 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Explicit representations of 2D contours.

Figure 3 . 2 :

 32 Figure 3.2: Level-sets: Topological changes are automatically handled.

Figure 3 . 3 :

 33 Figure 3.3: Top row: Evolution of the contour with the geodesic active contour model; Bottom row: Metric of the Riemanian space Eq. (3.4), and representation of the image and the contours in that space (the contour remains stuck into some kind of "saddle").

Figure 3 . 4 :

 34 Figure 3.4: Top row: Initial image; Edge map; Gradient vector flow & zoom; Bottom row: Evolution of the contour with the gradient vector flow model.

Figure 3 . 5 :

 35 Figure 3.5: Top : Evolution of the contour using the Gaussian model; Bottom : Gaussian distribution of the regions at initialization (left) and convergence (right).

Figure 3 . 6 :

 36 Figure 3.6: Top : Evolution of the contour using a non-parametric Parzen model; Bottom : Distribution of the regions at initialization (left) and convergence (right).

Figure 3 . 7 :

 37 Figure 3.7: 3D Segmentation of the inside of the skull (brain + cerebrospinal fluid) on a T2 MRI, using a non-parametric segmentation model.

Figure 3 . 8 : 2 Finite 4 . 3 Finite 4 . 4

 3824344 Figure3.8: Summary of models for image segmentation: there are two separates categories (discrete and continuous) having their own optimization methods (green squares). The continuous models can be seen in a hierarchical way. The blue curves correspond to the possible link between two models.

Figure 4 . 1 :

 41 Figure 4.1: ∂Ω defines the border of a domain Ω, and n the normal vector to Ω.

Figure 4 . 2 :

 42 Figure 4.2: Partition of a 2D domain Ω in triangles.

Figure 4 . 3 :

 43 Figure 4.3: Elements of type P n .

Si

Figure 4 . 5 :

 45 Figure 4.5: Basis functions: piecewise affine functions for triangular P 1 elements (left) and piecewise bilinear functions for quadrangular Q 1 elements (right).

Figure 4 . 6 :

 46 Figure 4.6: Basis functions of Q 1 elements.

  4.7).

Figure 4 . 7 :

 47 Figure 4.7: P 1 finite element representation of the level set function (top), and its mean curvature motion (bottom).

Figure 4 . 8 :

 48 Figure 4.8: top: mean curvature motion of a square ; bottom: anisotropic neighborhood of a vertex, and deviation of the solution from a circle (black discs).

Figure 4 . 10 :

 410 Figure 4.10: Basis functions of a vertex in 3D with Q 1 elements.

(u a + v a ) ∇ϕ a 2 =

 2

Algorithm 1

 1 Level-set evolution using Q 1 finite elements initialize u repeat for all elements e do for all edges i, j ∈ e do Compute contribution to matrices A 1 and A 2 Compute contribution to vectors b 1 and b 2 v ⇐ 0 repeat v ⇐ v + conjugate gradient step of system 4.25 until convergence u ⇐ u + v for all elements e do if e too far from the zero level-set then remove e from the mesh insert new elements using fast-marching extrapolation until convergence

Figure 4 . 11 :

 411 Figure 4.11: 2D mean curvature motion of a square with Q 1 finite-elements.

Figure 4 . 12 :

 412 Figure 4.12: 3D mean curvature motion of a cube with Q 1 finite-elements.

Q FE 1 FDFigure 4 . 13 :

 1413 Figure 4.13: Top row: mean curvature motion of a cube with finite difference implementation; Bottom row: mean curvature motion of a cube with Q 1 finite element implementation.

Figure 4 .

 4 Figure 4.14: Top: 3D mean curvature motion of a cube with Q 1 finiteelements;Bottom: result of the evolution at middle convergence, the distance function is well preserved.

Figure 5 . 1 :

 51 Figure 5.1: MRI difficulties: Intensity of white matter may vary across the image; and contour between white matter and grey matter may be locally diffuse.

Figure 5 . 2 :

 52 Figure 5.2: The Adelson illusion: the human visual system is more sensitive to local variations of the information than to global grey level values.

Figure 5 . 3 :

 53 Figure 5.3: Inside (yellow) and outside (blue) statistics: left computed in the whole region; right computed at the neighborhood of the active contour.

Figure 5 . 4 :

 54 Figure 5.4: Importance of the denominator in the formulation of regions local means.

Figure 5 . 5 :

 55 Figure 5.5: Extraction of the liver from 2D real data. Top: Global mean Chan & Vese model ; Bottom: Local means model computed with σ = 16.

Figure 5 . 6 :

 56 Figure 5.6: Local statistics computation: The regions delineated by the yellow contour (a) are first separated (b), blurred independently (c), and then normalized by a blurred characteristic function of the regions (d). This gives the local mean (e) and the local variance (f) using in this example with a Gaussian kernel of variance 16.

Figure 5 . 7 :

 57 Figure 5.7: Top: Initialization and convergence of the local model ; Bottom: evolution of the contour which lead to the creation of a singularity.

Figure 5 . 8 :

 58 Figure 5.8: (a) Initialization; (b) Convergence without the coherence constraint: the segmentation is "locally correct" but not globally; (c) Convergence with the coherence constraint: intensity change in the same way from the inside to the outside.

Figure 5 . 9 :

 59 Figure 5.9: Simulated annealing: the energy restricted to a point of the contour is locally increased in order to converge toward the desired local minimum.

Figure 5 . 10 :

 510 Figure 5.10: Top row: evolution without the coherence constraint; Middle row: evolution with the coherence constraint; Bottom: Energy of the two evolutions.

Figure 5 .

 5 Figure 5.11: Top: Initialization and convergence of the local model ; Bottom: absence of evolution in locally homogeneous areas which leads to the creation of "layers".

Figure 5 . 13 :

 513 Figure 5.13: Two pass filtering of a 2D image.

Figure 5 . 14 :

 514 Figure 5.14: Top row: Initialization and evolution of the contour with the global Gaussian model ; Middle row: Initialization and evolution of the contour with our model ; Bottom-row: Scales chosen a each point of the contour.

Figure 5 . 15 :

 515 Figure 5.15: Equivalence between the piecewise smooth model and the Mumford Shah functional.

Figure 5 . 16 :

 516 Figure 5.16: Segmentation of diffuse lesions in retina image.

6 . 2 6 . 3

 6263 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Anatomical representation of the head model . . . . . . . . . . Creation of Head Model from Bi-modal Anatomical MRIs . . . 6.3.1 MRIs Acquisition & Registration . . . . . . . . . . . . . . . . . 6.3.2 Scalp Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.3 CSF Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.4 WM Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.5 GM Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.6 Skull Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6 . 2 :

 62 Figure 6.2: Anatomical representation of a human skull.

Figure 6 . 3 :

 63 Figure 6.3: Anatomical representation of a human brain.

  Imaging parameters: TR = 2000 ms, TE = 355 ms, bandwidth = 355 Hz/pixel; image matrix = 256 x 256; FOV = 256 mm x 256 mm; slab thickness = 176 mm; 176 partitions; sagittal orientation; spatial resolution = 1 mm x 1 mm x 1 mm; 1 acquisition.

40 Figure 6 . 4 :

 4064 Figure 6.4: Left: T1-weighted anatomical MRI (top) and T2-weighted anatomical MRI (bottom); Right: Tissues response on T1-and T2-weighted anatomical MRI normalized between 0 and 255.

Figure 6 . 5 :

 65 Figure 6.5: Left : 3D evolution of the scalp segmentation, and corresponding slices in the T2 MRI of the evolution. Right : The nasal cavity and acoustic canals are segmented to prevent the skull segmentation from leaking in these areas.

Figure 6 . 6 :

 66 Figure 6.6: Left: Histogram of 4 anatomical structures in a bias corrected bimodal MRI: the skull (Grey), the CSF (Orange), the GM (Purple) and the WM (white); Right: 3D evolution of the CSF segmentation.

Figure 6 . 7 :

 67 Figure 6.7: Left: projection of the WM segmentation in a slice of the combined I 1 image (Eq.6.1); Right: 3D evolution of the WM segmentation, and corresponding slices in the T1 MRI of the evolution.

Figure 6 . 8 :

 68 Figure 6.8: Left: projection of the GM segmentation in a slice of the combined I 2 image (Eq.6.2); Right: 3D evolution of the GM segmentation, and corresponding slices in the T2 MRI of the evolution.

Figure 6 . 9 :

 69 Figure 6.9: Left: projection of the GM segmentation in a slice of the T1 MRI, without (first) and with (second) the topological constraint; Right: 3D segmentation of the GM, without (first) and with (second) the topological constraint.

Figure 6 .

 6 Figure 6.10: Top: GM segmentation (purple) with external skeleton (yellow); Middle: GM segmentation next to external skeleton; Bottom: Projection of the GM segmentation on : the T1 MRI; its corresponding level-set function (LS); gradient norm of the LS (inside narrow band); detected skeleton.

Figure 6 .

 6 Figure 6.11: Left : 3D view of the segmented skull ; Right : Projection of the skull on the T1 image (top row) and the T2 image (bottom row).

Figure 6 .

 6 Figure 6.15: Projections (coronal, frontal, sagittal) of the WM, GM, CSF, skull, skin interfaces inside the T1 (Left) and the T2 MRIs (Right).

Figure 6 .

 6 Figure 6.16: 3D Results of 6 different subjects. From left to right : skin, skull, CSF, GM and WM.

Figure 6 . 17 :

 617 Figure 6.17: Results of 2 different subjects (Top and Bottom), with their 3D representation (Left), and projections on their respective anatomical MRI (Right).

Figure 6 . 18 :

 618 Figure 6.18: Zoom on projections of the 5 segmentations in the T1 MRI (left) and the T2 MRI (right).

InitializationFigure 6 . 19 :

 619 Figure 6.19: Consecutive steps of our model, using different combinations of T1 and T2 MRIs (see text).

ϕ 1

 1 (x, y) e =

Fig. A. 1 .

 1 These basis functions are characterized by their particular shape, in form of "pagodas" (Fig 4.6).

Figure A. 1 :

 1 Figure A.1: Basis function of 2D regular Q 1 elements associated to vertex V .

Fig. A. 2 . 1 - 1 - 1 - 1 - 1 Figure A. 2 :

 2111112 Figure A.2: Basis function of 3D regular Q 1 elements associated to vertex V .

  , . . . , n} , (A.2) with Q i n × n matrix associated to the vertex i and defined by Q i ab = Ω ∇ϕ a , ∇ϕ b ϕ i .

2 u T Q 0 v = 1 12 u 01 u 02 u 03 u 13 Tv

 113 01 v 02 v 03 v 13

∇ϕ b ϕ 0 6

 6

= u 13 u 15 u 23 u 26 u 37 u 45 u 46 u 57 u 67 T v 13 v 15 v 23 v 26 v 37 v 45 v 46 v 57 v 67 •= u 03 u 05 u 06 u 12 u 14 u 24 T v 03 v 05 v 06 v 12 v 14 v 24 •= u 17 u 27 u 35 u 36 u 47 u 56 T v 17 v 27 v 35 v 36 v 47 v 56 •= u 01 u 02 u 04 u 07 u 16 u 25 u 34 Tv

 67672424565634 Diagonal edges of adjacent faces to vertex 0 Diagonal edges of non adjacent faces to vertex 0 Horizontal & vertical edges adjacent to vertex 0 and internal edges 01 v 02 v 04 v 07 v 16 v 25 v 34

Algorithm 3

 3 Conjugate Gradient Method initialize solution x 0 initialize direction d 0 = r 0 = b -Ax 0 for i = 0 to n do compute step α i = r T i r i d T i Ad i update solution x i+1 = x i + α i d i if solution x i+1 good enough then break • r i+1 = r iα i Ad i • β i+1 = r T i+1 r i+1 r T i r i update direction d i+1 = r i+1 + β i+1 d iIn our case, for the bi-partitioning problem, we have,E ′ (Ω, Ω), V = Ω f s (x, Ω, V) dx + Ω f s x, Ω, V dx -Γf (x, Ω)f x, Ω (V(x).N(x)) da(x)

G 2 (H

 2 1s (x, y, V)dy -Γ H 1 (x, y)(V(y).N(y))da(y) G ′ 2 (x, Ω).V = Ω H 2s (x, y, V)dy -Γ H 2 (x, y)(V(y).N(y))da(y) (B.4)Since H 1 and H 2 do not depend on Ω, we obtain H 1s = 0 and H 2s = 0.

= 2 ( 2 + 2 (

 222 I(x)µ ρ (x, Ω)) Ω g ρ (xy)dy Γ g ρ (xy) (I(y)µ ρ (x, Ω)) (V(y).N(y)) da(y)and finally we obtain, by changing the integration order,E ′ (Ω, Ω), V = Ω f s (x, Ω, V) dx + Ω f s x, Ω, V dx -Γ f (y, Ω)f y, Ω (V(y).N(y)) da(y) = Γq(y, Ω)q(y, Ω) -(I(y)µ ρ (y, Ω)) I(y)µ ρ (y, Ω) V(y).N(y)) da(y) with q(y, Ω) = Ω 2 (I(x)µ ρ (x, Ω)) g ρ (xy) (I(y)µ ρ (x, Ω))

  µ ρ (x, Ω) = (g ρ * H σ (Φ) I) (x) (g ρ * H σ (Φ)) (x) µ ρ (x, Ω) = (g ρ * (1 -H σ (Φ)) I) (x) (g ρ * (1 -H σ (Φ))) (x) ∂Φ ∂t (x) = (I(x)µ ρ (x, Ω)) 2 -(I(x)µ ρ (x, Ω)) 2 |∇Φ| .(B.7)

  Dice coefficient of GM and WM on 3 incrementally altered MRIs.In order to estimate a "ratio of accuracy", a common choice is to use the Dice coefficients related to the Jaccard index, defined as follows:

		GM Dice coefficient WM Dice coefficient
	0% noise & 0% bias	0.94	0.91
	5% noise & 20% bias	0.89	0.91
	9% noise & 40% bias	0.85	0.74
	Figure 6.13: Dice =	2 TP (FP + TP) + (TP + FN)	,
	where TP, FP and FN correspond respectively to the True Positive (result∈ground
	truth), False Positive (result / ∈ground truth), and False Negative (ground truth / ∈result) results (Fig. 6.14). A value of 0 indicates no overlap and a value of 1 indicates perfect agreement. Higher numbers indicate better agreement.

In the next sections, we will use the notion or "front" or "interface".

Courant-Friedrichs-Lewy conditions: conditions linking the time step to the spatial step

Path of minimal length in a Riemanian space.

The length of the front is expressed as the integral of the level-set Heaviside gradient

Due to the large amount of "smoothness" on the region histogram

under certain regularity conditions.

For example, if the order of derivation of the PDE is two (e.g. Laplacian), the polynomial order of the discretization basis has to be at least two (e.g. solution piecewise quadratic).

Using any linear system resolution scheme, such as Gauss, Choleski, or iterative methods like gradient descent.

We can use first order polynomial elements (P 1 , Q 1 ) on first order differential equations.

Either a contour strength or global distributions of pixel values inside regions.

Under some simplifications, see Appendix B.

See section 3.6.2 for more details.

One may add an auxiliary balloon force to make the contour evolve in homogeneous areas

The image gradient is globally oriented either toward the inside or the outside of the object, for all the points lying on its boundary.

www.openmp.org

As pointed out in[START_REF] Pons | Contributions méthodologiques et appliquées à la méthode des modèles déformables[END_REF] the authors introduced a variational formulation but made an abuse in the derivation of the energy with respect to the level-sets.

http://www.fmrib.ox.ac.uk/fsl

Only a few iterations are required, otherwise the surface will grow into deep structures of the skull.

http://www.nlm.nih.gov/research/visible/visible human.html

http://www.bic.mni.mcgill.ca/brainweb/

It is interesting to notice that this vein is hardly distinguishable in T1 MRIs, so that most segmentations based on such images wrongly aggregate this head area with the grey matter.

mean model. As the Chan & Vese model is not robust to bias, the liver is not correctly extracted, and the top "leak" in the part of the image where intensity is close to the global mean inside the front. It finally results in separating the bright area from the dark ones, with no preference to sharp edges. Note that non-parametric Bayesian models produce the same result, as they rely on a global region model as well.

On the contrary, our model behaves quite well as the liver is correctly delimited at the local level. The model is more robust to low frequency alterations, and is well suited for this kind of medical images.

Bayesian Formulation and Local Region Statistics

An alternative approach is to model the problem through a Bayesian formulation. It's a different kind of formulation for image segmentation, but this leads to the minimization of an energy similar to Eq. (5.2). In the following, we note K a partition of the image domain Ω, composed of several regions Ω i such that

We consider P (K|I), the a posteriori probability for having the segmentation K, given the observed image I. The optimal segmentation K is then obtained by maximizing this probability according to the associated hypothesis. Thanks to the Bayes rule, we can write: P (K|I) = P (I|K)

where:

• P (I) is the probability of having an image I. This term is constant during the segmentation, and can be neglected.

• P (K) is the probability of having a segmentation K. This probability stand for the desired properties of the region borders such as their smoothness, global shape, topology, etc ...

• P (I|K) is the a-posteriori probability for having the image I, given the segmentation K.

We choose to use the segmentation probability as a regularizer, by minimizing the total length of the border between regions. This corresponds to maximize the following probability: of Ω i at scale ρ j for all points in the active contour do Find optimal scale ρ opt using Eq. (5.18) with µ i,j and σ 2 i,j

Compute speed with µ i,opt and σ 2 i,opt

Update contour until convergence small scales are more important, so we decided to choose four scales: 2, 4, 8, 16.

The model is described by the following algorithm:

Finding the correct scale from the set of local statistics is very fast, and the bottleneck remains the creation of the scale space of statistics at each iteration. Although the recursive filter is fast, several blurs are needed at each iteration, depending on the chosen scales.

In Table 5.1, the execution time of 500 iterations of the algorithm is shown on a 2D and a 3D image. It is compared with the single scale model, and with its implementation on a GPU (graphics processing unit), using the NVidia's Cuda Library 7 . The GPU implementation leads to reasonable times for both 2D and 3D image. All the convolutions are parallelized by filtering each row/columns simultaneously. However note that modern GPU have a quite limited on-chip memory, and this prevent the parallelization of all the needed convolutions (a total of 6 and for each scale for the bi-partitioning case).

Images

single scale 4 scales 4 scales GPU 2D (300 2 ) 13.17 Sec 56.36 Sec 7.61 Sec 3D (256 3 ) 51Mn 20Sec 3h 4Mn 25Mn

Table 5.1: Time execution of 500 iterations. 'single scale' and '4 scales' corresponds to the CPU implementation of the method, and '4 scales GPU' corresponds to the multi-threaded algorithm on a GPU (NVidia's Quadro FX570M).

RESULTS

7 www.nvidia.com/cuda CHAPTER 7

CONCLUSION

In this thesis, we have proposed new models for image segmentation in the purpose of creating anatomical head models from Magnetic Resonance Imaging (MRIs).

We first presented an overview of different methods for image segmentation. Although discrete models (Graph-Cuts, etc), are now widely used in image segmentation, models based on hypersurface evolutions are still the best choice for incorporating some knowledge about the shape of the desired segmentation. In our case, preservation of the topology for the cortex segmentation is crucial, and cannot be done easily using a discrete segmentation model.

We then emphasized the fact that segmentation models based on hypersurfaces are classified in two categories (edge-based and region-based), and none of these categories is robust to all of the MRI alterations (noise, inhomogeneities, diffuse contours).

Then, we proposed a novel implementation of the level-set method using quadrilateral finite elements. This representation aims at improving the robustness and the accuracy of the model while providing a good sub-voxelic representation of the zero level-set. However, its implementation is quite difficult, and largely increases the complexity of the model. We did not use this implementation in the rest of the thesis due to this high complexity, and used instead a classical finite-difference based implementation, with some parallelization routines.

A new segmentation model based on local statistics that aimed in being robust to common MRIs alterations has been proposed. Its formulation unifies classical edge-based and region-based models, and its implementation is quite effective using fast recursive filters and modern GPUs. We then pointed out two major drawback of the model, and proposed two solutions. This model provides good results, and manages to perform both segmentation and restoration, by providing piecewise smooth partition of images.

Finally, we applied the previous segmentation model in the case of anatomical

Appendix DETAILS ABOUT FINITE-ELEMENT

LEVEL-SET IMPLEMENTATION

This appendix describes some details about the implementation of the level-set method using Q 1 finite elements, depicted in chapter 4.

A.1 CALCULUS OF BASIS FUNCTIONS

Basis functions are functions associated to each vertice of the mesh. (Fig. 4.5). The i th basis function ϕ i associated to the i th vertices V i is defined as follows:

Therefore, we can restrict the computation of ϕ i to the elements adjacent to vertex V i . This computation can be done using the following method:

• let e be an element of the mesh, defined by vertices V 1 , V 2 , ..., V n

• the restriction of the basis function ϕ 1 associated to vertex V 1 on the element e is defined as follows:

where the points r , V 1 , V 2 , ..., V n are expressed in their homogeneous notation.

A.4 CALCULUS OF EVOLUTION EQUATION SYS- TEM

The evolution equation:

gives us the following system to solve:

where P i , R i and B i are vectors of size n associated to the vertex i. In the case of only one regular Q 1 element with one of its faces belonging to the boundary ∂Ω of the domain Ω, the vectors P i , Q i and B i can be computed using the previous Q 1 basis functions. For the 2D and 3D case, the vectors associated to the vertex 0 are defined as follows:

The reformulation of this system is more trivial than the previous one since R i , B i and P i are no longer matrices but vectors. Therefore, we won't present in detail their computations. However, we proceed in the same way as for the distance constraint, such as each element in the mesh contributes to the equations induced by its vertices by looking at all the possible pairs of vertex.

However, note that the vectors B i depends on a vector -→ n normal to the boundary, and is thus only defined on the boundary of the mesh. Only the elements having a face included in the boundary of the domain will contribute this term. Therefore, the system A.4 can be computed very easily, by looking at all elements in the mesh, and for each of them adding its contribution to its vertices, by looking at all its possible "edges".

ENERGY DERIVATION OF THE PIECEWISE SMOOTH MODEL

This appendix describes the derivation of the local mean energy proposed in chapter 5, by using the shape gradient method.

B.1 DERIVATION OF THE ENERGY

We consider the following energy, for the bi-partitioning case:

We recall the main theorem of shape gradient method presented in [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients[END_REF]:

Theorem B.1.0.1. The G âteaux derivative of a functional J (Ω) = Ω g(x, Ω)dx in the direction of a vector field V is:

Where g s (x, Ω, V) is the shape derivative of g (x, Ω) in the direction of V, ∂Ω is the boundary of Ω, N is the unit inward normal to ∂Ω and da its area element.

We finally get the following gradient descent:

B.2 IMPLEMENTATION

Each term in the gradient descent can be seen as a convolution by a Gaussian kernel of variance ρ. Therefore, the computation of the evolution can be done in a very fast way by using efficient Gaussian filter scheme. By noting (. * .) | Ω the convolution operator on a domain Ω, we have, We can compute the domain convolution for the bi-partitioning case in a level-set formulation:

where Φ denotes the level set function separating the two regions Ω and Ω, and H σ is a regularized heavyside function.

Experimentations have shown that the auxiliary terms q(x, Ω) have indeed an influence in the evolution of the contour, especially concerning the energy