MAP–Based Framework for Segmentation of MR Brain Images Based on Visual Appearance and Prior Shape

,


Introduction
Accurate delineation of the brain tissues from Magnetic Resonance Imaging (MRI) is an essential step used in many clinical applications in diagnosis, therapy evaluation, human brain mapping, and neuroscience [1].However, segmentation approaches are subject to multiple challenges stemming from image noise, image inhomogeneities, image artifacts such as partial volume effect, and discontinuities of boundaries due to similar visual appearance of adjacent brain structures.A variety of segmentation techniques have been developed to address these challenges.Brain MR segmentation methods can be classified into three main categories: probabilistic and statistical-based, atlas-based, and deformable model-based techniques.
Statistical-based approaches depend on building a prior models to describe signal distributions of each brain sturcutre.Anbeek et al. [2] presented an automated method for the segmentation of different brain structures (White Matter (WM), Central Gray Matter (CEGM), Cortical Gray Matter (COGM), and CerebroSpinal Fluid (CSF)) using T2-weighted and Inversion Recovery (IR) MRI of the neonatal brains.Probability maps for each brain tissue were manually constructed to segment each tissue class.Then, a K Nearest Neighbor (KNN) classifier was employed using voxel intensity values and voxel coordinates as the classification features.The obtained classes were combined into a multi-label segmentation.A similar approach was proposed by Xue et al. [3].Their method employed a parametric (Gaussian) density estimation using an Expectation-Maximization (EM) algorithm.The information from a Markov random field prior was used to increase the spatial homogeneity limitations of the MR images.They developed a technique for eliminating partial volume averaging effects by predicting the misclassification (e.g.CSF and GM "average" into an intensity similar to WM).Song et al. [4] proposed a Probabilistic Neural Network (PNN) framework for MRI brain tissue segmentation.The probabilistic density functions of the brain tissues were estimated from reference vectors generated by the Self-Organizing Map (SOM) neural network.In order to reduce the partial volume effects, they used a Weighted Probabilistic Neural Network (WPNN) that added weighting factors in the pattern of the summation layer.Finally, a supervised probabilistic classification based on a Bayesian rule was applied for soft labeling.
To overcome signal inhomogeneity and overlaps between signal distributions of different brain structures, atlas-based approaches have emerged as powerful segmentation tools.These approaches are based on a priori knowledge about brain structures, and treat the segmentation problem as a registration task.Morin et al. [5] presented an atlas-based segmentation framework based on random walks that combines registration and labeling propagation steps.They used a generative model to provide pixel label probabilities to give a greater impact on the segmentation for high-confidence labels.To find the matches from target images to atlas images, they used the Affine-Scale Invariant Feature Transform (ASIFT) [6] and Speeded Up Robust Features (SURF) [7] registration techniques.In order to avoid segmentation errors produced by registration imperfec-tion, Lötjönen et al. [8] introduced an optimized pipeline for multi-atlas brain MRI segmentation.They introduced two approaches that combine multi-atlas segmentation and intensity modeling based on using EM and graph cuts for optimization.First, they register all atlases to the target data and a majority voting is applied to predict the segmentation of the target image.Then, the intensity modeling is used as a post-processing step to improve the segmentation.Artaechevarria et al. [9] proposed a generalized local weighting voting scheme in which the fusion weights are adapted for each voxel based on local estimation of the segmentation performance.The local weighting voting outperforms traditional global strategies that estimate a single value for the segmentation accuracy for the whole image.Lijn et al. [10] introduced a segmentation method based on the combination of structures' location information and appearance models.They generated a spatial probability map, obtained from multiple atlas-target image registrations, to implement the spatial model.The tissue appearance was modeled by a KNN classifier based on Gaussian scale-space features.Then, a Bayesian framework was used to combine both spatial and appearance models and graph-cut approach was used for optimization.Sabuncu et al. [11] proposed an automated, label fusion segmentation technique.In order to capture greater inter-subject anatomical variability, each training data set was individually coregistered to the test data set.Then, a nonparametric probabilistic model was employed to fuse the training labels to compute the final segmentation.
In order to obtain continuous segmentation of brain structures, deformable boundaries have also been recognized as more accurate segmentation techniques of MR brain tissues.Angelini et al. [12] introduced a multi-phase level set framework for the automated segmentation of brain MRIs.The segmentation of the brain tissues (WM, GM and CSF) was solely based on homogeneity (average grey level) measures.To avoid the need for any prior information and to speed up numerical calculation, a random seed for initialization of the deformable boundaries was used.Colliot et al. [13] proposed a deformable model-based approach that used spatial constraints, represented as fuzzy subsets of the 3D image space as an external force to control the boundary evolution.To avoid manual selection of the model parameters, a training step was required to estimate the spatial constraints parameters.Huang et al. [14] introduced an automated, hybrid deformable model framework that integrates both image edge geometry and voxel statistics features to regularize the convergence of the deformable contour.Wang et al. [15] proposed a muti-phase level set framework to segment brain MR images with intensity inhomogeneity.They modeled the local image intensities using Gaussian distributions with different means and variances.Then, a variational approach minimizes an energy function to compute the means and variances that will guide the contour evolution towards the target boundaries.Ciofolo et al. [16] developed an automated framework based on level sets for the simultaneous segmentation of multiple structures from brain MR images.The evolution of each level set was driven by a fuzzy decision system that combines three factors: intensity distribution of the 3D MR volume, the relative position of the evolving contours, and a priori knowledge provided by an anatomical atlas.
In summary, the segmentation of MR brain images has been the subject of extensive research.However, the above brief overview shows the following limitations of the existing literature approaches.While known statistical-based techniques are easy to implement, they depend only on predefined probability models that cannot fit all of the possible real data distributions.This is due to the fact that actual intensity distributions of brain structures are greatly affected by several factors, such as the unique patient, scanner, and scanning parameters.The problem of single atlas-based segmentation is that it is heavily registration dependent and does not cover all of the anatomical variability.While multi-atlas methods try to cover these problems, they are still challenged by atlas selection, combination, and the associated heavy computation time.The accuracy of deformable model-based techniques is based on the accurate design of the guiding forces (statistical, geometric, etc.).
To overcome these limitations, we propose an automated MAP-based approach aimed at unsupervised segmentation of different brain tissues from T1weighted MRI.The proposed segmentation approach is based on the integration of statistical approaches, namely a probabilistic shape prior, first-order intensity model, and second-order appearance model.These three features are integrated into a two-level joint Markov-Gibbs Random Field (MGRF) model of T1-MR brain images.In this paper, we focus on the accurate identification of the spatial interactions between the brain voxels and the first-order visual appearance descriptor for the brain tissues.

Joint MGRF Model of T1-weighted MR Brain Images
Let Q = {0, . . ., Q − 1} and L = {1, . . ., L} denote sets of gray levels q and region labels L, respectively.Let R denote a 3D arithmetic (x, y, z)-lattice that supports a given grayscale image g : R → Q to be segmented and its goal labled region map m : R → L. The 3D T1-weighted MR brain images, g, being coaligned to the 3D training data, and its map, m, are described with the following joint probability model P (g, m) = P (g|m)P (m) (1) The above joint model combines a 3D second-order MGRF, P (m), of region labels with the shape prior of each brain structure and a conditionally independent random field.The joint conditional distribution of image intensities given the map is P (g|m) = (x,y,z)∈R p(g x,y,z |m x,y,z ).The map model P (m) = P sp (m)P h (m) has two parts: (i) a shape prior probability P sp (m), and (ii) a 2 nd -order MGRF model P h (m) of a spatially homogeneous region map m for the image g, which is aligned to the training set.In this work we focus on accurate identification of the spatial interaction between the brain voxels (P (m)) and the 1 st -order visual appearance descriptor for the brain tissues (P (g|m)) as shown in Fig. 1.

Shape Model
To enhance the segmentation accuracy, expected shapes of each brain label are constrained with a probabilistic shape prior.A training set of images, collected for different subjects, are co-aligned by 3D affine transformations with 12 degrees of freedom in a way that maximizes their Mutual Information (MI) [17].The shape prior is a spatially variant independent random field of region labels P sp (m) = (x,y,z)∈R p sp:x,y,z (m x,y,z ) for the co-aligned, manually segmented five data sets that are provided to our team by the MRBrainS13 challenge team, specified by voxel-wise empirical probabilities for each brain label (p sp:x,y,z (l), l ∈ {1, . . ., L}).
Our framework exploits four prior shapes (built at the learning stage) for the WM, GM, CSF, and other brain structures (excluding the background).For the training phase, we use five manually segmented data sets provided to us by the MRBrainS13 challenge team to create the probabilistic map for each label.During the testing phase, each data to be segmented is globally registered with the set of training data that have been used to create the prior shape model for each brain label.

1 st -Order Visual Appearance Descriptor
In addition to the learned prior shape descriptor, we implement a 1 st -order visual appearance descriptor of each brain label.During segmentation of a data set, this visual appearance descriptor is roughly taken into account by approximating the 1D empirical marginal gray level distributions of the T1-weighted MR brain images with a Linear Combinations of Discrete Gaussians (LCDG).This LCDG model is a modified version of our previous linear combination of continuous Gaussian probabilistic model [18,19].This approximation adapts the segmentation to the changing appearance, such as non-linear intensity variations caused by patient weight and data acquisition systems (e.g.scanner types and scanning parameters).The LCDG models the empirical distribution of each brain label more accurately than a conventional mixture of only positive Gaussians.This yields a better initial region map that is formed by the voxel-wise classification of the gray values in the images.
Let F = f (q) : q ∈ Q; q∈Q f (q) = 1 be an empirical marginal probability distribution of gray levels, q, for the 3D T1-weighted MR brain image, g.Each distribution has a known number, L, of dominant modes related to the regions of interest (in our case, L = 4).To segment an image by separating these individual dominant modes, the individual gray level distributions that are associated with each mode must be estimated from F. Unlike conventional modeling with a mixture of Gaussians [20], or another simple distribution type [21] where there is only one distribution per dominant mode, we approximate F j much more closely using the LCDG.The LCDG is partitioned for the whole image into multiple LCDG submodels that each relate to a unique dominant mode (in our case, WM, GM, CSF, and other brain structures).The Discrete Gaussian (DG) components Ψ θ = (ψ(q|θ) : q ∈ Q) then integrate a continuous Gaussian with parameters θ = (µ, σ 2 ); where µ and σ 2 are the mean and the variance, respectively, over successive gray level intervals.The LCDG with C p positive and C n negative components, such that C p ≥ L, is defined as follows [18,19]: with the non-negative weights w = [w p,., w n,.] such that In order to precisely estimate the parameters of the LCDG model, including the numbers of positive and negative components, EM-based techniques, namely those introduced in [18] for approximation of a probability density model when using a linear combination of Gaussians, have been adapted to this LCDG model.
One of the advantages of the proposed approach is that it addresses the inhomogeneities in the transition regions of different brain structures.At these transition locations, there is typically a spatial inhomogeneity in the signal intensity [23].Importantly, in addition to signal intensity, our approach uses shape and spatial interaction components.This provides additional information that aids in overcoming the variability that often results when solely taking signal intensities into account.The proposed step-wise segmentation algorithm as demonstrated in Fig. 3 is summarized in Algorithm 1.
Algorithm 1: Key Steps for the Proposed Segmentation Approach 1. Use Brain Extraction Tool (BET) software [24,25] to remove the skull from T1-weighted MR brain images.2. Approximate the marginal intensity distribution P (g) of the T1-weighted MR brain image using LCDG with four dominant modes.3. Form an initial region map m using the marginal estimated density and prior shape of each brain label.4. Find the Gibbs potentials for the MGRF model from the initial map by using Eq. 3. 5. Improve the region map m using voxel-wise Bayes classifier.Fig. 2. Illustration of the proposed 3D neighborhood system.

Performance Evaluation of the Proposed Segmentation Approach
To evaluate the segmentation accuracy, we use three metrics, namely, the Dice Similarity Coefficient (DSC), the Hausdorff metric, and the percentage volume difference between the segmentation and the ground truth [26].The DSC characterizes the agreement between the segmented and ground truth regions [26]: ) where TP, FP, and FN denote the True Positive, False Positive, and False Negative segmentation results, respectively.For a segmented region, SR, and its ground truth, GT, TP = |SR ∩ GT| is the area of their overlap, i.e., the number of the common points in SR and GT; FP = |SR − SR ∩ GT| is the number of points in the difference between SR and TP, and FN = |GT − SR ∩ GT| is the number of points in the difference between GT and TP.The closer the DSC to "1", the better the segmentation.In addition to the DSC metric, a volumetric metric is used to assess the accuracy of the segmentation.The percentage volume difference is calculated between the segmentation and the ground truth and is used to test the segmentation accuracy.The smaller the volume difference, the better the segmentation.
Another metric, which is used to evaluate the segmentation accuracy, is the Hausdorff distance.Hausdorff distance from a set A 1 to a set A 2 is defined as the maximum distance of the set A 1 to the nearest point in the set A 2 [26]: where c and e are points of sets A 1 and A 2 , respectively, and d(c, e) is Euclidean distance between these points.The bidirectional Hausdorff distance, H Bi (GT, SR), between the segmented region SR and its ground truth GT is defined as: H Bi (GT, SR) = max{H(GT, SR), H(SR, GT)}.In this work, we use the 95-percentile bidirectional Hausdorff distance as a metric that measures the segmentation accuracy.The smaller the distance, the better the segmentation.The ideal case with perfect segmentation is when the bidirectional Hausdorff distance is equal to 0.

Experimental Results and Conclusions
To assess robustness and computational performance, the proposed segmentation techniques have been tested on twelve thick T1-weighted MR brain images.This data has been acquired at the UMC Utrecht (the Netherlands) from patients with diabetes and their matched controls (having an increased cardiovascular risk) with varying degrees of atrophy and white matter lesions (age>50).These 12 data sets were acquired with a 3T MRI scanner with voxel size 0.958mm × 0.958mm × 3.0mm.Our goal is to classify T1-weighted MR brain images to the following four classes: WM, GM, CSF, and other.By grayness, some brain tissues such as the brain stem and cerebellum are very close to WM and GM.Therefore, the segmentation cannot be based only on an image's 1 st -order visual appearance descriptor but has to also account for shape and 3D spatial relationships between the region labels in order to accurately segment each label.Empirical gray level density for the T1-weighted MR brain images, after removing the brain skull, and its estimated four dominant densities are shown in Fig. 4. Figure 5 demonstrates a step by step marginal density estimation for each class (1 st -order visual appearance descriptor) using the LCDG model.As demonstrated in Algorithm 1, the first step is to remove the skull using BET software [24,25] as shown in Fig. 6(b).The second step is to get the initial region map by using the estimated LCDG models for each class (see Fig. 5(f)) and our prior shape as demonstrated in Fig. 6(c).Finally, the initial map is refined by using the proposed three descriptors (1 st -order visual appearance, 3D spatial MGRF, and prior shape models) to get the final segmentation as shown in Fig. 6(d).Figure 7 shows the 3D results of our segmentation approach.In order to measure the accuracy of our segmentation approach, we provided our segmentation results of the 12 T1-weighted MR brain images to MRBrainS13 challenge team who have the ground truth that has been validated by three experts.The MR-BrainS13 challenge team evaluated the accuracy of our segmentation approach using the following three performance metrics: 1) DSC, 2) 95% Hausdorff metric, and 3) percentage absolute volume difference.As demonstrated in Table 7, the DSC for segmentation of the whole brain (GM + WM) from T1-weighted MR brain images only is 94.87% which confirms the high accuracy of the proposed segmentation techniques.Please see Table 7 for more details.
Our experiments show that the proposed accurate identification of the Joint MGRF model demonstrates promising results in segmenting the brain (GM + WM) from T1-weighted MR brain images.Our present implementation in the C++ programming language on a Dell precession T7500 workstation with an Intel quad-core processor (3.33 GHz each) with 48 GB of memory and a 1.5 TB hard drive with RAID technology takes about 5.78 ± 0.54 sec for processing 48 T1-MR images of size 240×240 pixels each, i.e about 0.12 sec per image.Table 1.Accuracy of our segmentation approach using DSC, 95% Hausdorff metric, and percentage absolute volume difference.Note that these evaluation results have been done by MRBrainS13 challenge team using the ground truth images that have been validated by three experts.

Fig. 1 .
Fig. 1.Illustration of the proposed joint Markov-Gibbs model of T1-weighted MR brain images.