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

Authors: Fuhua Chen; Yunmei Chen; Hemant D. Tagare
Journal: Applied Computational Intelligence and Soft Computing
 (2011)

---

## Abstract

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.


---

## Body

## 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 [5–7]. Let I denote an image containing K phases. The purpose of soft segmentation is to find a set of functions pi,(1≤i≤K) indicating the probabilities that point x belongs to ith 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(1)Pϵ(υ)=ϵ2∫Ω|∇υ|2dx+1ϵ∫ΩW(υ)dx,
where W:ℝp→ℝ∪+∞ is a scalar function with exactly two minimizers at 0 and 1 satisfying W(0)=W(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 [11–13]. 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, 15–17]. 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Ωk, the points follow a same Gaussian distribution with mean ck and variance σk. We assume that the clean true image u(x) is still piecewise constant (i.e., inside each phase, the intensities are always a constant equal to ck) but contaminated by a Gaussian noise n(x), that is, I(x)=u(x)+n(x). 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  fk(x), 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):(2)∑k=0K-1∫Ω(log⁡σk+(ck-I(x))22σk2)pk(x)dx,
where 0≤pk(x)≤1 is a smooth approximation of characteristic function of kth phase, and ck and σk are mean and standard variation, respectively. In application, pk(x) are usually chosen in such a way that ∑k=0K-1pk(x)=1 holds in the whole domain. So, the model becomes a standard soft segmentation with pk(x) as probability functions. Compared with those piecewise constant [1–7, 15, 17], where variants are not involved, for point x where (ck-I(x))2 is relatively larger than (ci-I(x))2 for some  i≠k 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 (ck-I(x))2 will be partly counteracted by the variance. On the other hand, based on probability theory, we know that the probability that S⊂Ωk for a connected area S with big difference |Ik(x)-ck| for all x∈Ωk is much smaller than the probability that x∈Ωk for one isolated point x with big difference |Ik(x)-ck|. 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 functionh(x) so that the composite function h∘(z(x)-k) has the property of probability function pk(x). 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 z(x), the probability pk(x) is followed by pk(x)=h∘(z(x)-k).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 Section2, 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. LetI∈L2(Ω) be an image defined in a bounded, smooth and open domain Ω⊂R2. Suppose the image contains K phases, and we take the image as a random field with the following assumptions:  (i)
at each pointx∈Ω, the intensity I(x) is a random variable;(ii)
all the random variables{I(x)∣x∈Ω} are independent;(iii)
in each phaseI∣Ωk,0≤k≤K-1, all the random variables {I(x):x∈Ωk} are identically distributed as a Gaussian distribution with same mean ck and same variance σk2 (which are to be determined).We want to maximize the likelihood, joint pdf of{I(x),x∈Ω}, which is equivalent to minimize the following energy (the details can be found in [16]):(3)-log⁡L(c,σ)=∑k=0K-1∫Ωk(log⁡σk+(ck-I)22σk2)=∑k=0K-1∫Ω(log⁡σk+(ck-I)22σk2)1Ωk.Second, we borrow the length term from [17] which is based on phase transition theory. Let I be the image defined above. The signature functionz(x) is defined by(4)z(x)=k,ifx∈Ωk,k=0,1,…,K-1.
Then the total variation of z(x) is ∫Ω|Dz|, where Dz denotes the vectorial Radon measure of the total variation (TV) of z. We have the following relation:(5)∑k=0K-1|∂Ωk|≤∫Ω|Dz(x)|≤K∑k=0K-1|∂Ωk|,
where ∂Ωk is the boundary of Ωk. Thus, ∫Ω|Dz(x)| works as the length of the edges of all phases. The ideal model is to minimize the sum E[z,c,σ∣I] of the length term and the negative log-likelihood(3)(6)E[z,c,σ∣I]=∫Ω|Dz(x)|+λ∑k=0K-1∫Ω(log⁡σk+(ck-I)22σk2)1z=k.Note that whenσ0=σ1=⋯=σK-1≠0, the model is equivalent to piecewise constant Mumford-Shah model. However, this model has intrinsic drawback due to the discreteness of the signature function z(x) and the characteristic function 1z=k, 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 [8–10, 17] for further understanding to phase transition theory.Letz̃(x) be a smoothed version of the signature function z(x), which is called phase fields. To be simple, we still use the same notation z(x) to denote the phase field. Modica and Mortola [9] established that(7)∫Ω[ϵ|∇z|2+1ϵsin2⁡πz]dxγ-converges to (4/π)∫Ω|Dz(x)| in L1(ℛ) for phase fields z(x) that ultimately only take integer values. Now we replace the length term in (6) by (7) when ϵ is small enough. For the characteristic function 1z=k, we can use any smooth function hk(x) as an approximation if only 0≤hk(x)≤1 and ∑k=0K-1 for all x∈Ω. The ideal smooth function hk(x) should have the following properties:(i)
at each pointx∈Ω where z(x)=k, hk(x) is close to 1;(ii)
at each pointx∈Ω where |z(x)-k|>0.5, hk(x) is close to 0.If we can chooseh(x) satisfying that h(x) is close to one at small neighborhood of 0 and close to zero elsewhere, then hk(x) can be denoted as h(z(x)-k) since z(x) is almost integer. Then, the fidelity term (3) becomes(8)∑k=0K-1∫Ω(log⁡σk+(ck-I)22σk2)h(z(x)-k).Now, as the final step, we integrate the relaxed length term (7) and the relaxed fidelity term (8). The new energy functional is(9)Eϵ[z,c,σ∣I]=∫Ω[ϵ|∇z|2+1ϵsin2⁡πz]dx+λ∑k=0K-1∫Ω(log⁡σk+(ck-I)22σk2)h(z-k):=Fϵ[z]+G[z,c,σ∣I].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 functionz(x) 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 c=(c0,c1,…,cK-1), the variances σ2=(σ02,σ12,…,σK-12), and the phase field z(x). We compute Eϵ[z,c,σ∣I] regarding z, c, and σ as independent variables. This allows the application of the alternating minimization (AM) scheme, that is, to alternatingly optimize the three conditional energies Eϵ[z∣c,σ,I],Eϵ[c∣σ,z,I], and Eϵ[σ∣c,z,I], under the iterations of z(n)→c(n)→σ(n)→z(n+1) given by(10)c(n)=argmin⁡Eϵ[c∣σ(n),z(n),I],(11)σ(n)=argmin⁡Eϵ[σc(n),z(n),I],(12)z(n+1)=argmin⁡Eϵ[z∣c(n),σ(n),I].It is well known (see Vogel [20]) that the AM scheme is monotone: (13)Eϵ[z(n+1),c(n+1),σ(n+1)∣I]≤Eϵ[z(n),c(n),σ(n)∣I].
For (10) and (11), one can simply have at pixel level,(14)ck(n)=∑i∑iIi,jh(zi,j(n)-k)∑i∑jh(zi,j(n)-k),k=0,1,…,K-1,(σk(n))2=∑i∑j(Ii,j-ck(n))2h(zi,j(n)-k)∑i∑jh(zi,j(n)-k),k=0,1,…,K-1,
where zi,j(n) denotes the computational phase field on the Cartesian image domain at step n.ForargminEϵ[z∣c(n),σ(n),I], 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 z(x) but the exact value of z(x) so that we can estimate pure matters form the partial volume in MRI brain images using the composition h∘(z(x)-k) for 0≤k≤K-1. Note that the Euler-Lagrange equation of Eϵ[z∣c(n),σ(n),I] with respect to z is(15)-2ϵΔz(x)+πϵsin⁡(2πz(x))+λ∑k=0K-1h′(z(x)-k)(log⁡σk+|ck-I|22σk2)=0.As for most multiphase segmentation models, the energy functionalEϵ[z,c,σ∣I] is not convex with respect to z. 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:(16)E(x⃗)=Econvex(x⃗)+Econcave(x⃗).
Then, the discrete iterative CCCP algorithm given by
(17)∇Econvex(x⃗(n+1))=-∇Econcave(x⃗(n)),n=0,1,…
is guaranteed to  monotonically decrease the energy E(x⃗) as a function of time and to converge to a local minimum or a saddle point of E(x⃗).Proposition 1.
LetF(u)=∫Ωf(u(x))dx, where f∈C2(R). If f′′≥-γ for some γ≥0. We define the splitting
(18)F(u)=∫Ω(γ2u2+f(u))dx-∫Ωγ2u2dx:=F1(u)-F2(u).
Then both F(u) and F2(u) are convex.The proof is trivial since(19)(F1)′′(u)=∫Ω(γ+f′′(u))dx≥0
if f′′≥-γ for some γ≥0.


### 3.2. Iteration Scheme for Phase Functionz(x)

To findargminEϵ[z|c(n),σ(n),I], we now apply the CCCP algorithm. In this section, we will temporally suppose the parameters ck,σk,1≤k≤K are all known for the purpose of statement.Note that(d2/dz2)sin2πz≥-2π2, we split Fϵ[z] as(20)Fϵ[z]={Fϵ[z]+π2ϵ∫Ωz2dx}-π2ϵ∫Ωz2dx:=Fϵ1[z]-Fϵ2[z].Similarly, suppose(d2/dz2)h(z-k)≥-γ, where γ depends on function h(·). We split G[z|I] as(21)G[zI]:=G1[zI]-G2[zI],
where(22)G1[z∣I]=∑k=0K-1∫Ω[log⁡σk+(ck-I)22σk2]h(z-k)dx+∑k=0K-1∫Ωγ2[log⁡σk+(ck-I)22σk2](z-k)2dx,G2[z∣I]=∑k=0K-1∫Ωγ2[log⁡σk+(ck-I)22σk2](z-k)2dx.
We simply denote it by(23)Eϵ[z∣I]={Fϵ1[z]+G1[z∣I]}-{Fϵ2[z]+G2[z∣I]}.By Proposition1, Fϵ1[z], Fϵ2[z], G1[z∣I] and G2[z∣I] 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 σk2≠0. Then(24)Min⁡{σk,0≤k≤K-1}>0.
Thus, the functional Eϵ[z∣c(n),σ(n),I] 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 Theorem1, we can now use the CCCP iteration scheme via Frechet derivative, that is,(25)(F1+λG1)′(z(n+1))=(F2+λG2)′(z(n)).
Under integration by parts (see (15)), the above equation is equivalent to the following PDE.(26)[-2ϵΔz(n+1)+πϵsin⁡2πz(n+1)]+2π2ϵz(n+1)+λ[∑k=0K-1(log⁡σk+|ck-I|22σk2)h′(z(n+1)-k)]+λγ∑k=0K-1(log⁡σk+|ck-I|22σk2)(z(n+1)-k)=2π2ϵz(n)+λγ∑k=0K-1(log⁡σk+|ck-I|22σk2)(z(n)-k),
where the terms in the square brackets come from the Euler-Lagrange equation of Eϵ. Then we can use any existing method to solve zn+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 Functionh(x)

We constructhδ(x) in such a way that hδ(z(x)-k) forms an ownership function of kth phase, where hδ(x) is defined as follows:(27)hδ(x)={1+δ2+x2-δπcos⁡⁡π(x+1)2δif|x+1|≤δ,1+xif-1+δ<x<-δ,2δπcos⁡⁡xπ2δ+1-δif|x|≤δ,1-xifδ<x<1-δ,1+δ2-x2-δπcos⁡⁡π(x-1)2δif|x-1|≤δ,0,otherwise,
where δ is a parameter theoretically satisfying 0<δ<0.5 and should be small enough. Figure 1 shows different branches of hδ∘(z(x)-k) for k=-1,0,1,2,3. In Figure 2, a three-phase example shows that(28)hδ(x+1)+hδ(x)+hδ(x-1)=1
holds only for x∈[-1-δ,1+δ]. In application, we define hδ(x):=1 for x∈(-∞,-δ)∪(K-1+δ,∞). Then ∑k=0K-1hδ(x-k)≡1, and so our model is a standard soft segmentation model.Different branches of the constructed functionhδ(x). (a) δ=0.2, (b) δ=0.02.
(a)(b)Figure 2
Functionhδ(x+1),hδ(x) and hδ(x-1) and their sum, where δ=0.2.Note that although the functionhδ(x) itself is not a good approximation of 1z=k based on its graph, the composition hδ∘(z(x)-k) is actually a very good approximation of 1z=k for the reason that, with the length term (7), the phase function z(x) will ultimately only take integer values which makes the composition hδ∘z(x) mostly evaluated only based on the piece (2δ/π)cos(xπ/2δ)+1-δ since |z(x)-k|≤δ for x belonging to phase k. Thus, hδ∘(z(x)-k) 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 hδ(x) 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 toz(x) the following hard thresholding decision rule:(29)x∈{Ω0ifz(x)<0.5,ΩK-1ifz(x)≥K-1.5,Ωk-1ifk-1.5≤z(x)<k-0.5,otherwise.When the choice ofh(x) makes hδ∘(z(x)-k) 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.


## 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:(16)E(x⃗)=Econvex(x⃗)+Econcave(x⃗).
Then, the discrete iterative CCCP algorithm given by
(17)∇Econvex(x⃗(n+1))=-∇Econcave(x⃗(n)),n=0,1,…
is guaranteed to  monotonically decrease the energy E(x⃗) as a function of time and to converge to a local minimum or a saddle point of E(x⃗).Proposition 1.
LetF(u)=∫Ωf(u(x))dx, where f∈C2(R). If f′′≥-γ for some γ≥0. We define the splitting
(18)F(u)=∫Ω(γ2u2+f(u))dx-∫Ωγ2u2dx:=F1(u)-F2(u).
Then both F(u) and F2(u) are convex.The proof is trivial since(19)(F1)′′(u)=∫Ω(γ+f′′(u))dx≥0
if f′′≥-γ for some γ≥0.


## 3.2. Iteration Scheme for Phase Functionz(x)

To findargminEϵ[z|c(n),σ(n),I], we now apply the CCCP algorithm. In this section, we will temporally suppose the parameters ck,σk,1≤k≤K are all known for the purpose of statement.Note that(d2/dz2)sin2πz≥-2π2, we split Fϵ[z] as(20)Fϵ[z]={Fϵ[z]+π2ϵ∫Ωz2dx}-π2ϵ∫Ωz2dx:=Fϵ1[z]-Fϵ2[z].Similarly, suppose(d2/dz2)h(z-k)≥-γ, where γ depends on function h(·). We split G[z|I] as(21)G[zI]:=G1[zI]-G2[zI],
where(22)G1[z∣I]=∑k=0K-1∫Ω[log⁡σk+(ck-I)22σk2]h(z-k)dx+∑k=0K-1∫Ωγ2[log⁡σk+(ck-I)22σk2](z-k)2dx,G2[z∣I]=∑k=0K-1∫Ωγ2[log⁡σk+(ck-I)22σk2](z-k)2dx.
We simply denote it by(23)Eϵ[z∣I]={Fϵ1[z]+G1[z∣I]}-{Fϵ2[z]+G2[z∣I]}.By Proposition1, Fϵ1[z], Fϵ2[z], G1[z∣I] and G2[z∣I] 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 σk2≠0. Then(24)Min⁡{σk,0≤k≤K-1}>0.
Thus, the functional Eϵ[z∣c(n),σ(n),I] 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 Theorem1, we can now use the CCCP iteration scheme via Frechet derivative, that is,(25)(F1+λG1)′(z(n+1))=(F2+λG2)′(z(n)).
Under integration by parts (see (15)), the above equation is equivalent to the following PDE.(26)[-2ϵΔz(n+1)+πϵsin⁡2πz(n+1)]+2π2ϵz(n+1)+λ[∑k=0K-1(log⁡σk+|ck-I|22σk2)h′(z(n+1)-k)]+λγ∑k=0K-1(log⁡σk+|ck-I|22σk2)(z(n+1)-k)=2π2ϵz(n)+λγ∑k=0K-1(log⁡σk+|ck-I|22σk2)(z(n)-k),
where the terms in the square brackets come from the Euler-Lagrange equation of Eϵ. Then we can use any existing method to solve zn+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 Functionh(x)

We constructhδ(x) in such a way that hδ(z(x)-k) forms an ownership function of kth phase, where hδ(x) is defined as follows:(27)hδ(x)={1+δ2+x2-δπcos⁡⁡π(x+1)2δif|x+1|≤δ,1+xif-1+δ<x<-δ,2δπcos⁡⁡xπ2δ+1-δif|x|≤δ,1-xifδ<x<1-δ,1+δ2-x2-δπcos⁡⁡π(x-1)2δif|x-1|≤δ,0,otherwise,
where δ is a parameter theoretically satisfying 0<δ<0.5 and should be small enough. Figure 1 shows different branches of hδ∘(z(x)-k) for k=-1,0,1,2,3. In Figure 2, a three-phase example shows that(28)hδ(x+1)+hδ(x)+hδ(x-1)=1
holds only for x∈[-1-δ,1+δ]. In application, we define hδ(x):=1 for x∈(-∞,-δ)∪(K-1+δ,∞). Then ∑k=0K-1hδ(x-k)≡1, and so our model is a standard soft segmentation model.Different branches of the constructed functionhδ(x). (a) δ=0.2, (b) δ=0.02.
(a)(b)Figure 2
Functionhδ(x+1),hδ(x) and hδ(x-1) and their sum, where δ=0.2.Note that although the functionhδ(x) itself is not a good approximation of 1z=k based on its graph, the composition hδ∘(z(x)-k) is actually a very good approximation of 1z=k for the reason that, with the length term (7), the phase function z(x) will ultimately only take integer values which makes the composition hδ∘z(x) mostly evaluated only based on the piece (2δ/π)cos(xπ/2δ)+1-δ since |z(x)-k|≤δ for x belonging to phase k. Thus, hδ∘(z(x)-k) 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 hδ(x) 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 toz(x) the following hard thresholding decision rule:(29)x∈{Ω0ifz(x)<0.5,ΩK-1ifz(x)≥K-1.5,Ωk-1ifk-1.5≤z(x)<k-0.5,otherwise.When the choice ofh(x) makes hδ∘(z(x)-k) 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). Figure3 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:(30)Volume(phase(k))=∫Ωpk(x)dx=∫Ωhδ(z(x)-(k))dx.


## 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). Figure3 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:(30)Volume(phase(k))=∫Ωpk(x)dx=∫Ωhδ(z(x)-(k))dx.


## 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.(a) Original Artificial synthetic image with Gaussian noise (m=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 (m=0,σ=0.005). (e) Segmentation using the Sine-Sinc model (λ=150). (f) Segmentation using the proposed model (λ=150,δ=0.2).
(a)(b)(c)(d)(e)(f)Soft segmentation. (a) Original image. (b)–(d) Membership functions, whose value is between 0 and 1.
(a)(b)(c)(d)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.Robustness to bias: synthetic image. First line: AFCM model; second line: the proposed framework. Middle: soft segmentations. Right: hard segmentations.
(a)(b)(c)(d)(e)(f)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.
(a)(b)(c)(d)(e)(f)(g)(h)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.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)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.Figure9 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),(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.
(a)(b)(c)(d)(e)(f)(g)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 Table1.Table 1
Error comparison.
MethodsWhite matterGray matterOld model6.83%7.22%New model4.68%2.73%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 Figure10(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.Error comparison.
(a)(b)We also compared the influence on the errors as the parameterδ in function hδ(x) changes. This is shown in Figure 10(b). From the graph, pure matter estimation based on phase function z(x) looks a little better than the estimation based on the ownership function hδ(z(x)-k) when δ is bigger, while the pure matter estimation based on the ownership function hδ(z(x)-k) is a little better than the estimation based on the phase function z(x) 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 h(x), and better discretization scheme and iteration scheme.



---

*Source: 786369-2011-04-26.xml*