Keywords: Markov Random Fields, Registration, Segmentation, Brain Atlas, Diffuse Low-Grade Gliomas Champs de Markov Aléatoires, Recalage, Segmentation, Atlas du Cerveau, Gliomes Diffus de Bas Grade 7

Quenette with whom I was always happy to cross path and have a chat.

Outside the lab, I would like to thank my dear friends who have supported me during all those years and made sure I had a social life,

First of all, I would like to express my sincere gratitude towards my advisor Nikos Paragios for his continuous support and advice during my whole PhD. It has been a long journey, and I certainly wouldn't be where I am now if it wasn't for him. Then, I would like to thank Intrasense and more specifically its CEO Stéphane Chemouny, as well as Hugues Duffau for trusting me with this project and for the freedom I was given during this thesis. I will always remember the incredible opportunity I was given to observe a brain surgery at the beginning of my thesis, which was a truly fascinating experience.

Second, I would like to express my deepest gratitude towards all the members of my thesis committee: the reviewers Alison Noble and Christos Davatzikos, the chairman Nicholas Ayache and the examiners William Wells III and Emmanuel Mandonnet, for taking time to read and evaluate my work and their insightful comments and questions.

I would like to thank William Wells for hosting me in the Surgical Planning Laboratory and the three great months I spent there. I was an incredibly interesting and fruitful experience that completely renewed my motivation for my work. I would also like to thank all the members of the lab for their warm welcome and the fun times we spent together.

I thank Amélie Darlix for the enormous work she carried out constructing the DLGG database and her help and advice on the clinical aspect of my work.

I would like to acknowledge Bo (aka Crazy Bobo) who has been a wonderful (the best!) colleague and friend throughout those years. Thank you for all the experiences we shared together, I am very glad that I had you by my side during the tough (such as our first MICCAI deadline) and the happy times. I would also like to thank Blandine for being such a great friend, I am thankful that we started our crazy PhD adventure almost simultaneously. To a larger extend, I thank all my great labmates for making this lab such a better place (and for eating my million cakes): Katerina, Aris, Loïc, Olivier, Chaohui, Panos, Pierre-Yves, Nicolas, Vivien, Sylvain, Fabrice, Haithem, the Stavroses, Nacho, Enzo, Eduard, Puneet, Krishna, Siddharta, Dimitris, Grigoris, Evgenios, Stefan, Mateuzs, Matthew, Iasonas and Pawan. Thank you for the countless times you have helped me with work, the Chinese and belated Greek lessons, the tea/coffee breaks, the sweets from all over the world, the dinner parties... I would like to thank Sylvie and Annie from MAS laboratory, as well as Carine and Natalia from CVN for their invaluable help and support as well as the great discussions we had together. I would also like to acknowledge my former ECP teachers Moncef Stambouli and Florence Mayo-Résumé L'objectif principal de cette thèse est la modélisation, compréhension et segmentation automatique de tumeurs diffuses et infiltrantes appelées Gliomes Diffus de Bas Grade. Deux approches exploitant des connaissances a priori de l'ordre spatial et anatomique ont été proposées. Dans un premier temps, la construction d'un atlas probabiliste qui illustre les positions préférentielles des tumeurs dans le cerveau est présentée. Cet atlas représente un excellent outil pour l'étude des mécanismes associés à la genèse des tumeurs et fournit des indications sur la position probable des tumeurs. Cette information est exploitée dans une méthode de segmentation basée sur des champs de Markov aléatoires, dans laquelle l'atlas guide la segmentation et caractérise la position préférentielle de la tumeur.

Dans un second temps, nous présentons une méthode pour la segmentation de tumeur et le recalage avec absence de correspondances simultanés. Le recalage introduit des informations anatomiques qui améliorent les résultats de segmentation tandis que la détection progressive de la tumeur permet de surmonter l'absence de correspondances sans l'introduction d'un à priori. La méthode est modélisée comme un champ de Markov aléatoire hiérarchique et à base de grille sur laquelle les paramètres de segmentation et recalage sont estimés simultanément. Notre dernière contribution est une méthode d'échantillonnage adaptatif guidé par les incertitudes pour de tels modèles discrets. Ceci permet d'avoir une grande précision tout en maintenant la robustesse et rapidité de la méthode. Le potentiel des deux méthodes est démontré sur de grandes bases de données de gliomes diffus de bas grade hétérogènes (en taille, forme et apparence). De par leur modularité, les méthodes proposées ne se limitent pas au contexte clinique présenté et pourraient facilement être adaptées à d'autres problèmes cliniques ou de vision par ordinateur. Each tag correspond to a functional area that is identified through cortical stimulation, (c) Postoperative MRI. Images extracted from [Duffau 2006] . . . . . . . . . 1.3 Modeling of the LGG's growth rate. (a) Linear growth rates of 27 patients, displaying a mean growth rate of 4mm/year. (b) Growth rate evolution of a patient after anaplastic transformation. Figures extracted from [Mandonnet 2003] and [Mandonnet 2008] respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1

List of Figures

Introduction

This thesis introduces methods for the study and monitoring of brain tumors and in particular the diffuse WHO grade II gliomas or Diffuse Low-Grade Gliomas (DLGG). This is carried out by focusing on two challenging and important problems in medical imaging: Image Segmentation and Image Registration. The segmentation task consists of contouring an object of interest in an image in order to provide its exact position. The task becomes increasingly difficult as the object of interest's appearance gets more heterogeneous and has irregular boundaries. The complemenary task of image registration aims at aligning two or more images by finding corresponding structures between them. The difficulty resides in the relative appearance of the images and fuzziness of the structures' edges. This chapter features a presentation of the clinical context associated to DLGGs and the motivation and main objectives of this thesis, followed by an overview of the thesis organization.

Context and Motivation

Primary brain tumors constitute a heterogeneous set of tumors associated with variable behaviors, symptoms, origins and malignancy. They are characterized by the type of cell from which they arise and are dominated by gliomas, arising from glial cells. Gliomas are classified by the World Heatlh Organisation (WHO) into 4 grades of malignancy based on their histological properties [START_REF] Louis | [END_REF]].

Grade I gliomas (mainly pilocytic astrocytomas) are benign tumors, as they are well delineated and non infiltrative. They have an excellent prognosis after complete surgical resection. Diffuse grade II gliomas are slow-growing infiltrative tumors with continuous growth [Mandonnet 2003]. They inevitably progress into a malignant tumor. Grade III gliomas (or anaplastic gliomas) are fast growing malignant tumors while grade IV gliomas (glioblastomas multiforme) are the most aggressive tumors with a low treatment success rate. They are unfortunately the most common type of glioma.

Diffuse Low-Grade Gliomas: Introduction

DLGGs are diffusely infiltrative primary brain tumors. They represent about 15 % of primary brain tumors [Rigau 2011, Dolecek 2012], and affect younger patients (slightly more males than females) than the most common glioblastomas (mean age about 38 years old) [Capelle 2013]. They are classified into 3 histological subtypes according to their predominant cell type: astrocytomas, oligodendrogliomas and oligoastrocytomas.

Most of the time, DLGGs are revealed by seizures that occur in about 80 % of the cases. Some mild symtoms such as headaches, nausea, speech disturbance can appear depending on the size and location of the tumor. This absence of symptoms rendered the treatment of the tumors very controversial as a surgical tumor resection could cause functional deficits that could be more important than the ones produced by the tumor itself. In the past years, a "wait and see" policy was predominantly considered, where the tumors' growth was simply followed by successive MRIs [START_REF] Whittle | What is the place of conservative management for adult supratentorial low-grade glioma? In Low-Grade Gliomas[END_REF]]. It has been observed however, that during their natural course, the tumors inevitably progress into higher grade fast growing tumors (anaplastic transformation into WHO grade III gliomas) which lead to strong neurological deficits and ultimately to the patient's death (median survival is about 10 years). Therefore, a LGG has to be monitored closely in order to adapt the treatment to its growth rate and corresponding prognosis.

Therapy and Functional Reshaping

Nowadays, a consensus has been reached that maximal safe surgical resection is the first therapeutic option for DLGG [START_REF] Soffietti | Guidelines on management of low-grade gliomas: report of an EFNS-EANO* Task Force[END_REF]]. Indeed, several studies reported an improvement in the prognosis of progression-free and overall survival of patients after surgery [Sanai 2008, Capelle 2013, Jakola 2012]. Tumors localized in or near an important functional area used to be considered inoperable. However, it has been shown that some patients could recover without permanent deficit after damages of eloquent areas [START_REF] Duffau | [END_REF], Duffau 2006]. Other studies show that functional recovery is better in cases of slow growing lesions (like the LGG) than in cases of fast growing lesions (such as strokes) [Desmurget 2007]. It is inferred that this is due to a functional reshaping of the brain induced by the slow growing lesions, which has no time to occur when the damage is rapid. This phenomenon is also known as brain plasticity. Different patterns of functional organization with respect to the tumor can be considered. Firstly, the function can persist within the tumor (36% of patients present intra tumoral activity in [Desmurget 2007]). In that case, it is quasi impossible to completely remove the tumor without sequelae. Another possibility is the redistribution of eloquent areas immediately around the tumor. A near total resection is possible then, but yields transient deficit after the surgery. Thirdly, regions implicated in the functional network in the same hemisphere but remote to the damaged area can be recruited. Finally, the other hemisphere can be solicited. This reshaping would explain the absence of symptoms even when the tumor is located in eloquent areas.

Functional imaging is routinely used for surgery planning in order to evaluate the relationship between the tumor's location and the brain's functional areas and determine the extent of resection Figure 1.1: Preoperative functional mapping. The functional areas are reorganized around the tumor's boundaries. [Duffau 2006] possible. However, it is not sufficient in order to evaluate if the eloquent area is essential or compensable, while precision decreases around the lesion [Holodny 2000]. A solution is to perform brain mapping during the tumor surgery. Direct stimuli are sent on the surface of the brain, which enables the mapping of the motor functions and the cognitive areas (if the patient is awake). If any functional disturbance is noticed, the resection is stopped immediately in that area. The tumor is then removed according to the functional boundaries. Radiotherapy and chemotherapy can be used after surgery for patients with large or fast growing residual tumors.

Follow-up: Tumor Segmentation

In addition to surgery planning, medical imaging plays a key role in DLGG therapy planning and evaluation. Since the tumors are mostly asymptomatic, their growth is monitored by serial MRI exams (usually every 3 months). The best way to follow the tumor's progression is to observe its mean diameter that is expected to follow a linear law [Mandonnet 2008]. Currently, two methods exist for evaluating the mean diameter. The simplest approach consists in manually measuring the tumor's three biggest diameters (in the sagittal, coronal and axial planes) and approximating the tumor's volume as an ellipsoid. The alternative is to compute the tumor's volume using manual segmentation. The mean diameter is finally obtained from the volume as Each tag correspond to a functional area that is identified through cortical stimulation, (c) Postoperative MRI. Images extracted from [Duffau 2006] Mandonnet 2003, Mandonnet 2008]. Despite being fast and convenient, the three diameters method often overestimates the volume of the tumor, especially in cases of irregular contours (which is often the case for low grade gliomas due to their infiltration properties). It can yield inexact results (especially during follow ups after resection, when the shape of the tumor is very different from an ellipse). Therefore, the 3D manual segmentation is currently the gold standard to evaluate the mean diameter. It is, however, tedious and time consuming and suffers from a high intra and inter operator variability [Kaus 2001].

D mean = (2V ) 1/3 [
Typically, the mean diameter progresses according to a linear law. The slope is constant (mean value 4 mm/year) until anaplastic transformation, where there is an abrupt change of speed (average of 8 mm/year). The evaluation of the growth rate is essential for LGG management as it is directly associated to the patient's prognosis and risk of anaplastic transformation. Pallud et al. [Pallud 2012] proposed to systematically evaluate the growth rate and adapt the treatment accordingly. Such evaluation was based on measuring the growth rate between two temporal MRIs separated by 4 months. A growth rate higher than 2mm/year triggers a surgical treatment. An alternative was to use Chang's score [Chang 2009] that is based on 4 risks factors: 1) tumor located in an eloquent area, 2) Karnosky Performance Status ≤ 80 (measure of functional impairment, a score higher than 80 implies no impairment), 3) age > 50 years old, 4) tumor's maximum diameter > 4cm. If the patient obtains a score of 3 -4 a second MRI is performed after 6 weeks. A mean tumor diameter growth larger than 1 mm suggests a growth rate greater than 8.7 mm/year implicating a behavior similar to the one of a high grade glioma. For a score lower than 3, the second MRI is to be performed after 3 months and the mean diameter is measured. The treatment is adapted accordingly.

This strong dependence on the growth diameter as well as the limitations of manual segmentation and ellipsoid approximation methods highlight the need for a fast and robust (i.e. independent of the observer) automatic segmentation method, that is the main objective of this thesis. [Mandonnet 2003] and [Mandonnet 2008] respectively.

Low-Grade Gliomas' Localization in the Brain

As stated earlier, the treatment of low-grade gliomas strongly depends on their localization. It was inferred in [Duffau 2004] that they tend to appear in preferential locations that correspond essentially to functional areas with a predominance in the sensory-motor area and the insula (that studies tend to associate with language functions [Ackermann 2004]). This observation is the starting point of the second objective of this thesis, which is the construction of a statistical atlas mapping the low-grade gliomas preferential locations. Such an atlas offers diverse possibilities of understanding the tumors better. First, the position of the tumors could be correlated with extrinsic factors in order to identify risk factors associated to the DLGG. Second, comparison with functional areas would allow to study the mechanisms of brain plasticity associated to the slow growing tumors. Eventually, we will show that mapping the tumors in the brain provides a useful tool for tumor detection and segmentation.

Thesis Overview

This thesis is divided into 3 main chapters. Chapter 2 presents an overview of the different brain tumor segmentation methods existing in the literature, coupled with an introduction to the image registration problem for atlas-based segmentation.

In chapter 3, the construction of a low-grade gliomas probabilistic atlas is presented. The atlas enables to identify the tumors' preferential locations in the brain while describing their spatial extension in each preferential location. The second part of this chapter presents an application of the atlas for brain tumor segmentation, using it as a spatial prior on where the tumor is expected to be.

Chapter 4 introduces a concurrent registration and segmentation method that simultaneously segments the brain tumor, and registers a healthy subject to the diseased subject with a tumor while taking into account the missing correspondences area represented by the tumor. The chapter presents an alternative method for tumor segmentation that exploits anatomical cues via the registration process, as well as a more precise registration scheme for a future instance of the probabilistic atlas. An extension of the method for the registration of preoperative and intraoperative MRIs during a brain tumor surgical resection is be presented in appendix A.

Chapter 5 concludes the thesis with a discussion on the presented work and its implications, and suggests future research directions with respect to the presented methods.

Chapter 2

Brain Tumor Segmentation

Introduction

Image segmentation is one of the most important and active research area in the medical imaging domain. It can be defined as the delineation of one or several structures of interest within the image. Automated methods are sought in order to avoid the time consuming burden of manually contouring the structures. The problem is particularly difficult in the context of brain tumors. Indeed, most tumors have heterogeneous appearances and their intensity range overlap with the healthy tissues'. The presence of a necrotic core is frequent (especially for glioblastomas, but it also occurs for DLGGs) resulting on a strong contrast with the "active" tumor. Prior information regarding the shape of the tumor cannot be used as they have variable sizes and shapes. DLGGs in particular, have very fuzzy and irregular boundaries due to their infiltrative nature. Edema (swelling of brain tissue around the tumor) and mass effect (tissue displacement induced by the tumor) are quite uncommon due to the slow-growing nature of the DLGGs [Sanai 2011]. In this context, the simplest segmentation methods such as thresholding or region growing are insufficient [Gibbs 1996]. Despite extensive and promising work in the tumor segmentation field, obtaining accurate and reliable segmentations of brain tumors remains a difficult task.

Segmentation methods can be grouped in two categories: surface and region-based approaches. The objective of surface based methods is to find the organ or tumor's boundary by propagation a curve/surface with a flow that is determined according to curvature and image constraints (generally the image gradient). Snakes and level sets are typically used in this context. The former defines the object's boundary explicitly as a parametric curve, while the latter defines the contour via an implicit function allowing for more complex geometries and topological changes.

Region based methods consider the segmentation problem from a different angle. Here, the goal is to identify all voxels belonging to the object and separate them from the rest of the image. Early work relied on fuzzy clustering to regroup similar voxels but quickly showed limitations. The remainder of this chapter is organized as follows: section 2.2 consists of an overview of the different medical imaging modalities used for brain imaging on which the segmentation algorithms are applied, followed by a review the different existing tumor segmentation methods (including, but not limited to DLGG segmentation methods) in the remaining sections. Section 2.3 will be dedicated to surface based methods while Section 2.4 will present the region based methods for tumor segmentation. Section 2.5 will focus on Atlas based tumor segmentation, a method that relies on the registration of a healthy anatomical or probabilistic reference image to incorporate prior anatomical knowledge in the segmentation process and has been getting more and more attention in the recent years.

Brain Imaging

There exists a variety of different imaging modalities that enable the study of the brain. This section presents a brief overview of the different imaging methods, then focus on Magnetic Resonance Imaging (MRI) which is the most common modality for brain tumor observation.

Imaging modalities

Brain Imaging modalities can be grouped in two categories according to the information being captured by them: structural and functional imaging. The structural or anatomical modalities aim at visualizing the different structures and tissues of the brain. Among them, the most popular for Neuroimaging studies are the Computed Tomography (CT) and MRI. CT imaging relies on X-ray technology that is based on the absorption of X-rays beams as they pass through the different tissues of a patient's body. CT scans are constructed by using a series of X-ray beams that rotate around the patient's head. Each beam yields a 2D image at a specific angle that are used to construct a 3D volume by tomographic reconstruction. CT scans are well contrasted and high resolution images. The main drawback of this technique is the the fact that it uses ionizing radiation.

MRI relies on the magnetic property of the hydrogen nuclei present in large quantity in the human body. It is the modality of choice for brain studies due to its high tissue contrast and details and non ionizing property. In the recent years, new techniques have developed such as Diffusion Tensor Imaging (DTI). This modality enables the reconstruction of the white matter tracts connecting the different parts of the neural networks in the brain by measuring the anisotropic diffusion of water inside the tissues. The presence of a tumor could have a direct impact on the fibers structures by causing disruption, displacement or infiltration of the fibers [Wei 2007]. Modalities such as MR spectroscopy (measurement of major metabolites in tumor tissue) and Perfusion MRI (measurement of the relative cerebral blood volume using a contrast agent) provide additional informations for the diagnosis and study of brain tumors [START_REF] Soffietti | Guidelines on management of low-grade gliomas: report of an EFNS-EANO* Task Force[END_REF]].

The goal of functional imaging is to study the human brain's function based on physiological changes caused by the brain's activity. Electroencephalography (EEG) and Magnetoencephalography (MEG) are techniques that offer means of directly measuring the brain activity. EEG detects the electrical impulses in the brain due to the neuronal activity via electrodes placed on the scalp. MEG measures the magnetic flux changes using sensors positioned closed to the scalp. The methods are popular due to their simplicity, the fact that they are non invasive and their very high temporal resolution. Determining precisely the spatial origin of the observed signal is however difficult.

Position Emission Tomography (PET) and Single Photon Emission Computed Tomography (SPECT) are nuclear imaging techniques that measure changes in the cerebral blood flow and tissue metabolism and how it is altered by brain disorders. A biological molecule is marked with a radioactive isotope and injected in the bloodstream and accumulated in areas where the molecule has affinity. The advantage of this method is that the radioactive tracer can be designed specifically to target specific organs or processes related to an illness. For instance, a radioactive isotope typically used is the Fluorodeoxyglucose (FDG) that behaves like glucose molecules and therefore is able to trace the brain's metabolic activity. However, the technique is invasive and potentially harmful due to the use of radioactive tracers while at the same time producing the isotopes is costly and difficult. Nowadays, CT and PET technologies are being combined in one device, allowing to combine the anatomical information recovered from the CT scans with the metabolic information given by the PET scans.

Functional MRI (fMRI) exploits the Blood Oxygen Level Dependent (BOLD) contrast [Ogawa 1990] to detect changes in the neuronal activity induced by sensorimotor or cognitive tasks. Neuronal activation causes an increase of the blood flow to compensate the oxygen consumption and therefore reduces the amount of deoxygenated hemoglobin molecules. The detection of changes in hemoglobin oxygenation relies on the paramagnetic properties of the deoxygenated hemoglobin which impacts the measured NMR signal. fMRI is often used for tumor surgery planning in order to identify the spatial relationship between the lesion and the functional area and evaluate the corresponding risks. Time constraints in the acquisition of the signal cause the fMRI images to be of lower quality than the structural MRI.

Functional imaging is particularly important for brain tumor surgery planning. It enables to establish a relationship between the lesion and the functional area and plan the surgery accordingly (extent of resection, awake craniotomy...).

Magnetic Resonance Imaging

Magnetic Resonance Imaging is one of the most popular medical imaging modalities. It has the advantage of being a non invasive and non ionizing technique that produces images with excellent tissue contrast, making it the modality of choice for brain (and more specifically brain tumors) study. MRI is based on the principle of Nuclear Magnetic Resonance (NMR), that is commonly used in spectroscopy to study the physical and chemical properties of molecules.

Physical principles

MRI relies on the magnetic properties (nuclear spin) of the hydrogen nuclei present in high quantities in the human body. The hydrogen nuclei spins around its axis, leading to the creation of an electromagnetic field. In natural condition, the spins' orientation is random and therefore, the total magnetic field is null. When placed into a large magnetic field B 0 , the nucleis' spins align along (positive spin) or against (negative spin) the external field. The spins precess around the magnetic field's axis with a precession frequency or Lamor frequency ω 0 proportional to the external field:

ω 0 = γB 0 (2.1)
The magnetic vector of spinning nuclei can be expressed through the longitudinal component (M z , parallel to B 0 ) and the transverse component (M xy , perpendicular to B 0 ). In practice, the number of positive spins is slightly higher than the number of negative spins, yielding a longitudinal magnetization M z . Since the spins do not precess in phase, the transverse magnetization is null. When a radiofrequency (RF) pulse B 1 is introduced at the resonance frequency ω 0 , there is an exchange of energy between the nuclei and the RF pulse which affects the spin equilibrium (phenomenon known as resonance). This causes the net magnetization vector to tilt away from the longitudinal axis and the creation of a transverse magnetization that can be measured by a receiver coil. In practice, the pulse B 1 is perpendicular to B 0 , which causes the longitudinal component to be null at resonance and the transverse signal to be the strongest possible.

When the RF pulse is turned off, the nuclei progressively returns to equilibrium (relaxation phenomenon). During relaxation, the absorbed electromagnetic energy is retransmitted and constitute the NMR signal. Two different kinds of relaxation processes are observed: longitudinal and transverse. The longitudinal relaxation corresponds to the spin-lattice (surrounding tissue) interaction. The spin returns to its equilibrium state (alignment with the magnetic field B 0 ) by giving energy to the lattice. The recovery of the longitudinal magnetism follows an exponential curve and is characterized by a tissue specific time constant T1 after which 63 % of the final value is recovered.

M z = M 0 (1 -exp(- t T 1 )) (2.2)
M 0 is the net magnetization at equilibrium that depends on the proton density and the strength of the external magnetic field. The transverse relaxation describes the spin-spin interaction and corresponds to the spins getting out of phase. This yields an exponential diminution of the transverse component characterized by a tissue specific time constant T2.

M xy = M xy0 exp(- t T 2 ) (2.3) M xy0
is the amplitude of the transverse signal following the RF pulse. When the RF pulse is perpendicular to B 0 , M xy0 is equal to M 0 . In practice, T2 is always smaller than T1. The two constants are the key for the image contrast and definition of the different MRI sequences.

Image Sequences and DLGGs Imaging

An MRI sequence is characterized by two parameters. The Echo Time (TE) is the time between the RF pulse and the signal measurement, and the Repetition Time (TR) is the time between two RF pulses. After each RF pulse, the received NMR signal shows the evolution of the transverse component during relaxation. The signal oscillates at the resonance frequency and its envelope decays following an exponential rule. The initial amplitude corresponds to the magnitude of the longitudinal component and impacts its recovery status. The signal will be strong if the longitudinal component has been well recovered. TR defines the time between two pulses and therefore the time allowed for longitudinal recovery. If TR is short, the longitudinal magnetization of tissues with a long T1 will not have time to recover completely and their initial signal will be weaker, while the signal of tissues with a small T1 will be higher (large magnetization). This defines the T1-weighted modality, where the image contrast between tissues is defined by their T1 constant value. The image contrast no longer depends on TR when TR is set long enough so that complete longitudinal recovery is achieved for all tissues (same initial magnitude M 0 ). A long TE renders the image dependent on the transverse relaxation time T2. The signal will be weaker for tissues that have faster transversal decays (small T2 value). The T2-weighted images are obtained in this manner. Eventually, T1 and T2 are no longer discriminant parameters when TR is long and TE is short. This yields the Proton Density weighted image, where the strength of the signal depends on the strength of the acquired magnetization M 0 (dependent on proton density). Brain tumor segmentation algorithms often exploit the different modalities to include maximal information. Visualization of brain tumor on T1 or T2 images alone is far from ideal. The tumors appear hypointense or normal on T1 images, and are therefore hard to detect and delineate. Better detection can be achieved by the injection of contrast agents, the most common being gadolinium contrast agents. It results in a tumor that is hyperintense in T1 images. However, complete en- hancement is far from systematic especially in the context of DLGGs. About 20 % of DLGGs display a patchy and faint contrast enhancement. Contrast enhancement can also be associated to anapestic transformation. Therefore, methods that depend on contrast enhancement can only segment a subset of tumors and very few DLGGs. Contrarily to T1 images, the tumor appears hyperintense on T2 images and can be well visualized. The contrast is poor though, due to the high signal intensity of fluids which can be similar to the tumor's intensity.

Another popular sequence for the observation of brain tumors is the T2 FLAIR (Fluid Attenuated Inversion Recovery) sequence. This sequence is constructed by initially flipping the longitudinal magnetic field M z in the opposite direction using an additional electromagnetic pulse oriented at 180 • from the original field B 0 . This pulse is followed by a longitudinal relaxation to return to the equilibrium, passing through the null value. An additional time constraint TI (Inversion Time) determines when the perpendicular RF pulse is added. It is possible to suppress certain tissues' signal by choosing TI so that it corresponds to the time when the signal magnitude is null. In FLAIR sequences, TI is set in order to suppress the high cerebrospinal fluid (CSF) signal in T2 images. The FLAIR sequence is one of the best contrasted modalities to visualize brain tumors and is the gold standard for DLGG observation [START_REF] Bynevelt | [END_REF]].

An additional difficulty inherent to MR images is the problem of intensity inhomogeneity. The problem, caused by imperfections in the image acquisition process (instrumentation limitations or patient motion), is observed from one acquisition to the next and even within the same acquisition (bias field corruption) [START_REF] Sled | [END_REF]]. This phenomenon can cause strong misclassifications when the segmentation technique relies on the observed intensity and/or on training data. Many algorithms have been developed in order to correct those inhomogeneities [Nyúl 1999[START_REF] Vovk | A review of methods for correction of intensity inhomogeneity in MRI[END_REF]], but the problem becomes even more difficult in the presence of a pathology since it can have similar visual behavior as bias field corruption.

Surface Based Segmentation: Deformable models

Deformable models have their origins on Snakes or Active Contours in [START_REF][END_REF]]. The object contour is represented as an explicit parametric surface S : R → R 3 that evolves from an initial configuration towards the optimal solution according to internal and external energy forces. The internal energy imposes smoothness and regularity of the curve while the external energy is based on the image features, guiding the curve towards locations of large gradient (object boundaries).

The optimal contour is obtained by minimizing the energy functional. Despite the appeal and popularity of the concept, it suffers from several drawbacks, notably its sensitivity to initialization and tendency to get trapped in local minima, and inability to handle topological changes between the initial and the final configuration due to the explicit representation.

The level set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]] was introduced in medical imaging in [Malladi 1995]. It constitutes an elegant way of dealing with the inherent limitations of point-based snakes representations and in particular the handling of topological changes. It consists of an implicit representation of the object contour as the zero level set of a function Φ(x, t) of higher dimension. The implicit function is designed so that it is positive outside the contour and negative inside, and commonly set as the signed distance from the contour.

S = {x|Φ(x, t = 0) = 0} (2.4)
This yields an attractive formulation for shape modeling [START_REF] Leventon | [END_REF], Chen 2002[START_REF] Rousson | Shape priors for level set representations[END_REF]] and image segmentation that handles topological changes as well as complex geometries and does not depend on the contour's parametrization. The evolution of the surface follows the surface's normal direction with a speed F, and is governed by a Partial Differential Equation: The design of the speed function F is an essential aspect of the framework and one of the strongest challenges of the method. Such a speed function can be either the natural outcome of a minimization process or can be directly introduced depending on the nature of problem being considered. In the most general case it is the sum of two terms: an image term F Im and a curvature term.

∂Φ ∂t = F (t)|∇Φ| (2.5)
F (t) = F im + γ∇ ∇Φ |∇Φ| (2.6)
The mean curvature term is a smoothness term which aim is to prevent leakage of the level set in noisy and indistinct areas. The parameter γ determines the amount of smoothing and is of key importance as a value too high would lead to under segmentation. Despite extensive work in the medical imaging domain, the use of level sets for the segmentation of brain tumors is limited. The main drawbacks of the level set framework is its sensitivity to initialization and parameter setting as well as the computational burden of solving the partial derivative equation. In order to enable interactive parameter tuning, Droske et al. [Droske 2001] proposed an adaptive multi grid framework allowing for almost real time feedback on the curvature terms and front evolution speed. They adopt the Fast Marching Algorithm [Sethian 1996]. Assuming the speed function is always positive and the surface is initialized inside the tumor, the evolution of the contour will always be outwards. This allows to compute the arrival time T (x) of the front on a voxel x. This converts the problem into solving the so-called Eikonal equation:

|∇T |F = 1 (2.7)
The method is accelerated by considering an adaptively sampled grid (the highest resolution corresponding to the voxel level) based on the local image homogeneity evaluated using the image gradient. Similarly, [START_REF] Lefohn | [END_REF]] used GPU programming to reduce the running time, enabling interactive handling of the free parameters. Both methods used simple image terms based on image gradient and contrast that are sensitive to image noise and unclear boundaries, which enforces the need for expert supervision. Ho et al. [Ho 2002] proposed an fully automatic 3D level sets framework that relies on region based methods [START_REF] Zhu | [END_REF], Chan 2001, Paragios 2002]to avoid manual initialization and parameter tuning. An initial separation of the (enhanced) tumor and background pixels is performed by fuzzy classification, used for initialization of the surface. Classification is performed by fitting two distributions to the histogram of the pre and post gadolinium enhanced T1 weighted MRI images difference.The derived probabilities are then used to guide the evolution of the level set function in a region competition manner:

F im = α(P tum (x) -P bg (x))
(2.8)

The method was applied to glioblastomas and meningiomas. The performance of the method depends heavily on the quality of the classification and while it also requires contrast enhancement.

A similar region competition level set framework was presented in [START_REF] Popuri | 3D variational brain tumor segmentation using Dirichlet priors on a clustered feature set[END_REF]]. A series of features are extracted from the images (multi modal intensities, Gabor features for texture information and brain hemispheres symmetry) and used to perform an initial classification of the image using K-means clustering [MacQueen 1967] (see section 3.2.2). The contour evolution is guided by the clusters value probabilities combined with a Dirichlet prior which role is to penalize clusters that are predominant in the healthy brain and wrongly assigned to the tumor region. A post processing step to remove CSF and skull misclassification was added. The method was tested on 15 multi modal MRIs (T1, T1 contrast enhanced (T1CE), and T2).

Taheri et al. [START_REF] Taheri | [END_REF]] combined adaptive image thresholding with level sets segmentation. The level set was initialized inside the tumor area, leading to an initial guess of the threshold value. The threshold value was then iteratively updated based on the current segmentation. The contour evolves based on the signed intensity difference with respect to the threshold value, and stops when the variation of the threshold from one iteration to the next is lower than a predefined parameter. This approach assumes relatively strong intensity difference between tumor and background.

Deformable models are attractive methods that have proven to handle complex geometries and produce smooth and precise segmentations. They suffer however from several drawbacks that reduce their interest for the difficult problem of tumor segmentation. Their main weakness is their sensitivity to initialization and parameter selection. Interactive parameter tuning and initialization is promising, but implies extensive user work and is dependent on the user's knowledge. The combination with region based approaches enables automation but renders the final segmentation strongly dependent on the quality of the classification scores.

Region-Based Segmentation

Region-based segmentation methods aim at directly identifying the tumor voxels instead of finding the tumor's surface. We distinguish unsupervised classification that is solely based on the image's appearance and properties, from supervised classification where the separation criteria are learned from a labeled training set. Both approaches have been considered for tumor segmentation and will be reviewed in this section.

Classification methods

Unsupervised Classification and Knowledge based Methods

Unsupervised classification relies on intrinsic properties of the samples. Those methods were among the first considered due to their simplicity and the fact that classification is performed solely based on the observed samples without the need for labeled training data. The concept of clustering is to regroup observations into groups (clusters) with similar properties. Early work in tumor segmentation relied on clustering to classify the image into different categories. We distinguish hard clustering where each sample belongs to only one cluster from fuzzy clustering, where each sample is assigned a membership score for each cluster.

Phillips et al. [START_REF] Phillips | [END_REF] proposed fuzzy clustering to perform automated glioblastoma segmentation. The method adopted is the Fuzzy-C-means (FCM) algorithm [Dunn 1973, Bezdek 1981]. This is an unsupervised clustering technique derived from the well known hard clustering K-means algorithm, with the difference that the cluster membership is probabilistic. Fuzzy partitioning is carried out by minimizing of the sum of weighted distances between the samples and the cluster centers:

J m = n 1 K 1 u m ij x i -c j 2 (2.9)
where K is the desired number of clusters (tissue classes), c j is the center of cluster j and u ij is the degree of membership of sample x i in cluster j that verifies:

C i=1 u m ij = 1 ∀j n j=1 u m ij > 0 ∀i (2.10)
After initialization of the memberships, the clustering is obtained via an iterative process, alternating between computation of the cluster centers (as the mean of all points weighted by their membership) and of the cluster memberships. The samples were clustered based on their multispectral intensities values (T1, T2 and Proton Density). Numerous misclassifications were observed, which revealed the strong intensity overlap between the healthy tissue and the tumor as well as the insufficiency of the algorithm alone to perform meaningful segmentation. Clark et al. [Clark 1998] proposed a knowledge based automated segmentation method on multispectral MRI data for the 2D segmentation of enhanced glioblastomas. After learning an ensemble of rules on a training set of 17 slices, an initial segmentation of each slice is performed using Fuzzy C-Means clustering. The predefined rules enable to identify the abnormal slices to be further processed. After adaptive histogram thresholding based on the pixels mutli spectral intensities, the obtained segmentations are post processed based on the different heuristics. A similar system was developed in [START_REF][END_REF] for non enhancing tumors by defining a different set of heuristics.

The quality of the segmentation depends strongly on the definition of the knowledge system, which should be general yet precise enough to obtain good segmentations for a broad range of tumors. The need for important expert guidance to define the system as well as the limited success on noisy or heterogeneous tumor have reduced the appeal of such approaches and highlighted the need for more elaborate methods.

Supervised Classification

Despite the theoretical and practical value of unsupervised classification, it has quickly shown limitations in the separation of object with complicated and heterogeneous appearances such as tumors. As a result, supervised classification techniques have become increasingly popular in the segmentation field. The idea is to learn a classifier from a set of pre labeled data and an ensemble of image features. Consider an ensemble of training pairs {x i , y i } where x i corresponds to an observation (typically an image voxel) and y i to a label ({-1, 1} in the context of binary classification). Each pair is associated to a feature vector f (x i ). Based on these training samples, a mapping function is learned and constitutes what is called a classifier. In its simplest form, the set of image features can the image intensities. Supervised classification involves a training phase during which the classifier is learned based on the labeled training set, followed by a testing phase during which a new unlabeled sample is assigned a label based on its features and the classifier's rules.

Vinitski et al. [Vinitski 1997] used multispectral (T1,T2 and Proton density) k Nearest Neighbor (kNN) classification for the segmentation of benign and malignant brain tumors. Each voxel is associated with a 3 dimensional vector corresponding to the intensities observed in the 3 different modalities. The training set consists of samples from the scan to be segmented (voxels) that are manually selected and labeled by an operator. A probability per class is then assigned to each unlabeled voxel based on the classes of its k closest samples, the k parameter being set manually. This relatively simple method requires manual selection of samples on every image to segment and The most popular supervised classification method for brain tumor segmentation has undoubtedly been the Support Vector Machines (SVM) [García 2004, Verma 2008[START_REF] Zhang | [END_REF]]. Introduced in [Boser 1992] for linear binary classification, it has become main stream in machine learning and pattern recognition. Given a set of training samples {x i , y i }, the SVM method seeks a hyperplane {x| < w, x > +b = 0} defined by its normal vector w and the offset b, that separates the two different classes. The optimal hyperplane is chosen so that it maximizes the margin between the plane and the closest samples. The Support Vectors are the samples that are the closest to the hyperplane and therefore determine its position. In practice, the data is rarely linearly separable. Non linear separation is performed by mapping the data into a higher dimensional feature space x → Φ(x). The optimal hyperplane is then found by minimizing the following objective function:

min w,ξ,b 1 2 w 2 + C n 1 ξ i s.t. y i (w • Φ(x i ) -b) ≥ 1 -ξ i (2.11)
where ξ i are slack variables introduced to allow mislabeled samples, perfect separation being rarely possible when dealing with noisy data such as medical images. This optimization problem is solved using the Lagrangian formulation, where unknown variables α i (lagrange multipliers) are introduced for each constraint. The support vectors are the samples that have non zero weights α i .

L n = i α i - 1 2 i j α i α j y i y j Φ(x i )Φ(x j ) (2.12)
Finding the appropriate mapping can be difficult, while the complexity of the algorithm augments with the size of the feature space. It becomes unnecessary when defining a kernel function

Margin

Separating hyperplane Support Vectors

Figure 2.9: Illustration of the SVM approach.

K(x i , x j ) = Φ(x i )Φ(x j ) that directly computes the dot product of the feature vectors. Using such a formulation, the problem becomes solely dependent on the number of samples, it is therefore possible to operate with a feature space of infinite dimension. Eventually, the classification decision of a new sample is recovered as:

h(x) = i α i y i K(x i , x) + b (2.13)
The SVMs being by definition binary classifiers, extension to multiclass classification is not straightforward. Current methods either use a hierarchical approach or a one vs all approach [START_REF] Duan | [END_REF]].

Another popular classification algorithm is the AdaBoost algorithm. Boosting is a general machine learning method which goal is to improve the accuracy of any learning algorithm. Ad-aBoost (short for Adaptive Boosting) is the most popular boosting algorithm and was introduced in [START_REF] Freund | [END_REF]]. It relies on the concept that a classifier with high discriminative performance (strong classifier) can be constructed as a weighted linear combination of weak classifiers. The most commonly used weak classifiers' structures are decision stumps or tree structured classifiers. The simplest classifier is preferred to maintain the speed of the algorithm. The weak prediction rules are learned iteratively by minimizing the classification error with an increased focus on misclassified samples at each new iteration. Several variants of the Adaboost algorithm have been developed, among them the Gentle Adaboost algorithm [Friedman 2000] that will be further detailed in section 3.3.1. The popularity of the algorithm can be explained by its simplicity, flexibility and speed. An example of AdaBoost being applied for tumor segmentation can be found in [Xuan 2007].

Recent work on Decision Forests (or Randomized Forests) showed strong promises [Zikic 2012, Geremia 2012]. The concept was introduced in [START_REF] Breiman | Random forests[END_REF]] and corresponds to a collection of tree structured classifiers (decision trees). A decision tree consists of a set of binary rules that progressively split the data into smaller sets based on their features values. Starting at the root node, the complete training data X is split and assigned to either the left X L or right X R child node. The split is chosen among a randomly selected set of possible splits regarding the feature values so that it maximizes the information gain:

I = - X R X H(X R ) - X L X H(X L ) (2.14)
where H(.) is the entropy calculated from a histogram of class values. The same process is repeated until the tree reaches a predefined depth d. Each leaf node (terminal node of the tree) is associated with a probability that corresponds to the fraction of the data of class c that has been assigned to the node. Several decision trees are learned in a similar fashion and the forest consists of the ensemble of those trees. During testing, the test data is pushed through all the decision trees by applying the different separation rules. Each tree yields a probability per sample and class based on the reached leaf node. The forest probability in then computed as the average probability over all the trees. Decision forests are inherently multi-class classifiers, robust to overfitting and fast, making them an appreciated classification approach.

A key element of the SVM, Adaboost or Decision forest algorithms is the design of the feature vector that should discriminate between healthy and diseased samples. The most straightforward features are intensity based, such as the multi modal intensities, intensity patches around the sample (to incorporate neighborhood information) on which first order statistical moments are often computed (mean, variance or entropy to name a few). Features related to the image texture have also been popular, among them the Haralick features [Haralick 1973], or the Gabor features [Manjunath 1996]. Xuan et al. [Xuan 2007] exploit the symmetry between the brain hemispheres by computing the intensity difference between the left and right hemisphere. Zikic et al. [Zikic 2012] use preliminary classification results obtained by modeling the data as a mixture of Gaussians and context aware features such as the intensity range along a 3D line. In the specific case of SVMs, kernel mapping can be used instead of (or combined with) the feature vector.

Despite the promises of supervised pattern classification methods, they are constrained by the i.i.d (independent and identically distributed) assumption that treats all samples independently. This leads to noisy and irregular segmentations that could be barely compensated by the introduction of neighborhood or context aware information in the feature vector or post processing with morphological operators. Those limitations have motivated the use of more elaborated models that take into account the spatial interaction between neighboring pixels. Graph based approaches offer an elegant and flexible way of doing so, modeling the spatial interactions as edges connecting the pixels.

Graph based Segmentation

The main idea behind graph based segmentation is to consider the image as an undirected graph G = (V, E), that consists of a set of nodes V = {1, 2, .., N }, corresponding to each voxel of the image, and a set of edges E = (i, j) ∈ V, that represent the dependencies between the nodes. Typically, the edges only connect the adjacent nodes (4-neighborhood in 2D, and 6-neighborhood in 3D).

Affinity based Segmentation

The most natural way of exploiting such a graphical formulation is to partition the graph into connected components of similar properties. In such methods, the edges encode an affinity measure w i,j that represent the similarity between the connected nodes i and j. The problem of graph partitioning can be seen as the search for m disjoint sub graphs {V 1 , .., V m }, m being the number of sought classes, that verify V = ∪ m i=1 V i and ∩ m i=1 V i = ∅. The subgraphs are constructed by cutting the edges between them and constitute ensembles of vertices with strong affinities.

Considering the binary case, the cost of partitioning the graph into subsets A and B can be computed as the total affinity of the removed edges, called the cut of the graph:

cut(A, B) = i∈A,j∈B w i,j
(2.15)

The optimal bipartition of the graph minimizes this cut value. The minimum cut method was used in [START_REF] Wu | [END_REF]] for image segmentation as a clustering problem, seeking the smallest k -1 cuts among all possible cuts to form k disjoint subgraphs. Despite promising results, it highlighted the tendency of the method to favor smaller subgraphs leading to inadequate segmentations. This led to the introduction of Normalized Cuts [Shi 2000] that compute the cost of a cut as a fraction of the total edges connections to all nodes in the graph:

N cut(A, B) = cut(A, B) assoc(A, V) + cut(A, B) assoc(B, V) (2.16)
where assoc(A, V ) = i∈A,j∈V w i,j is the sum of the connection weights from all nodes in A to all the nodes in the graph. The optimal Ncut can be approximated by eigenspectrum analysis [Shi 2000].

Affinity based segmentation has drawn limited interest for tumor segmentation. [Moonis 2002[START_REF] Liu | [END_REF] proposed interactive affinity based segmentation methods that rely on fuzzy theory. After manual selection and labeling of seeds, the affinity between neighboring voxels is computed. The smallest affinity value represents the strength of connectivity of a path from voxel i to voxel j. Eventually, the strength of connectivity between both voxels is set as the largest path among all possible paths from i to j. The tumor is segmented as the voxels that are connected to the tumor seeds with a strength larger than an empirically chosen threshold.

Corso et al. [START_REF] Corso | [END_REF]] proposed a hybrid approach between Bayesian classification and graph partitioning for segmentation of glioblastomas. The affinities are defined in a probabilistic fashion so that they are class dependent:

P (X ij |i, j) = c i c j P (X ij |i, j, c i , c j )P (i|c i )P (j|c j )P (i, j)
(2.17) where X ij is the binary event describing if nodes i and j are assigned to the same class. The likelihood probabilities P (i|c i ) are learned as mixtures of Gaussians and lower the affinity if the class assignment is unlikely. The first term P (X ij |i, j, c i , c j ) is the affinity function weighted by a class dependent parameter learned over a training set. Graph partitioning relies on the Segmentation by Weighted Aggregation (SWA) algorithm, a multilevel hierarchical approach proposed in [Sharon 2006] for Normalized cuts partitioning. An initial high resolution graph G 0 (one node per voxel) is progressively coarsened by aggregating nodes that have strong affinities. Weights are updated by interpolation and edges are constructed if the weights are positive. The algorithm stops when the coarsened graph is made of one single node or if there are no connections between the remaining nodes. At each level, structures of different scales are extracted. In [Sharon 2006] they are recovered by associating each voxel with the most salient segment (lowest normalized cut) of which it is a part in the pyramid. This segment recovery method has been that successful in the context of medical images. Corso et al. [START_REF] Corso | [END_REF]] proposed computing the maximum likelihood at each level weighted by the inter-level interpolation weights. The method was tested for the segmentation of multi spectral glioblastoma images (T1 pre and post contrast, T2 and FLAIR) into 4 classes (tumor, edema, brain and background).

Markov Random Fields

MRFs where first introduced into vision in [START_REF] Geman | [END_REF]] for image restoration, and have since been one of the most popular undirected graphical model in the domain. They have been successfully applied to a broad range of applications [Kolmogorov 2002, Boykov 2006, Komodakis 2007b, Glocker 2008a]. Let X = {X 1 , ...., X n } be a set of random variables defined on V, each variable X i taking a value l i in a label space L. A joint event {X 1 = l 1 , ..., X n = l n } is called a configuration and abbreviated as X = l, where l = {l 1 , ..., l n }. We denote P (X = l), , abbreviated as P (l), the joint probability of a configuration. A random field is a Markov Random Field on a graph G if the two following conditions are satisfied:

Positivity: P (l) ≥ 0 ∀l ∈ L Markov Property: P (l i |l V-{i} ) = P (l i |l N (i) )
(2.18)

N (i) is the neighborhood system defined as: N (i) = {j|i, j ∈ E}. The Markov property imposes that the conditional distribution of a variable is only dependent on its neighbors. The problem of image segmentation is reformulated as a labeling problem, where the goal is to assign a label to each node (voxel). The labels correspond to the different image classes. For binary classification, the label set is defined as L = {0, 1}. The optimal labeling should maximize the joint configuration probability:

l = argmax l P (l) (2.19)
The Hammersley-Clifford theorem [Besag 1974] establishes the equivalence between an MRF and a Gibbs distribution, allowing to link the local conditional probability to a global probability on the whole distribution:

P (l) = 1 Z exp - 1 T C V C (l C ) (2.20)
Where T is a constant called the temperature (holdover from the formulation's origins in statistical physics) typically set to 1 and Z is a normalizing constant called the partition function:

Z = L exp - 1 T C V C (l C ) (2.21)
C denotes a clique that is a fully connected subgraph and V C (.) is an energy function defined on the subset defined by the clique, called the clique potential. The estimation of the partition function is computationally intractable. It is not critical for inference, but becomes important for learning since it depends on the model's parameters. It can be avoided for inference by working in the energy domain and minimizing the energy function:

E(l) = C V C (l C ) (2.22)
Considering cliques of size up to two leads to the well known pairwise MRF energy:

E(l) = i∈V V i (l i ) + i∈V j∈N (i) V i,j (l i , l j ) (2.23)

MAP-MRF Framework

The MAP-MRF framework relates to Bayesian statistics, where one seeks an optimal labeling given an observed image I by maximizing the Maximum a Posteriori probability P (l|I). Following Bayes rule, the posterior probability is linked to the joint probability as:

l = arg max l P (l|I) ∝ P (I, l) = P (I|l)P (l) (2.24)
where the prior probability P (l) is modeled as an MRF. Typically, only pairwise cliques are considered, but higher order cliques can also be considered to incorporate more elaborate prior knowledge. Conditional independence of the voxels is assumed given their label class ( P (I|l) = i (P (I(i)|l i )) [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]]. The MAP estimate of the posterior probability is equivalent to minimizing the corresponding energy function:

E(l) = i -log(P (I(i)|l i )) + i j∈N (i) V i,j (l i , l j ) (2.25)
The singleton potential is the data term or likelihood potential that measures the cost of assigning a label. The pairwise potential is a regularization prior ensuring local consistency and smoothness of the labeling. The MRF formulation allows flexibility in the design of the energy potentials, which are the key aspects of the model. One of the most simple and popular spatial interaction model is the homogeneous and isotropic Potts model (or Ising model for binary labelings) that penalizes different segmentation labels between neighbors:

V i,j (l i , l j ) = β(1 -δ(l i , l j )) (2.26)
where δ(., .) is the Kronecker delta function and the parameter β determines the amount of regularization. The likelihood potential is usually modeled using pattern classification techniques, the role of the MRF prior being to counterpart the i.i.d assumption of those approaches. Conditional Random Fields (CRF) [Lafferty 2001] differ from MRFs by the fact that they are globally conditioned on the observation I. They are thus discriminative models that directly model the conditional posterior probability, allowing to incorporate information from the observed variables with more flexibility. Notably, it allows the definition of a pairwise regularization term that is dependent on the observation.

MRF Optimization

Inference of MRFs/CRFs is a challenging task and many algorithms have been developed to tackle the problem. Among the earliest methods, we can cite the Iterated Conditional Modes (ICM) [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]] and simulated annealing (SA) [Kirkpatrick 1983]. ICM performs iterative local optimization of the posterior probability. This makes a relatively efficient algorithm that is however strongly sensitive to initialization. Simulated annealing methods are, on the contrary, more likely to reach global optima but computationally intractable methods. Nowadays, efficient methods have been developed that can be separated in two categories: graph cuts based approaches and message passing methods.

Graph cuts convert the problem of pairwise MRF optimization for binary labeling into a graph partitioning one by the construction of a directed graph (the s-t graph) where two terminal nodes (source s and sink t) corresponding to the two different labels have been introduced. An s-t cut separates the graph into two disjoint sets S and T, each connected to one terminal node. Nodes are assigned the label associated to the terminal node they are connected to. The graph is constructed so that the cost of an s-t cut (sum of removed edges' weights) is equal to the MRF energy given the corresponding configuration l. The problem is rendered equivalent to the min-cut problem introduced in section 2.4.2.1. The theorem of Ford and Fulkerson [Ford 1962] is an important result that states that finding the min-cut is equivalent to finding the maximum flow from the source s to the sink t (max-flow problem). The max-flow saturates a set of edges that corresponds to the min-cut. Low complexity methods have been proposed to solve the min-cut/max-flow problem, among which augmenting path methods [Ford 1962] and push-relabel methods [Goldberg 1988]. A review and comparison of both approaches can be found in [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF]]. Two limitations of the method are the restriction to binary labels and the requirement that the pairwise potentials are submodular in order to reach a global optimum:

V pq (0, 0) + V pq (1, 1) ≤ V pq (0, 1) + V pq (1, 0) (2.27)
For non submodular functions, a method referred to as Quadratic Pseudo-Boolean Optimization (QPBO) was considered [Kolmogorov 2007]. QPBO constructs an s-t graph where each node in the original MRF is represented by two complementary nodes x i and xi = 1x i . Solving the graph's min-cut/max-flow problem yields three types of label: 0,1 or unknown if the label of nodes x i and xi is inconsistent. QPBO provides a partial solution with guarantee that the labeled variables correspond to the optimal solution. Extensions to the multi label case were proposed in [START_REF] Boykov | [END_REF]]. They proposed the αexpansion and αβ swap methods that iteratively optimize the MRF energy by defining a set of possible moves based on the initial MRF configuration at each iteration, and select the configuration that locally minimizes the energy. α-expansion allows any node to change its label to α and only α. αβ swap only allows nodes labeled α to change their label to β and vice versa. Each iteration corresponds to a binary s-t cut computation. The local minimum recovered by α-expansion is guaranteed to be within a known factor of the global optimum.

In [Komodakis 2007b, Komodakis 2008b] a linear programming based method called Fast-PD was proposed for optimization of MRFs. It relies on linear programming relaxation and primaldual gap minimization. The method can be seen as a generalization of the α-expansion method and is substantially faster. The algorithm is used for MRF optimization throughout this thesis and a short review of the method will be provided in the next chapter.

Message passing techniques propagate messages across the nodes to evaluate the optimal configuration. Belief Propagation (BP) [START_REF] Pearl | Judea Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF]] is an iterative algorithm where at each iteration, a node p sends a message m pq (l) to node q ∈ N (p) describing how likely the assignment of label l to node q is. The message is computed based on the messages node p previously received and the local energy costs. Two kinds of BP algorithms exist: the max-product that estimates the MAP configuration and the sum-product that computes the marginal probabilities. BP is exact for tree structured graph where messages are passed from the leafs to the roofs. It is approximate where the graph contains loops (loopy belief propagation) and there is no guarantee of convergence.

Sequential Tree Reweighted message passing (TRW-S) [Kolmogorov 2006] is a message passing algorithm strongly related to linear programming relaxation. Here, the solution is sought by maximizing the dual problem (lower bound) of the relaxed linear integer program. This is done by splitting the problem into convex subproblems defined as trees. Then, the algorithm considers each node sequentially and performs belief propagation in each tree it belongs to followed by an averaging step. In [Komodakis 2007a] a new message passing method was proposed based on dual decomposition. It relies on decoupling the Lagrangian dual problem into several subproblems through Lagrangian relaxation. The decoupled subproblems (the "slaves", defined as subtrees of the MRF graph) are easily solved using the max-product BP algorithm. The method iterates between solving the simple subproblems and sending the obtained minimizers to the "master" which role is to update the slaves parameters (using the projected subgradient method) so that the dual is maximized.

Applications to Tumor Segmentation

An interactive hierarchical MRF model was presented in [START_REF] Gering | [END_REF]]. An initial intensity based segmentation at the voxel level is progressively refined by incorporating higher level knowledge. The voxel based segmentation relies on estimation of the MRI bias field using the Expectation-Maximization (EM) algorithm [Wells III 1996b]. The EM algorithm [START_REF] Dempster | [END_REF]] is an iterative procedure for the estimation of the maximum likelihood of probabilistic models in the presence of missing data (also called latent variables). In this work, the missing data are the tissue classes. The healthy classes intensities (likelihood) are modeled as Gaussian distributions, while the tumor class is modeled as a deviation from normalcy (minimum intensity distance from all other classes). The bias field is modeled as a zero-mean additive Gaussian distribution. The EM algorithm iteratively alternates between two steps: the E-step corresponds to the estimation of the posterior probabilities given the currently estimated parameters (bias field and Gaussian parameters here) while the M-step is a MAP estimation of the bias field given the estimated tissue probabilities. Initialization of the Gaussian parameters is done by the user by selecting representative voxels of each class.

The next level refines the obtained classification using an MRF model, whose parameters are learned from training data. Additional multi level MRF models are then used to refine the classification results based on structure size (expected distance from the structures boundaries) and inter structure distances.

In [Wels 2008], the likelihood term was constructed by Probabilistic Boosting Trees classification (PBT) [START_REF] Tu | Probabilistic boosting-tree: Learning discriminative models for classification, recognition, and clustering[END_REF]] and smoothed with an Ising prior weighted by the relative voxels' intensities. The idea of PBT is to progressively diminish the complexity of the learning process through a recursive construction of the classification tree. At each tree node, a strong classifier is learned using the AdaBoost algorithm. The data is then split and assigned to two child nodes based on their classification results. This progressively reduces the heterogeneity of the data at each node and simplifies the classification task. The method was optimized using Graph Cuts and tested for 3D segmentation of 6 multi spectral (T1, T1CE, T2) pediatric MR images.

Lee et al. [Lee 2005] proposed to combine SVM classification with CRF smoothing for tumor segmentation. The smoothing term was modeled as a penalizing constant between neighboring nodes weighted by the nodes' feature vector computed for the SVMs:

V i,j (l i , l j , x)) = l i l j ν T γ i,j (x)
(2.28)

where γ i,j (.) is the weighting function describing the similarity between nodes i and j, and ν T is a node dependent parameter that has to be inferred. The idea is to encourage discontinuities between neighboring nodes if they display different properties. After learning of the CRF parameters, the optimization was performed using the ICM approach. The method was tested on 7 patients using patient specific training for the classifier (the training set consists of extracted slices of the image to be segmented). To cope with the difficult and time consuming task of learning a CRF's parameters (that requires the computation of the partition function), the same group proposed the Pseudo Conditional Random Fields (PCRFs) [START_REF] Lee | [END_REF]]. A logistic regression classifier's output is smoothed using two parameter free pairwise potential functions (one for label interaction and the second for voxel interaction). The learning task is thus reduced to learning the logistic regression's parameter. Inference was performed using graph cuts. The system was tested for multispectral (T1, T1CE, T2) 2D segmentation of 11 contrast enhancing tumors and showed similar results with their previous CRF formulation with the benefit of substantial decrease on the learning time. Bauer et al. [Bauer 2011] also coupled SVM classification with CRF smoothing. They use a two stage smoothing, first a 3D regularization that is orientation dependent (weaker in the z axis due to anisotropic data) to separate tumor from background and second a stronger 2D regularization on a highly connected neighborhood (8 neighbors in 2D) to perform a 6 classes classification (CSF, White Matter, Grey Matter, Necrosis, Active Tumor and Edema). They exploited recent advances in MRF/CRF optimization using the Fast-PD algorithm, yielding a very fast algorithm. The method was tested on 10 brain tumor patients using multi spectral data (T1, T1CE, T2 and FLAIR).

To conclude, graph based segmentation methods are attractive due to their ability to recover good optima as well as their modularity and scalability. Markov Random Field models in particular enable to couple inter node interaction with individual node classification properties. Further-more, they are computationally efficient thanks to the recent advances in discrete MRF optimization techniques. Although prior knowledge can easily be introduced in an MRF formulation, current methods are limited to local priors that model the dependencies between neighboring nodes. Hierarchical models and position dependent smoothing (CRFs) offer better flexibility and stronger knowledge but still lack global structural information on the brain.

Atlas based Segmentation

Atlas based segmentation methods have gained a lot of attention in the recent years for automatic segmentation of medical images mostly due to the progress made in the field of deformable registration [START_REF] Sotiras | [END_REF]]. They allow the introduction of prior anatomical knowledge in the segmentation process by registering (i.e. align) an annotated anatomical template or a probabilistic atlas to the volume to be segmented (target volume). The method is particularly interesting for the segmentation of brains due to their relatively limited inter patient variability compared to other organs. When the structures to be segmented are present both in the atlas and target volume, the task of segmentation is converted into a registration task (the segmentation by registration method) [START_REF] Bondiau | [END_REF]]. A different approach has to be adopted when the object to segment is a tumor that is not present in the atlas. This section will first introduce the problem of medical image registration that is the basis of atlas construction and atlas based segmentation, followed by a review of the different reference atlases used for brain segmentation and their construction. Eventually, the existing applications of atlas registration to tumor segmentation will be reviewed. A more detailed review of some key deformable registration methods, applied to the context of brain tumors, will be presented in chapter 4.

Medical Image Registration

Medical Image Registration is, concurrently to image segmentation, among the most active research areas in the domain. It is defined as the process of aligning two (or more) images by establishing point to point correspondences. Given a floating image J and a fixed image I defined on a domain Ω, the aim is to find the transformation T that will map the floating image to the fixed image:

I(x) = J • T (x) (2.29)
We differentiate monomodal (same image modality) and multimodal (two different modalities) registration, as well as intra-patient (the two images belong to the same patient) and inter-patient registration. Intra-patient monomodal registration is used for temporal follow up. The objective is to compare the patient's image taken at different times and evaluate changes, particularly the evolution of a pathology with or without treatment. Registration can be necessary due to the patient's motion or changes in the image acquisition process. The different imaging modalities display complementary informations but are obtained in different reference systems. Aligning multimodal images renders direct comparison possible such as the association of structural with functional information. Inter-patient registration refers to atlas construction and alignment to a new patient for segmentation or population studies. The registration process is typically carried out by minimizing a similarity criterion that evaluates the differences between the images to match:

T = argmin T Sim(I, J • T ) (2.30)
The problem is characterized by 3 key elements: the image features to match and corresponding matching criteria Sim(., .), the transformation model that describes how the image will be deformed, and the optimization technique used to recover the optimal transformation.

Matching Criteria

When trying to align two images, we rely on the concept that both images have the same structures and try to find their correspondences. Two kinds of approaches can be considered to match the structures: geometric and iconic registration. Geometric approaches define a set of key points in the image corresponding to structural or geometrical landmarks, or extract surfaces from the pre segmented structures' boundaries. The matching criteria to be minimized is defined as the distance between corresponding points or surfaces in both images. Geometric methods have the advantage of relying on the anatomical structures and yield transformations that are more faithful and preserving of the anatomy. On the other hand, the registration is as good as the landmarks whose automatic selection is not trivial. Iconic (or intensity based) approaches seek to optimize a similarity measure based on the voxels' intensity values. The measure is dependent on the intensity relationship between the two images The underlying principle is that when the images are matched, their intensities should verify a certain relationship. The simplest relationship is the identity. The images' intensities are assumed to have the same or similar intensities in corresponding structures. This relationship is considered when dealing with monomodal images and allows direct comparisons of the images intensities. The most popular measures are the Sum of Absolute Differences (SAD) and the Sum of Squared Differences (SSD).

A more sophisticated approach for monomodal registration is to assume an affine relationship between the images intensities to acknowledge the potential relative intensity changes between the images to be registered. In this setting, the criterion of choice is the correlation coefficient [Brown 1992, Kim 2004]:

CC = Cov(I, J • T ) V ar(I) V ar(J • T ) (2.31)
Those monomodal similarity criteria directly compare the intensities of the two images. However, as it is the case for image classification, intensity alone is not that discriminative and can lead to ambiguous matchings and convergence to bad local minima solutions. This issue has motivated the work of Shen et al. [Shen 2002], where they proposed to use feature or attribute vectors that encompass neighborhood information and local anatomical information based on pre established tissue segmentations.

The definition of a matching criterion becomes more difficult when dealing with multimodal images. Given the impossibility to estimate the relationship between multimodal intensities, it was suggested in [START_REF] Roche | [END_REF]] to assume an unknown functional relationship I ≈ f (J • T ). Given this functional assumption, they proposed to use the Correlation Ratio as a similarity measure:

CR = V ar(E(I|J • T )) V ar(I) (2.32)
Information Theoretic approaches assume no existing relationship between the intensities, making them more general and very popular. Since its introduction, the Mutual Information (MI) [Wells III 1996a, Viola 1997[START_REF] Collignon | [END_REF][START_REF] Maes | [END_REF]] has been widely used (see [Pluim 2003] for a review of MI based registration methods) and is probably the most commonly used multimodal registration criterion. The objective of this measure is to maximize the mutual information of two images which is a combination of their joint and individual entropies. The Shannon entropy of a distribution p can be interpreted as a measure of the distribution's dispersion and is defined as:

H(x) = - x p(x)log(p(x))
(2.33)

The joint probability distribution p(I, J • T ) of two images is estimated from the normalized joint histogram of the images intensity values, while the marginal distributions p(I), p(J • T ) are obtained by summing over the rows and columns of the joint histogram. When the images are aligned, the joint histogram will be sharp and focused, while its dispersion augments as they get misaligned. As a result, minimizing the joint entropy corresponds to aligning the images. This is the idea behind the mutual information defined as:

M I(I, J • T ) = H(I) + H(J • T ) -H(I, J • T ) (2.34)
Given this definition, the term -H(I, J • T ) implies that maximizing the mutual information is related to minimizing the joint entropy. The presence of the individual entropies makes it more robust than direct minimization of the joint entropy: transformations that reduce the amount of information in an image are penalized. The Normalized Mutual Information (NMI) was later introduced as a measure less sensitive to the amount of overlap between the images:

N M I(I, J • T ) = H(I) + H(J • T ) H(I, J • T ) (2.35)
Multimodal registration remains a difficult problem that is actively studied. Aside from information theoretic methods, a second kind of approach consists of transforming the problem into the simpler monomodal problem. This can be achieved by simulating one modality into the other [START_REF] Roche | [END_REF], Wein 2008], or mapping both images into a new space where they are directly comparable [Maintz 2001[START_REF] Lee | [END_REF], Bronstein 2010].

Transformation Model

The transformation model describes how the image will be deformed and the amount of deformation allowed. The computational complexity as well as the precision of the deformation are dictated by the transformation model. The rigid transformation is the simplest transformation model. It only allows rotation and translation that, in the 3D case, amounts to 6 parameters (degrees of freedom) that need to be estimated. Similarity/Affine registration allow scaling and shearing additionally to translation and rotation (parallel lines are preserved), which corresponds to 12 degrees of freedom in 3D. Rigid and affine transformations are essentially used for same subject registration of organs that display limited movements or deformation in time such as the healthy brain. When stronger deformations are observed, they are often used as preprocessing prior to a dense registration.

If no constraints are imposed on the transformation, the registration problem is referred to as non-linear, dense, curved or deformable registration. The problem is ill-posed with unstable and non smooth solutions. To tackle this issue, a priori knowledge is introduced in the minimization problem in the form of a regularization term:

E(T (x) = Sim(F, G • T (x)) + λReg(T (x))
(2.36)

The parameter λ is a constant that weights the influence of the regularization over the similarity measure. The regularization term ensures smoothness of the deformation field and preservation of the anatomical properties. Two kinds of deformation models can be considered in this setting: parametric transformations and physics based transformations. In both cases, the transformation is usually defined at each pixel position as an additive displacement:

T (x) = x + u(x).
Parametric models represent the transformation as a linear combination of basis functions. Those methods are interpolation based: the transformation is calculated on an ensemble of control points that are then used to produce a dense interpolated flow for the entire image. Their interest reside in the fact that the model's complexity (correlated to the number of control points) can be considerably reduced while the richness of the transformation is preserved. The most popular parametric transformation model is the Free Form Deformation model (FFD) [Sederberg 1986[START_REF] Rueckert | [END_REF]]. Here, the control points are placed uniformly across the image with uniform spacing δ and the basis functions used are cubic B-splines. The deformation field is computed as:

u(x) = 3 l=0 3 m=0 3 n=0 B l (u)B m (v)B n (w)d i+l,j+m,k+n
(2.37)

where i = [x/δ x ] -1, j = [y/δ y ] -1, k = [z/δ z ] -1, u = x/δ x -[x/δ x ], v = y/δ y -[y/δ y ] and w = z/δ z -[z/δ z ].
B l represents the lth basis function of the B-spline and d i+l,j+m,k+n is the displacement at the corresponding control point. It has been proven [START_REF] Choi | [END_REF]] that the obtained transformation is diffeomorphic if the maximum displacement is limited to 0.4δ.

Physics based deformation models rely on a physical process guided by a force that is associated to the image matching criterion. The deformation model is described by a partial differential equation that is associated to the kind of physical model employed. Among the different physical models used for deformable registration, we can cite elastic registration (stretching of an elastic material) [Bajcsy 1989], viscous fluid models [START_REF] Christensen | Deformable shape models for anatomy[END_REF]] and more famously diffusion models governed by: ∆u + F = 0 (2.38)

The use of registration through diffusion models was introduced in [Thirion 1998] as an analogy to a thermodynamic concept (Maxwell's demons). Demons are mathematical entities located at the boundary of an object in the image that is modeled as a membrane. Demons apply forces to push the deformable model in the normal direction of the membrane according to local image characteristics. For instance, a demon can be placed at each location in the image where it is not constant. In this case, the corresponding membrane is an isosurface. The demons forces are inspired from optical flow equations but are normalized in order to handle instabilities due to small images gradients. This yields the displacement field gradient or velocity as:

v = (I -J • T )∇I ∇I 2 + α(I -J • T ) 2 (2.39)
The method is iterative and alternates between computation of the forces with regularization by Gaussian smoothing. This resulted in a computationally efficient algorithm that lacked a proper theoretical formalism. The first successful attempt to provide theoretical grounds of support for demons can be found in [START_REF] Pennec | [END_REF]], casting the problem as a gradient descent based energy minimization process. In order to recover the algorithm's efficient two step optimization process, it was suggested in [Cachier 2003] to introduce hidden variables that decouple the problem in a matching and regularization step. An extension of Thirion's algorithm to diffeomorphism using this approach was later introduced in [Vercauteren 2008].

Another physics based registration approach that has proven relevant in the presence of brain tumors is the use of biomechanical models to describe the transformation. Such models describe the elastic deformation properties of the different tissues using Finite Element Methods. The main difficulty resides in finding the appropriate model for the anatomical tissues which is far from trivial. It should also be noted that inter subject registration (and by extension, population based atlas construction) is not possible since the deformations from one subject to the next are not natural. Registration methods in the presence of tumors using biomechanical models will be reviewed in chapter 4.

Optimization procedure

The optimization procedure is closely related to the design of the transformation model and similarity metric. Continuous methods seek real valued registration parameters Θ. The optimal set of parameters is determined in an iterative fashion given an initial estimate as:

Θ (t+1) = Θ (t) + α (t) d (t)
(2.40) t) is the search direction at iteration t, and α (t) is the step size along the search direction. Plethora of optimization methods have been proposed that differ in the definition of α (t) and d. Notably, a comparative study of continuous registration strategies can be found in [Klein 2007]. The most famous and simple is the Gradient Descent method that defines the search direction as the derivative of the energy evaluated given the currently estimated parameters. More elaborated strategies involve a combination of the gradient and search direction with the ones at the previous iteration (Conjugate Gradient methods) [Hestenes 1952], computing an approximation of the gradient (Stochastic Gradient Descent) [Kushner 2003] or using higher order derivatives in the search direction (Quasi-Newton methods) [START_REF] Dennis | [END_REF]]. Continuous methods are tuned to the specific registration method and similarity measure for which they are used. The use of discrete optimization methods has gained a lot of attention in the past years. Contrarily to continuous methods, they perform a global search which makes them robust to initialization while benefiting from fast convergence rates. Another strong advantage of using such an approach is the modularity with respect to the similarity criterion and regularization term. The use of discrete methods is associated to the modeling of the registration problem as a discrete pairwise MRF. Given a predefined discrete displacement set, the registration problem is reformulated as a labeling problem where each node (control point given a parametric transformation model) is assigned a label that corresponds to a displacement. The MRF energy is defined so that the unary terms correspond to the similarity criterion, while the pairwise term plays the part of the regularization. The optimization problem then boils down to an MRF optimization. Among the many MRF optimization methods developed, Fast-PD and TRW-S have proven successful in the context of image registration [Shekhovtsov 2008, Glocker 2008a]. The main limitation of discrete methods is the trade off between a precise sampling of the search space and control points and computational efficiency. The issue will be addressed in chapter 4.

d (

Brain Atlases

The investigation of the brain's structural and functional organization has motivated the construction of numerous maps of the brain. Brain atlases have been built to create standardized reference coordinate systems in which individual brain maps can be mapped and compared. One of the first reference system (the stereotaxic or Talairach coordinates) was proposed by [Talairach 1988] and became the standard reference in brain mapping. The stereotaxic coordinates are defined based on two landmarks in the brain: the Anterior (origin of the coordinates) and Posterior commissures (AC and PC respectively). They were defined in order to horizontally orient the line connecting the two landmarks and vertically orient the interhemispheric plane. Alongside with the definition of the atlas, a piecewise affine transformation method (the brain is divided into 6 associated blocks) has been proposed to register a volume to the Talairach atlas. Despite the success of the Talairach atlas, it suffered from severe limitations. It was constructed from the single post mortem brain of a 60 year old French woman that was far from representative of the population, had a variable slice thickness and inconsistent data from orthogonal planes. The developments in medical imaging paved the way for the construction of more representative and detailed atlases. Among them, a 3D T1-weighted MR based atlas (the Harvard Brain Atlas) was constructed by the Surgical Planning Laboratory, using the scan of a 25 years old healthy subject that was enhanced by anisotropic filtering [Kikinis 1996]. It was hand-labeled into 150 different structures. To cope with the lack of contrast and high signal to noise ratio, The Montreal Neurological Institute (MNI) constructed a single subject atlas by averaging 27 high resolution T1-weighted scans of the same individual. [Holmes 1998] Despite offering precise structural information, single subject atlases are not always representative of the population and can introduce a bias. This has motivated the construction of population atlases that are based on an ensemble of scans that are registered and averaged. The reference population atlases have been constructed by the MNI and Laboratory of Neuro Imaging of UCLA (members of the International Consortium for Brain Mapping (ICBM) project). Their first population atlas was the MNI-305 [Evans 1993, Collins 1994], constructed from 239 males and 66 females around 23 years old. All images where registered to the Talairach space and their intensities were normalized. The atlas was a voxel wise average of the intensities over all volumes. Two improved atlases where constructed using the same procedure with 152 (MNI-ICBM152) [Fonov 2009, Fonov 2011] and 452 (ICBM 452) [Mazziotta 2001] subjects. The resulting volumes are blurred where structures vary among subjects and are therefore not suited for structural stud- ies but offer good templates for registration with limited bias. Similar atlases where constructed for T2-weighted, PD-weighted and Diffusion weighted images. Last but not least, probabilistic interpretation of those atlases were constructed by segmenting the different subjects tissues (white matter, grey matter and CSF) after registration. The probabilistic maps are constructed by averaging the labels over all the subjects. Those atlases describe the spatial distribution of tissue classes over large populations and can be exploited as global anatomical priors. Extended reviews of brain mapping and atlases can be found in [START_REF] Toga | [END_REF], Cabezas 2011].

Atlas based Tumor Segmentation

The segmentation of brain tumors cannot be mapped into a registration process due to the absence of a tumor in the atlas to be registered. Existing methods exploit the information brought by the atlas as a spatial prior on healthy brain structures. Kaus et al. [Kaus 1999, Kaus 2001] extended an atlas based segmentation method (the ATM SVC algorithm) [Warfield 2000] for meningioma and low-grade glioma segmentation. The method alternates between supervised kNN classification and registration of a manually segmented anatomical template (single subject) based on the current segmentation. The atlas is initially rigidly registered to the volume to be segmented. The features for kNN classification are the voxel intensities as well as spatial features obtained from the aligned atlas. The spatial features are computed by converting each of the tissue classes into a distance map. Using distance maps allows to have a non optimal registration result since it doesn't rely on hard boundaries but on the distance from where the labels are quite certain. The classification is then followed by an elastic deformable registration on the basis of the segmented image [Dengler 1988]. During the registration steps, the tumor is labeled as healthy tissue. A spatial template for tumor is defined at the next iteration based on the tumor segmentation obtained at the previous iteration. The different brain structures (background, skin, brain, ventricles and tumor) were segmented hier- archically. The method was tested on 20 T1CE images: 7 meningiomas and 14 low grade gliomas. The method inherits the drawbacks of kNN classification and implies strong homogeneity of the tumor appearance and limited spectral overlap with the healthy tissue. The successive applications of registration and classification also make the algorithm slow.

Cobzas et al. [Cobzas 2007] integrated atlas priors in a combined statistical learning/level set framework. The level set is defined in a region competition manner and evolves based on pre estimated probabilities. After registration, the atlas prior was integrated in a high dimensional feature set for supervised statistical learning of the two classes probabilities.

The classification decisions were smoothed by an MRF prior designed as a Potts model. The optimal parameters (distributions mean and covariance, bias field parameters and MRF parameters) and tissue classes were estimated using the EM algorithm. Automation of the method was achieved using a probabilistic brain atlas (ICBM152). After affine registration of the atlas to the subject, the probabilities were used as a rough prior estimation of the classification at every voxel. Furthermore, the atlas probabilities values were used as a spatial prior to increase the quality of the classification. [Moon 2002, Prastawa 2003] built on this work and adapted it to the presence of tumor and edema. For each patient, the probabilistic atlas is modified to include prior probabilities on the tumor and edema positions. The tumor prior probability map is constructed based on contrast enhancement. A difference image between T1 and T1CE images is constructed. A mixture of distributions is fitted to the difference image's histogram: two Gaussian distributions that correspond to noise and errors in registration and a Gamma distribution for contrast enhancement. The probability that corresponds to contrast enhancement is used to convert the difference image into a tumor prior probability.To remove contrast enhanced areas that do not correspond to the tumor, a smoothing regularization process is included using the region competition level set method developed by [Ho 2002].The edema prior class is introduced similarly as a fraction (20 %) of the white matter prior probability. The method was applied to 5 glioblastoma and meningioma cases and used T1 and T2 images for classification.

Another approach, similar to the work of [START_REF] Gering | [END_REF]], was to detect the tumor voxels as outliers with respect to the normal voxels. The concept was proposed by Van Leemput et al. [Van Leemput 2001] as an extension of their previous work by including an outlier class for multiple sclerosis segmentation. Each class probability is here defined as a mixture of a normal distribution and a uniform distribution δ corresponding to the outlier class:

p k (I(x i )|Φ i ) = (1 -ǫ)N (µ k , Σ k ) + ǫδ (2.41)
Where ǫ ∈ [0, 1] represents the portion of outliers. After computation of the parameters and class memberships, a voxel will be assigned to the outlier class if its membership score is below a threshold for all classes. The intensity model employed here for multiple sclerosis is too simple for a direct application to tumor segmentation.

The same idea (model lesions as outliers) guided the work of Prastawa et al. [Prastawa 2004].

They presented a more general framework in comparison to their previous work [Prastawa 2003] that was limited to enhanced tumors featuring edema. The method consists of three steps: first the tumor and edema are detected as outliers from the healthy classes given their intensity profile, then the presence of an edema is detected while the last step refines the segmentation output using spatial and geometric constraints. Affine registration of a probabilistic atlas (ICBM atlas) allows the random selection of a set of training samples representative of each healthy class (high atlas probability for White Matter, Grey Matter and CSF) but that are contaminated with tumor and edema samples (the abnormal class). The healthy multimodal (T1 and T2) intensity distributions are assumed to be Gaussian. The minimum covariance determinant [Rousseeuw 1999] is used to estimate the classes mean and covariances due to its robustness to the presence of outliers. The abnormal class is then formed by rejecting samples that are three standard deviation away from their class distribution. Once the samples are correctly classified, the problem is defined as a MAP problem. The likelihoods are estimated using Parzen kernel density estimation [START_REF] Duda | [END_REF]]:

P (I(x)|C) = 1 N N i=1 K λ (I(x) -I(p i )) (2.42)
where K λ is the Gaussian kernel, λ is its standard deviation or kernel bandwidth that determines the amount of smoothing on the estimated probability and p i is a training sample. The spatial priors for healthy classes are extracted from the atlas, and a fraction of the sum of the white and grey matter probabilities are used as prior for the abnormal class. In an EM style approach, this classification process was jointly considered with the bias field correction method of [START_REF] Van Leemput | [END_REF].

Tumor and edema classes are separated by clustering based on their T2 intensity and spatial constraints are enforced using a region competition level set. Menze et al. [Menze 2010] proposed an algorithm for tumor segmentation that strongly exploits the fact that the tumor has different visibility in multi-channel images. They propose a generative model in which each modality is modeled and segmented independently. Each class (including the tumor class) and modality is modeled as a Gaussian distribution combined with a smoothing MRF prior. A vector of tumor appearance describes the observation of a tumor at each voxel for each modality. An EM style approach is adopted to alternatively estimate the classification probabilities from which are derived the multi channel classifications, and the Gaussian distribution's parameters. The global segmentation is obtained by averaging the classification decisions over all channels. The method is initialized by segmenting the image into 4 classes (White Matter, Grey Matter, CSF and outlier) using the method of Van Leemput et al. [Van Leemput 2001]. The outlier class is defined as being more than three standard deviations away from the centroid of any normal class. The method was tested on 25 glioma patients using T1, T1CE, T2 and FLAIR volumes.

Atlas based tumor segmentation methods have a strong dependency on the quality of the registration. Despite numerous efforts to make the segmentation robust to errors in registration, the methods will fail if the brain is deformed by the tumor. Detection as an outlier also implies that the tumor's appearance varies significantly from the one of the healthy tissue, an assumption that does not always hold given the frequent intensity overlaps. Nevertheless, the use of anatomical and probabilistic atlases is a powerful tool for incorporation of prior spatial information on the expected localization of the tumors.

Summary

Despite an extensive amount of work dedicated to brain tumor segmentation, the problem remains difficult and segmentation results are often not satisfactory. Generalization of the existing methods to several kinds of tumors (or pathologies) is often difficult. Furthermore, few methods have been adapted to heterogeneous DLGG segmentation. The focus of early methods was on contrast enhanced and homogeneous tumors and while glioblastomas have been mostly studied in the recent years. Additionally, many methods are based on the use of different image modalities which are not systematically available in a clinical setting. Despite the need of a prior training phase, machine learning methods offer a flexible and efficient way of capturing the tumor's properties that are not limited to intensity values. Coupled with MRF models, the i.i.d assumption is relaxed and smooth solutions can be obtained very fast due to recent advances in MRF optimization techniques. An additional advantage of MRF modeling is the natural inclusion of prior knowledge in the formulation. Such prior knowledge can be obtained via the use of anatomical or probabilistic atlases.

Next chapter will present a DLGG segmentation method that introduces prior atlas based knowledge in an MRF framework. Instead of detecting the tumor as an outlier with respect to the healthy tissue distribution, a tumor specific atlas is constructed describing the most likely position of tumor voxels.

Chapter 3

A Probabilistic Atlas of Diffuse WHO Grade II Gliomas' Preferential Locations

Introduction

Brain atlases are excellent tools for the study of the brain's anatomy and the variations within populations. Futhermore, they offer a reference template for registration and atlas based segmentation tasks. They are however, not representative of diseased populations since the pathology tends to alter or deform the brain's structures. Several atlases are available with repositories of annotated images for a diversity of pathologies. Such images are a rich source of information but not in the same coordinate space and therefore not directly comparable.

This has motivated the construction of disease specific atlases that reflect the particular clinical subpopulation affected by the disease and how it alters the brain. Such atlases enable to study the evolution and origin of the disease as well as the variations from one patient to another. Aside from aligning all volumes in the same coordinates system, the difficulty resides in obtaining a statistical interpretation of the inter subject variability. The construction of disease specific atlases has notably been one of the focuses of UCLA's Laboratory of Neuro Imaging. Among notable disease specific atlases, Thompson et al. [Thompson 2000] have constructed atlases dedicated to the impact of schizophrenia and Alzheimer's disease on the brain. Maps of anatomical variability are constructed by comparing the deformation fields obtained by deformable registration of all volumes.

In the context of Diffuse Low Grade Gliomas (DLGG), the construction of an atlas mapping the tumors appearances in the brain was motivated by the fact that the management of brain tumors is highly dependent on their location. The impact of the tumor is dependent on the eloquence of the region in which it is. Several studies have reported a functional reorganization of the brain induced by the slow growing tumors [Desmurget 2007, Duffau 2005] while it was inferred in [START_REF] Ius | Evidence for potentials and limitations of brain plasticity using an atlas of functional resectability of WHO grade II gliomas: towards a ¡Èminimal common brain¡É[END_REF]] that certain areas have low compensatory abilities. These observations are directly linked to surgery planning for tumor resection, where the eloquence of the region and its compensatory ability can determine the safety of the surgery and amount of resection possible. Furthermore, tumors could have a location dependent behavior due to local differences in their molecular biology [Gozé 2009, Ren 2012].

An interesting observation was made in [Duffau 2004], inferring that DLGGs tend to appear in preferential locations. This assumption was explored for all gliomas in [START_REF] Larjavaara | [END_REF], Gooya 2012] where it was observed that the distribution of gliomas in the brain is non-uniform. Following this observation, this chapter introduces a probabilistic atlas that maps the DLGGs appearances in the brain and provides their preferential locations. This is carried out via the construction of a complete graph, where each node represents a tumor and the edges represent the spatial proximity between them. The tumors' preferential locations are then identified by extraction of a sparse subgraph, where each cluster corresponds to a preferential location. This results in a DLGG atlas where the central position (cluster center) and spatial extension (position of tumor voxels in the cluster) of each preferential location is represented. An application of this atlas is proposed in the second part of this chapter for the segmentation and characterization (integration in the atlas) of a new patient's tumor with spatial priors. This is carried out via a two-level Markov Random Field formulation based on position priors inferred from the atlas. The tumor is characterized and segmented based on a prior tumor detection through statistical classification, while the knowledge on the tumors' spatial extension in the clusters is used as a prior knowledge on the position of the voxels of the tumor to be segmented.

The remainder of this chapter is organized as follows: first the methodology for the construction of the low-grade glioma atlas is presented in section 3.2, while section 3.3 presents a direct application of the atlas for tumor segmentation. The experimental validation is part of section 3.4 and discussion concludes the chapter.

Spatial Position Mapping of Diffuse Low-Grade Gliomas

To construct the atlas, we assume the existence of a database of patients MRIs featuring a DLGG. The position of each tumor is known through binary maps obtained by experts' manual segmentation. To allow position comparison, all volumes are to be aligned to the same reference template of a healthy brain. This section presents the construction of the complete graph (registration and proximity measure between the tumors) followed by the estimation of its optimal clustering.

Graph Construction

Registration to a Common Reference Space

The mapping of all annotated volumes to the same reference space could be done through deformable registration. However, the presence of tumors violates the correspondence hypothesis, which can cause displacements errors and a shift in the position of the tumor. Discarding the pathology during registration was considered in [Parisot 2011] but as will be shown in chapter 4, such an approach can still lead to errors in registration at the vicinity of the tumor. Despite being less precise, affine registration is the most common choice in atlas construction, because it handles properly linear variations and therefore roughly conserves the anatomy so that the form and location of the tumor is preserved and comparable. The use of affine registration reflects the diversity of the anatomy, which makes the atlas more adapted for the integration and tumor segmentation of a new, unsegmented patient.

We adopt the block matching method of [Ourselin 2000] to affinely register our data set to the same reference template using the MedInria software (http://med.inria.fr/). Instead of globally minimizing an objective function, the idea of block matching is to perform registration in a sequential manner by minimizing a local objective function between image regions (blocks).

For rigid or affine registration, the registration is performed in a multiscale iterative manner where two successive tasks are carried out at each iteration. First, the displacement field is computed between the images by maximizing a local similarity measure. This is done by moving a block of the floating image I and comparing it to blocks in the reference target image J that have coincident positions. The size of the blocks and extent of the search space are constant at a given scale level. The Normalized Mutual Information similarity metric is used in MedInria.

The block that maximizes the similarity measure is chosen, which yields a correspondence between the centers of the blocks (a i , b i ). Inadequate displacement fields that correspond to a similarity measure below a given threshold (uniform blocks or absence of correspondences) are not considered in the subsequent step. The transformation T that best fit the obtained displacement field is computed by least trimmed square regression:

argmin T h i=1 a i -T (b i ) 2 h = n 2 (3.1)
where n is the number of displacement vectors. An illustration of the block matching search process is shown in Fig. [3.1] All volumes in the data set are registered to the same reference template that is an arbitrarily chosen healthy subject, following this scheme. The obtained displacement fields are then applied to the segmentation binary maps. The next step consists in measuring the proximity between the different aligned tumor segmentation maps.

Proximity Between Tumors

When measuring the proximity between tumors, there are two issues to be considered: (i) the registration errors and (ii) the different sizes of the tumors. Here, we are not interested in the size and extent of the tumor, but solely on its position in the brain reference space. Measuring the distance between tumors as surfaces or point clouds could penalize small tumors that would be further away from its neighbors than a large tumor. Furthermore, considering surfaces makes the distance measure more sensitive to the registration and miscorrespondences that are bound to happen with the simple affine registration scheme. This has motivated the use of the tumors' center of gravity as a reference position, computed by averaging all tumor pixels coordinates. The proximity between the tumors is then computed as the Euclidean distance between the centers1 .

The distance is computed between each pair of images. This constitutes an adjacency matrix that is used to construct a complete graph G LGG = (V LGG , E LGG ) where each node corresponds to a patient/tumor and is spatially associated to the center of gravity of the tumor, while the edges connect all pairs of nodes (p,q) ∈ V LGG and are weighted by the euclidean distance between tumors centers d(p, q). The obtained complete graph can be observed in Fig. [3.2], where the color and thickness of the edges represent their strength (from red and thick to blue and thin). From this complete representation, the existence of high and low density areas for tumors can be observed, most notably around the temporal lobes.

Graph Clustering

The next step towards identifying the DLGGs' preferential locations is to regroup nearby nodes based on their proximity. This can be done by performing clustering on the graph. The concept of clustering, as introduced in section 2.4, is to regroup the data into segments so that samples within the same cluster are similar to each other and dissimilar with respect to the elements associated with other clusters. In this section, we will only consider "hard" clustering, where each sample is assigned to a single cluster. In our context, nodes in the same cluster should be close and present in the same anatomical area. The clusters made of the highest number of nodes would then correspond to the preferential locations, while an absence of cluster in a region suggests an extremely rare location for tumor appearance. Obtaining compact and well separated clusters, as well as cluster free areas would confirm the existence of preferential locations.

Clustering has been a well studied problem over the past decades, in a number of fields. Early popular clustering algorithms refer to hierarchical methods [Johnson 1967]. Hierarchical or agglomerative clustering initially considers all samples as clusters. The clusters are then progressively obtained by merging neighboring clusters together. The process continues until there is only one cluster left. This results in a dendogram, a tree that represents the progressive clusters merging. The clustering is chosen by selecting a number of clusters and cutting the grouping tree accordingly. There are two issues with such a method, the first refers to the absence of "cluster center" representing each population while the second to the choise of the number of clusters being retained. More adapted to our problem, Centroid based clustering methods group the data by selecting cluster centers and assign the nodes to their closest center. Each cluster center is then representative of the cluster's properties. The most popular centroid based clustering method is the K-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]].

K-means clustering is an instance of the EM algorithm where the objective is to minimize the square distance between the samples and the cluster center:

J = K k=1 p∈V LGG p -ck 2 (3.2)
Given the expected number of clusters K, and an initialization of the centers { c1 , .., c K }, the Kmeans algorithm alternates between assigning the samples to the closest center and reevaluating the centers as the mean coordinates of all samples within the cluster. The K-means algorithm is sensitive to initialization which can be a drawback when there is no pre assumption on the position of the clusters. Furthermore, the cluster centers are constructed as the mean value of the samples in the clusters, which can be sensitive to the presence of outliers and not represent the cluster appropriately since there is no guarantee that the mean lives on the same manifold as the one of the samples. These issues are often addressed by considering different initializations both in terms of the number of populations as well as their centers (this is done by randomly sampling the population) and then adopting criteria from information theory to determine the optimal result.

We adopt a more robust method [Komodakis 2008a] that has been shown to outperform the Kmeans algorithm. The method simultaneously estimates (i) the number of clusters K (ii) the cluster centers { c1 , .., c K } ∈ V LGG , and (iii) the remaining nodes' assignments to the different clusters A p ∈ {C 1 , ..., C K }. The objective function aims at minimizing the distance between nodes and their assigned cluster center:

min K min ck min Ap K k=1 p∈V LGG δ(C k , A p )d(p, ck ) + α K k=1 V pen ( ck ) (3.3)
where δ(.) is the Kronecker delta function and α is a constant parameter that balances the importance of the two terms. The first term of the objective function assigns nodes to the cluster of which center they are the closest according to a given distance on the manifold of the samples' space. The second term is introduced in order to avoid the trivial solution of defining each node as a cluster center. This can be seen as a regularization term that measures the stability of a node in order to determine which nodes should be selected as centers. In the context of seeking preferential locations, we aim at selecting nodes that are at the center of high density areas. We define the penalty term as the mean distance over all samples in the graph in order to identify such nodes:

V pen (p) = 1 N q∈V LGG d(p, q) (3.4)
The idea behind this term is that nodes should be selected as center if they are close to many other nodes, i.e. if its distance value with an important set of nodes is low. Contrarily to the K-means algorithm, it should be noted that the centers are among the samples of the population on which clustering is to be performed and that they are chosen as the ones situated on dense (statistical interpretation) manifold areas. This makes the center selection more robust to outliers and more faithful to the actual cluster structure. The additional advantages of this clustering approach is that it is modular with respect to the distance metric between samples while the optimization also recovers the optimal number of clusters.

Similarly to the Fast-PD algorithm for MRF optimization, the near optimal clustering is recovered using linear programming by minimizing the gap between the primal and dual problems. The problem is reformulated as the following relaxed linear integer program:

PRIMAL-IP ≡ min x i V pen (i)x i + j,i d(i, j)x ji (3.5) s.t. i x ji = 1 (3.6) x ji ≤ x i (3.7) x ji , x i ∈ {0, 1} (3.8)
The binary variable x ji indicates whether the sample j has been assigned to node i or not, while the binary variable x i indicates if node i has been selected as a cluster center. The constraints of equation 3.6 simply express the fact that each observation should be assigned to one central node, while the constraints of equation 3.7 impose that if observation j has been assigned to node i then node i must obviously be chosen as a cluster center. The most important issue for such an integer program optimization is the correct assignment of the variables x i , i.e., deciding which nodes will be chosen as the clusters centers. The center selection process is guided by the nodes' stability.

Once the central nodes are computed then each observation can be trivially assigned to its nearest central node. This process repeats until no more stable nodes can be found. The resulting method does not require any initialization and provides near-optimal clustering results both in terms of number of clusters as well as in terms of central nodes and remaining nodes assignments.

Clustering Selection

The parameter α weighting the cost of introducing new clusters versus the overall assignment distance error plays a critical role in the obtained clustering. We seek to find the clustering that is the most representative of the data. To this end, we rely on cluster validity indices that provide a quantitative measure of how good the obtained result is. Two kinds of validity criteria exist: external and internal indices. The external indices imply prior knowledge about the data and the expected clustering. They are used to compare the obtained clustering structure with the a priori knowledge. In practice, external criteria are mostly used to evaluate a clustering method as the partition is rarely known and internal criteria are the most useful. In this work, we make use of three widely used internal cluster validity measures: the Dunn index [Dunn 1974], the Davies-Bouldin index [Davies 1979] and the Silhouette index [Rousseeuw 1987].

The Dunn index is among the simplest clustering quality measures. Its goal is to determine if the clusters are compact and well separated. The index, introduced in [Dunn 1974] is defined as the ratio between the minimum distance between clusters d C i ,C j and the maximum cluster diameter diam(C k ):

D = min i=1,...,K min j=1,...,K d C i ,C j max k=1,...,K diam(C k ) (3.9)
In [Dunn 1974], the distance inter cluster is defined as the minimum distance between any samples in the two different clusters while the cluster diameter is the maximum distance between two samples of the same cluster:

d C i ,C j = min p∈C i ,q∈C j d(p, q) (3.10) diam(C k ) = max p,q∈C k d(p, q) (3.11)
Compact and well separated clusters yield a large Dunn index (large distance inter cluster and small cluster diameter). The Dunn index suffers from several shortcomings, the most important being its severe sensitivity to noise and outliers as they can strongly impact the diameter and inter cluster distance. In order to obtain more robust measurements, several Dunn-like indices can be considered with different definitions of the diameter and cluster distance.

A graph theoretic formulation was notably proposed in [Pal 1997]. Each cluster was represented using graphs of known structure (Minimum Spanning Tree (MST), Gabriel Graph or Relative Neighborhood Graph) and the cluster diameter was defined as the maximum edge value in the graph. We will review the MST based Dunn index that is used in this chapter. Given the set of nodes V C k that belong to cluster C k and the corresponding complete graph G C K that connects all pairs of nodes, a spanning tree is the minimal subgraph that connects all nodes together. The weight of the spanning tree is defined as the sum of the weights (i.e. distances between nodes) of all the edges in the tree. The Minimal Spanning Tree G

M ST C K = (V C k , E M ST C k
) is the weighted graph that has the smallest weight. Having constructed such a graph per cluster, the diameter is then defined as the edge of the MST that has the maximum weight diam(

C k ) = max p,q∈E M ST C k d(p, q).
For increased robustness, the inter cluster distance is modified and defined as the distance between the clusters centers d C i ,C j = d( ci , cj ).

Similarly to the Dunn index, the Davies-Bouldin index aims at identifying compact and well separated clusters. A similarity measure R ij between the clusters C i and C j is computed based on a cluster dispersion measure S i and compared to an inter cluster dissimilarity measure D ij . The similarity measure should verify the following conditions:

R ij > 0 R ij = R ji if S i = 0 andS j = 0 then R ij = 0 if S j > S k andD ij = D ik then R ij > R ik if S j = S k andD ij < D ik then R ij > R ik (3.12)
x Usually, the similarity, dispersion and dissimilarity measures are defined as follows:

S i = 1 |C i | p∈C i d(p, ci ) D ij = d( ci , cj ) R ij = S i + S j D ij (3.13)
The dissimilarity measure corresponds therefore to the ratio of the mean intra-cluster distance (distance of all nodes to their cluster center) for clusters C i and C j with the inter-cluster distance (distance between cluster centers). The Davies-Bouldin index is then calculated as the average maximum similarity measure over all clusters:

DB = 1 K K j=1 max i=1,...,K,i =j R ij (3.14)
The optimal clustering corresponds to a low value of the DB index, i.e. to a low similarity between the clusters. This index is more robust than the Dunn index due to the fact that it relies on average values and distances to the center of the cluster. The last index considered, the Silhouette index, takes a different approach. Computed for each sample p, it evaluates how adequate the sample's cluster assignation is via a score of confidence s(p):

s(p) = b(p) -a(p) max(a(p), b(p)) (3.15)
a(p) is the average distance between sample p and all the remaining elements assigned to the same cluster (dissimilarity to all other samples in the same cluster). For each cluster C k the node hasn't been assigned to, the same dissimilarity measure d p,C k is computed between node p and all the nodes in the cluster. We define b(p) = min k=1,...,K d p,C k , which can be seen as the average distance between p and all the elements in the closest cluster. The Silhouette index takes values between -1 and 1. If the value is close to 1, the node has been assigned correctly. A value close to zero suggests that the node is equally far away from 2 clusters, while a silhouette index close to -1 infers that the node has been misclassified. To evaluate the quality of the clustering, we compute the global silhouette index:

GS = 1 K K j=1 1 n j p∈C j s(p) (3.16)
where n j is the number of samples in cluster C j . The highest global silhouette index corresponds to the best clustering, where the individual Silhouette indexes are closest to 1. The silhouette index gets its name from the fact that it offers a visualization tool of the quality of the clustering. The so-called silhouette of the cluster can be plotted, displaying which samples are well classified and the relative quality of clusters, narrow silhouettes indicating better clusters.

In order to select the best fitting clustering, the three different indices are computed for several values of α and the parameter that yields the maximum combined index is selected. In the remainder of this chapter, we will consider the optimal clustered graph made of K clusters.

Tumor Characterization and Segmentation with Spatial Position Priors

Given the possibility to recover properties of the tumors with respect to their localization, it is of high interest to integrate a new, non segmented patient to the atlas and determine which preferential location the tumor belongs to. In parallel, the atlas offers powerful insight on the expected positions and expansion of the tumors and can be used as a spatial prior for segmentation. This section presents an application of the DLGG atlas for characterization of the tumor's location with respect to the atlas and for segmentation via the use of spatial position priors. For each cluster C i of the atlas, we build a probability map P (x|C i ), i ∈ [1 : K] that describes the rate of tumor appearance per voxel x. This is carried out in a similar fashion as the construction of the probabilistic ICBM atlases by averaging the binary segmentation maps of the different tumors in the same cluster. In order to ensure smoothness of the probability maps for clusters with a small number of samples, all probability maps are smoothed by Gaussian filtering. Fig. [3.3] shows examples of probability maps for different clusters, superimposed to the average volume.

We consider an image I defined on a domain Ω featuring a DLGG. An initial tumor detection, performed via statistical classification, is used within a coupled MRF framework for simultaneous characterization and segmentation of the tumor. The approach is guided by the appropriate cluster's prior probability map. The method differs from existing atlas-based tumor segmentation methods since our atlas describes the expected position of the tumor (instead of detecting the tumor with respect to the healthy structures [START_REF] Gering | [END_REF], Prastawa 2004, Menze 2010]) and independent from the observed image (contrarily to the work of [Prastawa 2003] that modifies a healthy atlas based on the observation).

Tumor Detection

The first step towards tumor segmentation and characterization is an initial detection of the tumor to roughly estimate its position. We rely on boosting due to the simplicity of the Adaboost algorithm [START_REF] Freund | [END_REF]], and its computational efficiency once the classifier is learned.

Adaboost Algorithm

Consider an ensemble of labeled volumes X featuring a brain tumor. N voxels are extracted from the different volumes to construct a training set that consists of pairs {x i , y i }, i ∈ [1 : N ] where x i ∈ X is an observed image voxel and y i is its corresponding label (-1 for background, 1 for tumor). Each pair is associated to a feature vector f (x i ) and a weight D i . The idea of boosting is to construct a strong classifier H as a linear combination weakly discriminative classifiers (the weak learners). It is an iterative process, where a new weak learner is constructed at each iteration. The weights D i determine the importance of their corresponding sample assignment during classification. In each step of the process, a each weak learner is built, and the samples' weights are updated so that the focus is on the misclassified and more difficult samples at the next iteration.

Initially, all weights are set uniformly: D 0 i = 1/N while the user specified the number of iterations T and the kind of weak learner used. At each iteration t, a new weak learner h t : X → {-1, 1} is constructed based on a value in the feature vector in order to minimize the classification error

ǫ t = N i=1 D t i |h t (x i ) -y i | (3.17)
Once the classifier has been chosen, its importance with respect to the strong classifier is computed via the weighting parameter α t and the weights are updated in order to enhance the focus on misclassified samples:

α t = 1 N log( 1 -ǫ t ǫ t ) (3.18) D t i = 1 Z i D t-1 i exp(-α t y i h t (x i )) (3.19)
where Z i is a normalizing constant. Eventually, the strong classifier is computed as:

H(x) = T t=1 α t h t (x) (3.20)
Among nice properties of the AdaBoost algorithm, it was shown in [START_REF] Freund | [END_REF]] that the classification error rate decreases exponentially fast while [Friedman 2000] showed that the boosting score could be converted to probabilities as:

P (y = 1|x) = e H(x) e -H(x) + e H(x) P (y = -1|x) = e -H(x) e -H(x) + e H(x) (3.21)
In [Friedman 2000], a variant of the AdaBoost algorithm called Gentle Adaboost was introduced. The idea was to provide a method that is less sensitive to outliers and noise. Indeed, the use of log-ratios in the evaluation of α t can be numerically unstable and lead to an excessive focus on outliers. The Gentle Adaboost algorithm has been shown to yield more robust results and proceeds as follows:

1. Initialize weights D i = 1/N, i ∈ [1 : N ] and classifier H(x) = 0 2. Repeat for t = 1, ..., T (a) Find a weak classifier that minimizes

ǫ t = N i=1 D i (y i -h t (x i )) 2 (b) Update H(x) ← H(x) + h t (x) (c) Update D i ← 1 Z i D i exp(-y i h t (x i ))
3. Output classification result: sign

T t=1 h t (x)
In this thesis, the Gentle Adaboost algorithm is used while the weak learners are defined as decision stumps.

Feature Vector

One of the key elements for a good separation between the tumor and background classes is the design of the feature vector that has to be discriminative. We consider iconic, texture and symmetry information.

• Iconic features rely on the intensity values of the image and the fact that the tumor has a different appearance than the background. In order to integrate neighborhood information, 9 × 9 × 5 intensity patches centered on the voxel of interest are added to the feature vector.

The patches are constructed by selecting one voxel out of two in the image for increased neighborhood exploration. Rotation invariant statistics (median, standard deviation and entropy) on patches of size k × k × 3 where k = [3, 5, 7] are also included. The smaller patch size in the z axis is due to the lower resolution in this direction.

• Texture information is added via the computation of Gabor filters [Manjunath 1996]. Gabor filters have been used in various medical image analysis problems in order to detect the texture patterns in the image. A 2-D Gabor function is a sinusoidal plane wave of a certain frequency λ and orientation θ (the carrier) modulated by a Gaussian envelope of scale σ = (σ x , σ y ):

G(x, y) = 1 2πσ x σ y exp - 1 2 x ′2 σ 2 x + y ′2 σ 2 y + i 2π λ (x ′ + y ′ )
where x ′ = x cos θ and y ′ = y cos θ

(3.22)
Convolution of the Gabor function to an image yields a high response when the local image structure is similar to the Gabor function based on its scale, wavelength and orientation. A set of Gabor functions G σ,θn (x, y) (filter bank) can be constructed by selecting a number of scales S and orientations K so that it thoroughly paves the frequency space. In practice, 3D computation of the Gabor filters is very slow. We adopt the method proposed in [Zhan 2003] that approximates the 3D filter by computing two orthogonal 2D filter banks. We use the fast rotation and scale invariant implementation used in [START_REF] Michel | [END_REF]].

It was observed that texture features have little discriminative power with respect to the tumor. However, they are very efficient towards determining the brain's main structures and consequently reducing the risk of false detections. The filters are computed on 2 scales and 10 orientations. The image response to each obtained filter at voxel x constitutes one element of the feature vector f (x).

• The last feature used is specific to the brain and exploits the approximative symmetry between the two hemispheres. The idea is that the presence of the tumor will introduce a notion of dissymmetry. Assuming that an approximate symmetry plane Π is known, the symmetry measure evaluates the difference between the intensity of a voxel x in the left part of the image and its symmetric x Π . To compensate the approximation in the selection of the symmetry plane, the mean intensity of a patch around x and x Π is computed instead of using the single intensity value. The symmetry feature is thus computed as follows:

S(x) = 1 N y∈P(x) I(y) - 1 N y Π ∈P(x Π ) I(y Π ) S(x Π ) = -S(x)
(3.23)

where P(x) corresponds to the patch centered around x and N is the number of voxels in P(x). Given that the volumes to be segmented are registered to the same reference space where the atlas was built, estimating an approximate symmetry plane is straighforward.

Visual examples of the different elements constituting the feature vector can be seen in Fig.

Coupled Markov Random Field Model

We will now proceed with the presentation of the Markov Random Field model that couples tumor characterization and segmentation with spatial priors. We formulate the problem of tumor segmentation and characterization as a labeling problem where we want to assign a global volume wise characterization label c(x) and a voxel wise segmentation label s(x) to each image voxel.

The characterization label corresponds to assigning the patient/tumor to one of the atlas' K clusters c(x) ∈ C 1 , ..., C K , i.e. determining which preferential location the tumor is situated in. Given the assumption that the patient has only one tumor, the same label has to be assigned to all the voxels. The segmentation label s(x) ∈ {0, 1} simply characterizes the voxel as being part of the tumor (s(x) = 1) or the background (s(x) = 0). We adopt a pairwise MRF model on a graph G = (V, E) where each node x ∈ V corresponds to an image voxel while the edges constitute the neighborhood N (x) by connecting node x to its 6 immediate neighbors in 3D. A node x is to be assigned a label l x ∈ L = {1, • • • , K} × {0, 1} that jointly characterizes and segments the tumor. We map this label space into one dimension as follows:

c(x) = c lx = C k where k = l x (mod K) s(x) = s lx = l x K (3.24)
The optimal labeling l = {l x |x ∈ V} is recovered by minimizing the MRF energy defined as:

E(l) = x∈V V x (l x ) + x∈V y∈N (x)
V xy (l x , l y ) (3.25)

Pairwise Term

The pairwise term acts as a smoothing prior. Its role is two-fold: it imposes local spatial consistency of the voxel wise segmentation and ensures that the same characterization label is assigned to the whole volume. It is designed following the Potts model, by penalizing neighboring nodes that have been assigned different labels:

V xy (l x , l y ) =      0, if l x = l y β, if s lx = s ly and c lx = c ly γ, otherwise (3.26)
Where β is a constant parameter that describes the strength of the penalization imposed on the segmentation field and β ≪ γ.

Unary Terms

We define the MRF unary term as a sum of three potentials :

V x (l x ) = V S (s lx ) + γV C (c lx ) + µV S,C (s lx , c lx ) (3.27)
The first two terms V S and V C are observation dependent and focus on the segmentation and characterization space respectively, the third term V S,C couples the two, while γ and µ are constant parameters that balance the relative importance of the different potentials.

For the segmentation term V S (s lx ), we make use of the classification likelihoods P tm (I(x)) and P bg (I(x)) recovered through Gentle AdaBoost learning. This terms seeks the label configuration that yields the maximum likelihood probability. Inspired by the MAP-MRF formulation and following the approach of [Boykov 2006], we use the log probabilities and define the term as follows:

V S (s lx ) = -s lx log P tm (I(x)) -(1s lx ) log P bg (I(x))

(3.28)

The characterization term's role is to identify which cluster the tumor to be segmented belongs to. Even though the precise position of the tumor is unknown, the boosting classification score H(x) gives a rough estimate of the position of the tumor pixels. The closest cluster is identify by measuring the amount of non overlap between the binary map associated to the center of each cluster and the voxels with a high boosting classification score.

V C (c lx ) = max(0, H T (I(x)) -Mc lx (x)) (3.29)
Where MC i (.) is the binary segmentation map of the tumor that is the center of the cluster C i , and H T (.) is a binary map representing voxels with a high classification score, obtained by thresholding (with threshold T) of the boosting classifier output:

H T (I(x)) 1, if H(I(x)) > T 0, otherwise (3.30)
This potential counts the number of voxels that are classified as tumor by boosting and do not overlap with the tumor volume of the cluster center. The closest cluster corresponds to the one that has the highest overlap with the most probable tumor voxels.

Last but not least, the coupling term V S,C (s lx , c lx ) links the segmentation and characterization labels via the DLGG atlas and plays an essential role in the segmentation process. Let us recall that each cluster C i is associated to a probability map P (x|C i ) that describes the spatial distribution of tumor voxels in the cluster. Given a characterization label (cluster assignment), the corresponding probability map gives insight on where the tumor is expected to be and not to be. This is used as a spatial prior that will encourage or discourage segmentation decisions based on the voxel's localization:

V S,C (s lx , c lx ) = -s lx log P (x|c lx ) -(1 -s lx ) log 1 -P (x|c lx ) (3.31)
Spatial consistency is enforced between the segmentation label and the expected spatial distribution of the tumors of the corresponding cluster.

MRF Optimization Using Linear Programming

The MRF energy is minimized using Fast-PD [Komodakis 2007b, Komodakis 2008b]. The algorithm relies on linear programming and primal-dual strategies. Considering a pair of primal and dual linear programs:

PRIMAL : min c T x DUAL : max b T y s.t. Ax = b, x ≥ 0 s.t. A T y ≤ c (3.32)
It is possible to obtain an f-approximation of the solution to the Primal problem by minimizing the cost gap between the primal and dual problems so that their ratio is inferior to f . This is carried out by iteratively generating pairs of primal-dual solutions {x t , y t } until the elements of the last generated pair are both feasible and have a primal dual gap inferior to f :

c T x t ≤ f.b T y t (3.33)
This strategy is applied to the MRF energy by reformulating the problem as the following integer program:

min p∈V l∈L V p (l)x p (l) + p,q∈E l,l ′ ∈L V pq (l, l ′ )x pq (l, l ′ ) (3.34) s.t. l x p (l) = 1 ∀p ∈ V (3.35) l x pq (l, l ′ ) = x q (l ′ ) ∀l ′ ∈ L, (p, q) ∈ E (3.36) l ′ x pq (l, l ′ ) = x p (l) ∀l ∈ L, (p, q) ∈ E (3.37) x pq , x p ∈ {0, 1} (3.38)
The binary variables x p (l) and x pq (l, l ′ ) characterize the label assignments. x p (l) = 1 implies that label l has been assigned to node p while x pq (l, l ′ ) = 1 indicates that nodes p and q are assigned label l and l ′ respectively. Constraint 3.35 imposes that the node is only assigned one label and constraints 3.36 and 3.37 impose consistency between the unary and pairwise binary variables.

The linear programming relaxation of the integer program consists in relaxing constraint 3.38 to the constraints x p (.) ≥ 0 and x pq (., .) ≥ 0. The primal-dual scheme is then applied to the relaxed integer program and its corresponding Dual problem:

max p∈V y p (3.39) s.t. y p ≤ min l∈L h p (l) ∀p ∈ V (3.40) y pq (l) + y qp (l ′ ) ≤ V pq (l, l ′ ∀(l, l ′ ) ∈ L, (p, q) ∈ E (3.41) (3.42)
where h p (.) = V p (.) + q,(p,q)∈E y p,q (.) is called the height variable. For MRF optimization, each iteration (generation of primal-dual pairs) corresponds to a max-flow optimization of a capacitated graph constructed based on the current pair of primal-dual variables.

The method has several advantages among which its generality (the pairwise potentials only have to be positive), optimality guarantee (f-approximation of the optimum) and computational efficiency.

Experimental Validation

Data-set and Preprocessing

The database used to build the statistical atlas consists of 210 carefully selected patients. All of them were over 18 years old and treated in the same institution. A histological diagnosis has been performed and FLAIR volumes as well as molecular analysis where available prior to any treatment. The tumor size was limited to 180cm 3 so that the center of the tumor remains representative of where the tumor initially appeared: the tumor expansion would shift its center of gravity. The distribution of tumors' sizes is shown in Fig. [3.6]. The database consisted of 97 females and 113 males of mean age 38 years old and mean tumor volume 45 cm 3 . The tumor size ranged from 0.1 to 180 cm 3 .

The patients' MRIs considered for the construction of the atlas were 3D FLAIR volumes that were all manually segmented. The image size varied from 144 to 576 in the (x,y) plane et 12 to 213 in the z plane, and the voxel resolution from 0.4x0.4 to 1x1 mm 2 in the (x,y) plane and 0.9 to 10 mm in the z plane.

The reference pose used for registration is a 256x256x24 FLAIR image of a healthy brain, with resolution 0.9x0.9x5.5. Because of the strong slice thickness, the segmentation maps were no longer binary after registration due to interpolation tasks. The same threshold of 0.2 determined experimentally as the most faithful segmentation of the registered tumor is used on all the data and in the remainder of this thesis to convert the registered maps to binary maps.

The boosting features rely heavily on the images intensities that need to be comparable between volumes. The intervolume intensity inhomogeneity inherent to the MRI modality is dealt with using a simple method, in order to maintain the contrast between the tumor and the background. All images are set to a similar intensity range by setting them to the same median and interquartile range:

I reg = I -med(I) IQR(I) (3.43)
where med(I) and IQR(I) are respectively the median and interquartile range of image I, computed without taking into account the background area (lowest intensity values). The median was used instead of the mean value due to the presence of the tumor that can impact the mean intensity value, especially if the tumor's size is significant. No correction was done on the image's bias field.

Sparse Graph Validation and Robustness

In order to select the best parameter α for clustering, we clustered the complete graph for different values of α ranging from 0.3 to 7 which yields 5 to 49 clusters. For each experiment, we computed the value of the three indices (Dunn, Davies-Bouldin and Silhouette). All indices showed a peak corresponding to 11 clusters. The indices have been normalized and combined using their product:

I c = (D.
(1 -DB).GS) 1/3 with the underlying idea that if an index has a low value, the combined value should also be low. The combined index value showed a strong maximum for α = 2.1 which corresponds to a 11 clusters graph. The different indices values as well as their combined value with respect to parameter α and the number of clusters can be seen in Next, we evaluated the robustness of the graph by performing Leave One Out Cross Validation. We carried out 210 experiments where one sample p was removed from the data set. For each experiment, we learned the graph topology on the 209 remaining samples by selecting the clustering that yielded the maximum value of the three indices. An 11 clusters graph was obtained for 207 experiments (98.5 % of the data). A 10 clusters graph was obtained for 3 samples. Those samples were the center or close to the center of the same loose cluster that was merged with its closest neighbor when the sample was removed.

We then studied the robustness of the graph by analyzing the influence of each sample. For 1. Distance from the center of the cluster. The sample has to be assigned to its closest cluster and its distance to the center should follow the distance distribution within the cluster. We determine if the sample is an outlier based on the distribution's quartiles :

d(p, ck ) = min j∈[1:K] d(p, cj ) and d(p, ck ) ≤ Q 3 (C k ) + 1.5(IQR(C k )) (3.44)
where Q 3 (C k ) and IQR(C k ) are respectively the third quartile and interquartile range of the distance distribution in cluster C k .

2. Individual Silhouette index. The global silhouette index is indicative of the quality of the clustering. The individual index tells how well the sample fits in its cluster. For each cluster k, we compute the corresponding silhouette index. If the index is close to 1, there is no doubt about the fact that the sample belongs to this cluster. A value close to zero does not necessarily mean that the clustering is bad. Indeed, the main drawback of the silhouette is that it can penalize loose clusters (large diameter) when compared to smaller ones. The sample should be assigned to a cluster if it maximizes the silhouette index and fits in the silhouette of the cluster. Similarly to the first criterion, we make use of the quartile values:

s k (p) = max j∈[1:K] s j (p) and s k (p) ≥ Q 1 (C k ) -1.5(IQR(C k )) (3.45)
where Q 1 (C k ) and IQR(C k ) are respectively the first quartile and interquartile range of the silhouette index distribution in cluster C k and s k (p) is the silhouette index of sample p with respect to cluster C k .

3. The last criterion verifies that the same cluster is selected by criteria 1 and 2: Using those criteria, we successfully classify 88% of the whole data-set. Last, we evaluated topology correspondences between the 207 eleven-components clusters graphs and the graph obtained from the whole database. To this end, we use the graph matching algorithm presented in [Torresani 2008]. Complete match of the graphs was obtained in 70% of the cases and only one node differed in 26% of the cases. The rest of the samples had 2 (9 samples) or 4 (1 sample) different nodes. For the cases where the cluster centers were different, they were still situated in the same area in the brain and the centers were close. The maximum distance between two different centers is 12.6 which corresponds to the mean intra cluster distance of tight clusters. Complete and partial matching are illustrated in Fig.

Segmentation Results

Segmentation was performed on the 210 volumes via Leave One Out Cross Validation. The atlas was constructed on the 209 remaining samples according to the optimal value of α. For boosting learning, 210 classifiers were learned by randomly selecting 35 volumes among the 209 samples used for graph construction. For the symmetry feature, we rely on the fact that all images are registered and use an approximate symmetry plane estimated on the reference pose. Estimating an exact symmetry plane is out of the scope of this work and has been the subject of elaborated methods [START_REF] Sylvain Prima | Computation of the midsagittal plane in 3-D brain images[END_REF].

We evaluated the segmentation quality with respect to the manual segmentations by computing the Dice Score (D), False (FP) and True positive (TP) rates and Mean Absolute Distance between contours (MAD). Given a manual segmentation M and an automatic segmentation A, the different scores are defined as follows:

D = 2 A ∩ M A ∪ M + A ∩ M T P = A ∩ M M F P = A -A ∩ M A M AD = 1 2 1 n A x∈A d min (x, M ) + 1 n M x∈M d min (A, x) (3.47)
where n A and n M are the number of voxels in A and M respectively and d min (x, M ) = inf(d(x, y)|y ∈ M ) is the smallest distance of a point x to any point in M . The distance d(., .) used here is the euclidean distance. Basically, the MAD score characterizes how close the two contours are.

The method was compared to sole boosting classification and pairwise MAP-MRF segmentation using the boosting probabilities as likelihood and an Ising model as pairwise cost. The parameter β determining the amount of smoothing of the segmentation field was set to 1 with prior and 3 without prior (MAP-MRF). Parameters γ and µ were set to 1 and 10000 respectively. The high value of µ was necessary for small tumors that caused very little variation in the value of the potential V C (.). The threshold parameter T driving the construction of the binary boosting score is set to 1.5. The parameters µ and T were determined so that the maximum number of images were assigned to their right cluster (inferred from the complete graph construction).

Using our framework, we observed an increase in the quality of the segmentation and a significant drop of false detections. The different mean and median scores for the three methods are shown in table 3.2 and boxplots of the scores over all images are shown in Fig. [3.12]. We observe an increase of the Dice score (+ 3% mean and median) and strong reduction of the false positive rate (-9-10 % mean and median). The increase in quality is highlighted by the MAD score, as illus- [3.15] also highlight the impact of the spatial prior, where erroneous boosting detections are corrected. A reduction of the true positive rate is observed (-4 %), mainly due to mischaracterized tumors and outliers as well as the limited number of samples in certain clusters. This phenomenon could be reduced by an increase in the size of the data set and detecting outlier tumors (as tumors equally distant from two cluster for instance) that should not be assigned to any cluster. Wrong clusters assignments are directly connected to the quality of the detection step and principally concern small tumors. Increasing the characterization quality would mostly depend on an increase of the quality of the detection. The discrepancy between the mean and median Dice values (69% and 77%) was due to those small tumors, that were poorly detected by the classifier and therefore poorly segmented.

Discussion

In this chapter, the construction of a probabilistic atlas mapping the preferential locations of DLGGs has been presented. The atlas was constructed through a complete graph where each node represents a tumor and the edges represent the proximity between the tumors, evaluated as the Euclidean distance between their center of gravity. The DLGGs' preferential locations were then identified by regrouping nearby tumors in clusters where each cluster corresponds to a preferential location. The optimal clustering was determined using well-established cluster validity indices. To each cluster is associated a probability map representing the spatial extension of tumors within the cluster.

A natural application of this tumor specific atlas was in the context of tumor segmentation as a spatial prior. Such an application was presented in the second part of this chapter. The atlas was integrated in a coupled MRF framework where the tumor is concurrently characterized (assigned to one of the atlas' clusters) and segmented. An initial statistical classification based detection enables assignment to the cluster as well as segmentation with spatial smoothing (similarly to MAP-MRF approaches). The two problems (characterization and segmentation) are coupled via a spatial prior term dependent on the characterization and segmentation labels. This term makes use of the cluster's probability maps to enforce spatial consistency with the segmentation and the expected localization of tumor pixels within a given cluster.

Aside from the obvious appeal of tumor segmentation, the atlas enables the study of the tumor behavior with respect to its localization which could be correlated with the patient's characteristics (age, sex...) or the tumor's characteristics (histological type, molecular biology...). This could potentially unravel behaviors (growth rate, anaplastic transformation) dependent on specific parameters and enable to adapt treatments to the patient. Going further, the interaction between the tumors and their surrounding tissue and functional areas as well as their growth patterns could be studied with respect to the cluster they're in. This could lead to models of tumor growth and their impact on brain functional reorganization.

The latter applications are however limited by the use of affine registration. Consequently, the interaction between the tumors and surrounding structures would strongly lack precision. A deformable registration could handle these issues, where the presence of the tumor has to be taken into account in order to obtain realistic transformations. Integration of anatomical information would also increase the quality of segmentation results that, in our approach, is strongly dependent on the adequacy of the selected cluster. An outlier tumor or a poorly detected tumor can be assigned to a cluster that does not correspond to its location. Anatomical information would reduce the impact of bad assignments as it offers a second means of integrating spatial prior knowledge.

The development of such a deformable registration method, coupled with tumor segmentation, will be the subject of the next chapter. Chapter 4

Joint Tumor Segmentation and Dense Deformable Registration

Introduction

In the previous chapter, a DLGG specific atlas was constructed in order to study the tumor appearances in the brain and their behavior. The methodology was designed to be robust to errors in registration, relying on the tumors center of gravity instead of surfaces. In general, the presence of tumors violates common registration procedures and therefore, a global affine registration was used to map all images in the same reference coordinate system. Interesting and promising results where obtained with respect to the mapping of DLGGs in the brain and their preferential locations.

Nevertheless, an improvement on the precision is expected to be obtained via the use of deformable registration that would also provide the possibility of associating the presence of the tumor with the surrounding structures and in particular measure the impact of tumors on them. Generally speaking, when constructing pathology related atlases, the registration method should take into account the presence of the pathology and the absence of correspondence in the atlas.

In parallel to atlas construction, the task of atlas based segmentation of tumor bearing images has gained a lot of interest. A healthy reference atlas is registered to a pathological subject in order to segment its healthy structures and/or the tumor. For tumor segmentation (see chapter 2), a simple registration scheme is usually proposed (affine registration [Prastawa 2004], or segmentation based registration with the tumor labeled healthy [Kaus 2001]) and subsequently, a segmentation method is designed in order to cope with registration errors. However, such methods fail when the registration error is significant (i.e. when the brain is deformed by the tumor). If the healthy structures are to be segmented based on the atlas' structures segmentation, affine registration cannot match the different anatomies. Those issues highlight the need for a non-rigid registration scheme for tumor bearing images in order to obtain robust and precise segmentations.

Most of the existing methods for registration with missing correspondences consider that the tumor has been determined by a prior segmentation. The coupled segmentation/registration task is thus treated sequentially. Such methods do not exploit the inherent link between segmentation and registration, where the segmentation benefits from a good registration (atlas anatomical information) and the registration benefits from the segmentation by taking into account the area without correspondences. This chapter presents a concurrent registration and segmentation framework where the segmentation map and deformation field are estimated simultaneously leading to better solutions for both problems. A Markov Random Field model is adopted, seeking to assign to each graph node a label that corresponds to a displacement and a segmentation label. Statistical classification is used to provide tumor vs background prior likelihoods along with conventional registration similarity metrics. In the missing correspondences area (tumor), the registration term is relaxed and the area is implicitly deformed by interpolation. The remainder of this chapter is organized as follows: section 4.2 present the related work on brain registration in the presence of a pathology and concurrent registration and segmentation. Section 4.3 introduces the joint MRF framework while section 4.4 present a stochastic approach to account for the main drawback of discrete MRF based registration, which is the trade off between complexity and sampling precision. Experimental validation will be part of section 4.4 and a discussion concludes this chapter. A short graphical overview of the method is shown in Fig.

Related work

The following notations will be used in this section: I is the fixed or target image and J is the moving or source image defined on a domain Ω. T is the transformation that maps the source image to the target image.

Registration with Missing Correspondences

For the task of non rigid registration in the presence of brain tumors, we can divide prior art into two groups of methods. The first class of methods masks the pathology and does not take it into account during the estimation of the matching criterion. The deformation in and around the tumor is carried out by interpolation via a regularization term. The second group simulate a tumor in the healthy atlas to create correspondences. This is typically designed as an atlas seeding process where the origin of the tumor in the atlas (the seed) is determined, followed by the growth of the artificial tumor carried out using growth models.

Pathology Masking

A simple approach for pathology masking was proposed in [Brett 2001] and tested on a variety of brain lesions (infarction, metastasis and aneurysm). They used the registration method proposed in [Friston 1995] (distributed in the Statistical Parametric Mapping (SPM) package). An affine registration is followed by a parametric non linear registration. The transformation model is based on a discrete cosine transform (the basis functions are cosine function) and while the Sum of Square Differences was used as matching criterion. The lesion is manually segmented and its corresponding area is masked so that it is not taken into account for the evaluation of the optimal displacement.

In the work of Stefanescu et al. [Stefanescu 2004], the registration problem was designed via two regularization processes. The following objective function is to be minimized:

E = SSD(I, J • T ) + β Ω k(x) ∇ ∂T ∂t 2 + γ Ω D(x) ∇T (x) 2 dx (4.1)
The first term is the similarity criterion, defined as the Sum of Square Differences (SSD). The second term is a fluid regularization of the deformation field that allows large local displacements, while the third term is an elastic regularization which role is to preserve anatomical coherence. The space varying weighting parameter D(x) varies with respect to the anatomical tissues deformability, simulating the behavior of a biomechanical model. The parameter that accounts for the presence of the tumor is the fluid regularization weighting parameter k(x) ∈ [0, 1]. This parameter is a confidence score on the reliability of the matching criterion. If k(x) is close to zero, the smoothness of the deformation is heavily imposed and almost interpolated from neighboring values. In healthy structures, the value of k depends on the image gradient. It is set to zero in the tumor area, and therefore the deformation field is interpolated. The optimal displacement is recovered through a four steps iterative process. First the velocity field v(x) is recovered by gradient descent of the SSD criterion (which amounts to approximating the demons algorithm forces [START_REF] Pennec | [END_REF]). Second, a weighted fluid regularization of this field is performed as ∂v ∂t (x) = k(x)∆v(x). Eventually, the deformation field is updated T t+1 = T t • v and smoothed by elastic regularization ∂T ∂t (x) = ∇(D(x)∇T (x)). An extension of the method was presented in [Commowick 2005] to construct an anisotropic weighting field D(x) based on image statistics.

The method was more recently adapted to the diffeomorphic demons framework [Vercauteren 2008] in [START_REF] Lamecker | Atlas to image-with-tumor registration based on demons and deformation inpainting[END_REF]]. Thirion's algorithm was cast as a well posed energy minimization problem via the introduction of hidden variables called correspondences. The idea is that instead of requiring point to point correspondences to be exact, local errors are allowed [Cachier 2003]. This translates in the following energy:

E(C, T ) = 1 σ 2 i Sim(I, J • C) + 1 σ 2 c dist(T , C) + 1 σ 2 r Reg(T ) (4.2)
Where σ c and σ r are weights depending on the spatial uncertainty of the correspondences and on the irregularity of the deformation respectively. Typically, the different terms are defined as

Sim(I, J •C) = I -J •C 2 , dist(T , C) = C -T 2
and Reg(T ) = ∇T 2 . This formulation enables to decouple the problem into two iterative optimization steps:

Correspondence step:

C t = argmin C 1 σ 2 i Sim(I, J • C) + 1 σ 2 c dist(T t , C)
Regularization step:

T t+1 = argmin T 1 σ 2 c dist(T , C t ) + 1 σ 2 r Reg(T ) (4.3)
In order to deal with the presence of the pathology, [START_REF] Lamecker | Atlas to image-with-tumor registration based on demons and deformation inpainting[END_REF]] introduced a third step called the inpainting step. This step overwrites the computed displacement field in the tumor area with the solution of the Laplace equation ∆T = 0. Since no information is given about the image in this area, the deformation is simply interpolated from the boundary of the tumor.

Atlas Seeding

Aside from pathology masking, a series of methods specifically designed for brain tumors have been developed. The common idea is to simulate the presence and growth of the tumor in the healthy atlas. A simple method was proposed by Dawant et al. [Dawant 1999]. First, a strongly regularized demons registration brings the atlas and patient images into approximate alignment. The strong regularization limits the amount of displacements allowed which preserves the brain anatomy and prevents errors due to the presence of the tumor. Tumor voxels in the patient's image are manually segmented and assigned a constant intensity value. A subset of the synthetic lesion voxels (selected by erosion of the segmentation map) are then artificially introduced in the healthy atlas (the atlas is said to be seeded). A second weakly regularized registration is carried out where the seed will trigger demon forces that will simulate the tumor growth in the atlas to match the patient's tumor. This implies however a sizable seed that can mask key anatomical information in the atlas. The work of [Bach Cuadra 2004, Pollo 2005] was a natural extension of this seeding approach. Automatic segmentation of the tumor is carried out using the statistical classification based ATM-SVC algorithm [Warfield 2000, Kaus 2001] (see section 2.5.3 for an overview of the method). Here, the seed (a single voxel in the atlas) is manually chosen by an expert. The registration method follows the demons algorithm but treats the voxels differently if they are inside or outside the tumor.

The displacement field inside the tumor follow a growth model that assumes a radial growth from the seed. The classical demon forces are computed outside the tumor. The displacement field is adaptively regularized using a Gaussian filter of variable standard deviation σ. The tumor area is not regularized (σ = 0), a small value is used at the tumor boundary (σ = 0.5) and smaller deformations are allowed in the rest of the brain (σ = 0.8). A variant was proposed later in [Bach Cuadra 2006] using a Mutual Information flow based registration scheme in the healthy region. The drawbacks of the method include expert seeding and a limitation to homogeneous and uniformly growing meningiomas in terms of growth model and automatic segmentation.

A second group of methods relied on biomechanical models to simulate the deformations induced by the tumor. The work of Kyriacou et al. [START_REF] Kyriacou | [END_REF]] was the first using biomechanical modeling for atlas to patient registration in the presence of a tumor. A Finite Element (FE) model is used to model the brain's deformations. The FE method is a numerical approach to approximate the solution of partial differential equations with boundary conditions by discrete meshing of the domain. Here, the brain boundaries are defined manually and used to construct a triangular mesh of the brain. The different brain structures are modeled as elastic materials with tissue specific parameters. Similarly to [Bach Cuadra 2004], the tumor growth is assumed to be radial and uniform. Given a circular seed of diameter D s , the grown tumor would have a diameter of D t = D s (e 0 + 1), where e 0 is a parameter that specifies the growth rate. Because of boundary constraints defined from the surrounding tissues, the resulting growth of the tumor will tend not to be uniform. First the model shrinks the tumor to an almost inexistent tumor. A normal registration is then carried out between the healthy atlas and simulated healthy patient. The third step simulates the growth of the tumor in the registered atlas via a regression procedure to estimate the seed location in the atlas. This methods inherits a high computational cost and was limited to 2D models.

Mohamed et al. [START_REF][END_REF]] proposed a registration framework that decomposed the deformation field between the healthy atlas and diseased patient into a tumor induced deformation (the mass-effect) and an inter-subject anatomical differences deformation. They used a more sophisticated finite elements biomechanical model of tumor growth [START_REF][END_REF]]. Tumor and edema are modeled as concentric spheres following the idea of [START_REF] Kyriacou | [END_REF]]. The mass effect induced by the tumor is modeled by a constant outwards pressure P from the tumor's boundary, while an isotropic expansion in white matter is assumed for the edema. In order to reduce the computational complexity inherited from the biomechanical model, a PCA based statistical model was used to learn the deformations induced by the tumors for different subjects and parameters (tumor and edema center and radii, and pressure P). Given such a model, an initial guess of the deformation field obtained by normal registration is projected onto the tumor induced deformation space. This enables the estimation of the optimal tumor growth parameters via a maximum likelihood framework and the simulation of the tumor in the atlas. The final step registers the tumor bearing atlas to the patient.

Zacharaki et al. [Zacharaki 2008] built upon this work using a more elaborated registration scheme and method for the evaluation of the tumor growth parameter. Their framework, called ORBIT, adapts the HAMMER algorithm [Shen 2002] to the presence of tumors. The HAMMER algorithm is a hierarchical registration framework that was proposed in order to reduce the ambi-guities of intensity matching driving the registration process to local minima. The main idea was to compute a similarity metric between attribute vectors instead of relying solely on the intensity values. The attribute vectors are computed at each image voxel and consist of edge information, normalized intensities and tissue specific geometric moments invariants (computed on a spherical neighborhood at different scales). In ORBIT, two tumor based attributes were introduced: the signed distance from the tumor boundary and the angular location with respect to its center. The similarity measure between two voxels is defined as a weighted sum of a criterion matching the healthy structures Sim B and a tumor criterion Sim T that matches the tumors' geometry:

Sim(x, y) = (1 -w(x, y))Sim B (x, y) + w(x, y)Sim T (x, y) (4.4)
w(x, y) is a weighting factor that is equal to one in the tumor area and decreases proportionally with respect to the distances to both tumors (one simulated in the atlas and the target's image tumor). Sim B matches the healthy attributes (edges, intensity, geometric moments) while the tumor attributes are matched by Sim T . The registration process starts by selecting the voxels with the most salient features (the driving voxels) that are then hierarchically enriched by new ones to increase the registration accuracy. The biomechanical model parameters are estimated by minimizing an objective function based on the assumption that a poor modeling of the tumor will produce important and unrealistic deformations. The objective function evaluates the quality of the obtained registration using 3 terms: (i) the segmentation overlap of the registered images, (ii) the similarity metric, and (iii) the laplacian of the deformation field which measures its smoothness. A new registration process is then carried out with the optimally estimated parameters. A summary of the iterative approach can be seen in Fig. [4.2]. The method was improved in [Zacharaki 2009] with a computationally efficient biomechanical model [START_REF] Hogea | [END_REF]] that does not require statistical training.

The appeal of growth models is undeniable, as they offer means to obtain anatomically realistic deformations induced by the tumor. However, the quality of the registration depends on the accuracy of the model and of its parameters. Simple radial growths are far from adapted to the many types of tumors while the biomechanical models suffer from their computational cost in modeling and seeking the optimal parameters. They are also prone to errors when dealing with infiltrative tumors, mistaking infiltrated tissue with new tissue. Despite its simplicity, pathology masking is more adapted to the diffusively infiltrative tumors that are DLGGs and could also be used in other clinical contexts since no assumption is made on the pathology nor its progression. A common drawback for both approaches is their dependency on an accurate pre-segmentation of the tumor.

Joint Registration and Segmentation

Healthy Volumes

Coupling registration and segmentation is not a new problem and several attempts have been made in the context of healthy structures. The key idea of these methods is that the structures to be segmented are present in both images. [START_REF] Yezzi | [END_REF]] proposed an active contour based framework where the goal is to segment the same object in two different images I and J. The goal is to find the contour C I of the object in I and the transformation T that maps I to J so that C J = T (C I ). An active contour energy is defined on each image as a sum of region terms.

E I = InCurve f in (C I ) + OutsideCurve f out (C I ) (4.5)
where the functional f in and f out are the difference between the observed intensity and mean intensity in their corresponding region. A rather similar concept was used in [Paragios 2005] where the aim was joint segmentation and sense optical flow estimation. The link between contours C I and C J enables to define a coupled objective function as the sum of both energies E(T , C) = E I (C I )+E J (T (C I )). The type of transformations where limited to rigid and affine transformations that have a finite number of parameters to evaluate. This idea was extended in the context of dense optical flow in [START_REF] Karantzalos | Implicit free-formdeformations for multi-frame segmentation and tracking[END_REF]].

Wyatt and Noble [Wyatt 2003] cast the problem of joint segmentation and rigid registration as a MAP estimation one regarding the optimal segmentation maps S I , S J and deformation field T . Assuming independence between observations I and J, the joint probability is written as: P (S I , S J , T |I, J) ∝ P (I|S I )P (S I )P (J|S J )P (S J )P (T ) (4.6)

A MAP-MRF formulation is used for the segmentation terms. The pairwise term follows the Potts model while the likelihood terms are mixtures of Gaussians. The rigid registration uses the mutual information similarity determined according to the joint class histogram obtained from the segmentation. The method alternates between estimating the optimal segmentation fields for both images by maximizing the segmentation posterior probability and estimating the registration parameters by minimizing the similarity measure.

Ashburner et al. [Ashburner 2005] presented a probabilistic framework for simultaneous tissue classification, atlas registration and bias field correction. The model is based on classifying the image using a mixture of Gaussians. Bias field correction is included in the model by locally adjusting the mean and standard deviation values of the Gaussian distributions. The local variations are governed by an additional parameter β that accounts for smooth intensity variations across the image. Spatial priors are introduced as the ICBM probabilistic atlas values. After an initial affine registration, the matching between the atlas and the image to be segmented is refined by the introduction of a parameter α that allows local deformation of the atlas. Instead of weighting the Gaussians by a constant parameter P (S(x) = k|γ) = γ k , they are weighted by a spatially varying prior:

P (S(x) = k|γ, α) = γ k Π k (α(x)) K i=1 γ i Π i (α(x)) (4.7)
where Π k is the atlas probability value for class k and K the number of classes. The transformation is parametrized as a linear combination of cosine transforms. Smoothing priors for the bias field vector and deformation field are added as multivariate Gaussians of zero mean. The optimal parameters are obtained by minimizing the following objective function using the ICM algorithm [START_REF] Besag | On the statistical analysis of dirty pictures[END_REF]]:

E(Φ, β, α) = - K i=1 log N µ k ρ(β) , σ k ρ(β) P (S(x) = k|γ, α) -log P (α) -log P (β) (4.8)
where Φ corresponds to the set of Gaussian parameters. An Expectation Maximization framework was proposed by Pohl et al. [Pohl 2006] for the estimation of the transformation, structures segmentation and bias field correction parameters. The method follows the segmentation and bias field estimation method of Wells et al. [Wells III 1996b] but integrates affine registration and atlas prior information. The E-step estimates the classes posterior probabilities based on the current bias field and aligned spatial priors (probabilistic atlas maps). The M-step updates the registration and bias field parameters based on the classes probabilities. The registration parameters are computed based on the current segmentation and minimize the Kullback Leibler divergence between the computed segmentation posterior probabilities and the probabilistic atlas prior probabilities. In other words, we seek the best possible agreement between both probability distributions.

Last but not least, Mahapatra and Sun [START_REF] Mahapatra | Integrating segmentation information for improved MRF-based elastic image registration[END_REF]] introduced a 2D MRF based concurrent registration and segmentation method. Each voxel is assigned a label that corresponds to a displacement vector and segmentation label. The unary term is the sum of a conventional registration similarity measure and a term that computes the agreement between the source and target images intensity based classes likelihoods. The pairwise term ensures smoothness of the deformation field by forcing pixels of the same class to have similar displacements and limiting the maximum displacement between voxels of different classes. The MRF is optimized using graph cuts.

Tumor Bearing Images

Joint segmentation and registration methods for healthy volumes exploit the fact that the same objects are present in both volumes. The problem becomes far more difficult in the presence of a tumor where the object sought to segment is absent in one of the volumes. Gooya et al. [Gooya 2012] proposed a Bayesian approach to tackle it, inspired from the work of [Zacharaki 2009] and [Pohl 2006]. They simulate a tumor in the atlas using a diffusion-reaction growth model coupled with an elastic modeling of the brain to simulate the mass effect [Hogea 2008]. An example of modified atlas is shown in Fig. 4.3. The optimal registration and growth model parameters as well as segmentation classes are inferred using the EM algorithm guided by the modified atlas. This joint segmentation and registration work is an extension of their growth model-EM registration framework [Gooya 2011], where only the registration and growth model parameters were inferred through SVM classification. Here, they work on multispectral images (T1, T1CE, T2, and FLAIR) and model the different classes as multivariate Gaussian distributions weighted by the atlas probability maps. The EM E-step computes the classification posterior probability given the current estimates of the Gaussian, growth model and registration parameters that are updated in the M-step.The method suffers from the drawbacks of growth model approaches, that is computational cost and lack of adaptability.

Related work in the context of surgical resection will be introduced in appendix A.

Concurrent Deformable Registration and Tumor Segmentation

Efficient coupled registration and segmentation frameworks, should inherit reasonable computational complexity while being modular with respect to the different components of the algorithm (similarity criterion for registration, tumor detection method, and even clinical context). In such a context, the model of choice was an MRF model that coupled with discrete optimization algorithms inherits such properties. In order to account for the presence of the tumor, we selected pathology masking instead of growth models since it didn't require prior knowledge on the evolution of the tumor and could easily be introduced in an MRF model. Furthermore, the deformation induced by DLGGs tends to be limited due to their infiltrative nature. Contrarily to pathology masking methods that require an initial segmentation and can introduce a bias in the registration process, in our approach the mask is progressively determined during registration and becomes more and more precise as the images get registered.

Discrete MRF Model for Deformable Registration

The formulation of the registration problem as an MRF optimization problem was introduced in [Glocker 2008a]. The problem's complexity was remarkably reduced by the use of a parametric deformation model as well as a discretized solution space. Consider a moving image A defined on a domain Ω and a fixed image I. We seek the transformation T that will map the moving image to the fixed image. In the context of parametric transformation models, Free Form Deformation models [Sederberg 1986] assume a uniform sparse grid superimposed to the moving volume [START_REF] Rueckert | [END_REF]]. The grid's control points constitute a graph G = (V, E), where the nodes correspond to the control points and the edges constitute a neighborhood system N . We define a predefined discrete displacement set ∆ = {d 1 , ..., d n } and a set of labels L = {1, ..., n}. Each graph node p is to be assigned a label l p ∈ L that is associated to a displacement d lp ∈ ∆. The complete deformation field is then estimated by interpolation:

T (x) = x + p∈G η( x -p )d lp (4.9)
where η( xp ) is the projection function that describes the influence of control point p over image voxel x. When the FFD model is adopted, the projection function is defined as cubic Bsplines.

Let us recall that the form of a deformable registration energy to be minimized is the following:

E(T ) = Sim(I, A • T ) + λReg(T ) (4.10)
The optimal labeling l = {l p , p ∈ G} is recovered by minimizing the MRF energy

E(l) = p∈G V p (l p ) + p∈G q∈N (p)
V pq (l p , l q ) (4.11)

The parallel between an MRF energy and the deformable registration's energy is clearly visible.

The pairwise costs play correspond to registration regularization inherent to non-rigid registration methods. The unary costs correspond to the similarity criterion that is to be minimized to match the images. In the work of Glocker et al., it is defined as a point wise criterion:

Sim(I, A • T ) = Ω Sim(I(x), A(T (x))) dx (4.12)
Because of the parametric formulation, the energy could be expressed through terms acting on the grid's control points. To this end, we need to back project the voxel wise similarity measure on the control points:

Sim(I, A • T ) = 1 |G| p∈G Ω η( x -p )Sim(I(x), A(T (x)) dx (4.13)
where η( xp ) is the back projection function that describes the influence of voxel x on control point p. For a voxel wise similarity measure, it is defined as:

η( x -p ) = η( x -p ) Ω η( y -p )dy (4.14)
The unary costs are defined according to the similarity criterion as:

V p (l p ) ≈ 1 |G| Ω η( x -p )Sim(I(x), A(T (x))) dx (4.15)
A key hypothesis of the MRF formulation is that the unary term V p (l p ) depends only on its corresponding node p. However, the formulation of the unary cost proposed in equation 4.15 does depend on the neighboring nodes because of the interpolation term in T (x). It was shown in [Glocker 2008b] that to preserve independence, the unary costs can approximate the similarity value as:

V p (l p ) = 1 |G| Ω η( x -p )Sim(I(x), A(x + d lp )) dx (4.16)
Considering that the fixed volume I features a pathology, the MRF registration method has to be adapted. Following the pathology masking methods, a natural approach would be to ignore the similarity criterion in the tumor area, which will be deformed by interpolation due to the regularization term. Two issues arise from such an approach. First, it is required that the tumor's position is known which implies the need for an initial segmentation. Second, masking the pathology can introduce an initial bias and prevent the correct registration of healthy structures that are close to the tumor. Errors get even more likely when the automatic segmentation of the tumor is not that precise.

Coupled Tumor Segmentation and Registration

Consider a moving image or reference brain atlas A defined on a domain Ω and a fixed tumor bearing brain image I. The objective here is twofold: register the healthy atlas to the diseased subject and segment the tumor. Given the strong dependencies between the two problems, we aim at simultaneously finding the deformation field T and the segmentation map S. Aside from increasing the registration quality, coupling segmentation with registration also benefits the segmentation output as it introduces anatomical knowledge in the segmentation process and reduces the risks of false positives.

We consider a sparse grid G superimposed to the atlas, a predefined set of n discrete displacements ∆ = {d 1 , .., d n } and a discrete label set L c = {1, .., 2n}. The problem of concurrent tumor segmentation and registration is converted into a discrete labeling problem. Assigning a label l p to a grid node p corresponds to assigning a pair segmentation/displacement {d lp , s lp } ∈ ∆ × {0, 1} where :

d lp = d k where k = l p (mod n) s lp = l p n (4.17)
As a result, when a grid node is assigned a label, it is simultaneously assigned a displacement (selected from the predefined discrete set) and labeled as a tumor or background node. The displacements and classification decisions are then propagated on the whole volume by interpolation:

T (x) = x + p∈G η( x -p )d lp S(x) = p∈G η( x -p )s lp (4.18)
The MRF energy to minimize is designed as an association of segmentation and registration terms:

E c = p∈G αV R (l p ) + (1 -α)V S (l p )
Registration and segmentation unary terms

+ p∈G q∈N (p) α V pq,R (l p , l q ) Registration pairwise cost +(1 -α) V pq,S (l p , l q ) Segmentation pairwise cost (4.19)
Where α is a parameter that balances the importance of the registration and segmentation term. The focus is on registration when α is high, and on segmentation when α is low. Linear programming is used to recovered the optimal labeling via the Fast-PD algorithm [Komodakis 2008b].

Unary terms

The registration term seeks for optimal correspondences between the fixed and moving image. Correspondences are sought in the healthy areas by computing a similarity measure Sim(.) between the fixed and deformed image. However, they do not exist in the tumor area (characterized by s lp = 1) where the similarity measure becomes unreliable. We use a constant cost C tm instead whose role is two-fold: first it provides a unary cost independent of the displacement so that the deformation is solely guided by the pairwise regularization term. Second, It identifies regions that are potentially tumors due to a high dissimilarity between the two images: The segmentation term corresponds to the prior likelihood probability of the fixed volume I: p(I(x)|s lp ). It refers to voxel wise classes probabilities computed using standard classification techniques. This term ensures that the nodes with a high probability of being tumor are labeled accordingly. Since the node is segmented and displaced simultaneously, the segmentation term is dependent on the node's displacement d lp :

V R (l p ) = Ω η( x -p )     s lp C tm
V S (l p ) = Ω η( x-p ) -s lp log p I(x + d lp )|s lp = 1 -(1 -s lp ) log p I(x + d lp )|s lp = 0 dx
(4.21) For both unary terms, the back projection function η(.) is used to back project the voxel wise information (class probabilities and similarity criterion) on the control points.

Pairwise terms

The role of the pairwise terms is to enforce smoothness of the segmentation and deformation field. If the resolution of the grid is low (large spacing between the nodes), it makes little sense to impose similar labels to the neighboring nodes. Furthermore, cubic B-spline interpolation on a large resolution grid produces a satisfactory smoothing effect. The pairwise costs are therefore strengthened as the resolution decreases. The segmentation pairwise cost is designed as a simple Potts model weighted by the distance between the nodes:

V pq,S (l p , l q ) = |s lp -s lq | p -q (4.22)
The registration pairwise cost is defined in a similar fashion and imposes local consistency of the displacement field:

V pq,reg (l p , l q ) = λ d lp -d lq 2 p -q (4.23)
where λ describes the influence of the registration regularization. The registration pairwise cost plays a key role in the tumor area as it drives the deformation (due to the constant registration unary cost) and determines the amount of admissible deformation. In healthy areas, λ should be high enough to preserve the anatomical structures. It should however be weakened around the tumor area where strong deformations can occur due to the tumor's growth.

Incremental implementation

The main limitation of such discrete models resides in the definition of the displacement search space and of the grid's resolution. On one hand, we want the algorithm to handle large and yet precise displacements which requires a large and thorough displacement set as well as a important grid resolution. On the other hand, they are both to be constrained in order to maintain the computational burden manageable. Adding segmentation labels in the process increases the difficulty as almost voxel wise labelings are required, especially given the arbitrary shapes of the tumor.

Those drawbacks can be handled by adopting a coarse to fine hierarchical approach. This enables to recover large deformation (low grid resolution) as well as precise displacements and segmentation at high resolution while defining small discrete displacement sets [Glocker 2008a]. Additionally, the propagation of segmentation decisions from low resolution to high resolution makes the segmentation less sensitive to noise and the tumor's heterogeneous appearance. In this setting, grids of increasing resolutions are considered: {G r , r ∈ [1, L]}, where L corresponds to the number of resolution levels. For increased precision, several iterations (MRF optimizations) are performed at each grid level r during which the displacement set ∆ t is progressively refined around the optimal displacement obtained at iteration t -1. At each iteration, the updated MRF energy is computed between the fixed image and the deformed image A • T t-1 on the basis of the deformation field estimated at the previous iteration:

l t opt = argmin l E t c (l |I, A • T t-1 ) (4.24)
The newly obtained deformation field after optimization is then composed with T t-1 . Three challenges arise from such a formulation for concurrent registration and segmentation. First, the segmentation decisions obtained at resolution r -1 should be propagated to the next level. The idea is that images at a low resolution should not suffer from local non homogeneity in tumor appearances resulting in a smoother segmentation than what would be obtained for a voxel wise segmentation. The finer the resolution becomes, the more focus should be paid on the tumor's boundaries. In [Parisot 2012b], we introduced into the MRF energy a cost penalizing classification decisions that are not in accordance with the segmentation map of the previous level:

V t pen (l p ) = Ω η( x -p )s lp S t-1 (x) 2 exp(r) dx (4.25)
However, such a cost is manually set and has to be adapted to the clinical context and potentially, to each patient. Furthermore, erroneous segmentation can be propagated from one level to the next that is not a desired behavior since only the most certain labels should be propagated. Second, the grid resolution levels have to be specified: the first and last resolution as well as the refinement process from one level to the next. This is the key aspect of the incremental approach in terms of efficiency and quality of registration/segmentation. While fine resolutions could lead to precise detection of the tumor's boundary and displacements, they also suffer from a decrease of robustness and increase substantially the computational burden. Typically, the grid is refined by doubling the number of nodes at each level, which limits the final level's resolution. Ideally, high resolutions should only be considered in areas with discontinuities. This is the idea behind the work of Shi et al. [Shi 2012], who proposed the Sparse Free Form Deformations, based on the assumption that it is possible to embed a realistic transformation in a sparse representation using a multi level FFD framework. Indeed, the B-spline coefficients (control points displacements) are often close to zero for the majority of control points. Given a set of m FFD grids of different resolutions, they aim at simultaneously estimating the control points' displacements of the m levels. This is achieved by adding a sparsity constraint on the estimated solution to the registration objective function. This results in a sparse representation of the deformation field across multiple resolutions, where importance is given to control points at different resolutions according to the observed images properties.

The third challenge is the definition and resampling of the displacement set at each iteration. In [Glocker 2008a], the displacement set is defined as an isotropic regular sampling around the control point's initial position. It is refined at each iteration by simply diminishing the maximum displacement allowed. While it allows additional precision, it fails to take into account the local anisotropy of the regions to be registered. An uncertainty based sampling was proposed in [Glocker 2008b] for MRF based registration. The impact of a label swap (different displacement) on the overall energy is evaluated by computing the min-marginal energies [Kohli 2008]. The uncertain directions are estimated as directions where there are small energy variations and the displacements are then sampled according to those directions.

Inspired by this work, we consider local uncertainties to design a content driven adaptive framework that handles the three main issues of the hierarchical approach. The resampling of the displacement sets and grid resolution is guided by the segmentation and registration uncertainties that are computed from the min marginal energies. A segmentation propagation term is also derived from the uncertainties, assuring propagation of reliable segmentations from the lower resolution levels.

Uncertainty Driven Resampling

Local Registration and Segmentation Uncertainties

The local segmentation and registration uncertainties are estimated using the min-marginal energies. The use of the min-marginal energies of an MRF to compute the likelihood of a label assignment was suggested in [Kohli 2008] for segmentation and adapted later in [Glocker 2008b] for registration. More recently, new methods for estimation of such min-marginals/uncertainties have been proposed [START_REF] Tarlow | [END_REF]. By definition, a min-marginal is a function that evaluates the minimum value of the MRF energy under different constraints. In order to evaluate the local uncertainties, the influence of a label selection on the overall energy has to be evaluated. At a given iteration t (associated with resolution r), the min marginal marginal energy associated to a control point p ∈ G r and a label k ∈ L computes the minimum value of the energy E t c with the value of l p fixed to k:

Ψ p,k,t = min l,lp=k E t c (l |I, A • T t-1 ) (4.26)
Recall that an MRF energy is associated with the probability of a configuration P (l) as:

P t (l) = 1 Z exp(-E t c (l)) (4.27)
Where Z is the partition function. The max marginals µ p,k,t (evaluated via message passing algorithms) give the value of the maximum probability given the constraint l p = k:

µ p,k,t = max l,lp=k P t (l) = 1 Z exp(-Ψ p,k,t ) (4.28)
By computing the min marginal energies for all possible labels, we can evaluate the influence of a label swap on the energy with respect to the optimal minimum energy. A measure of confidence can be derived from the max marginals by normalizing over all possible labels:

C p,k,t (l) = µ p,k,t i∈L t µ p,i,t = exp(-Ψ p,k,t ) i∈L exp(-Ψ p,i,t ) (4.29)
Note that the computation of the max marginals' partition function Z is unnecessary since it cancels out when the confidence measure is computed.

In the context of concurrent segmentation and registration, the label information is the product of two label spaces, the displacement and the segmentation space. Considering that label k is different from the optimal label l opt p and that the segmentation label is constant (s k = s l opt p ), the label swap correspond to a local perturbation from the optimal displacement. A substantial energy increase means that the chosen displacement is quite certain in the direction d k -d l opt p d k -d l opt p while a weak energy variation implies that the direction is uncertain and that both labels are almost as likely. A registration confidence measure is computed by normalizing the max marginals over all possible displacements, for a given segmentation label:

C t R (p, k) = exp(-Ψ p,k,t )
i∈L,s i =s l opt p exp(-Ψ p,i,t ) (4.30)

A segmentation confidence measure can be similarly derived from the max marginals by measuring the variations of the energy when the segmentation label changes s k = s l opt p . The confidence of a segmentation label S can be measured by normalizing over all labels:

C t S (p, S) = i∈L,s i =S exp(-Ψ p,i,t ) i∈L exp(-Ψ p,i,t ) (4.31)

Min Marginals Computation using Dynamic Graph Cuts

While using Fast-PD to compute the min-marginals would be the most straightforward method, it would imply performing an optimization for each fixed label which would considerably slow down the algorithm. We compute an approximation of the min marginals using dynamic graph cuts, following the method proposed by [Kohli 2007]. We recall here that graph cuts rewrite the labeling problem as a graph partitioning one. Each node is connected to a source and sink node, each terminal node corresponding to a label. The optimal labeling is found by solving the mincut/max-flow problem, that is, finding a partition of the graph by removing the edges of minimal weights. This is equivalent to finding the maximum flow from the source edge to the terminal edge with respect to the edges capacities. Dynamic graph cuts minimize similar functions but with slightly different energy terms. The maximum flow computed for the original problem is recorded and used as initialization for the new, slightly modified problem.

The min-marginals Ψ p,k,t are computed from the graph's source/sink (depending on the fixed label) flow potentials that compute the maximum amount of flow that can arrive/leave the node p. [Kohli 2008] show that they are computed as:

Ψ p,k,t = f max + f T (k) p (4.32)
where T (k) is the terminal node corresponding to label k, f max is the maximal flow computed by minimizing the complete MRF energy E t c and f

T (k) p
is the flow potential. The flow potentials are computed using dynamic graph cuts by adding an edge of infinite cost between node p and the terminal node T (k) to the residual graph to force the label assignment, and then computing the updated graph's min-cut. While exact solutions can be obtained for binary label sets, the min marginals can only be approximated for mutli-label configurations. They are approximated in the multi-label setting by defining 2n binary energies (one per label) where the energy k is computed between the optimal labeling and label k.

Content Driven Adaptive Resampling

Given the segmentation and registration confidence scores, we use this information to design an image driven displacement set and grid resolution dependent on the local image properties. The displacement set construction is guided by the local registration uncertainties, while the grid resolution depends on the homogeneity of the regions covered by the nodes and the associated segmentations confidences.

Displacement Sampling

Consider that the current displacement set ∆ t is a dense discretization of the displacement search space defined by a maximum displacement and a sampling rate. Each min-marginal energy Ψ p,k,t (and its corresponding confidence measure) is associated with a spatial location in the Euclidean space centered at node p's initial position. This spatial location is defined by the position of the node after a displacement of d k . This yields a confidence map describing which displacements are most likely. We follow the methodology of [Glocker 2008b] and approximate this map as a Gaussian distribution centered at the node's optimal displacement. The distribution's covariance matrix illustrates the directions towards which the displacements are uncertain (small energy variation) and is used to intelligently sample the displacement set ∆ t+1 according to the covariance's scale and main axes. Note that we do not carry out this process for the tumor label and when the focus is on segmentation (low α) since the registration is mostly driven by the pairwise cost that yields a spherical covariance matrix. An example of adaptive displacement sampling can be seen in Fig. [4.4].

Grid Refinement and Segmentation Propagation

We define a uniform grid G r,max of resolution M ×N ×P and an activation function A r : G r,max → {0, 1} that describes the resolution of the current grid G r ⊂ G r,max . A node c r ∈ G r,max is a part of G r if A r (c r ) = 1. In the remainder of this chapter, the nodes characterized by A r (c r ) = 1 will be referred to as activated nodes. After optimization, the min-marginal energies are estimated for the activated nodes in G j .

At the next resolution level r + 1, G r,max is refined as a uniform grid G r+1,max of resolution 2M -1 × 2N -1 × 2P -1 by splitting the edges of G r,max in two. The new level's grid resolution G r+1 is determined by activating selected nodes that are likely to increase the precision of the segmentation and registration. There are three possible ways of activating a node at level r + 1. First, a node p ∈ G r+1,max will be activated if it has the same spatial coordinates with a node c r ∈ G r,max that is activated (A r (c r ) = 1) at the previous level. Because of the way G r+1,max is defined with respect to G r,max , all low resolution nodes have a direct spatial correspondent in the higher resolution grid. The activation of the remaining nodes depends on the min-marginal energies evaluated at the low resolution level. A node p can either be activated for registration or for increased precision of the segmentation. The registration and segmentation activations are governed by two activation terms A R (p) and A S (p) respectively. In order to propagate the max marginals confidence scores and activation information between levels, we construct an inter-level 3D neighborhood N i (.) that connects the node c r ∈ G r,max to its 27 closest neighbors in G r+1,max , where the proximity is evaluated based on the nodes' spatial coordinates. A 2D representation of the neighborhood can be seen in Fig. [4.5].

The registration activation term (A R (p)) aims at identifying the presence of structures and discontinuities in the region surrounding the node. The idea is that precise displacements (and therefore a higher resolution) are necessary at boundaries while no increased focus is necessary in homogeneous regions. As a result, a node should be activated if it is at the boundary between different structures. If discontinuities are present, strong variations of the min-marginal energies are expected with respect to the displacement label, due to the different structures' properties. The registration activation criterion of a node p is defined based on the min-marginal energy range of its neighbors in G j given their optimal segmentation label:

A R (p) = H   cr∈N i (p) 1 N max k∈L,s k =s l opt cr Ψ cr,k,t - min k∈L,s k =s l opt cr Ψ cr,k,t -µ r   (4.33)
where H(.) is the Heaviside step function, N is the number of nodes in the neighborhood of p and µ r is the mean energy range over all nodes in G r . In other words, a node p will be activated if the mean min-marginal energy range between its neighbors in G r is larger than the overall mean energy range µ r . Similarly, the segmentation activation term (A S (p)) exploits the information of the min-marginal energies. A confidence score C t S (c r , S) is computed for each node in G r given a label S by normalizing the min marginals. This allows to evaluate how reliable a node's segmentation is. In the case of a binary segmentation, the confidence score can be converted into a node wise uncertainty measure independent of the label:

U t (c r ) = 1 -|C t S (c r , S) -0.5| (4.34)
This term measures the degree of uncertainty associated to the chosen segmentation label. A low value of U t (c r ) is desirable as it corresponds to highly reliable segmentation decisions. We use this information to propagate the recovered optimal segmentation labels to next levels based on their reliability. This is done by adding an inter-level potential to the MRF energy defined as:

V ucy (l opt cr , l p ) = 0 if s lp = s l opt cr -log(U (c r )) Otherwise (4.35)
This potential is a penalty cost that penalizes nodes whose the segmentation labels do not agree with their neighbors'. The strength of the penalty depends on the confidence score of the estimated segmentations.

The segmentation activation criterion can easily be derived from this propagation potential: nodes should be activated if their inter-level penalty is low. It is defined as:

A S (p) = 1 -H   cr∈N i (p) (-1) s lp V ucy (l opt cr , l p ) -tsh   (4.36)
Where N is the number of nodes in N (p), tsh is a constant parameter that evaluates the strength of the penalty and H(.) is the Heaviside step function. The neighborhood system is designed so that a node in G r+1,max is connected to several nodes in G r,max . As a results, the penalty of a node is low in two cases. The first and most straightforward case is when the segmentation labels of all the node's neighbors are uncertain. The second case occurs when the neighbors have highly reliable but different labels. This can occur if the node is at the tumor's boundary. In that case, the node will be equally penalized for both labels and the penalty will have almost no impact on the labeling. The segmentation activation criterion measures the strength of the penalty imposed on node p. The case where the labels are equally confident is taken into account by subtracting the penalties associated to different labels. A node will be activated if its penalty strength is smaller than a predefined threshold tsh.

Eventually, the complete MRF energy computed at resolution r and iteration t is written as follows:

E t c (l|V, A • T t-1 , G r ) = 1 |G r | p∈Gr V p (l p |V, A • T t-1 ) + p∈Gr q∈N (p) V pq, (l p , l q ) + p∈Gr c∈N i (p) 1 N V ucy (l opt c , l p ) (4.37)
where N is the number of nodes in G r-1 that are connected to node p ∈ G r .

To summarize, at each iteration t, a new grid G r and energy are defined based on the deformation field and min-marginals estimated at iteration t -1. The first two terms of the energy correspond to the MRF unary and pairwise term. The third term is an inter-level pairwise potential that propagates the previous level's segmentation decisions based on their confidence determined according to their min-marginals. The energy is then minimized and the min-marginals are computed for the next iteration t + 1.

Experimental Validation

Data-Set Preprocessing

The healthy reference subject was a 3D FLAIR volume of size 256 × 256 × 24 and resolution 0.9 × 0.9×5.5mm 3 . The method was tested on a data-set of 110 3D FLAIR volumes of patients featuring DLGGs. All tumors have been manually segmented by experts. The preprocessing step follows the methodology of chapter 3 by regularizing the volumes intensities to the same median and interquartile range as the reference pose. An initial rigid registration of all volumes to the reference pose is performed using the block matching method described in section 3.2.1.1. Eventually, all volumes where skullstripped prior to the registration.

The data-set was split in a training set of 36 volumes randomly selected and a testing set of 74 volumes. The tumor likelihood probabilities are obtained by boosting learning on the 36 volumes. The same features as the ones used in chapter 3 are used, that is, intensity patches and statistics, Gabor features and brain symmetry. The rigid registration step provides a common approximate symmetry plane for all volumes based on the reference pose. To evaluate the potential of our method, the joint framework was compared to a sequential approach. The tumor is first segmented using a pairwise MAP-MRF model, where the likelihood corresponds to the boosting probabilities and an Ising model is used as pairwise cost. Second, a discrete MRF deformable registration is performed by masking the pathology using the obtained segmentation map. Additionally, the potential of the uncertainty driven method is tested by comparing the segmentation and registration results to the ones obtained without uncertainties using a user defined segmentation propagation penalty cost. The results obtained with a maximum resolution and a uniform grid of the same final resolution as the adaptively sampled grid's final resolution are considered.

Implementation

All methods were compared using the set of parameters that yielded the best results. The coarse to fine hierarchical approach consists of 4 grid levels and 3 image levels. The image resolution increases as the grid resolution increases via a Gaussian pyramid, and remains constant at the fourth grid level (see Fig. [4.6]. The resolution increases from 9 × 9 × 5 to 65 × 65 × 37 (maximum). B-splines were used as a projection function, while the inverse projection η was defined as a linear one as suggested in [Glocker 2008b] to limit the impact of the approximation of the unary term. Because of the interpolation, the obtained segmentation field is not binary and was thresholded at 0.5. The registration similarity criterion was the Sum of Absolute Differences. The parameter α is initially set to 1 and is progressively diminished to 0.015. The idea is to give an increasing importance of the presence of the tumor on the registration to avoid an initial bias, while the focus is on segmentation at the highest resolution level that is not necessary for increased registration quality. Only the segmentation labels are considered at this level. The constant cost for registration increases with the grid and image resolution from 5 and 6 times the mean similarity cost value with and without uncertainties respectively. The parameter λ describing the influence of the registration regularization was set to 20 and relaxed in the tumor area to allow for substantial deformations due to the tumor presence. Lastly, the threshold tsh for node activation was set to 1.6, corresponding to a mean uncertainty measure of 0.2 (if neighboring labels have the same value). The same parameters where used for all volumes.

At each grid level, 3 iterations were performed. If the min-marginals are not computed, the displacement sampling is sparse and amounts to 31 labels sampled along the main axes. The sampling is refined at each iteration by reducing the maximum displacement. A dense sampling is adopted to compute the uncertainties at the first iteration, which correspond to 1331 labels in 3D uniformly sampled. Once computed, a sparse sampling is used for the last two iterations. The labels are sampled along the min marginals covariance matrix main axes. The registration uncertainties are not computed at the last two grid levels when α is low and the anisotropy cannot be captured. At these levels, a sparse sampling is used for all iterations. To prevent inadequate displacements towards the tumor area, a small pairwise penalty cost is added to the corresponding labels of nodes at the tumor boundary. The two sampling strategies can be visualized in Fig. [4.7] Given the boosting score likelihoods, the selected parameters yield a run time of about 40 seconds with uniform and sparse sampling and about 3-4 minutes with dense/sparse sampling.

Results

In order to evaluate the use of uncertainties for adaptive sampling, we computed the percentage of activated nodes with respect to the maximal resolution. The mean percentage over all nodes for all levels can be seen in Fig. complexity of the framework is considerably reduced, notably at the last, most computationally expensive levels where only 32% and 20 % of the nodes are activated. The relevance of the selected nodes is evaluated visually and several examples can be seen in Fig. [4.9] for different resolution levels. Nodes are activated around the brain's structures and the tumor, where the registration and segmentation have to be precise. This demonstrates the adequacy of the registration and segmentation activation terms. The impact on the run time and memory consumption can be evaluated through a complexity analysis (excluding the time required to build the data term potentials). We consider an MRF with L, E and N being respectively the number of labels, edges and nodes.

• Computational cost: O L × O E × N 2 ∼ O L × E × N 2 per iteration • Memory cost: O (L × O (N + E))
Reducing the number of nodes to approximately 20% at the finer resolution scale leads to:

• Computational cost: O L × O 0.2E × (0.2N ) 2 ∼ O 0.008L × E × N 2 . When tak-
ing into account the number of iterations, we can obtain a complexity that is approximately 3-4 orders of magnitude lower.

• Memory cost: O (L × O (0.2(N + E)))
, approximately one order of magnitude lower.

The quality of the segmentation was evaluated using the Dice Score, False Positive Rate (FP), True Positive Rate (TP) and Mean Absolute Distance between contours (MAD) with respect to the manual segmentation. We compared the segmentations obtained by simple MAP-MRF segmentation, concurrent segmentation and registration at high and low resolution, and uncertainty driven concurrent segmentation and registration. The joint framework shows increased performance in The registration quality is evaluated mostly qualitatively. The individual registration method with masked pathology (MAP-MRF segmentation) is compared to the joint framework without uncertainties at different resolutions as well as with the graphical model endowed with uncertainties. For quantitative evaluation, the ventricles where segmented on 33 volumes where the tumor was close to the ventricles. The Dice and MAD scores, FP and TP rates were computed between the segmented and registered ventricles using the different methods outside the tumor area. Quantitative results show equivalent performance between the individual and joint framework. The performance is preserved with the adaptive sampling while the quality diminishes when an equivalent (in term of nodes) uniform low resolution sampling. The quantitative results are shown in Fig.

Discussion

This chapter introduced a concurrent tumor segmentation and registration framework that directly exploits the dependencies between both problems. We adopt a discrete Markov Random field model, where the tumor is detected by boosting classification as well as a segmentation through registration term that detect areas with high dissimilarities. The registration similarity criterion is relaxed in the tumor area and as a byproduct of that, this area is deformed by interpolation with the displacement at the boundary. The principal shortcomings of discrete methods are addressed via the evaluation of the local segmentation and registration uncertainties and their exploitation for adaptive displacement set and sparse grid sampling. Segmentation results show a notable improvement towards eliminating false positives and an increased robustness due to the hierarchical approach, while the registration is smoother and preserves the deformed image's anatomy in the tumor area. The content driven grid sampling strongly reduces the complexity of the model with a maintained performance.

The choice of pathology masking for registration (justified by the infiltrative nature of DLGGs), along with the discrete formulation, results in a fast and versatile algorithm. Indeed, the method is modular with respect to the registration similarity criterion and the tumor detection method. No assumption is made on the pathology nor its progression, which makes the method easily adaptable to other pathologies and clinical contexts where correspondences are missing. It should be noted that in the context of fast growing, space occupying lesions, integrating tumor growth information would be necessary. Such tumors can impact the brain structures where a relaxation of the pairwise cost at the tumor boundary wouldn't be sufficient.

Several improvements of the algorithm could be considered. Our current implementation of the adaptive sampling simply discards the inactive nodes in the evaluation of the optimal labeling and associates them with infinite costs to changes of labels. This evaluates the impact of reducing the grid's complexity on the quality of the results but does not yet impact the run time. A direct construction of the non-uniform grid would strongly decrease the run time (both in terms of optimization and data term construction) and increase the robustness and quality of the results as active nodes would have direct edges connecting them.

The tumor segmentation results remain dependent on the initial boosting detection and partial segmentations will be obtained if the tumor is poorly detected. A first step in order to obtain better segmentation would be to focus on the detection step. The two most straightforward ways of improving the quality of the detection are the addition of multimodal information such as T1 and T2 volumes that can be easily introduced in the classification feature vector, and assigning the different parts of the tumor (notably necrosis and active tumor) to different classes. The discrete labeling formulation makes the extension to multi-class segmentation straightforward.

The algorithm presented is modular with respect to the registration similarity criterion, but optimal interdependency between segmentation and registration can only be obtained using intensity based similarity measures, that can roughly detect the tumor area via the sole registration unary term. The algorithm is therefore better suited to monomodal registration. An extension of the method to multimodal registration in the context of tumor resection is presented in appendix A. Chapter 5

Conclusion

The main objective of this thesis was the development of methods for tumor segmentation through atlas based prior knowledge. Two approaches have been proposed: one that exploits a tumor specific atlas as a spatial prior, and a second exploiting anatomical knowledge. The presented contributions are not limited to tumor segmentation as we proposed tools for the study of diffuse low-grade gliomas, registration with missing correspondences, concurrent registration and segmentation, and adaptive sampling for hierarchical grid based segmentation and registration.

Contributions

The first contribution of this thesis is the construction of a probabilistic atlas from 210 different patients mapping the preferential locations of DLGGs. This is carried out through a graph structure encoding the spatial proximity between tumors over 210 affinely registered volumes. The preferential locations are obtained by regrouping nearby tumors in clusters. This yields an elevencomponents graph where each cluster corresponds to a preferential location and potential position dependent behavior. This atlas paves the way for studies on the genesis and natural history of the tumors. Understanding the tumors better is key towards optimizing their treatment for a given patient.

A direct application of the atlas was proposed for tumor segmentation and characterization via a coupled Markov Random Field model. The method simultaneously assigns a new tumor to one of the atlas' clusters and labels the tumor voxels. The atlas plays the role of a prior on the tumor voxels' localization based on the spatial repartition of tumors voxels in the cluster. This can be seen as the association of local and global information where determining the tumor's global category (cluster) allows a more refined local segmentation based on the knowledge of its position and expected spatial extension directions. It is clear that this method is not limited to atlas based lowgrade glioma segmentation but can be considered in any segmentation process where the objects to be segmented tend to appear in predefined locations or follow a repetitive structural pattern.

The second important contribution of the thesis is the design of a concurrent segmentation and 119 registration method with missing correspondences. The method is among the first to tackle the problem of simultaneous registration and segmentation when the object to be segmented is absent from one of the images. The deformation and segmentation fields are estimated concurrently on a sparse grid using an MRF model where both solution spaces are interconnected. The proposed method is characterized by its modularity and scalability. It is indeed, designed in a very general way. Almost no assumptions are made on the pathology and on the object detection method, as well as the registration criterion. This modularity is notably highlighted in appendix A where a direct application of the concurrent registration and segmentation method is proposed in the context of surgical resection. A strong asset of the method is its speed, which justifies its extension to surgical interventions.

In the context of concurrent registration and segmentation, we proposed a technique in order to address one of the most important shortcomings of discrete and hierarchical registration and segmentation methods that is the trade off between precision and complexity. We proposed an adaptive sampling approach based on local uncertainties. The precision (higher resolution) is determined based on the reliability of the evaluated segmentation/registration labels. This leads to a considerable decrease of the complexity at high resolutions with a focus on areas with discontinuities and tumors boundaries. Furthermore, it provides a qualitative manner for the interpretation of the obtained results. The combined registration and segmentation method is designed in a way that allows its application to the sole problem of registration or segmentation.

Future Work

The most straightforward future development is the association of the two proposed segmentation methods in one single framework. Using manual segmentations as detection probabilities, the concurrent segmentation and registration (CSR) framework could be used for the construction of a deformable low-grade glioma atlas on the preferential locations of the tumors. The use of deformable registration would increase the precision as well as allow studies on the impact of the tumors on the brain structures and functional areas. This deformable atlas could then be incorporated as a spatial prior in the CSR method, while cluster assignment could be based on initial detections in the hierarchical framework.

Throughout this thesis, a single subject anatomical template has been used as a reference pose. This introduces a bias in the construction of the atlas as well as increases the difficulty of the registration problem when the patient and reference have strong anatomical differences. The construction and use of an averaged population template of the FLAIR modality would be more accurate and adapted to the latter problem. As it encompasses the inter subject variability, deformations would be more natural. Such an atlas would also reduce the bias on the former problem. Furthermore, atlas construction has raised a lot of interest as there has been many attempts at eliminating this bias [START_REF] Zollei | [END_REF][START_REF] Bhatia | [END_REF], Joshi 2004]. A notable approach for MRF based registration [Sotiras 2009] could be considered in our context.

The use of growth models [Jbabdi 2005, Konukoglu 2010, Rekik 2013] could be considered in order to increase the quality of both the atlas construction and CSR framework. Indeed, growth models would enable a more precise localization of the origin point of the tumor which could differ from the center of gravity [Rekik 2013], while including growth models in the registration/segmentation framework could prove useful in the context of heavily deformed brains. Multimodality would most likely bring useful additional information for both segmentation and registration. The different tumor appearances on multimodal images would make the tumor detection step more accurate and consequently increase the quality of the overall segmentation results. Diffusion Tensor Imaging is a very promising modality that could be relied upon for registration and segmentation. While providing structural information for increased robustness and anatomical structures preservation, it can constitute a strong cue on the tumor's (or resection cavity) position as the tumor tends to displace or disrupt the fibers' natural organization. Integration of this knowledge in the CSR framework could have a strong positive impact on the results.

Last but not least, the use of recent parallelized message passing optimization techniques [START_REF] Alchatzidis | Efficient parallel message computation for MAP inference[END_REF]] would further speed up the methods towards real time as well as provide minmarginals computation in a single optimization step. provide important information for neurosurgical guidance. Direct acquisition of such data during surgery is currently not an option due to long acquisition and processing times, and the presence of strong distortions.

For those reasons, aligning the preoperative data to the intraoperative image via non-rigid registration is highly desirable, while the segmentation of the resected area would enable direct evaluation of a potential residual tumor. Similarly to the lacking correspondences induced by the tumor with respect to a healthy subject, the registration task is here violated by the presence of a resected area that is not part of the preoperative image. Specific methods have to be designed in order to accommodate for the presence of this cavity.

Chapter 4 introduced a concurrent tumor segmentation and registration method in the presence of DLGGs that is easily adaptable to other clinical contexts. This chapter presents an extension of the method for the context of surgical brain tumor resection. We propose a concurrent preoperative to intraoperative registration method that is coupled with a segmentation of the resection cavity. A probabilistic formulation is proposed, where the registration task relies on the intensity joint histogram and the resected area is detected based on its intensity profile as a mixture of Gaussians. The knowledge on the position of the tumor in the preoperative image is exploited as a spatial prior on the position of the resection cavity.

The remainder of this chapter is organized as follows: section A.2 presents an overview of the non-rigid registration methods in the presence of a tumor resection, including methods that are coupled with a segmentation of the resection area. Section A.3 presents a MAP-MRF probabilistic interpretation of the concurrent segmentation and registration method while section A.4 describes the definition of the segmentation and matching terms. Experimental validation is part of section A.5 and Discussion concludes the chapter.

A.2 Related work

This section presents prior work on non-rigid registration with tumor resection. The general case of pathology masking presented in section 4.2 could also be applied in this context. Aside from the related work on registration with missing correspondences, there has been little focus on tumor resection and most specifically preoperative to intraoperative non-rigid registration.

Clatz et al. [Clatz 2005] presented a block matching method coupled with a biomechanical modeling of the brain. The block matching method, introduced in chapter 3, seeks correspondences between images blocks by minimizing local similarity measures. This results in a sparse estimation of the deformation field that is in a second step estimated on the whole volume. With the objective to have as little computational overhead as possible during surgery, several preprocessing steps are carried out based on the preoperative image. The brain, ventricles and tumor are segmented on the preoperative image and those segmentations are used to build a biomechanical model of the brain based on the Finite Elements model. In order to avoid comparing homogeneous and ambiguous blocks, a block rejection step is also carried out preoperatively, by measuring the block intensity variance and selecting only the blocks with the highest variance. A prohibited connectivity criterion is included to prevent selected blocks to have overlapping structures. Eventually, blocks are selected outside the tumor area to prevent matchings with the resection cavity. Last but not least, they construct a structure tensor that measures the smoothness of the intensity distribution and provide a measure of confidence of the block's displacement. During surgery, the matchings between blocks are estimated using the Cross Correlation. The dense deformation field is then estimated using the biomechanical model. The influence of a block on the deformation field depends on the value of the similarity metric and the direction of the displacement with respect to the tensor structure (the matching should be colinear with the intensity gradient direction in the block). Similarly to the method of [Ourselin 2000], an additional block rejection step using least trimmed squares is included to reject resection blocks. The method was tested "offline" on 6 pre and intraoperative T1 images, and more thoroughly in [START_REF] Archip | [END_REF]] on 11 patients during surgery. In [START_REF] Archip | [END_REF]], the data was also fused with fMRI and DTI data.

Risholm et al. [Risholm 2009[START_REF] Petter Risholm | Validation of a nonrigid registration framework that accommodates tissue resection[END_REF] proposed a registration framework based on the bijective demons algorithm. They adopt a fluid anisotropic regularization, where the amount of regularization differs from one direction to the next. The resected area was detected using levelsets. The speed function α defining the evolution of the front along its normal direction is defined based on the intensity disagreement between the source and target image:

α(x) = |J(x) -(I • T )(x)|e -|∇I| σ -L (A.1)
I and J are the preoperative and intraoperative images respectively, T is the transformation and L is a threshold parameter. Given this definition, the front expands when the intensity disagreement is strong and retracts when it is small, based on threshdold L, while the exponential term's role is to eliminate the propagation of the front at large gradients. The estimated resection area is exploited for registration as a diffusion sink. The idea is that instead of propagating the computed demon velocities over the whole image, the forces computed inside the resection area (that are not reliable) do not propagate beyond the boundary. In other words, diffusion is allowed towards the resection cavity's boundary, but not from it. This diffusion sink is simply defined by setting the anisotropic regularization parameters to 0 in the appropriate direction. The algorithm alternates between performing the registration given the segmented resected area and diffusion sink constraint and evaluating the resection area. Experiments were performed on 2D data in [Risholm 2009] and 3D in [START_REF] Petter Risholm | Validation of a nonrigid registration framework that accommodates tissue resection[END_REF]]. An Expectation-Maximization framework was introduced in [Liu 2010]. They proposed a point based registration framework inspired from the method of [START_REF] Chui | [END_REF]]. Given an initial set of points extracted from the source and target image (canny edge detection) and an estimation of the resected area M t , the method estimates the correspondences between the points as a normal distribution on the intensity difference between both sets of points. It also registers the images by minimizing the distance between the matched points (E-step). The points within the resected area are not taken into account (i.e. they are removed from the set of points in the source and target image). The deformation field on the whole image is recovered through biomechanical modeling and the Finite Element method. The M-step reestimates the resected area using the deformed mesh and a segmentation of the brain. Elements in the mesh that are contained in the background and constitute the largest connected mesh with the currently estimated resection mesh M t are added to the updated resection mesh M t+1 . Points of the source image that belong to this mesh are removed from the set.

Chitphakdithai et al. [START_REF] Chitphakdithai | Non-rigid registration with missing correspondences in preoperative and postresection brain images[END_REF]] adopted a Bayesian formulation inspired from the work of [Pohl 2006] for the registration of preoperative and intraoperative images, as well as the segmentation of the brain and the resection cavity. The registration problem is posed as a marginalized MAP estimation of the deformation field:

T opt = argmax T log S p(T , S|I, J) (A.2)
The problem is solved using the EM algorithm, where iteratively, the E-step estimates the segmentation map S given the current estimation of the transformation T t and the M-step updates the transformation T t+1 given the estimated segmentation map. The M-step is written as follows: where the label "l" can correspond to three classes: background, brain or resection. The term p(I(T (x))|J, S = l) corresponds to the similarity criterion for registration. In areas with correspondences (background and brain), it is modeled as a normal distribution with respect to the intensity difference between both images, with a class specific standard deviation. It is set as a constant in the resection area. The term p(J(x)|S(x) = l) is an intensity-based prior on the intraoperative image's intensity given a class. It is set as a uniform distribution in the classes that have correspondences and defined as a Maximum Likelihood estimate in the resection area based on a set of training data. Last, the transformation prior p(T ) is a smoothness term acting on the deformation field, expressed through Free Form Deformations. The E-step computes the segmentation posterior probabilities. The method was tested in synthetic data and 6 pairs of 3D clinical data using Leave One Out cross validation for resection probability learning. The method was extended in [START_REF] Chitphakdithai | [END_REF]] by the introduction of an MRF smoothing prior on the segmentation field.

T t+1 = argmax

A.3 Probabilistic Formulation

This section proposes a probabilistic adaptation of the concurrent segmentation and registration MRF framework. Let us consider the source preoperative image I of a subject with a tumor defined on a domain Ω, and the intraoperative target image J of the same subject during the resection of the tumor. We assume that the tumor in the preoperative image has been segmented prior to the surgery where S I is the segmentation map. We seek the transformation T that maps the preoperative image to the intraoperative image as well as the resection cavity segmentation map R. A sparse grid G is superimposed to the source volume I, and each grid node p is to be assigned a label

The classification can be performed in a supervised or unsupervised setting. In the supervised setting, a class dependent GMM is learned on the intensity distribution over all training samples of the same class. This yields a per class probability dependent on the image intensities. The unsupervised setting learns a GMM on the whole image's intensity distribution and a class can be associated to one or several learned Gaussians. In both cases, the Gaussians parameters are estimated from the set of samples X = {x 1 , ..., x n } by maximizing the log likelihood function using the Expectation-Maximization algorithm:

Expectation Maximization Algorithm

The Expectation-Maximization algorithm [START_REF] Dempster | [END_REF]] is an iterative procedure for maximum likelihood estimates in the presence of missing data. Consider X the observed data and Z the missing data. Their joint distribution p(X, Z|θ) is governed by a set of parameters θ and its maximization is supposed to be straightforward (complete log likelihood). Since Z is not observed, we need to maximize the incomplete log likelihood function:

log p(X|θ) = log Z p(X, Z|θ) (A.12)
Unfortunately, the optimization is not trivial due to the presence of the sum that results in complicated expressions. The problem can be approached by acknowledging that for any distribution q(Z) defined over the missing data, we have: log p(X|θ) = L(q, θ) + KL(q p) where L(q, θ) = Z q(Z) log p(X, Z|θ) q(Z) KL(q p) = -Z q(Z) log p(Z|X, θ) q(Z) (A.13) KL is the Kullback-Leibler divergence between distributions q(Z) and p(Z|X, θ) and is therefore always positive, and null when q(Z) = p(Z|X, θ). L(q, θ) constitutes thus a lower bound on log p(X|θ). The EM algorithm aims at maximizing this lower bound and consists of two steps that are repeated iteratively. Given a prior estimation of the parameters θ (t) , the Expectation step (E-step) maximizes the lower bound L(q, θ) with respect to q(Z) while holding θ fixed. Following equation A.13, this maximum is obtained when q(Z) = p(Z|X, θ (t) ). The Maximization step (M-step) maximizes L(q, θ) with respect to θ given the current estimation of q(Z): L(q, θ) = Z p(Z|X, θ (t) ) log p(X, Z|θ) -Z p(Z|X, θ (t) ) log p(X, Z|θ (t) ) = Q(θ, θ (t) ) + cst (A.14) Q(θ, θ (t) ) is defined as the expectation of the complete log likelihood and is what is maximized during the M-step. In summary, the EM algorithm alternates between computing a lower bound of the log likelihood for the current parameters, and maximizing this lower bound to obtain the new parameter values. E-step q (t+1) (Z) =p(Z|X, θ (t) ) M-step θ (t+1) = argmax θ Q(θ, θ (t) ) (A.15)

A.4.1.3 EM for Gaussian Mixture Models

For Gaussian Mixture Models, we can introduce a K-dimensional binary random variable z that corresponds to the hidden data [START_REF] Bishop | [END_REF]]. Each observed variable x n is associated with a hidden variable z n that defines which probabilty distribution the sample follows. Taking into account this variable, the mixture coefficients can be redefined as π k = p(z k = 1) and can be seen as the prior probability of z k = 1. Estimating the Gaussians parameters amounts to estimating the hidden variable's posterior probability or responsibility γ(z n k) during the E-step given the current The term n γ(z nk ) can be interpreted as the number of samples that have been assigned to the distribution k.

In this chapter, the approximate resection cavity area (delineated by the spatial prior S I,s ) is modeled very simply as a mixture of 4 Gaussians: two classes for residual tumor and brain, and two classes for the cavity. The cavity is then isolated by selecting the Gaussian with the brightest and darkest means. This results in a probability likelihood for the resection cavity defined as:

A.4.2 Background area

In the background area, the likelihood term plays the part of the similarity measure that will map the preoperative image to the intraoperative image. The likelihood registration term should take the form of a probabilistic measure that is directly comparable with the resection detection probability. Because of the strong differences in the acquisition process, images of the same modality still have very different intensity profiles. Intensity based similarity measures cannot be used in this situation. A natural probabilistic measure of the differences between the source and target image is to compute the joint intensity histogram between both images. The intensity values of the source and target image are partitioned in a number N of predefined bins B I = {B I (i)} i∈1,..,N , B J = {B J (i)} i∈1,..,N , each bin corresponding to an intensity range. The joint histogram H(B I (i), B J (j)) is computed by counting the number of voxels in I and J that have a value corresponding to the range defined by bins B I (i) and B J (j). When normalized, the joint histogram yields a probability distribution for pairs of intensities:

p(I(x), J(x)) = H(B I (i x ), B J (j x )) N IJ (A.21) where N IJ is the number of intensity pairs counted in the histogram, i x and j x are the bin indices associated to intensity values I(x) and J(x) respectively. A sharp histogram is characteristic of well registered images, as the matching intensities for a tissue remain concentrated in a small histogram location. Mutual information and joint entropy measures aim at directly minimizing the histogram dispersion. Maximizing the joint histogram yields the same effect. The source image is deformed so that its intensity pairs with the target image correspond to bin values that have a high probability. In other words, it reduces the histogram's dispersion.

A.5 Experiments

A.5.1 Implementation Details

We adopted a coarse to fine incremental hierarchical approach. We use 4 grid levels and 3 images levels, with 3 iterations at each level. The grid resolution ranged from 8 × 8 × 6 to 57 × 57 × 37. At each iteration, the joint histogram is re-estimated, exploiting the fact that the images should have stronger correspondences as they get registered. In order to avoid introducing a bias due to the high number of voxels in the background, we do not take the background into account when evaluating the joint histogram. The method was tested using the implementation proposed in chapter 4 that does not include computation of registration and segmentation uncertainties. We therefore use a sparse sampling for all iterations and resolution levels, as well as a manually set penalty term defined as V t pen = Ω η( xp )0.2R(x)s lp e -l where t is the iteration and l the resolution level. This penalty term had to be manually edited to the different context, which highlights the interest of the uncertainty based propagation.

The number of bins was set at N = 64 at the coarsest resolution level and doubled at each subsequent level of the image resolution, so that we have N = 256 at the last two grid levels (maximal image resolution). The last resolution level was dedicated to segmentation as the label set was limited to the segmentation labels. Such a high resolution was not necessary for increased quality of the registration. The smoothing prior weights λ s and λ r were set to 1 and 2 respectively. The projection function was defined as cubic B-splines.

Using those parameters, the algorithm runs for about 4 minutes.

A.5.2 Results

We present preliminary results on a pair of preoperative/intraoperative T2 images. The preoperative tumor, the resection cavity and residual tumor were manually segmented. Both volumes were skullstripped. In order to ensure high correspondences and overlap during the initial estimation of the joint histogram, the preoperative image was affinely registered to the intraoperative image prior to the non-rigid registration using block matching. The preoperative tumor, residual tumor and resection cavity were manually segmented. We compare the automatic resection cavity segmentation obtained with the manual segmentation, and the deformed preoperative tumor to the residual tumor plus resection cavity. Visually, we observe high correspondences between the segmented regions as shown in Fig.

A.6 Discussion

This chapter illustrates the adaptability and modularity of the Joint registration and segmentation method introduced in the previous chapter. A probabilistic interpretation of the method is proposed towards registration of preoperative and intraoperative images with tumor resection. The resection cavity is detected through modeling the intensity profile in the resection area as a mixture of Gaussians while a probabilistic similarity criterion for registration is designed based on the images' joint intensity histogram. The method registers the preoperative image to the intraoperative image while segmenting the resection cavity and preserving the surrounding tumor area. This scheme enables direct comparison of the preoperative surgery plan with the current surgery status and gives a direct estimate of the amount of tumor removed and the residual tumor. The speed of the method makes it easily exploitable in a surgery setting, especially given that many preprocessing steps (skullstripping, affine registration) could be carried out preoperatively or at the beginning of the surgery.

Being tailored for the specific application of tumor resection, the proposed method is less modular in the sense that the matching criterion has to be probabilistic. It is however more adapted to multimodal registration due to the use of the information theoretic based joint histogram as similarity criterion. Given the limited dataset, a training step was not included and the GMM was used to directly model the image's intensity profile and extract Gaussians corresponding to the resection cavity. Learning a GMM directly on a training set of resection cavity voxels would certainly yield 136APPENDIX A. CONCURRENT SEGMENTATION AND REGISTRATION APPLICATION more accurate results, especially in the areas filled with air. Additionally, the method would have to be tested on a larger dataset to thoroughly assess its potential and the possible improvements.

1. 1

 1 Preoperative functional mapping. The functional areas are reorganized around the tumor's boundaries. [Duffau 2006] . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Intraoperative brain mapping. (a) Preoperative MRI, (b) Intraoperative Mapping.

2. 1

 1 Examples of different DLGG appearances on FLAIR T2 MRIs. . . . . . . . . . . . 2.2 Impact of an magnetic field B 0 on the spins' magnetization. (a) Random spins orientation, (b) alignment with the external field, (c) spin precession around B 0 . Images taken from [Puddephat 2002] . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Impact of TE and TR on the NMR signal [Buxton 2009]. (a) NMR signal, amplitude dependent on T2, (b) Impact of TR on the initial amplitude of the signal, (c) longitudinal magnetization dependent on TR and T1. . . . . . . . . . . . . . . . . 2.4 Contrast enhancement of High grade gliomas. (a) T1 weighted image, (b) contrast enhanced T1 weighted image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Images sequences for the visualization of a Low grade glioma. (a) T1, (b) T1 contrast enhanced (no enhancement of the tumor), (c) T2, (d) FLAIR images. . . . 2.6 Level set of an implicit function u and corresponding curve C [Leventon 2000]. . . 2.7 Figure extracted from [Droske 2001]. Adaptive grid sampling, progressively constructed as the front progresses. The colors indicates the arrival time of the propatation (from blue (early) to red) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Illustration of the kNN clustering approach. We aim at classifying the black sample as either a blue square or red circle. Setting k=5 selects the circle class, while k=8 (dashed circle) selects the other class. . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Illustration of the SVM approach. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Illustration of the SWA algorithm extracted from [Corso 2008]. . . . . . . . . . . . 13 2.11 Evolution of the joint histogram between an image and iself as one image is progressively rotated (0,2, 5 and 10 degrees). Below the histogram are the joint entropy values. Example taken from [Pluim 2003] . . . . . . . . . . . . . . . . . . . . . . 45 2.12 Different types of transformation models. (a) original image, (b) rigid transformation, (c) affine transformation, (d) non-linear transformation. . . . . . . . . . . . . 46 2.13 Schematic illustrating the AC-PC line and axes of the Talairach space. Example taken from http://airto.ccn.ucla.edu/. . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.14 MNI-ICBM-152 probabilistic atlas. From left to right: T1, T2 and Proton Density weighted images ; Gray Matter, White Matter and CSF probability maps. Figure extracted from [Fonov 2011] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.15 Modified MNI-305 atlas for integration of tumor and edema probability [Prastawa 2003]. Top row: healthy probability (White Matter, Grey Matter, CSF). Bottom row: tumor(left) and edema (right) probabilities. . . . . . . . . . . . . . . . . . . . . . . 52 3.1 Illustration of the block matching process. The floating image is on the left and target on the right.Figure extracted from [Ourselin 2000] . . . . . . . . . . . . . . 58 3.2 Visualization of the complete graph superimposed to the mean registered image. For visibility reasons, edges of weight larger than 35 are not shown. Short distances appear red and thick while large distance values are thin and blue. . . . . . . . . . 59 3.3 Probability maps of three different clusters superimposed to the mean registered image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.4 Visual examples of features used for boosting training: (a) median filter, (b) standard deviation, (c) entropy, (d) symmetry, (e,f) Gabor features. . . . . . . . . . . . 68 3.5 Image courtesy of [Komodakis 2008b]. (a) Primal-Dual principle: the optimal cost always lies between b T y and c T x for any pair (x, y) of feasible solutions. If both values are close, the approximated solution to the primal problem is close to the optimal solution too. (b) Primal-Dual schema: iteratively generate pairs (x k , y k ) of feasible solution until they are close enough (their ratio is below f * ) to obtain an f * -approximation of the solution . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.6 Distribution of tumors sizes in the data set (number of tumors vs tumor size in cm 3 ) 73 3.7 Visualization of the complete clustered graph superimposed to the mean registered image. The numbers correspond to the number of nodes in each cluster. . . . . . . 74 3.8 Schematics of the brain structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.9 Cluster validity indices with respect to the value of α (a,b,c,d) and the number of clusters (e). Dunn index (a), Davies Bouldin index (b), Silhouette index (c) and combined indices (d,e). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.10 Errorbars of the indices computed for the cross validation experiments with respect to the value of α (a) and the number of clusters (b). . . . . . . . . . . . . . . . . . 77 3.11 Illustration of the graph matching results. (a) Complete match, (b,c) Partial match. Positive matches correspond to blue edges and mismatched samples to red edges. . 78 3.12 Boxplots of the Dice score (a), True Positive rate (b), False positive rate (c), and MAD score (d) between the automatic and manual tumor segmentation for the three different methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 3D visualization of an example segmentation result. (a) manual segmentation, (b) automatic segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Positions of the cluster centers on the ICBM atlas after registration. A strong symmetry can be observed between clusters in columns (a) and (b) of the same image row. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Visual Segmentation results. (a) boosting score, (b) Boosting classification (thresholding), (c) Pairwise MRF, (d) MRF with spatial prior. . . . . . . . . . . . . . . . 4.1 Overview of the Joint Registration and Segmentation method . . . . . . . . . . . . 4.2 Summary of different steps of the ORBIT algorithm. Figure extracted from [Zacharaki 2008]. 90 4.3 Modified ICBM probabilistic atlas by the tumor growth model. From left to right: tumor, edema, CSF, gray matter and white matter probabilities. Figure extracted from [Gooya 2012] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Registration uncertainty and displacement set resampling for one control point: (a) Min marginal values per displacement label (blue: low, red: high energy) associated covariance matrix centered at the optimal label, (b) Min marginals visualization on a 2D slice, (c) Original isotropic displacement set and (d) Uncertainty driven displacement set, following the brain boundaries. . . . . . . . . . . . . . . 4.5 Visual representation of the grid refinement from level j (left) to level j+1 (right). Grid resampling: the nodes that have direct correspondences appear in white, and the new nodes and edges are red. The edges connecting the 2 grids represent the nodes' neighborhood. The grid is shown in 2D for increased visibility. . . . . . . . 4.6 The 4 levels of the incremental approach. . . . . . . . . . . . . . . . . . . . . . . 4.7 The two displacement sampling strategies, centered around the node's initial position: (a) Dense displacement sampling, (b) Sparse displacement sampling. The red dots correspond to the possible displacements. . . . . . . . . . . . . . . . . . . . . 4.8 Mean percentage of activated nodes per resolution level. . . . . . . . . . . . . . . 4.9 Visual examples of the activated nodes for the 4 levels of the incremental approach. The nodes are superimposed to the target image. . . . . . . . . . . . . . . . . . . . 4.10 Quantitative Segmentation Results: Error bars (mean and standard deviation) of the Dice score, False Positive (FP) and True Positive (TP) rates and MAD score for the joint framework with low (JSRLow) and high resolution (JSRHigh), the MAP-MRF segmentation framework (SegMRF) and the uncertainty based approach (Ucy).111 4.11 Quantitative Registration Results: Error bar graphs of the Dice, True Positives (TP), False Positives (TP) and MAD scores obtained for the joint framework with low (JSRLow) and high resolution (JSRHigh), the individual registration framework with masked pathology (RegMask) and the uncertainty based approach (Ucy). 4.12 Visual examples of Registration results. (a) Target image, deformed image obtained using (b) the individual framework with masked pathology, (c) the Joint framework, high resolution, (d) the Joint framework, with adaptive sampling. . . . 4.13 Visual examples of Segmentation results. (a) Boosting Score, (b) MAP-MRF framework, (c) Joint framework, high resolution, (d) Joint framework, with adaptive sampling. Automatic segmentations (blue) are compared to the manual segmentation (red) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Gaussian Mixture Model of the resection cavity intensity. (a) Intensity histogram, (b) Modeled intensity distribution, (c) Gaussians constituting the GMM. . . . . . . A.2 Log probabilities associated to the resection cavity coupled with the spatial prior. . A.3 Difference Image after affine registration (a) and after non-rigid registration (b). . . A.4 Visual registration and segmentation results. (a) Preoperative image after affine registration, (b) non-rigidly deformed preoperative image, (c) intraoperative image segmented. Resection cavity segmentation: automatic (blue) and manual (red); Residual tumor segmentation: manual (magenta) and registered preoperative tumor (yellow). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Tables 3.1 Anatomical localization of the different clusters . . . . . . . . . . . . . . . . . . . 75 3.2 Mean and median (in parentheses) values of the different scores and methods for tumor segmentation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

Figure 1 . 2 :

 12 Figure 1.2: Intraoperative brain mapping. (a) Preoperative MRI, (b) Intraoperative Mapping.Each tag correspond to a functional area that is identified through cortical stimulation, (c) Postoperative MRI. Images extracted from[Duffau 2006]

  Figure 1.3: Modeling of the LGG's growth rate. (a) Linear growth rates of 27 patients, displaying a mean growth rate of 4mm/year. (b) Growth rate evolution of a patient after anaplastic transformation. Figures extracted from[Mandonnet 2003] and[Mandonnet 2008] respectively.

  Figure 2.1: Examples of different DLGG appearances on FLAIR T2 MRIs.

  Figure 2.2: Impact of an magnetic field B 0 on the spins' magnetization. (a) Random spins orientation, (b) alignment with the external field, (c) spin precession around B 0 . Images taken from [Puddephat 2002]

  Figure 2.3: Impact of TE and TR on the NMR signal [Buxton 2009]. (a) NMR signal, amplitude dependent on T2, (b) Impact of TR on the initial amplitude of the signal, (c) longitudinal magnetization dependent on TR and T1.

  Figure 2.4: Contrast enhancement of High grade gliomas. (a) T1 weighted image, (b) contrast enhanced T1 weighted image.

  Figure 2.5: Images sequences for the visualization of a Low grade glioma. (a) T1, (b) T1 contrast enhanced (no enhancement of the tumor), (c) T2, (d) FLAIR images.

Figure 2 . 6 :

 26 Figure 2.6: Level set of an implicit function u and corresponding curve C [Leventon 2000].

Figure 2 . 7 :

 27 Figure 2.7: Figure extracted from [Droske 2001]. Adaptive grid sampling, progressively constructed as the front progresses. The colors indicates the arrival time of the propatation (from blue (early) to red)

Figure 2 . 8 :

 28 Figure 2.8: Illustration of the kNN clustering approach. We aim at classifying the black sample as either a blue square or red circle. Setting k=5 selects the circle class, while k=8 (dashed circle) selects the other class.

Figure 2 .

 2 Figure 2.10: Illustration of the SWA algorithm extracted from [Corso 2008].

Figure 2 .

 2 Figure 2.11: Evolution of the joint histogram between an image and iself as one image is progressively rotated (0,2, 5 and 10 degrees). Below the histogram are the joint entropy values. Example taken from [Pluim 2003]

Figure 2

 2 Figure 2.13: Schematic illustrating the AC-PC line and axes of the Talairach space. Example taken from http://airto.ccn.ucla.edu/.

Figure 2

 2 Figure 2.14: MNI-ICBM-152 probabilistic atlas. From left to right: T1, T2 and Proton Density weighted images ; Gray Matter, White Matter and CSF probability maps. Figure extracted from [Fonov 2011]

Figure 2 .

 2 Figure 2.15: Modified MNI-305 atlas for integration of tumor and edema probability [Prastawa 2003]. Top row: healthy probability (White Matter, Grey Matter, CSF). Bottom row: tumor(left) and edema (right) probabilities.

Figure 3 . 1 :

 31 Figure 3.1: Illustration of the block matching process. The floating image is on the left and target on the right. Figure extracted from [Ourselin 2000]

  Figure 3.2: Visualization of the complete graph superimposed to the mean registered image. For visibility reasons, edges of weight larger than 35 are not shown. Short distances appear red and thick while large distance values are thin and blue.

Figure 3 . 3 :

 33 Figure 3.3: Probability maps of three different clusters superimposed to the mean registered image.

  [3.4] while examples of boosting classification results are shown in Fig.[3.15a]

  Figure 3.4: Visual examples of features used for boosting training: (a) median filter, (b) standard deviation, (c) entropy, (d) symmetry, (e,f) Gabor features.

  Figure 3.5: Image courtesy of [Komodakis 2008b]. (a) Primal-Dual principle: the optimal cost always lies between b T y and c T x for any pair (x, y) of feasible solutions. If both values are close, the approximated solution to the primal problem is close to the optimal solution too. (b) Primal-Dual schema: iteratively generate pairs (x k , y k ) of feasible solution until they are close enough (their ratio is below f * ) to obtain an f * -approximation of the solution

  Figure 3.6: Distribution of tumors sizes in the data set (number of tumors vs tumor size in cm 3 )

  Figure 3.7: Visualization of the complete clustered graph superimposed to the mean registered image. The numbers correspond to the number of nodes in each cluster.

Fig. [ 3

 3 Fig.[3.14], highlighting a strong symmetry between the hemispheres. We observe, as summarized in table3.1, a strong predominance of tumor appearances in the Insula (more than a third of the tumors) and towards the frontal lobe, while there are no clusters in the occipital lobe. Schematics of the brain structure are shown in Fig.[3.8]

  Figure 3.9: Cluster validity indices with respect to the value of α (a,b,c,d) and the number of clusters (e). Dunn index (a), Davies Bouldin index (b), Silhouette index (c) and combined indices (d,e).

  Figure 3.10: Errorbars of the indices computed for the cross validation experiments with respect to the value of α (a) and the number of clusters (b).

  [3.11] that highlights the proximity of mismatched centers. The worst case with 4 non matched samples is shown in Fig.[3.11c].

  Figure 3.11: Illustration of the graph matching results. (a) Complete match, (b,c) Partial match. Positive matches correspond to blue edges and mismatched samples to red edges.

  Figure 3.12: Boxplots of the Dice score (a), True Positive rate (b), False positive rate (c), and MAD score (d) between the automatic and manual tumor segmentation for the three different methods.

  Figure 3.13: 3D visualization of an example segmentation result. (a) manual segmentation, (b) automatic segmentation

  Figure 3.14: Positions of the cluster centers on the ICBM atlas after registration. A strong symmetry can be observed between clusters in columns (a) and (b) of the same image row.

  Figure 4.1: Overview of the Joint Registration and Segmentation method

  [4.1]

Figure 4 . 2 :

 42 Figure 4.2: Summary of different steps of the ORBIT algorithm. Figure extracted from [Zacharaki 2008].

Figure 4 . 3 :

 43 Figure 4.3: Modified ICBM probabilistic atlas by the tumor growth model. From left to right: tumor, edema, CSF, gray matter and white matter probabilities. Figure extracted from [Gooya 2012]

  s lp )Sim I(x), A(x + d lp )

  Figure 4.4: Registration uncertainty and displacement set resampling for one control point: (a) Min marginal values per displacement label (blue: low, red: high energy) associated covariance matrix centered at the optimal label, (b) Min marginals visualization on a 2D slice, (c) Original isotropic displacement set and (d) Uncertainty driven displacement set, following the brain boundaries.

Figure 4 . 5 :

 45 Figure 4.5: Visual representation of the grid refinement from level j (left) to level j+1 (right). Grid resampling: the nodes that have direct correspondences appear in white, and the new nodes and edges are red. The edges connecting the 2 grids represent the nodes' neighborhood. The grid is shown in 2D for increased visibility.

Figure 4 . 6 :

 46 Figure 4.6: The 4 levels of the incremental approach.

  Figure 4.7: The two displacement sampling strategies, centered around the node's initial position: (a) Dense displacement sampling, (b) Sparse displacement sampling. The red dots correspond to the possible displacements.

Figure 4 . 8 :

 48 Figure 4.8: Mean percentage of activated nodes per resolution level.

  [4.11], while qualitative results are shown in Fig.[4.12]. We observe higher performance around the tumor area where the individual framework can fail.

Figure 4 . 9 :

 49 Figure 4.9: Visual examples of the activated nodes for the 4 levels of the incremental approach. The nodes are superimposed to the target image.

Figure 4 .

 4 Figure 4.11: Quantitative Registration Results: Error bar graphs of the Dice, True Positives (TP), False Positives (TP) and MAD scores obtained for the joint framework with low (JSRLow) and high resolution (JSRHigh), the individual registration framework with masked pathology (RegMask) and the uncertainty based approach (Ucy).

Figure 4

 4 Figure 4.12: Visual examples of Registration results. (a) Target image, deformed image obtained using (b) the individual framework with masked pathology, (c) the Joint framework, high resolution, (d) the Joint framework, with adaptive sampling.

Figure 4

 4 Figure 4.13: Visual examples of Segmentation results. (a) Boosting Score, (b) MAP-MRF framework, (c) Joint framework, high resolution, (d) Joint framework, with adaptive sampling. Automatic segmentations (blue) are compared to the manual segmentation (red)

  p(S = l|I, J, T k ) (log p(I(T (x))|J, S = l, T ) + log p(J(x)|S(x) = l)) + log p(T ) (A.3)

  (x n |µ k , Σ k ) (A.11) A.4.1.

  Figure A.1: Gaussian Mixture Model of the resection cavity intensity. (a) Intensity histogram, (b) Modeled intensity distribution, (c) Gaussians constituting the GMM.

Figure A. 2 :=

 2 Figure A.2: Log probabilities associated to the resection cavity coupled with the spatial prior.

P

  (I(x), J(x)|l p , s lp = 1) = 1 2 N (J(x+d lp )|µ min , Σ min )+ 1 2 N (J(x+d lp )|µ max , Σ max ) (A.20)GMM modeling is illustrated in Fig.[A.1] and the corresponding probability score is shown in Fig.[A.2]

  [A.4]. Quantitatively, we obtain a dice score of 85 % for residual tumor and 80 % for resection cavity. Strong visual correspondences are observed between the registered preoperative image and the intraoperative image as shown in Fig.[A.3].

  Figure A.3: Difference Image after affine registration (a) and after non-rigid registration (b).

  Figure A.4: Visual registration and segmentation results. (a) Preoperative image after affine registration, (b) non-rigidly deformed preoperative image, (c) intraoperative image segmented. Resection cavity segmentation: automatic (blue) and manual (red); Residual tumor segmentation: manual (magenta) and registered preoperative tumor (yellow).

Table 3 . 2 :

 32 Mean and median (in parentheses) values of the different scores and methods for tumor segmentation.

	Method	Dice	True Positive False Positive MAD
	Boosting	50 (55) % 77 (83) %	59 (58) %	22 (22) mm
	MAP-MRF 66 (74) % 70 (78) %	29 (21) %	9.5 (6) mm
	Atlas Prior 69 (77) % 66 (74) %	20 (12) %	6 (3.5) mm

A.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

The aforementioned arguments are also supported from experimental evidence. Different distance metrics have been used, including for example the Mahalanobis distance between tumors. The distance between the tumors' centers of gravity has been proven to be more robust.

Acknowledgments

each learned clustered graph obtained by cross validation, we evaluated the impact of the 210th sample on the structure of the graph and the quality of the clusters. In order to determine how well the graph represents the data, we tried to assign the 210th sample to a cluster based on two measures: the distance from the center of the cluster and the silhouette index. To assign a sample to a cluster, three criteria have to be verified:

Publications of the Author

International Journals The treatment of DLGGs has long been controversial. Based on several retrospective studies, the positive impact of surgical resection on the patient's progression free and overall survival is undeniable and maximal safe tumor resection is now considered as the first therapeutic option for DLGGs. As a result, the surgical resection task is very frequent and an essential part of the treatment of the tumor. The process is very challenging as the surgeon seeks to remove as much tumor as possible while preserving the critical tissue associated to functional areas. Medical images are heavily used for surgery planning. Structural imaging enables to identify the position and extent of the tumor, while functional images allows mapping of the functional areas with respect to the tumor. Such data allows to estimate the volume of possible resection. During surgery, the preoperative data is often used in order to guide the surgeon and provide additional information on the brain's functional and structural organization. However, the preoperative image is generally aligned using rigid or affine registration which could be highly unreliable due to important deformations of the brain during surgery (mainly due to brain shift and tissue collapse). Those deformations are due to the resection as well as a brain shift during surgery caused notably by gravity, tumor mass effet, edema and CSF leakage [Hartkens 2003].

The development of intra operative MR imaging enabled to obtain anatomical data and increased the efficiency of the procedure through reliable anatomical guiding. However, surgical planning (mainly, the expected amount of tumor resection) as well as additional functional and structural information (fMRI, DTI for instance) acquired prior to the surgery. Such modalities l p ∈ L = {1, .., 2n}. Each label corresponds to a pair segmentation/displacement {d lp , s lp } where the displacements are selected among a predefined discrete set ∆ = {d 1 , ..., d n }. Once the grid nodes' labels are obtained, the transformation and segmentation fields are computed on the whole image by interpolation:

Where η(.) is the projection function.

The problem can be cast as a MAP-MRF estimation of the posterior probability P (l|I, J) on the labeling l = {l p } p∈V given the observed images : l = argmax l∈L P (l|I, J) ∝ P (I, J|l)P (l) (A.5)

The prior probability P (l) is modeled as a pairwise MRF defined on the sparse grid G:

Following the methodology of chapter 4, the pairwise term is defined as the sum of a registration and segmentation smoothing prior weighted by the distance between the nodes:

The unary term V p (l p ) is a spatial prior that exploits the fact that the localization of the tumor in the preoperative image gives insight on where the resection cavity is going to be. Indeed, the resection area is expected to be in the preoperative tumor area or within its close vicinity. The preoperative tumor segmentation is used to define the area where the resection cavity is expected to be and not to be. To allow segmentation outside the tumor area, the segmentation map is smoothed by Gaussian filtering.

where S I,s is the smoothed tumor segmentation map, and η(.) is the back projection function that describes the influence of voxel x on control point p. We recall that the transformation T (x) is approximated in the unary terms as x + d lp in order to preserve the independency asumption inherent to MRF models.

130APPENDIX A. CONCURRENT SEGMENTATION AND REGISTRATION APPLICATION

Assuming independency between the observations and the projection of the likelihood probability on the grid control points leads to a MAP estimation problem that aims at minimizing the following objective function, defined on the sparse grid G:

V p , q(l p , l q ) (A.9)

The obtained MRF energy is minimized using the Fast-PD algorithm.

A.4 Likelihood term

The next step consists in defining the likelihood term P (I(x), J(x)|l p ) that will play an essential part in the concurrent registration and segmentation framework. Two cases have to be distinguished: if we are in the background (s lp = 0) we want the joint probability to be maximized when the images are aligned while in the resection cavity area (s lp = 1) no correspondences should be sought. The likelihood must enable the detection of the resection cavity as well as the matching of both images.

A.4.1 Resection Cavity Detection

The intensity profile associated to a resection cavity is simpler than the one associated to brain tumors. The detection is therefore easier and simple classification methods can be used in order to separate it from the rest of the image. The resection cavity likelihood term is defined based on the observation that it is made of two possible elements: air and fluid (CSF or blood). When working on T2 images, the air appears as the darkest intensity class while fluid is the brightest. The resection area can be isolated from the rest of the image by modeling the intensity of the area delineated by the spatial prior (that includes the resected area) as a mixture of Gaussians and selecting the ones that correspond to its intensity profile.

A.4.1.1 Gaussian Mixture Models

The Gaussian Mixture Model (GMM) is a popular classification technique that models the intensity distribution as a mixture Gaussian distributions. By choosing a sufficient number of Gaussians, it is possible to capture the intensity distribution of different image classes. The Gaussian mixture distribution is written as follows:

where K is the number of Gaussians, and (µ k , Σ k ) are their respective mean and standard deviation. π k is a mixture coefficient describing the influence of the Gaussian on the overall distribution. The mixture coefficients must satisfy 0 ≤ π k ≤ 1 and K k=1 π k = 1.