Localized FCM Clustering with Spatial Information for Medical Image Segmentation and Bias Field Estimation

This paper presents a novel fuzzy energy minimization method for simultaneous segmentation and bias field estimation of medical images. We first define an objective function based on a localized fuzzy c-means (FCM) clustering for the image intensities in a neighborhood around each point. Then, this objective function is integrated with respect to the neighborhood center over the entire image domain to formulate a global fuzzy energy, which depends on membership functions, a bias field that accounts for the intensity inhomogeneity, and the constants that approximate the true intensities of the corresponding tissues. Therefore, segmentation and bias field estimation are simultaneously achieved by minimizing the global fuzzy energy. Besides, to reduce the impact of noise, the proposed algorithm incorporates spatial information into the membership function using the spatial function which is the summation of the membership functions in the neighborhood of each pixel under consideration. Experimental results on synthetic and real images are given to demonstrate the desirable performance of the proposed algorithm.


Introduction
Medical image segmentation plays an important role in a variety of biomedical-imaging applications, such as the quantification of tissue volumes, diagnosis, localization of pathology, study of anatomical structure, treatment planning, and computer-integrated surgery [1]. However, segmentation of medical images involves three main image-related problems [2]. First, images contain noise that can alter the intensity of a pixel such that its classification becomes uncertain. Second, images exhibit intensity inhomogeneity where the intensity level of a single tissue class varies gradually over the extent of the image. Third, images have finite pixel size and are subject to partial volume averaging where individual pixel volumes contain a mixture of tissue classes so that the intensity of a pixel in the image may not be consistent with any one class. To overcome these problems, many segmentation techniques have been proposed in the past decades, such as the expectation maximization (EM) algorithm [3][4][5], level set method [6][7][8][9], clustering [10][11][12][13][14][15][16][17], and so on.
Clustering for image segmentation usually classifies image pixels into -clusters such that members of the same cluster are more similar to one another than to members of other clusters, where the number, , of clusters is usually predefined or set by some validity criterion or a priori knowledge [18]. In the clustering methods, fuzzy -means (FCM) based algorithms have been widely used in medical image segmentation. Such a success chiefly attributes to the introduction of fuzziness for the belongingness of each image pixel. This enables the clustering methods to retain more information from the original image than the crisp or hard segmentation methods [10].
Pham and Prince proposed an adaptive FCM algorithm [10] and its extension to 3D data [11], which incorporated a spatial penalty term into the objective function to enable the estimated membership functions to be spatially smoothed. Ahmed et al. [12] 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 [13] used a B-spline surface to model the bias field and incorporated the spatial continuity constraints into fuzzy clustering algorithm. Zhang and Chen [14] 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. [15] proposed a modified FCM algorithm which was less sensitive to noise and yielded more homogeneous segmented regions. L. Szilágri et al. [16] proposed an efficient FCM clustering model for compensating intensity inhomogeneity and segmentation of magnetic resonance (MR) images, which drastically reduced the processing time without causing relevant change in terms of accuracy. Recently, local intensity information has been taken into account to deal with intensity inhomogeneity in fuzzy segmentation method. For example, Li et al. [17] proposed a new fuzzy energy minimization method based on coherent local intensity clustering (CLIC) for simultaneous tissue classification and bias field estimation of MR images. CLIC algorithm draws upon intensity information in local regions; therefore, it can be used to segment images with intensity inhomogeneity. However, spatial information is not taken into account in the CLIC algorithm; as a result, the CLIC algorithm is sensitive to noise.
Our proposed algorithm in this paper is motivated by the localized -means clustering model proposed by Chen et al. in [6]. By introducing the fuzzy belongingness of each pixel into Chen's model, we develop a localized FCM algorithm for image segmentation. We define a fuzzy energy that depends on membership functions, a bias field that accounts for the intensity inhomogeneity, and the constants that approximate the true intensities of the corresponding tissues. Hence, image segmentation and bias field estimation are simultaneously achieved as the result of minimizing this energy. Besides, we incorporate spatial information into the membership function to suppress noise. As an important application, our proposed algorithm can effectively segment medical images with intensity inhomogeneity and noise.
The remainder of this paper is organized as follows. Section 2 reviews a relevant method. In Section 3, we represent the definition and minimization of the proposed fuzzy energy in detail. We describe how to utilize neighborhood spatial information in Section 4. The algorithm 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.

Background
Chen et al. [6] applied a localized -means clustering to form an objective function modeling the problem of segmentation and bias field estimation for brain MR images. The image model of intensity inhomogeneity they used is defined as where is the measured image intensity, is the true image to be restored, and is an unknown bias field. Let̃,̃, and represent log , log , and log , respectively, then (1) can be rewritten as̃=̃+̃. (2) A generally accepted assumption on the bias field̃is that it is smooth or slowly varying [19]. Ideally, the intensitỹ belonging to the th tissue should take a specific value , which represents the measured physical property [6,8,17].
Chen's method is based on an observation that pixel intensities in a relatively small region are separable. Let = { : | − | ≤ } denote a circular neighborhood with a relatively small radius centered on each point in the image domain Ω. The partition {Ω } =1 ( is the total number of segmented regions) of the entire domain Ω induces a partition of the neighborhood ; that is, { ∩ Ω } =1 forms a partition of . For example, Figure 1 presents an image consisting of three disjoint regions: Ω 1 , Ω 2 , and Ω 3 , which divide the neighborhood into three subregions: ∩ Ω 1 , ∩ Ω 2 , and ∩ Ω 3 . Chen et al. defined an objective function to classify the datã( ) in the neighborhood into clusters using a -means clustering method: wherẽ( ) is the value of bias field̃at the center of , which is approximately equal to the valuẽ( ) for all ∈ on account of the smoothness of the bias field [6]; that is Thus, (̃( ) + ) ( = 1, . . . , ) are considered as the approximations of the cluster centers within the neighborhood , and ( − ) is a nonnegative weighting function such that ( − ) = 0 for | − | > and ∫ ( − ) = 1. Note that for each point , ( − ) has the nonzero value with respect to only in ∈ . Therefore, (3) can be rewritten as The ultimate goal is to find an optimal set of partitions for the entire image domain Ω, the bias field̃, and the constants . The minimization of a single criterion for a point does not accomplish this goal. The method minimizes for all ∈ Ω. This can be achieved by minimizing the integral of over Ω. Therefore, the energy is written as The above energy is expressed in terms of the regions Ω 1 , . . . , Ω . It is difficult to derive a solution to the energy minimization problem from this expression of . Alternatively, we can use one or multiple level set functions to represent the disjoint regions Ω 1 , . . . , Ω as in [20]. Thus, this energy can be converted into an equivalent level set formulation, which can be solved by using well-established variational methods [21]. centered on . The image domain Ω is divided into three disjoint regions Ω 1 , Ω 2 , and Ω 3 , which partition the neighborhood into three subregions ∩Ω 1 , ∩ Ω 2 , and ∩ Ω 3 .

Localized FCM Clustering
Chen's method can be considered as a hard segmentation method in which each pixel is assigned to an exclusive cluster. However, it is more suitable for medical images that each pixel is given a membership degree of belonging to each cluster, due to the impact of intensity inhomogeneity and noise. In this paper, we introduce the fuzzy belongingness of each pixel into Chen's model and, thus, propose a localized FCM clustering algorithm to implement the task of segmentation and bias field estimation.

Energy Formulation.
Similar to Chen's method, we first consider a task of classifying the datã( ) in the neighborhood into clusters. If the -means clustering in Chen's method is replaced by the FCM clustering, then the objective function in (3) can be converted to the following expression: where > 1 is the fuzzy coefficient, ( ) (0 ≤ ( ) ≤ 1) is the membership function of pixel belonging to the region Ω , and ( − ) is the same nonnegative weight function as in (3). 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 datã( ) for closer to the center of the neighborhood . In this paper, the weighting function is chosen as a truncated Gaussian kernel where is the standard deviation of the Gaussian kernel and is a constant to normalize the Gaussian kernel. The above objective function can be rewritten as follows: The desired clustering on the entire image domain Ω should have a good local performance in terms of the above objective function for every neighborhood . Therefore, we need to minimize for all ∈ Ω like Chen's method [6]. This can be achieved by minimizing the integral of over Ω. As a result, we define the following energy for our proposed localized FCM clustering: 3.2. Energy Minimization. The above energy 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 three necessary but not sufficient conditions for to be at a local extremum. In this subsection, we will derive these three conditions. (10) is subject to the constraint ∑ =1 ( ) = 1. Thus, this constrained optimization will be solved using one Lagrange multiplier

Membership Functions Evaluation. The energy in
Taking the derivative of with respect to ( ) and setting the result to zero, we have, for > 1 Solving for * ( ), we have * ( ) = ( 4 International Journal of Biomedical Imaging Since ∑ =1 ( ) = 1 for all , we have Substituting (15) into (13), the zero-gradient condition for the membership functions can be rewritten as * ( ) = 1

Bias Field Estimation.
In a similar way, taking the derivative of with respect tõ( ) and setting the result to zero, we have Solving for̃ * ( ), we havẽ

Constants Updating.
Likewise, taking the derivative of with respect to and setting the result to zero, we have Solving for * , we have * =

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 [15] and define a spatial function as follows: where NB( ) represents a square window centered on pixel in the spatial domain. A 5 × 5 window was 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 neighborhoods belong to the same cluster. The spatial function is incorporated into membership function as follows: To demonstrate the effect of removing noise by exploiting spatial information, we use a 3 × 3 neighborhood centered on a pixel under consideration. Without loss of generality, we assume that the image domain is divided into two regions; that is, = 2. Suppose that the values of membership functions of all neighborhood pixels belonging to the first cluster are shown in Figure 2. The upper row corresponds to the original values of membership functions, while the lower row shows the new values of membership functions by using (22). If we set the threshold to 0.50 for defuzzification, then the left and right columns of Figure 2 show the variations of the membership functions of a noisy pixel and a noisefree pixel, respectively. Obviously, the value of membership function of the noisy pixel has a desired correction from 0.60 to 0.33, while the noise-free pixel still belongs to the second cluster with a larger membership function.
In general, the spatial functions simply fortify the original membership in a homogeneous region, and the clustering result remains unchanged. However, for a noisy pixel, (22) reduces the weighting of a noisy cluster by the labels of its neighborhood pixels. As a result, misclassified pixels from noisy regions or spurious blobs can be easily corrected.

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 functions incorporated with the spatial information are updated, and the resulting new membership functions 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 functions ( ), constants , and bias field̃( ).
Step 5. Computing the new membership functions incorporated with spatial information using (22).
Repeat Steps 2-5 till termination. The iteration is stopped when the maximum difference between constants 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 function is maximal.
In this section, we apply the proposed algorithm to both synthetic and clinical medical images to demonstrate its effectiveness. The parameters used in our algorithm are as follows: fuzzy coefficient = 2, standard deviation of the Gaussian kernel = 4, and neighborhood radius of the Gaussian kernel = 15. To demonstrate the robustness, the initializations of the variables ( ),̃( ), and for the experiments in this paper are all generated as random fields or random numbers.

Segmentation of Synthetic
Images. The first experiment is performed in three synthetic images, which are displayed in the first column of Figure 3. In the first image, there is strong noise in both object and background regions. The images in the middle and bottom rows are corrupted by noise and intensity inhomogeneity. 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. It is revealed from Figure 3 that the result gradually improves 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.

Segmentation of Clinical Medical Images.
In this subsection, we compare the proposed algorithm with the biascorrected FCM (BCFCM) algorithm [12], the sFCM algorithm [15], and the CLIC algorithm [17] for clinical medical images.
The images in the first column of Figure 4 are two Xray vessel images with noise and intensity inhomogeneity. It can be seen that the upper parts of the images appear brighter, while the lower parts are darker due to the intensity inhomogeneity. As a result, the intensity values of the background in the upper region may be larger than the ones of vessels in the lower region. This phenomenon can cause serious misclassification for those clustering algorithms based on global region. Both of the BCFCM algorithm and the sFCM algorithm are based on global region clustering and hence cannot overcome this problem. This can be observed from the segmentation results which contain some parts of background in the upper brighter region while losing some vessel profiles in the lower darker region. However, the aforementioned phenomenon will become unobvious in local region because intensity inhomogeneity is slowing varying. Therefore, the local region clustering-based algorithm, namely the CLIC algorithm and the proposed algorithm can handle intensity inhomogeneity to obtain the complete vessel profile. Nevertheless, the CLIC algorithm has no step to resist noise so that its results contain some spurious blobs, due to the impact of noise. By contrast, our proposed algorithm utilizes spatial information to suppress noise and thus achieves the desirable segmentation results. The two images shown in the last column of Figure 4 are the estimated bias fields obtained by the proposed algorithm.  We also apply the aforementioned four algorithms to 3T brain MR images. The original images are also corrupted by intensity inhomogeneity and noise, which makes the images brighter in the middle than in the other regions. 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 5. Obviously, the BCFCM algorithm and the sFCM algorithm misclassify plentiful WM into GM in the vicinity of the skull because the WM in such region has approximate intensity values with the GM owing to the impact of the intensity inhomogeneity. The segmentation results of the CLIC algorithm show again that it is capable of dealing with the intensity inhomogeneity but unable to suppress the noise. However, 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 5.

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 include 40% image intensity inhomogeneity and 3% noise. The original images and the segmentation results are shown in Figure 6. We adopt Jaccard similarity (JS) [19] as a measurement of the segmentation accuracy. The JS between two regions 1 and 2 is defined as the ratio between the areas of the intersection and the union of them, namely, JS( 1, 2) = | 1 ∩ 2|/| 1 ∪ 2|. We compute the JS between the segmented region 1 obtained by the algorithm and the corresponding region 2 given by the ground truth. The closer the JS value 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 6 and Table 1 that the segmentation results of our proposed algorithm are more accurate than the other three algorithms.

Discussion
The proposed algorithm suffers from manually setting of two parameters: the neighborhood radius and the standard deviation of the truncated Gaussian kernel. Note that the radius should be selected appropriately according to the degree of the intensity inhomogeneity. For more localized intensity inhomogeneity, the bias field̃varies faster, and therefore the approximation in (4) is valid only in a smaller neighborhood. In this case, a smaller should be used as the radius of the neighborhood, and the standard deviation should be also selected a smaller value correspondingly. Figure 7 shows the JS values of the segmentation results with different parameters selection. The original image is obtained from Brain Web. The upper figure shows the influence of the radius, while the lower figure shows the influence International Journal of Biomedical Imaging 7 Figure 6: Comparison of segmentation results on two simulated brain MR images. Column 1: original images. Column 2: the BCFCM algorithm. Column 3: the sFCM algorithm. Column 4: the CLIC algorithm. Column 5: the proposed algorithm. Column 6: ground truth.  of the standard deviation. The accuracy of segmentations increases with the increasing of and . When > 10 or > 3, the JS values of WM and GM increase slightly, while the time consumption would have a significant increase. Considering the segmentation accuracy and the time consumption

Conclusion
In this paper, we have proposed a localized FCM clustering algorithm for simultaneous segmentation and bias field estimation of medical images. The proposed algorithm defines a fuzzy energy that depends on the bias field, membership functions, and the constants that approximate the true signal from the corresponding tissues. Bias field estimation and image segmentation are simultaneously achieved by minimizing this energy. Besides, we also utilize the neighborhood spatial information to resist the noise interference. Moreover, the proposed algorithm is robust to initialization, thereby allowing fully automatic applications. Comparisons with other approaches demonstrate the superior performance of the proposed algorithm.