PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Neuroimage. Author manuscript; available in PMC 2010 July 1.
Published in final edited form as:
PMCID: PMC2748869
NIHMSID: NIHMS145710

Automated Segmentation of Mouse Brain Images Using Extended MRF

Abstract

We introduce an automated segmentation method, extended Markov Random Field (eMRF) to classify 21 neuroanatomical structures of mouse brain based on three dimensional (3D) magnetic resonance imaging (MRI). The image data are multispectral: T2-weighted, proton density-weighted, diffusion x, y and z weighted. Earlier research (Ali et al., 2005) successfully explored the use of MRF for mouse brain segmentation. In this research, we study the use of information generated from Support Vector Machine (SVM) to represent the probabilistic information. Since SVM in general has a stronger discriminative power than the Gaussian likelihood method and is able to handle nonlinear classification problems, integrating SVM into MRF improved the classification accuracy. The eMRF employs the posterior probability distribution obtained from SVM to generate a classification based on the MR intensity. Secondly eMRF introduces a new potential function based on location information. Third, to maximize the classification performance eMRF uses the contribution weights optimally determined for each of the three potential functions: observation, location and contextual functions, which are traditionally equally weighted. We use the voxel overlap percentage and volume difference percentage to evaluate the accuracy of eMRF segmentation and compare the algorithm with three other segmentation methods – mixed ratio sampling SVM (MRS-SVM), atlas-based segmentation and MRF. Validation using classification accuracy indices between automatically segmented and manually traced data shows that eMRF outperforms other methods.

Keywords: Automated segmentation, Data mining, Magnetic resonance microscopy, Markov Random Field, Mouse brain, Support Vector Machine

Introduction

MRI is often the imaging modality of choice for noninvasive characterization of brain anatomy because of its excellent soft tissue contrast. Its detailed resolution allows the investigation of normal anatomical variability bounds, as well as the quantization of volumetric changes accompanying neurological conditions. A prerequisite for such studies is an accurate segmentation of the brain; therefore many studies have focused on tissue classification into white matter, gray matter and cerebrospinal fluid (CSF), as well as anatomical structure segmentation. Several successful methods for tissue segmentation include statistical classification and geometry-driven segmentation (Ballester et al., 2000), statistical pattern recognition methods based on a finite mixture model (Andersen et al., 2002), expectation maximization algorithm (Kovacevic et al., 2002), artificial neural network (Reddick et al., 1997), hidden Markov Random Field (Zhang et al., 2001), and region-based level-set and fuzzy classification (Suri, 2001). Compared to anatomical structure segmentation, tissue segmentation is a relatively easier task (Heckmann et al., 2006). This is because the MR signals are in general distinctive enough to have the brain tissues classified into white matter, gray matter and CSF while the anatomical structures can be composed of more than one tissue type, and have ambiguous contours. Nevertheless, segmenting neuroanatomical structures has recently attracted considerable attention since it can provide stronger foundation to help in the early diagnosis of a variety of neurodegenerative disorders (Fischl et al., 2002).

Automated methods for segmenting the brain in anatomically distinct regions can rely on either a single imaging protocol, or multispectral data. For example, Duta et al. (1998) have segmented T1-weighted MR images of the human brain into 10 neuroanatomical structures using active shape models. Fischl et al., (2002) accomplished an automated labeling of each voxel in the MR human brain images into 37 neuroanatomical structures using MRF. Heckemann et al. (2006) segmented T1-weighted human MR images into 67 neuroanatomical structures using nonrigid registration with free-form deformations, and combined label propagation and decision fusion for classification. Multispectral imaging was used as the basis of segmentation for example by Zavaljevski et al. (2000) and Amato et al. (2003). Zavaljevski et al. (2000) used Gaussian MRF and maximum likelihood method on multi-parameter MR images (including T1, T2, proton density, Gd+ T1 and perfusion imaging) to segment the human brain into 15 neuroanatomical structures. Amato et al. (2003) used independent component analysis (ICA) and nonparametric discriminant analysis to segment the human brain into nine classes.

Developments in small animal imaging capabilities have led to increased attention being given to the problem of mouse brain segmentation, since mice provide excellent models for the anatomy, physiology, and neurological conditions in humans, with whom they share more than 90 percent of the gene structure. Unfortunately, the methods for human brain anatomical structure segmentation cannot directly be applied with the same success to the mouse brain. This is probably due to the facts that (1) a normal adult mouse brain is approximately 3000 times smaller than an adult human brains – thus the spatial resolution needs to be scaled accordingly from the 1-mm3 voxels commonly used in the study for human brains to voxel volumes < 0.001 mm3; (2) the inherent MR contrast of the mouse brain is relatively low compared to the human brain.

Research on mouse brain image segmentation is still limited thus far and the major directions of research include the atlas based approach and the probabilistic information based approach. The atlas based approach is to create a neuroanatomical atlas, also called reference atlas, using training brains and nonlinearly register the atlas to test brains with the scope of labeling each voxel of the test brains (Mazoyer et al., 2002; Kovacevic et al., 2005; Ma et al., 2005; Bock et al., 2006). For example, Ma et al. (2005) segmented T2*-weighted magnetic resonance microscopy (MRM) images of 10 adult male formalin-fixed, excised C57BL/6J mouse brains into 20 anatomical structures. They chose one mouse brain image out of the 10 mouse brains as a reference atlas. The rest nine testing images were aligned to the reference atlas using a six-parameter rigid-body transformation and then the reference atlas were mapped to the nine testing brain images using a nonlinear registration algorithm. The registration step yields a vector field, called a deformation field, which contains information on the magnitude and direction required to deform a point in the reference atlas to the appropriate point in the testing brain images. Using the deformation fields, the labeling of the reference atlas is transformed to the testing images to predict the labeling of the testing images. Kovacevic et al. (2005) built a variational atlas using the MR images of nine 129S1/SvImJ mouse brains. The MR images of the genetically identical mouse brains are aligned and nonlinearly registered to create the variational atlas comprised of an unbiased average brain. The probabilistic information based approach uses probabilistic information extracted from training datasets of MR images for the segmentation. Ali et al. (2005) incorporates multi-spectral MR intensity information, location information and contextual information of neuroanatomical structures for segmentation of MRM images of the C57BL/6J mouse brain into 21 neuroanatomical structures using the MRF method. They modeled the location information as a prior probability of occurrence of a neuroanatomical structure at a location in the 3D MRM images and the contextual information as a prior probability of pairing of two neuroanatomical structures.

MRF has been widely used in many image segmentation problems and image reconstruction problems such as denoising and deconvolution (Geman and Geman, 1984). It provides a mathematical formulation for modeling local spatial relationships between classes. In the brain image segmentation problem, the probability of a label for a voxel is expressed as a combination of two potential functions: one is based on the MR intensity information and the other is based on contextual information, such as the labels of voxels in a predefined neighborhood around the voxel under study. Bae et al. (2008) developed an enhanced SVM model, called Mix-Ratio sampling-based SVM (MRS-SVM), using the multispectral MR intensity information and voxel location information as input features. The MRS-SVM provided a comparable classification performance with the probabilistic information based approach developed in Ali et al. (2005). Furthermore, Bae’s study also suggested that compared to MRF, the MRS-SVM outperforms for larger structures, but underperforms for smaller structures.

Based on these studies which suggest that integrating a powerful SVM classifier into the probabilistic information based approach may improve the overall accuracy of mouse brain segmentation, we introduce a novel automated method for brain segmentation, called extended MRF (eMRF). In this method (eMRF), we use the voxel location information by adding a location potential function. Other than using Gaussian probability distribution to construct a potential function based on MR intensity, the SVM classifier is used to model the potential function of MR intensity. We assess the accuracy of the automated brain segmentation method in a population of five adult C57BL6 mice, imaged using multiple (five) MR protocols. All 21 neuroanatomical structures are segmented, and the accuracy is evaluated using two metrics including the volume overlap percentage (VOP) and volume difference percentage (VDP). The use of these two metrics allows us to compare the accuracy of our method with the atlas-based segmentation, MRF and MRS-SVM methods.

Materials and Methods

Subjects and Magnetic Resonance Microscopy data

The MRM images used in this work were provided by the Center for In Vivo Microscopy in Duke University Medical Center and they were previously used in Ali et al. (2005) as well. Five formalin-fixed C57BL/6J male mice of approximately 9 weeks in age were used. The MRM image acquisition consisted of isotropic 3D T2-weighted, proton density-weighted, diffusion x, y and z weighted scans. Image acquisition parameters for all acquisition protocols include the field of view of 12×12×24 mm and matrix size of 128×128×256. All protocols used the same flip angle of 135 degrees and 2 excitations (NEX). Specific to the PD image, TE/TR was 5/400 ms and bandwidth was 62.5 KHz; for the T2 weighted image, TE/TR was 30/400 ms and bandwidth was 62.5 kHz bandwidth; and for the three diffusion scans, TE/TR was 15.52/400 ms, and bandwidth was 16.25 MHz. The Stejskal Tanner sequence was used for the acquisition of the diffusion-weighted scans. Bipolar diffusion gradients of 70 G/cm with pulse duration of 5 ms and inter-pulse interval of 8.56 ms were applied along the three axes and the effective b value of 2600 s/mm2 was used. A 9-parameter affine registration which accounts for scaling, rotation and translation was applied to each mouse brain, to bring it into a common space. Manual labeling of 21 neuroanatomical structures was done by two experts using T2-weighted datasets of the five mouse brains. These manual labelings were regarded as true labeling for each voxel. Table 1 presents the list of the 21 neuroanatomical structures and abbreviations to be segmented in this work.

Table 1
List of 21 segmented structures and their abbreviations.

Markov Random Field Theory

MRF theory is a class of probability theory for modeling the spatial or contextual dependencies of physical phenomena. It has been used for brain image segmentation by modeling probabilistic distribution of the labeling of a voxel jointly with the consideration of the labels of a neighborhood of the voxel (Fischl et al., 2002; Ali et al., 2005). For simplicity, let us consider a 2D image represented by an m×n matrix, let X be a vector of signal strength and Y be the associated labeling vector, that is, X=(x11, x12, …, x1n, …, xm1, … xmn), and Y=(y11, y12, …, y1n, …, ym1, … ymn). We will rewrite yab as yi where i=1, 2, …, S, S=mn, Y is said to be a MRF on S with respect to a neighborhood N if and only if the following two conditions are satisfied:

P(Y)>0
(Eq. 1)

P(yi|yS{i})=P(yi|yNi)
(Eq. 2)

where S-{i} denotes the set difference, and Ni denotes the set of sites neighboring site i. Eq. 1 is called the positivity property and Eq. 2 is the Markovianity property, which states only neighboring labels (or clique) have direct interaction with each other. If these conditions are satisfied, the joint probability P(Y) of any random field is uniquely determined by its local conditional probabilities.

The equivalence between the MRF and Gibbs distribution (Hammersley-Clifford theorem) provides a mathematically efficient way of specifying the joint probability P(Y) of an MRF. The theorem states the joint probability P(Y) can be specified by the clique potential function Vc(Y) which can be defined by any appropriate potential function based on a specific system’s behavior. For more information about potential function, please refer to Li (2001). That is, the probability P(Y) can be equivalently specified by a Gibbs distribution as follows:

P(Y)=1Zexp[U(Y)]
(Eq. 3)

where Z=YΩYexp[U(Y)] is a normalizing constant, ΩY is the set of the all possible Y on S, and U(Y) is an energy function which is defined as a sum of clique potential Vc(Y) over all cliques c [set membership] C:

U(Y)=cCVc(Y)
(Eq. 4)

where a clique, c, is a set of points that are all neighbors of each other and C is a set of cliques, or the neighborhood of the clique under study. The value of Vc(Y) depends on a certain configuration of labels on the clique c.

For the image segmentation problem, we try to maximize a label’s posterior probability for given specific features, that is P(Y|X). With the assumption of feature independency the posterior probability can be formulated using Bayesian theorem as:

P(YX)P(XY)P(Y)=P(Y)i=1SP(xiyi)
(Eq. 5)

Eq. 5 can be rewritten as:

P(YX)=1Zexp[iSlogP(xiyi)+cCVc(Y)]
(Eq. 6)

Typically in MRF, a multivariate Gaussian distribution is used for P(X|Y) and the maximum likelihood estimation is performed to find the labeling based on MR signals. This model assumes that the relationship between the features and labels follows the Gaussian distribution. In some of image segmentation problems, this assumption is too restrictive to model the complex dependencies between the features and the labels (Lee et al., 2005). If some data mining techniques, such as SVM, are integrated with MRF, the overall segmentation performance will be improved due to their powerful class discrimination abilities even for the cases with complex relationship between the features and labels.

Support Vector Machines Theory

SVM (Vapnik, 1995) was initially designed for binary classification by constructing an optimal hyperplane which gives the maximum separation margin between two classes. Considering a training set of m samples (xi, yi), i=1, 2, …, m where xi [set membership] Rn and yi [set membership] {+1, −1}. Samples with yi=+1 belong to positive class while those with yi=−1 belong to negative class. SVM training involves finding the optimal hyperplane by solving the following optimization problem:

MinQP(w,b,ξ)=12wTw+Ci=1mξis.t.yi(w·Φ(xi)+b)1ξi,ξi0,i=1,,m.
(Eq. 7)

where w is the n dimensional vector, b is a bias term, ξ= {ξ1, …, ξm} and QP is the objective function of the prime problem. The non-negative slack variable (ξi) allows Eq. 7 to always yield feasible solutions even in a non-separable case. The penalty parameter (C) controls the trade-off between maximizing the class separation margin and minimizing the classification error. A larger C usually leads to higher training accuracy, but may over-fit the training data and cause the classifier un-robust. To enhance the linear separability, the mapping function (Φ(xi)) projects the samples into a higher-dimensional dot-product space called the feature space. Figure 1 shows the optimal hyperplane in solid line, which can be obtained by solving Eq. 7. The squares represent the samples from the positive class and the circles represent the samples from the negative class. The samples which satisfy the equality are called support vectors. In Figure 1, the samples represented as the filled squares and the filled circles are the support vectors.

Figure 1
SVM binary classification problem (adopted from Vapnik, 1995). The solid line between the two dashed lines is the optimal hyperplane. The squares represent the samples from the positive class and the circles represent the samples from the negative class. ...

Eq. 7 presents a constrained optimization problem. By introducing the non-negative Lagrangian multiplier αi and βi, it can be converted to an unconstrained problem as shown below:

MinQP(w,b,ξ,α,β)=12wTw+Ci=1mξii=1mαi(yi(wT·Φ(xi)+b)1+ξi)i=1mβiξi
(Eq. 8)

where α = {α1, …,αm} and β = {β1, …, βm}. Furthermore, by differentiating with respect to w, b and ξi and introducing the Karush-Kuhn-Tucker (KKT) condition, Eq. 8 is converted to the following Lagrangian dual problem:

MaxQD(α)=i=1mαi12i,j=1mαiαjyiyjΦ(xi)TΦ(xj)s.t.0αiCi=1,,mandi=1maiyi=0
(Eq. 9)

The optimal solution αi* for the dual problem determines the parameters w* and b* of the following optimal hyperplane, also known as SVM decision function:

f(x)=sign(w·Φ(xi)+b)=sign(i=1myiαiΦ(x)TΦ(xi)+b)=sign(i=1myiαiK(x,xi)+b)
(Eq. 10)

where K(xi, xj) is a kernel function defined as K(xi,xj)=Φ(xi)TΦ(xj). The kernel function performs the nonlinear mapping implicitly. We chose a Radial Basis Function (RBF) kernel, defined as:

K(xi,xj)=exp(γ||xixj||2),γ>0
(Eq. 11)

where γ in Eq. 11 is a parameter related to the span of an RBF kernel. The smaller the value is, the wider the kernel spans.

To extend the application of SVM for multiclass classification, a number of methods have been developed, which mainly fall in three categories: One-Against-All (OAA), One-Against-One (OAO) and All-At-Once (AAO). In OAA method, one SVM is trained with the positive class representing one class and the negative class representing the others. Therefore, it builds n different SVM models where n is the number of the classes. The idea of AAO is similar to that of OAA, but it determines n decision functions at once in one model, where the kth decision function separates the kth class from the other classes. In OAO method, a SVM is trained to classify the kth class and the lth class. Therefore, it constructs n(n−1)/2 SVM models. Hsu and Lin (2002) reported that the training time of OAO method is less than that of OAA or AAO method. OAO is more efficient on large datasets than OAA and AAO method. Thus, OAO method is used in this study. In OAO method, the following problem is to be solved:

MinQP(wkl,bkl,ξkl)=12(wkl)T(wkl)+Ci=1mξikls.t.(wkl)T·Φ(xi)+bkl1ξikl,ifyi=k(wkl)T·Φ(xi)+bkl1+ξikl,ifyi=lξik0,i=1,,mandk,l=1,,n.
(Eq. 12)

eMRF

As discussed earlier in Eq. 6, traditional MRF mainly focuses on MR intensity as well as the contextual relationship. Research on brain segmentation has concluded that beside of MR intensity and contextual relationship with neighboring structures, location within the brain also plays an important role (Fischl et al., 2002; Ali et al., 2005, Bae et al., 2008). To incorporate the three different types of information into an image segmentation model, Bayesian theorem is used. Consider a 3D MR image represented by an m×n×o matrix, with its associated multivariate intensity vector X=(x111, x112, …, x11o, …, x1n1, … xm11, …, xmno), a location vector L=(l111, l112, …, l11o, …, l1n1, … lm11, …, lmno), and class label configuration Y=(y111, y112, …, y11o, …, yl1n1, … ym11, …, ymno). We rewrite yabc as yi where i=1, 2, …, S, and S=mno. The posterior probability of having a label configuration Y given a multivariate intensity vector X and a location vector L is formulated as follows:

P(YX,L)=P(X,LY)P(Y)P(X,L)
(Eq. 13)

Assuming X and L independent from each other yields following expression:

P(YX,L)=P(XY)P(LY)P(Y)P(X)P(L)P(YX)P(X)P(YL)P(L)P(Y)P(X)P(L)=P(YX)P(YL)P(Y)
(Eq. 14)

If we make logarithmic transformation to Eq. 14, we obtain

logP(YX,L)logP(YX)+logP(YL)+logP(Y)
(Eq. 15)

Each term of the right hand side of Eq. 15 can be regarded as the contribution to labeling from MR signal intensities, voxel coordinates, and prior belief of label, which incorporates the contextual information into the labeling decision. The function can be modified as following:

P(YX,L)exp[w1iSlogAi(yi,xi)+w2iSlogBi(yi,li)+w3iSVi(yi,yNi)]s.t.w1+w2+w3=1
(Eq. 16)

where w1, w2 and w3 are model parameters which indicate the degree of contribution of each term to the posterior probability P(Y|X,L), and Ni denotes the neighboring sites of site i.

The first term, Ai(yi, xi), in Eq. 16 is the observation potential function that models the MR intensity information. EMRF model employs SVM in order to take advantage of the excellent discriminative ability of SVM. Since SVM decision function gives as output the distance from an instance to the optimal separation hyperplane, Platt (2000) proposed a method for mapping the SVM outputs into posterior probability by applying a sigmoid function whose parameters are estimated from the training process. The observation potential function is defined as follows, for voxel i:

Ai(yi,xi)=P(yi=kxi)=11+exp(αfk(xi)+β)
(Eq. 17)

where fk(xi) is the SVM decision function for class k, α and β are the parameters determined from the training data. P(yi=k|xi) is the likelihood that the label yi is labeled as class k given the MR intensity xi. When applying SVM to MR image segmentation, the overlapping of MR signals is a critical issue. Adding the three coordinates (x, y and z coordinates of each data point in a 3D MR image) as features into SVM classifier helps in classification by relieving the class overlapping problem of MR intensities between different structures (Bae et al., 2008). In one of our preliminary experiments, it is interesting to find that using two MR protocols (T2-weighted and proton density-weighted) along with coordinate features yield better segmentation performance than using all the five MR protocols 1. Our preliminary result also indicates that adding more MR contrasts does not aid the classifier because it increases the overlapping among MR signals and makes the dataset noisier. Thus, the intensity feature vector X forms a five dimensional vector, consisting of the two MR intensities (T2-weighted and proton-density (PD) weighted acquisition) and the three coordinate features.

The second term Bi(yi, li) in Eq. 16, is the location potential function. Fischl et al. (2002) point out that the number of possible neuroanatomical structures at a given location in a brain atlas becomes small as registration become accurate. Therefore, the location of a voxel in a 3D image after registration is critical for classification of the voxel into the neuroanatomical classes. The location potential function is computed as follows:

Bi(yi,li)=P(yi=kli=r)=#ofvoxelslabeledaskatlocationr#ofvoxelsatlocationr
(Eq. 18)

where P(yi=k|li=r) is the probability that a voxel’s label yi is predicted as class k given that the location li of the voxel is r. The denominator in Eq. 18 is equal to the number of mice used in the training set. This location information is similar to the apriori probability used in human brain segmentation study (Ashburner and Friston, 1997) as implemented in SPM.

The third term Vi(yi,yNi) in Eq. 16 is the contextual potential function which models the contextual information using MRF. Based on the MRF theory, the prior probability of having a label at a given site i is determined by the label configuration of the neighborhood of the site i. The contextual potential function for site i will have a higher value as the number of neighbors that have the same label increases. It is defined as

Vi(yi,yNi)=jNiδ(yi,yj)n(Ni)whereδ(yi,yj)={1ifyi=yj0ifyiyj
(Eq. 19)

where n(Ni) is the number of voxels in a neighborhood of site i. We use a first order neighborhood system as a clique, which consists of the adjacent six voxels in the four cardinal directions in a plane and the front and back directions through the plane.

In fact, other than observation and contextual potential function, eMRF introduces location potential function to the problem. Secondly, noticing that each potential function can independently classify a voxel, eMRF model integrates them together and assigns weight optimally to each of them based on their classification performance from the training set. The weight parameter indicates the contributions from the three different potential functions to predict the labeling. The optimal weight parameter is determined by grid search on several possible sets of w1, w2 and w3 using cross-validation. The Y that maximizes the posterior probability P(Y|X,L) corresponds to the most likely label given the information of X and L. This is known as the maximum a posterior (MAP) solution which is well known in the machine vision literatures (Li, 2001). Finding the optimal solution of the joint probability of a MRF Y=maxYP(YX,L) is very difficult because of the complexity caused by the interaction among multiple labels. A local search method called iterated conditional modes (ICM) proposed by Besag (1986) is used in this study to locate the optimal solutions. The ICM algorithm is to maximize local conditional probabilities sequentially by using the greedy search in the iterative local maximization. It is expressed as

yi=argmaxyiYP(yixi,li)
(Eq. 20)

Given the data xi and li, the algorithm sequentially updates yi(t) into yi(t+1) by switching the different labels to maximizing P(yi|xi, li) for every site i in turn. The algorithm terminates when no more labels are changed. In the ICM algorithm, how to set the initial estimator y(0) is very important. We use the MAP solution based on only the location information as the initial estimator of the ICM algorithm, i.e.,

yi(0)=argmaxyiYP(yili)
(Eq. 21)

Performance Measurements

The VOP and VDP are used to measure the performance of the proposed automated segmentation procedure (Fischl et al., 2002; Ali et al., 2005, Hsu and Lin, 2002). The VOP and VDP are calculated by comparing the automated segmentation with the true voxel labeling (from the manual segmentation). Denote LA and LM as labeling of the structure i by automated segmentation and manual segmentation respectively, and V(L) as the volume of the labeling. The volume overlap percentage for class i is defined as

VOPi(LA,LM)=V(LALM)(V(LA)+V(LM))/2×100%

This performance index is the larger the better. When all the labels from the automated and manual segmentation coincide, VOPi(LA, LM) is 100%. VOP is very sensitive to the spatial difference of the two labelings, because a slight difference in spatial location of the two labelings can cause significant decreases in the numerator of VOPi(LA, LM).

The VDP for class i is used for quantifying the volume difference of the structures delineated by the two segmentations, and it is defined as

VDPi(LA,LM)=V(LA)V(LM)(V(LA)+V(LM))/2×100%

This performance index is smaller the better. When all the labelings from the two segmentations are identical, VDPi(LA, LM) is 0.

Results

Segmentation and Model Validation

The eMRF algorithm was implemented and tested on a group of 90 μm isotropic resolution MR images of five C57/BL6 mice. To assure the validity of all results from this experiment, a five-fold cross validation was used in the every step of the experiment. Each mouse was used as a test brain, while the remaining four mice were used as the training set. Hence we tested the algorithm on five distinct training and testing sets. To calculate the observation potential function in Eq. 17, we built five SVM models for the five different training sets using the MRS-SVM procedure (Bae et al., 2008). Building a SVM model for mouse brain segmentation is challenging because the number of the classes is large (>20) and the classes are highly imbalanced. The MRS-SVM procedure enabled us to build SVM model for the brain segmentation efficiently and effectively. Testing each SVM model yielded the posterior probability P(Y|X) in Eq. 17 for each tested mouse. The location potential function was calculated for the five mouse based on Eq. 18. After the observation and location potential functions were calculated for each mouse, the ICM algorithm was used to provide the contextual potential function in Eq. 19 and the posterior probability P(Y|X,L) in Eq. 16. To find the best weight for each potential function we conducted a grid search over the range of each wi= {0.1, 0.2, …, 0.9} (i=1, 2, 3), where w1 + w2 + w3= 1. The best model parameters determined were w1=.1, w2=.6 and w3=0.3 for observation, location and contextual functions, respectively.

We also implemented the atlas-based segmentation for the same mouse brain images to compare its performance with the eMRF method. One mouse brain image was chosen out of the five mouse brains as a reference brain. First, the rest four testing images were realigned to the reference brain using a 12-parameter affine registration which adjusts the global position and size differences between the individual brain images. Next, the reference brain was mapped to the testing brain images using a nonlinear registration which deforms each neuroanatomical structure into a common space. The registrations were done using Image Registration Toolkit (ITK) software and the registration parameters were optimized by the software. Each image of the five mouse brain was used as the reference brain, while the remaining four images were used as the testing brains. Therefore, each brain image was segmented four times with different reference brain each time.

Table 2 presents estimates of the segmentation performance based on two metrics: VOP and VDP between the automated and manual label images. These estimates are based on five mouse brains, and with the use of five-fold cross validation. The results for the eMRF method are compared with three other segmentation methods – the aforementioned atlas-based segmentation method, the MRS-SVM method (Bae et al., 2008) and the MRF method (Ali et al., 2005). Overall eMRF outperforms all the three segmentation methods. Compared to MRS-SVM the average VOP of eMRF is improved by 8.68%, 10.05% compared to the atlas-based segmentation, and 2.79% compared to MRF. Corresponding to an increase in the average VOP, the average VOD of eMRF is decreased. It is improved by 42.04% compared to MRS-SVM, 23.84% compared the atlas-based segmentation, and 12.71% compared to MRF. Figure 2 illustrates the VOP (top) and the VDP (bottom) of the four different segmentation methods.

Figure 2
Comparison of the relative performances of the eMRF, MRS-SVM, atlas-based segmentation and MRF methods based on the voxel overlap percent - VOP (top) and the volume difference percent -VDP (bottom) indices
Table 2
Accuracy of the four segmentation methods (eMRF, MRF-SVM, atlas-based segmentation and MRF) for each of the 21 segmented brain structures is evaluated using the VOP (top row) for each segmentation and VDP (second row) metrics. The percentage improvement ...

The eMRF method outperforms the MRS-SVM method in most structures. The improvement in segmentation performance as assessed by VOP is greatest for white matter structures like the anterior commissure (99.29%), and cerebral peduncle (15.55%) but also for ventricles (35.91%). Exceptions where eMRF underperforms to MRS-SVM are three large structures: CORT, CBLM and OLFB. The MRS-SVM method performs better for large structures because the classifier tends to generate the separation hyperplane between a large class and a small class such that it is biased towards the large class. Therefore, a voxel that is located at the boundary between the large class and the small class is more likely predicted as the large class and this reduces the misclassification of a large class voxel to a small class voxel. However, the eMRF method balances the contributions from the MR intensity information, the location information and the contextual information by assigning potential functions with different weights. In our study, the weight assigned to observation potential function (w1) is smaller than the other two. The averages of the VOP and the VDP of the eMRF method are improved by 8.68%, and 42.04% respectively, compared to the MRS-SVM method, which indicates that the overall segmentation performance has been improved by the use of eMRF.

The eMRF method provides better performances in most structures than the atlas-based segmentation method. Based on the VOP, the eMRF method outperforms the atlas-based segmentation method in all the structures except the two structures, AC and INTP. The performance differences in these two structures are only 0.45% in AC and 5.49% in INTP. The eMRF method has a better VDP performance than the atlas-based segmentation method in 17 structures. There are four small structures (AC, VEN, INTP and CC) that the atlas-based segmentation method has better VDP performance than eMRF. Since each labeling of the voxels in the reference brain is mapped to the testing brain in the atlas-based segmentation process, the number of labeling for a structure in the reference brain is well maintained in the number of labeling for the structure in the testing brain. Therefore, the VDP values of the atlas-based segmentation method in all the structures are quite consistent (11.14–15.23). However, this method has the worst VOP performance among the four different segmentation methods, indicating that the accuracy of atlas-based segmentation is poor. The eMRF method can improve the averages of the VOP and the VDP by 10.05% and 23.84%, respectively, which shows that overall the eMRF method outperforms the atlas-based segmentation method.

Based on the voxel overlap metric, the eMRF method shows better performance in 13 out of 21 structures compared to the MRF method. The eMRF method provides more than 5% performance improvement in 7 structures (CPED, PAG, MED, PON, SNR, OPT and TRI), the largest improvements are seen for optic tract (29.39%) and trigeminal tract (15.63%). Based on the volume difference metric, the eMRF method shows better performance in 16 structures out of 21 structures. VDP improvement ranged from 61.68% for optic tract and to 3.61% for globus pallidus. In four small structures (PAG, INFC, AC and INTP), which take less than 1% volume compared to the whole brain volume, the MRF method shows better performance. This is due to the fact that in the MRF method labeling a voxel depends on the labelings of the neighborhood of the voxel, thus the identification of small structures is enhanced. Nevertheless, the averages of the VOP and the VDP of the eMRF method are improved by 2.79%, and 12.71% respectively, compared to the MRF method. The eMRF method can make a balance between the overfitting of the SVM method to the larger classes and the overfitting of the MRF method to the smaller classes so that the overall segmentation performance is improved.

A visual comparison of the three segmentation methods, manual labeling (the gold standard), eMRF and MRS-SVM, is shown for two specific coronal levels (one at the level close to the hippocampal commissure, and one at the level of the pons) in Figure 3. On the left column of Figure 3 are displayed the manual segmentations (used as the gold standard), on the middle column are the automated segmented images by the eMRF method and the right column, those produced by the MRS-SVM method. The eMRF method seems to deviate less than the MRS-SVM method from the manual labeling, and to respect better the topology of the structures. The testing time of the MRS-SVM algorithm for a mouse brain dataset (472,100 voxels) was 289.4 minutes with a 3.4-GHz PC. The testing experiment was run using MATLAB, and executed the eMRF segmentation of one mouse brain in 75 minutes.

Figure 3
Coronal slices through the labeled brain at the level of anterior hippocampus and third ventricle (upper row), and pons and substantia nigra (lower row) show in a qualitative manner the relative superiority of eMRF compared to MRS-SVM. Note that eMRF ...

Discussion and Conclusion

In this paper, we developed an eMRF method, which first adds a potential function based on location information, then integrates the SVM and MRF for a more accurate and robust segmentation of mouse brain MR images by taking the advantages of the two methods. MRF has been used in the human brain image segmentation task because it utilizes the contextual information of a voxel, which describes the correlation of the voxel and its neighbors for segmentation (Fischl et al., 2002). A similar algorithm has been successfully implemented for the mouse brain (Ali et al., 2005). Other than contextual information and MR intensity information, Bae et al. (2008) recently studied the use of SVM based on MR intensity information and location information for mouse brain segmentation and a novel automated segmentation method, termed MRS-SVM, was proposed. Experimental results indicate that the MRF method outperforms the MRS-SVM method for smaller structures while the MRS-SVM method outperforms the MRF method for larger structures. The complementary nature of the two methods directs the development of eMRF. Specifically, other than using Gaussian probability distribution to model the MR intensity signal, the eMRF method employs the posterior probability distribution obtained from the SVM to generate classification based on the MR intensity information. Secondly, eMRF introduces a new potential function based on the location information. Instead of considering the contributions from the three potential functions – observation, location and contextual functions – equally, eMRF further applies ICM to optimally determine the contribution weights for each function.

To validate the proposed method, we conducted a comparison of the four different algorithms: MRF (Ali et al., 2005), MRS-SVM (Bae et al., 2008), atlas-based segmentation and eMRF for the automate segmentation of mouse brain into 21 structures, using the same dataset. Our test results show that the overall performance of the eMRF method is better than all the three segmentation methods in terms of both the average VOP and average VDP, even though the MRS-SVM method is slightly better on a small number of large structures, and the atlas-based segmentation and MRF methods are slightly better on a few small structures.

In the future, we will extend and adapt the eMRF method to human brain segmentation and compare the results with the existing study in literature (Powell et al., 2008). Further developments of this work would also include obtaining higher resolution images and more labels of neuroanatomical structures, which could provide more information to be incorporated in atlases of the normal mouse brain. Moreover, we foresee that the automated segmentation method described in the paper will accelerate the study of brain images of large quantity, thus help developing small animal disease models for many neurodegenerating disorders, such as Alzheimer disease and Parkinson disease.

Acknowledgments

The authors would like to thank the Duke Center for In Vivo Microscopy, an NCRR/NCI National Biomedical Technology Research Resource (P41 RR005959/U24 CA092656), for providing the images. They would also like to thank the Mouse Bioinformatics Research Network (MBIRN) (U24 RR021760) for providing support for this imaging study.

Footnotes

1Using the two indices VOP (the larger the better), VDP (the smaller the better) introduced in the next section, the experiment on two MR protocols with coordinates yields 72.42% VOP and 19.21% VDP in average. The experiment on all five MR protocols with coordinates generates 59.38% VOP and 24.16% VDP.

References

  • Ali AA, Dale AM, Badea A, Johnson GA. Automated segmentation of neuroanatomical structures in multispectral MR microscopy of the mouse brain. NeuroImage. 2005;27 (2):425–435. [PubMed]
  • Amato U, Larobina M, Antoniadis A, Alfano B. Segmentation of magnetic resonance images through discriminant analysis. Journal of Neuroscience Methods. 2003;131:65–74. [PubMed]
  • Andersen AH, Zhang Z, Avison MJ, Gash DM. Automated segmentation of multispectral brain MR images. Journal of Neuroscience Methods. 2002;122:13–23. [PubMed]
  • Ashburner J, Friston K. Multimodal image coregistration and partitioning – a unified framework. NeuroImage. 1997;6:209–217. [PubMed]
  • Bae MH, Wu T, Pan R. Mix-Ratio Sampling: Classifying Multiclass Imbalanced Mouse Brain Images Using Support Vector Machine. 2008. Technical Report available at http://swag.eas.asu.edu/vcie/
  • Besag J. On the statistical analysis of dirty pictures. Journal of Royal Statistical Society, Series B. 1986;48(3):259–302.
  • Bock NA, Kovacevic N, Lipina TV, Roder JC, Ackerman SL, Henkelman RM. In vivo magnetic resonance imaging and semiautomated image analysis extend the brain phenotype for cdf/cdf mice. The Journal of Neuroscience. 2006;26 (17):4455–4459. [PubMed]
  • Duta N, Sonka M. Segmentation and interpretation of MR brain images: an improved active shape model. IEEE Transactions on Medical Imaging. 1998;16:1049–1062. [PubMed]
  • Fischl B, Salat DH, Busa E, Albert M, Dieterich M, Haselgrove C, Van der Kouwe A, Killiany R, Kennedy D, Klaveness S, et al. Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain. Neuron. 2002;33:341–355. [PubMed]
  • Geman S, Geman D. Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1984;6:721–741. [PubMed]
  • Gonzalez Ballester MA, Zisserman A, Brady M. Segmentation and measurement of brain structures in MRI including confidence bounds. Medical Image Analysis. 2000;4:189–200. [PubMed]
  • Heckemann RA, Hajnal JV, Aljabar P, Rueckert D, Hammers A. Automatic anatomical brain MRI segmentation combining label propagation and decision fusion. NeuroImage. 2006;33:115–126. [PubMed]
  • Hsu CW, Lin CJ. A comparison of methods for multi-class support vector machines. IEEE Transactions on Neural Networks. 2002;13(2):415–425. [PubMed]
  • Kovacevic N, Lobaugh NJ, Bronskill MJ, Levine B, Feinstein A, Black SE. A robust method for extraction and automatic segmentation of brain images. NeuroImage. 2002;17:1087–1100. [PubMed]
  • Kovacevic N, Henderson JT, Chan E, Lifshitz N, Bishop J, Evans AC, Henkelman RM, Chen XJ. A three-dimensional MRI atlas of the mouse brain with estimates of the average and variability. Cerebral Cortex. 2005;15 (5):639–645. [PubMed]
  • Lee CH, Schmidt M, Murtha A, Bistritz A, Sander J, Greiner R. Segmenting brain tumors with conditional random fields and support vector machines. Lecture Notes in Computer Science. 2005;3765:469–478.
  • Li SZ. Markov Random Field Modeling in Image Analysis. Springer-Verlag; Tokyo: 2001.
  • Ma Y, Hof PR, Grant SC, Blackband SJ, Bennett R, Slatest L, Mcguigan MD, Benveniste H. A three-dimensional digital atlas database of the adult C57BL/6J mouse brain by magnetic resonance microscopy. Neuroscience. 2005;135 (4):1203–1215. [PubMed]
  • Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, Mazoyer B, Joliot M. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage. 2002;15:273–289. [PubMed]
  • Platt J. Advances in Large Margin Classifiers. MIT Press; Cambridge, MA: 2000. Probabilistic outputs for support vector machines and comparison to regularized likelihood methods.
  • Powell S, Magnotta VA, Johnson H, Jammalamadaka VK, Pierson R, Andreasen NC. Registration and machine learning-based automated segmentation of subcortical and cerebellar brain structures. NeuroImage. 2008;39:238–247. [PMC free article] [PubMed]
  • Reddick WE, Glass JO, Cook EN, Elkin TD, Deaton RJ. Automated segmentation and classification of multispectral magnetic resonance images of brain using artificial neural networks. IEEE Transactions on Medical Imaging. 1997;16:911–918. [PubMed]
  • Suri J. Two-dimensional fast Magnetic Resonance brain segmentation. IEEE Engineering in Medicine and Biology. 2001 July/August;:84–95. [PubMed]
  • Vapnik VN. The Nature of Statistical Learning Theory. Springer; Verlag: 1995.
  • Zavaljevski A, Dhawan AP, Gaskil M, Ball W, Johnson JD. Multi-level adaptive segmentation of multi-parameter MR brain images. Computerized Medical Imaging and Graphics. 2000;24:87–98. [PubMed]
  • Zhang Y, Smith S, Brady M. Segmentation of brain MR images through a hidden Markov random field model and the expectation–maximization algorithm. IEEE Transactions on Medical Imaging. 2001;20:45–57. [PubMed]