New Expectation Maximization segmentation pipeline in Slicer 3

Many neuroanatomy studies rely on brain tissues segmentation of magnetic resonance (MR) images. Expectation-maximization (EM) algorithm are very popular framework for this task. The Surgical Planning Laboratory (SPL) of Harvard Medical School developped its own segmentation algorithm. Nevertheless, the segmentation is not accurate if the image exhibits intensity inhomo-geneity. Moreover the optimum parameters are challenging to estimate. In this context, this document ﬁrst proposes a full description of the EM algorithms. Then solutions are presented to enhance the segmentation. Contributions are very various, from the algorithm to the graphic user interface (GUI). We add a bias ﬁeld correction to correct the intensity inhomogeinities in the MR image, as a preprocessing step in the algorithm. We also propose new procedures in the graphic user interface to deﬁne the distribution of the tissues to be segmented. Some tools to make the segmentation process easier and more accurate are also presented. Finally, the results obtain with the new segmentation pipeline are evaluated by an expert.


Acknowledgment
First of all, I would like to thank Ron KIKINIS who gave me the opportunity to carry out this internship at the Surgical Planning Laboratory of Harvard Medical School.He has always took the time discuss with me and to give me precious advices in my work.I would also like to thank Katie MASTROGIACOMO for her help on the administrative side, especially to obtain the visa.
Of course, my thoughts also go to Sylvain JAUME.We worked in close collaboration.He was very patient and comprehensive.Regarding my student status, he was very pedagogic to improve my knowledge in a lot of different ways and I am very grateful for that.He spent a lot of time to guide me and to help to achieve my internship.
I wish to express my gratitude to Steve PIEPER, Andriy FEDOROV and Daniel HAEHN for their assistance .They helped a lot me to get familiar with the whole working environment.I must also acknowledge Christine CAVARO-MENARD who has accepted to supervise this internship and Dominique MARATRAY for her assistance.
Finally, I would like to thank all the people working in the laboratory for their reception and their friendliness.Thanks to them, my time in the Surgical Planning Laboratory was a very enriching experience and it was a real pleasure to work with them.

Context and motivation
Nowadays, medical image processing is becoming a major field of research in most of the laboratories.Indeed, because of the increasing complexity of the data they have to deal with, physicists need something to help.Help must be provided in many different ways.Before the surgery, to etablish a fast and accurate diagnosis.During the surgery to prevent physicians from errors and to help them to proceed to more precise moves.After the surgery, to see if it succeeded, or to follow the pathology of a patient.The informations brought to the physicians by the tools must be accurate, robust and provide a fast feedback.
In this context of pre and post operation, plenty of work has already been done.Nevertheless, there is still a lot of work to achieve.Regarding data storage and exchange, the increasing among of information leads us to find more approriate methods for the same purpose.Another interesting contribution of computer science to medicine is image segmentation.New methods have to be developped for a better diagnosis, or to detect new pathologies.A lot of segmentation techniques appeared like level-set segmentation, region growing or texture based segmentation.Each one is adapted to a specific problem like vessel segmentation or tumor detection.Another remarkable contribution is the segmentation based on expectation maximization (EM).It is very well suited for brain tissue segmentation.
For segmentation of brain MR images, the Surgical Planning Laboratory (SPL) at the Brigham and Women's Hospital, affiliated to Harvard Medical School, has developped an EM algorithm for segmentation.This model assumes that the MR volumes does not exhibit an intensity inhomogeinity.Moreover, the implementation is not widely used so far, because the complexity of the segmentation process.Finally, some steps have to be enhanced.In this report we will present an approach to enhance the segmentation for inhomogeneous MR images, correcting the intensity bias field.We will also provide tools to the end-user for an easier and more accurate segmentation process.

Contents
The main body of this report is divided as follows.Chapter 2 focuses on the EM segmentation.Background is reminded and the algorithm used in the SPL is fully described.We also present in this chapter the segmentation pipeline developped in the SPL and the limitations of the current implementation.Chapter 3 describes our contributions.It explains the solutions we proposed to enhance the segmentation and to improve the usability of the current framework.Chapter4 shows the results obtained.It aslo discuss about the results obtained.An expert radiologist evaluates the segmentation.Finally, some perspectives about the future work to be done are presented.

Expectation-maximization applied to brain segmentation
Here we get going with theory of the expectation-maximization (EM) applied to brain segmentation.We show first a simple approach of the problem, in the particular case of Gaussian mixture models followed by a more generalistic approach.The simple approach will give the reader an intuitive understanding of the problem then the general approach will formalize it in order to adapt it to most of the segmentation problems.Finally, there will be a presentation of the algorithm implemented in the software suit Slicer 3 [22].

Presentation of the EM segmentation
The EM algorithm was originally described in 1977 by Arthur Dempster, Nan Laird, and Donald Rubin [1].They generalized and developped a method used previously by authors, for specific applications.It is widely used to solve problems where data are "missing".The EM algorithm is an iterative algorithm which works in two steps: Expectation and Maximization.It can be use to solve a lot of image processing's problems like classification, restoration [3] and motion estimation [2].The EM segmentation has been applied to segmentation Nowadays, EM algorithms are become a popular tool for classification problems.It is particulary well suited for brain MR images segmentation.A lot of algorithms already exist.They present complex frameworks using spatial information, neighborhood or intensity inhomogeneities to enhance the classification.
In the SPL, the algorithm developped uses spatial, tree structure arborescence and information about intensity inhomogeneity to segment the brain.

Background
We start with a presentation of all the fundamentals the reader needs to have a good understanding of EM segmentation.We begin with a description of the statistical model used for the brain.Then we describe briefly the widely used GM model.Finally we will present what the mamixum likelihood function is.This part is mainly based on [4], [5] and [6]

Statistical model used for the brain
We define the voxel intensities of a volume as Y = {y 1 , ..., y n } when the image counts n voxels.Y is called observed data because this is the the data we see when we observe the image.Each y is a realization of the random variable Y .The actual labelling of the image is Z = {z 1 , ..., z n }.Z is called hidden data because we do not know the value of each label.This is precisely the purpose of the segmentation: estimating the hidden data from the observed data.We assume that the observed data is generated from the hidden data and a parameter Φ.The parameter Φ can either be a probability density function, noise,or bias field, among others, depending on the model.
Y and Z can be viewed as n-dimensional random variables Y = {Y 1 , ..., Y n } and Z = {Z 1 , ..., Z n } then each y i is a realization of Y i and each z i is a realization of Z i .The conditional probability function describing The easiest model assumes that all intensities in one class are the same, but these intensities are corrupted by factors like noise modelled by a Gaussian distribution.We can describe the relationship as below: where µ k is the mean intensity of the k th tissue and n i a random sample generated by the corrupting factor(s).Let's say that n i is generated by a Gaussian probability distribution function G(., 0, σ), with 0 mean and σ variance.That means that y i is a random sample generated by a Gaussian probability density function G(., µ k , σ).Let's assume that each class has a different variance, G(., µ k , σ) becomes G(., µ k , σ k ) and it leads to: As the labelling is not known, it is usefull to express the probability density function (PDF) of Y i only depending on parameter Φ with the total probability theorem: p(Z i = k|Φ) is the prior probability.It expresses the probability that a voxel i belongs to a class k. p(y i |Z i = k, Φ) is the likelihood.In our case, we will assume that the prior probability is unique for each each label.The new model we obtain for the labelling is a widely used one: the gaussian mixture model.

Gaussian mixture model
Let's remind the first hypotesis: the conditional probability function for each tissue to segment is defined as in Equation 2.1.Moreover, we will assume that prior probability is a constant c k for each class k. c k is the weight of the class k.
The last assumption will be that Φ contains unknown means, variances and weights for each tissue.Then we can express Φ as a set of parameters such as Φ Using Equations 2.1 and 2.3, Equation 2.2 becomes: In the case of GM model, each voxel is considered to be independent.That means that each voxel will have its own probability density function.Consequently, the normalized histogram of the whole volume can be interpreted as an approximation of the sum of all the probability density functions.

Maximum likelihood
In our case, we know the intensity of each observed voxel y i .Φ are to be found.The best estimation of Φ will be obtain using the maximum likelihood principle.p(y i |Φ) is called likelihood function.For each value of Φ it returns the value of the likelihood of y i , given Φ.
The voxels are considered to be independent.Through the whole volume, it leads to: The objective is to find the parameter Φ which will maximize the likelihood of the observed volume.We can note this parameter: Therefore, it is more convenient to work with logarithm because the product from Equation 2.5 will be converted into a summation.Equation 2.6 becomes: With Equation 2.2, let us denote: Finally, in case of GM model, with Equation 2.4, L(Φ) becomes: The log likelihood can be maximized by finding partial derivatives for each parameter.When the partial derivative of L(Φ) is 0 for a parameter, we have the maximum likelihood for the parameter Φ.
For example, to find the maximum likelihood regarding µ k we first have to find when: Then we compute the partial derivative of L(Φ) over µ k : Using Bayes' theorem (see App. A, Sec.A.2), we notice that: Thus, setting the denominator to 0 in Equation 2.7 and using Equation 2.8 yields: Let us denote Equation 2.9 leads to: Proceeding the same way as we did for Equation 2.9, we can get similar equations for variance σ k and weight c k .We find that: (2.12) Equations 2.11,2.12 and 2.13 provides us This equation expresses that a voxel i belongs to the class k.It is formulation for soft segmentation.p ik is called soft assignment.p ik is used to fill a "map" of soft segmentation.At the end of the segmentation process, this map contains the probability that the voxel i belongs to class 1, 2, ..., K.We determine the class of the voxel i looking at the class which has the highest probability in the map for this given voxel i.
The segmentation can now be calculated following an iterative process called expectation maximization.

Expectation maximization algorithm
The EM algorithm is a method to find the maximum likelihood for a given set of parameter (Φ in our case).Here we start with an intuitive description of the algorithm in the particular case of Gaussian mixture model then we will present a more general definition.

Algorithm in case of Gaussian mixture data model
Let's assume that we can find the maximum likelihood of the hidden data by a direct differentiation (because of the Gaussian mixture model).The EM algorithm is an iterative process of two steps: the expectation step (E-Step) and the maximization step (M-Step).At each iteration, the maximum likelihood will be increased until convergence is reached.

• E-step
In this step, we calculate an estimation of soft segmentation p (m+1) with Equation 2.14 as below.We know all the variables needed for the calculation from the observed data and the current parameter estimate Φ (m) .Note that an initialization is necessary for the first iteration.
In this step, we estimate the maximum likelihood for parameter Φ (m+1) .We do it with Equations 2.11,2.12 and 2.13 as below.We know all the variables needed for the calculation from the observed data and the current estimate p (m+1) of hidden data.
The problem is simple as long as we are working with Gaussian mixture model.In the other case, the log-likelihood can not be maximized by direct differentiation and a generalized approach must be used.

Generalized algorithm
Now we assume that we are no longer working with Gaussian mixture model.Thus, we must use a more general algorithm.To explain the general algorithm, we will start from the log-likelihood L(Φ).As presented in the previous subsection, we wish to find Φ such as p(Y |Φ) is a maximum.We have: Since log is a strictly increasing function, the value of Φ which maximizes p(Y |Φ) also maximizes L(Φ).Assume that after the m th iteration, the current estimate of Φ is given by Φ (m) .Since the objective is to maximize L(Φ) we want an updated estimate Φ such that In other words, we want to maximize the difference L(Φ) − L(Φ (m) ).For convenience, we introduce a new variable z ik which means that Z i = k.Using the new notation, we can transform this difference as below: We can deduce this inequality from Jensen's inequality (see App. A, Sec.A.3) since p(z ik |y i , Φ (m) ) is a probability mesure and log a concave function ( [5]).
We will then use the fact that k p(z ).This allows us to bring log p(Y |Φ (m) ) into the summation.
We will also use a new variable e k to express the difference through the whole volume without a summation.e k is defined as e k = {z 1k , ..., z nk } when the image consists of n voxels.For example, We can then conclude that: We have now a function l(Φ|Φ (m) ) which is bounded above by L(Φ).Additionnaly, we observe that (l(Φ (m) , Φ (m) )) = L(Φ (m) ).
Our objective is to choose values of Φ which will maximize L(Φ).We have shown that the function l(Φ, Φ (m) ) is bounded above by the likelihooh function L(Φ).Moreover, the value of Φ for which the function (l(Φ, Φ (m) )) and L(Φ) are egals is Φ = Φ (m) .Therefore, any Φ which will increase l(Φ, Φ (m) ) will increase L(Φ).In order to maximize L(Φ) as much as possible, we maximize L(Φ) such as: )....As soon as L(Φ(m + 1)) < L(Φ(m)), the convergence has been reached and the maximum likelihood found.We can formalize this search.Φ (m+1) is the updated value which is found after maximization of Φ using Φ (m) .
We notice from Equation 2.15, that for a given voxel i: Now, both expectation and maximization steps are apparents.
• E-step This is the expectation step.During this step, we estimate the probability that the pixel i belongs to class k regarding Φ (m) .This equation is obtained from Equation 2.16 with Bayes formula.
Using this probability, E Z|Y,Φ (m) returns the expected value of the parameter Φ regarding Φ (m) .
• M-step This is the maximization step.During this step arg max Φ maximizes E Z|Y,Φ (m) for the parameter Φ.It returns Φ (m+1)   arg max The EM algorithm iterates until convergence is reached.The condition of convergence can differ from an algorithm to another.A possibility is to fix the number of iterations, C, of the algorithm.Most of the time, the following approach is used; convergence is reached when the difference between the estimation of the parameter Φ, at step m and step m + 1 is smaller than ε: If after the C th iteration, this condition is not satisfied, the EM algorithm is stopped.
We are now familiar with the theory behind the EM algorithm.The logic appears and we can summarize the basic algorithm with the

Expectation maximization algorithm used in Slicer 3
In this part, we present the EM algorithm which has been developped in Slicer 3 and the pipeline in which it has been integrated as well.Finally, we briefly discuss the limitations of this one.
The EM algorithm which is used in the SPL is derived from the original one.It enhances the original algorithm by the addition of informations like a probalistic atlas, a multichannel segmentation, a bias correction and a structure information.In Slicer 3, the Gaussian mixture model is used to describe the tissues to segment.Thus it will simplify the problem and the notations.

Probabilistic atlas
The EM algorithm is very sensitive to initialization since it only finds local extremums during the maximization step.A solution to enhance the initialization is to use atlases.An atlas is needed for each tissue the user wants to segment.For each voxel i of the volume, the atlas returns the probability that this voxel belongs to class k.This probability can be used as initialization.
From this value, it estimates Φ (1) and the algorithm iterates until convergence is reached.The probabilistic atlases are not only used to inialize the process.It is also use to get a more robust algorithm.Indeed, we can use the spatial information given by the atlases.Voxels will be classified not only based on intensity but regarding spatial position too.Van Lemput et al. ([8] and [9]) used the spatial prior at each iteration.It is constant and we then have a spatial information.
The probability that a pixel i belongs to class k, in the E-Step changes.From Equation 2.17, it becomes:

Multichannel segmentation
Most of the time, several sequences of a given modality are used to process brain segmentation.Indeed, the best suited sequence to use depends on the tissue you want to segment.For example, MR images provided by T1 sequences [12] are well suited to segment white matter (WM) but are not accurate for cerebrospinal fluid (CSF).On the contrary, T2 sequences [12] MR images which are well suited for CSF and not for WM.To formalize the utilization of different MR images sequence (i.e.different channels)s during the segmentation, we will change the definition of y i , µ k and σ k we did at the beginning.Let y i = {y i1 , y i2 , ..., y iR }, µ k = {µ k1 , µ k2 , ..., µ kR } and σ k = {σ k1 , σ k2 , ..., σ kR } when we use R channels, from different modalities to do the segmentation.The equations for the E-Step and the M-Step will remain the same.

Bias field correction
A major issue in MR modality is that the images can be corrupted by a low field bias field.It is mainly due to equipment limitations and to patient induced electrodynamic interactions ( [12]).We will now present how this bias field can be estimated and corrected in the EM algorithm.

Principle
Let I = (I 1 , ..., I n ) the observed intensities in an image, I * = (I * 1 , ..., I * n ) the ideal intensities and F = (F 1 , ..., F n ) the bias field.We use the assumption that the bias field is only multiplicative.Degradation of each voxel can then be expressed as: .., Y * n ) be the log-transformed observed and ideal intensities.B = (B 1 , ..., B n ) the log-bias field.This transform makes the bias field become additive instead of multiplicative without the log-approach.
We can model the PDF of the voxel intensity with a Gaussian distribution The low frequency characteristic of the bias field B can be modeled by a linear combination of smooth basic functions Ψ l (x) ( [13]).Let b i be the realization of the random variable pos(i) returns the 3D position (x, y, z) of the voxel i. a i is the i (th) value of the vector A = (a 1 , ..., a L ).A represents the bias field parameters.
In the Gaussian model, bias field can then be estimated using EM framework.The bias field parameter A will be used during the E-Step , through b i to estimate the soft segmentation.A will be re-estimated during the M-Step, after the maximization of the tissue class parameters (mean, variance and weight).Van Leemput formalised the two steps as below ( [8] and [9]): (2. 19) with: The bias field correction can be interpreted as follows: the estimated soft segmentation (p (m) ik ) and tissue class parameters are used to reconstruct the image Ỹ = ( ỹ1 , ..., ỹn ).This new image is supposed to be not corrupted by the bias field.We then substract the reconstructed image Ỹ from the observed image Y .We obtain the residual image R. From R, we estimate the bias field.F represents the discretized geometry of the bias field.W is an inverse covariance matrix.It returns information about the possible error for each voxel.The covariance matrix is described in details in [17].
The approach used in Slicer 3 is based on the same principle but differs regarding the maximization method and the parameter which is maximized.

Variation used in Slicer 3
In this method, we are working with a Gaussian mixture model.Moreover the parametres of this Gaussian distribution are assumed to be known.The idea of estimating the field in EM framework was originally proposed by Wells et al. ([10]).He proposed to only use maximization to re-estimate the bias field.The maximum a posteriori principle (MAP) instead of the maximum likelihood principle (Equation 2.6) is used to find the lower bound, during the maximization: As Bayes's theorem can be applied Proceeding the same way as we did in section 2.3, the new E-Step becomes: In Wells' method, the only parameter to estimate is the bias field.We assume that the noise has a Gaussian distribution: The equation for the bias field will change.We add the smoothness constraint in the Gaussian distribution: Σ −1 B .We also set F to a matrix filled with ones as no parametric model for the bias field is assumed.Finally, we define the mean residual image M (m+1) The Equation 2.19 for the bias field will then be remplaced by (B (m+1) ) T :

Hierarchical information
The last modification of the original EM algorithm is the addition of hierarchical information in the iterative segmentation process.The algorithm was described by Pohl et al. [11] The idea was to describe the structures we want to segment as a tree.It allows us to subdivide the segmentation process into subproblems, that are easier to solve.
Here we continue with an intuitive description of the process.It is a brief explanation of how the image would be segmented using the hierarchical information 2.2.At the first iteration, the MR image will be segmented into the background (BG) and the intracranial cavity (ICC) with the EM algorithm.At the second iteration, the BG will be segmented into the air and the skull.Finally, at the last iteration, the ICC will be segmented into white matter (WM), grey matter (GM) and cerebrospinal fluid (CSF).
To formalize it, we incorporate H, a set of structure-specific information in Equation 2.20.H contains a lot of information like the structures of the tree which have to be segmented, an approximative size of the structure to be segmented (the Global Prior Weight) and information about which modality is the best suited to segment this structure.

Summary
We have shown that the EM algorithm is very flexible and can be transformed to solve a lot of segmentation problems.It is very well suited for segmentation of MR brain images and we can add a lot of informations through this algorithm to enhance the segmentation.The iterative general process is divided in two steps: the expectation (E-Step) and the maximization (the M-Step).

• E-Step
Estimates a soft segmentation (p ik ), given parameter Φ (m) .The soft segmentation creates a map of probability.Each voxel contains the probabilities that it belongs to each class.It is used for the final segmentation.

• M-Step
Estimates Φ (m+1) , using the soft segmentation done in the E-Step.
-Estimates the intensity distribution for each tissue to be segmented.
-Estimates the bias field, Φ (m+1) , using the soft segmentation and the intensity class distribution The image will be segmented until the whole tree has been processed as described in section 2.4.4.
We can also describe the segmentation process in the Figure 2.3.The EM segmentation box represents the figure 2.1.To describe how it works, we use the tree structure presented in Figure 2.2.It first segments the node n = 0, i.e., BG and ICC.Once the EM Segmentation has converged, BG and ICC have been segmented and we move to the next node.Air and Skull will then be segmented.Following this process, all the structures of the tree will be segmented.

Workflow in Slicer 3
We will now present the whole segmentation workflow used in Slicer 3. It will describe all the initialisation steps done by the user, via the graphical user interface (GUI).It will also present how the whole algorithm works.We will explain why each initialisation step is important and where this information is used.

User interface
It consists of a manual initialization of the parameters which are required for the segmentation on Slicer 3. The user chooses the optimal values via the GUI.

• Step 1: Tree structure creation
We first create our problem specific tree structure for the segmentation.It will be used to define H (section 2.4.4).

• Step 2: Atlas assignment
We assign to each node of the tree, i.e. each tissue to segment, the related atlas.It will be used for the spatial information (section (2.4.1)).It implies that an atlas for each structure to segment is needed.

• Step 3: Mutlichannel segmentation
We choose the images we want to use for the segmentation.As discussed in section (2.4.2), it is usefull for the multichannel segmentation since some tissues are not visible in some sequence, in some modalities.

• Step 4: Intensity normalization
We choose the value for the intensity normalization.We normalize the intensity of the images to be segmented, regarding the related atlas.The utility of this normalization step is presented in the next section.

• Step 5: Class definition
We define mean value, variance and covariance for each class and modality.It is an accurate way to initialiaze the class tissues distribution for the algorithm.It it useful because the EM algorithm only estimates the bias field but still requires information about the tissues to be segmented.Using these values, the algorithm is initialized and estimates the bias field.

• Step 6: Hierarchical parameters
Here we set some parameters for the hierachical segmentation.We define the utility of each target image, of the atlas for each tissue to segment and approximated the size of the tissue to be segmented.This will be stored in H (section 2.4.4).For each tissue, H knows now useful informations like which volume is the more relevant for the class and the size that the class is supposed to be.
Let's take the example of the CSF.We assume a multichannel segmentation where T1 and T2 MR images are available.We give a weight of one (maximum) to the T2 target volume and zero (minimum) to the T1 target volume.It means that the only relevant information for CSF will be in the T2 volume.The algorithm will act accordingly and only use the information from T2 to perform the segmentation.It is the same for the atlas.If we set the weight of the atlas to one, the algorithm will use the spatial information.If we set it to zero, it will ignore the spatial information.

• Step 7: Registration method
We choose the type of registration we want.Different kind of registrations are available.The default registration method is a non-rigid registration.The rationale behind selecting a specific registration method is presented in the next section.

Algorithm
After all the initialization steps done via the GUI (section 2.5.1),we will now present the segmentation pipeline.

• Step 1: Intensity normalization
The intensity of the target volumes are normalized to the value that the user chose in the GUI.The normalization works in two steps: the background is first detected in the histogram, using its low intensity and the significant number of pixels of which it is created.Then it estimates the mean values of the pixels, background excluded.Then this mean value is normalized the normalization value defined by the user via the GUI.Thus, atlas and target images have the same mean intensity (background expected).This normalization is useful if we want to use the command-line interface (i.e. if we want to run a lot of segmentations without the GUI).We just choose one atlas and all the volumes to be segmented.By normalization, the mean intensity of the target images to be segmented will be the same as in the atlas.This step is usefull because now all the target images should have the same initial parameters.

• Step 2: Image registration
In order for the atlas to guide the segmentation (spatial information), it has to be aligned to the target images.The transformation between the atlas and the first related target volume is evaluated during this step, using the registration method defined by the user in the GUI.

• Step 3: Spatial prior alignement
The transformation computed during the image registration step (step 2) is applied to all the structure-specific atlases.We finally obtain new atlases which are aligned to the images to be segmented.

• Step 4: EM Algorithm in the tree structure
Now, all the required initializations and preprocessing steps have been done.The whole segmentation workflow is then applied (section 2.4.5).A typical multichannel segmentation, using two 200x300x300 targets volumes (T1 and T2 MR images) and segmenting five tissues, lasts around twenty minutes.
The whole algorithm pipeline in summarized in Figure (2.4).[18].The first step of the segmentation process consists of intensity normalization.The second consists of an estimation of the transformation to apply for the registration.Third the atlases are aligned using the transformation.Finally, the EM segmentation algorithm is performed.

Summary
In Slicer 3, the whole segmentation pipeline can be described as in Figure 2.5.There is first an initialization step, done by user via the GUI.Then, some pre-processing steps are applied in order to enhance the segmentation.Finally runs the EM pipeline to segment the MR images.

Limitations
Even if the EM segmentation pipeline is robust in Slicer 3, some limitations appear.The first problem the user has to face appears during the intensity normalization step.Because of the lack of vizual feedback, it is challenging for the user to find the good normalization value.This problem will be discussed in section 3.4.
The second problem is directly linked to the EM algorithm.As the maximization method is a local one, the class distribution, which are used for the initialization of the algorithm has to be well defined.So far, we have no possibilities to know how accurate our definition of the class is.This problem will be discussed in details in section 3.2.
Another problem appears at the same steps.The actual method for defining means and variances for each class appears not as efficient as we want it to be.This problem will be discussed in details in section 3.1.
Moreover, if the user is not familiar with the EM segment algorithm, defining good hierarchical parameters can be challenging.In section 3.5 we will provide a tool to help the user.
The last problem we encounter appears after the initialization, during the preprocessing.Only one preprocessing step (intensity normalization) is done on the target images before the atlases are registered to it.Since the EM algorithm is targeted to work on MR images, bias field is a recurrent problem.Performing the registration without bias correction might affect the results of the registration.This problem will be discussed in details in section 3.3.
In this chapter, we have described the EM segmnetation algorithm and its implementation in Slicer 3. We have concluded this the limitations of the current implementation.In the next chapter we will propose some solutions and describe them in details.

Contributions
In this chapter we will present all the contribution realized to enhance the segmentation workflow in Slicer 3. We propose solutions to enhance the class selection method and to allow the user to evaluate the selection.We will address the registration problem in a third step.Finally, we present the tools that we have added to help the user to find the optimum intensity normalization value and to estimate the hierarchical parameters.

Class Distribution Selection
During parameters initialization, the user has to define each class distribution.The previous method of selection presents some limitations and we propose a new approach.

Motivation
So far, the user has two choices to define each class distribution.The first possiblity consists of entering manually the intensities mean value and variance for each class, for each volume to be processed.This way, the user is precise to define each class.However it is challenging to find the optimum mean value and variance for each class for each volume.Morever, each time the user wants to process a new volumen, mean values and variances have to be redefined.It is not convenient and it can be time-consuming to find accurate values for the parameter initialization.
The next approach consists of defining a class model by manual sampling.For each class, the user clicks in the related part of the volume.The problem with this method is that the computation of the mean value and variance only uses a few samples.Because manual sampling is time consuming, the user will only select a reduced number of sample compared to the total number of voxels in the volume.Consequently mean values and variances are not accurate.Besides, reproducibility is an issue.Indeed it is challenging for an user to perform an identical sampling on two volumes.
Because of all these limitations, we proposed a new approach using a labelmap, to estimate each class model.

Solution
Our solution is to create a labelmap where each tissue class is given a label and a color.The relation color/class is stored in parameter H, in the EM algorithm.This relation color/class is set up during the tree creation step (section 2.5.1).
The user creates a labelmap by coloring characterisc regions for each tissue to segment in the appropriate color.This gives a spatial information to the algorithm.The labelmap allows to compute automatically the mean value and covariance of each class, for each tissue.In the next section we explain why we compute the covariance matrix instead of the variance.The labelmap sampling is a convenient solution because a large number of sample can be drawn for each class.Besides the sampling is now reproducible since we can store then re-use it later.
Figure (3.1)illustrates an example of a labelmap.Left image (A) represents an axial view of the T1 volume to be segmented and right image (B) the labelmap created for this volume.Each color represents a tissue to be segmented.

Evaluation
To estimate the importance of this contribution, we perform a simple comparison.We select manually ten points of the white matter as accurately as possible.With this manual sampling method, we estimate mean value and covariance matrix, within two volumes.The second step of the evaluation consists of estimating the same values for the same tissue using a labelmap (more than 200 sampling points).
Results: The mean values (µ M anual and µ Label ) differ slightly and covariance matrices (Σ M anual and Σ Label ) differ significantly.It means that our approach is useful and it shows the importance of having a large sample to evaluate means and covariance values.The square variance of the class for the first volume, is in position (1,1) in the covariance matrix.For the second volume, it is in position (2,2).Variance expresses the range through which the class is expected to be.The estimation of this interval is significantly more accurate with labelmap sampling than manual sampling.We can explain it with the important number of samples used for the estimation.

Class Distribution visualization
An important contribution is a tool which allows to visualize the distribution of the classes to be segmented.

Motivation
As discussed before, the algorithm is sensible to the initalization.In other words, the quality of the initialization is important.Once the parameters are chosen, the user has no tools to know if his selection is accurate.Two classes to segment can not have too close means and variances and even if the user sees the values obtained after sampling, it is challenging to know if two classes to be segmented are too similar or not.

Solution
The objective is to provide the user the most accurate and usefull vizualisation tool as possible.We first assume that each class has a normal distribution.This assumption is correct since the EM algorithm present in Slicer 3 is based on Gaussian mixture models.We decided to plot the Gaussian in the 3D space, using the multivariate normal distribution.In the 2D non-singular case, the probability density function is (see [15]).x and y are the positions of the pixel in the plane.f (x, y) will return the value (height) of the (x, y) pixel.The X axis represents a volume while the Y axis represents another volume.Let's first say that the range of the X and Y axes are the intensity range of the X and Y volumes.µ x is the mean value of the class in the X volume.µ y is the mean value of the class in the Y volume.σ x and σ y are the variance of the tissue in its respective volume.ρ is the correlation between X and Y .It indicates the strength and direction of a linear relationship between two random variables ( see [16]).
We can easily deduce ρ from the covariance matrix Σ (see [17]).Indeed, in the 2 dimensional case, the covariance matrix can be expressed as The covariance matrix and the mean values for each class to segment for each image are computed during the labelmap sampling (section 3.1).
Please note that if ρ = 1, it means that the two radom variables are the same and we can not use the same density probability function anymore.We express the classic normal distribution A is the amplitude of the Gaussian.We said that the range of the X and Y axes are the intensity range of the X and Y volumes.The problem with this approach appears if the classes to segment are not spread over all the intensities.Indeed, the vizualisation is poor then: the Gaussian is only localized in a small portion of the 2D plane.We want to "zoom" on the region of interest.We decided to change the range of the two axes.Then range will be redefined for each image.Let's present it for a given image X.The range is now defined as the difference between M ax, the maximum value extracted with the labelmap, between all the samples between all the classes for the image X, and M in, its opposite.For our particular purpose of tissues distributions vizualisation we did not use exactly the formulas (3.1) and (3.2).Indeed, we do not normalize the curves: we set 1 2πσyσy

√
1−ρ 2 to 1 in Equations 3.1 and A to 1 in Equation 3.2.We do it because we do not know the importance of each class so we do not want to "disavantage" any one in the visualization.A compromise could be to have an amplitude factor proportional to the number of pixels which constitute the class.

Evaluation
We finally obtain the results Figure 3.2 for different sampling methods.The left Figure (A) illustrates the visualization we obtain after a manual sampling and (B) the results after a labelmap sampling.The distributions slightly differ.The dark blue point, on the left corner represents the air.The skull is represented in purple, the white matter in yellow, the gray matter in green and the cerebrospinal fluid (CSF) in red.The X axis represents the T1 volume and the Y axis the T2 volume.In left image (A), the skull is significantly more spread than in right image (B), especially along the X axis (T1).According to a clinician expert, the labelmap sampling offers better results.Indeed, the skull should almost have the same variance in the T1 and T2 images, which is illustrated in Figure (B).According to the same expert, the intensities of CSF and grey matter, in the T1 volume can not be as distant as these are in image (A).The distributions of these two tissues should be close, as it is in right image (B).These two observations validates the utility of the labelmap sampling in order to have accurate tissue distributions.

MRI Bias Field Correction
The registration step could present some problems if the image to segment exhibits intensity inhomogineity.Moreover, bias field can also be a problem regarding the accuracy of the class distributions.We will first remind the problems and then propose some solutions.

Motivation
In the segmentation process, a registration step is required.Registration consists of finding a transformation to fit two images as well as possible.The main methods are described in [21].Only one pre-processing (intensity normalization) is done before the registration.The problem is that the algorithm is designed to process MR images.MR images are often corrupted by a bias field.Thus, the image to register exhibits intensity inhomogeneity.These inhomogeinities significantly affect the registration.
On Figure 3.3, we present the result of the registration between an atlas (left image (A)) and a biased MR image (right image (B)).Note that the target MR image has been normalized to have the same mean value as the atlas.The image (C) shows that the registration is poor.Moreover this bias field leads to another important issue.Indeed, if the MR image intensities are biased, the intensity of a given tissue varies significantly in the volume, even if it should not.Then the algorithm will be initialized with wrong mean and variance values and the segmentation will not be as good as it should be.

Solution
The idea simply consists of correcting the bias field of the MR image before the intensity normalization step.Thus, the registration and the class distributions will be significantly enhanced.Since the registrationand class distributions is better, it should also improve the segmentation.
To correct the bias field, we used the non-parametric approach presented by Sled et al in [19].We choose a non-parametric approach because it does not require prior information like the number of tissues to correct or the mean value of each tissue to be corrected.We implement an ITK filter [13] and [14] in Slicer 3.
We can describe the new segmentation workflow in Slicer 3 as we do in Figure 3.4.We choose not to implement it in Slicer 3 as part as the EM Segment module.Indeed, users may want to correct the bias field in MR images for other purposes.Moreover, it is possible because it would be the first pre-processing step.The user will first have to correct the intensity inhomogineity via the module developped for this particular purpose then use the corrected images in the EM segmentation module.After the bias correction, we obtain interesting results presented in Figure3.5.(A) represents the atlas to be registered.(B) represents the target volume for the registration.(C) represents the atlas after the registration.The result of the registration visually appears to be better but we can not be satisfied of this visual assessment.We need a qualitative evaluation method.

Evaluation
We evaluated accuracy of the registration using the joint histograms method.The joint histogram evaluation method is basic comparison between two images.Let A be a matrix of size W xL. W will be the intensity range of the first image used for the comparison.L will be the intensity range of the second image to be compared.The matrix is initialized to zero.Each time that in the same position, there is the same intensities in the two images, we add one in the corresponding cell in the matrix.Thus, a perfect registration, would lead to an array of zeros, expect on the diagonal.After the joint histogram creation, the value at the coordinate i, j in the matrix is the number of pixel pairs having intensity i in image onee and intensity j in image two, at the same position in the volume.
In Figure 3.6, the left image (A) and the right image (B) compares the joint-histogram of the biased image and its atlas, respectively before and after registration.The color scale used is the following one: if there are a lot of pixels in a cell of the array, the cell will be displayed in red.We the number of points decreases continously.The color scale goes from red to blue.The color fo medium intensities is yellow.
In Figure 3.7, left image (A) and right image (B) compare the joint-histograms of the image after a bias field correction and its atlas, before and after registration.
It appears that the difference between before and after registration is more significant if the bias field is corrected.Indeed, we compare the joint histograms for volumes which exhibits intensity inhomogeneity: there is no significative difference between before and after registration.That means that the registration did not improve the similarity between the images.On the contrary, if we compare joint histogram for an intensity inhomogeneity free vollume, the number of points around the matrix diagonal increases.It means the images are more similiar after than before registration.It shows the utility of the contribution.The influence of these results after the whole segmentation process will be presented in the next chapter.
Regarding the effect of the intensity inhomogeneity on tissue intensity distribution, we will present results using the tool we developped in the previous section, to evaluate the classe distribution.Using the same labelmap for sampling, we obtain two significantly different distributions for a image corrupted with a bias field and the same image corrected.In the left image (A), the tissues variances are huge, especially for the CSF (red) in the T2 volume.According to an expert, it is not a good assumption.CSF should on appear in the high intensities of the T2 volume, as presented in image (B).Image (A) means that the CSF class contents pixels of high and low intensities in the T2 volume.It is a bad assumption and we can explain it because of the intensity inhomogeneity.
Intuitively, we understand that the segmentation process should be a lot affected because of the two issues we presented (registration and distribution).Concrete results on the whole segmentation  process will be presented in the next chapter.

Registration parameters
Even if we are performing a non-parametric registration, some parameters have to be defined."Non-parametric" means no information about the volume and classes to correct.We will first present and explain the parameters.In a second step, we will propose some parameters adjusted to different problems.

• Shrink factor
It is a factor which is used to reduce the size of the image to be processed.A downsampling is done by the bias correction filter.

• Maximum number of iterations
Optimization of the bias field occurs iteratively until the number of iterations exceeds the maximum specified by this variable.

• Bias field full width at maximum iteration
The bias field is modelled with a Gaussian.This variable characterizes this Gaussian (see [20]) and can be presented as a parameter which defines the strength of the bias.
From the understanding of the parameters, it makes sense that if the time of processing has tp be reduced, the shrink factor must be increased and the maximum number of iteration must be reduced.The limitation is that it can deteriore the bias field correction.
The bias field full width at maximum iteration (BFFWMI) can also be reduced or increased, depending on the importance of the bias.If the bias field deteriorates significantly the image, an important BFFWMI must be used.On the contrary, if the intensity inhomogeneities are not significant, BFFWMI can be small.

Intensity Normalization
Another useful contribution is a tool which helps the user to determine the good normalization value.

Motivation
As discussed in section (2.5.1), at the step 4, an intensity normalization is done.We have already presented the utility of an intensity normalization in the same section.The problem is that the user has no tool to find the optimum values for the segmentation.The user has to find the mean intensity of the voxels in the MR image, excluding the background.In practise it is challenging to estimate the good value since no tool exist for this particular purpose.

Solution
We implemented a tool, to allow the user to find easily and accurately this normalization value.The first step of the work consited in creating the histogram corresponding to the image.The Y axis which presents the number of pixels for each intensity in the volume uses a log-scale because the range is large.The log scale significantly reduces the range.We then added a cursor in this histogram.With this feature, the user can choose the intensity which will be the threshold between the background and structure of interest.Finally, while the cursor is moving, our algorithm computes automatically the mean value of the voxels in the volume, from this intensity range to the highest intensity range.This is the normalization value.We present in Figure 3.9 the tool we developped.A T1 volume has been loaded.The user can move the cursor in the histogram.While moving, it returns in real-time the normalization value for the given position of the cursor in the lower frame.

Global Prior Weights Estimation
The last contribution to the EM segmentation is a tool which provides the user an easy and fast way to estimate the approximate size (also called Global Prior Weights) of each class to be segmented.

Motivation
During the segmentation process,at the 6 th step of the intialization (section 2.5.1), the user has to provide to the algorithm an estimation of the size of each tissue to be segmented.First of all if there are a lot of structures to segment, this step is time consuming.Moreover, the user may not know at all which size to choose.A tool to estimate of the good sizes to choose is needed.

Solution
Our solution is to divide the problem in two parts.We provide the user with a real-time feedback regarding the global prior weights estimation.The second part consists of developping an algorithm which fills automatically the tree.
The new tool is presented in Figure 3.10, in the left image (A).When the user moves a cursor, it changes the intensity range for each class.The tool provides a real time feedback to the user.The user sees Figure 3.10) (B) which is updated in real time, regarding the position of all the cursors consequence.Clicking on "update tree" (in Figure 3.10) (A)), the estimated size of a class is computed and the information is stored into the tree structure (parameter H).In this chapter we have presented our contribution to the EM segmentation algorithm.The illustrations shows that the selection of the parameter has been improved: in the next chapter we present the impact on the segmentation.

Results and discussion
We will now discuss the impact of the enhanced usability on the segmentation results.Firstly we present the influence of each contribution on the final segmentation.Secondly suggest some potential improvements for the EM segment module in Slicer.

Results
Here we start with a presentation of different results, using the different contributions.The results obtained will then be reviewed by a expert radiologist.The expert who accepted to evaluate the results is Dr. KIKINIS, professor of radiology at Harvard Medical School.We report his assessments for each segmentation in the "expert's validation" section.We work with the same dataset for all the segmentation.We chose the datasets with significant intensity inhomogeneity in order to show the importance of our contributions.

Original segmentation
First we present the segmentation obtained without any contribution.Based on this segmentation, the importance of the contributions will be presented in the next sections.First we describe the testing process.Second we present the results of the segmentation.Then the clinical expert's assessment will be reported.Finally, we discuss about the expert's observations.

Method
Thanks to Sonia Pujol, a tutorial for the EM segmentation in Slicer 3 is available at http://www.na-mic.org/Wiki/images/2/2f/AutomaticSegmentation_SoniaPujol_Munich2008.ppt.We follow it with an image which exhibits intensity inhomogeneity.The atlases we use for the segmentation are the one available at http://www.na-mic.org/Wiki/images/b/b7/AutomaticSegmentation.tar.gz.

Expert's validation
Based on the assessment of the expert, the grey matter is significantly over estimated in the whole volume.The segmentation is poor and can not be use for any medical application purpose.

Discussion
We speculate that the reasons of the poor segmentation are the incorrect registration and the class ditributions which is poorly defined.It leads to inaccurate means and excessively large variances which deteriorate a lot the segmentation.This is also illustrated Figure 3.8 where the right image shows that grey matter (light blue) is overestimated in the right image (B) in comparison to the left image (A).

Bias corrected segmentation
Here we present the results of a segmentation using the intensity inhomogeinity correction tool developped.First we describe the testing process.Second we present the results of the segmentation.Then the clinical expert's assessment will be reported.Finally, we discuss about the expert's observations.

Method
Similarly to the previous section, we follow the EM segmentation tutorial with one extra-step: the bias field correction.

Results
The results of the final segmentation are presented in Figure 4.2.The left image (A) presents an axial view of the T1 target volume to be segmented.The right image (B) presents and axial view of the segmentation's results.

Expert's validation
The expert pointed misegmentation in the skull area.Indeed, some voxels in the skull are incorrectly classified as grey matter and white matter.The expert explained that this misclassification is due to partial volume effect.The partial volume effect is caused by imaging voxels containing two different tissues (skull and air in our case).The mesured intensity for that voxel will be an average of the two tissues at that location.Another misclassification occurs in areas where the skull is classified as air.The reason is that in some regions, the skull is porous and contains air.Even is pores are not significant, it leads to partial volume effect the misclassification.

Discussion
Based on our understanding of the algorithm this better segmentation can be explained by a better registration as shown section (3.3).Moroever, as demonstrated it in Figure 3.8, the tissues have a better distribution, especially the grey matter (light blue) in this case.The misclassifications in the skull could be corrected by a pre-processing step using skull stripping algorithm.

Label map sampling segmentation
This segmentation will now let us evaluate the influence of the sampling method on the segmentation.First we describe the testing process.Second we present the results of the segmentation.Then the clinical expert's assessment will be reported.Finally, we discuss about the expert's observations.

Method
Similarly to the previous section, we follow the EM segmentation tutorial .The extra step that we take are the bias field correction and the initialization of the class distribution using a labelmap sampling.

Results
The results of the final segmentation are presented in

Expert's validation
The skull is better segmented than in the manual sampling segmentation.Moreover, there is less misclassification than before.Nevertheless, in cerebulum, the segmentation differs from the previous one.It is underestimated with the labelmap sampling whereas before it was over estimated.The cerebulum is a specific area of the brain because there is a lot of white and grey matter in a small area.The segmentation is challenging and similarly as the point raised during the validation of the bias field correction segmentation, the results are due to partial volume effects.According to our expert segmentation of this area is a complex task.A solution to improve the segmention the cerebrulum properly, would be to extract the region of interest and to apply a finer algorithm to this area.

Discussion
We explain the enhanced classification because the tissues have more located distributions than in the previous segmentation, as we se in Figure 3.2.More screenshots with a better resolution are presented in App.B to have a better understanding of the interest of of a manual sampling instead of a labelmap sampling.

Perspectives
Our work focused on improving the usability of EM segment and the selection of the optimum parameters.As a future direction, some effort could spend on the optimization of the algorithm.Indeed, the long processing time is a major drawback of the current implementation.
On one side the processing time both for the bias field correction and for the EM segmentation needs to be improved.A number of users reported this issue as the main limitation of the current implementation.We observe that the processing for 255x255x100 volume is approximatly 30 minutes.Towards this goal the code of the bias field correction algorithm needs to be optimized.
Another issue occurs with images where the distributions of two classes have similar means.In this case, the exploration of the parameters is challenging, especially for the global prior weights.
A useful contribution could be to represent each class as a gaussian instead of a finite portion of the histogram.

Conclusion
We have presented an improved algorithm and implementation of the EM segmentation.Our main contributions are the addition of a bias field correction step and the creation of tools to select the optimum parameters.First we implemented a tool to apply the non parametric method for automatic correction of intensity uniformity in MR images of Sled et al. [19].Our second contribution enables the visualization of the tissue classes as Gaussian distributions and the selection of the optimum statistical parameters.
We have reported our segmentation results and our qualitative validation performed by a clinical expert.The expert assessed our interface as a significant improvement in usability and validated our segmentation improved results.
This work opens new perspectives for improving the impact of medical imaging in a clinical environment.We suggest that the future effort should focus on reducing the processing time.The interaction with clinicians at Harvard Medical School make us believe that the outcome of this internship will be developped further and ultimatly make a difference for the patients.

Figure 2 . 1 :
Figure 2.1: Illustration of the basic EM algorithm

Figure 2 . 2 :
Figure 2.2: A simple tree structure of the brain The brain is first divided into the background (BG) and the intracranial cavity (ICC).Second the BG is dived into Air and Skull.Finally, the ICC is divided into white matter (WM), grey matter (GM) and cerebrospinal fluid (CSF)

Figure 2 . 3 :
Figure 2.3: Illustration of the EM segment algorithm in Slicer

Figure 2 . 4 :
Figure 2.4: Algorithm segmentation pipeline in Slicer 3 Courtesy Pohl et al.[18].The first step of the segmentation process consists of intensity normalization.The second consists of an estimation of the transformation to apply for the registration.Third the atlases are aligned using the transformation.Finally, the EM segmentation algorithm is performed.

Figure 3 . 1 :
Figure 3.1: Axial view of the labelmap Left image (A) represents an axial view of the T1 volume to be segmented.Right image (B) represents the volume overlays by the labelmap.Pink represents the air.White represents the skull.Yellow represents the white matter (WM).Blue represents the grey matter (GM).Red represents the cerebrospnial fluid(CSF).

Figure 3 . 2 :
Figure 3.2: Distribution of the class to be segmented in an homogeneous volume Left image (A) represents the distribution of the class after a manual sampling.Right (B) image represents the class distribution after a labelmap sampling.Air is plotted in dark blue, skull in blue, WM in yellow, GM in green and CSF in red.X axis represents a T1 volume and Y axis represents a T2 volume.

Figure 3 . 3 :
Figure 3.3: Registration of a MR image which exhibits intensity inhomogeneity Image (A) is the atlas to be registered.Image (B) is the target of the registration.Image (C) is the image (A) after registration to image (B)

Figure 3 . 4 :
Figure 3.4: New algorithm segmentation pipeline in Slicer 3 First the intensity inhomogeneity are corrected.Second the intensity is normalized.Third an estimation of the transformation to apply for the registration is done.Fouth the atlases are aligned using the transformation.Finally, the EM segmentation algorithm is performed.

Figure 3 . 5 :
Figure 3.5: Registration of a MR image without intensity inhomogeneity Image (A) is the atlas to be registered.Image (B) is the target of the registration.Image (C) is the image (A) after registration to image (B)

Figure 3 . 6 :
Figure 3.6: Joint histograms to evaluate the registration of the inhomogeneous MR image Figure (A) represents the joint histogram of the inhomogeneous MR image and its atlas before registration.Figure (B) represents the joint histogram of the inhomogeneous MR image and its atlas after registration.

Figure 3 .
8 shows the two different distributions.The relation tissue/color is the same as the one in Figure 3.2.The left image (A) presents the distribution before bias field correction.The right (B) presents the distribution in the corrected volume.

Figure 3 . 7 :
Figure 3.7: Joint histograms to evaluate the registration of the homogeneous MR image Figure (A) represents the joint histogram of the homogeneous MR image and its atlas before registration.Figure (B) represents the joint histogram of the homogeneous MR image and its atlas after registration.

Figure 3 . 8 :
Figure 3.8: Distribution of the class to be segmented in an inhomogeneous volume Left image (A) represents the class distribution in an homogeneous volume.Right image (B)represents the distribution of the class in a volume which exhibits intensity inhomogeneity.Air is plotted in dark blue, skull in blue, WM in yellow, GM in green and CSF in red.X axis represents a T1 volume and Y axis represents a T2 volume.

Figure 3 . 9 :
Figure 3.9: Tool developped for the intensity normalization parameter estimation The user moves the cursor and get an immediate feedback about the normalization value to use in the "recommended normal" box

Figure 3 . 10 :
Figure 3.10: Tool developped for the global prior weights estimation Left image (A) presents the tool we developped for the global prior weights estimation.Right image (B) shows the visual feedback provided to the user given the tool presented in (A)

Figure 4 .
Figure 4.1 shows the result of the segmentation.The left image (A) presents an axial view of the T1 target volume to be segmented.The right image (B) presents and axial view of the segmentation's results.

Figure 4 . 1 :
Figure 4.1: Results of the segmentation without bias field correction Left figure (A) is an axial view of the T1 the target volume to be segmented.Right figure (B) is an axial view of the segmentation's result.Air is represented in pink, skull in white, white matter in yellow, grey matter in blue and cerebrospinal fluid in red.

Figure 4 . 2 :
Figure 4.2: Results of the segmentation with bias correction Left figure (A) is an axial view of the T1 the target volume to be segmented.Right figure (B) is an axial view of the segmentation's result.Air is represented in pink, skull in white, white matter in yellow, grey matter in blue and cerebrospinal fluid in red.
Figure B.3.The left image (A) presents an axial view of the T1 target volume to be segmented.The right image (B) presents and axial view of the segmentation's results.

Figure 4 . 3 :
Figure 4.3: Results of the segmentation with labelmap sampling Left figure (A) is an axial view of the T1 the target volume to be segmented.Right figure (B) is an axial view of the segmentation's result.Air is represented in pink, skull in white, white matter in yellow, grey matter in blue and cerebrospinal fluid in red.

Figure B. 2 :
Figure B.2: Coronal view of the segmentation with different sampling methods Left figure (A) is an axial view of the T1 the target volume to be segmented.Right figure (B) is an axial view of the segmentation's result.Air is represented in pink, skull in white, white matter in yellow, grey matter in blue and cerebrospinal fluid in red.

Figure B. 3 :
Figure B.3: Sagittal view of the segmentation with different sampling methods