Gaussian Intensity Model with Neighborhood Cues for Fluid-Tissue Categorization of Multi-Sequence MR Brain Images

. This work presents an automatic brain MRI segmentation method which can classify brain voxels into one of three main tissue types: gray matter (GM), white matter (WM) and Cerebro-spinal (cid:13)uid (CSF). Intensity-model based classi(cid:12)cation of MR images has proven problematic. The statistical approach does not carry any spatial, textural and neighborhood information in it. We propose to use a computationally fast and novel feature-set to facilitate voxel wise classi(cid:12)cation based on regional intensity, texture, spatial location of voxels in addition to posterior probability estimates. Information available through T1-weighted (T1), T1-weighted inversion recovery (IR) and T2-weighted FLAIR (FLAIR) MRI sequences was also leveraged. An aggregate overlap of 90.21% for all intracranial structures was reported between the automatic classi(cid:12)cation and available expert annotation as measured by the DICE coe(cid:14)cient.


Introduction
Magnetic Resonance Imaging (MRI) has become important in the diagnosis and monitoring of brain anatomical structures and has emerged as a key supportive therapeutic outcome measure in clinical diagnosis and biomedical research [15].In clinical and medical studies on brain anatomical structures, a successful partition of images into Gray Matter (GM), White Matter (WM) and Cerebrospinal Fluid (CSF) is often an important first step.This work presents an automatic brain MRI segmentation method which can classify brain voxels into one of three main tissue types.
Volumetric analysis of different parts of the brain are useful in assessing the progress or remission of various other diseases, such as Alzheimers disease, epilepsy, multiple sclerosis, and schizophrenia.In neurodegeneration studies of schizophrenic patients, the volume changes of gray matter in thalamus, frontal and temporal lobes as well as CSF in ventricles are of relevance [9] [13].It is also of considerable interest to study regional volumes of GM and WM across different developmental stages of the human brain [5].Lesion segmentation is also another important application in clinical studies [16].
Amongst majority of biomedical research community, the image segmentation is still supervised by professional experts on a slice-by-slice basis.This process is not only labor intensive but also introduces large intra-and interobserver variability on account of partial volume effect, variability of scanning procedures and human factors.Further, with advancement of technology, multiple MRI modalities are available viz.T1-weighted (T1), T1-weighted inversion recovery (IR) and T2-weighted FLAIR (FLAIR).For a human expert it may not be practicable to observe multiple modalities simultaneous and draw their inferences.Therefore, in order to reduce the subjectivity in analysis and to make use of information available through multiple modalities it is highly desirable to apply automatic image analysis methods for biomedical studies involving large image data.
This work presents an automatic brain MRI segmentation method which can classify brain voxels into one of three main tissue types: Gray matter (GM), White matter (WM), and Cerebro-spinal fluid (CSF).Gaussian distributions are used to model voxel intensities from multi-sequence MRI scan (T1, IR, FLAIR).Since statistical models do not take into account spatial information, we propose to use neighborhood intensity cues, texture cues and positional information along with cues obtained from statistical modeling to classify the voxels using support vector machines.This work was done as part of a contest titled MR Brain Image Segmentation 20131 to be held during 16th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI-2013).
The contest dataset consists of 20 fully annotated multi-sequence (T1, IR and FLAIR) 3T MRI brain scans.These scans have been acquired at the UMC Utrecht (The Netherlands).All scans are bias corrected, and the thick-slice scans are registered.Each thick sliced sequence contains 48 image slices, each of size 240X240 voxels.Out of the 20 scans, 5 scans were provided with annotations to be used as training set.12 scans were provided without annotations to be used as test set.Remaining 3 scans will be used for onsite evaluation of the algorithms.Further details of the dataset can be found on the contest website 1 .
The organization of the rest of the paper is as follows.A concise literature review is presented in Section 2. Following that, Section 3 describes in detail the procedure that was used to solve the problem at hand.Results were presented in Section 4. Section 5 concludes the paper and gives direction for future work.

Literature Review
In recent years, researchers have proposed various approaches for brain MRI classification.A wide range of statistical methods are proposed.These approaches assume that mixed voxel intensities reflect distinct tissue groups, and individual voxels are assigned to different groups through modeling the intensity histogram as a mixture of probability distributions.Ashburner et al. [1] have proposed a modified version of the mixture model algorithm for analyzing spatial distributions of tissue groups.Rajapakse et al. [17] have proposed a Bayesian maximum a posteriori (MAP) model for image segmentation by incorporating Markov random fields (MRFs) as spatial priors for tissue distributions.Expectation Maximization (EM) is yet another popular approach used by various medical imaging researchers in varied domains.Zhang et al. [19] improved EM by using hidden Markov Random Field (HMRF) model in which state sequence is estimated through observations.The advantage of the HMRF model is that the spatial information is encoded through the influences of neighboring sites.
Yet another class of approaches includes clustering.Fuzzy-based approaches generalize the K-means algorithm to allow for soft segmentations such that each voxel can be assigned to more than one class of tissue [6] [12].On the other hand, the adaptive mean shift methods enable the integration of intensity and spatial features [14][7].Glass et al. have used the Kohonen Map and multi-layered backpropagation neural networks for segmentation in inversion recovery (IR) images [4].
In addition to these, some miscellaneous approaches have also been proposed.Li et al. [10] have proposed an interactively guided method utilizing dual-front active contours which minimizes image-based energies.In another of their work Li et al. have combined the watershed transform and region-based level set method [11].Tian and Fan [18], have proposed a feature set which uses pixel intensity of the voxel and of its spatial neighbors.Classification of voxels is done with self-organizing map (SOM) neural network.
In such a well-established research area, there is a tremendous need for fair comparison of these methods with respect to accuracy and robustness.Many of the previous algorithms make use of publicly available databases like the Alzheimers Disease Neuroimaging Initiative (ADNI) and the Internet Brain Segmentation Repository (IBSR).However, these datasets either lack full manual segmentations (ADNI) or provide low-field (1.5T) single sequence (T1-weighted) MRI data (IBSR).Most of the algorithms have been tested on privately held datasets.The aim of this challenge is to compare algorithms for segmentation of gray matter, white matter and cerebrospinal fluid on a common dataset consisting of 3T MRI multi-sequence scans of the brain.

Methodology
The main stages of the proposed segmentation scheme are : 1) Preprocessing, 2) Statistical Modeling, 3) Feature Extraction, 4) Supervised Classification and 5) Post Processing. Figure 1 shows the overview diagram of the proposed segmentation scheme.
where, I is the original image, M ax(I) and M in(I) denote the maximum and minimum intensities of I respectively and I N is the intensity normalized image.

Region of Interest (ROI) selection
As explained in Section 1, for each patient there are three thick sliced sequences viz.T1, IR and FLAIR.Each thick sliced sequence contains 48 image slices, each of size 240X240 voxels.Since, the proposed segmentation scheme is voxel based, it presents a huge amount of data to be processed.Hence there is a need to select a Region of Interest (ROI).An edge map of each slice in a given volume in T1 sequence was computed using Canny Edge Detection [3] as shown in Fig. 2a.N-Radially outward directed lines were constructed from the center of the slice as shown in Fig. 2b.Each line was drawn using Eq. 3.
Where N = 60, k ∈ {0 to N − 1}, r = 1 to 115 and C x , C y denote the center of image slice.Intersection points of these lines with the outermost edge of the edge map were computed as shown in Fig. 2c.A closed cubic spline was fitted through these computed intersection points as shown in Fig. 2d.Area enclosed in this closed cubic spline was used as ROI.T1, IR and FLAIR scans are aligned to eliminate the influence of different registrations.Hence, the ROI obtained from volumes of T1 sequence can be used for IR and FLAIR sequences as well.

Statistical Modeling
MRI provides good contrast between the different soft tissues and fluids in the brain.Scans obtained with different MRI modalities (T1, IR, FLAIR) help in analysis of different parts of the brain.For example, using T1 scan, the CSF appears darker than WM and GM.CSF is bright, but fat (and thus white matter) is darker in IR.The FLAIR modality suppresses CSF and white matter is highlighted.All these observations can be verified from histogram plots for each of the labels (CSF, GM, WM) for each of the modalities (see Fig. 3).
Let, l ∈ {CSF, GM, WM} at a given voxel s ∈ {T1, IR, FLAIR} sequence X = Intensity tuple at a given voxel representing intensities of three sequences viz.T1, IR, FLAIR X s = Intensity at a given voxel for a sequence s.
Intensity distribution of each label (CSF, GM, WM) in each of the sequences can be approximated as a Gaussian distribution with mean µ and variance σ 2 .Fig. 3 shows the histogram of all the three region viz.CSF, GM, WM of an image volume in T1, IR and FLAIR sequences.Let Θ s,l = {µ s,l , σ 2 s,l } denote the mean Gaussian parameters for all the volumes in the given training dataset.Hence underlying likelihood of each sequence were modeled as three independent Gaussian distributions.
A stationary prior probability distribution on label class was used.This was assigned as the fraction of occurrence of label l in the ground truth and is denoted by p(l).Probability density of voxel intensity conditioned on a given label for a sequence is defined as We assume statistical independence of the voxel intensities in each of the sequences (T1, IR, FLAIR).We may then write the combined probability density for all the sequences as, Posterior probability for the labels, given observed intensity data was obtained using Bayes' rule as follows, where P(X) is an unimportant normalizing constant and hence was ignored.

Feature Extraction
Empirical results indicate that label based (CSF, GM, WM) intensity distributions in image volume sequences(T1, IR, FLAIR) exhibit a partial Gaussian shape but were modeled exactly as Gaussian distributions in the Section 3.2.The statistical approach also does not carry any spatial, textural and neighborhood information in it.Hence it presents a need for a feature set through which an accurate and robust segmentation can be achieved.A novel feature set was developed to facilitate voxel wise classification based on regional intensity, texture, spatial location of voxels in addition to posterior probability estimates as developed in Section 3.2.Following features were computed for each of the voxel in ROI : To evaluate the mean of far neighborhood efficiently, mean filtering was employed in which each voxel value in an image slice is replaced with the mean value of its neighbors.This was achieved using convolution of the image slice with a 7X7 square kernel (k 7 ).
where J 7 is 7x7 ones matrix.
where * stands for convolution and I N is the intensity normalized image.
An edge map was computed using canny edge detection as discussed in Section 3.1.Number of edges in the near neighborhood were computed by convolution of edge map with a 3X3 square kernel J 3 , where J 3 is 3x3 ones matrix.Hence, each voxel value was replaced with the number of edges in its near neighborhood.
where * stands for convolution and I ed is the edge map image.The Power Spectral Density in the Far neighborhood is essentially mean of squared neighborhood values.It was computed by taking the point-wise product of the normalized image slice with itself followed by mean filtering which was achieved by convolution of the obtained point-wise product with a 7X7 square kernel (k 7 ) as depicted in following equation.
It was observed that initial slices of MR Image data consist of higher number of anatomical structures of brain while they attain regularity over the higher slice numbers and are constituted mainly by CSF, GM and WM.For the available thick slice dataset containing 48 slice MRI data three logical partitions based on the observed anatomy of the brain were made.Partition 1 contained data from slice number 1 to 14, partition 2 contained data from slice number 14 to 24 and partition 3 contained data from slice number 25 to 48.Above described feature set was computed separately for each of the partition.Thus in all 24 features were extracted for each voxel in ROI.

Classification
The feature set described in Section 3.3 was used for classification of voxels into one of the four class viz.CSF, GM, WM and Background.Supervised classification was performed using Multi Category Support Vector Machine.Since features were extracted for each voxel, the number of samples were prohibitively large to enable SVM to learn a model.Random sampling of the data provided a simple solution to this problem.Separate classifiers were learned for each of the three logical partitions described in Section 3.3 using the computed training feature data for each partition.

Post Processing
The resultant voxel wise classification in label masks contains small noisy blobs.Connected components analysis (CCA) labels the blobs in a label image, as per its connectivity.The labels thus formed were used to iterate through each of the blob to extract the blob area.Isolated blobs with very small areas were removed.

Results and Discussion
The Gaussian intensity models treat each voxel independently and thus suffer from an intrinsic limitation that such models do not take into account spatial information.This section visually demonstrates the limitation of using only an intensity based model.Further it shows qualitatively the advantage of including spatial features as proposed in Section 3.3 in addition to the features derived from intensity model.Finally we also present quantitative evaluation of our segmentation using metrics as proposed by Babalola et al. [2].
It can be observed from the ground truth (Fig. 4b) that in partition 1 there are regions in brain scan which cannot be classified only on the basis of statistical approach.Such regions also require edge information, textural cues, and spatial information for accurate classification.It is evident from the Fig. 4b and Fig. 4d that ground truth and segmentation result from our algorithm are comparable.Due to similar intensity composition of the skull and white matter, parts of the skull were accidentally classified as white matter when using only the statistical model.However, such misclassifications are eliminated by use of proposed feature set.This can be observed in Fig. 4. The classification scheme described in Section 3 was applied on the contest testing set which consisted of 12 MRI brain scans.The automatic segmentation results were compared to manual annotated scans.For each tissue type (gray matter, white matter and cerebrospinal fluid), the Dice coefficient (DC) , the Modified Hausdorff distance (MHD) measure and absolute brain (gray matter + white matter) volume difference (ABVD) were calculated.Dice coefficient is a measure of volume overlap between the manually annotated results and automatic results.
Hausdorff distance is a boundary based measure.It measures the degree of separation of the manually marked contour and automatically derived contour.Modified Hausdorff Distance (MHD) is defined as the 95th-percentile Hausdorff distance(HD).
Where C is the total number of pixels in the image and  1 shows the results summarized by the organizers using the above mentioned metrics.It can be noted that there was an overlap of 74.22% between the manually marked and automatic results for gray matter as measured with Dice coefficient.Similarly, an overlap of 78.72% and 70.23% was reported for white matter and cerebrospinal fluid respectively.An aggregate overlap of 90.21% of all intracranial structures was reported.Since the previous methods have used different data-sets and often privately held data-sets, it was not possible to compare our results with those results.The results from other competitors of this contest were not available as of writing of this paper.The comparison results of this algorithm with other competitors will be available on the contest site after the onsite contest.
The prototype for the proposed method was built using MATLAB on Windows 7 on an Intel i5 machine with 4GB RAM.Windows version of SVM Light multi-category [8] implementation was used.Approximately 5 minutes are required on the above described machine to process each volume.

Conclusion and Future Work
In this paper we proposed to use information from multiple MRI sequences to automatically classify brain voxels into one of three main tissue types: gray matter (GM), white matter (WM), and Cerebro-spinal fluid (CSF).Since the intensity model assumes voxel independence, it fails to capture the spatial information.To incorporate the spatial information, we proposed to make use of regional intensity, texture, and location cues of voxels in addition to posterior probability estimates.As of now, the system does not require user interaction and an aggregate overlap of 90.21% for all intracranial structures was achieved between the automatic classification and available expert annotation as measured by the Dice coefficient.
Use of 3-D neighborhood may prove effective in further boosting the classification rates.In statistical models such as the hidden Markov model spatial information is encoded in the model.Use of such a model can be a subject of future research.Further, with advent of faster GPUs it may be possible to reduce computational time for our system.Parallel deployment of the algorithm may be required for processing massive amount of data that is needed to be processed in a typical clinical setting.A multi-site, multi-discipline collaborative study may be a vehicle to develop practical computerized analysis tools to study the brain and to understand and combat various brain ailments.

Table 1 :
|x| represents cardinality of any set x. S and G represent segmentation result and ground truth.A(S) and A(G) are the areas of the close boundary of segmentation results and manual delineation respectively.For boundary based measures S and G are closed boundaries of segmentation results and manual delineations.V is defined as the number of absolute brain (gray matter + white matter) labeled voxels multiplied by the voxel dimensions.Overall Dice Coefficient(%), Modified Hausdorff Distance (mm) and Absolute Brain Volume Difference (%) for the testing set.