Brain Tissue Segmentation and Bias Field Correction of MR Image Based on Spatially Coherent FCM with Nonlocal Constraints

Influenced by poor radio frequency field uniformity and gradient-driven eddy currents, intensity inhomogeneity (or bias field) and noise appear in brain magnetic resonance (MR) image. However, some traditional fuzzy c-means clustering algorithms with local spatial constraints often cannot obtain satisfactory segmentation performance. Therefore, an objective function based on spatial coherence for brain MR image segmentation and intensity inhomogeneity correction simultaneously is constructed in this paper. First, a novel similarity measure including local neighboring information is designed to improve the separability of MR data in Gaussian kernel mapping space without image smoothing, and the similarity measure incorporates the spatial distance and grayscale difference between cluster centroid and its neighborhood pixels. Second, the objective function with an adaptive nonlocal spatial regularization term is drawn upon to compensate the drawback of the local spatial information. Meanwhile, bias field information is also embedded into the similarity measure of clustering algorithm. From the comparison between the proposed algorithm and the state-of-the-art methods, our model is more robust to noise in the brain magnetic resonance image, and the bias field is also effectively estimated.


Introduction
Magnetic resonance image has been widely used in diagnostic imaging for general check-up in clinical application, especially the detection and diagnosis of brain diseases. e volume change for brain tissues often indicates various diseases [1], such as brain tumor, leukoencephalopathy, olivopontocerebellar atrophy (OPCA), etc. erefore, brain tissue segmentation of MR image has become a very important medical treatment step. However, brain MR image has some lacks such as noise, intensity inhomogeneity, low contrast, the partial volume effect, and so on, which brings serious obstacles to segment the brain MR images. To this end, the multitudinous brain MR image segmentation methods have been put forward by using the theory of fuzzy set, random field, and level set.
Currently, there are two popular methods-based models for medical image segmentation: the random field theory [2][3][4] and the fuzzy c-means (FCM) algorithm. Random field is density-based unsupervised method where it finds the maximum likelihood estimate of the parameters from a given dataset. However random field algorithm has the disadvantages in high complexity and slow convergence and will drop into local optimization. FCM clustering is another efficient method used in image segmentation because it has robust characteristics for ambiguity and can retain much more information than random field algorithm [5].
e neighboring pixels in an image are highly correlated, i.e., the pixels in the immediate neighborhood possess nearly the same feature data. erefore, the spatial relationship of neighboring pixels is an important characteristic that can be of great aid in imaging segmentation. However, the conventional FCM algorithm does not fully utilize this spatial information. Pedrycz and Waletzky [9] took advantage of the available classified information and actively applied it as part of their optimization procedures. Szilagyi et al. [10] proposed the enhanced FCM (EnFCM) algorithm to accelerate the image segmentation process in which the pixels of an image are replaced the gray-level histogram and the statistical number and calculation are much smaller than FCM. In order to further reduce the computation time and improve the parameter inflexibility, Cai et al. [11] presented a fast generalized FCM (FGFCM) method, and FGFCM introduced a flexible locality factor S ij incorporating simultaneously both the gray-level difference and spatial distance in a local window. Ji et al. [12] proposed a robust spatially constrained fuzzy c-means (RSCFCM) algorithm for brain MR image segmentation. First, a spatial factor is constructed based on the posterior probabilities and prior probabilities and takes the spatial direction into account. Second, the negative log-posterior is utilized as dissimilarity function by taking the prior probabilities into account.
FCM with spatial constraint and its variants greatly improved the antinoise performance compared with FCM, but when the noise is very serious in the image, the performance of the algorithm may be worse. erefore, the nonlocal spatial information was often used and incorporated into the distance metric of FCM in recent years [13][14][15][16]. Zhao [14] brought in a nonlocal adaptive regularization term in its energy function, and the control factor is adaptive determined to adjust the balance of the objective function. Feng et al. [15] proposed a FCM method with specific nonlocal information for the segmentation of synthetic aperture radar (SAR) image. Ma et al. [16] proposed a modified FGFCM approach by introducing nonlocal constraint term, and local distance metric and nonlocal distance metric are used, respectively, in its objective function. By introducing nonlocal constraint term, the features of image can be used more comprehensively and effectively, and the robustness to noise of FCM-based algorithm is significantly improved. However, there generally exists intensity inhomogeneity in brain MR images. erefore, it is necessary to further design relevant algorithms to correct the intensity inhomogeneity. Sled et al. designed a set of software package for the estimation of bias field [17], and the characteristic of the method is nonparametric nonuniform intensity normalization or N3 for short; the distribution of the true tissue intensities can be achieved by an iterative method. Tustison et al. [18] improved the N3 algorithm based on modified B-spline approximation and hierarchical optimization algorithm (called N4ITK); N4ITK can also automatically perform without the priori knowledge. Liew and Hong Yan [19] introduced a spatial constraint to a fuzzy cluster method where the inhomogeneity field was modeled by a B-spline surface. e spatial pixel connectivity was implemented by a dissimilarity index, which enforced the connectivity constraint only in the homogeneous areas. Ji et al. proposed the modified possibilistic FCM (MPFCM) algorithm for bias field [20], generalized rough fuzzy c-means (GRFCM) algorithm, [21] and fuzzy local Gaussian mixture model (FLGMM) for brain MR image segmentation [22], respectively.
ose methods can estimate bias field and segment the MR images simultaneously.
In this paper, a brain tissue classification and intensity inhomogeneity correction model of MR image based on spatially coherent FCM with nonlocal constraints is proposed. In this model, firstly, both the local constraint term and nonlocal regularization term about brain MR image are incorporated into the objective function, and an adaptive control factor is used to maintain the balance between them. Secondly, the similarity measure is designed in Gaussian kernel mapping space without image filtering, and the detail information and the edge of the image can be preserved well. Meanwhile, bias field model is also embedded into the objective function of clustering algorithm. erefore, after the intensity inhomogeneity of the MR image is corrected, the segmentation accuracy is improved significantly. Experiments demonstrate that this algorithm can not only effectively estimate the bias field of MR image but also has stronger antinoise ability.

Preliminary Theory
2.1. Fuzzy Clustering Algorithms. Let X � x 1 , x 2 , . . . , x N denote an image with N pixels, where x k is gray value of the kth pixel in the image. FCM clustering aims at partitioning X into c clusters by minimizing the following objective function: where v i denotes the ith cluster prototype, u ik denotes the membership degree of x k belonging to ith cluster and follows c i�1 u ik � 1, c denotes the number of centers, ‖ · ‖ denotes the Euclidean norm, and the parameter m is a weight exponent on each fuzzy membership that determines the amount of fuzziness of the resulting partition.
Ahmed et al. proposed a modification to FCM objective function by introducing a term that allows the labeling of a pixel to be influenced by the labels in its immediate neighborhood [23]. is effect acts as a regularizer and biases the solution toward piecewise homogeneous labeling. It proved useful in segmenting images corrupted by noise. e modified objective function is given by where x r represents the neighbor voxels of x k and N R and N k stand for the number of voxels in the neighborhood of the kth voxel. e parameter α controls the intensity of the neighboring effect. One disadvantage of BCFCM is that the neighborhood labeling is computed in each iteration step, something that is very time consuming.

Spatially Coherent Fuzzy c-Means Clustering (SCFCM).
In view of some drawbacks of standard FCM algorithm, a modified scheme is proposed by Despotović et al. [24]. e similarity measure D ik � ‖x k − v i ‖ 2 is replaced by D * ik � (1 − S ik )‖x k − v i ‖ 2 introducing a weight factor S ik , and the objective function is where S ik is a weight factor including both local spatial information and grayscale difference and is designed as follows: where Ω k denotes a local neighboring window around x k , u ir denotes the membership degree of neighborhood pixels belonging to the ith cluster, a kr � |x k − x r | is the absolute intensity difference between the study pixel and its neighbor, is the Manhattan distance between the coordinates of pixel x k and x r , and (p k , q k ) and (p r , q r ) denote the coordinates x k and x r in the image, respectively. By minimizing equation (3) by Lagrangian multiplier approach, u ik and v i can be derived as shown in the following equation: Compared with the FCM, this algorithm has two advantages: firstly, it enhances the robustness to all kinds of noise. e constraint item of neighborhood information is included into the similarity measure, so as to effectively utilize the local information of the image. Secondly, it considers anisotropic neighborhood information, and more details and edges information can be preserved. However, the influence of bias field for MR images in segmentation algorithm is not mentioned.

Coherent Local Intensity Clustering Model.
Bias field of the MR image usually embodies slowly and smoothly varying for the pixel grayscale of the local region across an image. Meanwhile, in a neighboring local window of the image, bias field can be approximately considered as a constant.
erefore, the most popular model can be described in equation (6) [25]; let Y denote the observed image, b denote the unknown bias field, X denote the true image to be restored, and n denote the additive noise.
In the observed image, noise n is often assumed to obey Gaussian distribution with zero mean and variance σ 2 n , and the gray level value of the true image approximately takes a constant in a local window. erefore, the gray level of the observed image can be approximated to obey Gaussian distribution with mean bX and variance σ 2 n . In coherent local intensity clustering (CLIC) model [26], a novel metric introducing spatially coherent local intensity convergence criterion for bias field estimation and image segmentation simultaneously is proposed. A Gaussian kernel weight parameter K(r − k) is introduced into the similarity measure of each pixel gray level x k and its neighbor pixels x r , and the objective function of CLIC is where b r v i is the clustering prototype with bias field in the selective region Ω k , K(r − k) is the weight of a truncated Gaussian kernel allocated for the intensity x k , and the weight parameter can be defined as where σ denotes the standard deviation, a denotes a normalization factor to standardize Gaussian kernel, and ρ denotes a radius to measure the size of the local region. In the CLIC model, intensity inhomogeneity of the MR image can be effectively corrected and can reduce the misclassification rate, but there are some drawbacks in CLIC. When computing the distance metric between the central pixel and its surrounding pixels in a local region, the model only used the local neighborhood information of the pixel without considering the global structure information of the entire image. As a result, the antinoise performance of this algorithm is unsatisfactory.

Proposed Method
e standard FCM algorithm has the shortcoming of being sensitive to noise. ough, the modified FCM algorithms are improved by adding the spatial information, it is difficult to get a satisfied segmentation result for noise robustness. erefore, an improved FCM approach based on CLIC and SCFCM is proposed; its objective function is constructed according to the local constraint term and global regularization term; the similarity measure including local neighboring information is designed in Gaussian kernel mapping space, and brain tissue classification and intensity inhomogeneity correction can be realized simultaneously.

Nonlocal Weighted Constraint.
In a discrete noisy image X � x 1 , x 2 , . . . , x N , x k is the kth pixel and y k is its nonlocal weighted average value. e derivation method of the nonlocal weighted average can be acquired in [27], and the mathematical expression of y k is where w u k indicates a search region of radius u around x k , w kl denotes nonlocal weight coefficient depending on similarity measure between x k and its neighboring pixels x l in window w u k , and w kl satisfies the constraint conditions 0 ≤ w kl ≤ 1 and l w kl � 1. e weight coefficient w kl is computed as follows: where x(N k ) and x(N l ) denote the grayscale vectors in the square neighborhood N k and N l of radius s around x k and Computational and Mathematical Methods in Medicine x l , respectively, and ‖x(N k ) − x(N l )‖ 2 2,σ is the weighted Euclidean distance between x(N k ) and x(N l ); its expression is defined in equation (11). σ is the same as in equation (8), and h denotes a control factor to adjust the variation of the similarity measure w kl .
where x p (N k ) is the pth pixel in the grayscale vectors x(N k ) and σ p is defined as follows: where y p � mod(p, (2s + 1)) and z p � floor(p, (2s + 1))+ 1, (y p , z p ) denote the coordinates of the pth component in a preselected region and d � max(|y p − s − 1|, |z p − s − 1|).

Objective Function.
In order to correct bias field and classify the brain tissues simultaneously, the modified objective function-incorporated local constraint term and nonlocal regularization term is as follows: where u ik is the membership degree of x k belonging to the ith cluster, Ω k denotes a local square region of the radius s around the center x k , b k v i is the ith cluster center in Ω k , y k denotes the nonlocal weighted average value of x k , S ik denotes the weight factor of local neighborhood information in equation (4), and the definitions of K(r − k) and b k are the same as equation (7). α k is a trade-off weight factor to adjust the balance of local neighborhood information and nonlocal constraints information, and the definition of parameter α k is where x denotes the gray level mean of all pixels in the local region Ω k , and var(x) denotes the variance of pixel gray values in the same window. e larger the α k value is, the smaller the influence of the noise is. e factor α k can be obtained adaptively with the change of the local window Ω k without being given in advance.
Theorem. Assume c i�1 u ik � 1, 0 ≤ u ik ≤ 1, and m > 1. On the basis of Lagrange multiplier approach, equation (13) is minimized with respect to u ik , v i and b k can be derived as shown in the following equation: Proof. According the method of Lagrange multiplier, equation (13) can be converted to unconstrained optimization problem: where λ k is the Lagrange multiplier of the constraint condition c i�1 u ik � 1, by computing the partial derivatives of polynomial L in regard to u ik and λ k , respectively, and set zL/zu ik � 0, zL/zλ k � 0, as shown in the following equation: e following equation can be obtained by mathematical derivation of equation (17): Substituting equation (19) into equation (17), we obtain the following equation: .

(20)
And then substituting equation (20) into equation (19), the following equation can be obtained: Similarly, set zL/zv i � 0, that is e following equation can be obtained from equation (22) by mathematical derivation: We adopt the same mathematical derivation process to estimate bias field b k , for fixed u ik and λ k , and computing the partial derivative of L with respect to b k , set zL/zb k � 0, that is b k can be obtained from equation (24).
e theorem proves to be completed. Finally, the framework of the proposed algorithm can be summarized in Table 1.

Experimental Results and Analysis
In this section, several classical algorithms for intensity inhomogeneity correction and brain image segmentation are selected as the reference for comparison; bias field estimation and antinoise performance analysis for the brain MR images are the main experimental contents. For the experiments in the following sections, the related parameter values are fuzzy exponential m � 2, the stop criterion ε � 0.001, 3 × 3 neighborhood window, and the radius u � 10 of the search window.

MR Image Database.
Intensity inhomogeneity is one of the problems in interfering brain MR image segmentation; in the experiments of bias field correction, the dataset is from a simulated brain database (SBD) : BrainWeb [28] in which the brain MR images have three types: T1-, T2-, and proton density-(PD-) weighted 3D data volumes. In Figure 1, the T1-weighted normal brain MR images with 181 × 217 × 181 cubic voxels, 1 mm slice thickness, 40% intensity nonuniformity, and 3% noise are used to test; all the skulls and blood vessels are already stripped ahead of image processing, and the image is segmented into four regions: white matter (WM), gray matter (GM), cerebrospinal fluid (CSF), and background. Figure 1 shows the results of bias field correction and segmentation for three brain MR images. Figure 1(a) shows the brain slice images from three different directions: transaxial mode, sagittal mode, and coronal mode. Figure 1(b) shows the estimated bias field, Figure 1(c) shows the final segmentation results, and Figure 1(d) shows the corrected images after removed bias field. Figure 2 shows the histogram comparison of original MR image and bias corrected images corresponding to three images in Figure 2. From Figures 1 and 2, three brain tissues are more homogeneous after bias field correction; each brain tissue approximately obeys Gaussian distribution with different mean and variance, and WM, GM, and CSF can be clearly distinguished. In addition, the histogram distribution of corrected MR image is more reasonable, from which we can see three approximative peaks representing three brain tissues.

Computational and Mathematical Methods in Medicine
To further validate the performance of bias field correction, three bias correction algorithms including BCFCM [23], CLIC [26], and N4ITK [18] are chosen as comparative methods, as shown in Figure 3. Figure 3(a) is a T1-weighted transaxial slice of normal brain MR image with 40% spatial inhomogeneity; its slice thickness is 1 mm and the noise level is about 2%. Figures 3(b) and 3(c) are the obtained bias inhomogeneity and the amendatory MR images by BCFCM, CLIC, N4ITK, and our method, respectively. Figure (4) shows the histograms of the image with spatial inhomogeneity and the corrected images in Figure 3. It is can be seen that the distribution of pixel gray level of blue line is more accurate than red dotted line; it indicates that all the algorithms of bias field estimation are more or less effective. However, our method is more reasonable than other three algorithms from Figure 4. Because the histogram normally should have three peaks corresponding to WM, GM, and CSF; the peak value of CSF is minimal according to tissue volume, and the gray level's mean and variance of each tissue are also obviously different.

Antinoise Performance.
In the third experiment, first of all, the 104th transaxial slice of simulated T1-weighted brain MR image with 1 mm slice thickness and 7% Gaussian white noise is used to analyze the robustness to the noise, and we Table 1: e basic flow of image segmentation and bias field estimation.
Step 1. Set the number of cluster c, the exponent of fuzziness m, stopping condition of the iteration ε > 0, the radius s of local neighborhood window, and the radius u of search window in equation (11).
Step 2. Set the iteration initial value t � 0 and initialize bias field b 0 � 1 and the center of cluster Step 3. Calculate nonlocal weight coefficient w t kl in equation (8) and then obtain the nonlocal weighted value y t k by equation (7). Step 4. Update the membership degree u t ik by equation (14).
Step 5. Update the center of clustering v t i by equation (14). Step 6. Calculate the estimated bias field b t k by equation (15).
Step 7. If satisfying max‖V t+1 − V t ‖ < ε, then terminate iteration; otherwise, go to step 3 and set t � t + 1. select four algorithms: standard FCM, BCFCM [23], CLIC [26], and SCFCM [24] as the compared algorithms. e segmentation results are exhibited in Figure 5; Figure 5 because many pixels are misclassified; the segmentation result by SCFCM is better than FCM, BCFCM, and CLIC; however, there are still some noisy points that need to be removed. It can be observed that the presented method has more superior segmentation effect than four classical algorithms and can effectively eliminate the influence of the noise. Furthermore, the ability of detail and edge preservation is also compared and analyzed for five algorithms; we select a local region of the MR image to observe by enlarging 3 times, and the detailed images are presented in Figure 6; it is clearly seen that Figure 6(f ) is most similar subimage with the ground truth in Figure 6(g), and the vast majority of image details and edge are completely preserved.
In order to further evaluate and compare the antinoise ability of aforementioned five fuzzy clustering algorithms, we choose a brain slice image with 14% additive Gaussian white noise as the test object, as shown in Figure 7. Computational and Mathematical Methods in Medicine 9 [29] is adopted for comparison and quantitative analysis on the different level noise, JS is given as where S 1 represents a set of pixels of the segmented region by a clustering algorithm, S 2 denotes the set of pixels of the corresponding region acquired from the ground truth, ∩ denotes the intersection operation, and ∪ denotes the union operation. As a quantitative evaluation index, the values of JS belong to the interval of [0, 1], and the higher the JS value is, the better the segmentation performance is. We selected 15 noisy brain MR images with 20% intensity nonuniformity as the experimental objects and the noise level from 5% to 30%. ese images are segmented three regions: WM, GM, and CSF by FCM, BCFCM, CLIC, SCFCM, and our method, respectively; JS values comparison results are shown in Figure 8(a)-8(c). It is clearly shown that the presented approach has better matching degree with the ground truth and higher accuracy rate than other four clustering methods. en, we evaluate the effect of the search window radius u and the square neighborhood radius s on the performance of the proposed method. Here, we test u and s on the sets {4, 6, 8, 10, 12} and {1, 3, 5, 7}, respectively. In this experiment, the tested images perform 8 independent runs of the algorithm under each pair (u, s), and the noise level is 1% and 3%, respectively. Under each s value, the average JS curve of the algorithm with the increase of u value is shown in Figure 9. It can be found from Figure 9 that the algorithm under s � 3 and u � 10 can obtain satisfying performance on the noisy images.
In the aforementioned sections, the model is applied in the synthetic brain MR images. Next, this model is also applied to the real clinical images with noise. We selected three normal MR slice images from transaxial, coronal, and sagittal views, and these MR images are obtained from the Whole Brain Atlas clinical MR image database by the Harvard medical school [30]. Figure 10(a) shows three 2D T1-weighted brain MR slice images; the left image is a transaxial slice, the right image is a coronal slice, and the middle image is a sagittal slice. e segmentation results of brain slice images are given in Figure 10(b)-10(e) by BCFCM, HMRF-EM, SCFCM, and our method, respectively. From the experimental results, it is obvious that the proposed method can effectively segment each brain tissues as well as preserve more detail information of the original MR image. Furthermore, the experimental results of brain tissues in real MR images also further prove the robustness to noise of the proposed method.

Conclusion
Brain MR imaging has wide clinical application as an effective medical imaging diagnostic technique; however, the real brain MR images often suffer from some interference such as noise, intensity inhomogeneity, and low contrast. erefore, a brain tissue classification and nonuniformity field correction scheme in MR images based on spatially coherent FCM with nonlocal constraints is proposed in our study. e available information including local adjacent constraint and nonlocal global information of brain MR image is fully used in our model, and the similarity measure is designed in Gaussian kernel mapping space. Furthermore, the algorithm corrects the bias field of the MR image and improves its antinoise performance. Several experiments on the simulated brain MR images and real brain MR images have demonstrated that the proposed model can effectually overcome the effects of the noise while estimating the bias field existing in brain MR images.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.