Table of Contents Author Guidelines Submit a Manuscript
Applied Computational Intelligence and Soft Computing
Volume 2011, Article ID 786369, 11 pages
Research Article

A New Framework of Multiphase Segmentation and Its Application to Partial Volume Segmentation

1Department of Mathematics, University of Florida, Gainesville, FL 32611-8105, USA
2Department of Diagnostic Radiology, Yale University, New Haven, CT 06520-8042, USA

Received 14 October 2010; Revised 24 January 2011; Accepted 16 February 2011

Academic Editor: Antonio Di Nola

Copyright © 2011 Fuhua Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We proposed a novel framework of multiphase segmentation based on stochastic theory and phase transition theory. Our main contribution lies in the introduction of a constructed function so that its composition with phase function forms membership functions. In this way, it saves memory space and also avoids the general simplex constraint problem for soft segmentations. The framework is then applied to partial volume segmentation. Although the partial volume segmentation in this paper is focused on brain MR image, the proposed framework can be applied to any segmentation containing partial volume caused by limited resolution and overlapping.

1. Introduction

Roughly speaking, multiphase segmentation is more challenging than two-phase segmentation. There are a lot of explanations for this phenomenon. But for energy-based models, one reason is that it is harder to control the convergence in the implementation for multiphase segmentation than for two-phase segmentation. One extreme case is that the implementation may lead to less phases (local minimum) than expected. Recently, multiphase segmentation has been studied extensively. One way is to employ level set method. For example, papers in [1, 2] use multiple level set functions, and papers in [3, 4] use multiple layers for each level set function. Soft segmentation is another way developed for multiphase segmentation [57]. Let 𝐼 denote an image containing 𝐾 phases. The purpose of soft segmentation is to find a set of functions 𝑝𝑖,(1𝑖𝐾) indicating the probabilities that point 𝑥 belongs to 𝑖th phase.

The similarities between image segmentation and the phase transition theory in material sciences and fluid mechanics have inspired people to borrow some ideas in contemporary material sciences, for example, the diffuse interface model of Cahn-Hilliard [8], and its rigorous mathematical analysis in the framework of Γ convergence approximation by Modica and Mortola [9, 10]. The phase-field relaxation consists in approximating the perimeter using a Cahn-Hilliard type penalization functional [8] with the form𝑃𝜖(𝜖𝜐)=2Ω||||𝜐21𝑑𝑥+𝜖Ω𝑊(𝜐)𝑑𝑥,(1) where 𝑊𝑝+ is a scalar function with exactly two minimizers at 0 and 1 satisfying 𝑊(0)=𝑊(1)=0. The second term of the penalty functional ensures that the values of the material density 𝜐 converges to 0 or 1 as 𝜖0, while the first term controls the perimeter. The parameter 𝜖 can be interpreted as the width of the diffused edge representation in 𝜐. The phase-field approach has been used in topological optimization problems [1113]. In [14], the authors used the phase field to approximate sharp edges and a variational phase field model is derived to compute a shape average of a given number of shapes. In [12], the authors used the phase transition theory in a Cahn-Hilliard inpainting model.

Since applying phase-transition theory into a model can enhance pattern separation and make boundaries smooth, the theory has been widely applied to image segmentation [6, 1517]. The authors in [17] presented a model for image segmentation based upon phase transition theory of Modica and Mortola and discussed its connections to the Mumford-Shah model and some related works. Then Chen et al. extended it under stochastic theory [16]. In [6], Shen proposed a general multiphase stochastic variational fuzzy segmentation model combining stochastic principle and the Modica-Mortola’s phase-transition theory.

There is a common point in the models mentioned above. They all contain a regularity term and a fidelity term and assume that the image is piecewise constant or smooth enough, which makes them vulnerable to noise. Although we can choose large weight for the fidelity term to restrain noise, the segmentation result is sensitive to the choice of the weights between the regularity term and the fidelity term. If the weight of the regularity term is much bigger than that of the fidelity term, then the edge and fine structure can be preserved very well but the noise may not be removed ideally. On the other hand, if we choose the contrary, that is, the weight of fidelity term is much bigger than that of the regularity term, noise will be removed very well but edge may be damaged and some fine structure may be lost. To avoid this dilemma, there have been some studies in recent years aiming at employing nonlocal information of images, such as graph-cut-based method (discrete case) [18] and nonlocal-variation-based method (continuous case) [19]. In these methods, whether a point is an edge depends not only on the local intensity difference, but also on finding how often the similar features of the point have been repeated in the whole domain. By taking nonlocal information, the edge can be well preserved while the noise is smoothed.

Different from using nonlocal information, the framework in this paper uses stochastic theory to restrain noise and improve segmentation result. It can be thought an extension of piecewise constant Mumford-Shah kind models mentioned above. More precisely, we assume that the intensity of each point is a Gaussian distributed random sample. In each phase Ω𝑘, the points follow a same Gaussian distribution with mean 𝑐𝑘 and variance 𝜎𝑘. We assume that the clean true image 𝑢(𝑥) is still piecewise constant (i.e., inside each phase, the intensities are always a constant equal to 𝑐𝑘) but contaminated by a Gaussian noise 𝑛(𝑥), that is, 𝐼(𝑥)=𝑢(𝑥)+𝑛(𝑥). As a result, the intensities of points in a same phase will not be a constant, but a family of samples from a same Gaussian distribution 𝑓𝑘(𝑥), that is, a Gaussian distribution PDF (Probability Density Function). By maximizing the likelihood (joint PDF) of all random samples, the fidelity term becomes the following form (see detail in Section 3):𝐾1𝑘=0Ωlog𝜎𝑘+𝑐𝑘𝐼(𝑥)22𝜎2𝑘𝑝𝑘(𝑥)𝑑𝑥,(2) where 0𝑝𝑘(𝑥)1 is a smooth approximation of characteristic function of 𝑘th phase, and 𝑐𝑘 and 𝜎𝑘 are mean and standard variation, respectively. In application, 𝑝𝑘(𝑥) are usually chosen in such a way that 𝐾1𝑘=0𝑝𝑘(𝑥)=1 holds in the whole domain. So, the model becomes a standard soft segmentation with 𝑝𝑘(𝑥) as probability functions. Compared with those piecewise constant [17, 15, 17], where variants are not involved, for point 𝑥 where (𝑐𝑘𝐼(𝑥))2 is relatively larger than (𝑐𝑖𝐼(𝑥))2 for some 𝑖𝑘 due to noise, the model with fidelity term (2) can still classify it to kth phase (the correct phase) and the noise can therefore be restrained. The reason is that the effect of large (𝑐𝑘𝐼(𝑥))2 will be partly counteracted by the variance. On the other hand, based on probability theory, we know that the probability that 𝑆Ω𝑘 for a connected area 𝑆 with big difference |𝐼𝑘(𝑥)𝑐𝑘| for all 𝑥Ω𝑘 is much smaller than the probability that 𝑥Ω𝑘 for one isolated point 𝑥 with big difference |𝐼𝑘(𝑥)𝑐𝑘|. This fact can guarantee that the model based on (2) can preserve small structure while removing isolated noise.

In this paper, we did not introduce membership functions as an approximation of characteristic function. Instead, we introduce a constructed function (𝑥) so that the composite function (𝑧(𝑥)𝑘) has the property of probability function 𝑝𝑘(𝑥). As a result, the model itself is still a soft segmentation. The advantage of applying constructed function in the model lies in the fact that there will be less variables introduced in the model which makes the discussion and the implementation easier. For example, as long as we know 𝑧(𝑥), the probability 𝑝𝑘(𝑥) is followed by 𝑝𝑘(𝑥)=(𝑧(𝑥)𝑘).

In MRI brain image, there are three tissues, that is, white matter, gray matter, and cerebrospinal fluid (CSF). Due to limited resolution and nonregularity of the boundaries of pure matters, there are some partial volumes formed by overlaps of pure matters. Those partial volumes are hard to be classified to pure matters. We assume each partial volume follows a mixed Gaussian distribution generated by two Gaussian distributions which correspond to pure matters. The proposed frame work can then be applied to partial volume segmentation. The rest of the paper is organized as follows. In Section 2, we develop the framework of multiphase segmentation by combining phase transition theory and Gaussian distribution. Section 3 is the implementation and some considerations. Then in Section 4, we apply the frame work to partial volume segmentation. Experiments are carried out in Section 5. We show by examples the advantage of the proposed model by comparing with other multiphase segmentation models. Finally, we close the paper with a short conclusion.

2. Framework Development

In this section, we develop a framework of multiphase segmentation by integrating phase transition theory and statistics method. Let 𝐼𝐿2(Ω) be an image defined in a bounded, smooth and open domain Ω𝑅2. Suppose the image contains 𝐾 phases, and we take the image as a random field with the following assumptions: (i)at each point 𝑥Ω, the intensity 𝐼(𝑥) is a random variable;(ii)all the random variables {𝐼(𝑥)𝑥Ω} are independent;(iii)in each phase 𝐼Ω𝑘,0𝑘𝐾1, all the random variables {𝐼(𝑥)𝑥Ω𝑘} are identically distributed as a Gaussian distribution with same mean 𝑐𝑘 and same variance 𝜎2𝑘 (which are to be determined).

We want to maximize the likelihood, joint pdf of {𝐼(𝑥),𝑥Ω}, which is equivalent to minimize the following energy (the details can be found in [16]):log𝐿(𝑐,𝜎)=𝐾1𝑘=0Ω𝑘log𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘=𝐾1𝑘=0Ωlog𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘1Ω𝑘.(3)

Second, we borrow the length term from [17] which is based on phase transition theory. Let 𝐼 be the image defined above. The signature function 𝑧(𝑥) is defined by𝑧(𝑥)=𝑘,if𝑥Ω𝑘,𝑘=0,1,,𝐾1.(4) Then the total variation of 𝑧(𝑥) is Ω|𝐷𝑧|, where 𝐷𝑧 denotes the vectorial Radon measure of the total variation (TV) of 𝑧. We have the following relation:𝐾1𝑘=0||𝜕Ω𝑘||Ω||||𝐷𝑧(𝑥)𝐾𝐾1𝑘=0||𝜕Ω𝑘||,(5) where 𝜕Ω𝑘 is the boundary of Ω𝑘. Thus, Ω|𝐷𝑧(𝑥)| works as the length of the edges of all phases. The ideal model is to minimize the sum 𝐸[𝑧,𝑐,𝜎𝐼] of the length term and the negative log-likelihood(3)𝐸[]=𝑧,𝑐,𝜎𝐼Ω||||𝐷𝑧(𝑥)+𝜆𝐾1𝑘=0Ωlog𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘1𝑧=𝑘.(6)

Note that when 𝜎0=𝜎1==𝜎𝐾10, the model is equivalent to piecewise constant Mumford-Shah model. However, this model has intrinsic drawback due to the discreteness of the signature function 𝑧(𝑥) and the characteristic function 1𝑧=𝑘, which will impede the application of PDE-based method. So, we need to use some relaxed version. For the signature function, we introduce its relaxed version via the celebrated model of Modica and Mortola on phase transitions in material science and fluid mechanics. We refer the authors to papers [810, 17] for further understanding to phase transition theory.

Let ̃𝑧(𝑥) be a smoothed version of the signature function 𝑧(𝑥), which is called phase fields. To be simple, we still use the same notation 𝑧(𝑥) to denote the phase field. Modica and Mortola [9] established thatΩ𝜖||||𝑧2+1𝜖sin2𝜋𝑧𝑑𝑥(7)𝛾-converges to (4/𝜋)Ω|𝐷𝑧(𝑥)| in 𝐿1() for phase fields 𝑧(𝑥) that ultimately only take integer values. Now we replace the length term in (6) by (7) when 𝜖 is small enough. For the characteristic function 1𝑧=𝑘, we can use any smooth function 𝑘(𝑥) as an approximation if only 0𝑘(𝑥)1 and 𝐾1𝑘=0 for all 𝑥Ω. The ideal smooth function 𝑘(𝑥) should have the following properties:(i)at each point 𝑥Ω where 𝑧(𝑥)=𝑘, 𝑘(𝑥) is close to 1;(ii)at each point 𝑥Ω where |𝑧(𝑥)𝑘|>0.5, 𝑘(𝑥) is close to 0.

If we can choose (𝑥) satisfying that (𝑥) is close to one at small neighborhood of 0 and close to zero elsewhere, then 𝑘(𝑥) can be denoted as (𝑧(𝑥)𝑘) since 𝑧(𝑥) is almost integer. Then, the fidelity term (3) becomes𝐾1𝑘=0Ωlog𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘(𝑧(𝑥)𝑘).(8)

Now, as the final step, we integrate the relaxed length term (7) and the relaxed fidelity term (8). The new energy functional is𝐸𝜖[]=𝑧,𝑐,𝜎𝐼Ω𝜖||||𝑧2+1𝜖sin2𝜋𝑧𝑑𝑥+𝜆𝐾1𝑘=0Ωlog𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘(𝑧𝑘)=𝐹𝜖[𝑧][].+𝐺𝑧,𝐜,𝝈𝐼(9)

This is the proposed framework of multiphase segmentation. Compared with most existing models, the proposed model is robust to noise and more applicable; moreover, since the phase function 𝑧(𝑥) is involved in the model, segmentation can be derived directly from the phase function.

3. Implementation and Considerations

In our framework (9), there are three groups of parameters to be determined, the means 𝐜=(𝑐0,𝑐1,,𝑐𝐾1), the variances 𝝈𝟐=(𝜎20,𝜎21,,𝜎2𝐾1), and the phase field 𝑧(𝑥). We compute 𝐸𝜖[𝑧,𝐜,𝝈𝐼] regarding 𝑧, 𝐜, and 𝝈 as independent variables. This allows the application of the alternating minimization (AM) scheme, that is, to alternatingly optimize the three conditional energies 𝐸𝜖[𝑧𝐜,𝝈,𝐼],𝐸𝜖[𝐜𝝈,𝑧,𝐼], and 𝐸𝜖[𝝈𝐜,𝑧,𝐼], under the iterations of 𝑧(𝑛)𝐜(𝑛)𝝈(𝑛)𝑧(𝑛+1) given by𝐜(𝑛)=argmin𝐸𝜖𝐜𝝈(𝑛),𝑧(𝑛)𝝈,𝐼,(10)(𝑛)=argmin𝐸𝜖𝝈𝐜(𝑛),𝑧(𝑛)𝑧,𝐼,(11)(𝑛+1)=argmin𝐸𝜖𝑧𝐜(𝑛),𝝈(𝑛).,𝐼(12)

It is well known (see Vogel [20]) that the AM scheme is monotone: 𝐸𝜖𝑧(𝑛+1),𝐜(𝑛+1),𝜎(𝑛+1)𝐼𝐸𝜖𝑧(𝑛),𝐜(𝑛),𝜎(𝑛).𝐼(13) For (10) and (11), one can simply have at pixel level,𝑐𝑘(𝑛)=𝑖𝑖𝐼𝑖,𝑗𝑧(𝑛)𝑖,𝑗𝑘𝑖𝑗𝑧(𝑛)𝑖,𝑗𝜎𝑘,𝑘=0,1,,𝐾1,𝑘(𝑛)2=𝑖𝑗𝐼𝑖,𝑗𝑐𝑘(𝑛)2𝑧(𝑛)𝑖,𝑗𝑘𝑖𝑗𝑧(𝑛)𝑖,𝑗,𝑘𝑘=0,1,,𝐾1,(14) where 𝑧(𝑛)𝑖,𝑗 denotes the computational phase field on the Cartesian image domain at step 𝑛.

For argmin𝐸𝜖[𝑧𝑐(𝑛),𝜎(𝑛),𝐼], since the model is not an active contour model (the integral of the fidelity term cannot be separated by boundaries), it is not proper to use level set method for implementation. Besides, our purpose is not to find the rough segmentation of 𝑧(𝑥) but the exact value of 𝑧(𝑥) so that we can estimate pure matters form the partial volume in MRI brain images using the composition (𝑧(𝑥)𝑘) for 0𝑘𝐾1. Note that the Euler-Lagrange equation of 𝐸𝜖[𝑧𝑐(𝑛),𝜎(𝑛),𝐼] with respect to 𝑧 is𝜋2𝜖Δ𝑧(𝑥)+𝜖sin(2𝜋𝑧(𝑥))+𝜆𝐾1𝑘=0(𝑧(𝑥)𝑘)log𝜎𝑘+||𝑐𝑘||𝐼22𝜎2𝑘=0.(15)

As for most multiphase segmentation models, the energy functional 𝐸𝜖[𝑧,𝑐,𝜎𝐼] is not convex with respect to 𝑧. Thus, simple iteration scheme obtained directly from the Euler-Lagrange equation may not converge. In our application, we adopt the convex-concave procedure (CCCP). Before applying CCCP to our model, we first give a short review on CCCP. We only describe the basic result that is necessary to understand this paper. For more details, we refer the readers to [6, 21]. We also recommend the reader to use selected initial value to help converge (but not guaranteed).

3.1. Convex-Concave Procedure (CCCP)

The convex-concave procedure is a convex splitting method in optimization which was explored by Yuille and Rangarajan [21]. The following theorem is the foundation of CCCP.

Theorem 1. Consider an energy function which is bounded below and is an addition of convex and concave functions: 𝐸𝑥=𝐸convex𝑥+𝐸concave𝑥.(16) Then, the discrete iterative CCCP algorithm given by 𝐸convex𝑥(𝑛+1)=𝐸concave𝑥(𝑛),𝑛=0,1,(17) is guaranteed to monotonically decrease the energy 𝐸(𝑥) as a function of time and to converge to a local minimum or a saddle point of 𝐸(𝑥).

Proposition 1. Let 𝐹(𝑢)=Ω𝑓(𝑢(𝑥))𝑑𝑥, where 𝑓𝐶2(𝑅). If 𝑓𝛾 for some 𝛾0. We define the splitting 𝐹(𝑢)=Ω𝛾2𝑢2+𝑓(𝑢)𝑑𝑥Ω𝛾2𝑢2𝑑𝑥=𝐹1(𝑢)𝐹2(𝑢).(18) Then both 𝐹(𝑢) and 𝐹2(𝑢) are convex.

The proof is trivial since 𝐹1(𝑢)=Ω𝛾+𝑓(𝑢)𝑑𝑥0(19) if 𝑓𝛾 for some 𝛾0.

3.2. Iteration Scheme for Phase Function 𝑧(𝑥)

To find argmin𝐸𝜖[𝑧|𝑐(𝑛),𝜎(𝑛),𝐼], we now apply the CCCP algorithm. In this section, we will temporally suppose the parameters 𝑐𝑘,𝜎𝑘,1𝑘𝐾 are all known for the purpose of statement.

Note that (𝑑2/𝑑𝑧2)sin2𝜋𝑧2𝜋2, we split 𝐹𝜖[𝑧] as𝐹𝜖[𝑧]=𝐹𝜖[𝑧]+𝜋2𝜖Ω𝑧2𝜋𝑑𝑥2𝜖Ω𝑧2𝑑𝑥=𝐹1𝜖[𝑧]𝐹2𝜖[𝑧].(20)

Similarly, suppose (𝑑2/𝑑𝑧2)(𝑧𝑘)𝛾, where 𝛾 depends on function (). We split 𝐺[𝑧|𝐼] as𝐺[]𝑧𝐼=𝐺1[]𝑧𝐼𝐺2[],𝑧𝐼(21) where𝐺1[]=𝑧𝐼𝐾1𝑘=0Ωlog𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘+(𝑧𝑘)𝑑𝑥𝐾1𝑘=0Ω𝛾2log𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘(𝑧𝑘)2𝐺𝑑𝑥,2[]=𝑧𝐼𝐾1𝑘=0Ω𝛾2log𝜎𝑘+𝑐𝑘𝐼22𝜎2𝑘(𝑧𝑘)2𝑑𝑥.(22) We simply denote it by𝐸𝜖[]=𝐹𝑧𝐼1𝜖[𝑧]+𝐺1[]𝐹𝑧𝐼2𝜖[𝑧]+𝐺2[].𝑧𝐼(23)

By Proposition 1, 𝐹1𝜖[𝑧], 𝐹2𝜖[𝑧], 𝐺1[𝑧𝐼] and 𝐺2[𝑧𝐼] are all convex. In order to apply CCCP, we still need the functional to be bounded below. To do that, we must assume that the image is not constant for any phase. Or in detail, we suppose each variance 𝜎2𝑘0. Then𝜎Min𝑘,0𝑘𝐾1>0.(24) Thus, the functional 𝐸𝜖[𝑧𝑐(𝑛),𝜎(𝑛),𝐼] is lower bounded.

Remark 1. This assumption is true for most applications since if the image has been piecewise constant, we do not need to do any segmentation. Even if some phases are exactly constant, the assumption can still be true by adding some noise or spurious dots in the phases.

By Theorem 1, we can now use the CCCP iteration scheme via Frechet derivative, that is,𝐹1+𝜆𝐺1𝑧(𝑛+1)=𝐹2+𝜆𝐺2𝑧(𝑛).(25) Under integration by parts (see (15)), the above equation is equivalent to the following PDE.2𝜖Δ𝑧(𝑛+1)+𝜋𝜖sin2𝜋𝑧(𝑛+1)+2𝜋2𝜖𝑧(𝑛+1)+𝜆𝐾1𝑘=0log𝜎𝑘+||𝑐𝑘||𝐼22𝜎2𝑘𝑧(𝑛+1)𝑘+𝜆𝛾𝐾1𝑘=0log𝜎𝑘+||𝑐𝑘||𝐼22𝜎2𝑘𝑧(𝑛+1)=𝑘2𝜋2𝜖𝑧(𝑛)+𝜆𝛾𝐾1𝑘=0log𝜎𝑘+||𝑐𝑘||𝐼22𝜎2𝑘𝑧(𝑛),𝑘(26) where the terms in the square brackets come from the Euler-Lagrange equation of 𝐸𝜖. Then we can use any existing method to solve 𝑧𝑛+1 such as Gauss-Seidel method. Now we have all the three minimizations (10), (11), and (12). We can use the alternating minimization scheme as discussed at the beginning of this section.

3.3. Construct Approximation Function (𝑥)

We construct 𝛿(𝑥) in such a way that 𝛿(𝑧(𝑥)𝑘) forms an ownership function of 𝑘th phase, where 𝛿(𝑥) is defined as follows:𝛿=(𝑥)1+𝛿2+𝑥2𝛿𝜋cos𝜋(𝑥+1)||||2𝛿if𝑥+1𝛿,1+𝑥if1+𝛿<𝑥<𝛿,2𝛿𝜋cos𝑥𝜋2𝛿+1𝛿if|𝑥|𝛿,1𝑥if𝛿<𝑥<1𝛿,1+𝛿2𝑥2𝛿𝜋𝜋cos(𝑥1)||||2𝛿if𝑥1𝛿,0,otherwise,(27) where 𝛿 is a parameter theoretically satisfying 0<𝛿<0.5 and should be small enough. Figure 1 shows different branches of 𝛿(𝑧(𝑥)𝑘) for 𝑘=1,0,1,2,3. In Figure 2, a three-phase example shows that𝛿(𝑥+1)+𝛿(𝑥)+𝛿(𝑥1)=1(28) holds only for 𝑥[1𝛿,1+𝛿]. In application, we define 𝛿(𝑥)=1 for 𝑥(,𝛿)(𝐾1+𝛿,). Then 𝐾1𝑘=0𝛿(𝑥𝑘)1, and so our model is a standard soft segmentation model.

Figure 1: Different branches of the constructed function 𝛿(𝑥). (a) 𝛿=0.2, (b) 𝛿=0.02.
Figure 2: Function 𝛿(𝑥+1),𝛿(𝑥) and 𝛿(𝑥1) and their sum, where 𝛿=0.2.

Note that although the function 𝛿(𝑥) itself is not a good approximation of 1𝑧=𝑘 based on its graph, the composition 𝛿(𝑧(𝑥)𝑘) is actually a very good approximation of 1𝑧=𝑘 for the reason that, with the length term (7), the phase function 𝑧(𝑥) will ultimately only take integer values which makes the composition 𝛿𝑧(𝑥) mostly evaluated only based on the piece (2𝛿/𝜋)cos(𝑥𝜋/2𝛿)+1𝛿 since |𝑧(𝑥)𝑘|𝛿 for 𝑥 belonging to phase 𝑘. Thus, 𝛿(𝑧(𝑥)𝑘) will ultimately take the value (2𝛿/𝜋)+1𝛿 which approximates 1 very well when 𝛿 is small enough. At the end, we want to mention that as 𝛿(𝑥) is defined in an explicit way, its derivative can be calculated easily. So we can finally use the iteration scheme CCCP developed in the above section.

3.4. Segmentation Decision

Finally, once the iterations are stable, we have two ways to determine the segmentation. One way is hard segmentation, that is, we apply to 𝑧(𝑥) the following hard thresholding decision rule:Ω𝑥0Ωif𝑧(𝑥)<0.5,𝐾1Ωif𝑧(𝑥)𝐾1.5,𝑘1if𝑘1.5𝑧(𝑥)<𝑘0.5,otherwise.(29)

When the choice of (𝑥) makes 𝛿(𝑧(𝑥)𝑘) an ownership function, the model is a soft segmentation. Which scheme (hard segmentation or soft segmentation) should be chosen depends on the purpose of application. If we only want to know the pieces of segments, we should use the hard segmentation. However, in the application to partial volume estimation, we should choose the soft segmentation scheme.

4. Applied to Partial Volume Segmentation

Since our model is a soft segmentation and the membership functions can be viewed as ownership, we can apply our framework to partial volume estimation.

4.1. Introduction to Partial Volume Segmentation

One important application of multiphase segmentation is MRI brain image segmentation. There are three different tissues in human brains, that is, the white matter, the gray matter, and the cerebrospinal fluid (CSF). It is well known that the volume of gray matter has a close relation to some intelligence diseases. Precisely computing the volume of white matter and gray matter can help to diagnose those diseases earlier. On the other hand, due to the limited spatial resolution of imaging equipment, not all voxels in the image contain the same type of tissue, especially the voxels near the tissue borders, which highly likely contain more than one tissue types, called partial volume (PV). Figure 3 shows the principle of partial volume formulation (which is originally used by [22]). The left image contains two phases with higher resolution. Due to lower resolution (half of the left image in each dimension), four squares in the left image contribute to one square in the right image. As a result, the right image contains more phases (for this example, it contains four different phases). When all four subsquares with the same phase in the left image contribute to one square in the right image, the resulted square with lower resolution will be still the same phase as original one, called pure matter. However, when the four subsquares contain different phases, the resulted square will present a phase between the original two phases. In this case, the resulted square is called partial volume. The intensity of the partial volume is a weighted average of original pure matters with combination ratios depending on the number of subsquares of each pure matter in the original image forming the partial volume.

Figure 3: Formulation of partial volume. One partial volume square contains half white matter and half black matter; another partial volume square contains 25 percent white matter and 75 percent of black matter.

This kind of partial volume may cause a big error in the estimation of pure tissue volumes. The error is sometimes as big as 40–60 percent [23]. Thus, partial volume segmentation of MRI images has now received considerable attention. The ideal partial volume segmentation should contain two aspects: finding partial volume and deciding its combination ratio of different pure matters. Most recent work on partial volume segmentation based on statistical principal, for example, the expectation-maximization (EM) method [22, 24]. These methods improved the precision of pure matter estimation. However, they do not contain length term, which makes them sensitive to noise. In [25], the author applied a reparameterized level set algorithm to partial volume segmentation. The method does include the length term. It takes the partial volume part as separated classes that are composition of pure tissues. The drawback of the paper is to use fixed ratios (e.g., 50%) of combinations for partial volumes.

4.2. Apply the framework to Partial Volume Segmentation

We can now apply the proposed framework to partial volume segmentation for brain MR images with mixed Gaussian distribution. In order to apply the proposed framework to partial volume segmentation, we treat the brain image as three different phases: white matter, gray matter, and pure CSF. We calculate the pure matter volumes in a natural way:Volume(phase(𝑘))=Ω𝑝𝑘(𝑥)𝑑𝑥=Ω𝛿(𝑧(𝑥)(𝑘))𝑑𝑥.(30)

5. Experiment and Discussion

Since this paper focus on soft segmentation and applications to partial volume segmentation, our experiments are mostly applied to MRI brain images although a natural image is also presented which is not very close to piecewise constant. Because this work is related to the Sine-Sinc model [17] (where the author proposed a multiphase segmentation model based on phase-transition theory but it is not a soft model), we first give a comparison between our model and the Sine-Sinc model in Figure 4 with synthetic image and MRI brain image. Then in Figure 5, we show three membership functions of a natural image.

Figure 4: (a) Original Artificial synthetic image with Gaussian noise (𝑚=0,𝜎=0.02). (b) Segmentation using the Sine-Sinc model (𝜆=150). (c) Segmentation using the proposed model (𝜆=150,𝛿=0.2). (d) Original brain MR image with Gaussian noise (𝑚=0,𝜎=0.005). (e) Segmentation using the Sine-Sinc model (𝜆=150). (f) Segmentation using the proposed model (𝜆=150,𝛿=0.2).
Figure 5: Soft segmentation. (a) Original image. (b)–(d) Membership functions, whose value is between 0 and 1.

The rest of the examples give a comparison between adaptive fuzzy c-mean method [26] and the proposed model. Figure 6 is the segmentations for a synthetic biased image which is first used by Bresson and Chan [19]. The contour of segmentations in the third column is obtained by thresholding the membership functions (same thing is true for next two figures). Figures 7 and 8 are for MRI brain images.

Figure 6: Robustness to bias: synthetic image. First line: AFCM model; second line: the proposed framework. Middle: soft segmentations. Right: hard segmentations.
Figure 7: Comparison of MRI brain image segmentation. First line: AFCM model. Second line: proposed model. From left to right are original image, gray matter, white matter, and white matter contour by thresholding. The first line contains too much partial volume.
Figure 8: Comparison of MRI brain image segmentation. First line: AFCM model. Second line: proposed model. From left to right are original image, CSF, gray matter, white matter, and white matter contour by thresholding. The first line contains too much partial volume.

Finally, we apply our model to partial volume segmentation using simulated brain images. Then we compare the ground truth of pure matters and our segmentation results. We calculate the errors of partial volume estimation in two ways. One way is based on hard segmentation. Another way is based on soft segmentation.

Figure 9 is a comparison between the ground truth of the original simulated brain MR image and the membership functions obtained using the proposed framework. Figure 9(a) is the original simulated noised image. The corresponding ground truth of white matter, gray matter, and CSF are shown in Figures 9(b), 9(c), and 9(d), respectively. Phase membership functions are shown in Figures 9(e)9(g),

Figure 9: (a) Original emulated brain image. (b) Ground truth of white matter. (c) Ground truth of gray matter. (d) Ground truth of CSF (e) Segmentation of white matter. (f) Segmentation of gray matter. (g) Segmentation of CSF.

We carried out the experiment with 35 consecutive 2D slices of a 3D simulated brain MR image. Then we compared the errors between the Sine-Sinc model and the proposed model. As an average, the errors are shown in Table 1.

Table 1: Error comparison.

Finally, we add series of Gaussian noises to the images with zero mean and different variances. Then compare their errors among the AFCM model, the proposed model with hard segmentation by thresholding, and the proposed model with soft segmentation. The errors are shown in Figure 10(a). From the graph, we can see that as the variance of the noise rises, the error will also rise for all the three cases. However, compared with the AFCM model, the errors using the proposed model rises much more slowly as the variance of the noise rises.

Figure 10: Error comparison.

We also compared the influence on the errors as the parameter 𝛿 in function 𝛿(𝑥) changes. This is shown in Figure 10(b). From the graph, pure matter estimation based on phase function 𝑧(𝑥) looks a little better than the estimation based on the ownership function 𝛿(𝑧(𝑥)𝑘) when 𝛿 is bigger, while the pure matter estimation based on the ownership function 𝛿(𝑧(𝑥)𝑘) is a little better than the estimation based on the phase function 𝑧(𝑥) when 𝛿 is smaller.

6. Conclusions

A critical problem in energy-based multiphase segmentation is the nonconvexity of the energy functional. The problem is more difficult to handle than two-phase case. Level set-based method is proved very successful in multiphase segmentation based on hard segmentation, but they do not work for soft segmentation since the area of different phases may be overlapped and so no clear boundaries between different phases. The soft segmentation is more useful in partial volume segmentation for brain MR images. This paper borrowed the ideas from phase-transition-based methods [6, 15, 17]. The key point of this paper is to construct an approximation function so that the membership functions can be obtained by its composition with phase function. With this constructed function, we can avoid adding new variables for membership functions and so it saves memory space and promotes efficiency. Moreover, since the composition of the constructed function and phase functions forms membership functions, we also avoid the general constraint problem for soft segmentation in implementation. The framework is then applied to partial volume segmentation. The future contains choosing better constructed function (𝑥), and better discretization scheme and iteration scheme.


Publication of this article was funded in part by the University of Florida Open-Access Publishing Fund.


  1. H. K. Zhao, T. Chan, B. Merriman, and S. Osher, “A variational level set approach to multiphase motion,” Journal of Computational Physics, vol. 127, no. 1, pp. 179–195, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  2. L. A. Vese and T. F. Chan, “A multiphase level set framework for image segmentation using the Mumford and Shah model,” International Journal of Computer Vision, vol. 50, no. 3, pp. 271–293, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  3. J. Lie, M. Lysaker, and X. C. Tai, “A variant of the level set method and applications to image segmentation,” Mathematics of Computation, vol. 75, no. 255, pp. 1155–1174, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  4. G. Chung and L. A. Vese, “Energy minimization based segmentation and denoising using a multilayer level set approach,” in Energy Minimization Methodsin Computer Vision and Pattern Recognition, vol. 3757 of Lecture Notes in Computer Science, pp. 439–455, 2005.
  5. S. R. Thiruvenkadam, S. Arcot, and Y. Chen, “A PDE based method for fuzzy classification of medical images,” in Proceedings of the International Conference on Image Processing (ICIP '06), vol. 481, pp. 1805–1808, 2006. View at Publisher · View at Google Scholar
  6. J. Shen, “A stochastic-variational model for soft Mumford-Shah segmentation,” International Journal of Biomedical Imaging, vol. 2006, Article ID 92329, 14 pages, 2006. View at Publisher · View at Google Scholar
  7. S. Esedoglu and Y. H. R. Tsai, “Threshold dynamics for the piecewise constant Mumford-Shah functional,” Journal of Computational Physics, vol. 211, no. 1, pp. 367–384, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  8. J. W. Cahn and J. E. Hilliard, “Free energy of a nonuniform system. I. Interfacial free energy,” The Journal of Chemical Physics, vol. 28, no. 2, pp. 258–267, 1958. View at Google Scholar
  9. L. Modica and S. Mortola, “Un esempio di Γ-convergenza,” Bollettino della Unione Matematica Italiana, vol. 14, pp. 285–299, 1977. View at Google Scholar
  10. L. Modica, “The gradient theory of phase transitions and the minimal interface criterion,” Archive for Rational Mechanics and Analysis, vol. 98, no. 2, pp. 123–142, 1987. View at Publisher · View at Google Scholar · View at MathSciNet
  11. B. Bourdin and A. Chambolle, “Design-dependent loads in topology optimization,” ESAIM Control, Optimisation and Calculus of Variations, no. 9, pp. 19–48, 2003. View at Publisher · View at Google Scholar
  12. M. Burger and R. Stainko, “Phase-field relaxation of topology optimization with local stress constraints,” SIAM Journal on Control and Optimization, vol. 45, no. 4, pp. 1447–1466, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  13. M. Y. Wang and S. Zhou, “A variational method for structural topology optimization,” Computer Modeling in Engineering and Sciences, vol. 6, no. 6, pp. 547–566, 2004. View at Google Scholar
  14. M. Rumpf and B. Wirth, “A nonlinear elastic shape averaging approach,” SIAM Journal on Imaging Sciences, vol. 2, pp. 800–833, 2009. View at Google Scholar
  15. J. H. An and Y. Chen, “Region based image segmentation using modified Mumford-Shah algorithm,” in Proceedings of the 1st International Conference on Scale Space and Variation Variational Methods in Computer Vision, vol. 4485, pp. 733–742, 2007.
  16. F. Chen, Y. Chen, and H. D. Tagare, “An extension of sine-sinc model based on logarithm of likelihood,” in Proceedings of the International Conference on Image Processing, Computer Vision, and Pattern Recognition (IPCV '08), vol. 1, pp. 222–227, 2008.
  17. Y. M. Jung, S. H. Kang, and J. Shen, “Multiphase image segmentation via Modica-Mortola phase transition,” SIAM Journal on Applied Mathematics, vol. 67, no. 5, pp. 1213–1232, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  18. Y. Boykov and G. Funka-Lea, “Graph cuts and efficient N-D image segmentation,” International Journal of Computer Vision, vol. 70, no. 2, pp. 109–131, 2006. View at Publisher · View at Google Scholar
  19. X. Bresson and T. F. Chan, “Non-local unsupervised variational image segmentation model,” UCLA CAM Report 08-67, 2008. View at Google Scholar
  20. C. Vogel, Computational Methods for Inverse Problems, SIAM, Philadelphia, Pa, USA, 2002.
  21. A. L. Yuille and A. Rangarajan, “The concave-convex procedure,” Neural Computation, vol. 15, no. 4, pp. 915–936, 2003. View at Publisher · View at Google Scholar · View at PubMed
  22. K. V. Leeput, F. Maes, D. Vandermeulen, and P. Suetens, “A unifying frame-work for partial volume segmentation of brain MR images,” IEEE Transactions on Medical Imaging, vol. 22, pp. 105–119, 2003. View at Google Scholar
  23. W. J. Niessen, K. L. Vincken, J. Weickert, B. M. Ter Haar Romeny, and M. A. Viergever, “Multiscale segmentation of three-dimensional MR brain images,” International Journal of Computer Vision, vol. 31, no. 2, pp. 185–202, 1999. View at Google Scholar · View at Scopus
  24. D. Eremina, X. Li, W. Zhu, J. Wang, and Z. Liang, “Investigation on an EM framework for partial volume image segmentation,” in Medical Imaging: Image Processing, vol. 6144 of Proceedings of SPIE, pp. 1–9, February 2006. View at Publisher · View at Google Scholar
  25. H. D. Targare, Y. Chen, and R. K. Fulbright, “A reparameterized level set algorithm and its comparison with EM-based partial volume segmentation of MR brain images,” SPIE 2008.
  26. D. L. Pham and J. L. Prince, “An adaptive fuzzy C-means algorithm for image segmentation in the presence of intensity inhomogeneities,” Pattern Recognition Letters, vol. 20, no. 1, pp. 57–68, 1999. View at Publisher · View at Google Scholar