Segmentation of magnetic resonance brain images through discriminant analysis

Journal of Neuroscience Methods 131 (2003) 65–74

Segmentation of magnetic resonance brain images through discriminant analysis Umberto Amato a,∗ , Michele Larobina b , Anestis Antoniadis c , Bruno Alfano b a

b

Istituto per le Applicazioni del Calcolo ‘Mauro Picone’ CNR—Sezione di Napoli, Consiglio Nazionale delle Ricerche, Via Pietro Castellino 111, Napoli 80131, Italy Istituto di Biostrutture e Bioimmagini, Consiglio Nazionale delle Ricerche, Via Pansini 5—Edificio 10, Napoli 80131, Italy c Institut d’Informatique et Mathématiques Appliquées de Grenoble, Université J. Fourier, 51 Rue des Mathematiques, B.P. 53, 38041 Grenoble Cedex 09, France Received 2 April 2003; received in revised form 25 June 2003; accepted 21 July 2003

Abstract Segmentation (tissue classification) of medical images obtained from a magnetic resonance (MR) system is a primary step in most applications of medical image post-processing. This paper describes nonparametric discriminant analysis methods to segment multispectral MR images of the brain. Starting from routinely available spin-lattice relaxation time, spin–spin relaxation time, and proton density weighted images (T1 w, T2 w, PD w), the proposed family of statistical methods is based on: (i) a transform of the images into components that are statistically independent from each other; (ii) a nonparametric estimate of probability density functions of each tissue starting from a training set; (iii) a classic Bayes 0–1 classification rule. Experiments based on a computer built brain phantom (brainweb) and on eight real patient data sets are shown. A comparison with parametric discriminant analysis is also reported. The capability of nonparametric discriminant analysis in improving brain tissue classification of parametric methods is demonstrated. Finally, an assessment of the role of multispectrality in classifying brain tissues is discussed. © 2003 Elsevier B.V. All rights reserved. Keywords: Magnetic resonance imaging; Brain; Segmentation; Discriminant analysis; Principal component analysis; Independent component analysis; Kernel density estimation

1. Introduction Nuclear magnetic resonance is recognized as a very important tool for clinical diagnosis. The main task in MR imaging is to infer useful information for diagnostic purposes from the images. From the methodological point of view this problem is called classification. From the voxel information one wants to understand the type of tissue associated to that voxel, both for making morphometric analyses, and for picking pathological conditions somewhere. Traditionally, this task was accomplished by radiologists, however there is a strong tendency to complement this activity by computer-based methodologies and quality and speed of the algorithms that actually process MR images are being quickly improved. The result of these efforts, combined with

∗

Corresponding author. Tel.: +39 0816132377; fax: +39 0816132597. E-mail address: [email protected] (U. Amato).

0165-0270/$ – see front matter © 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0165-0270(03)00237-1

an appropriate choice of the measurement protocol, is an increased reliability of the classification methodology in terms of success of pathological conditions retrieved. Among the technological progresses candidate to improve classification of MR images, the most noteworthy is the availability of so-called multispectral images: several independent images of the same anatomical slice are provided that increase the capability of discrimination between different tissues as a consequence of different responses of the tissues to particular pulse sequences (Fletcher et al., 1993; Jackson et al., 1994; Vannier et al., 1985). At present there are several methodologies for tissue classification of multispectral MR images: discriminant analysis (e.g. DeCarli et al., 1995; Zavaljevski et al., 2000), fuzzy methods (e.g. Suckling et al., 1999), neural networks (e.g. Alirezaire et al., 1998; Lin et al., 1996; Reddick et al., 1997), atlas methods (e.g. Nakazawa and Saito, 1994), knowledge-based techniques (e.g. Clark et al., 1998), shape-based methods (e.g. Duta and Sonka, 1998; Suri,

66

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

2001), variational methods (e.g. Pien and Gauch, 1994; Tang et al., 2000), morphological (texture) based techniques (Chen et al., 1993), multivariate principal component analysis methods (Antalek et al., 1998), statistical pattern recognition (e.g. Andersen et al., 2002; Bezdek et al., 1993; see also Clarke et al., 1995 and Niessen et al., 1999 for a review). The present paper deals with the discriminant analysis that is well-known in the statistics research and that has been proven successful in many applications arising in most research fields, including medical diagnostics and brain MR imaging. The main ingredient of discriminant analysis is the estimation of the probability density function of the MR signals for each normal brain tissue that can be detected in a MR scan. Most discriminant analysis methods are based on a parametric estimate of the density functions, where an analytic shape is given for the densities and a few parameters are sought (Student (Holden et al., 1995), Gaussian (Zavaljevski et al., 2000); see also Gudbjartsson and Patz, 1995 for Rician distributions). As well known departure from this assumption can dramatically degrade accuracy of the classification, an alternative is given by nonparametric methods, where the analytic expression of the density is unknown. In this case, an important issue is the so called curse of dimensionality that degrades the convergence rate of the density estimator to the actual density in the presence of multivariate data, i.e. multispectral data. This point is strictly related to how the multivariate character of the data is faced: a simple yet effective rule would be to treat each component separately; however this does not take account of the intrinsic lack of independence of the MR images (e.g. T1 w, T2 w, PD w) that would justify this approach rigorously. Therefore, a correct way to face the problem is to perform a priori transform of the T1 w, T2 w, PD w triplets into a target domain where new coordinates can be considered independent from a statistical point of view. In this respect, Amato et al. (2003) propose a discriminant analysis methodology aimed at solving all these issues based on (i) a (linear) transform of the original components (i.e. images) into independent components by use of the independent component analysis (ICA) tools; (ii) univariate nonparametric estimate of the density function for each component separately by means of a Kernel estimator; (iii) a classical optimal Bayes decision rule for 0–1 loss. We have to notice that the idea of transforming the original signals has quite a history in the signal processing literature and also in brain MR imaging. Wang et al. (2001) consider a Generalized Orthogonal Subspace Projection (GOSP) methodology to estimate the linear transform. Under some respects, the approach proposed by Wang et al. (2001) is opposite to the one followed in Amato et al. (2003): a dimensionality expansion is pursued in Wang et al. (2001) aimed at generating nonlinearly correlated images from the input ones, whereas independency is the key word of Amato et al. (2003). Other transforms were considered in Soltanian-Zadeh et al. (1992). We also mention that a nonlinear transform is considered in Alfano et al. (1997) aimed at yielding the apparent relaxation times

of the tissues (R1 = 1/T1 , R2 = 1/T2 ) starting from the MR signals, T1 w, T2 w. The rationale for this last transform is the independence of the derived quantities on the experimental protocol. The aim of the present paper is to follow a nonparametric approach to the discriminant analysis for the classification of MR brain images. Furthermore, the role of transforming images prior to classification to yield independent images will be investigated with an emphasis on nonlinear transforms.

2. Materials 2.1. Discriminant analysis The purpose of any discriminant analysis method is to assign an observation x (voxel) to one class from a set of K classes (brain tissues) with the lowest possible error rate. In the standard setting, observations are described by multivariate random vectors coming from a certain class k (k = 1, . . . , K) characterized by a density function fk (x). An observation is assumed to be drawn from one and only one class and error is incurred if it is assigned to a different one. The cost or loss associated with such an error is usually defined ˜ 1 ≤ k, k˜ ≤ K, where k is the correct class asby L(k, k), signment and k˜ is the assignment that was actually made. A special but commonly occurring loss L is the 0–1 loss defined by
˜ 0, if k = k, ˜ L(k, k) = 1, otherwise. In this case, the Bayes decision rule actually used in the ˜ such that fk (x)πk is maxalgorithm allocates x to class k, imum, where πk are unconditional class prior probabilities approximated by tissue volume fractions in the present paper. Discriminant analysis requires a training data set that can be considered as a sample of feature vectors from each class used to learn the density functions of the classes. 2.2. Probability density function The mostly applied classification rules are based on the normal-theory classification, which assumes that the densities fk are Gaussian. Such standard parametric rules include linear and quadratic discriminant analysis (e.g. Anderson, 1984), which have been shown to be quite useful in a wide variety of problems. However, in practice, the form of the class-conditional densities is seldom known. A deep analysis of the density functions of the brain tissues is very important for the correct application of the discriminant analysis. First of all, overall feasibility of the discriminant analysis in classifying voxels can be inferred: the more density functions are separated for different tissues, the more the discriminant analysis will be able to classify voxels correctly. Then, if density functions are recognized as

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

belonging to classic known types, simple parametric discriminant analysis methods can be applied; in the opposite case discriminant analysis naturally extends to the situation where nothing is known about the densities fk except possibly for some assumptions about their general behavior. The suggested approach is to estimate the densities fk from a training set using nonparametric density estimates and to substitute these estimates into the Bayes decision rule to give a nonparametric discriminant rule. The most used procedure for nonparametric density estimation is kernel density estimation with appropriate smoothing parameter selection: N

k 1  fk (x) = H(x − z ), Nk

=1

where Nk is the size of the training set for class k, z ≡ (p) (1) (z )=1,... ,Nk ≡ (z , . . . , z )=1,... ,Nk (dropping the dependence on the class for the multispectral components for simplicity’s sake). A popular choice of the Kernel function H is the Gaussian one:   p  1 (x(j) )2 H(x) = , exp − √ 2σj2 σ1 · · · σp 2π j=1 where σj = 1.09σˆ j Nk−0.2 and σˆ j2 is an estimate of the variance of the multispectral component j for the (dropped) class k. In the present paper, density estimation is computed by the algorithm of Beardah (1995). Most often the p-variate density function is taken to be a product of univariate functions. However, it is known that while in one-dimensional density estimation it is not crucial to estimate the tails accurately, this is no longer true in high dimensional spaces where regions of relatively low density can still be extremely important parts of the multidimensional density. 2.3. Transforms To estimate a density function as the product of univariate functions requires an assumption of independence of data that does not always hold in practice. In the case of Gaussian distributions, this can be fixed by transforming original multivariate data into principal components by principal component analysis (PCA) and then proceeding with classification. For general distributions, application of PCA to the data only decorrelates them, without yielding full independence. A possible remedy is to seek for a transform that makes data mutually independent irrespective of their distribution. ICA achieves such a task. It is a statistical method for linearly transforming an observed multidimensional random vector into a random vector whose components are stochastically as independent from each other as possible. Several procedures to find such transformations

67

have been recently developed in the signal processing literature relying either on Comon’s information-theoretic approach (Comon, 1994) or Hyvärinen’s maximum negentropy approach (Hyvärinen, 1997). The present paper considers this latter approach and relies on the Matlab fastica package (Hyvärinen, 1999), available at the Website http://www.cis.hut.fi/projects/ica/fastica/, for its implementation. 2.4. Classification Three nonparametric and one parametric discriminant analysis methods for classifying multispectral MR brain images have been considered: • Linear discriminant analysis (LDA), based on Gaussian density functions with common variance among classes; • Nonparametric discriminant analysis (NPDA), where a nonparametric estimate of the density functions is made for each component separately; • Principal component discriminant analysis (PCDA), where original components are transformed into principal components prior to nonparametric density estimation; • Independent component discriminant analysis (ICDA), where original components are transformed into independent components prior to nonparametric density estimation. Two kinds of experiments have been conducted: synthetic and real data. Discriminant analysis requires two data sets for each experiment: a training set that is fed into the classification method to learn about tissue density functions and the test set on which classification is actually performed. The following indicators have been used to rank the effectiveness of the discriminant analysis methods to classify MR images; all of them are evaluated for the whole volume of data: • S: global agreement percentage defined as the percentage of correctly classified voxels; index S can refer to each tissue or cumulatively to the whole volume; • F+ : false positive rate for each tissue, defined as the percentage of voxels erroneously classified as belonging to that tissue; • F− : false negative rate for each tissue, defined as the percentage of voxels erroneously not classified as belonging to that tissue; • κ: kappa coefficient (Zijdenbos et al., 1994). It is a chance-corrected measure of agreement, defined by κ = (p0 − pc )/(1 − pc ), where p0 is the observed percentage of agreement (the percentage of correctly classified voxels) and pc is the percentage agreement that would occur by chance alone; values of κ > 0.7 are claimed to indicate a good classification; • WM : volume percentage of white matter; • GM : volume percentage of gray matter; • CSF: volume percentage of cerebro spinal fluid.

68

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

2.5. MR synthetic data A digital phantom to test, validate, and estimate performance of the classification methods was used as described by Kwan et al. (1996, 1999), Cocosco et al. (1997), Collins et al. (1998). It is based on a high-resolution (1 mm3 ) 3D low-noise data set created by registering 27 scans (T1 -weighted gradient-echo acquisitions with TR = 18 ms, TE = 10 ms, and flip angle 30◦ ) of the same normal volunteer in stereotaxic space where they were subsampled and intensity averaged. The volume contains 181 × 217 × 181 voxels and covers the brain completely. Full details of the digital brain phantom are given in Collins et al. (1998). The digital brain phantom is available at the website http://www.bic.mni.mcgill.ca/brainweb. Nine classes are defined in brainweb: cerebro spinal fluid, gray matter, white matter, fat, muscle/skin, skin, skull, glial matter, connective; tissue assignment for each voxel of the data set (needed for the training data set) is described in Kwan et al. (1999). Images are provided normalized to the value 255 for each signal. 2.6. MR real data MR studies were performed at 1.5 T (Siemens) with sequential acquisition of two interleaved sets of oblique axial slices (15 slices per set, thickness 4 mm, slice interval 4 mm). Two spin-echo sequences (TR = 600 ms, TE = 15 ms; TR = 2300 ms, TE = 15/90 ms; both with 90◦ flip angle, 256 × 256 matrix, 0.98 mm pixel size) for each set of slices were acquired sequentially in the same spatial position to obtain T1 w, T2 w, and PD w images. Slice positioning and angle were chosen on a sagittal localizer T1 w scan to cover the entire brain, including cerebellum. Eight studies from patients with no evidence of pathology have been considered. In these studies, we verified the absence of motion artifacts caused by head movement either during the acquisition or in the time interval between the two sequences. Nine classes are defined: partial volume bone-low density tissues (PV bone), fat, muscle, gray matter, white matter, pallidus, putamen, cerebro spinal fluid (CSF), and partial volume gray matter-CSF (PV gray-CSF). Tissue assignment for each voxel of the data set (needed for the training data set) is described in Alfano et al. (1997).

3. Results 3.1. Synthetic data We considered the T1 w, T2 w, and PD w images available in the framework of the brainweb project for a gradient-echo configuration, no noise, 1 mm vertical thickness (over 8 million voxels) to reduce the effect of mixed voxels. Fig. 1 shows the scatterplot of the triplets (T1 w, T2 w, PD w) projected along the (T1 w, T2 w), (T1 w, PD w), and (T2 w, PD w)

axes for the whole volume, colored according to the nine different tissues defined in the brainweb project. The voxel classification is available in the brainweb database. This kind of plot gives a quite good indication, even though not exhaustive, of the full multidimensional structure of the data. Tissues that show a strong overlap in all the three scatterplots are prone to poor classification by discriminant analysis. Fig. 2 shows the density functions of each tissue and each component, grouped by T1 w, T2 w, and PD w components. Analysis of density functions enforces the conclusions that can be drawn by scatterplots and allows one to investigate in more detail the overlap of density functions. An experiment of tissue classification was performed, where the brainweb data set has been considered both as a training and test data set and LDA, NPDA, and PCDA methods have been used (results for ICDA are not shown for lack of convergence of the fastica package). Table 1 shows the performance indicators for each tissue and the whole brain volume as resulting from the experiment. 3.2. Real data Estimates of voxel density functions and results of classification experiments are presented. Differently from the synthetic case a transform of T1 w, T2 w and PD w to apparent R1 , R2 , and PD has been done as in Alfano et al. (1997); due to the nonlinear character of the transform, it is expected that the density functions will be very different from the ones shown for the synthetic case (brainweb) and generally more complicated. Fig. 3 shows the scatterplots of (R1 , R2 , PD ) triplets (for the whole volume) projected along the (R1 , R2 ), (R1 , PD ), and (R2 , PD ) axes for the pat3 patient; scatterplots for the other patients are similar. Class assignment of voxels among the defined tissues relies on the method fully described in Alfano et al. (1997). Fig. 4 shows the density functions estimated for the defined tissues grouped together by R1 , R2 , and PD components for the whole volume (over 500,000 voxels). Three classification experiments were performed. First, we considered each of the eight available patients as a training set and then each of the eight patients as a test set to be classified, yielding 64 classification runs globally. The results are summarized in Table 2, where for each patient considered as a training set the average agreement percentage, S, of the four classifiers over each patient considered as a test set is shown together with its standard deviation. The percentage agreement when the test set coincides with the training set is also reported for each patient; this corresponds to the best percentage agreement achieved, with the notable exceptions of pat2 and pat6 patients as a training set, where top performance of the nonparametric discriminant analysis is reached for the pat1 and pat5 as test data, respectively. As a second experiment we consider in more detail the case of pat3 as a training and test data set. Table 3 shows the resulting performance indicators. Fig. 5 shows T1 w, T2 w, and PD w multispectral images for one brain axial slice at the level of

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

Fig. 1. Scatterplot of the T1 w, T2 w, PD w brainweb multispectral MR images projected along the couples of axis (T1 w, T2 w), (T1 w, PD w), (T2 w, PD w), with different colors assigned to the pixels according to the type of tissues.

69

Fig. 3. Scatterplot of the R1 , R2 , PD multispectral MR images for patient pat3 projected along the couples of axis (R1 , R2 ), (R1 , PD ), (R2 , PD ), with different colors assigned to the pixels according to the type of tissues.

the basal ganglia from the same pat3 data set. Fig. 6 shows their segmented images obtained by the method of Alfano et al. (1997), LDA, NPDA, and ICDA. Finally we consider an experiment where classification is accomplished basing on single multispectral components or their couples from the data set. The purpose of this analysis is to assess the performance of each component in the classification analysis and to rank them accordingly. Table 4 shows the resulting tissue agreement percentage, S, in the same case previously analyzed (pat3 as training and test data set and NPDA classifier).

4. Discussion 4.1. Synthetic data

Fig. 2. Density functions of the tissues in brainweb MR images grouped by T1 w, T2 w, and PD w.

Analysis of Fig. 1 shows that not all tissues are well separated in the scatterplot of the brainweb data set. This is the case for the glial matter that strongly overlaps with the gray matter and for the connective tissue that overlaps with the white matter in the multiparametric space. Fig. 2 confirms these overlaps, which causes tissues with small volume fractions (glial matter and connective tissue) to be severely underestimated in the classification by discriminant analysis. Analysis of Fig. 2 also shows that even for synthetic data not all density functions have a regular shape and that Gaussianity cannot be taken as a general type of density. The experiment on tissue classification (Table 1) shows that performance of the LDA is poor with respect to NPDA and PCDA, as expected from the analysis of the tissue density functions of the voxels. The most improvement is due to the choice of a nonparametric method for estimating density functions,

70

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

Table 1 Performance indicators for the LDA, NPDA, and PCDA methods applied to the brainweb data set: agreement percentage S for the nine tissues and for the whole volume; kappa coefficient κ; white and gray matter, and CSF volume percentage LDA

CSF Gray matter White matter Fat Muscle/Skin Skin Skull Glial matter Connective

All

NPDA

S

κ

77.2 80.6 34.2 95.0 92.5 74.1 94.2 77.2 89.2

0.65 0.85 0.45 0.93 0.90 0.73 0.95 0.06 0.48

46.6 3.2 2.7 7.8 9.8 16.4 3.4 2200 146

22.8 19.4 65.8 5.0 7.5 25.9 5.8 22.8 10.8

69.8 95.3 94.1 91.8 86.7 81.3 97.0 0.0 59.3

0.66 0.92 0.88 0.89 0.88 0.77 0.97 0.00 0.65

S

κ

GM

WM

CSF

S

κ

GM

WM

75.7

0.73

51.5

24.2

86.3

0.85

45.6

35.8

16.9

F−

S

PCDA

F+

κ

F+

F−

S

κ

F+

F−

32.3 8.3 14.7 13.3 6.7 20.4 2.9 0.0 16.9

30.2 4.7 5.9 8.2 13.3 18.7 3.0 100 40.7

65.8 93.4 97.9 94.7 92.6 78.7 97.6 0.0 71.4

0.65 0.93 0.90 0.95 0.90 0.75 0.95 0.0 0.77

28.5 4.5 15.2 4.5 8.8 19.5 7.6 0.5 10.4

34.2 6.6 2.1 5.3 7.4 21.3 2.4 100 28.6

CSF

S

κ

GM

WM

CSF

87.6

0.86

44.1

38.1

18.7

17.7

True values of gray and white matter and CSF percentage are 46.2, 34.6 and 19.2%, respectively.

whereas improvement due to prior transforms into principal components yields a minor increase of performance. Estimate of the gray and white matter and CSF percentage is quite good, in particular with NPDA, where the error is about 1% (true values are 46.2, 34.6, and 19.2%, respectively). Connective tissue and glial matter, especially, have

Fig. 4. Density functions of the tissues in pat3 MR images grouped by R1 , R2 , and PD .

a poor agreement percentage and transformation into principal components improves the performance of the former.

4.2. Real data Analysis of Fig. 4 shows much less overlap of voxel distributions among different tissues than in the case of brainweb tissue definition, although different tissue classes are defined. This suggests a better link of the definition of tissue classes to the actual capabilities of MR images to discriminate among different tissues. Less overlap is to be intended in the sense that distributions overlap little with respect to at least one of the three possible couples of the R1 , R2 , and PD coordinates. More than the case of synthetic data, density functions show a quite complicated behavior, frequently with bimodal distributions, well far from classic regular types. This makes parametric discriminant analysis potentially not well suited for classifying MR images. Also note that some distributions totally overlap in some (R1 , R2 , PD ) projections, but are separated in others, confirming the role of multispectrality in classifying MR images. Analysis of the tissue classification experiments of Table 2 confirms that nonparametric discriminant analysis methods significantly outperform parametric LDA, but the transformation into new coordinates is not able to improve the performance indicators significantly. Classification of a test set starting from a different training set degrades the percentage agreement by around 5% on average (7% for LDA) with respect to the case of coinciding training and test data sets. Therefore, the method is quite robust with respect to the choice of the training set. From Table 3, we see that global percentage agreement of nonparametric discriminant analysis is significantly better than LDA and that the percentage agreement is more homogeneous among tissues for the latter (it ranges between 74.7 and 90.7% for LDA and, e.g. between 61.8 and 93.1% for NPDA). Interestingly, even though the percentage agreement, S, decreases when moving to

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

71

Table 2 Average percentage agreement (with standard deviation) of the four classification methods when one patient is considered as training set and each of the eight patients is classified Training patient

LDA

NPDA

Average Pat0 Pat1 Pat2 Pat3 Pat4 Pat5 Pat6 Pat7

79 73 65 78 78 75 74 79

± ± ± ± ± ± ± ±

4 8 9 3 4 7 7 3

Test = Train

Average

83.4 82.2 78.2 83.1 82.6 83.8 83.3 82.4

87 84 80 87 87 85 84 87

± ± ± ± ± ± ± ±

PCDA

3 5 7 3 3 5 5 3

Test = Train

Average

90.4 89.1 87.2 89.9 89.6 90.7 89.6 90.0

86 83 80 87 87 85 85 87

± ± ± ± ± ± ± ±

ICDA

4 5 7 3 3 5 5 3

Test = Train

Average

Test = Train

90.5 89.4 87.8 89.8 90.0 90.8 89.4 90.0

86 ± N.A. 80 ± 87 ± N.A. 85 ± 85 ± 87 ±

90.3 N.A. 87.7 90.3 N.A. 91.1 89.7 90.3

3 7 3 5 5 3

Percentage agreement of the case when training and test sets coincide is also shown. N.A. (not available) means that the fastica package was not able to reach convergence for the corresponding training patient. Table 3 Performance indicators for the LDA, NPDA, PCDA, and ICDA methods applied to the pat3 patient data as training and test set: agreement percentage S for the nine tissues and for the whole volume; kappa coefficient κ; white and gray matter, and CSF volume percentage LDA

PV bone Fat Muscle Gray White Pallidus Putamen CSF PV gray-CSF

All

NPDA

PCDA

S

κ

F+

90.7 89.0 77.5 81.5 81.6 79.0 74.7 80.1 85.8

0.76 0.83 0.80 0.82 0.85 0.43 0.50 0.70 0.30

35.6 19.9 11.0 4.7 4.0 185 122 46.0 364

9.3 11.0 22.5 18.5 18.4 21.0 25.3 19.9 14.2

83.9 87.6 89.3 93.1 91.8 73.2 77.6 91.7 61.8

0.87 0.89 0.87 0.88 0.90 0.54 0.58 0.80 0.57

S

κ

GM

WM

CSF

S

κ

GM

WM

CSF

83.1

0.81

49.8

11.6

89.9

0.89

54.0

36.4

6.5

34.3

F−

S

κ

F+

F−

S

5.7 7.1 12.2 8.8 7.3 97.2 88.4 36.7 53.0

16.1 12.4 10.7 6.9 8.2 26.5 22.4 8.3 38.2

ICDA κ

F+

F−

S

κ

F+

F−

82.3 85.8 90.0 93.6 93.8 63.1 44.3 87.4 60.0

0.84 0.88 0.84 0.90 0.90 0.56 0.50 0.82 0.65

9.0 6.3 17.8 7.1 9.3 59.9 32.9 25.5 23.2

17.7 14.2 10.0 6.4 6.2 36.9 55.7 12.6 40.0

84.1 86.5 91.0 94.5 93.2 66.4 52.5 81.3 64.7

0.85 0.88 0.86 0.90 0.90 0.59 0.54 0.82 0.68

8.7 6.7 16.2 7.4 8.4 59.7 40.1 16.6 25.3

15.9 13.5 9.0 5.5 6.8 33.6 47.5 18.7 35.3

S

κ

GM

WM

CSF

S

κ

GM

WM

CSF

89.8

0.89

54.1

38.4

5.7

90.3

0.89

54.7

38.1

5.2

True values of gray and white matter and CSF percentage are 54.7, 37.9, and 5.3%, respectively.

nonparametric methods for selected tissues (Pallidus, Putamen, and PV gray in particular), the κ coefficients always assume better values. The reason of this behavior is a better balance of misclassifications of tissues between false positive and negative occurrences, especially for ICDA.

Finally we discuss the positive effect of multispectrality as resulting from Table 4. In particular R1 data strongly contribute to discrimination of gray and white matter and fat; R2 to muscle and CSF, besides gray and white matter; PD to PV bone mainly. It is noteworthy that for tissues whose

Fig. 5. MR multispectral images T1 w (left), T2 w (center), and PD w (right) for one brain axial slice (pat3) at the level of the basal ganglia.

72

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

Table 4 Tissue agreement percentage S for the patient pat3 as training and test data set when single multispectral components or their couples are considered in the NPDA classifier R1

R2

PD

R1 −R2

R1 −PD

R2 −PD

R1 −R2 −PD

PV bone Fat Muscle Gray White Pallidus Putamen CSF PV gray-CSF

10.4 81.1 6.6 95.6 91.5 0.0 0.0 0.0 0.0

5.0 24.7 85.4 74.9 67.8 0.0 0.0 89.9 0.0

86.3 54.6 16.1 78.7 57.1 0.0 0.0 0.0 0.0

22.3 80.7 89.9 92.1 93.9 28.0 2.4 86.2 42.7

73.1 85.4 44.8 94.7 89.2 0.0 26.4 54.2 31.6

80.6 60.7 75.2 88.8 72.8 9.9 34.0 84.7 0.2

83.9 87.6 89.3 93.1 91.8 73.5 77.6 91.7 61.8

All

64.2

57.9

59.9

79.5

80.1

77.4

89.9

As a reference in the last column the values of the full multispectral data are also given.

monospectral image fails in providing acceptable success percentages, use of multispectral images significantly improves performances, in particular pallidus and putamen. 4.3. Perspectives The present paper demonstrated feasibility of nonparametric discriminant analysis in classifying brain MR multispectral images and, in particular, their superiority with

respect to parametric methods. The ability to deal with density functions of general shapes also permits it to take partial account of mixed voxels that tend to spread and make tissue distributions more irregular. In addition, estimates of gray and white matter and CSF percentage by nonparametric discriminant analysis discussed in this paper are already quite good. This makes such crispy methods competitive with the fuzzy alternatives. Robustness of the method with respect to the choice of the training data set has also been shown. No

Fig. 6. Segmented images of the MR multispectral images of Fig. 5. Upper left: classification resulting from the method of Alfano et al. (1997); upper right: LDA; lower left: NPDA; lower right: ICDA. Color map is the same as Figs. 3 and 4.

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

clear evidence was demonstrated about the overall improvement of the classification by linearly decorrelating or making data as independent as possible; rather than improving the percentage agreement, transform into independent components has the effect of better balancing misclassifications among false positive and negative occurrences. The methods presented in this paper can be improved with the twofold objective of: (1) increasing the overall performance of nonparametric discriminant analysis and, (2) better fitting the methods to the brain segmentation problem. The first objective will be reached through nonlinear transforms into independent components. The second one will be based on a better estimate of the priori probabilities πk to be fed into the Bayes rule; they will take account of local textures that will force contiguity of voxels, thus mitigating the effect of isolated misclassifications, and of the brain morphology that will limit the number of local possible tissues. The result will be a fully local classification method. In addition, some recent powerful iterative techniques for the reduction of misclassifications will be considered (e.g. bagging and boosting; Breiman, 1996; Freund and Schapire, 1997).

Acknowledgements Part of this work has been supported by Grant #RF01.167 from the Italian Ministero della Sanità, “Programmi speciali” D. Lgs 502/92 e D. Lgs 229/99. References Alfano B, Brunetti A, Covelli EM, Quarantelli M, Panico MR, Ciarmiello A, et al. Unsupervised, automated segmentation of the normal brain using a multispectral relaxometric magnetic resonance approach. Magn Reson Med 1997;37:84–93. Alirezaire J, Jernigan ME, Nahmias C. Automatic segmentation of cerebral MR images using artificial neural networks. IEEE Trans Nucl Sci 1998;45:2174–82. Amato U, Antoniadis A, Grégoire G. Independent component discriminant analysis. Int J Math 2003;3:735–53. Anderson TW. An introduction to multivariate statistical analysis. 2nd ed. New York: Wiley, 1984. Andersen AH, Zhang Z, Avison MJ, Gash DM. Automated segmentation of multispectral brain MR images. J Neurosci Methods 2002;122:13– 23. Antalek B, Hornak JP, Windig W. Multivariate image analysis of magnetic resonance images with the direct exponential curve resolution algorithm (DECRA). Part 2. Application to human brain images. J Magn Reson 1998;132:307–15. Beardah CC. The Kernel Density Estimation toolbox for Matlab, 1995. Available at http://euler.ntu.ac.uk/maths.html. Bezdek JC, Hall LO, Clarke LP. Review of MR image segmentation techniques using pattern recognition. Med Phys 1993;20:1033–48. Breiman L. Bagging predictors. J Mach Learn 1996;24:123–40. Chen Y, Dougherty ER, Totterman SM, Hornak JP. Classification of trabecular structure in magnetic resonance images based on morphological granulometries. Magn Reson Med 1993;29:358–70. Clark MC, Hall LC, Goldgof DB, Velthuizen R, Murtagh FR, Silbiger S. Automatic tumor segmentation using knowledge-based techniques. IEEE Trans Med Imag 1998;17:187–201.

73

Clarke LP, Velthuizen RP, Camacho MA, Heine JJ, Vaidyanathan M, Hall LO, et al. MRI segmentation: methods and applications. Magn Reson Imag 1995;13:343–68. Cocosco C, Kollokian V, Kwan R-S, Evans A. Brainweb: Online interface to a 3-D MRI simulated brain database. Neuroimage 1997;5:S425. Collins DL, Zijdenbos AP, Kollokian V, Sled JG, Kabani NJ, Holmes CJ, et al. Design and construction of a realistic digital brain phantom. IEEE Trans Med Imag 1998;17:463–8. Comon P. Independent component analysis, a new concept? Signal Process 1994;36:287–314. DeCarli C, Murphy DGM, McIntosh AR, Teichberg D, Shapiro MB, Horwitz B. Discriminant analysis of MRI measures as a method to determine the presence of dementia of the Alzheimer type. Psychiatry Res 1995;57:119–30. Duta N, Sonka M. Segmentation and interpretation of MR brain images: an improved active shape model. IEEE Trans Med Imag 1998;16:1049– 62. Fletcher LM, Barsotti JB, Hornak JP. A multispectral analysis of brain images. Magn Reson Med 1993;29:623–30. Freund Y, Schapire RE. A decision-theoretic generalization of online learning and an application to boosting. J Comp Syst Sci 1997;55:119– 39. Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med 1995;34:910–4. Holden M, Steen E, Lundervold A. Segmentation and visualization of brain lesions in multispectral magnetic resonance images. Comput Med Imag Graphics 1995;19:171–83. Hyvärinen A. Independent component analysis by minimization of mutual information. In: Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP’97); 1997. p. 3917–20. Hyvärinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans Neural Networks 1999;10:626–34. Jackson EF, Narayana PA, Falconer JC. Reproducibility of nonparametric feature map for determination of normal human intracranial volumes with MR imaging data. J Magn Reson Imag 1994;4:692–700. Kwan RK-S, Evans AC, Pike GB. An extensible MRI simulator for post-processing evaluation. In: Proceedings of the Fourth International Conference on Visualization in Biomedical Computing, VBC’96, Hamburg; 1996. p. 135–40. Kwan RK-S, Evans AC, Pike GB. MRI simulation-based evaluation of image-processing and classification methods. IEEE Trans Med Imag 1999;18:1085–97. Lin JS, Cheng KS, Mao CW. Multispectral magnetic resonance images segmentation using fuzzy Hopfield neural network. Int J Biomed Comput 1996;46:205–14. Nakazawa Y, Saito T. Region extraction with standard brain atlas for analysis of MRI brain images. Image Process 1994;1:387–91. Proceedings of the IEEE International Conference on ICIP-94. Niessen WJ, Vincken KL, Weickert J, ter Haar Romeny BM, Viergever MA. Multiscale segmentation of three-dimensional MR brain images. Int J Comput Vision 1999;31:185–202. Pien HH, Gauch JM. Variational segmentation of multi-channel MRI images. In: Proceedings of the IEEE International Conference on Image Processing, Austin, Texas, November; 1994. 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 Trans Med Imag 1997;16:911–8. Soltanian-Zadeh H, Windham JP, Peck DJ, Yagle AE. A comparative analysis of several transformations for enhancement and segmentation of Magnetic Resonance scene sequences. IEEE Trans Med Imag 1992;11:302–18. Suckling J, Sdsson T, Greenwood K, Bullmore ET. A modified fuzzy clustering: algorithm for operator independent brain tissue classification of dual echo MR images. Magn Reson Imag 1999;17:1065–76.

74

U. Amato et al. / Journal of Neuroscience Methods 131 (2003) 65–74

Suri J. Two-dimensional fast Magnetic Resonance brain segmentation. IEEE Eng Med Biol 2001;July/August:84–95. Tang H, Wu EX, Ma QY, Gallagher D, Perera GM, Zhuang T. MRI brain image segmentation by multi-resolution edge detection and region selection. Comput Med Imag Graphics 2000;24:349– 57. Vannier MW, Butterfield RL, Jordan D, Murphy WA, Levitt RG, Gado M. Multispectral analysis of magnetic resonance images. Radiology 1985;154:221–4.

Wang CM, Yang SC, Chung PC, Chang CI, Lo CS, Chen CC, et al. Orthogonal subspace projection-based approaches to classification of MR image sequences. Comput Med Imag Graphics 2001;25:475–6. Zavaljevski A, Dhawan AP, Gaskil M, Ball W, Johnson JD. Multi-level adaptive segmentation of multi-parameter MR brain images. Comput Med Imag Graphics 2000;24:87–98. Zijdenbos AP, Dawant BM, Margolin RA, Palmer AC. Morphometric analysis of white matter lesions in MR images: methods and validation. IEEE Trans Med Imag 1994;13:716–24.