Modiﬁed Expectation Maximization Method for Automatic Segmentation of MR Brain Images

. An automated method of MR brain image segmentation is presented. A block based Expectation Maximization Algorithm is proposed for the tissue classiﬁcation of MR brain images. The standard Gaussian Mixture Model is the most widely used method for MR Brain image segmentation and Expectation Maximization algorithm is used to estimate the model parameters. The Gaussian Mixture Model considers each pixel as independent and does not take into account the spatial correlation between the neighboring pixels. Hence the segmentation result obtained using standard GMM is highly sensitive to Intensity Non-Uniformity and Noise. The image is divided into blocks before applying EM since the GMM is preserved in the local image blocks. Also, Nonsub-sampled Contourlet Transform is employed to incorportate the spatial correlation among the neighbouring pixels. The method is applied to the 12 MR brain volumes of MRBRAINS13 test data and the White Matter, Gray Matter and CSF structures were segmented.


Introduction
Segmentation of MR Brain images into White Matter, Gray Matter and Cerebrospinal Fluid(CSF) is an important step in quantitative morphology of the brain.The segmentation of MR Brain images aids in detection of various brain diseases like Multiple Sclerosis, Alzheimer's disease, epilepsy and others.Manual Segmentation is time consuming and also it leads to intra and inter observer variability.Hence automatic segmentation methods are in research focus and wide range of methods are proposed in the literature.The standard GMM and EM for estimating the model parameters is widespread method for image segmentation.The drawbacks of the method are that it is sensitive to noise and it does not take into account the spatial correlation between the pixels.Koen Van Leemput et al. [3][4][5] developed the model based method for automated classification of MR brain images.The method applies an iterative Expectation Maximization algorithm that interleaves pixel classification with estimation of class distribution and bias field parameters, improving the likelihood of the parameters at each iteration.Nathalie Richard et al. proposed distributed Markovian segmentation performed under a collabarative and decentralized strategy to ensure the consistency of segmentation over neighbouring zones.[6].Thanh Minh Nguyen and Q. M. Jonathan Wu proposed an algorithm which incorporates the spatial relationship between neighboring pixels by replacing each pixel value in an image with the average value of its neighbors including itself.[1].Zexuan Ji, Yong Xia et al. proposed the fuzzy local gaussian mixture model in which a truncated gaussian kernel function is used to impose the spatial constraint and fuzzy memberships are employed to balance the contribution of each GMM.[2].To overcome the effects of bias field and noise in segmentation accuracy, a wide range of techniques are being referenced.Since the GMM is preserved in the local image blocks, block based EM segmentation is suggested to remove the effect of bias field in segmentation.Also, to incorporate the spatial correlation among the neighboring pixels, Nonsubsampled contourlet transform low pass filter is employed.

Gaussian Mixture Model and EM Segmentation
The EM algorithm is a general technique for finding Maximum Likelihood parameter estimates in problems with hidden data.[3][4][5].It first estimates the hidden data based on the current parameter estimates.The initial parameters mean and variance of the classes are calculated from the histogram of the image.The estimated completed data consisting of both observed and hidden data are then used to estimate the parameters through maximizing the likelihood of the complete data.It consists of two steps E step and M step.The E step computes expectation of the unobserved data.The M step computes Maximum Likelihood estimates of the unknown parameters.The process is repeated until it converges.The likelihood increases with each iteration.For segmentaiton problem, the observed data are the signal intensities, the missing data are the classification of images, and the parameters are class-conditional intensity distribution parameters.Each voxel value is selected at random from one of K classes.The EM algorithm then interleaves class-conditional parameter estimation and it performs statistical classification of the image voxel into K classes.Each class is modeled by a normal distribution G σ (y − µ) with mean µ and variance σ 2 .The probability density that class j has generated the voxel value y i at position i is with Γ i {j/j = 1...K} the tissue class at position i and θ j = {µ j , σ j }, the distribution parameters for class j.Defining θ = {θ j } as the model parameters, the overall probability density for y i is which is a mixture of normal distributions.Since all the voxel intensities are assumed to be statistically independent, the probability density for the image y given the model is The maximum likelihood estimates for the parameters µ j and σ j can be found by maximization of p(y/θ), equivalent to minimization of − i log e (p(y i /θ)).
The expression for µ j is given by the condition that Differentiating and substituting p(y i /Γ i = j), θ j by the Gaussian distribution 1 yields where Bayes' rule was used.
Rearranging the terms yields the expression for µ j .
Equation 6 performs classification while 7 and 8 are parameter estimates.The algorithm fills in the missing data during step 1 and then finds the parameters that maximize the likelihood for the complete data during step 2. The likelihood is guaranteed to increase at each iteration.The method is simple when compared to other segmentation methods like fuzzy classification, deformable model etc.,.
The drawbacks observed in the method are that it treats each pixel as independent and it does not take into account the spatial relationship between pixels.This leads to misclassification of the pixels to different classes and hence the observed accuracy is low.Also if the noise level and intensity in-homogeneity increase, the accuracy decreases rapidly.

Nonsubsampled Contourlet Transform
The NSCT is constructed by combining non sub sampled pyramids and non sub sampled directional filter banks as shown in fig (1).The Nonsubsampled pyramid structure results the multi scale property.The Nonsubsampled directional filter bank results the directional property.The NSCT has similar subband decomposition as that of contourlets but without up samplers and down samplers.The standard GMM considers each pixel as independent and does not take

Method
The algorithm is summarized in the following steps.1. Skull Stripping is done to extract the brain region using threshold and morphological operations 2. The Nonsubsampled Contourlet Transform is applied to decompose the image into low pass sub band and band pass directional subbands.
3. The low pass sub band image is divided into 16 blocks each of 60 x 60 pixels.

Then Expectaion Maximization
Segmentation is applied to each block.4. The number of classes for EM segmentation is selected as 5.
5. The tissues are labeled as Gray Matter, White Matter and Cerebro Spinal Fluid(CSF).

Results and Discussion
In this section, the performance of the algorithm on 12 MR brain volumes of MRBRAINS13 test data is presented.The system configuration used is Intel Core 2 Duo CPU @2.53GHz with 1.98GB of RAM.The algorithm is executed using Matlab.Only T1 weighted images are used and no training data is used.

Performance Metrics
Dice Metrics The Dice Similarity (DS) is a measure of the similarity between manual segmentation (X) and the automatic segmentation (Y).The equation for calculation can be written as Absolute Volume Difference The total volume of the segmentation is divided by the total volume of the reference.From this number 1 is subtracted, the absolute value is taken and the result is multiplied by 100.This value is 0 for a perfect segmentation and larger than zero otherwise.

Discussion
In order to overcome the drawbacks of the standard EM for MR brain image segmentation, a modified EM is proposed and implemented.The experimental results prove that the method achieves overall Dice Similarity Index of 0.79, 0.64 for GM and 0.77 for WM.The results achieved are comparable to that of the state-of-the art algorithms.The computational complexity is very much reduced and only T1 weighted images are used.The method is a fully unsupervised technique and as such no training data are required.The time required for segmentation of a single volume is around 300 to 350 seconds which include the time for brain extraction also.The method may be extended for pathology identification and segmentation like Multiple Sclerosis lesion segmentation.The method also can be used to segment MR brain images with more bias field and noise.