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

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.

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 𝐼 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_Ω||∇𝜐||²𝑑𝑥+1𝜖_Ω𝑊(𝜐)𝑑𝑥,(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 [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 Ω_𝑘, 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):^𝐾⁻¹_𝑘=0_Ωlog𝜎^𝑘+𝑐^𝑘−𝐼(𝑥)²2𝜎^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 [1–7, 15, 17], where variants are not involved, for point 𝑥 where (𝑐_𝑘−𝐼(𝑥))² is relatively larger than (𝑐_𝑖−𝐼(𝑥))² 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 (𝑐_𝑘−𝐼(𝑥))² 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 𝐼∈𝐿²(Ω) be an image defined in a bounded, smooth and open domain Ω⊂𝑅². 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𝜎^𝑘+𝑐^𝑘−𝐼²=_𝐾−12𝜎^2𝑘_𝑘=0_Ωlog𝜎^𝑘+𝑐^𝑘−𝐼²2𝜎^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:^𝐾⁻¹_𝑘=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𝜎^𝑘+𝑐^𝑘−𝐼²2𝜎^2𝑘1^𝑧=𝑘.(6)

Note that when 𝜎₀=𝜎₁=⋯=𝜎_𝐾−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 𝑧(𝑥) 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 [8–10, 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_Ω𝜖||∇𝑧||²+1𝜖sin²𝜋𝑧𝑑𝑥(7)𝛾-converges to (4/𝜋)∫_Ω|𝐷𝑧(𝑥)| in 𝐿¹(ℛ) 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^𝐾⁻¹_𝑘=0_Ωlog𝜎^𝑘+𝑐^𝑘−𝐼²2𝜎^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𝐸_𝜖[𝑧,𝑐,𝜎∣𝐼]=_Ω𝜖||∇𝑧||²+1𝜖sin²+𝜆_𝐾−1𝜋𝑧𝑑𝑥_𝑘=0_Ωlog𝜎^𝑘+𝑐^𝑘−𝐼²∶=𝐹_𝜖2𝜎^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 𝐜=(𝑐₀,𝑐₁,…,𝑐_𝐾−1), the variances 𝝈^𝟐=(𝜎²⁰,𝜎²¹,…,𝜎^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,_𝑘(𝑛)²=∑^𝑖∑_𝑗𝐼^𝑖,𝑗−𝑐^𝑘(𝑛)²ℎ𝑧^{(𝑛)𝑖,𝑗}−𝑘∑_𝑖∑𝑘_𝑗=ℎ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𝜖Δ𝑧(𝑥)+𝜋+𝜆_𝐾−1𝜖sin(2𝜋𝑧(𝑥))_𝑘=0ℎ^(𝑧(𝑥)−𝑘)log𝜎_𝑘+||𝑐^𝑘−𝐼||²2𝜎^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 𝑓∈𝐶²(𝑅). If 𝑓^≥−𝛾 for some 𝛾≥0. We define the splitting 𝐹(𝑢)=_Ω𝛾2𝑢²+𝑓(𝑢)𝑑𝑥−_Ω𝛾2𝑢²∶=𝐹₁𝑑𝑥(𝑢)−𝐹²(𝑢).(18) Then both 𝐹(𝑢) and 𝐹²(𝑢) are convex.

The proof is trivial since 𝐹¹^(𝑢)=_Ω𝛾+𝑓^(𝑢)𝑑𝑥≥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 (𝑑²/𝑑𝑧²)sin²𝜋𝑧≥−2𝜋², we split 𝐹_𝜖[𝑧] as𝐹_𝜖[𝑧]=𝐹_𝜖[𝑧]+𝜋²𝜖_Ω𝑧²𝑑𝑥−𝜋²𝜖_Ω∶=𝐹_1𝜖𝑧²𝑑𝑥[𝑧]−𝐹^2𝜖[𝑧].(20)

Similarly, suppose (𝑑²/𝑑𝑧²)ℎ(𝑧−𝑘)≥−𝛾, where 𝛾 depends on function ℎ(⋅). We split 𝐺[𝑧|𝐼] as𝐺[𝑧𝐼]∶=𝐺¹[𝑧𝐼]−𝐺²[𝑧𝐼],(21) where𝐺¹[𝑧∣𝐼]=^𝐾−1_𝑘=0_Ωlog𝜎^𝑘+𝑐^𝑘−𝐼²+2𝜎^2𝑘ℎ(𝑧−𝑘)𝑑𝑥_𝐾−1_𝑘=0_Ω𝛾2log𝜎^𝑘+𝑐^𝑘−𝐼²2𝜎^2𝑘(𝑧−𝑘)²𝐺𝑑𝑥,₂[𝑧∣𝐼]=^𝐾−1_𝑘=0_Ω𝛾2log𝜎^𝑘+𝑐^𝑘−𝐼²2𝜎^2𝑘(𝑧−𝑘)²𝑑𝑥.(22) We simply denote it by𝐸_𝜖[𝑧∣𝐼]=𝐹^1𝜖[𝑧]+𝐺¹[𝑧∣𝐼]−𝐹^2𝜖[𝑧]+𝐺²[𝑧∣𝐼].(23)

By Proposition 1, 𝐹^1𝜖[𝑧], 𝐹^2𝜖[𝑧], 𝐺¹[𝑧∣𝐼] and 𝐺²[𝑧∣𝐼] 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. ThenMin𝜎_𝑘,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)=𝐹²+𝜆𝐺²^𝑧^(𝑛).(25) Under integration by parts (see (15)), the above equation is equivalent to the following PDE.−2𝜖Δ𝑧^(𝑛+1)+𝜋𝜖sin2𝜋𝑧^(𝑛+1)+2𝜋𝜖²𝑧^(𝑛+1)+𝜆^𝐾−1_𝑘=0log𝜎^𝑘+||𝑐^𝑘−𝐼||²2𝜎^2𝑘ℎ^+𝜆𝛾_𝐾−1𝑧^(𝑛+1)−𝑘_𝑘=0log𝜎^𝑘+||𝑐^𝑘−𝐼||²2𝜎^2𝑘𝑧^(𝑛+1)=2𝜋₂−𝑘𝜖𝑧^(𝑛)+𝜆𝛾^𝐾−1_𝑘=0log𝜎^𝑘+||𝑐^𝑘−𝐼||²2𝜎^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+𝛿12𝛿+2𝑥+𝑥2−𝜋𝛿cos𝜋(𝑥2𝛿+1)iiff||−𝑥1++1||𝛿≤<𝛿,𝑥<−𝛿,11𝜋−+c𝑥𝛿os𝑥2𝜋𝛿+1−𝛿iiff𝛿|𝑥<|≤𝑥𝛿<,1−𝛿,0,2−𝑥2−𝜋𝛿cos𝜋(𝑥2𝛿−1)iofth||e𝑥rw−i1se||,≤𝛿,(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.

(a)

(b)

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:𝑥∈⎧⎪⎨⎪⎩Ω⁰Ω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.

(a)

(b)

(c)

(d)

(e)

(f)

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

(a)

(b)

(c)

(d)

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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

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)

(g)

(h)

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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

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),

(a)

(b)

(c)

(d)

(e)

(f)

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

(a)

(b)

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.

Acknowledgment

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