The automatic segmentation of three dimensional (3-D) magnetic resonance images (MRI) of the brain is increasingly important for neuroscientific studies. Studies concerning morphological measurements of different brain structures can benefit from the cost-effectiveness, speed and good reproducibility of automatic image segmentation [1
]. The basic form of the 3-D brain MRI segmentation involves the classification of the voxels in the MR image into the three primary tissue types: white matter (WM), gray matter (GM) and cerebro-spinal fluid (CSF). This tissue classification has various important applications. For example, several cortical surface extraction algorithms depend on tissue classification [2
The research problems in tissue classification methodology include how to best cope with the intensity non-uniformity, partial volume effect, and scanner specificity of MRI as well as how to best apply the existing neuroanatomical knowledge and the spatial dependency of the labels of close-by image voxels to increase the quality of tissue classification. Recently, information and/or classifier fusion for the tissue classification, or more generally for MRI segmentation, have drawn an increasing amount of attention. For example, an interesting line of research has focused on the simultaneous registration and tissue classification [8
]. These approaches can be viewed as information fusion methods where the information from the brain atlases is combined with the information on the tissue composition of the individual image in a non-trivial manner. Similarly, label propagation techniques for the brain MRI segmentation can use information from multiple templates to better adapt to the anatomical variability of human brain leading to decision fusion approaches [12
In this work, we propose a different kind of the information/classifier fusion approach targeted to the tissue classification problem. Below, we argue that the MR brain images are more accurately modeled locally than globally. We therefore divide the image domain into local subdomains that are described by their own local image model and derive a framework in which image models with local spatial support can be combined to yield a global, statistically based image model. Both local and global image models will be posed in a Markov Random Field (MRF) framework allowing the incorporation of the spatial information to improve the quality of the voxel classification. We apply a Sub Volume Probabilistic Atlas (SVPA) to define the local support regions so that the division of the image domain into sub-domains is anatomically informed. The model-fusion framework is complemented by computational algorithms for solving the parameters of the model and performing the tissue classification in a completely automated manner. Particularly, we apply a global optimization approach based on a genetic algorithm for unsupervised estimation of the tissue type parameters for each local model [13
]. We note that the approach presented here is in many senses complementary to the above referenced frameworks for joint image segmentation and registration. This is because, although we start from local image models, the final image model is a global generative model that could then be unified with a joint segmentation and registration approach.
As already mentioned, we argue that local image models can be more accurate than global image models. There are at least two important reasons for this. First, MR images contain low frequency intensity variations that are often described by the terms intensity non-uniformities or bias fields. In this work, we refer to this phenomenon as instrumentation non-uniformity. These low frequency variations cause non-stationarity of the parameters for the class conditional probability density functions (pdfs) of individual tissue classes. This has a negative effect on the quality of automatic tissue classifications that are based on the assumption that tissue classes have the same statistical properties throughout the image. Several methods for the correction of this low frequency noise in images have been proposed and a recent review can be found in [14
]. The methods for the retrospective bias field correction typically assume that the bias field is multiplicative and spatially smooth. Some of them, including [15
] are stand-alone methods, while others perform the tissue classification jointly with non-uniformity correction [19
]. With modern MR imaging devices, state-of-the-art methods for non-uniformity correction can be considered to be largely successful. However, these methods are not able to remove all of the low frequency noise in all parts of the image. In addition, there is evidence that the assumption of the spatial smoothness of the bias field may not be completely valid for T1 and T2 weighted scans [22
]; furthermore Marroquin [23
] presented a case against multiplicative nature of the bias fields. The framework we propose in this paper can reduce the unfavorable effects of unsuccessful intensity non-uniformity correction.
Second, the intrinsic properties of tissue types, such as T1 and T2 relaxation times and proton density (PD), vary across an individual brain. Since the MRI signal is a function of these intrinsic tissue characteristics, their variation across the brain causes variation in the MRI signal values as demonstrated in [24
]. This phenomenon - termed here as biological non-uniformity - has been confirmed in various studies. The variability of GM T1 relaxation time across the cortex is well documented [9
] as is the T1 relaxation time difference between cortical GM and various sub-cortical GM structures such as caudate, putamen, thalamus, amygdala and hippocampus [9
]. Differences in T2 between various GM and WM structures have been reported [27
]. There is evidence demonstrating a slight hemispheric asymmetry of PD as well as variability of PD in different WM structures [30
]. All the above findings are about healthy adult populations with the exception of [29
] where also schizophrenic subjects were studied. The biological non-uniformity can be even greater in pathology or in aging or d Since image intensities from a tissue class have changing statistical characteristics across the brain due to biological non-uniformity that may not be accounted for by methods designed to correct for instrumentation non-uniformity, the assumption of global models for tissue types has a negative effect on the quality of tissue classification and subsequent analysis.
] and Kovacevic [32
] have proposed tissue classification strategies based on similar arguments to those we present here, considering both biological and instrumentation non-uniformity. They first estimated parameters of a global mixture of Gaussians (MoG) based intensity model. Thereafter, the image domain was divided into smaller sub-domains whose intensities were modeled by individually by local MoGs. The parameters for these were estimated by a local optimization method starting from the global image parameter estimates. These previous approaches differ from our method in three important aspects. As opposed to our framework, these approaches did not include a spatial MRF prior and they did not lead to a generative, global image model. Moreover, partial volume effect (PVE) was not modeled in [32
] and in [31
] it was assumed that the PVE could be modeled using 50:50 mixtures, e.g., voxels containing 50% WM and 50% GM were included in the intensity model. That assumption is clearly overly simplistic.
Many methods for simultaneous instrumentation non-uniformity correction and tissue classification are similar to our method in that they consider tissue class parameters to be spatially varying over the image domain. These methods can include also an MRF prior for tissue classifications, and some of them incorporate statistically based modeling for PVE, e.g.,[33
]. The general difference between these methods and our framework is that in these the spatial variation of tissue class parameters is strongly parametrized whereas our method is nonparametric in this respect. Note that we do not argue that our method should be used to completely replace existing non-uniformity correction methods, but rather to complement them. The adaptive K-means algorithm [34
], modified and extended to MRI application in [35
] and the Adaptive MAP algorithm [37
] are perhaps closest to our framework. The adaptive K-means algorithm accounts for local intensity variations by an iterative procedure involving averaging over a sliding window whose size decreases as the algorithm progresses. These variations may be tissue class specific. The algorithm also includes an MRF prior. The adaptive MAP algorithm can be considered as an extension of the adaptive K-means, where also the variance of the tissue class distribution is allowed to vary over the image. Both of these methods rely on the standard K-means clustering for their initialization, and as local optimization methods, their performance is dependent on this initialization. In contrast to our method, these methods also assume a single, global MRF prior and do not account for PVE. We will also provide experimental results clearly demonstrating the superiority of the proposed approach to the Adaptive MAP algorithm [37
The organization of the paper is as follows. Section 2 describes the combination of the local MRFs defined on image sub-domains into a single global MRF-based image model. Section 3 describes the segmentation method rooted in the image model of Section 2. In Section 4, we present experimental validation of the method by comparing it quantitatively to other recent tissue classification schemes as well as the corresponding tissue classification method based on a traditional MRF-based image model. This validation is performed using both synthetic and real MR brain images. Section 5 summarizes the approach and discusses it in relation to the other approaches.