Multimodal MR Brain Segmentation Using Bayesian-based Adaptive Mean-Shift (BAMS)

. In this paper, we validate our proposed segmentation algo-rithm called Bayesian-based adaptive mean-shift (BAMS) on real multimodal MR images provided by the MRBrainS challenge. BAMS is a fully automatic unsupervised segmentation algorithm. It is based on the adaptive mean shift wherein the adaptive bandwidth of the kernel for each feature point is estimated using our proposed Bayesian approach [1]. BAMS is designed to segment the brain into three tissues; white matter (WM), gray matter (GM) and cerebrospinal ﬂuid (CSF). The performance of the algorithm is evaluated relative to the manual segmentation (ground truth). The results of our proposed algorithm show the average Dice index 0 . 8377 ± 0 . 036 for the WM, 0 . 7637 ± 0 . 038 for the GM and 0 . 6835 ± 0 . 023 for the CSF.


Introduction
Automated segmentation of MR images of the brain is an active area of research in the field of neuro imaging.The resulting segmentation yields a patient-specific labeling of individual tissues which makes possible for the quantitative characterization of these tissues (for example, in the study of Alzheimer's disease and multiple sclerosis).The majority of unsupervised methods that have been proposed for automated segmentation of brain tissues are based on statistical parametric models.These methods assume some distributional form for the underlying probability distribution of the data and seek to estimate its parameters.Some of these [2,3] are purely voxel wise intensity based clustering methods.A downside of these is that they may give poor tissue classifications in the presence of additive noise and the multiplicative bias field [4].Some of the parametric methods [5,6] using a Markov random filed (MRF) statistical spatial model to improve the smoothness of segmentation.A drawback with these approaches is that MRF algorithm is computationally expensive and requires critical parameters settings at higher dimensional feature space [4].
An alternative unsupervised approach that doesn't require many parameters, incorporates the spatial information easily into a higher dimensional feature space (multimodal MR images) is mean shift (MS) clustering.In MS, only the parameter that influences the clustering is called bandwidth of the kernel.A couple of MS methods [4,7] based on the adaptive bandwidth have been proposed for brain tissue segmentation in MR images.The adaptive bandwidth estimator [8] used in [4], is based on the k nearest neighbour (kNN) distance.A downside is that this approach is known to be biased by outliers for Euclidean distance [9].Drawback of adaptive bandwidth estimator used in [7] include that it requires an initial density estimate (called the pilot).These collective limitations motivated the development of new algorithm presented in [1].In [1], we applied our proposed algorithm only on single modality (T1) real MR data sets.In this paper we evaluate the performance of our proposed algorithm on multimodal MR data sets obtained from the MRBrainS challenge.The rest of the paper is organized as follows.In section 2 we summarize the adaptive mean shift algorithm.In section 3 we present the Bayesian approach for adaptive bandwidth estimation.The proposed segmentation algorithm is presented in section 4. Finally, the experimental results are presented in section 5.

Adaptive mean shift
The adaptive kernel density estimate of the underlying multivariate probability function at point x is given by fK where h i > 0 is the adaptive bandwidth and k : [0, 1] → R is the kernel profile of the spherically symmetrical kernel K with bounded support defined as and c k,d is a normalizing constant ensuring that the kernel K integrates to 1.In this work, the Gaussian kernel is applied which is defined as The derivative of the adaptive density estimator in eq. 1 leads to where g(x) = −k (x).The right-most factor (in square brackets) in eq. 4 is called the mean shift vector.It points toward the direction of maximum increase in density and is defined as where G represents the kernel and defined as Adaptive mean shift is an iterative hill climbing procedure in which the kernel G starts from an initial position y 1 and moves towards the position closer to the higher dense region.This process is continued until the position in higher dense region is achieved which represents a mode of the density.The update rule of kernel position is given by where {y j } j=1,2,... represents the successive locations of the kernel G.The feature points that converge to the same mode constitute a cluster.In order to perform the clustering using spatial information of voxel with its range (intensity) value, a joint spatial-range domain kernel K hs,hi (x) is used.It is defined as a product of spatial and range domain kernels and is given by where x s represents a vector of voxels spatial coordinates, x r represents a vector of voxels range (intensity) values and h s and h i are their corresponding kernel bandwidths, and C is a normalization constant.

Bayesian approach for adaptive bandwidth estimation
Our proposed approach is based on a novel variation on the Bayesian approach initially proposed in [9] for the global bandwidth estimation of the kernel.In our variation we used this idea locally for adaptive bandwidth estimation.The bandwidth is modeled by the a posteriori probability density function p(s|x) of local data spread s given the feature point x.Let M < n (total number of feature points) be the number of nearest neighborhoods to a feature point x.We can then define the pseudolikelihood where P (s|x Mj ) is the probability of local data spread s depending on the nearest neighborhood samples to x Mj and {M j | j = 1, 2, ...N } is the set of N such neighborhoods of various sizes.The evaluation of probabilities over the entire range of number of neighborhoods M j is given as Bayes rule yields where P (x Mj |M j ) is the probability of the feature point x Mj given the M j nearest neighborhood.Hereinafter P (M j ) is considered to have uniform distribution on the interval [M 1 M 2 ].Several values are selected for M j in this interval according to For a given M j the local variance s j is computed as where x (i,l) is the l-th nearest neighbor to the feature point x i .The distribution of variances is modeled as the Gamma distribution defined as where α and β define the shape and the scale of the Gamma distribution, respectively.These parameters are estimated using the maximum likelihood approach.Finally the adaptive bandwidth is estimated by the product of these parameters, identically the mean of Gamma distribution, i.e. ĥ( 4 Proposed segmentation algorithm The pre-processing steps includes: (1) extraction of brain mask which is done by employing the BET2 [2] tool on a T1 weighted image with parameter settings; fractional intensity threshold=0.2and threshold gradient=0.05.These parameters optimal values are determined by using the 5 training data sets.The resultant brain mask is then used to extract the brain of multimodal MR images.(2) Normalization of intensity values of each modality to the interval [0 1] which is done by applying the linear histogram stretching.The proposed algorithm BAMS involves the following steps for segmenting the brain into WM, GM and CSF. 1. BAMS takes feature vectors x i as input, which represent the thick slice multimodal MR brain images.The modes or clusters of the feature points are determined using the eq.7.The clustering is done in the joint spatial-range domain using the joint-domain kernel defined in eq.8.The spatial bandwidth is set as h s = 9 and the adaptive range bandwidth h i is estimated by employing our proposed Bayesian approach described in section 3 with parameters settings M 1 = 100, M 2 = 750 and N = 10.These parameters optimal values and the spatial bandwidth optimal value are determined empirically by using the 5 training datasets.The output of BAMS is a large number of modes or clusters.2. Finally the three tissues are obtained by merging the clusters (output from step 1) using the voxel weighted k-means clustering algorithm [4].Cluster initialization is performed making use of prior information of tissue intensity ranges in the multimodal MR images.The schematic procedure of the proposed algorithm is shown in Fig. 1.

Results
Our proposed unsupervised algorithm BAMS is applied on 12 test data sets, provided from the MRBrainS challenge.Each data set consists of thick slice multimodal (T1-weighted, T1-weighted inversion recovery (T1 IR)and FLAIR) MR images with voxel resolution 0.958mm×0.958mm×3mm.The segmentation performance of BAMS for each data set is evaluated relative to manual segmentation.The evaluation is performed by the MRBrainS challenge organizers.The average segmentation performance of proposed method BAMS is presented in

Table 1 .
The proposed algorithm BAMS is implemented in MATLAB (R2010a) and it takes approximately one and a half hours for segmenting the multimodal MR images on a single desktop PC running 64-bit Ubuntu 10.04 with 8 GiB ram memory and 8 Intel Core i7 CPU 870 at 2.93 GHz.

Table 1 .
Average segmentation performance of proposed algorithm BAMS in terms of Dice, modified Hausdorff distance (Mod.HD) and absolute volume difference (Abs.volume diff.)