Advanced Segmentation Techniques PDF

Chapter 9 Advanced Segmentation Techniques Aly A. Farag Mohamed N.AhmedAyman El-Baz and Hossam Hassan 9.1 Introduction The principal goal of the segmentation process is to partition an image into regions that are homogeneous with respect to one or more characteristics or features. Segmentation...

Chapter 9 Advanced Segmentation Techniques Aly A. Farag Mohamed N.AhmedAyman El-Baz and Hossam Hassan 9.1 Introduction The principal goal of the segmentation process is to partition an image into regions that are homogeneous with respect to one or more characteristics or features. Segmentation is an important tool in medical image processing and it has been useful in many applications including lesion quatification, surgery simulations,surgical planning, multiple scleroris,functional mapping,computer assisted diagnosis, image registration and matching, etc. A wide varity of segmentation techniques has been proposed.However,there is no one standard segmentation technique that can produce satisfactory results for all imaging applications.Quite often,methods are optimized to deal with spe cific imaging modalities such as magnetic resonance (MR) imaging and X-ray computed tomography CT),or modeled to segment specific anatomic struc- tures such as the brain, the lungs and the vascular system. Recent research has demonstrated that the segmentation of anatomical structures from MRI and CT will benefit from the exploitation of three different types of knowledge: intensity models that describe the gray-level appearance of individual structures, shape models that descibe the shape of different struc- tures as well as imaging models that capture the characteristics of the imaging process. 1Computer Vision and Image Processing Laboratory, Department of Electrical and Com- puter EngineeringUniversity of Louisville LouisvilleKY 40292USA 2Software ResearchC19LLexmark International Inc.LexingtonKY 40550USAE-mail: [email protected] 479 480 Farag,Ahmed, El-Baz and Hassan Stochastic image models are useful in quantitatively specifying natural con- straints and general assumption about the physical world and the imaging pro- cess. Random field models permit the introduction of spatial context into pixel labeling problem. An introduction to random fields and its application in lung CT segmentation will be presented in Section 9.2. Crisp segmentation by which a pixel is assigned to a one particular region, often presents problems. In many situations it is not easy to determine if a pixel should belong to a region or not. This is because the features used to determine homogeneity may not have sharp transitions at region boundaries.To alleviate this situation, we can inset fuzzy set concepts into the segmentation process. In Section 9.4,we will present an algorithm for fuzzy segmentation of MRI data and estimation of intensity inhomogeneities using fuzzy logic. MRI intensity inhomogeneities can be attributedto imperfectionsinthe RF coils ortoproblems associated with the acquisition sequences.The result is a slowly varying shading artifact over the image that can produce errors with conventional intensity- based classification. The algorithm is formulated by modifying the objective function of the standard fuzzy c-means (FCM) algorithm to compensate for such inhomogeneities and to allow the labeling of a pixel voxel) to be influenced by the labels in its immediate neighborhood.The neighborhood effect acts as a regularizer and biases the solution toward piecewise-homogeneous labelings. Such aregularization is useful in segmenting scans corrupted by salt and pepper noise. Section 9.5 is devoted to the description of geometrical methods and their application in image segmentation.Among many methods used for shape recov- ery,the level sets has proven to be a successful tool.The level set is a method for capturing moving fronts introduced by Osher and Sethian in 1987.It was used in many applications like fluid dynamics,graphics, visualization, image process- ing and computer vision. In this chapter, we introduce an overview of the level set and its use in image segmentation with application in vascular segmentation. The human cerebrovascular system is a complex three-dimensional anatomical structure.Serious types of vascular diseases such as carotid stenosis, aneurysm, and vascular malformation may lead to brain stroke, which are the third leading cause of death and the main cause of disability. An accurate model of the vas- cular system from MRA data volume is needed to detect these diseases at early stages and hence may prevent invasive treatments.In this section, we will use Advanced Segmentation Techniques 481 a method based on level sets and statistical models to improve the accuracy of the vascular segmentation. 9.2Stochastic Image Models The objective of modeling in image analysis is to capture the intrinsic character ofimages in a few parameters so as to understand the nature of the phenomenon generating the images. Image models are also useful in quantitatively specifying natural constraints and general assumptions about the physical world and the imaging process.The introduction of stochastic models in image analysis has led to the development of many practical algorithms that would not have been realized with ad hoc processing.Approaching problems in image analysis from the modeling viewpoint, we focus on the key issuesof model selection,sampling, parameter estimation, and goodness-of-fit. Formal mathematical image models have long been used in the design of image algorithms for applications such as compression, restoration, and en- hancement. Such models are traditionally low stochastic models of limited complexity. In recent years, however, important theoretical advances and in- creasingly powerful computers have led to more complex and sophisticated image models.Depending on the application, researchers have proposed both low-level and high-level models. Low-level image models describe the behavior of individual image pixels rel- ative to one another. Markov random fields and other spatial interaction models have proven useful for a variety of applications, including image segmentation and restoration [2,3].Bouman et al. ,along with Willsky and Benvensite [5,6], have developed multiscale stochastic models for image data. High-level models are generally used to describe a more restrictive class of images.These models explicitly describe larger structures in the image, rather than describing individual pixel interactions.Grenander et al., for example,pro- pose a model based on deformable templates to describe images of nonrigid ob- jects , while Kopec and his colleagues model document images using a Markov source model for symbol generation in conjunction with a noisy channel [8,9]. The following part of this chapter is organized as follows: First a short introduction about Gibbs random field (GRF) and Markov random field (MRF) 482 Farag,Ahmed, El-Baz and Hassan is given.A detailed description of our proposed approach to get an accurate image model is then presented.Finally, we will apply the proposed model in the segmentation of lung CT. 9.2.1Statistical Framework The observed image is assumed to be a composites of two random process: a high-level process X, which represents the classes that form the observed image; and a low-level process Y,which describes the statistical characteristics of each class. The high-level process X is a random field defined on a rectangular grid S of N2 points and the value of X will be written as Xs.Points in X will take values in the set (T1.Im), where m is the number of classes in the given image. Given x, the conditional density function of y is assumed to exist and to be strictly positive and is denoted by pY=y|X=x or py|x. Finally,an image is a square grid of pixels, or sites,{iji=1 to Nj= 1 to N}. We adopt a simple numbering of sites by assigning sequence number t=j+Ni-1) to site s.This scheme numbers the sites row by row from 1 to N2,starting in the upper left. 9.2.2Gibbs Random Fields In 1987,Boltzmann investigated the distribution of energy states in molecules of an ideal gas.According to the Boltzmann distribution the probability of a molecule to be in a state with energyis 1 ps= 7 9.1) where Z is a normalization constant,that makes the sum of probabilities equal to 1.T is the absolute temperature, and K is Boltzmann's constant.For simplicity we assume that the temperature is measured in energy units, hence K T will be replaced by T. Gibbs used a similar distribution in 1901 to express the probability of a whole system with many degrees of freedom to be in a state with a certain energy.A discrete GRF provides a global model for an image by specifying a probability Advanced Segmentation Techniques 483 mass function in the following form px= Z 9.2) where Z=xeand the function Ex is called energy function. 9.2.3Markov Random Fields Hassner and Sklansky introduced Markov random fields to image analysis and throughout the last decade Markov random fields have been used extensively as representations of visual phenomena.A Gibbs random filed describes the global properties of an image in terms of the joint distributions of colors for all pixels.An MRF is defined in terms of local properties.Before we show the basic properties of MRF, we will show some definitions related to Gibbs and Markov random fields [10-15]. Definition 1: A clique A is a subset of S for which every pair of sites is a neighbor. Single pixels are also considered cliques.The set of all cliques on a grid is called A. Definition 2: A random field X is an MRF with respect to the neighborhood system n={nss ES} if and only if p(X=x>0 for all xwhereis the set ofall possible configurations on the given grid pX=x|Xs=xs=pX=xXs=xs)wheresr refers to allN2 sites excluding site r,and as refer to the neighborhood of site s; pX=xsXas=xas is the same for all sites s. The structure of the neighborhood system determines the order of the MRF. For a first-order MRF the neighborhood of a pixel consists of its four nearest neighbors.In asecond-order MRF the neighborhood consists ofthe eight nearest neighbors. The cliques structure are illustrated in Figs 9.1 and 9.2. Consider a graph tn) as shown in Fig.9.3 having a set of N sites. The energy function for a pairwise interaction model can be written in the following form: Ex=Fx+Hxxr, t=1r=1 9.3 484 Farag,Ahmed, El-Baz and Hassan Figure 9.1: Cliques for a first-order neighborhood, where d1,and are the cliques coefficients for first-order neighborhood system. Figure 9.2: Cliques for a second-order neighborhood, where 01.09 are the cliques coefficients for second-order neighborhood system. t-11 t:-7 t-6 t+8 t:+12 5 4 3 4 5 t-9 t:-3 t-2 t+4 t+10 4 2 2 4 t-5 t-1 t t+1 t+5 3 t 1 3 t-10 t+3 t:-4 t:+2 t+9 4 2 2 4 t:-12 t:-8 t+6 t+7 t+11 5 4 3 4 S (a) (b) Figure 9.3: Numbering and order coding of the neighborhood structure. where F is the potential function for single-pixel cliques and H is the poten tial function for all cliques of size 2.The parameter w depends on the size of the neighborhood around each site.For example,w is 24,6,10, and 12 for neighborhoods of orders 1,2,3,4,5,respectively. Using the Derin-Elliott model to compute F and H, we have Fx=axandHxxt+r=0Ixxt+r, where I(a,b) is called indicator function where Iab=-1ifa=b =lifa b. 9.2.4Image Models As mentioned before, the observed image is modeled as a composite of two random processes, a high-level process X and a low-level process Y [16-20] Advanced Segmentation Techniques 485 The maximum aposterioriparameters estimation involves the determination of x that maximizes p(x|y) with respect to x.By Bayes'rule, 9.4 p(y) Since the denominator of Eq.9.4 does not affect the optimization the MAP parameters estimation can be obtained, equivalently, by maximizing the numer- ator ofEq.9.4 or its natural logarithm;that is,we need to find x which maximizes the following criterion: Lxy=lnpyx+ln px. 9.5) The first term in Eq.9.5 is the likelihood due to the low-level process and the second term is due to the high-level process.Based on the models of the high- level and low-level processes, the MAP estimate can be obtained. In order to carry out the MAP parameters estimation in Eq.9.5, one needs to specify the parameters of the two processes.A popular model for the high-level process is the Gibbs Markov model. In the following sections we introduce a new accurate model to model the low-level process. In this model we will assume that each class consists of a mixture of normal distributions which follow the following equation: py|= pyC for i=1,2,...m (9.6) l=1 where n is the number of normal components that formed class T is the corresponding mixing proportion, and {Cis the number of Gaussian com- ponents that formed class I'.So the overall model for the low-level process can be expressed as follows: Pesy=pTpyi 9.7) i=1 In our proposed algorithm the priori probability p ) is included in the mixing proportion for each class. 9.2.5Parameter Estimation for Low-Level Process In order to estimate the parameters for low-level process, we need to esti- mate the number of Gaussian components that formed the distribution for each class, their means, the variances, and mixing proportions for each Gaussian 486 Farag,Ahmed, El-Baz and Hassan component. To estimate the distribution for each class, we use the expecta- tion maximization algorithm.The first step to estimate the distribution for each class is to estimate the dominant Gaussian components in the given empirical distribution. 9.2.5.1 Dominant Gaussian Components Extracting Algorithm 1.Assume the number of Gaussian components that represent the classes Ti= 1,..m.Initialize the parameters of each distribution randomly 2.The E-stepCompute that represent responsibility that the given pixel value is extracted from certain distribution as py1 fort=1toN, 9.8) pyr where yt is the gray level at location t in the given image, is the mix ing proportion of Gaussian component i at step k and is estimated parameter for Gaussian component i at step k. 3.The M-step:we compute the new mean the new variance,and the new proportion from the following equations: N2 J8it, 9.9) t=1 9.10 t+y=n-) 9.11) 4.Repeat steps 1 and 2 until the relative difference of the subsequent values of Eqs.9.9,9.10and 9.11 are sufficiently small. Let pir1yP1 2yP1my be the dominant Gaussian components that are estimated from the above algorithm. Then the initial estimated density p1y) for the given image can be defined as follows: p1y=pI1y+Pr2y++mPImy. 9.12) Because the empirical data does not exactly follow mixture of normal distri- bution there will be error between pr(y and pem(y.So we suggestthe following Advanced Segmentation Techniques 487 models for the empirical data: Pemy=p1y+y) (9.13) where y represent the error between pemyand pry.From Eq.9.13,y can be rewritten as y=|Pemy-p1ysignPemy-p1y 9.14) We assume that the absolute value of (y) is another density which consists of a mixture of normal distributions and we will use the following EM algorithm to estimate the number of Gaussian components in y and the mean, the variance, and mixing proportion. 9.2.5.2 Sequential EM Algorithm 1.Assume the number of Gaussian components (n) in (y is 2. 2.The E-stepGiven the current value of the number of Gaussian compo nents in ycompute as py fori 1ton andt=1toN. 9.15) py 3.The M-step:We compute the new mean, the new variance, and the new proportion from the following equations: N2 J (9.16) 9.17) (9.18) 4. Repeat steps 1 and 2 until the relative differences ofthe subsequent values of Eqs.9.16,9.17, and 9.18 are sufficiently small. 5.Compute the conditional expectation and the error between |(y)| and the estimated density py)for |y|from the following equations: Qn=uln pyi, (9.19) t=1 i=1 n=5y|- piy 9.20) 488 Farag,Ahmed, El-Baz and Hassan 6.Repeat steps 2,3,4,and 5, and increase the number of Gaussian com- ponents n by 1 if the conditional expectation Qn) is still increasing and n is still decreasing, otherwise stop and select the parameters which correspond to maximum Q(n) and minimum (n. Since EM algorithm can be trapped in a local minimum, we run the above algorithm several times and select the number of Gaussian components and their parameters that give maximum Q(n and minimum e(n. After we determined the number ofGaussian components that formed |(y)l, we need to determine which components belong to class T1,and belong to class T2,and so on. In this model we classify these components based on the mini- mization of risk function under 0-1 loss. In order to minimize the risk function, we can use the following algorithm.Note that the following algorithm is writen for two classes but it is easy to generalize to n classes. 9.2.5.3Components Classification Algorithm 1. All Gaussian components that have mean less than the estimated mean for PIr1(y) belong to the first class. 2.All Gaussian components that have mean greater than the estimated mean for pir(y belong to the second class. 3.For the components which have mean greater than the estimated mean for Pry and less than the estimated mean for pr(y,do the following a Assume that the first component belongs to the first class and the other components belong to the second class. Compute the risk value from the following equation TI RTh= JTh p(y|l1dy+ py|2dy, 9.21) where Th is the threshold that separates class from class 2.The above integration can be done using a second-order spline. b) Assume that the first and second components belong to the first class and other components belong to the second class, and from Eq.9.21 compute R(Th).Continue this process as R(Th) decreases, and stop when R(Th) starts to increase. Advanced Segmentation Techniques 489 Finally, to show the convergence of the proposed model, we will show ex- perimentally, when we use this model, the Levy distance will decrease between the estimated distribution Pes(y) and empirical distribution Pem(y).The Levy distancePemPes) is defined as PemPes=inf{>0VyPemy--PesyPemy++}. 9.22 As pPemPes approach zeroPesy converge weakly to Pemy 9.2.6 Parameter Estimation for High-Level Process In order to carry out the MAP parameters estimation in Eq.9.5, one needs to specify the parameters of high-level process.A popular model for the high-level process is the Gibbs Markov model which follows Eq.9.2. In order to estimate the parameters of GMRF, we will find the parameters that maximize Eq.9.2, and we will use the Metropolis algorithm and genetic algorithm (GA). The Metropolis algorithm is a relaxation algorithm to find a global maximum. The algorithm assumes thatthe classes of all neighbors of yare known.The high- level process is assumed to be formed of m-independent processes; each of the m processes is modeled by Gibbs Markov random which follow Eq.9.2.Then y can be classified using the fact that pxyis proportional to py|xPxns, where s is the neighbor set to site S belonging to class xtpxtns) is computed from Eq.9.2,and p(y x is computed from the estimated density for each class. By using the Bayes classifier, we get initial labeling image.In order to run the Metropolis algorithm, first we must know the coefficients of potential function Ex),so we will use GA to estimate the coefficient of E(x and evaluate these coefficients through the fitness function. 9.2.6.1 Maximization Using Genetic Algorithm To build the genetic algorithm, we define the following parameters Chromosome:A chromosome is represented in binary digits and consists of representations for model order and clique coefficients.Each chromosome has 41 bits.The first bit represent the order of the system (we use digit 0 for first- order and digit 1for second-order-GMRF).The remaining bits represent the 490 Farag,Ahmed, El-Baz and Hassan clique coefficients, where each clique coefficient is represented by 4 bits (note that for first-order system, we estimate only five parameters, and the remaining clique's coefficient will be zero,but for the second-order system we will estimate ten parameters). Fitness Function:Since our goal is to select the high-level process X that maximize Eq.9.5, we can use Eq.9.5 as the fitness function. High-level parameters estimation algorithm: 1.Generate the first generation which consists of 30 chromosomes. 2. Apply the Metropolis algorithm for each chromosome on each image and then compute the fitness function as shown in Eq.9.5. 3. If the fitness values for all chromosomes do not change from one popula- tion to another population, then stop and select the chromosome, which gives maximum fitness value.(If there are two chromosomes that give the same fitness value,we select the chromosome which represents lower order system.) Otherwise go to step2. Using the results obtained by this algorithm, we will repeat the estimation of low-level process and high-level process.We will stop when the difference between the current parameters and previous parameters is small. 9.3Applications Lung Cancer remains the leading cause of mortality cancer.In 1999, there were approximately 170 000 new cases of lung cancer. The 5-year survival rate from the diseases is 14% and has increased only slightly since the early 1970s despite extensive and expensive research work to find effective therapy. The disparity in survival between early and late-stage lung cancer is substantial, with a 5-year survival rate of approximately 70% in stage 1A disease compared to less than 5% in stage IV disease according to the recently revised lung cancer staging criteria. The disproportionately high prevalence of and mortality from lung cancer has encouraged attempts to detect early lung cancer with screening programs aimed at smokers.Smokers have an incidence rate of lung Advanced Segmentation Techniques 491 cancer that is ten times that of nonsmokers and accounts for greater than 80% of lung cancer cases in the United States. One in every 18 women and every 12 men develop lung cancer, making it the leading cause of cancer deaths. Early detection of lung tumors (visible on the chest film as nodules may increase the patient's chance of survival. For this reason the Jewish Hospital designed a program for early detection with the following specific aims:A number of lung cancer screening trials have been conducted in the United States,Japan, and Europe for the purpose of developing an automatic approach of tummor detection. At the University of Louisville CVIP Lab a long-term effort has been ensued to develop a comprehensive image analysis system to detect and recognize lung nodules in low dose chest CT (LDCT) scans.The LDCT scanning was performed with the following parameters:slice thickness of 8 mm reconstructed every 4 mm and scanning pitch of 1.5.In the following section we highlight our approach for automatic detection and recognition of lung nodules; further details can be found in. 9.3.1Lung Extraction The goal of lung extractionis to separate the voxels corresponding to lung tissue from those belonging to the surrounding anatomical structures. We assume that each slice consists of two types of pixels: lung and other tissues (e.g.chest, ribs, and liver).The problem in lung segmentation is that there are some tissues in the lung such as arteries, veins, bronchi and bronchioles having gray level close to the gray level of the chest.Therefore,in this application if we depend only on the gray level we lose some of the lung tissues during the segmentation process.Our proposed model which depends on estimating parameters for two processes (high-level process and low-level process) is suitable for this appli- cation because the proposed model not only depend on the gray level but also takes into consideration the characterization of spatial clustering of pixels into regions. We will apply the approach that was described in Section 9.2.4 on lung CT. Figure 9.4 shows a typical CT slice for the chest.We assume that each slice consists of two types of tissues:lung and other tissues e.g.chest, ribs, and liver).As discussed above, we need to estimate parameters for both low-level process and high-level process.Table 9.1 presents the results of applying the 492 Farag,Ahmed, El-Baz and Hassan Table 9.1: Estimated using dominant Gaussian components extracting algorithm Parameter 1 112 111 12 Value 59.29 139.97 177.15 344.29 0.25 0.758 dominant Gaussian components extractingalgorithm describedin9.2.5.1.Figure 9.5 shows the empirical density for the CT slice shown in Fig.9.4 and the initial estimated density (which represented the two dominant Gaussian components in the given CT). The Levy distance between the two distribution functions which represented the densities shown in Fig.9.5 is 0.09.This value is large and this means there is a mismatch between empirical pem(y and pi(y.Figure 9.6 shows the error and absolute error between pem(y) and p1y After we apply sequential EM algorithm to |(yl,we get that the number of normal components that represent y) is 10 as shown in Fig.9.7.Figure 9.8 Figure 9.4:A typical slice form of a chest spiral CT scan. Advanced Segmentation Techniques 493 Figure 9.5:Empirical density for given CT slice and initial estimated density. Figure 9.6:Error and absolute error between PemY yand piY=y 494 Farag,Ahmed, El-Baz and Hassan Figure 9.7:Conditional expectation Q(n and the error function en)ver sus the number of Gaussians approximating the scaled absolute deviation in Fig.9.6. Figure 9.8:Estimated density for nY=y. Advanced Segmentation Techniques 495 Figure 9.9:12 Gaussian components which are used in density estimation. shows the estimated density for |y.Figure 9.9 shows all Gaussian compo nents which are estimated after using dominant Gaussian components extract- ing algorithm and sequential EM algorithms. Figure 9.10 shows the estimated density for the CT slices shown in Figure 9.4.The Levy distance between the distributions Pes(y and Pem(y) is 0.0021 which is smaller compared to the Levy distance between the distributions Pem(y and PI(y) Now we apply components classification algorithm on the ten Gaussian com- ponents that are estimated using sequential EM algorithm in order to determine which components belong to lung tissues and which components belong to chest tissues.The results of components classification algorithm show that the minimum risk equal to 0.00448 occurs at threshold Th108 when Gaussian components 1,2,3,and 4 belong to lung tissues and component 5,6,7,8,9, and 10 belong to chest tissues.Figure 9.11 shows the estimated density for lung tissues and estimated density for chest and other tissues that may appear in CT. Thenext step ofour algorithm istoestimatethe parametersforhigh-level pro- cess.A popular model for the high-level process is the Gibbs Markov mode, and we use the Bayes classifier to get initial labeling image.After we run Metropolis algorithm and GA to determine the coefficients of potential function E(x),we get 496 Farag,Ahmed, El-Baz and Hassan Figure 9.10:Estimated density for lung tissues and chest tissues. Figure 9.11:Empirical density and estimated density for CT slice shown in Fig.9.4. Advanced Segmentation Techniques 497 Figure 9.12(a Segmented lung using the proposed algorithm,error 1.09%. b) Output of segmentation algorithm by selecting parameters for high-level process randomly, error =1.86%.(c Segmented lung by radiologist. the following results:=10=0.890=0.803=0.7804=0.6905=0.54, 06=0.6107=0.890s=0.56and0g=0.99. The result of segmentation for the image shown in Fig.9.4 using these pa- rameters is shown in Fig.9.12.Figure 9.12(a) shows the results of proposed algo rithm.Figure 9.12(b) shows output of the Metropolis algorithm by selecting pa rameters randomly.Figure 9.12(c) shows the segmentation done by a radiologist. As shown in Fig.9.12(a) the accuracy of our algorithm seems good if it is compared with the segmentation of the radiologist.Figure 9.13 shows compari- son between our results and the results obtained by iterative threshold method which was proposed by Hu and Hoffman.It is clear from Fig.9.13 that the Figure 9.13:(aOriginal CT,(b) segmented lung using the proposed model,(c) segmented lung using the iterative threshold method, and (d) segmented lung by radiologist.The errors with respect to this ground truth are highlighted by red color. 498 FaragAhmed,El-Baz,and Hassan Figure 9.14:a Generated Phantom, b ground truth image (black pixel rep- resent lung area,and gray pixels represent the chest area) and (c) segmented lung using the proposed approach (error 0.091).The errors with respect to this ground truth are highlighted by red color. proposed algorithm segments the lung without causing any loss of abnormality tissues if it is compared with the iterative threshold method. Also, in order to validate our results we create aphantom which has the same distribution as lung and chest tissues.This phantom is shown in Fig.9.14.One of the advantages of this phantom is that we know its ground truth.It is clear from Fig.9.14 that the error between segmented lung and ground truth is small and this shows that the proposed model is accurate and suitable for this application. 9.4FuzzySegmentation As mentioned before, the objective of image segmentation is to divide an image into meaningful regions. Errors made at this stage would affect all higher level activities.Therefore, methods that incorporate the uncertainty of object and region definitionsand the faithfulness ofthe features torepresent variousobjects are desirable. In an ideally segmented image, each region should be homogeneous with respect to some predicate such as gray level or texture, and adjacent regions should have significantly different characteristics or features.More formally, segmentation is the process of partitioning the entire image into c crisp maxi- mally connected regions Ri} such that each R is homogeneous with respect to some criteria. In many situations, it is not easy to determine if a pixel should belong to a region or not. This is because the features used to determine homo- geneity may not have sharp transitions at region boundaries. To alleviate this situation, we can inset fuzzy set concepts into the segmentation process. Advanced Segmentation Techniques 499 In fuzzy segmentation, each pixel is assigned a membership value in each of the c regions. If the memberships are taken into account while computing properties of regions, we oftain obtain more accurate estimates of region prop erties.One of the known techniques to obtain such a classification is the FCM algorithm [40,41].The FCM algorithm is an unsupervised technique that clus- ters data by iteratively computing a fuzzy membership function and mean value estimates for each class.The fuzzy membership function,constrained to be be tween O and 1,reflects the degree of similarity between the data value at that location and the prototypical data value, or centroid, ot its class. Thus, a high membership value near unity signifies that the data value at that location is close to the centroid of that particular class. FCM has been used with some success in image segmentation in general [45,46], however, since it is a point operation, it does not preserve connectivity among regions.Furthermore,FCM is highly sensitive to noise.In the following sections, we will present anew systemto segment digital imagesusingamodified Fuzzy c-means algorithm.Our algorithm is formulated by modifying the objec tive function of the standard FCM algorithm to allow the labeling of a pixel to be influenced by the labels in its immediate neighborhood. The neighborhood ef fect acts as a regularizer and biases the solution toward piecewise-homogeneous labelings.Such a regularization is useful in segmenting scans corrupted by scan- ner noise. In this paper, we will present the results of applying this algorithm to segment MRI data corrupted with a multiplicative gain field and salt and pepper noise. 9.4.1Standard Fuzzy-C-Means The standard FCM objective function for partitioning =1into c clusters is given by J=ull-vill 9.23) where are the feature vectors for each pixel, [i are the prototypes of the clusters and the array [u]=U represents a partition matrix,U U,namely u{uik[0,1uk=1 k i=1 500 Farag,Ahmed, El-Baz and Hassan and Uik1, [8Fm P. =0. 9.28) Suik NR ux=u Solving for uwe have oDik (9.29) Sincej=1j=1 Vk, Yi (9.30) or (9.31) Substituting into Eq.9.29, the zero-gradient condition for the membership esti mator can be rewritten as 1 (9.32) Dik+ Dj+Yj 502 FaragAhmed,El-Baz and Hassan 9.4.3.2Cluster Prototype Updating Using the standard Eucledian distance and taking the derivative of Fm w.r.t.vi and setting the result to zero, we have N x-vi =0. 9.33) k-1 yEN v=v Solving for viwe have x+xNx 9.34) 1+a 9.4.4Application:Adaptive MRI Segmentation In this section, we describe the application of the MFCM segmentation on MRI images having intensity inhomogeneity.Spatial intensity inhomogeneity induced by the radio frequency (RF) coil in magnetic resonance imaging (MRI) is amajor problem in the computer analysis of MRI data [24-27].Such inhomogeneities have rendered conventional intensity-based classification of MR images very difficult, even with advanced techniques such as nonparametric, multichannel methods [28-30].This is due to the fact that the intensity inhomogeneities ap- pearing in MR images produce spatial changes in tissue statistics, i.e. mean and variance. In addition, the degradation on the images obstructs the physician's diagnoses because the physician has to ignore the inhomogeneity artifact in the corrupted images The removal of the spatial intensity inhomogeneity from MR images is diffi- cult because the inhomogeneities could change with different MRI acquisition parameters from patient to patient and from slice to slice.Therefore,the correc- tion of intensity inhomogeneities is usually required for each new image.In the last decade, a number of algorithms have been proposed for the intensity inho- mogeneity correction.Meyer et al. presented an edge-based segmentation scheme to find uniform regions in the image followed by a polynomial surface fit to those regions. The result of their correction is, however, very dependent on the quality of the segmentation step. Several authors have reported methods based on the use of phantoms for intensity calibration.Wicks et al. proposed methods based on the signal Advanced Segmentation Techniques 503 produced by a uniform phantom to correct for MRI images of any orienta- tion. Similarly,Tincher et al. modeled the inhomogeneity function by a second-order polynomial and fitted it to a uniform phantom-scanned MR image. These phantom approaches, however, have the drawback that the geometry relationship of the coils and the image data is typically not available with the image data. They also require the same acquisition parameters for the phan- tom scan and the patient. In addition, these approaches assume the intensity corruption effects are the same for different patients, which is not valid in general. The homomorphic filtering approach to remove the multiplicative effect of the inhomogeneity has been commonly used due to its easy and efficient im- plementation [29,34]. This method, however, is effective only on images with relatively low contrast. Some researchers [33,35] reported undesirable artifacts with this approach. Dawant et al. used operator-selected reference points in the image to guide the construction ofathin-plate spline correctionsurface.The performance of this method depends substantially on the labeling of the reference points. Considerable user interactions are usually required to obtain good correction results. More recently, Gilles et al. proposed an automatic and iterative B- spline fitting algorithm for the intensity inhomogeneity correction of breast MR images.The application ofthis algorithm is restricted to MR images with a single dominant tissue class, such as breast MR images.Another polynomial surface fitting method was proposed based on the assumption that the number of tissue classes, the true means, and standard deviations of all the tissue classes in the image are given.Unfortunately,the required statistical information is usually not available. A different approach used to segment images with intensity inhomogeneities is to simultaneously compensate for the shading effect while segmenting the image.This approach has the advantage of being able to use intermediate infor- mation from the segmentation while performing the correction.Recently, Wells et al. developed a new statistical approach based on the EM algorithm to solve the bias field correction problem and the tissue classification problem. Guillemaud et al. further refined this technique by introducing the extra class "other."There are two main disadvantages of this EM approach.First, the EM algorithm is extremely computationally intensive, especially for large 504 FaragAhmedEl-Bazand Hassan problems.Second, the EM algorithm requires a good initial guess for either the bias field or the classification estimate.Otherwise, the EM algorithm could be easily trapped in a local minimum, resulting in an unsatisfactory solution. Another approach based on the FCM [40,41] clustering technique has been introduced lately [42-44].FCM has been used with some success in image seg- mentation in segmenting MR images [42,47,50].Xu et al. proposed a new adaptive FCM technique to produce fuzzy segmentation while compensating for intensity inhomogeneities. Their method, however, is also computationally intensive. They reduced the computational complexity by iterating on a coarse grid rather than the fine grid containing the image. This introduced some er- rors in the classification results and was found to be sensitive to a considerable amount of salt and pepper noise To solve the problem of noise sensitivity and computational complexity of the Pham and Prince method, we will generalize the MFCM algorithm to segment MRI data in the presence of intensity inhomogeneities. 9.4.4.1 Signal Modeling The observed MRI signal is modeled as a product of the true signal generated by the underlying anatomy and a spatially varying factor called the gain field: Y=XGK VkE [1,N] 9.35) where X and Y are the true and observed intensities at the kth voxel,respec tively, G is the gain field at the kth voxel, and N is the total number of voxels in the MRI volume. The application of a logarithmic transformation to the intensities allows the artifact to be modeled as an additive bias field yk=Xk+Bk VkE[1N] (9.36) where x and yk are the true and observed log-transformed intensities at the kth voxel,respectively, and is the bias field at the kth voxel.If the gain field is known,it is relatively easy to estimate the tissue class by applying a conventional intensity-based segmenter to the corrected data.Similarly, if the tissue classes are known, we can estimate the gain field,but it may be problematic to estimate Advanced Segmentation Techniques 505 either without the knowledge ofthe other.We will show that by using an iterative algorithm based on fuzzy logic, we can estimate both. 9.4.4.2Bias Corrected Fuzzy C-means (BCFCM) Objective Function Substituting Eq.9.36 into Eq.9.25, we have N N y-Br i=1k=1 yN 9.37) Formally, the optimization problem comes in the form min Jm subject to UEU. (9.38) U. 9.4.4.3 BCFCM Parameter Estimation The objective function Jm can be minimized in a fashion similar to the MFCM algorithm.Taking the first derivatives of Jm with respect to uikViand k and setting them to zero results in three necessary but not sufficient conditions for Jm to be at a local extrema.In the following subsections, we will derive these three conditions. 9.4.4.4 Membership Evaluation Similar to the MFCM algorithm, the constrained optimization in Eq.9.38 will be solved using one Lagrange multiplier Fm= 9.39) where Di=1-B-vand =|--v.The zero gradient condition for the membership estimator can be written as 1 uik (9.40) Dik+ ik+ 506 Farag,Ahmed, El-Baz and Hassan 9.4.4.5Cluster Prototype Updating Taking the derivative of Fm w.r.t.v and setting the result to zero, we have N N =0. 9.41 k=1 Solving for vi,we have -B+Ny-B 9.42) 1+= 9.4.4.6Bias Field Estimation In a similar fashion taking the derivative of Fm w.r.t.k and setting the result to zero we have Ca 21 0 9.43) =B Since only the kth term in the second summation depends on k,we have uyk-Bk-vi2 =0. 9.44) =P Differentiating the distance expression, we obtain yu-Bu i- 0. 9.45 i=1 i=1 =P Thus, the zero-gradient condition for the bias field estimator is expressed as B*=yk iunui 9.46) iuk 9.4.4.7BCFCM Algorithm The BCFCM algorithm for correcting the bias field and segmenting the image into different clusters can be summarized in the following steps: Step 1.Select initial class prototypes {vi}1SetB1 to equal and very small values (e.g. 0.01). Step 2.Update the partition matrix using Eq.9.40. Advanced Segmentation Techniques 507 Step 3.The prototypes of the clusters are obtained in the form of weighted averages of the patterns using Eq.9.42 Step 4.Estimate the bias term using Eq.9.46. Repeat steps 2-4 till termination. The termination criterion is as follows |Vnew-Vold| 9.47) where ||| is the Euclidean norm,V is a vector of cluster centers, and e is a small number that can be set by the user. 9.4.4.8BCFCM Results In this section, we describe the application of the BCFCM segmentation to syn- thetic images corrupted with multiplicative gain, as well as digital MR phan toms [51 and real brain MR images. The MR phantoms simulated the appear- ance and image characteristics of the T1 weighted images. There are many advantages of using digital phantoms rather than real image data for validating segmentation methods.These advantages include prior knowledge of the true tissue types and control over image parameters such as mean intensity values, noise, and intensity inhomogeneities. We used a high-resolution T1 weighted phantom with in-plane resolution of0.94 mm, Gaussian noise with = 6.0, and 3D linear shading of 7% in each direction.All of the real MR images shown in this section were obtained using a General Electric Signa 1.5 T clinical MR imager with the same in-plane resolution as the phantom. In all the exam- ples,we set the parameter the neighbors effect) to be 0.7,p2,NR=9 a 3 x3 window centered around each pixel),and =0.01.For low SNR im ages, we set a 0.85.The choice of these parameters seems to give the best results. Figure 9.15(a) shows a synthetic test image.This image contains a two-class pattern corrupted by a sinusoidal gain field of higher spatial frequency.The test image is intended to represent two tissue classes,while the sinusoid represents an intensity inhomogeneity.This image was constructed so that it would be dif- ficult to correct using homomorphic filtering or traditional FCM approaches. As shown in Fig.9.15b),FCM algorithm was unable to separate the two classes, while the BCFCM and EM algorithms have succeeded in correcting and classi fying the data as shown in Fig.9.15(c).The estimate of the multiplicative gain 508 Farag,Ahmed, El-Baz and Hassan Figure 9.15:Comparison of segmentation results on a synthetic image cor- rupted by a sinusoidal bias field.(a) The original image,(b) FCM results,c) BCFCM and EM results, and (d) bias field estimations using BCFCM and EM algorithms:this was obtained by scaling the bias field values from 1 to 255. using either BCFCM or EM is presented in Fig.9.15(d).This image was obtained by scaling the values of the bias field from 1 to 255.Although the BCFCM and EM algorithms produced similar results BCFCM was faster to converge to the correct classification, as shown in Fig.9.16. Figures 9.17 and 9.18 present a comparison of segmentation results between FCM,EM, and BCFCM, when applied on T1 weighted MR phantom corrupted with intensity inhomogeneity and noise.From these images, we can see that Advanced Segmentation Techniques 509 Figure9.16: Comparisonofthe performanceofthe proposed BCFCMalgorithm with EM and FCM segmentation when applied to the synthetic two-class image shown in Fig.9.15a traditional FCM was unable to correctly classify the images.Both BCFCM and EM segmented the image into three classes corresponding to background, gray matter (GM), and white matter (WM).BCFCM produced slightly better results than EM due to its ability to cope with noise. Moreover, BCFCM requires far less number of iterations to converge compared to the EM algorithm.Table 9.2 depicts the segmentation accuracy (SA) of the three mentioned method when applied to the MR phantom.SA was measured as follows: Number of correctly classified pixels SA= 100% 9.48) Total number of pixels SA was calculated for different SNR.From the results, we can see that the three methods produced almost similar results for high SNR.BCFCM method, however, was found to be more accurate for lower SNR. 510 FaragAhmed,El-Baz,and Hassan Figure 9.17:Comparison of segmentation results on a MR phantom cor- rupted with 5% Gaussian noise and 20% intensity inhomogeneity:(a) original T1 weighted image, (b) using FCM,c) using EM, and (d using the proposed BCFCM. Advanced Segmentation Techniques 511 Figure 9.18:Comparison of segmentation results on an MR phantom cor- rupted with 5% Gaussian noise and 20% intensity inhomogeneity: (a) original T1 weighted image,b) using FCM, c) using EM, and (d using the proposed BCFCM. 512 Farag,Ahmed, El-Baz and Hassan Table 9.2:Segmentation accuracy of different methods when applied on MR simulated data SNR Segmentation Method 13 db 10 db 8 db FCM 98.92 86.24 78.9 EM 99.12 93.53 85.11 BCFCM 99.25 97.3 93.7 Figure 9.19 shows the results of applying the BCFCM algorithm to segment a real axial-sectioned T1 MR brain.Strong inhomogeneities are apparent in the im- age.The BCFCM algorithm segmented the image into three classes correspond- ing to background, GM, and WM.The bottom right image shows the estimate of the multiplicative gain scaled from 1 to 255. Figure 9.20 shows the results of applying the BCFCM for the segmentation of noisy brain images. The results using traditional FCM without considering the neighborhood field effect and the BCFCM are presented. Notice that the BCFCM segmentation, which uses the the neighborhood field effect, is much less fragmented than the traditional FCM approach.As mentioned before, the relative importance of the regularizing term is inversely proportional to the SNR of MRI signal. It is important to note however, that the incorporation of spatial constraints into the classification has the disadvantage of blurring some fine details.There are current efforts to solve this problem by including contrast information into the classification.High contrast pixels,which usually represent boundaries between objects, should not be included in the neighbors. 9.5Level Sets The mathematical foundation of deformable models represents the confluence of physics and geometry.Geometry serves to represent object shape and physics puts some constrains on how it may vary over space and time.Deformable mod- els have had great success in imaging and computer graphics.Deformable mod- els include snakes and active contours.Snakes are used based on the geometric properties in image data to extract objects and anatomical structures in medi- cal imaging.After initialization,snakes evolve to get the object.The change of Advanced Segmentation Techniques 513 Figure 9.19:Brain MRI example: (upper left) the original MR image corrupted with intensity inhomogeneities. (Upper right) crisp gray matter membership using traditional FCM. (Middle left) crisp gray matter membership using the proposed BCFCM algorithm. (Middle right) the bias-field corrected image using BCFCM.The segmented image and bias field estimate using BCFCM are shown in bottom left and bottom right,respectively. 514 FaragAhmed,El-Baz,and Hassan Figure 9.20:Brain tumor MRI examples.Upper row:Original MR images cor- rupted with salt and pepper noise.Middle row: the segmented images using FCM without any neighborhood consideration.Bottom row:The segmented images using BCFCM 0.85). Advanced Segmentation Techniques 515 snakes with time is guided by differential equations.These equations are de rived from the energy minimization concept to describe the change of snakes with time. The output obtained using snakes depends highly on the initializa- tion.It was found that initial curve has to lie close to the final solution to obtain required results.The initialization is relatively easy in the case of 2D images but in the 3D case it is very difficult.Also the topology change of the solution needs a special regulation to the model. Level sets were invented to handle the problem of changing topology of curves. The level sets has had great success in computer graphics and vision. Also, it was used widely in medical imaging for segmentation and shape re- covery. It proved to have advantages over statistical approaches followed by mathematical morphology. In the following section we will give a brief overview on level sets and its application in image segmentation. 9.5.1 Level Set Function Representation Level sets wasinvented by Osher and Sethian to handle thetopology changes of curves.A simple representation is that a surface intersects with the zero plane to give the curve. When this surfaces changes the curve changes.The surface can be described by the following equation: xt>0ifxxt