Elsevier

NeuroImage

Volume 23, Issue 1, September 2004, Pages 84-97
NeuroImage

Fast and robust parameter estimation for statistical partial volume models in brain MRI

https://doi.org/10.1016/j.neuroimage.2004.05.007Get rights and content

Abstract

Due to the finite spatial resolution of imaging devices, a single voxel in a medical image may be composed of mixture of tissue types, an effect known as partial volume effect (PVE). Partial volume estimation, that is, the estimation of the amount of each tissue type within each voxel, has received considerable interest in recent years. Much of this work has been focused on the mixel model, a statistical model of PVE. We propose a novel trimmed minimum covariance determinant (TMCD) method for the estimation of the parameters of the mixel PVE model. In this method, each voxel is first labeled according to the most dominant tissue type. Voxels that are prone to PVE are removed from this labeled set, following which robust location estimators with high breakdown points are used to estimate the mean and the covariance of each tissue class. Comparisons between different methods for parameter estimation based on classified images as well as expectation–maximization-like (EM-like) procedure for simultaneous parameter and partial volume estimation are reported. The robust estimators based on a pruned classification as presented here are shown to perform well even if the initial classification is of poor quality. The results obtained are comparable to those obtained using the EM-like procedure, but require considerably less computation time. Segmentation results of real data based on partial volume estimation are also reported. In addition to considering the parameter estimation problem, we discuss differences between different approximations to the complete mixel model. In summary, the proposed TMCD method allows for the accurate, robust, and efficient estimation of partial volume model parameters, which is crucial to a variety of brain MRI data analysis procedures such as the accurate estimation of tissue volumes and the accurate delineation of the cortical surface.

Introduction

The quantitative analysis of magnetic resonance (MR) images in the study of human brain anatomy is becoming more and more important. For example, a range of brain disorders as well as brain development and healthy aging can cause structural changes in the brain. These changes can be quantified by measuring volumes or other properties of anatomical structures of interest providing information, for example, on disease severity. Before measurements can be performed, the structures of interest must be extracted from the image data. This often includes the labeling of voxels according to their tissue type. This labeling, or classification, can be performed based on a single MR image or based on a multispectral image constructed by combining series images of the same subject acquired with different pulse sequence parameters. Typically, the tissue types of interest are white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF), but also more anatomical labels can be assigned to image voxels Collins et al., 1999, Fischl et al., 2002. However, voxel classification with three basic tissue types has itself rather direct applications such as quantifying disease burden in multiple sclerosis by estimation of the amount of brain atrophy (Collins et al., 2001). Moreover, many procedures aiming at the extraction of particular brain structures, such as cerebral cortex, can gain from the initial tissue classification MacDonald et al., 2000, Xu et al., 1999, Zeng et al., 1999 as can partial volume correction in positron mission tomography (Rousset et al., 1998).

Because of the finite resolution of the imaging devices, a single voxel may contain several tissue types. This is known as partial volume effect (PVE). Due to PVE, the classification of a voxel reflecting the dominant tissue type (WM, GM, or CSF), does not reveal all possible information about the tissue content of that voxel. This can be problematic in small structures or highly convoluted areas of the brain. For example, algorithms aimed at extraction of the cortical surface often omit entire sulci due to the PVE between the thin ribbon of sulcal CSF and the surrounding gray matter. These problems are especially serious when pediatric subjects are considered as illustrated in Fig. 1. Estimation of the amount of each tissue type present in each voxel, that is, partial volume estimation, provides an interesting possibility to improve the accuracy of cortical surface extraction (cf. Fig. 1). Other applications that gain from modeling of PVE have been as well considered within brain MRI. Santago and Gage (1993) apply information about partial volume voxels to improve tissue quantification. González Ballester et al. (2002) study the asymmetry of temporal horns taking PVE into account, and in an earlier work (González Ballester et al., 2000), they suggest that PVE and discrete sampling at boundary locations can lead to volume measurement errors in the range 20–60%.

Partial volume effect and PV estimation have been addressed in various ways in MR imaging literature. For example, Pham and Prince (1999) have proposed a fuzzy C-means algorithm. They have also studied the relationship between the fuzzy C-means objective function and statistical models of PVE in the simplified case (only two type is of tissues and single-spectral data), showing that in this case, these two approaches can be considered equivalent (Pham and Prince, 1998). However, these considerations do not extend to more realistic situations where one would have more than two tissue types and where the data would be multispectral.

Wang et al. (2001) propose to use a Bayesian classifier with a variable number of tissue classes, including classes of mixed tissue types. However, as the authors model the image histogram with a finite mixture of normal distributions and parameters (means and variances) for mixed tissue classes, which are independent of the parameters of the related pure tissue classes, there is no explicit model of the PVE. While the procedure may be reasonable for detection of voxels containing PVE, estimation of the mixing proportions can be challenging.

The most commonly used, statistically based model of PVE is the mixel model proposed by Choi et al. (1991). (A similar model was proposed earlier by Kent and Mardia, 1988, but without consideration of medical imaging applications.) This approach 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. We call these weighting factors partial volume coefficients (PVCs). The method involves maximum-likelihood estimation of the PVCs for each voxel that model PV fractions of pure tissue types. Some authors have studied the identification of voxels containing PVE based on the mixel or a closely related model without trying to estimate the PVCs for each voxel Laidlaw et al., 1998, Ruan et al., 2000, Santago and Gage, 1993. Our interest in this study is in estimating PVCs and not in merely identifying voxels containing PVE.

Before statistical PV estimation can be performed, the probability density functions (pdfs) of the RVs describing pure tissue types must be specified. In practice, this typically means that the parameters of the pdfs—usually normal distributions—need to be estimated. Unfortunately, errors in the parameter estimation often have a major impact on the quality of PV estimation. In the case of statistical “hard” classification, where each voxel is classified according to the most dominant tissue type, the parameter estimation problem differs from the one studied here as each intensity is a realization of a single RV. On the other hand, the mixel model assumes that each voxel's intensity represents a weighted sum of several RVs and aims to estimate these unknown weight parameters.

In general, there are three approaches to the parameter estimation problem: histogram analysis (Santago and Gage, 1993), simultaneous parameter, and partial volume estimation by expectation maximization (EM)-like algorithms (Noe and Gee, 2001), and estimation based on a hard segmentation of the image (Shattuck et al., 2001). These three approaches, however, each have their drawbacks. Histogram analysis requires a mixture probability density to be fit to an image histogram by parameter optimization. This involves finding the minimizer of a multimodal objective function and therefore reliability of histogram analysis for parameter estimation depends heavily on the optimization algorithm used for the fitting task. If a standard nonlinear optimization algorithm aimed at local minimization (e.g., Levenberg–Marquee algorithm used commonly for curve fitting) is used, the initialization for the algorithm has to be chosen carefully to avoid convergence to a poor local minimum. These considerations call for the use of advanced global optimization algorithms, for example, Santago and Gage (1993) propose to use the tree-annealing method (Bilbro and Snyder, 1991). The problem with global optimization methods is that they are usually far more time consuming than local optimization methods.

Like histogram analysis, expectation–maximization schemes for parameter estimation are time consuming. Besides, the use of spatial information in the form of Markov Random Fields Besag, 1974, Besag, 1986 causes practical problems with the E-step of the EM-algorithm (Van Leemput et al., 2003). To solve this, statistical dependency between voxel labels can be ignored during the E-step as in Noe and Gee (2001). However, this leads to an algorithm that is merely a heuristic and does not necessarily share the convergence properties of the original EM-algorithm (Dempster et al., 1977). An elegant solution to this problem based on a Monte Carlo EM-algorithm (Wei and Tanner, 1990) was recently proposed by Van Leemput et al. (2003). The algorithm was implemented only for 2D image slices and the authors reported that partial volume estimation typically requires about 20 min for a single slice of an MR image with a very fast 1.7 GHz processor. Hence, the total time consumption for an (moderately sized) image of 100 slices would be over 30 h. The authors claimed that the computation time can be significantly reduced but did not offer any figures to support this claim.

Parameter estimation based on a hard labeling can be computationally efficient due to the prior knowledge of labels of voxels that can be utilized in parameter estimation. However, due to the PVE and classification errors in the hard labeling, each class in the hard segmented image contains a large number of outliers. This fact is taken into the account in Shattuck et al. (2001), but their approach for the parameter estimation applies only to single-spectral images and assumes that each tissue type has a Gaussian noise distribution with the same variance. Also, since this method involves detecting modes of histograms, it is sensitive to noise and the estimates may not be unique.

In this study, using minimum volume ellipsoid and minimum covariance determinant estimators Rousseeuw, 1984, Rousseeuw and Leroy, 1987, we propose routines for parameter estimation based on segmented images that are well-defined and can be used in single-spectral as well as multispectral cases. The studied estimators are robust; in other words, they tolerate deviations from the parametric form of the distribution assumed for the data. This is important because in our case, part of the data is good (correctly classified voxels of pure tissue), but the data contain also outliers (e.g. PV voxels). There are fast algorithms for computing the estimators and hence increasing the robustness of the parameter estimation does not lead to significantly increased running time for the PV estimation. We compare different techniques of parameter estimation using simulated MR data (Kwan et al., 1999) and demonstrate how the errors in the parameter estimation affect the final PV estimation results. The results are also compared against those obtained from an EM-like method, similar to Noe and Gee (2001), but with an advanced initialization technique. Furthermore, we compare our results with a fast PV estimation technique proposed by Shattuck et al. (2001) and demonstrate that our technique can yield clear improvements in the accuracy of PV estimates without considerable loss in time efficiency. As a secondary contribution, we consider differences between material and sampling noise models for PVE also on a more theoretical level (cf. Santago and Gage, 1995).

Section snippets

Statistical model for the partial volume effect

In this section, we state the PVE parameter estimation problem and describe the mixel model more formally. In the following, random variables (RVs) are denoted by boldface letters, while both scalars and vectors are shown in italics. Let us denote the observed image by X = {xi: i = 1,…, N}, with xiRK, and K the number of data channels. Let the set of possible tissue types present in the image be L = {1,…, M}. Moreover, lj is the RV describing the tissue type j and the pdf of lj is Gaussian g(·∣

Simulated data

Different methods for parameter estimation and the influence of the quality of parameter estimates to PVC estimates were studied using the BrainWeb Simulated Brain Database of Montreal Neurological Institute (Cocosco et al., 1997) that is available at http://www.bic.mni.mcgill.ca/brainweb. The images in the database are generated by an MRI simulator (Kwan et al., 1999), which models the MRI data acquisition process starting from the Bloch equation. The input for simulations is a fuzzy realistic

Parameter estimation

Results for the single-spectral case are presented in Table 2. Robust estimators combined with trimming were the most reliable among the segmentation-based parameter estimation techniques. They were also better than parameter estimates that were used to initialize EM-ICM and usually better than the final results of EM-ICM. For example, with image set 1 and 5% noise, TMCD produced a Mahalanobis error of 0.06, whereas the error resulting from MCD was 0.19 and the error resulting from EM-ICM was

Discussion

In this paper, various methods for the parameter estimation for a statistical model of the partial volume effect have been studied. It has been shown that it is possible to estimate parameters in a reliable and fast way based on the initial hard labeling of the image before the actual partial volume estimation procedure. For this, outliers of each class in the classified image were eliminated using a simple morphological rule, and thereafter, parameters for model were computed by using robust

Acknowledgements

Fig. 1 was generously provided by Jason Lerch, McConnell Brain Imaging Centre, Montreal Neurological Institute. Thanks to Chris Cocosco, Jason Lerch, and Steve Robbins for help in implementing the algorithms for this paper. J. Tohka acknowledges financial support from the Tampere Graduate School in Information Science and Engineering, the Academy of Finland, the KAUTE foundation, and the Jenny and Antti Wihuri fund.

References (46)

  • C. Cocosco et al.

    Brainweb: online interface to a 3D MRI simulated brain database

  • C. Cocosco et al.

    Automatic generation of training data for brain tissue classification from MRI

  • L. Collins et al.

    Design and construction of a realistic digital brain phantom

    IEEE Trans. Med. Imag

    (1998)
  • D.L. Collins et al.

    ANIMAL+INSECT: improved cortical structure segmentation

  • D.L. Collins et al.

    Automated estimation of brain volume in multiple sclerosis with BICCR

  • A. Dempster et al.

    Maximum likelihood from incomplete data via the EM algorithm

    J. R. Stat. Soc., Ser. B Methodol

    (1977)
  • E.R. Dougherty

    Probability and Statistics for the Engineering, Computing and Physical Sciences

    (1990)
  • S. Geman et al.

    Stochastic relaxation, gibbs distributions and the Bayesin restoration of images

    IEEE Trans. Pattern Anal. Mach. Intell

    (1984)
  • F.R. Hampel et al.

    Robust Statistics. An Approach Based on Influence Functions

    (1985)
  • M. Kamber et al.

    Model-based 3-D segmentation of multiple sclerosis lesions in magnetic resonance brain images

    IEEE Trans. Med. Imag

    (1995)
  • J.T. Kent et al.

    Spatial classification using fuzzy membership models

    IEEE Trans. Pattern Anal. Mach. Intell

    (1988)
  • Kollokian, V., 1996. Performance analysis of automatic techniques for tissue classification in magnetic resonance...
  • R.-S. Kwan et al.

    MRI simulation-based evaluation of image-processing and classification methods

    IEEE Trans. Med. Imag

    (1999)
  • Cited by (568)

    View all citing articles on Scopus
    View full text