Partial volume effect modeling for segmentation and tissue classification of brain magnetic resonance images: A review

Abstract

Quantitative analysis of magnetic resonance (MR) brain images are facilitated by the development of automated segmentation algorithms. A single image voxel may contain of several types of tissues due to the finite spatial resolution of the imaging device. This phenomenon, termed partial volume effect (PVE), complicates the segmentation process, and, due to the complexity of human brain anatomy, the PVE is an important factor for accurate brain structure quantification. Partial volume estimation refers to a generalized segmentation task where the amount of each tissue type within each voxel is solved. This review aims to provide a systematic, tutorial-like overview and categorization of methods for partial volume estimation in brain MRI. The review concentrates on the statistically based approaches for partial volume estimation and also explains differences to other, similar image segmentation approaches.

Key Words: Magnetic resonance imaging, Segmentation, Tissue classification, White matter, Gray matter, Image processing, Brain imaging, Image analysis

Core tip: Each voxel in a brain magnetic resonance imaging (MRI) may contain multiple types of tissue. Partial volume estimation refers to a generalized image segmentation task where the amount of each tissue type within each image voxel of brain MRI is solved. This is important for volume quantification and cortical thickness analysis due to the geometrical complexity of human brain structure. This review aims to provide a systematic, tutorial-like overview of methods for partial volume estimation in brain MRI.

INTRODUCTION

Quantitative analysis of magnetic resonance (MR) brain images to gain knowledge about human brain structure is increasingly important. For example, various neuropsychiatric and neurodegenerative diseases, such as schizophrenia[1] and Alzheimer’s disease[2], alter the brain structure. By analyzing these alterations, a better understanding of the underlying disease mechanisms could be gained and diseases could potentially be diagnosed more rapidly and accurately[3]. This is important since brain diseases represent a major source of the overall disease burden[4] and are often associated with heavy impact to informal caregivers.

The typical quantitative analyses to detect and quantify differences in brain structure between two or more subject groups include voxel based morphometry[5] and cortical thickness analysis[6]. These analyses are facilitated by the development of automated MR image (MRI) segmentation algorithms, which are standard tools in modern neuroscience. The image processing chain leading to MRI segmentation and, finally, to statistical analyses, comprises of a long pipeline of different operations including skull stripping, intensity non-uniformity correction, tissue classification, registration to the stereotactic space and cortical surfaces extraction. The point of interest in this review is the tissue classification. This refers to assigning a tissue type label to each voxel of a brain image. Typically, the three main tissue types, white matter (WM), gray matter (GM), and cerebro-spinal fluid (CSF), are considered.

A single voxel may contain of several types of tissues due to the finite spatial resolution of the imaging device. This phenomenon, termed partial volume effect (PVE), complicates the segmentation process, and, due to the complexity of human brain anatomy, the PVE is an important factor when accurate brain structure quantification is needed; Figure 1 for a schematic explanation of the PVE in the context of brain MRI. González Ballester et al[7,8] reported that ignoring the PVE can lead to volume measurement errors in the range of 20%-60%. Most widely used MRI segmentation algorithms account for PVE, for example by incorporating extra tissue classes[9-11]. Ruan et al[12] demonstrated that the intensity distributions of the partial volume voxels can be approximated using Gaussian distributions and an early work attributed the non-normality of the intensity distributions of the tissue classes to partial volume artefact[13]. However, some algorithms take a step further and try to solve an extended version of the tissue classification problem, where the amount of each tissue type within each voxel is solved. For example, hard or crisp tissue classification provides information whether a particular voxel is WM, GM, or CSF. In the extended problem, one wants to know that a voxel contains 20% GM, 80% WM and 0% of CSF and we say that the partial volume coefficients (PVCs) are 20% for GM, 80% for WM and 0% for CSF. The extended problem has various names. It has been referred to as fuzzy segmentation, partial volume segmentation, partial volume estimation, and tissue fraction estimation. It will be referred to as partial volume estimation in the remainder of this paper. In order for the partial volume estimation problem to be solvable, the intensity of a partial volume voxel has to be expressed with a model that depends on the parameters of image intensity distributions of pure tissue classes. Figure 2 exemplifies partial volume estimation as compared to hard tissue classification and also points out a specific problem of hard tissue classification particularly important to cortical thickness computations. Namely, insufficient image resolution may lead to hard tissue classification miss sulcal CSF and this may subsequently lead to incorrect cortical thickness computation if hard tissue classification is used as a preprocessing operation to the cortical thickness computation.

Figure 1 A schematic explanation of the partial volume effect in the context of brain magnetic resonance imaging. Voxels composed of purely gray matter (GM) are colored in black color while voxels composed of cerebro-spinal fluid (CSF) or white matter (WM) are in white color. These are termed pure tissue voxels or pure voxels. Voxels composed of multiple tissue types, termed mixed voxels, are colored in gray. In the figure, these can be either voxels containing both CSF and GM tissue types or voxels containing both WM and GM tissue types. The actual anatomical boundaries between tissue types are shown in blue and red color is used to indicate voxel boundaries.

Figure 2 Example of partial volume estimation. Top row, from left: A coronal section of T1 weighted MR image; A skull stripped version of the coronal section; A manual labeling into gray matter (GM) (gray color), white matter (WM) (white color), and cerebro-spinal fluid (CSF) (dark gray color). Bottom row: Estimates of partial volume coefficients (PVCs) for CSF, GM, and WM. The color bar refers to the PVC estimates in the bottom row. The image is obtained from the IBSR2 dataset provided by the Center for Morphometric Analysis at Massachusetts General Hospital and PVCs were computed as described in the ref. [28]. Note how the manual hard labeling completely misses the CSF in the interhemispheric fissure as well as in the superior frontal sulcus pointed by red arrows. Instead PVC estimates of CSF in the bottom row capture well the sulcal CSF.

This review aims to provide a systematic, tutorial-like overview and categorization for different approaches for partial volume estimation in brain MRI. In addition of the author’s knowledge about existing literature, the articles to be included in this review were searched on Pubmed: Search term: [(magnetic resonance [Title/Abstract] OR MRI [Title/Abstract]) AND brain [Title/Abstract] AND partial volume [Title/Abstract] AND (segmentation [Title/Abstract] OR tissue classification [Title/Abstract] OR partial volume coefficient estimation [Title/Abstract])] NOT (PET [Title/Abstract] OR emission tomography [Title/Abstract]). The search yielded 80 articles, majority of which were found relevant to this review.

IMAGE PRE-PROCESSING

The algorithms introduced in next sections require various image pre-processing steps to be performed before the partial volume estimation can take place. The pre-processing pipeline can include intensity non-uniformity correction, brain extraction (or skull stripping) and registration to a sterotactic space.

Intensity non-uniformity correction is required because MR images are known to contain low frequency spatial intensity variations often referred to as radio frequency inhomogeneity or shading artifact[14]. All segmentation algorithms in brain MRI must account for this artifact to produce accurate segmentations. There are several ways to correct for the shading artifact[14]. This can be assumed to be an image pre-processing step or to be performed jointly with the PV estimation, interleaving PV estimation (segmentation) and non-uniformity correction steps. In what follows, we will assume that the images have been corrected for this artifact.

Although we are interested in segmentation of the brain tissues, brain MR images contain signal from other, extracerebral tissue types, such as skull or scalp. Because these extracerebral tissue types are often irrelevant for brain image quantification, it is useful to mask out the voxels outside the brain out before the PV estimation. This is termed skull stripping or brain extraction and the reference[15] provides a comparison of skull stripping algorithms.

The registration to stereotactic space is usually carried out to be able to utilize information of the tissue type probability maps, which, for each voxel, give a prior probability that the voxel is of certain issue type[16]. It should be noted that this is not as useful for partial volume estimation as it can be for hard segmentation, because tissue probability maps provide no information on tissue fractions[17]. Moreover, if the registered images are resampled to the stereotactic space, this amplifies the partial volume effect and may not be a recommended action.

MIXEL MODEL

Definition and approximations

The most commonly used model of PVE in brain MRI is the mixel model[18]. The mixel model assumes that each intensity value in the image is a realization of a weighted sum of random variables (RVs), each of which characterizes a pure tissue type. The original formulation[18] requires images to be multispectral, i.e., that image data from multiple pulse sequences are available (for example, T1, T2,and proton density weighted images). However, there are approaches to overcome this problem by utilizing clever approximations as we shall see in Section 3.2.

We now proceed to a more formal description of the mixel model. For this, we need to establish some notation. The observed image is X = {x_i: i = 1,…,N}, with the voxel intensity x_i∈ R^K, and K the number of data channels in the multispectral case. For example, if we have T1-, T2-, and proton density-weighted images, then K = 3. N denotes the number of brain voxels in the image and i is the voxel index. The voxel index has three components that correspond to the position of the voxel in the left-right, anterior-posterior, and inferior-superior axes. There are M tissue types in the image. Typically, M is equal to 3, and the tissue types are WM, GM, and CSF. The mixel model is statistically based. Thus, a voxel intensity x_i is considered to be a realization of random variable x_i. (We use bold-face symbols to refer to random variables and the corresponding normal-face symbols denote their realizations.) Similarly, each tissue type j is described by a random variable l_j , which is assumed to be distributed according to the multivariate normal distribution with the mean μ_j and covariance Σ_j. Random variable x_i is written as a weighted sum

Math 1

Math 1 Math(A1).

where n represents measurement noise, typically assumed to be Gaussian (with a covariance matrix Σ^*) and partial volume coefficients (PVCs) w_ij∈ [0, 1] for all i,j and

Math 2

Math 2 Math(A1).

for all i. The PVCs model the fraction of each tissue type in the voxel, for example, if w_iGM has a value of 0.8 then the voxel contains 80% of the GM tissue type. This is similar to the fuzzy classification/segmentation problem, but in the mixel model the coefficients w_ij specifically model the fraction of tissue type j present in the voxel i. We will return to connections of the mixel model and the Fuzzy C-means algorithm in Section 5.

In practice, the mixel model has to be simplified because it is impossible to distinguish between measurement noise and variability within tissue types. Various simplifications have been studied by Santago et al[19,20]. They identified two possible types of simplification, namely, the sampling noise model and material dependent noise model as depicted in Figure 3. The sampling noise model assumes that all the randomness in the model is due to measurement noise. This leads to a model, where the tissue types are described by mean intensities of tissue types:

Figure 3 Sampling and material dependent noise models. Sampling noise model assumes that each tissue type is represented by a single average value and Gaussian-distributed noise is then added. Material dependent noise model assumes that the tissue types are represented by random variables. CSF: Cerebro-spinal fluid; GM: Gray matter; WM: White matter.

Math 3

Math 3 Math(A1).

The material dependent noise model is obtained by embedding the measurement noise into material noise components, i.e., n is dropped from Eq. (1)

Math 4

Math 4 Math(A1).

This model is more complex than the sampling noise model, but it is probably more realistic.

Solving the mixel model

Direct solutionviapenalized least squares: Assuming the sampling noise model, the PVCs can be solved directly from Eq. (2) if enough data channels are available[18]. Denoting a matrix of all PVCs by w, the least squares criterion to minimize for solving Eq. (2) is written as

Math 5

Math 5 Math(A1).

with constraints that

Math 6

Math 6 Math(A1).

and 0 ≤w_ij≤ 1. Note that this equation can be solved individually for each voxel. In the case of single image channel and two tissue types, the solution is particularly simple:

Math 7

Math 7 Math(A1).

and the function r limits the solution to the interval from 0 to 1, i.e., r (y) = 0 when y < 0, r(y) = y when 0 ≤y≤ 1, and r (y) = 1 when y > 1. This solution is also the maximum likelihood solution and it accounts to a simple scaling of the image intensities to the interval from 0 to 1. For this reason, the solution is also very noisy and Choi et al[18] suggested to regularize it with a Markov Random Field (MRF) prior (see also Li et al[21]). The idea is that PVCs of neighboring voxels should have similar values. This leads to a modified criterion to minimize, with the same constraints as above,

Math 8

Math 8 Math(A1).

where the term P(w) penalizes differences between w_i = [w_i1,…,w_iM] and w_k = [w_k1,…,w_kM] if the voxels i and k are neighbours. Unfortunately, this objective cannot be anymore minimized separately for each voxel, but all the voxels must be taken into account. Besag et al[22] used Iterative Conditional Modes algorithm to minimize the penalized least squares criterion in Eq.(6).

Two step algorithms: The simple two-class, one-channel solution above motivates a set of techniques allowing the standard PVC estimation for three tissue types even if just data from just a single image (usually T1-weighted) is available. The idea is that since the combination of more than two tissue types in a voxel is very rare, we can estimate which two tissue types are present in a voxel before the PVC estimation; Alike idea was already mentioned for multichannel data in[18,23]. The steps of the two step algorithm can be given as follows, and they are schematically represented in Figure 4: (1) Partial volume classification: Estimate which is most likely tissue type configuration containing at most two tissue types in each voxel; and (2) PVC estimation: Solve the partial volume estimation problem limited to tissue types found in Step 1 for all the voxels.

Figure 4 Direct vs two step procedure for partial volume coefficient estimation. CSF: Cerebro-spinal fluid; GM: Gray matter; WM: White matter; PVC: Partial volume coefficient.

There are at least three different approaches to solve the task in the step 1. In the simplest approach, used for example in the reference[24], the tissue classes are ordered based on their mean values so that μ₁< μ₂< ..... < μ_M. Then, if the intensity value x_i lies in the interval [μ_k,μ_k+1], it is assumed that the voxel i is a mixture of tissue types k and k + 1. This simple model does not account for the noise in the images and is not applicable for multichannel data because it assumes that the mean intensity values of tissue types can be ordered. The second approach is to detect most likely pure tissue types within the voxel based on the Bayes classifier[18,23]. This is done computing the two most probable tissue types within a voxel. However, this approach, as the first one, ignores the possibility that voxels may be composed of a single tissue type. The third and preferred approach, which is we term as probabilistic partial volume classification, fixes the just mentioned problem. The probabilistic partial volume classification approach is to compute the probability of each possible tissue type mixture appearing in the voxel[19,20,25-27]. For example, if the tissue types of interest are WM, GM, and CSF, the following 6 probabilities are computed: (1) Voxel is solely CSF; (2) Voxel is solely GM; (3) Voxel is solely WM; (4) voxel is a mixture of background and CSF; (5) voxel is a mixture of CSF and GM; and (6) voxel is a mixture of GM and WM. (Some tissue type combinations are not considered due to their rarity in the brain.) The technical problem in the probabilistic partial volume classification approach is the construction of the probability models for mixed tissue classes; the class conditional densities for pure tissue classes are modelled by the normal density. The probability densities for the mixed tissue types can be constructed based on a marginalization technique developed originally in references[19,20] and further applied in references[25-28]. The idea is to integrate out the variable w_i1describing the percentage of tissue type 1 in a voxel by numerical integration. Note that with current computers the numerical integration does not present computational problem and can be solved very fast[28]. Advantages of this more complicated probabilistic approach over the two simple approaches include possibility to include spatial regularization in the form of MRFs to the step 1[25,26] and the applicability to multispectral images[26]. Additionally, it is often expected that the number of the pure tissue voxels should be greater than the number of mixed tissue voxels. The probabilistic partial volume classification includes automatic and elegant control for this issue that has been solved elsewhere by using Bayesian methods at the expense of introducing extra user-defined parameters[24,29].

Once the tissue types that are probable to appear in a voxel are determined, then the PVCs can be estimated using Eq. (5) if the sampling noise model is assumed. Note that if voxel i is determined to be a voxel of pure tissue type k, then w_ik = 1 and w_ij = 0 for other tissue types j ≠ k. One can also adopt the material dependent noise model leading to a maximum likelihood criterion. If i is a mixed voxel of tissue types j and k, the maximum-likelihood solution is

Math 9

Math 9 Math(A1).

where μ(w) = wμ_j + (1-w)μ_k; Σ(w) = w²Σ_j + (1-w)²Σ_k or Σ(w) = wΣ_j + (1-w)Σ_k. Furthermore, w_ik^* = 1 - w_ij^* and all the other PVCs are zero. The correct model for Σ(w) has caused some controversy (see the references[30,31] for details). The difference in the two models is that the first one (Σ(w) = w²Σ_j + (1 - w)²Σ_k) results in a more regularized solution of Eq. (7) while the second one (Σ(w) = wΣ_j + (1 - w)Σ_k ) is conceptually more pleasing. The maximum-likelihood PVC-estimate in Eq. (7) is solved by a simple grid search. Extensions to the maximum likelihood principle of Eq. (7) include Bayesian methods[24].

As mentioned above, the two-step algorithms can use the MRF prior to regularize the partial volume classification and this has been demonstrated to lead to more accurate partial volume estimates when the images are noisy[25]. The use of the MRF requires the user to set a proper weighting parameter for the prior which may be considered as a disadvantage[8]. However, often quoted disadvantage of the added computational cost (e.g., the reference[8]) of the MRF, can be overcome by new rapid algorithms capable of performing MRF based segmentation of the typical 3-D MR images within few seconds[28]. While the two-step algorithms often use spatial MRF prior during the partial volume classification step, they typically do not utilize spatial information during the second, PVC estimation, step. Manjón et al[27] introduced an MRF for modelling of the spatial information during the PVC estimation step and compared it to the usage of prefiltering the images with a non-local means filter. The results suggested that using spatial information improved the PVC estimates and non-local means filtering performed better than the MRF-based approach.

Discretization approaches: An alternative to try to find real-valued PVC estimates is to discretize the PVC estimation problem[32-34]. This means that instead of letting each PVC w_ij lie freely in the interval from zero to one, the discretization-based methods restrict the PVCs to have only a discrete set of values. For example, w_ij can be 0, 0.1, 0.2,…,1.0. The discretization-based methods then try to solve maximally probable PVCs from this discretized set resorting MRF approaches to model spatial interaction between adjacent voxels[32-34]. While the restriction to a discrete set of PVC values is perfectly reasonable given the noisiness of the images, the discretization approaches are usually very time consuming, especially when compared to fast two step approaches[25,28].

Parameter estimation

The necessary model parameters μ_j,j = 1,…,M and Σ^* or Σ_j,j = 1,…,M must be estimated before or during the solution of the mixel model. Correct estimation of these parameters is essential for partial volume estimation[35]. Tohka et al[26] identified three potential approaches to the parameter estimation problem: (1) histogram analysis; (2) simultaneous parameter, and partial volume estimation by expectation maximization (EM)-like algorithms; and (3) the estimation based on a hard segmentation of the image.

The conceptually simplest alternative is to fit a parametric model (a mixture model of pure and mixed tissue intensity densities) to an image histogram. The objective function can be based on the maximum likelihood or least squares criterion. The disadvantage of parametric model fitting is that the formulated minimization problem is complex and non-convex rendering the standard optimization algorithms useless. Various global optimization algorithms, including genetic algorithms and tree annealing, have been used for the task[19,36]. The EM-like algorithms start from an initial rough parameter estimates and refine the estimates jointly with the partial volume estimation[32,34] or classification[37] through alternating expectation and maximization steps. This can guarantee accurate parameter estimates, but the estimates depend strongly on the initial guess and the convergence of the process can be slow. The third alternative is to generate an initial rough segmentation of the image, and thereafter use outlier detection techniques based on the mathematical morphology, robust point estimates, or image gradient values to prune the set of voxels belonging to a certain tissue class[25,26,35,38,39]. Comparisons of these three techniques have been reported in the references[26,35]. The main result of these comparisons has been that the parameter estimation based on the hard segmentation of the image is fast and usually, but not always, works as well or better than the other two approaches.

RELATED METHODS

Fuzzy C-means

The standard Fuzzy C-means (FCM) algorithm optimizes a cost function

Math 10

Math 10 Math(A1).

where u_ij are the fuzzy membership values μ_k are the class centroids, and q is the fuzzification parameter. This objective function and its modifications have been widely and successfully used for brain MRI tissue classification[40-43]. As shown in the reference[29], if q = 3, M = 2, and K = 1, optimizing the objective J_FCM for fixed centroids leads to the identical PVCs as PVCs derived based on Eq. (5). However, with more than two tissue types or multispectral data, fuzzy segmentations by FCM and mixel model are different.

Bayesian tissue classifiers

Often the tissue classification is casted as the Bayesian decision problem[9,16,17,44,45]. In that, one tries to estimate the posterior probability map that the tissue type is c given the image intensities. Often approaches use prior information from tissue probability maps[9,16] or MRFs[44,45] or both[17]. It should be noted that the tissue type probabilities are different from the partial volume coefficients. The exact difference of the segmentation results depends on the probability model selected, but usually these Bayesian tissue classifiers produce more crisp tissue type maps than the partial volume estimation algorithms. This issue and its ramifications are considered in a more detail by Manjón et al[27].

APPLICATIONS OF PARTIAL VOLUME ESTIMATION

Voxel based morphometry

Voxel-based morphometry (VBM) involves a voxel-wise comparison of the local concentration of gray matter between two groups of subjects. The procedure consists of segmenting the gray matter from the MR images and spatially normalizing these gray matter images from all the subjects in the study into the same stereotactic space[5]. These gray matter images can either represent GM tissue probabilities, for example, as in the reference[46] or GM tissue fractions resulting partial volume estimation, for example, as in the reference[47]. While it seems clear that the PVCs are better representations of gray matter density than gray matter probabilities, it is not clear whether this particular modelling choice has a major effect on the accuracy of the results. To author’s knowledge, gray matter probability and gray matter PV-coefficient based VBM methods have not been directly compared. Tardif et al[48] examined two pipelines resulting in GM probability based VBM and PVC based VBM but the main focus of the work was on a comparison of 1.5T and 3T imaging protocols. The VBM8 software package (http://dbm.neuro.uni-_jena.de/vbm/) offers possibility to VBM using PVCs[49].

Cortical thickness

Cortical thickness is a quantitative measure describing the combined thickness of the layers of the cerebral cortex that can be measured using MRI either using mesh based[6,50,51] or voxel based techniques[52]. The thickness of the cortex, and its local variations, are of great interest in both normal development as well as a wide variety of neurodegenerative and psychiatric disorders[6]. Cortex is a highly folded structure with an approximate average thickness of 2.5 mm[53] and hence it is not difficult to appreciate that the partial volume effect has been an important consideration when measuring cortical thickness. Both surface mesh based[54] and voxel based[55-57] cortical thickness measures can be shown to be improved if the partial volume effect is taken into account. Especially, as demonstrated in Figure 2 and discussed further in the references[26,54], hard tissue classifications may miss some of the sulcal CSF because of an insufficient image resolution. This causes incorrect reconstruction of the GM/CSF boundary, which, in turn, leads to errors in the cortical thickness computation.

Other applications

Other applications of segmentation with the PVE modeling identified during the literature review were segmentation of the brain images of the neonates[58-61], hemisphere segmentation and related shape analysis[62,63], EEG source localization[64], and lesion load computations based on MRI[65-68]. Especially, in the case of the Multiple Sclerosis (MS) lesion volumetry, the correction for the partial volume effects has a large positive effect on the reproducibility and accuracy of the analysis[69]. In particular, it was found to be important in avoiding of misclassification of some non-lesion voxels (between CSF and brain tissue) into lesion voxels[69].

CONCLUSION

An interesting recent development in MRI segmentation and partial volume estimation is the use of quantitative tissue type maps for the purpose[70-72]. For example, Ahlgren et al[70] utilized the signal of a spoiled gradient-recalled echo (SPGR) sequence acquired with multiple flip angles to map T1, and subsequently to fit of a multi-compartment model yielding parametric maps of partial volume estimates of the different compartments. West et al[71] used quantitative MRI values of the longitudinal relaxation rate, the transverse relaxation rate and the proton density to define tissues (WM,GM,CSF) and constructed a lookup table for partial volume estimation. These quantitative approaches show good potential to improve the partial volume estimation accuracy. Another recent development is the use of high-field MRI to map smaller and smaller brain structures[73], such cortical layers or hippocampal subfields[74]. These efforts will benefit from automated segmentation. Despite of improved image resolution provided by higher field strengths the problems related to partial volume effect will remain as the structures of interest will become smaller at the same time. For example, while the improved image resolution will diminish (but not completely erase) the challenges related to partial volume effect in the cortical thickness computation, it will also possibly allow studies concerning individual cortical layers requiring a higher image resolution, where partial volume effect is again an important consideration.