A modified possibilistic fuzzy c-means clustering algorithm for bias field estimation and segmentation of brain MR image

Computerized Medical Imaging and Graphics 35 (2011) 383–397

Contents lists available at ScienceDirect

Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag

A modified possibilistic fuzzy c-means clustering algorithm for bias field estimation and segmentation of brain MR image Ze-Xuan Ji, Quan-Sen Sun ∗ , De-Shen Xia The School of Computer Science and Technology, Nanjing University of Science and Technology, No. 200, Xiao Ling Wei Street, Nanjing 210094, China

a r t i c l e

i n f o

Article history: Received 7 July 2010 Received in revised form 27 October 2010 Accepted 9 December 2010 Keywords: Brain MR image segmentation Intensity inhomogeneity Possibilistic fuzzy c-means algorithm Anisotropic weight Coherent local and global intensity

a b s t r a c t A modified possibilistic fuzzy c-means clustering algorithm is presented for fuzzy segmentation of magnetic resonance (MR) images that have been corrupted by intensity inhomogeneities and noise. By introducing a novel adaptive method to compute the weights of local spatial in the objective function, the new adaptive fuzzy clustering algorithm is capable of utilizing local contextual information to impose local spatial continuity, thus allowing the suppression of noise and helping to resolve classification ambiguity. To estimate the intensity inhomogeneity, the global intensity is introduced into the coherent local intensity clustering algorithm and takes the local and global intensity information into account. The segmentation target therefore is driven by two forces to smooth the derived optimal bias field and improve the accuracy of the segmentation task. The proposed method has been successfully applied to 3 T, 7 T, synthetic and real MR images with desirable results. Comparisons with other approaches demonstrate the superior performance of the proposed algorithm. Moreover, the proposed algorithm is robust to initialization, thereby allowing fully automatic applications. © 2010 Elsevier Ltd. All rights reserved.

1. Introduction Magnetic resonance (MR) imaging has several advantages over other medical imaging modalities, including high contrast among different soft tissues, relatively high spatial resolution across the entire field of view and multi-spectral characteristics. Therefore, it has been widely used in quantitative brain imaging studies. Quantitative volumetric measurement and three-dimensional (3D) visualization of brain tissues are helpful for pathological evolution analyses, where image segmentation plays an important role. The size alterations in brain tissues often accompany various diseases, such as schizophrenia [1]. Thus, estimation of tissue sizes has become an extremely important aspect of treatment which should be accomplished as precisely as possible. This creates the need to properly segment the brain MR images into gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) and also to identify tumors or lesions, if present [2]. The main difficulties in brain segmentation are the intensity inhomogeneities and noise. In fact, intensity inhomogeneity occurs in many real-world images from different modalities [3,4]. In particular, it is often seen in medical images, such as X-ray radiography/tomography and MR images. For example, the intensity inhomogeneity in MR images often appears as an intensity varia-

∗ Corresponding author. Tel.: +86 13611507508. E-mail addresses: [email protected] (Z.-X. Ji), [email protected] (Q.-S. Sun). 0895-6111/$ – see front matter © 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2010.12.001

tion across the image, which arises from radio-frequency (RF) coils or acquisition sequences. Thus the resultant intensities of the same tissue vary with the locations in the image. The noise in MR images is Rician distributed and can affect significantly the performances of classification methods. The best solutions consist of either filtering the image prior to classification or embedding spatial regularization inside the classifier itself. Bias field correction refers to a procedure to estimate the bias field from the measured image so that its effect can be eliminated. Methods for intensity inhomogeneity correction in MR images have been reviewed [5,6]. Currently, there are two popular segmentation based models for intensity inhomogeneity correction: the expectation-maximization (EM) algorithm [7,8] and the fuzzy Cmean (FCM) algorithm [9–12]. In these methods, the tasks of bias field correction and segmentation are interleaved in an iterative process such that they benefit from each other to yield better results. Wells et al. [7] developed an approach based on EM algorithm for interleaved bias field correction and image segmentation. However, there are two main disadvantages of this EM approach. First, the EM algorithm is extremely computationally intensive. Second, the EM algorithm requires a good initial guess for either the bias field or for the classification estimate. Otherwise, the EM algorithm could be easily trapped in a local minimum, resulting in an unsatisfactory solution. Based on the EM framework in, Leemput et al. [8] proposed an explicit parametric model of the bias field. Instead of manual intervention, their method used a digital brain atlas that provides a priori probability maps for WM, GM, and CSF. Although this method is claimed to be more robust

384

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

than the method of Wells et al., the initialization of the parameters remains critical. Recently, Li et al. proposed a new energy minimization method based on coherent local intensity clustering (CLIC) for simultaneous tissue classification and bias field estimation of MR images [9]. The CLIC energy is convex in each of its variables, which renders the proposed method robust to initialization, thereby allowing fully automatic applications. Meanwhile, no extra effort is needed for the bias field smoothing because the smoothness of the derived estimate of the bias field is intrinsically ensured by the spatial coherent nature of the CLIC criterion function. CLIC model draws upon intensity information in local regions, therefore, can be used to segment images with intensity inhomogeneity. However, the CLIC model is not robust to the parameters. Moreover, CLIC would lead the misclassification especially for the pixels on the boundaries. On the other hand, segmentation that utilizes local scale as homogeneous criteria has been presented and applied to intensity inhomogeneity correction in [13,14]. In [13] local scale was defined as the radius of the largest ball centered at every voxel, such that all voxels within the ball satisfy a predefined homogeneity criterion. In [14], an alternative definition of local scale, referred to as generalized scale was proposed. The g-scale of any voxel in the scene was defined as the largest connected set of voxels associated with such that all voxels within this set satisfy a predefined homogeneity criterion. The generalized scale can be used for noise reduction, inhomogeneity correction, registration, and segmentation [15–17]. Standardization is a pre-processing technique which maps nonlinearly image intensity gray scale into a standard intensity gray scale through a training and a transformation step. Standardized intensity level helps segmentation and registration, also inhomogeneity correction. Once the global intensity patterns are fixed to a certain scale, then the segmentation algorithms will work better [18]. Based on the traditional FCM objective function, most improved approaches embodied regularization terms to show increased robustness of the classification to noisy images. Pham and Prince [10] modified the FCM objective function by introducing a spatial penalty, enabling the iterative algorithm to estimate spatially smooth membership functions. Meanwhile, Ahmed et al. [11] introduced a neighborhood averaging additive term into the objective function of FCM, named the algorithm bias corrected FCM (BCFCM). Liew and Yan [12] introduced the spatial constraint to a fuzzy cluster method where the inhomogeneity field was modeled by a B-spline surface, while the spatial voxel connectivity was implemented by a dissimilarity index, which enforced the connectivity constraint only in the homogeneous areas. In this way the tissue boundaries were better preserved. Although suppressing the impact of noise and intensity inhomogeneity to some extent, this method still produces misclassified small regions, especially in case of heavy noise. These approaches can overcome the impact of noises and intensity inhomogeneity, but it computes the neighborhood term in every iteration step, giving the algorithm a computational load. In order to overcome these defects, Cai et al. [19] introduced a new local similarity measure, combining spatial and gray level distances, applied it as an alternative pre-filtering to EnFCM, and named this approach fast generalized FCM (FGFCM). This method is able to extract local information causing less blur than averaging filter. But it still has an experimentally adjusted parameter and the precision of the segmentation is not good enough. Szilágyi et al. [20] modified the FGFCM (MFGFCM) to improve the precision of segmentation. Similar to this anisotropic weight construction method, Kang et al. [21] improved FCM with adaptive weighted averaging filter (FCM AWA), and Kang et al. [22] proposed a spatial homogeneity-based FCM (SHFCM). Wang et al. [23] incorporated both the local spatial context and the non-local information into the standard FCM cluster algorithm (LNLFCM), using a novel dissimilarity in place of the usual distance metric. These approaches

could overcome the noise impact, but the intensity homogeneity cannot be handled at the same time. While stable, FCM-based algorithms are known to be vulnerable to outliers/noise. To address this, possibilistic clustering, which is pioneered by the possibilistic c-means (PCM) algorithm [24], is developed and has been shown to be more robust to outliers as compared with FCM. However, the robustness of PCM comes at the expense of the algorithm’s stability [25]. Early PCM-based algorithms are known to suffer from the coincident cluster problem, which makes them too sensitive to initialization [25,26]. There have been many efforts to improve the stability of possibilistic clustering [27,28]. Among the recent ones is the possibilistic fuzzy c-means (PFCM) algorithm [29]. PFCM provides an easily expanded hybrid model of FCM and PCM. By combining FCM and PCM, PFCM can simultaneously exhibit the FCM’s stability (i.e., invulnerability to the coincident-cluster problem) and the PCM’s robustness to outliers. In addition, the possibilistic-fuzzy memberships produced by PFCM can provide a better insight into how objects are distributed. These desirable properties of PFCM make it suitable to be a basic model of the proposed algorithm. However, PFCM is only robust to outliers when estimating the centroids. This kind of algorithm cannot label the outliers accurately. Therefore, this drawback will reduce the clustering performance in real applications. In this paper, we modify the PFCM algorithm to segment MR images with intensity inhomogeneities and noise. Firstly, we propose an adaptive method to compute the weights for the neighborhood of each pixel in the image. This adaptive method can not only overcome the effect of the noise effectively, but also prevent the edge from blurring. Secondly, we propose a novel energy minimization method that combines both local and global intensity information to address intensity inhomogeneity. The segmentation target therefore is driven by two forces: one induced by the coherent local intensity and the other by the coherent global intensity. We define an energy that depends on the bias field, membership functions, typicality, label function and centroids. Bias field estimation and image segmentation are simultaneously achieved as the result of minimizing this energy. A salient advantage of our method is that its result is independent of initialization and robust to the parameter setting, which allows more accuracy intensity inhomogeneity corrections and image segmentations. The remainder of this paper is organized as follows. Section 2 provides background information on bias field formulation and clustering algorithms. The proposed method is introduced in Section 3. The implementation and results are given in Section 4, followed by some discussions in Section 5. This paper is summarized in Section 6. Note that part of results in this paper was reported in our recent conference paper [30]. 2. Related work 2.1. Bias field formulation The intensity inhomogeneity in a MR image can be modeled as a multiplicative component of an observed image described by I = bJ + n

(1)

where I is the measured image intensity, J is the true signal to be restored, b is an unknown bias field, and n is additive noise. The goal of bias correction is to estimate the bias field b from the measured intensity I. Without loss of generality, the bias field b is slowly varying in the entire image domain. Ideally, the intensity J in each tissue should take a specific value vi of the physical property being measured (e.g. the proton density for MR images). This property, in conjunction with the spatially coherent nature of each tissue’s dis-

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

tribution, implies that the true signal J is approximately a piecewise constant map. In addition, the additive noise n can be assumed to be zero-mean Gaussian noise.

The  i can be defined as

i =

n m u D k=1 ik ik  n m

The FCM clustering algorithm was first proposed by Dunn [31] and promoted as the general fuzzy c-means clustering algorithm by Bezdek [32]. The main purpose of FCM algorithm is to make the vector space of a sample point be divided into a number of sub-spaces in accordance with distance measure [33]. However, FCM algorithm fails to deal with significant properties of images, since neighbor pixels are strongly correlated, which leads to strong noise sensitivity. To overcome this weakness, Krishnapuram and Keller [26] proposed a new clustering algorithm named PCM. PCM relaxes the column sum constraint of fuzzy membership matrix in FCM and introduces a possibilistic partition matrix, so that possibilistic memberships may reflect the typicalities of data points to their clusters. Compared with FCM, PCM can effectively eliminate the influence of noise and outliers on clustering results. However, firstly the price PCM pays for its freedom to ignore noisy points is that PCM is very sensitive to initializations, and often results in the coincident cluster problem [25,28]. Secondly, typicalities, i.e., possibilistic memberships, are very sensitive to the choices of the additional parameters of PCM, which directly decide the clustering accuracy. To overcome the weaknesses of the original PCM algorithm, Pal et al. [29] combined the objective functions of PCM and FCM into a new objective function and presented an improved version, called PFCM, which can be interpreted as PCM and FCM, respectively, in some special cases where some proper parameters were adopted. So, PFCM can inherit the merits of both clustering algorithms. The algorithm divides the data set I ={I1 , I2 , . . ., In } into c clusters and n is the number of all the pixels in the image. Let the membership function uik , uik ∈ [0, 1] show the degree of the pixel Ik , k = 1, 2, . . ., n belonging to cluster i(1 ≤ i ≤ c). Then the result can be denoted by a matrix of fuzzy membership function matrix U = [uik ]c×n . We represent typicality by tik , tik ∈ [0, 1] and the typicality matrix by T = [tik ]c×n . According to the definition of the theory, we have  c u = 1 for every pixel in the image. The objective function i=1 ik to be minimized is: n c  



(CF um + CT tik ) × ||Ik − i ||2 + ik

i=1 k=1

c n  

i

i=1

(1 − tik ) (2)

k=1

where V ={1 , 2 , . . ., c } is the characterized intensity center. The parameters CF > 0, CT > 0, m > 1,  > 1, The  i > 0 are user defined constants. The constants CF and CT define the relative importance of fuzzy membership and typicality values in the objective function. Note that uik has the same meaning of membership as that in FCM. Similarly, tik has the same interpretation of typicality as in PCM. Let, the objective function of PFCM can get the minimum by updating the membership U, the typicality T and the cluster centers V as follows:

uik =

tik =

⎧  1/(m−1) ⎫−1 c ⎨ ⎬ D ik

⎩

j=1

⎭

Djk 1

1 + ((CT /i )Dik )

1/(−1)

n  (CF um + CT tik )Ik ik i = k=1 n  m k=1

(6)

u

k=1 ik

2.2. PFCM

JPFCM =

385

(CF uik + CT tik )

(3)

(4)

(5)

2.3. CLIC Li et al. proposed a new energy minimization method based on coherent local intensity clustering (CLIC) for simultaneous tissue classification and bias field estimation of MR images [9]. For the intensity Ik in the neighborhood ok , a clustering criterion function is defined, with the clustering centers i replaced by bk i . In addition, a weight for each pixel is introduced to control its influence on the clustering criterion function. The intensity Ik at location k far away from the neighborhood center should have less influence in the clustering criterion function than the locations close to the neighborhood center. The clustering criterion function is shown as follow

JCLIC =

n c   i=1 k=1

um ik



K(r − k)|Ik − b(r)i |2

(7)

r ∈ ok

where U ={uik } are the membership functions for c regions (tissues), and K(r − k) is the weight assigned to the target pixel and the weighting function K is a truncated Gaussian kernel

 K(u) =

1 −|u|2 /2 2 e , a 0,

for |u| < 

(8)

else

where  is the standard deviation of the Gaussian kernel,  is the radius of the neighborhood window and a is a constant to normalize the Gaussian kernel. Tissue classification and bias field estimation are simultaneously achieved by minimizing this energy. The smoothness of the derived optimal bias field is intrinsically ensured by the spatially coherent nature of the CLIC criterion function. As a result, no extra effort is needed to smooth the bias field in the method. Moreover, the CLIC algorithm is robust to the choice of initial conditions, thereby allowing fully automatic applications. However, the CLIC algorithm is not robust to the choice of the parameter in the algorithm. Moreover, from Eq. (7) we can observe that the CLIC model introduce the local intensity into the objective function with Gaussian kernel. The membership function cannot be estimated accurately because not all the pixels in the kernel belong to the same tissue, which would make the estimated membership functions very ‘fuzzy’, and finally cause misclassifications, especially for the pixels around the boundaries.

3. Methods In this section, we modify the PFCM algorithm to segment MR images with intensity inhomogeneity and noise. To segment the pixels with noise in the image, we introduce an adaptive method to compute the weights for the neighborhood of each pixel. Meanwhile, we define the global intensity and introduce this global intensity information into the objective function to make our method take the local and global intensity information into account. Therefore, a model based on coherent local and global intensity clustering is proposed to handle intensity inhomogeneity and segment the brain MR images, simultaneously.

386

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

3.1. Objective function

3.2. Incorporation of local spatial continuity

We define the energy minimization function as

JMPFCM =

n c  

⎡  (CF um +CT tik ) ik

+(1−˛)

 ωks

˛

s ∈ Nk

i=1 k=1



⎣



⎞⎤

K(r−k)||Iks −b(r)i ||2 ⎠⎦ +

r ∈ Ok



Li (r)||Iks −b(r)i ||2

r ∈˝ c n  

i

i=1

(1 − tik )

(9)

k=1

where U = [uik ] is the membership functions, T = [tik ] is the typicality matrix, V ={1 , 2 , . . ., c } is the characterized intensity center. ωks ∈ [0, 1] is the local spatial continuity weights detailed in the next subsection. Nk is the neighborhood of local spatial continuity for each pixel in image and the radius is r; Ok is the neighborhood of intensities Ik to estimate the bias field and the radius is . ˝ stands for the whole image. The parameters CF > 0, CT > 0, m > 1,  > 1. The constants CF and CT define the relative importance of fuzzy membership and typicality values in the objective function. Note that uik and tik have the same meaning of membership and typicality as those in PFCM. n is the number of all the pixels in the image and c is the number of clusters.Li (r), i = 1, . . ., c ; r = 1, . . ., n is a label function defined as



The incorporation of local spatial continuity considers the influence of neighboring pixels on the center pixel of interest during classification [12]. Let r be the radius of the neighborhood Nk centered as pixel k. Denote s be the pixel in the neighborhood Nk . There are (2r + 1)2 pixels in each neighborhood obviously, but not all the pixels are useful for similarity computation. In other words, different pixels in the neighborhood should have different weights. Gaussian kernel is the most common method to describe the impact of different pixels. However, this method may blur the edges of the image. Szilágyi et al. [20] and Kang et al. [22] all concentrate on the weights computation in the neighborhood for each pixel. But these filter methods are complex and need at least one parameter to tune. In this section we propose a novel method to incorporate the local spatial continuity into the objective function, in order to improve the accuracy of segmentation and make the filter method easy and simple. In image processing, standard deviation is in common use. But it only describes the difference between the target pixel and mean values or central pixel. Motivated by standard deviation and FGFCM [19], we define the filter coefficients as (s)

(g)

ωks = ωks · ωks

(12)

All neighbor pixels will have coefficients ωks ∈ [0, 1], depend(s) ing on their space distance and gray level difference. ωks reflects (g)

i = max(uir ) (10)

the local spatial relationship, while ωks reflects the local gray-level relationship.

pixels in the ith tissue. The first term N i is the total number 2 stands for the global intensity; the second L (r)||I − b(r) || i i k r∈˝ term K(r − k)||Ik − b(r)i ||2 is the same as the term in Eq. r ∈o

(1) Construction of local spatial relationship. For each pixel k with coordinate (pk , qk ), the spatial component reflects the damping extent of the neighbors with the spatial distances from the central pixel and defined as

l (r) Li (r) = i ; Ni

li (r) =

1, 0,

i

i= / max(uir ) i

k

(3) which stands for the local intensity; and ˛ is a positive constant (0 ≤ ˛ ≤ 1) to balance the influence of local intensity and global intensity. In our model, the  i can be defined as i =

n m u D l (k) k=1 ik ik i  n m u l (k) k=1 ik i

(11)

For the existence of the second term, the value of b(x) depends on all the intensities Ik in the neighborhood Ok . This indicates the smoothness of estimated bias field. However, the membership function cannot be estimated accurately because not all the pixels in the neighborhood Ok belong to the same tissue, which would make the estimated membership functions very ‘fuzzy’, finally, cause misclassifications, especially for the pixels around the boundaries. Therefore, we utilize the first term to introduce the global intensity by computing the distance between gray value of pixel k and the centroids with the pixels belonging to the same tissue with the pixel k. It should be noted that the typicality and membership in the objective function play different roles during segmenting an image. Typicality is used to alleviate the undesirable effects of outliers/noise, which means that the centroids can be estimated accurately for the introduction of typicality. Membership function reflects the degrees of the pixel belonging to each cluster, which means that we can get the final segmentation with membership values. If CT = 0 and  i = 0, the proposed model degenerate to the FCM optimization problem; while CF = 0 converts it to the PCM based model. If ˛ = 0, our model can degenerate to the energy minimization method based on CLIC. The proposed model can be treated as a general framework of FCM, PCM, CLIC and PFCM to segment the MR images with intensity inhomogeneities and noise.

(s)

ωks = exp[−max(|pk − ps |, |qk − qs |)]

(13)

(s)

ωks makes the influence of the pixels within the local window change flexibly according to their distance from the central pixel and thus more local information can be used. This local spatial relationship reflects the damping extent of the neighbors with the spatial distances from the central pixel. (2) Construction of local gray-level relationship. We firstly get the root of mean-square deviation  ks for each pixel s: for each pixel s, we compute the total distance between s and other pixels in the neighborhood, square this distance, average the values and take the square root. This kind of distance computing can reflect the difference between the target pixel and the other pixels in the neighborhood. When the interested pixels in the neighborhood are superior in numbers, such as for the pixels in the homogenous areas, around the edges and around the corners with obtuse angles,  ks can be used to compute the weights accurately to get rid of the impact of noise and preserve the edges or corners with obtuse angles. Otherwise, for the pixels around the corners with acute angles and in the areas with micro-texture, this method will treat the target pixel as noise and lead misclassification. To overcome this defect, for each pixel in the neighborhood, the intensity difference from the central pixel is added to get the final  ks . Then we project  ks into kernel space. Finally, the weights are normalized. Due to the fast decay of the exponential kernel, large distance between  ks and the mean of these mean-square deviations will lead to nearly zero weights. The formulas show as follows

 ks =

s ∈ Nk \(s)

(Is − Is )2

nk − 1

1/2 + (Is − Ik )

2

(14)

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

387

Fig. 1. The weights of a synthetic image in different conditions. Part A is a synthetic image without noise and Part B is a synthetic image with noise. The upper number in each cell represents the intensity value, while the lower number shows the obtained weight.

⎛ ks = exp ⎝−

(g)

ωks =



ks −

ks

 s ∈ Nk ks





s ∈ Nk



ks /nk

⎞ ⎠

(15)

(16)

where  ks is the mean-square deviation of pixel r from other pixels in the neighborhood Nk centered as pixel k. Its value reflects gray value similarity degree for the target pixel in the local window. The smaller its value is, the higher similarity degree is, and vice versa. The motivation behind the use of pixeldependent parameter  ks is to exploit the local statistics varying with each image pixel in the neighborhood. can make  ks be the same if the image values are multiplied with a constant. In this paper, is the greatest common divisor of the gray values in (g) the neighborhood, then  ks is projected between 0 and 1. ωks is the weight of pixel s, xs is the intensity of pixel s, nk is the potency of Nk . Through the formulas above, the pixels with similar intensity can be separated as a class and then we can get the final gray level of the center pixel k by the weight ωks . Fig. 1 shows the weights of pixels in different conditions. In Part A of Fig. 1, there are four homogenous regions and the gray value are 25, 50, 100 and 200 from left to right, respectively. The coefficients for the neighborhood of three pixels around the edge are computed. For the introduction of parameter , the corresponding weights of these three regions are the same. Meanwhile, the proposed method can exactly restore the gray value of center pixel and makes the boundaries preserved perfectly for the image without noise. In Part B of Fig. 1, the synthetic image in Part A is added with Gaussian noise. The proposed method shields the impact of the edge and noise pixels and makes the estimate of center pixel much more accurate.

Fig. 2 shows the weights of a MR image with noise in different conditions. Pixel A, B and C locate in the homogenous area of WM, CSF and GM, respectively; Pixel D, E and F locate around the edge between WM and CSF, edge between WM and GM and edge between GM and CSF, respectively; Pixel G is around the corner of CSF; Pixel G is around the corner of GM; Pixel I locate around the micro-texture of CSF. The upper number in each cell represents the intensity value, while the lower number shows the obtained weight. From the lower number of each cell, we can find that the proposed method can get rid of noise effectively (the weight of the noisy pixel approaches zero), while preserve the edge and corner of the image and get a much more accurate weight for each pixel. MFGFCM [20], FCM AWA [21] and SHFCM [22] all concentrate on the construction of the weights. Unfortunately, Both MFGFCM and FCM AWA need two thresholds which can only get by experiment. On the other hand, FCM AWA and SHFCM would not get good results if the noise is heavy. LNLFCM [23] introduced non-local spatial information and incorporated both the local and non-local spatial context into the standard FCM cluster algorithm. However, this improved method is computationally expensive and can not segment the images accurately when the similar structure can not be found. The weight computation of our method is more effective and simple, and has a better stability and accuracy for different noise.

3.3. Parameter estimation By changing the order of integration, the above energy JMPFCM can be written in the form JMPFCM =

n c   i=1 k=1



(CF um + CT tik )Dik + ik

c n  

i

i=1

k=1

(1 − tik )

(17)

388

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

Fig. 2. The weights of a MR image with noise in different conditions. The upper number in each cell represents the intensity value, while the lower number shows the obtained weight. Pixel A: homogenous WM; Pixel B: homogenous CSF; Pixel C: homogenous GM; Pixel D: edge between WM and CSF; Pixel E: edge between WM and GM; Pixel F: edge between GM and CSF; Pixel G: corner around CSF; Pixel G: corner around GM; Pixel I: micro-texture around CSF.

where



Dik =

 ωks

s ∈ Nk



˛

2



Li (r)||Iks − b(r)i || + (1 − ˛)

r ∈˝

2

K(r − k)||Iks − b(r)i ||

r ∈ Ok

(18) The minimization of the energy JMPFCM is subject to the constraint c 

uik = 1,

uik ∈ [0, 1], tik ∈ [0, 1]

(19)

i=1

The four variables U, T, b and V in the energy function JMPFCM can be estimated by minimizing the energy function given the other three variables remain fixed. Thus, minimization of the energy JMPFCM can be performed by an interleaving process of minimization with respect to the variables U, T, b and V. Moreover, the energy minimization algorithm is robust to the initialization of the variables U, T, b and V. The segmentation results can be gotten when the convergence criterion reaches and the estimated bias field and true signal are given by b and I/b, respectively. The updating function of U can be shown as

⎧  1/(m−1) ⎫−1 c ⎨ ⎬ D ik

⎩

j=1

⎭

Djk

(20)

The updating function of T can be shown as tik =

1 1 + ((CT /i )Dik )

(21)

1/(−1)

The updating function of b can be shown that

b=

c  (2) Jg = (CF um + CT ti )i2 Li , i i=1 c  (2) m 2

(1) ˛˜I Jg + (1 − ˛) (2)

˛Jg + (1 − ˛)



˜I J (1) l



(2)

Jl





∗K

∗K





(22)

(1)

Jl

=

=

c

c

i=1

i=1



(CF um + CT ti )i Li , i 

(CF um + CT ti )i i

and

(C u + CT ti )i . Note that ui and ti are the memJl = i=1 F i bership and typicality of the ith tissue over the whole image, respectively. ˜I is the filtered image with the proposed method to compute the weights for each pixel. The updating function of V can be written as i =

n  (CF um + CT tik )˜Ik (˛bLi + (1 − ˛)(b ∗ K)) ik nk=1  m 2 2 k=1

uik =

(1)

where * is the convolution operation. Jg

(CF uik + CT tik )(˛b Li + (1 − ˛)(b ∗ K))

(23)

where ˜Ik is the gray value of kth pixel in the filtered image ˜I . The entire procedure of minimization of the energy JMPFCM is described as below Step 1. Initialization. Initialize U, T, b, V and the parameters used in the proposed algorithm. Step 2. Update U given by Eq. (20). Step 3. Update L given by Eq. (10). Step 4. Update  given by Eq. (11). Step 5. Update T given by Eq. (21). Step 6. Update b given by Eq. (22). Step 7. Update V given by Eq. (23). Step 8. Determine the termination of iteration. Check convergence criterion. If convergence has been reached, stop the iteration, otherwise, go to Step 2. Remark. In this work, the proposed algorithm and all other compared methods are all fed only with the intensity of the pixels. So the intensity value of pixels is the only feature vector that is employed for the clustering. Our work in this paper focuses on the improvement for the accuracy of segmentations and estimation of bias field, so the number of clusters is assigned manually. Moreover, since the number of cluster has been set, the initialization of each cluster is assigned randomly, e.g. if c = 4, we will randomly create 4 values that lie between the minimum and maximum intensity value of the whole image and assign these four values to be the initialization of each cluster. In Section 5,

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

389

Table 1 The description of the dataset used in the experiments. Data

Type/protocol

Scene domain

Noise level

Inhomogeneity level

D1 D2 D3

Synthetic Synthetic T1 Real T1

256 × 256 181 × 217 × 181 255 × 255 × 108

0–5% 0–9% Round 3%

NAN 0–100% Round 60–80%

Fig. 3. Segmentation results on phantom images.

we will show the the initialization robustness of the proposed algorithm. We will discuss the choice of the parameters used in our paper in Section 5. Unless otherwise specified, we use the following default setting of the parameters in the experiments: CF = 1, CT = 1, m = 2  = 2, ˛ = 0.5, the radius of neighborhood Nk is r = 1, the parameters of Gaussian kernel are  = 4 and  = 10 for all images in this work. 4. Experimental results In this section, we test and compare the proposed method (MPFCM) with some other reported algorithms on several synthetic images and synthetic brain MR images from two aspects (anti-noise ability and bias field estimation). Our model is also extended to segment 3D MR images.

noise levels, and levels of intensity non-uniformity. More importantly, BrainWeb provides the ‘ground truth’ for each simulated 3D brain MR image. In our experiments, we use the T1 brain MR images with different level of noise and intensity inhomogeneity to test the anti-noise ability of the algorithms and to estimate the bias field in the image. D3 is a series of real brain MR images from BRAINMAPS [36]. BRAINMAPS provides real adult T1 whole-brain MRI images with expert segmentations. The noise level is around 3% and the intensity inhomogeneity is about 60–80%. Moreover, we use some real 3 T and 7 T brain MR images with noise and intensity inhomogeneity to test the performance of the proposed algorithm. However, there do not exit the ‘ground truth’ for these images. We can only compare the segmentations visually. Fortunately, for heavy intensity inhomogeneity exists in these images, the visual comparisons are possible.

4.1. Describing of dataset In this paper, we mainly use three dataset for experiments, which have been listed in Table 1. D1 is a dataset with 256 × 256 synthetic images taken from IBSR [34]. D1 consists of 4 increasingly complex shapes with three different levels of signal-to-noise and three different levels of contrast-to-noise. There are 3 clusters for each image. In our work, we only use one kind of image with different level of noise to test the anti-noise ability of the algorithms. D2 is a dataset of simulated 3D brain MR images from BrainWeb [35]. BrainWeb provides full three-dimensional data volumes which have been simulated using three sequences (T1-, T2-, and proton-density- (PD-) weighted) and a variety of slice thicknesses,

4.2. Quantitative evaluation The experimental results on synthetic phantoms and MR images can be evaluated with the following equations. (1) The evaluation of synthetic phantoms A ground truth is used to quantify the segmentation quality of synthetic phantoms and the Kappa Index (KI) overlap measure is used KI =

2 · TP 2 · TP + FP + FN

(24)

Table 2 Quantitative accuracy (KI, expressed as %) of segmentations on phantom images with different methods. Bold value shows the highest quantitative accuracy in each line. Images

FCM

PFCM

FCMAWA

SHFCM

MFGFCM

LNLFCM

MPFCM

G1 G2 G3 Impulse Mixed1 Mixed2 Average

0.9938 0.8397 0.7946 0.9733 0.9672 0.8271 0.8993

0.9941 0.8403 0.7946 0.9733 0.9671 0.8273 0.8995

0.9986 0.9673 0.9056 0.9735 0.9602 0.9271 0.9553

0.9993 0.9733 0.9061 0.9591 0.9585 0.9331 0.9549

0.9997 0.9968 0.9812 0.9665 0.9664 0.9590 0.9783

0.9980 0.9504 0.8905 0.9733 0.9719 0.9531 0.9562

0.9999 0.9986 0.9985 0.9996 0.9991 0.9971 0.9988

390

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

Fig. 4. The segmentations of brain MR images by using the algorithms FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM. The brain images were simulated with 9% noise. Row 1: original images; Row 2: FCM AWA; Row 3: SHFCM; Row 4: MFGFCM; Row 5: LNLFCM; Row 6: MPFCM; Row 7: ground truth. The three figures on the right show JS values of these methods for MR images with increasing noise. From top to bottom, the figures show the JS value of GM, WM and CSF, respectively.

where TP is the number of true positives, FP is the number of false positives and FN is the number of false negatives. (2) The evaluation of synthetic brain MR images The common method, Jaccard similarity (JS), is used to compute the quantitative evaluation of the segmentations. It defines as the ratio between intersection and union of two sets S1 and

S2 , representing the obtained and gold standard segmentations, respectively. The closer the JS values to 1, the better the segmentation and bias correction. The equation is shown as follows

JS(S1 , S2 ) =

|S1 ∩ S2 | |S1 ∪ S2 |

(25)

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

391

4.3. Anti-noise ability Segmentation results on phantom images are shown in Fig. 3. The noise removal performances are compared on D1 and the experiments compare the accuracy of the algorithms FCM, PFCM, FCM AWA, SHFCM, MFGFCM, LNLFCM with MPFCM. In this experiment, we set the radius of neighborhood r = 1 for the algorithms FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM. And we use the suggested values for the tradeoff parameters in different algorithms. From Fig. 3, it is obvious that FCM is very sensitive to the noise, while the result and efficiency of PFCM are not satisfied. Though FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM all provide better segmentations. There exist obvious misclassification pixels in FCM AWA. The result of LNLFCM has obvious misclassification pixels around the corner, because the number of similar pixels around corners is small. This makes the ability of anti-noise limited for the pixels around corners. Visually, the proposed method achieves the best result, slightly over SHFCM, MFGFCM and LNLFCM in details, and significantly over all others. The quantitative accuracy of the results is shown in Table 2. Here, different types of noise (Gaussian noise, salt-and-pepper impulse noise, and mixtures of these) are used, where G1, G2, G3 are the Gaussian noise with 0 mean and increasing STD  = 5,  = 15 and  = 25, respectively. Impulse is a salt-and-pepper impulse noise. Mixed1 and Mixed2 are the mixed noise with increasing standard deviation (STD). Table 2 demonstrates that the proposed method is suitable for segmenting images corrupted with unknown noises and gets higher accuracy. Table 2 shows that the methods FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM all have the accuracy above 98%. Therefore, we only use these five algorithms to segment the MR images with different noise level and the MR images are obtained from the Brain Web Simulated Brain Database [35]. A detailed view containing five segmentations, FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM, is shown in Fig. 4. In Fig. 4, all the brain images are got from D2 with T1-weighted contrast, 1-mm cubic voxels, 0% image intensity inhomogeneity and 9% noise. From the second row to the sixth row, the figures show the corresponding results with FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM, respectively. The images in the last row show the anatomical model of the segmented brain image. Through the comparison between the five algorithms FCM AWA, SHFCM, MFGFCM, LNLFCM and MPFCM, the proposed algorithm performs the best in denoising ability especially in the area with abundant textures and details. The FCM AWA can not get an accurate segmentation. The segmentations of SHFCM, MFGFCM and LNLFCM have some misclassification pixels. In contrast, the proposed method obtains a higher accuracy in these areas. The three figures on the right of Fig. 4 show JS values of these methods for MR images with increasing noise. From top to bottom, the figures show the JS value of GM, WM and CSF, respectively. It can be seen that our method achieves more accurate results.

Fig. 5. Applications of our method to 3 T brain MR images. Column 1: original 3 T brain MR images. Column 2: estimated bias fields. Column 3: bias corrected images. Column 4: segmentation results.

images shown in the first column in Fig. 6 without stripping the skull. Desirable results of bias field image, bias corrected image and tissue segmentation are obtained as shown in the seconde, third and fourth column in Figs. 6, respectively. To compare the proposed algorithm with other algorithms, we use the dataset D2 with ground truth. We compare our method with the methods of Wells et al. [7], Leemput et al. [8], Li et al. [9] and SMH [18]. We tested these methods on three MR brain images with noise 3% and synthetic intensity non-uniformity (INU). The segmentation results for the three noisy images are shown in Fig. 7. The original noisy images are shown in the first column. The segmentations of Wells’ method, Leemput’s method, Li’s method, SMH and the proposed method are shown in the second, third, fourth fifth and sixth column, respectively. The last column shows the ground truth of the three MR images. To evaluate the accuracy of segmentation, we compute the JS between the segmented region S1 by the algorithm and the corresponding region S2 given by the ground truth. The typical noise levels corresponding to real MR data would be around 3%. So we tested these methods on 20 MR brain images with noise 3% and different level of intensity non-uniformity (INU), from 20% to 100%. The average JS values for WM, GM and CSF of the

4.4. Bias field estimation We firstly show the results of our method for the 3 T MR images in the first column of Fig. 5. The estimated bias fields, bias corrected images, and the segmentation results are shown in the second, third, and fourth columns, respectively. In the bias corrected images in the third column, the intensities within each tissue become quite homogeneous. From the segmentation results of the fourth column, we can observe that the tissue segmentation results are quite consistent with the desire. It is worth noting that the proposed algorithm does not rely on skull stripping. For example, we apply our model to 7 T brian MR

Fig. 6. Applications of our method to 7 T brain MR images. Column 1: original 7 T brain MR images. Column 2: estimated bias fields. Column 3: bias corrected images. Column 4: segmentation results.

392

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

Fig. 7. Comparison results for three synthetic MR images with noise 3% and synthetic INU. Column 1: original MR images. Column 2: segmentations of the Wells’ method. Column 3: segmentations of the Leemput’s method. Column 4: segmentations of the Li’s method. Column 5: segmentations of SMH. Column 6: segmentations of the proposed method. Column 7: ground truth. Table 3 The quantitative evaluation (JS) for Wells’ method, Leemput’s method, CLIC, SMH and the proposed method. Images

Class

Wells et al.

Leemput et al.

CLIC

SMH

MPFCM

20%

GM WM CSF

0.8300 0.8951 0.8792

0.8181 0.8741 0.8642

0.8314 0.8952 0.8855

0.8563 0.9180 0.8890

0.9043 0.9360 0.8912

40%

GM WM CSF

0.8729 0.9282 0.8809

0.8449 0.8847 0.8875

0.8968 0.9349 0.8878

0.8559 0.9178 0.8851

0.9048 0.9384 0.8962

60%

GM WM CSF

0.8233 0.8984 0.8638

0.8171 0.8695 0.8665

0.8416 0.8998 0.8722

0.8558 0.9176 0.8844

0.8644 0.9253 0.8920

80%

GM WM CSF

0.8108 0.8873 0.8445

0.8145 0.8662 0.8643

0.8380 0.8984 0.8680

0.8528 0.9148 0.8758

0.8618 0.9252 0.8867

100%

GM WM CSF

0.7553 0.8667 0.7876

0.8072 0.8633 0.8589

0.8327 0.8990 0.8677

0.8489 0.9112 0.8716

0.8513 0.9210 0.8862

Fig. 8. Comparison of our method with the methods of Wells et al., Leemput et al. on real brain MR images. Column 1: original images; column 2: Wells et al.; column 3: Leemput et al.; column 4: CLIC; column 5: SMH; column 6: our method; column 7: ground truth.

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

393

Fig. 9. JS values for WM (left) and GM (right). The x-axis represents four images in Fig. 8 in the same order.

Fig. 10. Comparison results for two synthetic images shown in the fist column in each part. The estimated bias fields, bias corrected images and segmentation results, are shown in the second, third, and fourth columns, respectively. The results of the methods of Wells et al., Leemput et al., Li et al., SMH and the proposed method are shown in rows 1, 2, 3, 4, respectively.

394

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

Fig. 11. Results for 3D brain MR images. For each column, the left part shows the surfaces of GM and WM, and the right part shows the corresponding 2D depth maps of GM and WM. Column 1: ground truth of GM and WM surfaces; column 2: GM and WM surfaces of brain MR images with 1% noise; column 3: GM and WM surfaces of brain MR images with 5% noise; column 4: GM and WM surfaces of brain MR images with 9% noise.

Table 4 Quantitative accuracy (JS) of the segmentations on brain MR images with different radius. Noise level

Class

r=1

r=2

r=3

r=4

r=5

1%

GM WM CSF

0.9301 0.9697 0.8948

0.8755 0.9488 0.8392

0.8410 0.9283 0.8050

0.8165 0.9126 0.7813

0.7952 0.8977 0.7595

3%

GM WM CSF

0.9044 0.9592 0.8897

0.8650 0.9439 0.8292

0.8341 0.9257 0.7965

0.8112 0.9093 0.7817

0.7932 0.8975 0.7610

5%

GM WM CSF

0.8838 0.9454 0.8788

0.8536 0.9349 0.8239

0.8267 0.9211 0.7929

0.8049 0.9074 0.7744

0.7853 0.8951 0.7546

7%

GM WM CSF

0.8621 0.9329 0.8599

0.8374 0.9256 0.8120

0.8156 0.9145 0.7873

0.7944 0.9023 0.7629

0.7763 0.8911 0.7514

9%

GM WM CSF

0.8553 0.9255 0.8430

0.8146 0.9090 0.8015

0.7911 0.8965 0.7742

0.7706 0.8830 0.7590

0.7529 0.8732 0.7463

11%

GM WM CSF

0.7151 0.8386 0.7244

0.7376 0.8613 0.6995

0.7004 0.8309 0.6463

0.6736 0.8182 0.6008

0.6520 0.7977 0.5741

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

395

Fig. 12. Axial view of the 3D segmentation and bias correction results. Row 1: original images; Row 2: ground truth; Row 3: bias corrected images; Row 4: segmentation results; Row 5: membership function of WM; Row 6: membership function of GM.

results obtained by the five methods are listed in Table 3 . It can be observed from both Fig. 7 and Table 3 that the segmentation results of our method are more accurate than the other four methods. Fig. 8 shows the experimental results of the images from D3 (see the first column). The segmentation results obtained by the methods of Wells et al., Leemput et al., CLIC, SMH and the proposed method are shown in columns 2, 3, 4, 5 and 6, respectively. The ground truth segmentations are obtained from expert manual segmentation, as shown in the sixth column. It can be seen that the results of our model and the methods of Wells et al., Leemput et al. CLIC and SMH look similar by visual comparison. However, we can show by quantitative comparison that our model produces more accurate results. Fig. 9 shows the JS values for WM and GM of these methods. It can be seen that our method achieves more accurate results. Fig. 10 shows the comparison results for two synthetic images with severe intensity inhomogeneities. The results of the methods by Wells et al., Leemput et al., Li et al., SMH and the proposed one are shown in the first, second, third, fourth and fifth row in each part, respectively. While it is difficult to visually compare the bias

corrected images, the segmentation results of our method are more accurate than the other three methods, especially in the lower part of the image.

4.5. 3D MR brain images Our model can be easily extended to segment 3D MR images with intensity inhomogeneities and noise. Fig. 11 shows the surfaces of the GM and WM segmentation of 3D brain MR images with 40% intensity inhomogeneity and increasing noise level (1%, 5% and 9% from Column 2 to Column 4) and all these 3D images are got from D2. It can be observed that our method can get accurate segmentations of 3D brain MR images. To demonstrate the advantage of our method clearly, we show five slices and the corresponding segmentations obtained by our method in Fig. 12. Fig. 12 also shows the membership function results of the proposed algorithm by segmenting five sequential brain MR images with 9% noise and 40% intensity inhomogeneity. It can be clearly seen that our method segments the WM, GM, and CSF accurately.

396

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

Fig. 13. JS results of the segmentation with different parameter selections of Gaussian kernel. Left: the influence of the standard deviation . Right: the influence of the radius .

5. Discussion

5.3. The parameter ˛

In this section, the influence of parameters in our algorithm is discussed: the selection for the radius r of Nk neighborhood, the two parameters selections of Gaussian kernel (the standard deviation  and the radius of the neighborhood window ) and the parameter ˛. Meanwhile, we will discuss the initialization robustness of the proposed model. It should be emphasized that we do not discuss the two parameters CF and CT , which define the relative importance of fuzzy membership and typicality values in the objective function. A similar discussion has been done by Pal et al. [29]. For simplicity, we set CF = 1, CT = 1, m = 2 and  = 2 for all the experiments.

In this paper, the parameter ˛ is a constant, which controls the influence of the global intensity force and local intensity force. When the intensity inhomogeneity is severe, the bias estimation relies on the local intensity force. In such case, we should choose small ˛ as the weight of the global intensity force. Otherwise, the bias field estimation may perform poorly. For images with minor inhomogeneity, the accuracy of segmentation relies on the global intensity force. In this case, we can use relatively larger ˛ as the weight of global intensity. Thus, the global intensity reduces the misclassification for the pixels around the edges. In the experiment, we need to choose appropriate value for ˛ according to the degree of inhomogeneity. We suggest using ˛ = 0.5 for general MR images.

5.1. The radius r of neighborhood Nk Table 4 shows the JS and DC results of the segmentations with different window radius. The original images are brain MR images from BrainWeb [35]. From Table 4, we can observe that the quantitative accuracy of the segmentation results is decreasing when the noise level in brain MR images is increasing. The segmentations with r = 1 get the largest accuracy when the noise level is less than 9%. Large radius will lead to the loss of texture, so CSF can not be segmented accurately with the radius bigger than 1. On the other hand, the computational complexity of the algorithm will increase with the increasing of the radius. Generally, the typical noise levels corresponding to real MR data would be around 3%. Therefore, we suggest using r = 1 for general images such as MR images and real images.

5.4. The initialization robustness A salient advantage of our method is that its result is robust of initialization, which allows fully automated application. To visualize this robustness, we tested our model on 20 brain images from Brain Web [35] with 40% intensity inhomogeneity and 9% noise. 20 different initializations are used. The JS values for WM, GM and CSF of the results obtained by our model are shown in Fig. 14. The JS values of our method for the 20 different initializations show no visible difference in Fig. 14, which demonstrate the robustness of initialization of our method.

5.2. The parameters of Gaussian kernel The Gaussian kernel influences the bias field estimation. Fig. 13 shows the JS results of the segmentation with different parameter selections. The original image is obtained from Brain Web [35] with 40% intensity inhomogeneity and 9% noise. The left figure shows the influence of the standard deviation , while the right figure shows the influence of the radius . From Fig. 13, we can observe that the accuracy of segmentations increases with the increasing of  and . When  > 5 or  > 10, the JS values of GM decrease slightly, while the JS values of WM and CSF still increase. It is very hard to find the best parameter setting for WM, GM and CSF. Considering the segmentation accuracy of WM, GM and CSF comprehensively, we suggest that 4 ≤  ≤ 7 and 8 ≤  ≤ 15. In our experiments, we set  = 4 and  = 10 for general MR images.

Fig. 14. Test of robustness to initialization for our method. The x-axis represents 20 images, and the y-axis represents the JS values for WM (red), GM (green) and CSF (blue) for 20 different initializations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Z.-X. Ji et al. / Computerized Medical Imaging and Graphics 35 (2011) 383–397

6. Conclusions A modified possibilistic fuzzy c-means clustering algorithm is presented for fuzzy segmentation of MR images that have been corrupted by intensity inhomogeneities and noise. We propose an adaptive method to compute the weights for the neighborhood of each pixel in the image. The proposed adaptive method can not only overcome the effect of the noise effectively, but also prevent the edge from blurring. To address intensity inhomogeneity, the proposed algorithm introduces the global intensity into the CLIC algorithm and combines the local and global intensity information into account to ensure the smoothness of the derived optimal bias field and improve the accuracy of the segmentations. The proposed model can be treated as a general framework of FCM, PCM, CLIC and PFCM to segment the MR images with intensity inhomogeneities and noise. Our model can segment a brain MR image in 10–20 iterations within no more than 2 min. With good initialization, the model may need less iteration and can obtain results in less time. A variety of images, including synthetic images, synthetic brain MR images and real brain MR images are used to compare the performance of the proposed algorithm. Experimental results show that the proposed method is effective, more robust for different level of noise and more accuracy for both 2D and 3D brain MR image segmentation. Moreover, it is robust to initialization, thereby allowing fully automated applications. References [1] Ho BC. MRI brain volume abnormalities in young, nonpsychotic relatives of schizophrenia probands are associated with subsequent prodromal symptoms. Schizophrenia Research 2007;96(1):1–13. [2] Sikka K, Sinha N, Singh PK, Mishra AK. A fully automated algorithm under modified FCM framework for improved brain MR image segmentation. Magnetic Resonance Imaging 2009;27(7):994–1004. [3] Awate S, Tasdizen T, Foster N, Whitaker R. Adaptive Markov modeling for mutual-information-based, unsupervised MRI brain-tissue classification. Medical Image Analysis 2007;10(5):726–39. [4] Wong W, Chung A. Bayesian image segmentation using local isointensity structural orientation. IEEE Transactions on Image Processing 2005;14(10):1512–23. [5] Vovk U, Pernus F, Likar B. A review of methods for correction of intensity inhomogeneity in MRI. IEEE Transactions on Medical Imaging 2007;26(3):405–21. [6] Hou Z. A review on MR image intensity inhomogeneity correction. International Journal of Biomedical Imaging 2006:1–11. [7] Wells III WM, Grimson W, Kikinis R, Jolesz F. Adaptive segmentation of MRI data. IEEE Transactions on Medical Imaging 1996;15(4):429–42. [8] Van Leemput K, Maes F, Vandermeulen D, Suetens P. Automated model-based tissue classification of MR images of the brain. IEEE Transactions on Medical Imaging 1999;18(10):897–908. [9] Li C, Xu C, Anderson A, Gore J. MRI tissue classification and bias field estimation based on coherent local intensity clustering: a unified energy minimization framework. In: IPMI’09. Berlin, Heidelberg: Springer-Verlag; 2009. p. 288–99. [10] Pham D, Prince J. Adaptive fuzzy segmentation of magnetic resonance images. IEEE Transactions on Medical Imaging 1999;18(9):737–52.

397

[11] Ahmed MN, Yamany SM, Mohamed N, Farag AA. A modified fuzzy C-mean algorithm for bias field estimation and segmentation of MRI data. IEEE Transactions on Medical Imaging 2002;21(3):193–9. [12] Liew AWC, Yan H. An adaptive spatial fuzzy clustering algorithm for 3-D MR image segmentation. IEEE Transactions on Medical Imaging 2003;22(9):1063–75. [13] Saha PK, Udupa JK, Odhner D. Scale-based fuzzy connected image segmentation: theory, algorithms, and validation. Computer Vision and Image Understanding 2000;77(2):145–74. [14] Madabhushi A, Udupa JK, Souza A. Generalized scale: theory, algorithms, and application to image inhomogeneity correction. Computer Vision and Image Understanding 2006;101:100–21. [15] Nyul L, Udupa JK, Saha P. Incorporating a measure of local scale in voxel-based 3-D image registration. IEEE Transactions on Medical Imaging 2003;22(2):228–37. [16] Souza A, Udupa JK, Madabushi A. Image filtering via generalized scale. Medical Image Analysis 2007;12(2):87–98. [17] Zhuge Y, Udupa JK, Liu L, Saha PK. A scale-based method for correcting background intensity variation in acquired images. Proceedings of SPIE Medical Imaging 2002;4684:1103–11. [18] Zhuge Y, Udupa JK. Intensity Standardization simplifies brain MR image segmentation. Computer Vision and Image Understanding 2009;113(10):1095–103. [19] Cai W, Chen S, Zhang DQ. Fast and robust fuzzy c-means algorithms incorporating local information for image segmentation. Pattern Recognition 2007;40(3):825–38. [20] Szilágyi L, Szilágyi SM, Benyó Z. A modified FCM algorithm for fast segmentation of brain MR images. In: ICIAR LNCS 2007, vol. 4633. 2007. p. 866–77. [21] Kang JY, Min LQ, Luan QX, Li X, Liu JZ. Novel modified fuzzy c-means algorithm with applications. Digital Signal Processing 2009;19(2):309–19. [22] Kang BY, Kim DW, Li Q. Spatial homogeneity-based fuzzy c-means algorithm for image segmentation. In: FSKD LNAI 2005, vol. 3613. 2005. p. 462–9. [23] Wang JZ, Kong J, Lu YH, Qi M, Zhang BX. A modified FCM algorithm for MRI brain image segmentation using both local and non-local spatial constrains. Computerized Medical Imaging and Graphics 2008;31(8):685–98. [24] Krishnapuram R, Keller JM. Possibilistic approach to clustering. IEEE Transactions on Fuzzy Systems 1993;1(2):98–100. [25] Barni M, Cappellini V, Mecocci A. Comments on a possibilistic approach to clustering. IEEE Transactions on Fuzzy Systems 1996;4(3):393–6. [26] Krishnapuram R, Keller JM. The possibilistic c-means algorithm: insights and recommendations. IEEE Transactions on Fuzzy Systems 1996;4(3):383–95. [27] Timm H, Borgelt C, Doring C, Kruse R. An extension to possibilistic fuzzy cluster analysis. Fuzzy Sets and Systems 2004;147(1):3–16. [28] Zhang JS, Leung YW. Improved possibilistic c-means clustering algorithms. IEEE Transactions on Fuzzy Systems 2004;12(2):209–17. [29] Pal NR, Pal K, Keller JM, Bezdek JC. A possibilistic fuzzy c-means clustering algorithm. IEEE Transactions on Fuzzy Systems 2005;13(4):517–30. [30] Ji ZX, Chen Q, Sun QS, Xia DS, Heng PA. MR image segmentation and bias field estimation using coherent local and global intensity clustering. In: FSKD 2010, proceedings 2010 seventh international conference on fuzzy systems and knowledge discovery, vol. 2. 2010. p. 578–82. [31] Dunn JC. A fuzzy relative of the ISODATA process and its use in detecting compact well separated cluster. Journal of Cybernet 1974;3(3):32–57. [32] Bezdek JC. Pattern recognition with fuzzy objective function algorithms. Norwell, MA, USA: Kluwer Academic Publishers; 1981. [33] Bezdek JC, Hall LO, Clarke LP. Review of MR image segmentation techniques using pattern recognition. Journal of Cybernet 1999;20(4):1033–48. [34] http://www.cma.mgh.harvard.edu/ibsr/. [35] http://www.bic.mni.mcgill.ca/brainweb/. [36] http://brainmaps.org/index.php.