# A New and Fast Multiphase Image Segmentation Model for Color Images

Authors: Yunyun Yang; Boying Wu
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

This paper presents a new and fast multiphase image segmentation model for color images. We propose our model by incorporating the globally convex image segmentation method and the split Bregman method into the piecewise constant multiphase Vese-Chan model for color images. We have applied our model to many synthetic and real color images. Numerical results show that our model can segment color images with multiple regions and represent boundaries with complex topologies, including triple junctions. Comparison with the Vese-Chan model demonstrates the efficiency of our model. Besides, our model does not require a priori denoising step and is robust with respect to noise.


---

## Body

## 1. Introduction

Image segmentation [1–8] is an important technique for detecting objects and analyzing images in computer vision and image processing. Image segmentation is concerned with the task of partitioning a given image into different classes or regions corresponding to different objects and the background in the image.Mumford and Shah [9] formulated the image segmentation problem and proposed the famous Mumford-Shah model. For a special case of the Mumford-Shah model, Chan and Vese [1] proposed the famous Chan-Vese model without using the image gradient. Chan and Vese gave the two-phase level set formulation for gray images in [1]. Then, the authors extended their scalar Chan-Vese algorithm to the vector-valued case in [2]. There are many other image segmentation models for color images, and we mention the works in [10, 11].Either the scalar Chan-Vese algorithm in [1] or the Chan-Vese model for vector-valued images in [2] is mainly for images with two phases. To deal with images with multiple regions, Vese and Chan extended their original two-phase Chan-Vese model [1, 2] to a multiphase model by using a multiphase level set formulation in [3]. In their multiphase Vese-Chan model [3], multiple regions can be represented by multiple level set functions. Vese and Chan described their model in two cases: piecewise constant case and piecewise smooth case. In this paper, we mainly focus on the piecewise constant multiphase Vese-Chan model.In fact, Zhao et al. [12], Samson et al. [13], and Paragios and Deriche [14] have already proposed several multiphase image segmentation models to segment images with multiple regions before the multiphase Vese-Chan model. However, these models have the natural problems of vacuum and overlap. Compared with these models, the multiphase Vese-Chan model has several advantages. Firstly, it automatically avoids the problems of vacuum and overlap by construction. Secondly, it needs fewer level set functions to represent the same number of phases in the piecewise constant case. Finally it can represent boundaries with complex topologies, including triple junctions.Recently, the split Bregman method has been applied to solve image segmentation problem more efficiently. The efficiency of the split Bregman method for image segmentation has been demonstrated in [15–18].In this paper, we propose a fast multiphase image segmentation model in a variational level set formulation for color images. Our model is an extension of our previous model for gray images proposed in [19]. Our model has been applied to many synthetic and real color images. Numerical results show that our model has the advantages of the original Vese-Chan model [3] for multiphase image segmentation, but our model is much more efficient. The robustness of our model to noise has also been demonstrated by the numerical results.The outline of the paper is as follows. We review our improved multiphase image segmentation model for gray images in Section2.1. In Section 2.2, we extend our model to work for color images. Then, we apply the split Bregman method to solve our proposed minimization problem more efficiently in Section 2.3. Section 3 gives the algorithm and numerical results of our model. We conclude this paper in Section 4.


## 2. Our Model

### 2.1. Our Model for Gray Images

In our previous work [19], we have proposed an improved active contour model for multiphase image segmentation based on the piecewise constant multiphase Vese-Chan model [3], the globally convex image segmentation method [20] and the split Bregman method [15–17]. The model we proposed in [19] is mainly for gray images. For simplicity, we only review the four-phase model in this section.LetΩ⊂ℜ2 be the image domain and u0:Ω→ℜ be a given gray level image. Our proposed model is given in the following level set formulation:
(2.1)min⁡0≤ϕ1,ϕ2≤1E(ϕ1,ϕ2)=min⁡0≤ϕ1,ϕ2≤1(|∇ϕ1|g+|∇ϕ2|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉),
where ϕ1 and ϕ2 are two level set functions, and μ is a positive constant. |∇ϕi|g and 〈ϕi,ri〉 are defined as
(2.2)|∇ϕi|g=∫g(|∇u0(x)|)|∇ϕi(x)|dx,〈ϕi,ri〉=∫ϕi(x)ri(x)dx,
where
(2.3)g(ξ)=11+β|ξ|2
is an edge detector function [17, 18, 21], and β is a parameter that determines the detail level of the segmentation. r1 and r2 are defined as
(2.4)r1=((u0-c11)2-(u0-c01)2)ϕ2+((u0-c10)2-(u0-c00)2)(1-ϕ2),r2=((u0-c11)2-(u0-c10)2)ϕ1+((u0-c01)2-(u0-c00)2)(1-ϕ1),
where c=(c11,c10,c01,c00) is a constant vector of averages in the four phases.Givenϕ1(0,x)=ϕ1,0(x),ϕ2(0,x)=ϕ2,0(x), once the minimization problem (2.1) is solved, the image domain Ω can be partitioned into four regions Ωi(i=1,2,3,4) by thresholding the level set functions as follows:(2.5)Ω1={x:ϕ1(x)>α,ϕ2(x)>α},Ω2={x:ϕ1(x)>α,ϕ2(x)<α},Ω3={x:ϕ1(x)<α,ϕ2(x)>α},Ω4={x:ϕ1(x)<α,ϕ2(x)<α},
for some α∈(0,1), we choose α=0.5 in this paper. This is different from [3] where the zero level set is used to identify the boundary. We use the α=0.5 level set to identify the boundary because in our paper we have restricted 0≤ϕi≤1. The average vector c=(c11,c10,c01,c00) can be obtained by(2.6)c11=mean(u0)inΩ1,c10=mean(u0)inΩ2,c01=mean(u0)inΩ3,c00=mean(u0)inΩ4.


### 2.2. Our Model for Color Images

In this section, we extend our model in Section2.1 from gray images to color images. Let u0:Ω→ℜ3 be a given color image and u0i(i=1,2,3) the three channels of u0. Let ci=(c1i,c2i,c3i,c4i)(i=1,2,3) be the average vector of the ith channel.The extension of our model to color images is still in the same formulation as (2.1):(2.7)min⁡0≤ϕ1,ϕ2≤1E(ϕ1,ϕ2)=min⁡0≤ϕ1,ϕ2≤1(|∇ϕ1|g+|∇ϕ2|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉).The difference between (2.7) and (2.1) lies in the definition of r1 and r2. In our model for color images (2.7), r1 and r2 are defined differently from (2.4):
(2.8)r1=(∑i=13(u0i-c1i)2-∑i=13(u0i-c3i)2)ϕ2+(∑i=13(u0i-c2i)2-∑i=13(u0i-c4i)2)(1-ϕ2),r2=(∑i=13(u0i-c1i)2-∑i=13(u0i-c2i)2)ϕ1+(∑i=13(u0i-c3i)2-∑i=13(u0i-c4i)2)(1-ϕ1).Besides, another difference is the computing of the edge detector functiong(|∇u0(x)|). According to the definition of g in (2.3), we need to compute |∇u0(x)|2 when computing g.Whenu0 is a gray image, for each pixel x=(x,y), |∇u0(x)|2 is defined as
(2.9)|∇u0(x)|2=[(∂u0∂x)2+(∂u0∂y)2](x).Whenu0 is a color image, |∇u0(x)|2 is defined as
(2.10)|∇u0(x)|2={[(∂u01∂x)2+(∂u02∂x)2+(∂u03∂x)2]+[(∂u01∂y)2+(∂u02∂y)2+(∂u03∂y)2]}(x).Similarly to the case of gray images, once the minimization problem (2.7) is solved, the image domain Ω can also be partitioned into four regions Ωi(i=1,2,3,4) by (2.5). Then, the average vector ci=(c1i,c2i,c3i,c4i)(i=1,2,3) for each channel can be updated by(2.11)c1i=mean(u0i)inΩ1,c2i=mean(u0i)inΩ2,c3i=mean(u0i)inΩ3,c4i=mean(u0i)inΩ4,(i=1,2,3).Letu be the fitting image of the given color image u0 and ui(i=1,2,3) the three channels of u. Then, each channel ui of the fitting image u can be expressed as
(2.12)ui=c1iϕ1ϕ2+c2iϕ1(1-ϕ2)+c3i(1-ϕ1)ϕ2+c4i(1-ϕ1)(1-ϕ2).Thus, we can obtain the fitting imageu by(2.13)u(:,:,1)=u1,u(:,:,2)=u2,u(:,:,3)=u3.We denote the four segmented phases ofu as phase1, phase2, phase3, and phase4, then the four phases can be given as(2.14)phase1(:,:,i)=c1iϕ1ϕ2,phase2(:,:,i)=c2iϕ1(1-ϕ2),phase3(:,:,i)=c3i(1-ϕ1)ϕ2,phase4(:,:,i)=c4i(1-ϕ1)(1-ϕ2),
where i=1,2,3 denotes the three channels.In fact,u can also be obtained by
(2.15)u=phase1+phase2+phase3+phase4.Remark 2.1.
The multiphase image segmentation model we proposed is different from the piecewise constant multiphase Vese-Chan model in [3] for the following reasons. First, our energy functional E(ϕ1,ϕ2) in (2.1) or (2.7) is different from [3]. Second, we incorporate information from the edge into the energy functional by using a nonnegative edge detector function g to detect boundaries more easily. Third, we apply the split Bregman method, and our model is much more efficient than the model in [3]. Last, we do not need the procedure to reinitialize ϕi to the signed distance function to its zero-level curve as [1–3], we simply initialize ϕi as a binary step function which takes 1 inside a region and 0 outside. We have restricted ϕi to [0,1], thus we do not need the Heaviside function H or its smooth versions Hε as [3].


### 2.3. Application of the Split Bregman Method to Updateϕ1 and ϕ2

To minimize the energy functionalE(ϕ1,ϕ2) in (2.1) or (2.7) with respect to ϕ1 and ϕ2, one traditional method is using the standard gradient descent method directly [1–3]. In this section, we apply the split Bregman method to solve our proposed minimization problem (2.7) more efficiently. First, we introduce two auxiliary variables: d1⃗←∇ϕ1 and d2⃗←∇ϕ2. Then, we add two quadratic penalty functions to weakly enforce the resulting equality constraints and get the unconstrained problem as follows:(2.16)(ϕ1*,ϕ2*,d1⃗*,d2⃗*)=arg⁡min⁡d1⃗,d2⃗0≤ϕ1,ϕ2≤1(|d1⃗|g+|d2⃗|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉+λ2‖d1⃗-∇ϕ1‖2+λ2‖d2⃗-∇ϕ2‖2),
where λ is a positive constant.We then strictly enforce the constraints by applying the Bregman iteration, and this results in the following optimization problem:(2.17)(ϕ1k+1,ϕ2k+1,d1⃗k+1,d2⃗k+1)=argmin⁡d1⃗,d2⃗0≤ϕ1,ϕ2≤1(|d1⃗|g+|d2⃗|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉+λ2‖d1⃗-∇ϕ1-b1⃗k‖2+λ2‖d2⃗-∇ϕ2-b2⃗k‖2),b1⃗k+1=b1⃗k+∇ϕ1k+1-d1⃗k+1,b2⃗k+1=b2⃗k+∇ϕ2k+1-d2⃗k+1.Keepd1⃗, d2⃗, and ϕ2 fixed, the Euler-Lagrange equation of the optimization problem (2.17) with respect to ϕ1 is(2.18)Δϕ1=μλr1+∇⋅(d1⃗-b1⃗),whenever0≤ϕ1≤1.Note here that from the definition ofr2 in (2.8), r2 is a function of ϕ1. Here, when we minimize (2.17) with respect to ϕ1, we just consider r2 as a constant. From the numerical results in the next section, we can see that this consideration does not affect the performance our model. Similarly, when we minimize (2.17) with respect to ϕ2 for fixed d1⃗, d2⃗, and ϕ1, we also consider r1 as a constant and get the following Euler-Lagrange equation:(2.19)Δϕ2=μλr2+∇⋅(d2⃗-b2⃗),whenever0≤ϕ2≤1.For (2.18) and (2.19), the central difference and the backward difference are used for the Laplace operator and the divergence operator, respectively. Then the resulting numerical scheme is
(2.20)αl,i,j=dl,i-1,jx-dl,i,jx+dl,i,j-1y-dl,i,jy-(bl,i-1,jx-bl,i,jx+bl,i,j-1y-bl,i,jy),βl,i,j=14(ϕl,i-1,j+ϕl,i+1,j+ϕl,i,j-1+ϕl,i,j+1-μλrl+αl,i,j),ϕl,i,j=max⁡⁡{min⁡⁡{βl,i,j,1},0},
where l=1or2 denotes the two schemes for (2.18) and (2.19).Then we minimize (2.17) with respect to d1⃗ (d2⃗) for fixed ϕ1, ϕ2, and d2⃗ (ϕ1,ϕ2 and d2⃗) and obtain
(2.21)d1⃗k+1=shrinkg(b1⃗k+∇ϕ1k+1,1λ)=shrink(b1⃗k+∇ϕ1k+1,gλ),d2⃗k+1=shrinkg(b2⃗k+∇ϕ2k+1,1λ)=shrink(b2⃗k+∇ϕ2k+1,gλ),
where
(2.22)shrink(x,γ)={x|x|max⁡(|x|-γ,0),x≠0,0,x=0.


## 2.1. Our Model for Gray Images

In our previous work [19], we have proposed an improved active contour model for multiphase image segmentation based on the piecewise constant multiphase Vese-Chan model [3], the globally convex image segmentation method [20] and the split Bregman method [15–17]. The model we proposed in [19] is mainly for gray images. For simplicity, we only review the four-phase model in this section.LetΩ⊂ℜ2 be the image domain and u0:Ω→ℜ be a given gray level image. Our proposed model is given in the following level set formulation:
(2.1)min⁡0≤ϕ1,ϕ2≤1E(ϕ1,ϕ2)=min⁡0≤ϕ1,ϕ2≤1(|∇ϕ1|g+|∇ϕ2|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉),
where ϕ1 and ϕ2 are two level set functions, and μ is a positive constant. |∇ϕi|g and 〈ϕi,ri〉 are defined as
(2.2)|∇ϕi|g=∫g(|∇u0(x)|)|∇ϕi(x)|dx,〈ϕi,ri〉=∫ϕi(x)ri(x)dx,
where
(2.3)g(ξ)=11+β|ξ|2
is an edge detector function [17, 18, 21], and β is a parameter that determines the detail level of the segmentation. r1 and r2 are defined as
(2.4)r1=((u0-c11)2-(u0-c01)2)ϕ2+((u0-c10)2-(u0-c00)2)(1-ϕ2),r2=((u0-c11)2-(u0-c10)2)ϕ1+((u0-c01)2-(u0-c00)2)(1-ϕ1),
where c=(c11,c10,c01,c00) is a constant vector of averages in the four phases.Givenϕ1(0,x)=ϕ1,0(x),ϕ2(0,x)=ϕ2,0(x), once the minimization problem (2.1) is solved, the image domain Ω can be partitioned into four regions Ωi(i=1,2,3,4) by thresholding the level set functions as follows:(2.5)Ω1={x:ϕ1(x)>α,ϕ2(x)>α},Ω2={x:ϕ1(x)>α,ϕ2(x)<α},Ω3={x:ϕ1(x)<α,ϕ2(x)>α},Ω4={x:ϕ1(x)<α,ϕ2(x)<α},
for some α∈(0,1), we choose α=0.5 in this paper. This is different from [3] where the zero level set is used to identify the boundary. We use the α=0.5 level set to identify the boundary because in our paper we have restricted 0≤ϕi≤1. The average vector c=(c11,c10,c01,c00) can be obtained by(2.6)c11=mean(u0)inΩ1,c10=mean(u0)inΩ2,c01=mean(u0)inΩ3,c00=mean(u0)inΩ4.


## 2.2. Our Model for Color Images

In this section, we extend our model in Section2.1 from gray images to color images. Let u0:Ω→ℜ3 be a given color image and u0i(i=1,2,3) the three channels of u0. Let ci=(c1i,c2i,c3i,c4i)(i=1,2,3) be the average vector of the ith channel.The extension of our model to color images is still in the same formulation as (2.1):(2.7)min⁡0≤ϕ1,ϕ2≤1E(ϕ1,ϕ2)=min⁡0≤ϕ1,ϕ2≤1(|∇ϕ1|g+|∇ϕ2|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉).The difference between (2.7) and (2.1) lies in the definition of r1 and r2. In our model for color images (2.7), r1 and r2 are defined differently from (2.4):
(2.8)r1=(∑i=13(u0i-c1i)2-∑i=13(u0i-c3i)2)ϕ2+(∑i=13(u0i-c2i)2-∑i=13(u0i-c4i)2)(1-ϕ2),r2=(∑i=13(u0i-c1i)2-∑i=13(u0i-c2i)2)ϕ1+(∑i=13(u0i-c3i)2-∑i=13(u0i-c4i)2)(1-ϕ1).Besides, another difference is the computing of the edge detector functiong(|∇u0(x)|). According to the definition of g in (2.3), we need to compute |∇u0(x)|2 when computing g.Whenu0 is a gray image, for each pixel x=(x,y), |∇u0(x)|2 is defined as
(2.9)|∇u0(x)|2=[(∂u0∂x)2+(∂u0∂y)2](x).Whenu0 is a color image, |∇u0(x)|2 is defined as
(2.10)|∇u0(x)|2={[(∂u01∂x)2+(∂u02∂x)2+(∂u03∂x)2]+[(∂u01∂y)2+(∂u02∂y)2+(∂u03∂y)2]}(x).Similarly to the case of gray images, once the minimization problem (2.7) is solved, the image domain Ω can also be partitioned into four regions Ωi(i=1,2,3,4) by (2.5). Then, the average vector ci=(c1i,c2i,c3i,c4i)(i=1,2,3) for each channel can be updated by(2.11)c1i=mean(u0i)inΩ1,c2i=mean(u0i)inΩ2,c3i=mean(u0i)inΩ3,c4i=mean(u0i)inΩ4,(i=1,2,3).Letu be the fitting image of the given color image u0 and ui(i=1,2,3) the three channels of u. Then, each channel ui of the fitting image u can be expressed as
(2.12)ui=c1iϕ1ϕ2+c2iϕ1(1-ϕ2)+c3i(1-ϕ1)ϕ2+c4i(1-ϕ1)(1-ϕ2).Thus, we can obtain the fitting imageu by(2.13)u(:,:,1)=u1,u(:,:,2)=u2,u(:,:,3)=u3.We denote the four segmented phases ofu as phase1, phase2, phase3, and phase4, then the four phases can be given as(2.14)phase1(:,:,i)=c1iϕ1ϕ2,phase2(:,:,i)=c2iϕ1(1-ϕ2),phase3(:,:,i)=c3i(1-ϕ1)ϕ2,phase4(:,:,i)=c4i(1-ϕ1)(1-ϕ2),
where i=1,2,3 denotes the three channels.In fact,u can also be obtained by
(2.15)u=phase1+phase2+phase3+phase4.Remark 2.1.
The multiphase image segmentation model we proposed is different from the piecewise constant multiphase Vese-Chan model in [3] for the following reasons. First, our energy functional E(ϕ1,ϕ2) in (2.1) or (2.7) is different from [3]. Second, we incorporate information from the edge into the energy functional by using a nonnegative edge detector function g to detect boundaries more easily. Third, we apply the split Bregman method, and our model is much more efficient than the model in [3]. Last, we do not need the procedure to reinitialize ϕi to the signed distance function to its zero-level curve as [1–3], we simply initialize ϕi as a binary step function which takes 1 inside a region and 0 outside. We have restricted ϕi to [0,1], thus we do not need the Heaviside function H or its smooth versions Hε as [3].


## 2.3. Application of the Split Bregman Method to Updateϕ1 and ϕ2

To minimize the energy functionalE(ϕ1,ϕ2) in (2.1) or (2.7) with respect to ϕ1 and ϕ2, one traditional method is using the standard gradient descent method directly [1–3]. In this section, we apply the split Bregman method to solve our proposed minimization problem (2.7) more efficiently. First, we introduce two auxiliary variables: d1⃗←∇ϕ1 and d2⃗←∇ϕ2. Then, we add two quadratic penalty functions to weakly enforce the resulting equality constraints and get the unconstrained problem as follows:(2.16)(ϕ1*,ϕ2*,d1⃗*,d2⃗*)=arg⁡min⁡d1⃗,d2⃗0≤ϕ1,ϕ2≤1(|d1⃗|g+|d2⃗|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉+λ2‖d1⃗-∇ϕ1‖2+λ2‖d2⃗-∇ϕ2‖2),
where λ is a positive constant.We then strictly enforce the constraints by applying the Bregman iteration, and this results in the following optimization problem:(2.17)(ϕ1k+1,ϕ2k+1,d1⃗k+1,d2⃗k+1)=argmin⁡d1⃗,d2⃗0≤ϕ1,ϕ2≤1(|d1⃗|g+|d2⃗|g+μ〈ϕ1,r1〉+μ〈ϕ2,r2〉+λ2‖d1⃗-∇ϕ1-b1⃗k‖2+λ2‖d2⃗-∇ϕ2-b2⃗k‖2),b1⃗k+1=b1⃗k+∇ϕ1k+1-d1⃗k+1,b2⃗k+1=b2⃗k+∇ϕ2k+1-d2⃗k+1.Keepd1⃗, d2⃗, and ϕ2 fixed, the Euler-Lagrange equation of the optimization problem (2.17) with respect to ϕ1 is(2.18)Δϕ1=μλr1+∇⋅(d1⃗-b1⃗),whenever0≤ϕ1≤1.Note here that from the definition ofr2 in (2.8), r2 is a function of ϕ1. Here, when we minimize (2.17) with respect to ϕ1, we just consider r2 as a constant. From the numerical results in the next section, we can see that this consideration does not affect the performance our model. Similarly, when we minimize (2.17) with respect to ϕ2 for fixed d1⃗, d2⃗, and ϕ1, we also consider r1 as a constant and get the following Euler-Lagrange equation:(2.19)Δϕ2=μλr2+∇⋅(d2⃗-b2⃗),whenever0≤ϕ2≤1.For (2.18) and (2.19), the central difference and the backward difference are used for the Laplace operator and the divergence operator, respectively. Then the resulting numerical scheme is
(2.20)αl,i,j=dl,i-1,jx-dl,i,jx+dl,i,j-1y-dl,i,jy-(bl,i-1,jx-bl,i,jx+bl,i,j-1y-bl,i,jy),βl,i,j=14(ϕl,i-1,j+ϕl,i+1,j+ϕl,i,j-1+ϕl,i,j+1-μλrl+αl,i,j),ϕl,i,j=max⁡⁡{min⁡⁡{βl,i,j,1},0},
where l=1or2 denotes the two schemes for (2.18) and (2.19).Then we minimize (2.17) with respect to d1⃗ (d2⃗) for fixed ϕ1, ϕ2, and d2⃗ (ϕ1,ϕ2 and d2⃗) and obtain
(2.21)d1⃗k+1=shrinkg(b1⃗k+∇ϕ1k+1,1λ)=shrink(b1⃗k+∇ϕ1k+1,gλ),d2⃗k+1=shrinkg(b2⃗k+∇ϕ2k+1,1λ)=shrink(b2⃗k+∇ϕ2k+1,gλ),
where
(2.22)shrink(x,γ)={x|x|max⁡(|x|-γ,0),x≠0,0,x=0.


## 3. The Algorithm and Numerical Results

### 3.1. The Algorithm

The split Bregman algorithm for the proposed minimization problem (2.7) in Section 2.3 can be summarized as in  Algorithm 1.Algorithm 1
1: Input u0 and ϕ10, ϕ202: Compute initial Ωi and c by (2.5) and  (2.11)3: while∥ϕ1k+1-ϕ1k∥>ϵ or ∥ϕ2k+1-ϕ2k∥>ϵ do4: Define r1k, r2k by  (2.8)5: Update ϕ1k+1, ϕ2k+1 and d1⃗k+1, d1⃗k+1 by (2.20) and  (2.21)6: Update bl⃗k+1=bl⃗k+∇ϕlk+1-dl⃗k+1 (l=1,2)7: Update Ωi and c by (2.5) and  (2.11)8: end while


### 3.2. Numerical Results

Our model has been applied to synthetic and real color images with multiple regions in this section. When working with the level set functionϕi, we simply initialize the level set function ϕi as a binary step function which takes 1 inside a region and 0 outside. The advantage of using binary step functions as the initial level set functions is that new contours can emerge easily and the curve evolution is significantly faster than the evolution from initial functions as signed distance maps. In this paper, λ=1 and β=0 are used for all images. Unless otherwise specified, we choose μ=50/2552 for all synthetic color images and μ=30/2552 for all real images. The size of each image is specified in each figure.Figure1 shows the application of our model to a synthetic color image with a triple junction. Most of the models [12, 13] need three level set functions to represent the triple junction. Here, we need only two level set functions. The α=0.5 level sets of the final ϕ1 and ϕ2 are shown in (d) and (h), which have to overlap on a segment of the triple junction.Results for a synthetic image with our model. (a)–(c): The active contour evolving process from the initial contour to the final contour; (e)–(g): The corresponding fitting imagesu at different iterations; (d) and (h): The α level sets of the final ϕ1 and ϕ2; (i)–(l): The final four phases; size = 256 × 256.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)We then apply our model to another synthetic image in Figures2 and 3. In Figure 2, we show the results for the clean image with three different initial conditions. When we seed with small initials in (a) and (b), our model can work well, both the final contours and the final fitting images are correct. However, when we use two disjoint circles as the initial curves, our model will be trapped in a local minimum, and the pink object will be missed as shown in (f) and (i). Then we add random noises with standard deviation 30.0 to the clean image. Figure 3 shows the results for the noisy version. We can observe that our model works well even if the noise level is high. In the following experiments, the noises we added are all random noise with standard deviation 30.0.Results for a synthetic image with different initial conditions; (a)–(c): The original image with three different initial contours; (d)–(f): The corresponding final contours; (g)–(i): The corresponding final fitting images; size = 256 × 256.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Application of our model to the noisy version of the image from Figure2. (a)–(c): The curve evolution process; (e)–(g): The corresponding fitting images u; (d) and (h): The α level sets of the final ϕ1 and ϕ2; (i)–(l): The final four phases.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)The results for another synthetic color image without or with noise are shown in Figures4 and 5. The results for the clean image with two different initial conditions are shown in Figure 4. When we use two disjoint initial curves as shown in (a), both the final contour and the final fitting image are good. When using two joint initial curves as shown in (g), we can get the right final contour shown in (h), but the final fitting image is incorrect shown in (i). The two bottom objects are considered to have the same average value. The results for the noisy version are shown in Figure 5.Results for a synthetic image with two different initial conditions. (a) and (g): The two different initial contours; (b) and (h): The corresponding final contours; (c) and (i): The corresponding final fitting images; (d) and (j): Theα level sets of the final ϕ1 and ϕ2 with the initial condition (a); (i)–(l): The final four phases with the initial condition (a); size = 240 × 110.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Application of our model to the noisy version of the image from Figure4. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting images u.
(a)(b)(c)(d)(e)(f)(g)(h)Figures6, 7, and 8 show the applications of our model to two synthetic images and their noisy versions. Both the two images have three objects with different shapes, which may be difficult to detect. However, our model can identify the objects correctly for both clean images and noisy versions.Results for a synthetic image without noise and with noise. (a)–(c) and (e)–(g): The initial contour, final contour and final fitting image for the clean and noisy images; (d), (h) and (i)–(l): Theα level sets of the final ϕ1, ϕ2 and the final four phases for the noisy version; size = 200 × 200.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Results for a clean synthetic image with our model. (a)–(c): The curve evolution process; (e)–(g): The corresponding fitting imagesu; (d) and (h): The α level sets of the final ϕ1 and ϕ2; (i)–(l): The final four phases; size = 256 × 256.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Application of our model to the noisy version of the image from Figure7. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting images u.
(a)(b)(c)(d)(e)(f)(g)(h)We have declared that the standard deviation of the random noise added to the above synthetic images is 30.0. Now, we increase the noise level and give the results with our model in Figures9 and 10. For the above five synthetic images, when we increase the noise level to 50.0, our model can give the correct final contours and final fittings, which can be clearly seen from Figure 9. Even the noise is increased to a much higher level 80.0, our model can still detect the object boundaries of the synthetic images and give good fitting images. The corresponding initial contours, final contours and final fitting images are shown in Figure 10. Here, we only need to choose μ=1/2552 to avoid detecting the high-level noise for the noisy versions in Figures 9 and 10. This demonstrates the robustness of our model to noise.Results for the synthetic images added by random noise with standard deviation 50.0. (a)–(e): The initial contours; (f)–(j): The final contours; (k)–(o): The final fitting images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)Results for the synthetic images added by random noise with standard deviation 80.0. (a)–(d): The initial contours; (e)–(h): The final contours; (i)–(l): The final fitting images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)We then apply our model to three real color images, and the results are shown in Figures11, 12, 13, 14, and 15. We seed with small initial curves for all of these real images. Figures 11 and 12 show the results for an image of flowers, while Figures 13 and 14 give the results for a 4-colors image. We can observe that our model can segment these color images well. In Figure 15 we apply our model to an image with contours without gradient (cognitive contours). From the final fitting image shown in (d), it can be seen that our model can also detect “contours without edges".Application of our model to a real color image of flowers. (a)–(d): The curve evolution process. (e)–(h): The corresponding fitting imagesu; size = 160 × 240.
(a)(b)(c)(d)(e)(f)(g)(h)Final four phases of the flowers image.
(a)(b)(c)(d)Application of our model to a 4-color image. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; size = 185 × 136.
(a)(b)(c)(d)(e)(f)(g)(h)Final four phases of the 4-colors image.
(a)(b)(c)(d)Application of our model to a real color image. (a): The original image; (b): The initial contour; (c): The final contour; (d): The final fitting image; size = 160 × 240.
(a)(b)(c)(d)Figures16 and 17 show the results of our model for two real color images from the Berkeley image set. In Figures 16 and 17, we give the curve evolution process in Row 1 and the corresponding fitting images in Row 2, while the final four phases are shown in Row 3. It can be observed that our model can handle these two real color images very well.Results for a real color image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 160 × 240.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Results for a real color image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 160 × 240.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Our model can be applied to segment real medical MRI images. We apply our model to two real brain MRI images in Figures18 and 19. Similarly, the curve evolution process, corresponding fitting images and final four phases are shown in Row 1, Row 2 and Row 3, respectively. We can see that our model can identify quite well the gray matter, the white matter, and so forth.Results for a brain MRI image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 120 × 117.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Results for a brain MRI image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 125 × 160.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)To demonstrate the efficiency of our model, we compare the computation (CPU) time for the above five synthetic images with the Vese-Chan model from [3] and our model in Tables 1 and 2. In Table 1, we test the CPU time for all the clean synthetic images used in this paper. Except the image from Figure 1, we use the noisy versions for other four synthetic images for Table 2. From Tables 1 and 2, we can see that our model is much more efficient than the Vese-Chan model by applying the split Bregman method. We record the CPU time from our experiments with Matlab codes run on an ACPI Multiprocessor PC, Intel(R) Core(TM)2 Quad CPU Q8200, 2.33 GHz, 2 GB RAM, with Matlab R2010a on Windows XP.Table 1
The CPU time (in second) for the Vese-Chan model and our model for the clean synthetic images from Figures1, 2, 4, 6, and 7. The sizes of images are 256 × 256, 256 × 256, 240 × 110, 200 × 200, and 256 × 256, respectively.
Figure1Figure2Figure4Figure6Figure7The Vese-Chan model107.0692106.894160.485861.9038108.2529Our model0.4887060.4879070.1983820.2030330.494109Table 2
The CPU time (in second) for the Vese-Chan model and our model for the noisy synthetic images from Figures3, 5, 6, and 8. The sizes of images are 256 × 256, 240 × 110, 200 × 200, and 256 × 256, respectively.
Figure3Figure5Figure6Figure8The Vese-Chan model228.9331131.935074.3890198.9657Our model1.0449410.4327220.2439820.908158


## 3.1. The Algorithm

The split Bregman algorithm for the proposed minimization problem (2.7) in Section 2.3 can be summarized as in  Algorithm 1.Algorithm 1
1: Input u0 and ϕ10, ϕ202: Compute initial Ωi and c by (2.5) and  (2.11)3: while∥ϕ1k+1-ϕ1k∥>ϵ or ∥ϕ2k+1-ϕ2k∥>ϵ do4: Define r1k, r2k by  (2.8)5: Update ϕ1k+1, ϕ2k+1 and d1⃗k+1, d1⃗k+1 by (2.20) and  (2.21)6: Update bl⃗k+1=bl⃗k+∇ϕlk+1-dl⃗k+1 (l=1,2)7: Update Ωi and c by (2.5) and  (2.11)8: end while


## 3.2. Numerical Results

Our model has been applied to synthetic and real color images with multiple regions in this section. When working with the level set functionϕi, we simply initialize the level set function ϕi as a binary step function which takes 1 inside a region and 0 outside. The advantage of using binary step functions as the initial level set functions is that new contours can emerge easily and the curve evolution is significantly faster than the evolution from initial functions as signed distance maps. In this paper, λ=1 and β=0 are used for all images. Unless otherwise specified, we choose μ=50/2552 for all synthetic color images and μ=30/2552 for all real images. The size of each image is specified in each figure.Figure1 shows the application of our model to a synthetic color image with a triple junction. Most of the models [12, 13] need three level set functions to represent the triple junction. Here, we need only two level set functions. The α=0.5 level sets of the final ϕ1 and ϕ2 are shown in (d) and (h), which have to overlap on a segment of the triple junction.Results for a synthetic image with our model. (a)–(c): The active contour evolving process from the initial contour to the final contour; (e)–(g): The corresponding fitting imagesu at different iterations; (d) and (h): The α level sets of the final ϕ1 and ϕ2; (i)–(l): The final four phases; size = 256 × 256.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)We then apply our model to another synthetic image in Figures2 and 3. In Figure 2, we show the results for the clean image with three different initial conditions. When we seed with small initials in (a) and (b), our model can work well, both the final contours and the final fitting images are correct. However, when we use two disjoint circles as the initial curves, our model will be trapped in a local minimum, and the pink object will be missed as shown in (f) and (i). Then we add random noises with standard deviation 30.0 to the clean image. Figure 3 shows the results for the noisy version. We can observe that our model works well even if the noise level is high. In the following experiments, the noises we added are all random noise with standard deviation 30.0.Results for a synthetic image with different initial conditions; (a)–(c): The original image with three different initial contours; (d)–(f): The corresponding final contours; (g)–(i): The corresponding final fitting images; size = 256 × 256.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Application of our model to the noisy version of the image from Figure2. (a)–(c): The curve evolution process; (e)–(g): The corresponding fitting images u; (d) and (h): The α level sets of the final ϕ1 and ϕ2; (i)–(l): The final four phases.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)The results for another synthetic color image without or with noise are shown in Figures4 and 5. The results for the clean image with two different initial conditions are shown in Figure 4. When we use two disjoint initial curves as shown in (a), both the final contour and the final fitting image are good. When using two joint initial curves as shown in (g), we can get the right final contour shown in (h), but the final fitting image is incorrect shown in (i). The two bottom objects are considered to have the same average value. The results for the noisy version are shown in Figure 5.Results for a synthetic image with two different initial conditions. (a) and (g): The two different initial contours; (b) and (h): The corresponding final contours; (c) and (i): The corresponding final fitting images; (d) and (j): Theα level sets of the final ϕ1 and ϕ2 with the initial condition (a); (i)–(l): The final four phases with the initial condition (a); size = 240 × 110.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Application of our model to the noisy version of the image from Figure4. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting images u.
(a)(b)(c)(d)(e)(f)(g)(h)Figures6, 7, and 8 show the applications of our model to two synthetic images and their noisy versions. Both the two images have three objects with different shapes, which may be difficult to detect. However, our model can identify the objects correctly for both clean images and noisy versions.Results for a synthetic image without noise and with noise. (a)–(c) and (e)–(g): The initial contour, final contour and final fitting image for the clean and noisy images; (d), (h) and (i)–(l): Theα level sets of the final ϕ1, ϕ2 and the final four phases for the noisy version; size = 200 × 200.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Results for a clean synthetic image with our model. (a)–(c): The curve evolution process; (e)–(g): The corresponding fitting imagesu; (d) and (h): The α level sets of the final ϕ1 and ϕ2; (i)–(l): The final four phases; size = 256 × 256.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Application of our model to the noisy version of the image from Figure7. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting images u.
(a)(b)(c)(d)(e)(f)(g)(h)We have declared that the standard deviation of the random noise added to the above synthetic images is 30.0. Now, we increase the noise level and give the results with our model in Figures9 and 10. For the above five synthetic images, when we increase the noise level to 50.0, our model can give the correct final contours and final fittings, which can be clearly seen from Figure 9. Even the noise is increased to a much higher level 80.0, our model can still detect the object boundaries of the synthetic images and give good fitting images. The corresponding initial contours, final contours and final fitting images are shown in Figure 10. Here, we only need to choose μ=1/2552 to avoid detecting the high-level noise for the noisy versions in Figures 9 and 10. This demonstrates the robustness of our model to noise.Results for the synthetic images added by random noise with standard deviation 50.0. (a)–(e): The initial contours; (f)–(j): The final contours; (k)–(o): The final fitting images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)Results for the synthetic images added by random noise with standard deviation 80.0. (a)–(d): The initial contours; (e)–(h): The final contours; (i)–(l): The final fitting images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)We then apply our model to three real color images, and the results are shown in Figures11, 12, 13, 14, and 15. We seed with small initial curves for all of these real images. Figures 11 and 12 show the results for an image of flowers, while Figures 13 and 14 give the results for a 4-colors image. We can observe that our model can segment these color images well. In Figure 15 we apply our model to an image with contours without gradient (cognitive contours). From the final fitting image shown in (d), it can be seen that our model can also detect “contours without edges".Application of our model to a real color image of flowers. (a)–(d): The curve evolution process. (e)–(h): The corresponding fitting imagesu; size = 160 × 240.
(a)(b)(c)(d)(e)(f)(g)(h)Final four phases of the flowers image.
(a)(b)(c)(d)Application of our model to a 4-color image. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; size = 185 × 136.
(a)(b)(c)(d)(e)(f)(g)(h)Final four phases of the 4-colors image.
(a)(b)(c)(d)Application of our model to a real color image. (a): The original image; (b): The initial contour; (c): The final contour; (d): The final fitting image; size = 160 × 240.
(a)(b)(c)(d)Figures16 and 17 show the results of our model for two real color images from the Berkeley image set. In Figures 16 and 17, we give the curve evolution process in Row 1 and the corresponding fitting images in Row 2, while the final four phases are shown in Row 3. It can be observed that our model can handle these two real color images very well.Results for a real color image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 160 × 240.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Results for a real color image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 160 × 240.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Our model can be applied to segment real medical MRI images. We apply our model to two real brain MRI images in Figures18 and 19. Similarly, the curve evolution process, corresponding fitting images and final four phases are shown in Row 1, Row 2 and Row 3, respectively. We can see that our model can identify quite well the gray matter, the white matter, and so forth.Results for a brain MRI image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 120 × 117.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Results for a brain MRI image with our model. (a)–(d): The curve evolution process; (e)–(h): The corresponding fitting imagesu; (i)–(l): The final four phases; size = 125 × 160.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)To demonstrate the efficiency of our model, we compare the computation (CPU) time for the above five synthetic images with the Vese-Chan model from [3] and our model in Tables 1 and 2. In Table 1, we test the CPU time for all the clean synthetic images used in this paper. Except the image from Figure 1, we use the noisy versions for other four synthetic images for Table 2. From Tables 1 and 2, we can see that our model is much more efficient than the Vese-Chan model by applying the split Bregman method. We record the CPU time from our experiments with Matlab codes run on an ACPI Multiprocessor PC, Intel(R) Core(TM)2 Quad CPU Q8200, 2.33 GHz, 2 GB RAM, with Matlab R2010a on Windows XP.Table 1
The CPU time (in second) for the Vese-Chan model and our model for the clean synthetic images from Figures1, 2, 4, 6, and 7. The sizes of images are 256 × 256, 256 × 256, 240 × 110, 200 × 200, and 256 × 256, respectively.
Figure1Figure2Figure4Figure6Figure7The Vese-Chan model107.0692106.894160.485861.9038108.2529Our model0.4887060.4879070.1983820.2030330.494109Table 2
The CPU time (in second) for the Vese-Chan model and our model for the noisy synthetic images from Figures3, 5, 6, and 8. The sizes of images are 256 × 256, 240 × 110, 200 × 200, and 256 × 256, respectively.
Figure3Figure5Figure6Figure8The Vese-Chan model228.9331131.935074.3890198.9657Our model1.0449410.4327220.2439820.908158


## 4. Conclusion

In this paper, we propose a new and fast multiphase image segmentation model for color images. Our model is an extension of the model for gray images we proposed in [19]. We have tested our model with many synthetic and real color images. Numerical results show that our model has all the benefits of the Vese-Chan model [3], including robustness even with noisy data and automatic detection of interior contours. However, our model is much more efficient by applying the split Bregman method. Of course, our model has its limitation. Our model is a piecewise constant multiphase segmentation model, and it mainly focuses on homogeneous multiphase images. It only considers the global information of the given image, thus it cannot deal with images with inhomogeneity. In future work, we will extend our model to handle inhomogeneous multiphase images by taking both local and global information into consideration.



---

*Source: 494761-2012-03-22.xml*