# Brain MR Image Segmentation Based on an Adaptive Combination of Global and Local Fuzzy Energy

Authors: Wenchao Cui; Yi Wang; Tao Lei; Yangyu Fan; Yan Feng
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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. Section2 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:
(1)I=bJ+n,
where I is the observed image, b is an unknown bias field, J is the inhomogeneity-free image or true image, and n is the additive noise. A generally accepted assumption on the bias field b is that it is smooth or slowly varying. Ideally, the intensity J belonging to the ith tissue should take a specific value ci, which represents the measured physical property. The additive noise n is assumed to be Gaussian noise with zero mean and variance σ2.


## 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 y in the image domain Ω, we consider a circular neighborhood with a small radius r, which is defined as Oy={x:|x-y|≤r}. The partition {Ωi}i=1N  (N is the total number of segmented tissues) of the entire domain Ω induces a partition of the neighborhood Oy; that is, {Oy∩Ωi}i=1N forms a partition of Oy. For a slowly varying bias field b, the values b(x) for all x in the circular neighborhood Oy can be well approximated by its value b(y) at the center of Oy; that is,
(2)b(x)≈b(y),x∈Oy.Together with the aforementioned assumption with respect to the true imageJ and the Gaussian noise n in (1), we can describe the intensities I(x) within subregion Oy∩Ωi as a Gaussian distribution with mean b(y)·ci and variance σ2. In other words, the probability density in subregion Oy∩Ωi is
(3)pi,y(I(x))=12πσexp(-(I(x)-b(y)ci)22σ2).Now, we consider the segmentation of the circular neighborhoodOy based on the maximum a posteriori probability (MAP). Let p(x∈Oy∩Ωi∣I(x)) be the a posteriori probability of the subregion Oy∩Ωi given the neighborhood gray values I(x). According to the Bayes rule,
(4)p(x∈Oy∩Ωi∣I(x))=p(I(x)∣x∈Oy∩Ωi)p(x∈Oy∩Ωi)p(I(x)),
where p(I(x)∣x∈Oy∩Ωi) is the probability density in region Oy∩Ωi, which is denoted by pi,y(I(x)) in (3). p(x∈Oy∩Ωi) is the a priori probability of the partition Oy∩Ωi among all partitions of Oy, which is usually set equal in all partitions [12]. The a priori probability p(I(x)) 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∏i=1N∏x∈Ωi∩Oypi,y(I(x)). Taking a logarithm, the maximization can be converted to the minimization of the following energy:
(5)Ey=∑i=1N∫Oy∩Ωi-logpi,y(I(x))dx.Now, we introduce the fuzzy membership functionui(x) into the above energy. Since ui(x) represents the probability of pixel x belonging to region Ωi such that ui(x)=0 for x∉Ωi, we can rewrite (5) as
(6)Ey=∑i=1N∫Oy-uim(x)logpi,y(I(x))dx,
where m 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:
(7)Ey=∑i=1N∫Oy-uim(x)ω(x-y)logpi,y(I(x))dx,
where ω(x-y) is a nonnegative weighting function such that ω(x-y)=0 for |x-y|>r and ∫Oyω(x-y)dx=1.Although the choice of the weighting function is flexible, it is preferable to use a weighting functionω(x-y) such that larger weights are assigned to the data I(x) for x closer to the center y of the neighborhood Oy. In this paper, the weighting function ω is chosen as a truncated Gaussian kernel as follows:
(8)ω(d)={1ae-|d|2/2τ2,|d|≤r,0,else,
where τ is the standard deviation of the Gaussian kernel and a is a constant to normalize the Gaussian kernel. The energy function Ey in (7) can be rewritten as
(9)Ey=∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx,
as ω(x-y)=0 for x∉Oy.The ultimate goal is to minimizeEy for all the center points y in the image domain Ω, which directs us to define the following double integral energy functional as our proposed local fuzzy energy:
(10)EL=∫ΩEydy=∫Ω(∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx)dy.


### 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
(11)EG=∫Ω(∑i=1Nuim(x)|I(x)-ci|2)dx,
where ci is the constant value of intensity J belonging to the ith tissue in the inhomogeneity-free image. Because ui(x)(0≤ui(x)≤1) represents the probability that pixel I(x) belongs to the ith tissue, the distance at location x between the original image and the corresponding inhomogeneity-free image should be a weighted summation, that is, d(x)=∑i=1Nuim(x)|I(x)-ci|2. Adding up the distance d(x) for all x 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 energyEL with the global fuzzy energy EG, and thus obtain the proposed objective function as follows:
(12)E=β∫Ω(∑i=1Nuim(x)|I(x)-ci|2)dx+∫Ω(∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx)dy,
where β(0≤β≤1) 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 densitypi,y(I(x))(x∈Oy∩Ωi) for different cluster i in local region Oy is almost the same, thus the estimated membership functions by minimizing the energy Ey 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 pi,y(I(x)) 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:
(13)β(x)=γ·mean(LCRW(x))·(1-LCRW(x))x∈Ω,
where γ is an adjustable parameter and LCRW denotes the local contrast ratio of the given image, which is defined as
(14)LCRW(x)=(ILmax-ILmin)Ig,
where W defines the size of the local square window centered on x. we set W=5 for all experiments in this paper. ILmax and ILmin are the maximum and minimum of intensities within this local window, respectively. Ig represents the total intensity levels of the image, for example, Ig = 255 for gray image.In the above weight function (13), mean (LCRW(x)) is the average value of LCRW(x) 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 (LCRW(x)) should be chosen for the global fuzzy term. (1-LCRW(x)) 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 functionE in (12) can be minimized in a fashion similar to the standard FCM algorithm. Taking the first derivatives of E with respect to ui(x), b(y), ci and σ, and setting them to zero results in four necessary but not sufficient conditions for E to be at a local extremum. In this subsection, we will derive these four conditions.The energyE is subject to the constraint ∑i=1Nui(x)=1. Thus, this constrained optimization will be solved using one Lagrange multiplier as follows:
(15)F=β∫Ω(∑i=1Nuim(x)|I(x)-ci|2)dx+∫Ω(∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx)dy+λ(1-∑i=1Nui(x)).Taking the derivative ofF with respect to ui(x) and setting the result to zero, we have, for m>1,
(16)[∂F∂ui(x)=βmuim-1(x)|I(x)-ci|2+∫Ω-muim-1(x)ω(x-y)logpi,y(I(x))dy-λ∂F∂ui(x)]ui(x)=ui*(x)=0.
Solving for ui*(x), we have
(17)ui*(x)=(λβm|I(x)-ci|2+∫Ω-mω(x-y)logpi,y(I(x))dy)1/(m-1).
Since ∑K=1Nuk(x)=1 for all x, we have
(18)∑k=1N(λβm|I(x)-ck|2+∫Ω-mω(x-y)logpk,y(I(x))dy)1/(m-1)=1,
or
(19)λ=m×((∫Ω-ω(x-y)logpk,y(I(x))dy∑k=1N({+∫Ω-ω(x-y)logpk,y(I(x))dy)-1}(∫Ω-ω(x-y)logpk,y(I(x))dyβ|I(x)-ck|2+∫Ω-ω(x-y)logpk,y(I(x))dy)-1)1/(m-1))m-1)-1.Substituting (3) and (19) into (17), the zero-gradient condition for the membership functions can be rewritten as(20)ui*(x)=1∑k=1N(β|I(x)-ci|2+∫Ωω(x-y)(log(2πσ)+(I(x)-b(y)ci)2/2σ2)dyβ|I(x)-ck|2+∫Ωω(x-y)(log(2πσ)+(I(x)-b(y)ck)2/2σ2)dy)1/(m-1).In a similar way, taking the derivative ofF with respect to b(y) and setting the result to zero, we have
(21)[∑i=1N∫Ωuim(x)ω(x-y)×(I(x)-b(y)ci)(-ci)dx∑i=1N∫Ωuim(x)ω(x-y)]b(y)=b*(y)=0.
Solving for b*(y), we have
(22)b*(y)=∑i=1N∫Ωuim(x)ω(x-y)I(x)cidx∑i=1N∫Ωuim(x)ω(x-y)ci2dxLikewise, taking the derivative ofF with respect to ci and σ and setting the results to zero, we have
(23)[2β∫Ωuim(x)(I(x)-ci)(-1)dx+∫∫Ωuim(x)ω(x-y)×(I(x)-b(y)ci)(-b(y))σ2dxdy∫Ωuim(x)(I(x)-ci)(-1)dx]ci=ci*=0,(24)[∑i=1N∫∫Ωuim(x)ω(x-y)×(σ2-(I(x)-b(y)ci)2)dxdy∑i=1N∫Ω∫Ωuim(x)ω(x-y)]σ=σ*=0.
Solving for ci and σ, we have(25)ci*=2β∫Ωuim(x)I(x)dx+∫∫Ωuim(x)ω(x-y)I(x)b(y)(1/σ2)dxdy2β∫Ωuim(x)dx+∫∫Ωuim(x)ω(x-y)b2(y)(1/σ2)dxdy,(26)σ*=∑i=1N∫∫Ωuim(x)ω(x-y)(I(x)-b(y)ci)2dxdy∑i=1N∫∫Ωuim(x)ω(x-y)dxdy.


## 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 y in the image domain Ω, we consider a circular neighborhood with a small radius r, which is defined as Oy={x:|x-y|≤r}. The partition {Ωi}i=1N  (N is the total number of segmented tissues) of the entire domain Ω induces a partition of the neighborhood Oy; that is, {Oy∩Ωi}i=1N forms a partition of Oy. For a slowly varying bias field b, the values b(x) for all x in the circular neighborhood Oy can be well approximated by its value b(y) at the center of Oy; that is,
(2)b(x)≈b(y),x∈Oy.Together with the aforementioned assumption with respect to the true imageJ and the Gaussian noise n in (1), we can describe the intensities I(x) within subregion Oy∩Ωi as a Gaussian distribution with mean b(y)·ci and variance σ2. In other words, the probability density in subregion Oy∩Ωi is
(3)pi,y(I(x))=12πσexp(-(I(x)-b(y)ci)22σ2).Now, we consider the segmentation of the circular neighborhoodOy based on the maximum a posteriori probability (MAP). Let p(x∈Oy∩Ωi∣I(x)) be the a posteriori probability of the subregion Oy∩Ωi given the neighborhood gray values I(x). According to the Bayes rule,
(4)p(x∈Oy∩Ωi∣I(x))=p(I(x)∣x∈Oy∩Ωi)p(x∈Oy∩Ωi)p(I(x)),
where p(I(x)∣x∈Oy∩Ωi) is the probability density in region Oy∩Ωi, which is denoted by pi,y(I(x)) in (3). p(x∈Oy∩Ωi) is the a priori probability of the partition Oy∩Ωi among all partitions of Oy, which is usually set equal in all partitions [12]. The a priori probability p(I(x)) 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∏i=1N∏x∈Ωi∩Oypi,y(I(x)). Taking a logarithm, the maximization can be converted to the minimization of the following energy:
(5)Ey=∑i=1N∫Oy∩Ωi-logpi,y(I(x))dx.Now, we introduce the fuzzy membership functionui(x) into the above energy. Since ui(x) represents the probability of pixel x belonging to region Ωi such that ui(x)=0 for x∉Ωi, we can rewrite (5) as
(6)Ey=∑i=1N∫Oy-uim(x)logpi,y(I(x))dx,
where m 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:
(7)Ey=∑i=1N∫Oy-uim(x)ω(x-y)logpi,y(I(x))dx,
where ω(x-y) is a nonnegative weighting function such that ω(x-y)=0 for |x-y|>r and ∫Oyω(x-y)dx=1.Although the choice of the weighting function is flexible, it is preferable to use a weighting functionω(x-y) such that larger weights are assigned to the data I(x) for x closer to the center y of the neighborhood Oy. In this paper, the weighting function ω is chosen as a truncated Gaussian kernel as follows:
(8)ω(d)={1ae-|d|2/2τ2,|d|≤r,0,else,
where τ is the standard deviation of the Gaussian kernel and a is a constant to normalize the Gaussian kernel. The energy function Ey in (7) can be rewritten as
(9)Ey=∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx,
as ω(x-y)=0 for x∉Oy.The ultimate goal is to minimizeEy for all the center points y in the image domain Ω, which directs us to define the following double integral energy functional as our proposed local fuzzy energy:
(10)EL=∫ΩEydy=∫Ω(∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx)dy.


## 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
(11)EG=∫Ω(∑i=1Nuim(x)|I(x)-ci|2)dx,
where ci is the constant value of intensity J belonging to the ith tissue in the inhomogeneity-free image. Because ui(x)(0≤ui(x)≤1) represents the probability that pixel I(x) belongs to the ith tissue, the distance at location x between the original image and the corresponding inhomogeneity-free image should be a weighted summation, that is, d(x)=∑i=1Nuim(x)|I(x)-ci|2. Adding up the distance d(x) for all x 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 energyEL with the global fuzzy energy EG, and thus obtain the proposed objective function as follows:
(12)E=β∫Ω(∑i=1Nuim(x)|I(x)-ci|2)dx+∫Ω(∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx)dy,
where β(0≤β≤1) 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 densitypi,y(I(x))(x∈Oy∩Ωi) for different cluster i in local region Oy is almost the same, thus the estimated membership functions by minimizing the energy Ey 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 pi,y(I(x)) 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:
(13)β(x)=γ·mean(LCRW(x))·(1-LCRW(x))x∈Ω,
where γ is an adjustable parameter and LCRW denotes the local contrast ratio of the given image, which is defined as
(14)LCRW(x)=(ILmax-ILmin)Ig,
where W defines the size of the local square window centered on x. we set W=5 for all experiments in this paper. ILmax and ILmin are the maximum and minimum of intensities within this local window, respectively. Ig represents the total intensity levels of the image, for example, Ig = 255 for gray image.In the above weight function (13), mean (LCRW(x)) is the average value of LCRW(x) 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 (LCRW(x)) should be chosen for the global fuzzy term. (1-LCRW(x)) 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 functionE in (12) can be minimized in a fashion similar to the standard FCM algorithm. Taking the first derivatives of E with respect to ui(x), b(y), ci and σ, and setting them to zero results in four necessary but not sufficient conditions for E to be at a local extremum. In this subsection, we will derive these four conditions.The energyE is subject to the constraint ∑i=1Nui(x)=1. Thus, this constrained optimization will be solved using one Lagrange multiplier as follows:
(15)F=β∫Ω(∑i=1Nuim(x)|I(x)-ci|2)dx+∫Ω(∑i=1N∫Ω-uim(x)ω(x-y)logpi,y(I(x))dx)dy+λ(1-∑i=1Nui(x)).Taking the derivative ofF with respect to ui(x) and setting the result to zero, we have, for m>1,
(16)[∂F∂ui(x)=βmuim-1(x)|I(x)-ci|2+∫Ω-muim-1(x)ω(x-y)logpi,y(I(x))dy-λ∂F∂ui(x)]ui(x)=ui*(x)=0.
Solving for ui*(x), we have
(17)ui*(x)=(λβm|I(x)-ci|2+∫Ω-mω(x-y)logpi,y(I(x))dy)1/(m-1).
Since ∑K=1Nuk(x)=1 for all x, we have
(18)∑k=1N(λβm|I(x)-ck|2+∫Ω-mω(x-y)logpk,y(I(x))dy)1/(m-1)=1,
or
(19)λ=m×((∫Ω-ω(x-y)logpk,y(I(x))dy∑k=1N({+∫Ω-ω(x-y)logpk,y(I(x))dy)-1}(∫Ω-ω(x-y)logpk,y(I(x))dyβ|I(x)-ck|2+∫Ω-ω(x-y)logpk,y(I(x))dy)-1)1/(m-1))m-1)-1.Substituting (3) and (19) into (17), the zero-gradient condition for the membership functions can be rewritten as(20)ui*(x)=1∑k=1N(β|I(x)-ci|2+∫Ωω(x-y)(log(2πσ)+(I(x)-b(y)ci)2/2σ2)dyβ|I(x)-ck|2+∫Ωω(x-y)(log(2πσ)+(I(x)-b(y)ck)2/2σ2)dy)1/(m-1).In a similar way, taking the derivative ofF with respect to b(y) and setting the result to zero, we have
(21)[∑i=1N∫Ωuim(x)ω(x-y)×(I(x)-b(y)ci)(-ci)dx∑i=1N∫Ωuim(x)ω(x-y)]b(y)=b*(y)=0.
Solving for b*(y), we have
(22)b*(y)=∑i=1N∫Ωuim(x)ω(x-y)I(x)cidx∑i=1N∫Ωuim(x)ω(x-y)ci2dxLikewise, taking the derivative ofF with respect to ci and σ and setting the results to zero, we have
(23)[2β∫Ωuim(x)(I(x)-ci)(-1)dx+∫∫Ωuim(x)ω(x-y)×(I(x)-b(y)ci)(-b(y))σ2dxdy∫Ωuim(x)(I(x)-ci)(-1)dx]ci=ci*=0,(24)[∑i=1N∫∫Ωuim(x)ω(x-y)×(σ2-(I(x)-b(y)ci)2)dxdy∑i=1N∫Ω∫Ωuim(x)ω(x-y)]σ=σ*=0.
Solving for ci and σ, we have(25)ci*=2β∫Ωuim(x)I(x)dx+∫∫Ωuim(x)ω(x-y)I(x)b(y)(1/σ2)dxdy2β∫Ωuim(x)dx+∫∫Ωuim(x)ω(x-y)b2(y)(1/σ2)dxdy,(26)σ*=∑i=1N∫∫Ωuim(x)ω(x-y)(I(x)-b(y)ci)2dxdy∑i=1N∫∫Ωuim(x)ω(x-y)dxdy.


## 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:
(27)hi(x)=∑y∈NB(x)ui(y),
where NB(x) represents a square window centered on pixel x in the spatial domain. A 5 × 5 window is used throughout this work. Just like the membership function, the spatial function hi(x) represents the probability that pixel x belongs to the ith 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:
(28)ui′(x)=ui(x)hi(x)∑k=1Nuk(x)hk(x).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 variablesui(x), b(y), ci, 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 clustersN, membership function ui(x), constant intensity value ci, and bias field b(y).Step 2.
Computing the weight functionβ using (13).Step 3.
Updatingσ, b(y), ui(x), and ci sequentially using (26), (22), (20), and (25), respectively.Step 4.
Computing the new membership functions incorporated with spatial information using (28).
Repeat Steps3 and 4 till termination. The iteration is stopped when the maximum difference between two constant intensity values ci 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 functionui(x) and constant intensity value ci are initialized randomly at the intervals [0,1] and [0,255], respectively. Bias field b(y) is initially equal to 1 for every pixel y. The parameters used in this paper are mostly set as the following default values: fuzzy exponent m=2, standard deviation of the Gaussian kernel τ = 4, and neighborhood radius of the Gaussian kernel r=14. 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 Figure1 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.Figure 1
Segmentation results of the proposed algorithm on three synthetic images. Column 1: original images. Columns 2 and 3: intermediate results. Column 4: final results. Column 5: the estimated bias fields.


### 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 Figure2 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.Figure 2
Comparison of segmentation results on two X-ray vessel 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: the estimated bias fields 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 Figure3. 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.Figure 3
Comparison of segmentation results on two 3T 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: the estimated bias fields by the proposed algorithm.


### 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 linkhttp://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 S1 and S2 is defined as the ratio between the areas of the intersection and the union of them, namely, JS(S1,S2)=|S1∩S2|/|S1∪S2|. We compute the JS between the segmented region S1 obtained by the algorithm and the corresponding region S2 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.Table 1
Comparison of the JS values of the segmentation results obtained by the four algorithms.
Image
Tissue
BCFCM
sFCM
CLIC
Proposed algorithm
Brain 1
WM
0.9042
0.9168
0.9361
0.9562
GM
0.8433
0.8629
0.8985
0.9133
CSF
0.8337
0.8491
0.8543
0.8958
Brain 2
WM
0.8561
0.8674
0.8683
0.9430
GM
0.8260
0.8387
0.8331
0.9098
CSF
0.7938
0.8035
0.8184
0.8274Figure 4
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.Meanwhile, we also compare the computation time of the four algorithms for the images in Figure4. 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.Table 2
CPU time (seconds) comparison for the four algorithms.
Image
BCFCM
sFCM
CLIC
Proposed algorithm
Brain 1
0.8268
0.6552
24.2894
11.0917
Brain 2
0.7176
0.4992
21.2317
9.9841


## 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 Figure1 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.Figure 1
Segmentation results of the proposed algorithm on three synthetic images. Column 1: original images. Columns 2 and 3: intermediate results. Column 4: final results. Column 5: the estimated bias fields.


## 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 Figure2 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.Figure 2
Comparison of segmentation results on two X-ray vessel 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: the estimated bias fields 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 Figure3. 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.Figure 3
Comparison of segmentation results on two 3T 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: the estimated bias fields by the proposed algorithm.


## 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 linkhttp://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 S1 and S2 is defined as the ratio between the areas of the intersection and the union of them, namely, JS(S1,S2)=|S1∩S2|/|S1∪S2|. We compute the JS between the segmented region S1 obtained by the algorithm and the corresponding region S2 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.Table 1
Comparison of the JS values of the segmentation results obtained by the four algorithms.
Image
Tissue
BCFCM
sFCM
CLIC
Proposed algorithm
Brain 1
WM
0.9042
0.9168
0.9361
0.9562
GM
0.8433
0.8629
0.8985
0.9133
CSF
0.8337
0.8491
0.8543
0.8958
Brain 2
WM
0.8561
0.8674
0.8683
0.9430
GM
0.8260
0.8387
0.8331
0.9098
CSF
0.7938
0.8035
0.8184
0.8274Figure 4
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.Meanwhile, we also compare the computation time of the four algorithms for the images in Figure4. 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.Table 2
CPU time (seconds) comparison for the four algorithms.
Image
BCFCM
sFCM
CLIC
Proposed algorithm
Brain 1
0.8268
0.6552
24.2894
11.0917
Brain 2
0.7176
0.4992
21.2317
9.9841


## 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 r, 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 r while the right figure shows the influence of the standard deviation τ. It is clear that when r>8 or τ>3, 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 9≤r≤17 and 3≤τ≤6. In our experiments, we set r=14 and τ = 4 for all test images.Figure 5
The JS values of the segmentation results obtained by using different parameter settings of the truncated Gaussian kernel.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.Figure 6
Segmentation results obtained by the proposed algorithm with different parameterγ for images with less intensity inhomogeneity. Column 1: original images. Column 2: γ = 0.001. Column 3: γ = 0.005. Column 4: γ = 0.009.Figure 7
Segmentation results obtained by the proposed algorithm with different parameterγ for images with strong intensity inhomogeneity. Column 1: original images. Column 2: γ = 0.001. Column 3: γ = 0.005. Column 4: γ = 0.009.


## 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.



---

*Source: 316546-2013-11-25.xml*