Abstract
This paper presents a novel fuzzy algorithm for segmentation of brain MR images and simultaneous estimation of intensity inhomogeneity. The proposed algorithm defines an objective function including a local fuzzy energy and a global fuzzy energy. Based on the assumption that the local image intensities belonging to each different tissue satisfy Gaussian distributions with different means, we derive the local fuzzy energy by utilizing maximum a posterior probability (MAP) and Bayes rule. The global fuzzy energy is defined by measuring the distance between the original image and the corresponding inhomogeneity-free image. We combine the global fuzzy energy with the local fuzzy energy using an adaptive weight function whose value varies with the local contrast of the image. This combination enables the proposed algorithm to address intensity inhomogeneity and to improve the accuracy of segmentation and its robustness to initialization. Besides, the proposed algorithm incorporates neighborhood spatial information into the membership function to reduce the impact of noise. Experimental results for synthetic and real images validate the desirable performances of the proposed algorithm.
1. Introduction
Magnetic resonance imaging (MRI), compared to other medical imaging modalities, has many significant advantages, including high contrast among different soft tissues, relatively high spatial resolution across the entire field of view, and multispectral characteristics [1]. Thus, it has been widely used in both clinical practice and neuroscience research. For example, quantitative analysis of brain tissues, that is, white matter (WM), gray matter (GM), and cerebral spinal fluid (CSF) from brain MR images provides very helpful information for the study and the treatment of various pathologies such as Alzheimer disease, Parkinson, or Parkinson-related syndrome. Therefore, accurate segmentation of brain MR images is an essential and vital step. However, due to nonuniform magnetic field or susceptibility effects, brain MR images are usually corrupted by intensity inhomogeneity which is also referred to as the bias field and manifests itself as a smooth intensity variation across the image [2]. As a result, the intensities of the same tissue vary with the locations of the tissue within the image. This may cause serious misclassification when intensity-based segmentation algorithms are used. Therefore, intensity inhomogeneity has been a challenging difficulty in the task of image segmentation.
Among many available segmentation algorithms, fuzzy segmentation methods, especially the fuzzy c-means (FCM)-based algorithms, have been broadly applied to brain MR image segmentation. Such a success chiefly attributes to the introduction of fuzziness for the belongingness of each image pixel [3]. This enables the clustering methods to retain more information from the original image than the crisp or hard segmentation methods [4].
Pham and Prince [4] proposed an adaptive FCM algorithm which incorporated a spatial penalty term into the objective function to enable the estimated membership functions to be spatially smoothed. Ahmed et al. [5] developed the bias-corrected FCM (BCFCM) algorithm which modified the objective function of the standard FCM algorithm to compensate for intensity inhomogeneity and to allow the labeling of a pixel to be influenced by the labels in its immediate neighborhood. Liew and Yan [6] used a B-spline surface to model the bias field and incorporated the spatial continuity constraints into fuzzy clustering algorithm. Zhang and Chen [7] replaced the original Euclidean distance with a kernel-induced distance and supplemented the objective function with a spatial penalty term, which modeled the spatial continuity compensation. Incorporating spatial information into the membership function, Chuang et al. [8] proposed a modified FCM algorithm which was less sensitive to noise and yielded more homogeneous segmented regions. Recently, local intensity information has been taken into account to deal with intensity inhomogeneity in fuzzy segmentation methods. For example, Li et al. [9] proposed a fuzzy energy minimization method based on coherent local intensity clustering (CLIC) for simultaneous tissue classification and bias field estimation of MR images. This CLIC algorithm draws upon intensity means in local regions; therefore, it can be used to segment images with intensity inhomogeneity. However, global information is not taken into account in this model; as a result, this model is not robust to the parameters and also apt to suffer from misclassification for the pixels on the boundary especially in the regions with heavy level of intensity inhomogeneity [1]. Ji et al. [10] assumed that the local image data within each pixel’s neighborhood satisfied the Gaussian mixture model (GMM), and thus proposed the fuzzy local GMM (FLGMM) algorithm. Since this algorithm exploits local intensity means and variances, it produces more accurate segmentation results than several state-of-the-art algorithms. Nevertheless, it has a nonconvex objective function, and hence could be easily trapped in a local minimum, resulting in an unsatisfactory segmentation [10]. In addition, for clustering-based image segmentation algorithm, spatial correlation is frequently used to reduce the impact of noise [3, 5–8]. However, both the CLIC algorithm and the FLGMM algorithm ignore the spatial information; thus, they are sensitive to noise to some extent.
In this paper, in order to handle intensity inhomogeneity, we assume that the local image intensities belonging to each different tissue satisfy Gaussian distributions with different means, and then apply maximum a posteriori probability (MAP) and Bayes rule to derive a local fuzzy energy function in which the bias field accounting for intensity inhomogeneity acts as a variable. For the purpose of improving the accuracy of segmentation and alleviating the problems induced by the nonconvexity of the proposed local fuzzy energy function, we introduce a global fuzzy energy which is defined by measuring the distance between the original image and the corresponding inhomogeneity-free image. Moreover, these two energy functions are combined using an adaptive weight function whose value varies with the local contrast of the image. Besides, the proposed algorithm incorporates neighborhood spatial information into the membership function to reduce the impact of noise. Therefore, our proposed algorithm can effectively segment brain MR images with severe intensity inhomogeneity and noise.
The remainder of this paper is organized as follows. Section 2 gives the image model for intensity inhomogeneity. In Section 3, we present the definition and minimization of the proposed fuzzy objective function in detail. We describe how to utilize neighborhood spatial information in Section 4. The implementation and experimental results are given in Section 5. The discussion on the setting of important parameters is given in Section 6. We end this paper by the conclusion in Section 7.
2. Image Model for Intensity Inhomogeneity
Intensity inhomogeneity can be modeled as a multiplicative component of an observed image [2], as shown in the following: where is the observed image, is an unknown bias field, is the inhomogeneity-free image or true image, and is the additive noise. A generally accepted assumption on the bias field is that it is smooth or slowly varying. Ideally, the intensity belonging to the th tissue should take a specific value , which represents the measured physical property. The additive noise is assumed to be Gaussian noise with zero mean and variance .
3. Our Proposed Objective Function
The proposed objective function includes a local fuzzy energy and a global fuzzy energy, which are formulated by exploiting local intensity information and global intensity information, respectively. In this section, we will give a detailed description with respect to the definition and minimization of the objective function.
3.1. Local Fuzzy Energy
To effectively exploit information on local intensities, we need to characterize the distribution of local intensities via partition of neighborhood as in [11]. For each point in the image domain , we consider a circular neighborhood with a small radius , which is defined as . The partition ( is the total number of segmented tissues) of the entire domain induces a partition of the neighborhood ; that is, forms a partition of . For a slowly varying bias field , the values for all in the circular neighborhood can be well approximated by its value at the center of ; that is,
Together with the aforementioned assumption with respect to the true image and the Gaussian noise in (1), we can describe the intensities within subregion as a Gaussian distribution with mean and variance . In other words, the probability density in subregion is
Now, we consider the segmentation of the circular neighborhood based on the maximum a posteriori probability (MAP). Let be the a posteriori probability of the subregion given the neighborhood gray values . According to the Bayes rule, where is the probability density in region , which is denoted by in (3). is the a priori probability of the partition among all partitions of , which is usually set equal in all partitions [12]. The a priori probability is independent of the choice of the region.
Assuming that the pixels within each region are independent, the MAP can be achieved by finding the maximum of . Taking a logarithm, the maximization can be converted to the minimization of the following energy:
Now, we introduce the fuzzy membership function into the above energy. Since represents the probability of pixel belonging to region such that for , we can rewrite (5) as where is a weighting exponent on each fuzzy membership and determines the amount of fuzziness of the resulting classification.
Next, we use a weighting function to impose spatial constraint as in [9, 13], we define the following objective function: where is a nonnegative weighting function such that for and
Although the choice of the weighting function is flexible, it is preferable to use a weighting function such that larger weights are assigned to the data for closer to the center of the neighborhood . In this paper, the weighting function is chosen as a truncated Gaussian kernel as follows: where is the standard deviation of the Gaussian kernel and is a constant to normalize the Gaussian kernel. The energy function in (7) can be rewritten as as for .
The ultimate goal is to minimize for all the center points in the image domain , which directs us to define the following double integral energy functional as our proposed local fuzzy energy:
3.2. Global Fuzzy Energy
Although many algorithms, such as fuzzy segmentation method [10] and active contour models [11, 13–16] have utilized only the local intensity information to successfully segment images with intensity inhomogeneity, a common drawback for these algorithms is that they have a nonconvex energy function and thus are easy to result in local optima. To avoid unsatisfying segmentation results led by bad initial conditions, a usual settlement is to incorporate the global intensity information [17–22]. In this paper, we measure the distance between the original image and the corresponding inhomogeneity-free image as the global energy function to be minimized. Our proposed global fuzzy energy is defined as where is the constant value of intensity belonging to the th tissue in the inhomogeneity-free image. Because represents the probability that pixel belongs to the th tissue, the distance at location between the original image and the corresponding inhomogeneity-free image should be a weighted summation, that is, . Adding up the distance for all in the whole image domain , we can get the expression as presented in (11).
3.3. Definition of the Objective Function
We combine the local fuzzy energy with the global fuzzy energy , and thus obtain the proposed objective function as follows: where is the weight of the global fuzzy term.
Generally speaking, in regions far away from the desired boundary of objects (i.e., tissues) where the intensities vary slowly and appear approximately homogeneous, the probability density for different cluster in local region is almost the same, thus the estimated membership functions by minimizing the energy in (7) will be approximately equal to each other; that is, the segmentation result is very “fuzzy”. This will result in inaccuracy after defuzzification. However, the global fuzzy energy in (11) is very appropriate to segment piecewise homogeneous regions because it is similar to the standard FCM objective function. Hence, we should increase the weight of the global fuzzy term to reduce the “fuzzy” degree of the estimated membership functions. On the contrary, in regions close to the desired boundary of objects, the probability density can precisely reflect the difference of distributions for different tissues and may have a large deviation if we introduce the global fuzzy term redundantly. Thus, we need to reduce the weight of the global fuzzy term in such regions.
Based on the above discussion, we should choose a weight function for that varies dynamically with the locations of the image instead of a constant value. Therefore, we refer to [19–22] and define the weight function as follows: where is an adjustable parameter and denotes the local contrast ratio of the given image, which is defined as where defines the size of the local square window centered on . we set for all experiments in this paper. and are the maximum and minimum of intensities within this local window, respectively. represents the total intensity levels of the image, for example, = 255 for gray image.
In the above weight function (13), mean () is the average value of over the whole image and it can reflect the overall contrast information of the image. For an image with a higher overall contrast, the background and foreground may be more clearly classified on the whole; thus, a bigger weight proportional to mean () should be chosen for the global fuzzy term. can adjust the weight of global fuzzy term dynamically in all regions, making it smaller in regions with high local contrast and larger in regions with low local contrast.
3.4. Minimization of the Objective Function
Our proposed objective function in (12) can be minimized in a fashion similar to the standard FCM algorithm. Taking the first derivatives of with respect to , , and , and setting them to zero results in four necessary but not sufficient conditions for to be at a local extremum. In this subsection, we will derive these four conditions.
The energy is subject to the constraint . Thus, this constrained optimization will be solved using one Lagrange multiplier as follows:
Taking the derivative of with respect to and setting the result to zero, we have, for , Solving for , we have Since for all , we have or
Substituting (3) and (19) into (17), the zero-gradient condition for the membership functions can be rewritten as
In a similar way, taking the derivative of with respect to and setting the result to zero, we have Solving for , we have
Likewise, taking the derivative of with respect to and and setting the results to zero, we have Solving for and , we have
4. Exploiting Spatial Information
One of the important characteristics of an image is that neighborhood pixels are highly correlated. In other words, these neighborhood pixels possess similar intensity, and the probability that they belong to the same cluster is great. This spatial relationship is important in clustering, but it is not utilized in a conventional FCM algorithm. To exploit the spatial information, we refer to [8] and define a spatial function as follows: where represents a square window centered on pixel in the spatial domain. A 5 × 5 window is used throughout this work. Just like the membership function, the spatial function represents the probability that pixel belongs to the th cluster. The spatial function of a pixel for a cluster is large if the majority of its neighborhood pixels belong to the same cluster. The spatial function is incorporated into membership function as follows:
In a homogeneous region, the spatial functions simply fortify the original membership, and the clustering result remains unchanged. However, for a noisy pixel, this formula reduces the weight of a noisy cluster by the labels of its neighborhood pixels. As a result, misclassified pixels from noisy regions or spurious blobs can easily be corrected.
5. Implementation and Experimental Results
The proposed algorithm is a two-pass process at each iteration. The first pass is to calculate the corresponding variables , , , and . In the second pass, the membership function incorporated with the spatial information is updated and the resulting new membership function will be inserted into the next iteration. The detailed procedures can be summarized in the following steps.
Step 1. Initialize the number of clusters , membership function , constant intensity value , and bias field .
Step 2. Computing the weight function using (13).
Step 3. Updating , , , and sequentially using (26), (22), (20), and (25), respectively.
Step 4. Computing the new membership functions incorporated with spatial information using (28).
Repeat Steps 3 and 4 till termination. The iteration is stopped when the maximum difference between two constant intensity values at two successive iterations is less than a threshold (e.g., 0.001). After the convergence, defuzzification is applied to assign each pixel to a specific cluster for which the membership is maximal.
In this paper, membership function and constant intensity value are initialized randomly at the intervals and , respectively. Bias field is initially equal to 1 for every pixel . The parameters used in this paper are mostly set as the following default values: fuzzy exponent , standard deviation of the Gaussian kernel = 4, and neighborhood radius of the Gaussian kernel . Unless otherwise specified, the parameter is set to 0.005.
5.1. Segmentation of Synthetic Images
We first apply the proposed algorithm to three synthetic images to demonstrate its effectiveness. The original images shown in the first column of Figure 1 are corrupted by strong noise and intensity inhomogeneity. The parameters for these three images are set to 0.001. The intermediate segmentation results obtained by running the proposed algorithm for different numbers of iterations are shown in the second and third columns, and the final results obtained after the convergence of our algorithm are shown in the fourth column. The three images in the fifth column are the estimated bias fields. It is revealed from Figure 1 that the results gradually improve during the iterative segmentation process. In the final segmentation results, objects and background can be clearly differentiated despite of the impact of noise and intensity inhomogeneity.
5.2. Segmentation of Blood Vessel Images
In this subsection, we compare the proposed algorithm with the bias-corrected FCM (BCFCM) algorithm [5], the sFCM algorithm [8], and the CLIC algorithm [9] for blood vessel images.
The images in the first column of Figure 2 are two X-ray vessel images with noise and intensity inhomogeneity. In addition, some parts of the boundaries are quite weak. Our goal in this experiment is to segment the vessels from the inhomogeneous and noisy background. From the segmentation results shown in Figure 2, we can observe that both the BCFCM algorithm and the sFCM algorithm are based on the global intensity information and neither can address intensity inhomogeneity as well as the other two algorithms, that is, the CLIC algorithm and the proposed algorithm, which are based on local intensity information. Thus, they get rather unsatisfactory vessel segmentation results with some noisy spots, spurious blobs, and incomplete boundaries. By contrast with the CLIC algorithm, the proposed algorithm integrates the global and local intensity information and tries to exploit the neighborhood spatial information, and hence obtains desirable segmentation results with less spurious blobs. The two images shown in the last column of Figure 2 are the estimated bias fields obtained by the proposed algorithm.
5.3. Segmentation of Brain MR Images
We also apply the aforementioned four algorithms to 3T brain MR images. The original images are also corrupted by intensity inhomogeneity and noise. The task of segmentation is to partition the brain MR images into four regions, that is, white matter (WM), gray matter (GM), cerebral spinal fluid (CSF), and background. The comparison of segmentation results obtained by these four algorithms is shown in Figure 3. Obviously, the BCFCM algorithm and the sFCM algorithm misclassify plentiful WM into GM in the vicinity of the skull, while the segmentation results of CLIC algorithm include a lot of noisy spots over the images. By contrast, our proposed algorithm gets fairly better segmentation with clear and correct classification of tissues. The estimated bias fields obtained by our proposed algorithm are shown in the sixth column of Figure 3.
5.4. Quantitative Comparison
To quantitatively compare the proposed algorithm with the above-mentioned other three algorithms, we use the T1-weighted simulated brain MR images with ground truth from Brain Web in the link http://www.bic.mni.mcgill.ca/brainweb/. The selected MR images contain 40% image intensity inhomogeneity and 3% noise. We set the parameter = 0.008 for these images. The original images and the segmentation results are shown in Figure 4. We adopt Jaccard similarity (JS) [2] as a measurement of the segmentation accuracy. The JS between two regions and is defined as the ratio between the areas of the intersection and the union of them, namely, . We compute the JS between the segmented region obtained by the algorithm and the corresponding region given by the ground truth. The closer the JS value is to 1, the better the segmentation result. The resulting JS values for the four algorithms are listed in Table 1. It can be observed from both Figure 4 and Table 1 that the segmentation results of our proposed algorithm are more accurate than the other three algorithms.
Meanwhile, we also compare the computation time of the four algorithms for the images in Figure 4. The CPU times are listed in Table 2, which are recorded from our experiments with Matlab programs run on a Dell Vostro 1450 notebook with Intel Core i5 CPU, 2.5?GHz, 4?GB RAM, with Matlab R2010a on Windows 7. It is obvious that the BCFCM and sFCM algorithms have less time consumption, but they have lower segmentation accuracy. The CLIC algorithm and our proposed algorithm improve the segmentation accuracy at the cost of higher computational complexity. However, our algorithm accelerates the convergence of the iteration process by introducing the global fuzzy energy term; hence, it takes less CPU time than the CLIC algorithm.
6. Discussion
In this paper, the truncated Gaussian kernel in (8) is used to determine the size of the local regions in which the proposed local fuzzy energy is defined. Its two parameters, neighborhood radius , and standard deviation , can influence the segmentation accuracy to some extent. Figure 5 shows the JS values of the segmentation results with different parameter settings. The original image is obtained from Brain Web. The left figure shows the influence of the radius while the right figure shows the influence of the standard deviation . It is clear that when or , the JS values of WM, GM, and CSF remain stable on the whole. However, the time consumption would have a significant increase with the increase of the size of neighborhood radius. Considering the segmentation accuracy and the time consumption of the algorithm, we suggest that and . In our experiments, we set and = 4 for all test images.
As an adjustable variable of the weight function , the parameter plays an important role to control the weight of the global and local fuzzy energy in our proposed objective function. Its value is typically selected from 0.001 to 0.01. If intensity inhomogeneity in the given image is severe, a value approximate to 0.001 should be chosen so that the local fuzzy energy is dominant. Otherwise, a bigger value close to 0.01 should be chosen. We first illustrate the effect of this parameter for the images with less intensity inhomogeneity. The original images are shown in the first column of Figure 6. The segmentation results obtained by our proposed algorithm with different are displayed from the second to fourth column of Figure 6. It can be seen that = 0.009 is more appropriate for this type of images. In addition, Figure 7 shows the segmentation results for images with strong intensity inhomogeneity using different . It can be observed from the second to fourth column of Figure 7 that the segmentation results become more unsatisfactory with the increase of . Therefore, we should choose the smaller for images with severe intensity inhomogeneity.
7. Conclusion
In this paper, we have proposed a global and local fuzzy energy based algorithm for simultaneous segmentation and bias field estimation of images with different imaging modalities. The global fuzzy energy aims to minimize the distance between the original image and the corresponding inhomogeneity-free image. The local fuzzy energy is defined by using local intensity statistics and MAP and Bayes rule. These two fuzzy energy terms are combined using an adaptive weight function which depends on the local contrast of the image. This enables the proposed algorithm not only to address intensity inhomogeneity effectively but also to improve the accuracy of segmentation and its robustness to initialization. Moreover, we utilize neighborhood spatial information to update the membership function in order to reduce the impact of noise. Comparisons with other approaches demonstrate the superior performance of the proposed algorithm.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (60903127, 61202314), the NPU Foundation for Fundamental Research (JCY20130130), “Top Paper” Project (13GH013308), and Ao-Xiang Star Project (11GH0315) at the Northwestern Polytechnical University, Xi’an, Shannxi, China.