# Fast Global Minimization of the Chan–Vese Model for Image Segmentation Problem

Authors: Ran Gao; Li-Zhen Guo
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

The segmentation of weak boundary is still a difficult problem, especially sensitive to noise, which leads to the failure of segmentation. Based on the previous works, by adding the boundary indicator function withL2,1 norm, a new convergent variational model is proposed. A novel strategy for the weak boundary image is presented. The existence of the minimizer for our model is given, by using the alternating direction method of multipliers (ADMMs) to solve the model. The experiments show that our new method is robust in segmentation of objects in a range of images with noise, low contrast, and direction.


---

## Body

## 1. Introduction

Image segmentation [1–6] is the process of separating objects of interest from each other or backwards to find boundaries of objects. It has become increasingly important in the last decade, due to being a fast-expanding field of applications in image analysis and computer vision. Not only that, it is a fundamental problem in the field of computer vision, because recognition and reconstruction often rely on this information [7, 8]. There are many circumstances affecting the segmentation results such as the noise, weak edge, and directional texture information, and these effects lead the segmentation problem to be a typical structural ill-posed problem. Many approaches exist to deal with this problem including histogram analysis, region growth, and edge detection [9, 10]. Up to now, although there are many existing segmentation methods, they still lack the universality. Over the last decade, researchers are constantly exploring new segmentation methods to make image segmentation results as good as possible and variational methods and partial differential equation (PDE) base techniques have been introduced to image segmentation. In this paper, we present a classical method of variation for these problems. Specially, these methods are divided into three methods: threshold-based segmentation [11–13], region-based segmentation [8, 14], and edge-based segmentation [7, 15, 16]. Here, we mainly consider the region-based method and edge-based method, which are uniformed with the most common segmentation method.The edge-based models such as the GAC model [15–18] mainly utilize the gradient information to direct the contours toward the boundaries of desired objects and then obtain the segmentation image. Therefore, these classes of models are very sensitive to the noise and also show difficulty to detect the weak boundaries. That is to say, when there is strong noise in the image, the boundary between the target and the background must not be obvious enough and the segmentation will fail. In order to overcome this drawback, one way to process a noisy image is by adding a smoothing step prior to segmentation, but doing this also smooths image edges. Moreover, the segmentation result generated by these models is highly dependent on the initial contour placement due to the nonconvexity of the proposed model. For the region-based models [8, 14], they mainly incorporate region information so that the image within each region has uniform characteristics such as intensities and textures. Then, these models are more suitable for segmentation of the image with weak boundaries or noise. Furthermore, they are also less sensitive to the location of the initial contour. The well-known Mumford–Shah (MS) variational model [7] can achieve both goals simultaneously by using a piecewise smooth representation of an image [19, 20] and separate different regions by using a set of close contours (curves in 2D and surfaces in 3D). However, due to an unknown set of boundaries and its nonconvexity, the numerical experiment of this model has made very hard. With the assumptions of the segmentation regions to be piecewise constant, the Mumford–Shah model can be simplified as the Chan–Vese (CV) model [8, 21, 22]:(1)minc1,c2,Γλ1∫insideCI0−c12dx+λ2∫outsideCI0−c22dx+ν⋅lengthΓ,where ν, λ1, and λ2 are positive parameters, c1 and c2 are the intensity averages of the original image I0 inside and outside the contour Γ, respectively, dx is the region element, and the term of length is to regularize the contour C. Original CV model (1) only considers to segment two phases as the foreground and the background and then it extends to the multiphase image segmentation problem [23]. Obviously, the CV model can generally get satisfactory results for the images with intensity homogeneity. Simultaneously, the level set method was also proposed to solve problem (1) in [8]. However, due to the nonconvexity of problem (1), we only get the local minima which may lead to wrong levels of detail and scale and then show the sensitivity to the placement of initial contour. Hence, some global minimization active contour models had been proposed to avoid this problem for obtaining the segmentation regions when fixing intensity averages c1 and c2. Chan et al. [14] proposed a convex relaxation method of the CV model. Bresson et al. [24] further established theorems to determine the existence of global minimization of the active contour/snake model. Wang et al. [25] proposed a novel global minimization hybrid active contour model. By designing a convex energy functional and the dual algorithm, Xu et al. [26] proposed a global and local active contour model.Motivated by the Euler–Lagrange equation of problem (1) by introducing the level set function, Chan et al. [14] transformed to consider the following global segmentation model as(2)min0≤u≤1∫Ω∇udx︸≔TVu+λ∫ΩI0−c12−I0−c22udx,where λ is an arbitrary positive parameter, Ω⊂RN is an open set representing the image domain, and TVu is the total variation norm of the function u. Although the model has the ability of global segmentation, the segmentation for gray-scale inhomogeneous images is invalid. Later, a large number of scholars made a number of improvements to the model, but there are still some problems in the segmentation of image gray-scale unevenness. Recently, Bresson et al. [24] proposed to determine a global minimum of the snake model by enhancing model (2). The enhancement is realized by unifying it with the classic GAC model [27], and then, they considered the following model:(3)minu∈0,1∫Ωgx∇udx︸≔TVgu+λ∫ΩI0−c12−I0−c22udx,where gx is the edge indicator function as defined in the above. Obviously, compared with model (2), model (3) segments the image more robustly due to the term gx. In order to calculate the problem easily, similar to the case in [28], Bresson et al. [24] transform the above constraint problem into an unconstrained problem at Theorem 3, and it needs to compute the minimization of the function:(4)Ec1,c2,u=∫Ωgx∇udx+λ∫ΩI0−c12−I0−c22u+ανudx,where νξ≔max0,2ξ−1/2−1 is an exact penalty function provided that the constant α is chosen large enough compared to λ such as α>λ/2I0−c12−I0−c22L∞Ω. Because problem (4) is unconstrained, it includes a nonsmoothing term and a nonlinear term which lead to numerical difficulties by directly employing numerical methods. One efficient method is to introduce a constrained variable v≔u and then get the following form based on the penalty method:(5)Ec1,c2,u,v=∫Ωgx∇udx+12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx,where θ>0 is the penalty parameter. Furthermore, in order to overcome numerical difficulties generated by the L1-norm term, an alternative minimization scheme was introduced by adding some auxiliary variables in [24]. However, the convergence analysis of this scheme is lost.In this section, based on (5), let Ω∈Rn be a bounded open set with Lipschitz boundary and u be a given image. We propose a new model based on Bresson’s method [24] as follows:(6)Ec1,c2,u,v=∫ΩG∇u2,1dx+12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx,where I0, c1, and c2 are defined in (1) and νv is defined in (4). Also,(7)G=gx00gy,ξ2,1=ξx2+ξy2,∇g=gx,gy.The minimizing energy functionEc1,c2,u,v is as follows:(8)minc1,c2,u,vEc1,c2,u,v.Our novel contribution is to put an edge indication function into the L2,1 norm, and a novel strategy for the image with weak boundaries is present. The existence of the minimizer for (8) is given as follows.The alterative scheme to solve the minimizer model (8) is as follows:(1)
Choose the original value forc10,c20, u0, and v0, and set n≔0.(2)
When the stop criterion is not satisfied, compute the following iteration:(9)argminc1λ∫ΩI0−c12vdx,argminc2λ∫Ω−I0−c22vdx,argminuG∇u22+12θu−vn−122,argminv12θun−v22+λ∫Ωr1x,c1,c2v+αvvdx.The rest of this paper is organized as follows: In Section2, the general framework of the proposed method and the concrete solution to each subproblem are presented. In Section 3, numerical algorithm and experimental results to illustrate the effectiveness of our model in image segmentation are given. Finally, we conclude the paper in Section 4.


## 2. The Weighted Chan–Vese Model

In this section, we mainly show the related definitions and properties and solve every subproblem in model (8).


### 2.1. Basic Notations

Without loss of generality, we present a gray image as anN×N matrix, that is, Ω=1,2,…,N×1,2,…,N. The Euclidean space RN×N is denoted as W, let H:W×W, and we define the gradient operator ∇:W⟶H, for ∀v∈W, where ∇v denotes the gradient operator in the discrete context.(10)u2,1=∑i=1N∑j=1Nui,j21/2,u2=∑i=1N∑j=1Nui,j21/2.LetX=RN×N be a finite-dimensional vector space equipped with a standard scalar product for u,v∈X:(11)u,vX=∫Ωuvdx,v,vX=v2,1.For discretization of the gradient operator∇u, we use standard the finite difference operators with periodic boundary condition as follows:(12)∇ui,j=∇x+ui,j∇y+ui,j,where(13)∇x+ui,j=ui+1,j−ui,j,if1≤i<N,1≤j≤N,u1,j−ui,j,ifi=N,1≤j≤N,∇y+ui,j=ui,j+1−ui,j,if1≤i≤N,1≤j<N,ui,1−ui,j,if1≤i≤N,j=N,∇x−ui,j=ui,j−ui−1,j,if1<i≤N,1≤j≤N,ui,j−uN,j,ifi=1,1≤j≤N,∇y−ui,j=ui,j+1−ui,j,if1≤i≤N,1<j≤N,ui,j−ui,N,if1≤i≤N,j=1.Furthermore, we define the discrete divergence operatordivp:H⟶W by using the divergence theorem [27].(14)∇u,pH=−u,divpW,for ∀u∈RN×N, where p∈RN×N×RN×N and div denotes the adjoint operator of ∇; then, we have(15)divpi,j=∇x−pi,j+∇y−pi,j,∇ui,j=∇x+i,j2+∇y+i,j2.


### 2.2. The Minimization Problem

In this section, we first review some definitions and facts which will be used in our convergence analysis Algorithm1. These related contents can be found in the books [29, 30]. We use ADMM to solve every subproblem.Algorithm 1: The alterative scheme to solve the minimizer models (19) and (20).
(1)Choose the original valueu0 and v0, and set n≔0.(2)
While the stop criterion is not satisfied, compute the following iteration:argminuG∇u22+1/2θu−vn−122argminv1/2θun−v22+λ∫Ωr1x,c1,c2v+αvvdx.


#### 2.2.1. Energy Minimization with respect toc1 and c2

We fixedu and v in order to determine an exact solution of c1 and c2; then, the minimization of Ec1,c2,u,v with respect to c1 and c2 is as follows:(16)argminc1λ∫ΩI0−c12vdx,argminc2λ∫Ω−I0−c22vdx.Then, we can obtain(17)c1=c2=∫I0∫I.


#### 2.2.2. Energy Minimization with respect tou and v

For any givenc1 and c2, we have the problem(18)argminu,v≔∫ΩG∇u2,1dx+12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx.Problem (18) is convex and is not strictly convex about u and v, so its local minimizer is also the global minimizer. One of the effective methods is to change it into an unconstrained minimization problem according to the theorems in [14, 24].Obviously, it is easy to find that the objective function in problem (18) is strictly convex and coercive, so it owns a unique solution u∗,v∗. Formally, problem (18) including two variables u and v as a coupled problem, and the classical method is to divide it into two subproblems. Thus, we alteratively consider the following minimization problem:(19)Sv≔u=minuG∇u22+12θu−v22,(20)Tu≔v=minv12θu−v22+λ∫ΩI0−c12−I0−c22v+αvvdx.Theorem 1.
The following assertions hold:(1)
The solution of problem (19) is given by(21)u=v−PθSv,whereP is the project operator and S is defined by(22)S=closuredivξx:ξx∈Cc1Ω,R2,ξx≤gx,∀x∈Ω.(2)
The solution of problem (17) satisfies(23)v=minmaxu−θλr1x,c1,c2,0,1,wherer1x,c1,c2=I0−c12−I0−c22.Remark 1.
Obviously, problem (19) is a generated ROF model by adding a weighted function gx to the regularization term and so we can get the solution based on the scheme in [27,31,32]. For problem (20), it is similar to what happens in [28] for this circumstance; we notice that it is equivalent to the following minimization problem:(24)minv12θu−v22+λ∫Ωr1x,c1,c2u+δCvdx,where δ is the indicator function defined by(25)δx=0,ifx∈C,+∞,ifx∉C,with C=0,1. Thus, we can easily find that the solution of problem (20) satisfies condition (23).Definition 1.
The proximal operator of a proper, convex, semicontinuous functionφ on Rn is defined for x∈Rn by(26)proxγφ∘Bx≔argminy∈Rn12γy−x22+φBy,where B:Rn⟶Rn is a bounded linear operator. Furthermore, when B=I, it is the classical proximal operator defined in [33].Definition 2.
Letφ be a real-value convex function on Rn. The subdifferential of φ at x∈Rn is defined by(27)∂φx≔y:y∈Rnandφz≥φx+y,z−x∀z∈Rn.
The elements in∂φx are called subgradients. Especially, ∂φx=∇φx when φx is differential.Definition 3.
An operatorQ:Rn⟶Rn is firmly nonexpansive if it satisfies one of the following equivalent conditions:(i)
Qx−Qy22≤Qx−Qy,x−y,∀x,y∈Rn×Rn(ii)
Qx−Qy22≤x−y22−I−Qx−I−Qy22,∀x,y∈Rn×RnDefinition 4.
An operatorQ:Rn⟶Rn is nonexpansive if for any x,y∈Rn×Rn, we have(28)Qx−Qy2≤x−y2.It follows that a firmly nonexpansive operatorP is nonexpansive. Based on the proof of the proximal operator in [34], we obtain the similar result of the proximal operator defined in Definition 1.Lemma 1.
The proximal operatorproxγφ∘B defined in Definition 1 is firmly nonexpansive.Proof.
Following the definition of the proximity in Definition1 and the chain rule, we can get(29)φB∘proxγφ∘By−φB∘proxγφ∘Bx≥1γx−proxγφ∘Bx,proxγφ∘By−proxγφ∘Bx,φB∘proxγφ∘Bx−φB∘proxγφ∘By≥1γy−proxγφ∘By,proxγφ∘Bx−proxγφ∘By.We obtain(30)proxγφ∘Bx−proxγφ∘By22≤proxγφ∘Bx−proxγφ∘By,x−y,which implies the nonexpansion of the operator proxγφ∘B.Corollary 1.
The operatorS defined in minimization problem (19) is nonexpansive.Proof.
Based on above definitions, setB=gx∇u; then, the operator T is the proximal operator of ∫Ωgx∇udx. Furthermore, we can get S=pro⋅L1∘B. According to Lemma 1, the operator S is nonexpansive.Lemma 2.
Assume thatC⊂Rn is a closed convex set and PC⋅:Rn⟶C is a projection operator; then, the projector PC⋅ is firmly nonexpansive.Proof.
Letx∈Rn and ∀y∈C; we have(31)y−PCx,x−PCx≤0.
SincePCy=y while y∈C, we can get(32)PCy−PCx,x−y+PCy−PCx≤0.
This yields(33)PCx−PCy22≤x−y,PCx−PCy≤0.Corollary 2.
The operatorT defined in minimization problem (20) is firmly nonexpansive.Proof.
Obviously, from Theorem1 and Remark 1, we derive that(34)Tu≔v=P0,1u−θλr1x,c1,c2.
It follows from Lemma2 that T is firmly nonexpansive.In view of the above results, we know thatℛ1⋅=ST⋅ and ℛ2⋅=TS⋅ are firmly nonexpansive.Lemma 3.
Letuk and vk be generated by (2) and (3), respectively. Then, ∑k=0∞uk−uk+1L22 and ∑k=0∞vk−vk+1L22 are convergent.Proof.
DenoteE1u,v=u−vL22 and E2v=∇xvL1+∇yvL1. We have(35)Eu,v=∫Ωu+ge−udxdy+β1E1u,v+β2E2v.
Hence,(36)Euk+1,vk−Euk+1,vk+1=β1E1uk+1,vk−E1uk+1,vk+1+β2E2vk−E2vk+1.
Consider the Taylor series expansion ofE1u,v in the second variable, i.e.,(37)E1uk+1,vk=E1uk+1,vk+1+vk−vk+1T∂E1∂vuk+1,vk+1+12vk−vk+1T∂2E1∂v2uk+1,vk+1vk−vk+1.
Here,xT denotes the transpose of x. We notice that E1 is quadratic about v. Then, ∂2E1/∂v2≡2I where I is the identity matrix. Moreover, since E2 is a convex function, we get(38)E2vk≥E2vk+1+vk−vk+1T∂E2∂vvk+1.
Combining (36), (37), and (38), we obtain(39)Euk+1,vk−Euk+1,vk+1≥uk−uk+1Tβ1∂E1∂vuk+1,vk+1+β2∂E2∂vvk+1+β1vk−vk+1L22.
The subdifferential ofE with respect to v is equivalent to the vector sum of the subdifferential of E1 and E2 about v, i.e.,(40)∂E∂v=β1∂E1∂v+β2∂E2∂v.
Sincevk+1 is the minimizer of Euk+1,v, we have(41)∂E∂vuk+1,vk+1=0,that is(42)β1∂E1∂vuk+1,vk+1+β2∂E2∂vzk+1,vk+1=0.
Therefore, the first term in the right hand side of (39) is zero. When we solve successive minimization problems (19) and (20), we note that Euk,vk≥Euk+1,vk. Hence, we get(43)Euk,vk−Euk+1,vk+1≥Euk+1,vk−Euk+1,vk+1≥β1vk−vk+1L22.
It follows that the partial sum of the sequencevk−vk+1L22k=0∞ is bounded. Thus, the infinite series ∑k=0∞vk−vk+1L22 is convergent.
LetE3u=∫Ωu+ge−udxdy. By considering Euk,vk−Euk+1,vk and using the similar method, we can prove that ∑k=0∞uk−uk+1L22 is convergent.Definition 5.
A operatorP:Rn⟶Rn is asymptotically regular if for any x in Rn, the sequence Pk+1x−Pkx tends to zero as k⟶∞.Based on Lemma3, we get the following result.Lemma 4.
For any initial valuesv0andu0, assume uk and vk are generated by (19) and (20), respectively; then, T1 and T2 are asymptotically regular.Proof.
According to Lemma3, we obtain(44)limk⟶∞vk−vk+1L22=0,limk⟶∞uk−uk+1L22=0.
Asvk=T1vk−1 and uk=T2uk−1, by using the recurrence method, we have(45)vk=T1kv0,uk=T2ku0.
Therefore, we get(46)limk⟶∞T1k+1v0−T1kv0L2=0,limk⟶∞T2k+1u0−T2ku0L2=0.
This indicates thatT1 and T2 are asymptotically regular.Lemma 5.
Suppose the unique minimizer ofEu,v is u∗,v∗. Then, u∗ and v∗ are the unique fixed points of T2 and T1, respectively.Proof.
SinceEu,v is differentiable with respect to u and v separately, we obtain(47)∂E∂uu∗,v∗∂E∂vu∗,v∗=00.
This implies that(48)u∗=ℛv∗=argminuEu,v∗,v∗=Su∗=argminvEu∗,v.
We easily getu∗=ℛv∗=T2u∗ and v∗=Su∗=T1v∗. Therefore, u∗ and v∗ are the corresponding fixed points of T1 and T2.On the other hand, we note thatEu,v is strictly convex and differentiable about u and v, respectively. Therefore, the fixed points of T1 and T2 are the minimizers of Eu,v. By virtue of the uniqueness of the minimizer of Eu,v, T1 and T2 show a unique fixed point separately. That is to say, that u∗ and v∗ are the unique fixed points of T2 and T1, respectively.According to Theorem 1 in [7], we get the following result.Theorem 2.
The sequenceun,vn generated by Algorithm 1 converges to the solution u∗,v∗ of problems (19) and (20).Proof.
We should first notice the fact that(49)vn=minvFun−1,vun=minuFu,vn,so we can get(50)Fun,vn≤Fun−1,vn≤Fun−1,vn−1.
In particular, the sequenceFun,vn is nonincreasing.Theorem 3.
For any initial valuesv0andu0, assume uk and vk are generated by (2) and (3), respectively; then ukandvk converge to the corresponding fixed points of T2 and T1, i.e., uk,vk converges to u∗,v∗ which is the unique minimizer of Eu,v, as k⟶∞.Proof.
In view of Lemmas1–5, we know that T1:Rn⟶RnandT2:Rn⟶Rn are nonexpansive asymptotically regular mappings and have fixed points. By using Theorem 1 in [7], we get that ukandvk converge, respectively, to the fixed points of T2 and T1, i.e., uk,vk converges to u∗,v∗ which is the unique minimizer of Eu,v, as k⟶∞.


## 2.1. Basic Notations

Without loss of generality, we present a gray image as anN×N matrix, that is, Ω=1,2,…,N×1,2,…,N. The Euclidean space RN×N is denoted as W, let H:W×W, and we define the gradient operator ∇:W⟶H, for ∀v∈W, where ∇v denotes the gradient operator in the discrete context.(10)u2,1=∑i=1N∑j=1Nui,j21/2,u2=∑i=1N∑j=1Nui,j21/2.LetX=RN×N be a finite-dimensional vector space equipped with a standard scalar product for u,v∈X:(11)u,vX=∫Ωuvdx,v,vX=v2,1.For discretization of the gradient operator∇u, we use standard the finite difference operators with periodic boundary condition as follows:(12)∇ui,j=∇x+ui,j∇y+ui,j,where(13)∇x+ui,j=ui+1,j−ui,j,if1≤i<N,1≤j≤N,u1,j−ui,j,ifi=N,1≤j≤N,∇y+ui,j=ui,j+1−ui,j,if1≤i≤N,1≤j<N,ui,1−ui,j,if1≤i≤N,j=N,∇x−ui,j=ui,j−ui−1,j,if1<i≤N,1≤j≤N,ui,j−uN,j,ifi=1,1≤j≤N,∇y−ui,j=ui,j+1−ui,j,if1≤i≤N,1<j≤N,ui,j−ui,N,if1≤i≤N,j=1.Furthermore, we define the discrete divergence operatordivp:H⟶W by using the divergence theorem [27].(14)∇u,pH=−u,divpW,for ∀u∈RN×N, where p∈RN×N×RN×N and div denotes the adjoint operator of ∇; then, we have(15)divpi,j=∇x−pi,j+∇y−pi,j,∇ui,j=∇x+i,j2+∇y+i,j2.


## 2.2. The Minimization Problem

In this section, we first review some definitions and facts which will be used in our convergence analysis Algorithm1. These related contents can be found in the books [29, 30]. We use ADMM to solve every subproblem.Algorithm 1: The alterative scheme to solve the minimizer models (19) and (20).
(1)Choose the original valueu0 and v0, and set n≔0.(2)
While the stop criterion is not satisfied, compute the following iteration:argminuG∇u22+1/2θu−vn−122argminv1/2θun−v22+λ∫Ωr1x,c1,c2v+αvvdx.


### 2.2.1. Energy Minimization with respect toc1 and c2

We fixedu and v in order to determine an exact solution of c1 and c2; then, the minimization of Ec1,c2,u,v with respect to c1 and c2 is as follows:(16)argminc1λ∫ΩI0−c12vdx,argminc2λ∫Ω−I0−c22vdx.Then, we can obtain(17)c1=c2=∫I0∫I.


### 2.2.2. Energy Minimization with respect tou and v

For any givenc1 and c2, we have the problem(18)argminu,v≔∫ΩG∇u2,1dx+12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx.Problem (18) is convex and is not strictly convex about u and v, so its local minimizer is also the global minimizer. One of the effective methods is to change it into an unconstrained minimization problem according to the theorems in [14, 24].Obviously, it is easy to find that the objective function in problem (18) is strictly convex and coercive, so it owns a unique solution u∗,v∗. Formally, problem (18) including two variables u and v as a coupled problem, and the classical method is to divide it into two subproblems. Thus, we alteratively consider the following minimization problem:(19)Sv≔u=minuG∇u22+12θu−v22,(20)Tu≔v=minv12θu−v22+λ∫ΩI0−c12−I0−c22v+αvvdx.Theorem 1.
The following assertions hold:(1)
The solution of problem (19) is given by(21)u=v−PθSv,whereP is the project operator and S is defined by(22)S=closuredivξx:ξx∈Cc1Ω,R2,ξx≤gx,∀x∈Ω.(2)
The solution of problem (17) satisfies(23)v=minmaxu−θλr1x,c1,c2,0,1,wherer1x,c1,c2=I0−c12−I0−c22.Remark 1.
Obviously, problem (19) is a generated ROF model by adding a weighted function gx to the regularization term and so we can get the solution based on the scheme in [27,31,32]. For problem (20), it is similar to what happens in [28] for this circumstance; we notice that it is equivalent to the following minimization problem:(24)minv12θu−v22+λ∫Ωr1x,c1,c2u+δCvdx,where δ is the indicator function defined by(25)δx=0,ifx∈C,+∞,ifx∉C,with C=0,1. Thus, we can easily find that the solution of problem (20) satisfies condition (23).Definition 1.
The proximal operator of a proper, convex, semicontinuous functionφ on Rn is defined for x∈Rn by(26)proxγφ∘Bx≔argminy∈Rn12γy−x22+φBy,where B:Rn⟶Rn is a bounded linear operator. Furthermore, when B=I, it is the classical proximal operator defined in [33].Definition 2.
Letφ be a real-value convex function on Rn. The subdifferential of φ at x∈Rn is defined by(27)∂φx≔y:y∈Rnandφz≥φx+y,z−x∀z∈Rn.
The elements in∂φx are called subgradients. Especially, ∂φx=∇φx when φx is differential.Definition 3.
An operatorQ:Rn⟶Rn is firmly nonexpansive if it satisfies one of the following equivalent conditions:(i)
Qx−Qy22≤Qx−Qy,x−y,∀x,y∈Rn×Rn(ii)
Qx−Qy22≤x−y22−I−Qx−I−Qy22,∀x,y∈Rn×RnDefinition 4.
An operatorQ:Rn⟶Rn is nonexpansive if for any x,y∈Rn×Rn, we have(28)Qx−Qy2≤x−y2.It follows that a firmly nonexpansive operatorP is nonexpansive. Based on the proof of the proximal operator in [34], we obtain the similar result of the proximal operator defined in Definition 1.Lemma 1.
The proximal operatorproxγφ∘B defined in Definition 1 is firmly nonexpansive.Proof.
Following the definition of the proximity in Definition1 and the chain rule, we can get(29)φB∘proxγφ∘By−φB∘proxγφ∘Bx≥1γx−proxγφ∘Bx,proxγφ∘By−proxγφ∘Bx,φB∘proxγφ∘Bx−φB∘proxγφ∘By≥1γy−proxγφ∘By,proxγφ∘Bx−proxγφ∘By.We obtain(30)proxγφ∘Bx−proxγφ∘By22≤proxγφ∘Bx−proxγφ∘By,x−y,which implies the nonexpansion of the operator proxγφ∘B.Corollary 1.
The operatorS defined in minimization problem (19) is nonexpansive.Proof.
Based on above definitions, setB=gx∇u; then, the operator T is the proximal operator of ∫Ωgx∇udx. Furthermore, we can get S=pro⋅L1∘B. According to Lemma 1, the operator S is nonexpansive.Lemma 2.
Assume thatC⊂Rn is a closed convex set and PC⋅:Rn⟶C is a projection operator; then, the projector PC⋅ is firmly nonexpansive.Proof.
Letx∈Rn and ∀y∈C; we have(31)y−PCx,x−PCx≤0.
SincePCy=y while y∈C, we can get(32)PCy−PCx,x−y+PCy−PCx≤0.
This yields(33)PCx−PCy22≤x−y,PCx−PCy≤0.Corollary 2.
The operatorT defined in minimization problem (20) is firmly nonexpansive.Proof.
Obviously, from Theorem1 and Remark 1, we derive that(34)Tu≔v=P0,1u−θλr1x,c1,c2.
It follows from Lemma2 that T is firmly nonexpansive.In view of the above results, we know thatℛ1⋅=ST⋅ and ℛ2⋅=TS⋅ are firmly nonexpansive.Lemma 3.
Letuk and vk be generated by (2) and (3), respectively. Then, ∑k=0∞uk−uk+1L22 and ∑k=0∞vk−vk+1L22 are convergent.Proof.
DenoteE1u,v=u−vL22 and E2v=∇xvL1+∇yvL1. We have(35)Eu,v=∫Ωu+ge−udxdy+β1E1u,v+β2E2v.
Hence,(36)Euk+1,vk−Euk+1,vk+1=β1E1uk+1,vk−E1uk+1,vk+1+β2E2vk−E2vk+1.
Consider the Taylor series expansion ofE1u,v in the second variable, i.e.,(37)E1uk+1,vk=E1uk+1,vk+1+vk−vk+1T∂E1∂vuk+1,vk+1+12vk−vk+1T∂2E1∂v2uk+1,vk+1vk−vk+1.
Here,xT denotes the transpose of x. We notice that E1 is quadratic about v. Then, ∂2E1/∂v2≡2I where I is the identity matrix. Moreover, since E2 is a convex function, we get(38)E2vk≥E2vk+1+vk−vk+1T∂E2∂vvk+1.
Combining (36), (37), and (38), we obtain(39)Euk+1,vk−Euk+1,vk+1≥uk−uk+1Tβ1∂E1∂vuk+1,vk+1+β2∂E2∂vvk+1+β1vk−vk+1L22.
The subdifferential ofE with respect to v is equivalent to the vector sum of the subdifferential of E1 and E2 about v, i.e.,(40)∂E∂v=β1∂E1∂v+β2∂E2∂v.
Sincevk+1 is the minimizer of Euk+1,v, we have(41)∂E∂vuk+1,vk+1=0,that is(42)β1∂E1∂vuk+1,vk+1+β2∂E2∂vzk+1,vk+1=0.
Therefore, the first term in the right hand side of (39) is zero. When we solve successive minimization problems (19) and (20), we note that Euk,vk≥Euk+1,vk. Hence, we get(43)Euk,vk−Euk+1,vk+1≥Euk+1,vk−Euk+1,vk+1≥β1vk−vk+1L22.
It follows that the partial sum of the sequencevk−vk+1L22k=0∞ is bounded. Thus, the infinite series ∑k=0∞vk−vk+1L22 is convergent.
LetE3u=∫Ωu+ge−udxdy. By considering Euk,vk−Euk+1,vk and using the similar method, we can prove that ∑k=0∞uk−uk+1L22 is convergent.Definition 5.
A operatorP:Rn⟶Rn is asymptotically regular if for any x in Rn, the sequence Pk+1x−Pkx tends to zero as k⟶∞.Based on Lemma3, we get the following result.Lemma 4.
For any initial valuesv0andu0, assume uk and vk are generated by (19) and (20), respectively; then, T1 and T2 are asymptotically regular.Proof.
According to Lemma3, we obtain(44)limk⟶∞vk−vk+1L22=0,limk⟶∞uk−uk+1L22=0.
Asvk=T1vk−1 and uk=T2uk−1, by using the recurrence method, we have(45)vk=T1kv0,uk=T2ku0.
Therefore, we get(46)limk⟶∞T1k+1v0−T1kv0L2=0,limk⟶∞T2k+1u0−T2ku0L2=0.
This indicates thatT1 and T2 are asymptotically regular.Lemma 5.
Suppose the unique minimizer ofEu,v is u∗,v∗. Then, u∗ and v∗ are the unique fixed points of T2 and T1, respectively.Proof.
SinceEu,v is differentiable with respect to u and v separately, we obtain(47)∂E∂uu∗,v∗∂E∂vu∗,v∗=00.
This implies that(48)u∗=ℛv∗=argminuEu,v∗,v∗=Su∗=argminvEu∗,v.
We easily getu∗=ℛv∗=T2u∗ and v∗=Su∗=T1v∗. Therefore, u∗ and v∗ are the corresponding fixed points of T1 and T2.On the other hand, we note thatEu,v is strictly convex and differentiable about u and v, respectively. Therefore, the fixed points of T1 and T2 are the minimizers of Eu,v. By virtue of the uniqueness of the minimizer of Eu,v, T1 and T2 show a unique fixed point separately. That is to say, that u∗ and v∗ are the unique fixed points of T2 and T1, respectively.According to Theorem 1 in [7], we get the following result.Theorem 2.
The sequenceun,vn generated by Algorithm 1 converges to the solution u∗,v∗ of problems (19) and (20).Proof.
We should first notice the fact that(49)vn=minvFun−1,vun=minuFu,vn,so we can get(50)Fun,vn≤Fun−1,vn≤Fun−1,vn−1.
In particular, the sequenceFun,vn is nonincreasing.Theorem 3.
For any initial valuesv0andu0, assume uk and vk are generated by (2) and (3), respectively; then ukandvk converge to the corresponding fixed points of T2 and T1, i.e., uk,vk converges to u∗,v∗ which is the unique minimizer of Eu,v, as k⟶∞.Proof.
In view of Lemmas1–5, we know that T1:Rn⟶RnandT2:Rn⟶Rn are nonexpansive asymptotically regular mappings and have fixed points. By using Theorem 1 in [7], we get that ukandvk converge, respectively, to the fixed points of T2 and T1, i.e., uk,vk converges to u∗,v∗ which is the unique minimizer of Eu,v, as k⟶∞.


## 2.2.1. Energy Minimization with respect toc1 and c2

We fixedu and v in order to determine an exact solution of c1 and c2; then, the minimization of Ec1,c2,u,v with respect to c1 and c2 is as follows:(16)argminc1λ∫ΩI0−c12vdx,argminc2λ∫Ω−I0−c22vdx.Then, we can obtain(17)c1=c2=∫I0∫I.


## 2.2.2. Energy Minimization with respect tou and v

For any givenc1 and c2, we have the problem(18)argminu,v≔∫ΩG∇u2,1dx+12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx.Problem (18) is convex and is not strictly convex about u and v, so its local minimizer is also the global minimizer. One of the effective methods is to change it into an unconstrained minimization problem according to the theorems in [14, 24].Obviously, it is easy to find that the objective function in problem (18) is strictly convex and coercive, so it owns a unique solution u∗,v∗. Formally, problem (18) including two variables u and v as a coupled problem, and the classical method is to divide it into two subproblems. Thus, we alteratively consider the following minimization problem:(19)Sv≔u=minuG∇u22+12θu−v22,(20)Tu≔v=minv12θu−v22+λ∫ΩI0−c12−I0−c22v+αvvdx.Theorem 1.
The following assertions hold:(1)
The solution of problem (19) is given by(21)u=v−PθSv,whereP is the project operator and S is defined by(22)S=closuredivξx:ξx∈Cc1Ω,R2,ξx≤gx,∀x∈Ω.(2)
The solution of problem (17) satisfies(23)v=minmaxu−θλr1x,c1,c2,0,1,wherer1x,c1,c2=I0−c12−I0−c22.Remark 1.
Obviously, problem (19) is a generated ROF model by adding a weighted function gx to the regularization term and so we can get the solution based on the scheme in [27,31,32]. For problem (20), it is similar to what happens in [28] for this circumstance; we notice that it is equivalent to the following minimization problem:(24)minv12θu−v22+λ∫Ωr1x,c1,c2u+δCvdx,where δ is the indicator function defined by(25)δx=0,ifx∈C,+∞,ifx∉C,with C=0,1. Thus, we can easily find that the solution of problem (20) satisfies condition (23).Definition 1.
The proximal operator of a proper, convex, semicontinuous functionφ on Rn is defined for x∈Rn by(26)proxγφ∘Bx≔argminy∈Rn12γy−x22+φBy,where B:Rn⟶Rn is a bounded linear operator. Furthermore, when B=I, it is the classical proximal operator defined in [33].Definition 2.
Letφ be a real-value convex function on Rn. The subdifferential of φ at x∈Rn is defined by(27)∂φx≔y:y∈Rnandφz≥φx+y,z−x∀z∈Rn.
The elements in∂φx are called subgradients. Especially, ∂φx=∇φx when φx is differential.Definition 3.
An operatorQ:Rn⟶Rn is firmly nonexpansive if it satisfies one of the following equivalent conditions:(i)
Qx−Qy22≤Qx−Qy,x−y,∀x,y∈Rn×Rn(ii)
Qx−Qy22≤x−y22−I−Qx−I−Qy22,∀x,y∈Rn×RnDefinition 4.
An operatorQ:Rn⟶Rn is nonexpansive if for any x,y∈Rn×Rn, we have(28)Qx−Qy2≤x−y2.It follows that a firmly nonexpansive operatorP is nonexpansive. Based on the proof of the proximal operator in [34], we obtain the similar result of the proximal operator defined in Definition 1.Lemma 1.
The proximal operatorproxγφ∘B defined in Definition 1 is firmly nonexpansive.Proof.
Following the definition of the proximity in Definition1 and the chain rule, we can get(29)φB∘proxγφ∘By−φB∘proxγφ∘Bx≥1γx−proxγφ∘Bx,proxγφ∘By−proxγφ∘Bx,φB∘proxγφ∘Bx−φB∘proxγφ∘By≥1γy−proxγφ∘By,proxγφ∘Bx−proxγφ∘By.We obtain(30)proxγφ∘Bx−proxγφ∘By22≤proxγφ∘Bx−proxγφ∘By,x−y,which implies the nonexpansion of the operator proxγφ∘B.Corollary 1.
The operatorS defined in minimization problem (19) is nonexpansive.Proof.
Based on above definitions, setB=gx∇u; then, the operator T is the proximal operator of ∫Ωgx∇udx. Furthermore, we can get S=pro⋅L1∘B. According to Lemma 1, the operator S is nonexpansive.Lemma 2.
Assume thatC⊂Rn is a closed convex set and PC⋅:Rn⟶C is a projection operator; then, the projector PC⋅ is firmly nonexpansive.Proof.
Letx∈Rn and ∀y∈C; we have(31)y−PCx,x−PCx≤0.
SincePCy=y while y∈C, we can get(32)PCy−PCx,x−y+PCy−PCx≤0.
This yields(33)PCx−PCy22≤x−y,PCx−PCy≤0.Corollary 2.
The operatorT defined in minimization problem (20) is firmly nonexpansive.Proof.
Obviously, from Theorem1 and Remark 1, we derive that(34)Tu≔v=P0,1u−θλr1x,c1,c2.
It follows from Lemma2 that T is firmly nonexpansive.In view of the above results, we know thatℛ1⋅=ST⋅ and ℛ2⋅=TS⋅ are firmly nonexpansive.Lemma 3.
Letuk and vk be generated by (2) and (3), respectively. Then, ∑k=0∞uk−uk+1L22 and ∑k=0∞vk−vk+1L22 are convergent.Proof.
DenoteE1u,v=u−vL22 and E2v=∇xvL1+∇yvL1. We have(35)Eu,v=∫Ωu+ge−udxdy+β1E1u,v+β2E2v.
Hence,(36)Euk+1,vk−Euk+1,vk+1=β1E1uk+1,vk−E1uk+1,vk+1+β2E2vk−E2vk+1.
Consider the Taylor series expansion ofE1u,v in the second variable, i.e.,(37)E1uk+1,vk=E1uk+1,vk+1+vk−vk+1T∂E1∂vuk+1,vk+1+12vk−vk+1T∂2E1∂v2uk+1,vk+1vk−vk+1.
Here,xT denotes the transpose of x. We notice that E1 is quadratic about v. Then, ∂2E1/∂v2≡2I where I is the identity matrix. Moreover, since E2 is a convex function, we get(38)E2vk≥E2vk+1+vk−vk+1T∂E2∂vvk+1.
Combining (36), (37), and (38), we obtain(39)Euk+1,vk−Euk+1,vk+1≥uk−uk+1Tβ1∂E1∂vuk+1,vk+1+β2∂E2∂vvk+1+β1vk−vk+1L22.
The subdifferential ofE with respect to v is equivalent to the vector sum of the subdifferential of E1 and E2 about v, i.e.,(40)∂E∂v=β1∂E1∂v+β2∂E2∂v.
Sincevk+1 is the minimizer of Euk+1,v, we have(41)∂E∂vuk+1,vk+1=0,that is(42)β1∂E1∂vuk+1,vk+1+β2∂E2∂vzk+1,vk+1=0.
Therefore, the first term in the right hand side of (39) is zero. When we solve successive minimization problems (19) and (20), we note that Euk,vk≥Euk+1,vk. Hence, we get(43)Euk,vk−Euk+1,vk+1≥Euk+1,vk−Euk+1,vk+1≥β1vk−vk+1L22.
It follows that the partial sum of the sequencevk−vk+1L22k=0∞ is bounded. Thus, the infinite series ∑k=0∞vk−vk+1L22 is convergent.
LetE3u=∫Ωu+ge−udxdy. By considering Euk,vk−Euk+1,vk and using the similar method, we can prove that ∑k=0∞uk−uk+1L22 is convergent.Definition 5.
A operatorP:Rn⟶Rn is asymptotically regular if for any x in Rn, the sequence Pk+1x−Pkx tends to zero as k⟶∞.Based on Lemma3, we get the following result.Lemma 4.
For any initial valuesv0andu0, assume uk and vk are generated by (19) and (20), respectively; then, T1 and T2 are asymptotically regular.Proof.
According to Lemma3, we obtain(44)limk⟶∞vk−vk+1L22=0,limk⟶∞uk−uk+1L22=0.
Asvk=T1vk−1 and uk=T2uk−1, by using the recurrence method, we have(45)vk=T1kv0,uk=T2ku0.
Therefore, we get(46)limk⟶∞T1k+1v0−T1kv0L2=0,limk⟶∞T2k+1u0−T2ku0L2=0.
This indicates thatT1 and T2 are asymptotically regular.Lemma 5.
Suppose the unique minimizer ofEu,v is u∗,v∗. Then, u∗ and v∗ are the unique fixed points of T2 and T1, respectively.Proof.
SinceEu,v is differentiable with respect to u and v separately, we obtain(47)∂E∂uu∗,v∗∂E∂vu∗,v∗=00.
This implies that(48)u∗=ℛv∗=argminuEu,v∗,v∗=Su∗=argminvEu∗,v.
We easily getu∗=ℛv∗=T2u∗ and v∗=Su∗=T1v∗. Therefore, u∗ and v∗ are the corresponding fixed points of T1 and T2.On the other hand, we note thatEu,v is strictly convex and differentiable about u and v, respectively. Therefore, the fixed points of T1 and T2 are the minimizers of Eu,v. By virtue of the uniqueness of the minimizer of Eu,v, T1 and T2 show a unique fixed point separately. That is to say, that u∗ and v∗ are the unique fixed points of T2 and T1, respectively.According to Theorem 1 in [7], we get the following result.Theorem 2.
The sequenceun,vn generated by Algorithm 1 converges to the solution u∗,v∗ of problems (19) and (20).Proof.
We should first notice the fact that(49)vn=minvFun−1,vun=minuFu,vn,so we can get(50)Fun,vn≤Fun−1,vn≤Fun−1,vn−1.
In particular, the sequenceFun,vn is nonincreasing.Theorem 3.
For any initial valuesv0andu0, assume uk and vk are generated by (2) and (3), respectively; then ukandvk converge to the corresponding fixed points of T2 and T1, i.e., uk,vk converges to u∗,v∗ which is the unique minimizer of Eu,v, as k⟶∞.Proof.
In view of Lemmas1–5, we know that T1:Rn⟶RnandT2:Rn⟶Rn are nonexpansive asymptotically regular mappings and have fixed points. By using Theorem 1 in [7], we get that ukandvk converge, respectively, to the fixed points of T2 and T1, i.e., uk,vk converges to u∗,v∗ which is the unique minimizer of Eu,v, as k⟶∞.


## 3. General Framework of the Proposed Model

Now, we discuss the solution to each variable separately; without loss of generality, we use the ADMM method to solve minimization problem (18).(51)u≔minu∫ΩG∇u2,1dx+12θu−v22,(52)v≔minv12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx.Then, in the following sections, we obtain the solutions to problems (51) and (52).


### 3.1. Algorithm for Solving (51)

In this section, we solve model (51) about u; because the model is not easy to solve, we use ADMM to solve problem (51). We can take w=∇u and T=Gw in our model, and by introducing two variables, we can obtain the following constrained problem:(53)minu,w,T∫ΩT2,1dx+12θu−vL22s.t.w=∇u,T=Gw,where(54)w≔wxwy=uxuy,T≔TxTy=gx00gywxwy.By using the augmented Lagrangian method to turn the constrained problem into the unconstrained problem, we have the augmented Lagrangian function as follows:(55)ℒu,w,T;β1,β2=∫ΩT2,1dx+12θu−vL22+β1,w−∇u+γ12w−∇uL22+β2,T−Gw+γ22T−GwL22,where β1 and β2 are Lagrangian multiplier that can alternatively be regarded as the dual variables of problem w and T and γ1>0 and γ2>0 are penalty parameters which are constant that should be chosen properly. We convert it into the maximum minimization problem and simplify(56)minu,w,Tmaxβ1,β2∫ΩT2,1dx+12θu−vL22+γ12w−∇u+β1γ1L22+γ22T−Gw+β2γ2L22.In the following, we use ADMM to solve each subproblem about (56). First, the minimum energy function Eu,w,T is as follows:(57)minu,w,TEu,w,T.The maximum energy functionEβ1,β2 is as follows:(58)maxβ1,β2Eβ1,β2.Then, we solve minimization and maximization problems by using ADMM.(1) Subminimization with respect to. T. We fixed uandw in order to determine an exact solution of T; then, we have the minimization of T with respect to uandw as follows:(59)limT∫ΩT21dx+γ22T−Gw+β2γ2L22.In order to solve T, we use the soft threshing operator defined in [35] and obtain(60)T=maxGw−β2γ21−1γ2,0Gw−β2/γ2Gw−β2/γ21.For the corresponding Lagrangian multiplierβ2, we update it by using the gradient ascent method:(61)β2=β2+γ2T−Gw.(2) Subminimization with respect tow. We fixed uandT in order to determine an exact solution of w; the minimization of T with respect to uandw is as follows:(62)minwγ12w−∇u+β1γ1L22+γ22T−Gw+β2γ2L22,and then we calculate the equation about the w and obtain(63)γ1w−∇u+β1γ1−γ2GTT−Gw+β2γ2=0.Simplifying equation (63), we can obtain(64)w=γ2GT+β2/γ2+γ1∇u−β1/γ1γ1+γ2G2.For the corresponding Lagrangian multiplierβ1, we update it by using the gradient ascent method:(65)β1=β1+γ1w−∇u.(3) Subminimization with respect tou. We fixed Tandw in order to determine an exact solution of u; we begin by considering the minimization of (55) with respect to u, and then,(66)minu12θu−vL22+γ12w−∇u+β1γ1L22.Then, using the divergence theorem in the equation about u, we obtain(67)1θu−v+γ1divw−∇u+β1γ1=0.Similarly, due to periodic boundary condition being imposed, we can solveu efficiently by the fast Fourier transform (FFT). This yields the following solution:(68)u=ℱ−1ℱv−θγ1ℱdivw+β1/γ11−θγ1ℱΔ.with ℱ being the FFT operator and ℱ−1 being the complex conjugate. In (68), ℱdivw+β1/γ1=ℱ∇1+∘ℱwx+β1x/γ1+ℱ∇2+∘ℱwy+β1y/γ1, where “∘” denotes the componentwise multiplication and the division is also componentwise.


### 3.2. Algorithm for Solving (52)

In this section, we solve model (52) about v:(69)v≔minv12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx.The details are as follows: as introduced in (23) and since the constrained variable v=u, we calculate the minimization problem and get(70)θλI0−c12−I0−c22+v−u=0,that is,(71)vx=ux−θλI0−c12−I0−c22.Thus, the v-minimization can be achieved through the following update:(72)vx=minmaxux−θλI0−c12−I0−c22,0,1.Then, the framework of the algorithm which solves problem (52) can be illustrated as follows.In this section, we give experimental justification of our proposed model and compare it with ACWE model [8] and Bresson method [24]. These experiments show that our method is robust in segmentation of objects in a range of images that have noise, have low contrast, and are directional. The proposed model was implemented by Matlab7 on a computer with Intel Core2 Duo2.2 GHz CPU, 2GB RAM, and Windows XP operating system. In numerical implementations, we only consider to use proposed model (52) for the basic image segmentation problems. In fact, our model has been applied to synthetic and real images in this section.


### 3.3. Experimental Results of Real Images

In this section, the effectiveness of the segmentation models and algorithms is verified by some synthetic images. We first segment every original image into four images by adding 0.1, 0.15, and 0.2 random noise. The computer is Windows 7, a 64-bit operating system. The four original images with their three-level noise have been processed, and the results are illustrated in Figures1–4. Slight tuning of model parameters is necessary between images (but not between models for a single image).Figure 1
Illustration of without noisy image and three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure 2
Illustration of without noisy natural image and three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure 3
Illustration of without noisy natural image as well as three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure 4
Illustration of without noisy natural image as well as three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure5 gives four original images which will be segmented. Figures 1–4 show the segmentation results of the proposed method and the related methods. The first line of Figure 1 shows the segmentation of the original pictures, and the second and third lines are the segmentation results by adding 0.1 noise and 0.15 noise to the original image, respectively; we can see that our method can be used to successfully segment for the bottom part of the image, and the other two methods all failed; from the fourth line, although we have not fully segmented, we can see that our segmentation is more full at several locations below the image.Figure 5
The given four original images.Figure2 shows that our method is successful when the image has no noise and 0.1 random noise. When the image is added 0.15 random noise and 0.2 random noise, the CV method failed in the three methods, and compared with the Bresson method, we can see l and p provide more details in the lion’s neck. In Figure 3, compared with the CV method, we have a more obvious advantage for the arc of the aircraft’s tail and the head of the aircraft. And compared with the Bresson method, for the letter A on the tail, we are more clear. From Figure 4, we can see that our method is better.


### 3.4. Experimental Results of Synthetic Images

We give three synthetic images and add three different degrees of random noise in this part, and the effectiveness of the segmentation models and algorithms proposed in the following is verified by comparison with related models. We test four typical real images and tested each image with three different degrees of random noise. Figure6 shows that the size of the synthetic images is 128×128 and the pixel values of these images are only 0.4 and 0.8.Figure 6
Sizes of the synthetic images is128×128; the pixel values of these images are only 0.4 and 0.8.
(a)(b)(c)Before doing these experiments, we make some explanations as follows:(1)
For the above three synthetic images, we have two pixel values of 0.4 and 0.8.(2)
We use the index error rate (err) to measure the similarity of the initial image and the segmentation, which are defined as follows, whereP represent the test image, I represent the initial image:(73)errP,I=P≠II.From the definition, we can obtain that the closer the value of err is to 0, the closer the imageP is to the image I (Figures 7–9).Figure 7
Segmentation results to the synthetic image by using Bresson method [24], Chan–Vese [8], and our Algorithm 2 for different levels of noises. First row: the result of three methods for the original image segmentation with 0.1 noise. Second row: the result of three methods for the original image segmentation with 0.15 noise. Third row: the result of three methods for the original image segmentation with 0.2 noise.Figure 8
Segmentation results to the synthetic image by using Bresson method [24], Chan–Vese [8], and our Algorithm 2 for different levels of noises. First row: the result of three methods for the original image segmentation with 0.1 noise. Second row: the result of three methods for the original image segmentation with 0.15 noise. Third row: the result of three methods for the original image segmentation with 0.2 noise.Figure 9
Segmentation results to the synthetic image by using Bresson method [24], Chan–Vese [8], and our Algorithm 2 for different levels of noises. First row: the result of three methods for the original image segmentation with 0.1 noise. Second row: the result of three methods for the original image segmentation with 0.15 noise. Third row: the result of three methods for the original image segmentation with 0.2 noise.Algorithm 2: Solving the model (52).
Inputθ,γ1,γ1>0, set the starting values w=T, γ1=γ2, the constants c1 and c2 are updated periodically every 10 iterations, let k=0 and start k−th iteration, which includes the following steps, until a stopping criterion maxun+1−un,vn+1−vn≤ε is satisfied:(i)
UpdateT by using Tk+1=maxGwk−β2k/γ21−1/γ2,0Gwk−β2k/γ2/Gwk−β2k/γ21,(ii)
Updateβ2 by using β2k+1=β2k+γ2Tk+1−Gwk,(iii)
Update thew, solving wk by using wk+1=γ2GTk+1+β2k+1/γ2+γ1uk−β1k/γ1/γ1+γ2G2,(iv)
Updateβ1 by using β1k+1=β1k+γ1wk−∇uk,(v)
Updateu by using u=ℱ−1ℱvk−θγ1ℱdivwk+1+β1k+1/γ1/1−θγ1ℱΔ,(vi)
Updatev by using v=minmaxuk+1−θλI0−c12−I0−c22,0,1,(vii)
Letk=k+1 go to the k+1 iteration until converge.The error rate for Figure7 by using Algorithm 2 is shown in Table 1.Table 1
The err value in Figure7 by using the Algorithm 2.
Image with 0.1 noiseImage with 0.15 noiseImage with 0.2 noiseChan−Vese method0.00850.01500.0356Bresson method8.5 449e − 040.00310.0143Our method7.3 242e − 040.00290.0120Bold indicates that the err rate is the smallest compared with other two methods.The error rate for the synthetic image is shown in Table2.Table 2
The err value in Figure8 through Algorithm 2.
Image with 0.1 noiseImage with 0.15 noiseImage with 0.2 noiseChan−Vese method0.01390.03310.0598Bresson method0.00310.01140.0199Our method0.00140.00930.0190The error rate for the synthetic image is shown in Table3.Table 3
The err value in Figure9 by Algorithm 2.
Image with 0.1 noiseImage with 0.15 noiseImage with 0.2 noiseChan−Vese method0.00900.02440.0422Bresson method0.00130.00450.0114Our method0.00100.00400.0109The bold values illustrate that our method is the most effective since it is smallest.The closer the value of err is to 0, the better the quality of the segmentation is. The err tests can be seen in Tables1 to 3 where all the three models are tested on simple synthetic images, and these results validate quantitatively the outstanding performance of the Chan–Vese model and the Bresson method in comparison with the competing models.


## 3.1. Algorithm for Solving (51)

In this section, we solve model (51) about u; because the model is not easy to solve, we use ADMM to solve problem (51). We can take w=∇u and T=Gw in our model, and by introducing two variables, we can obtain the following constrained problem:(53)minu,w,T∫ΩT2,1dx+12θu−vL22s.t.w=∇u,T=Gw,where(54)w≔wxwy=uxuy,T≔TxTy=gx00gywxwy.By using the augmented Lagrangian method to turn the constrained problem into the unconstrained problem, we have the augmented Lagrangian function as follows:(55)ℒu,w,T;β1,β2=∫ΩT2,1dx+12θu−vL22+β1,w−∇u+γ12w−∇uL22+β2,T−Gw+γ22T−GwL22,where β1 and β2 are Lagrangian multiplier that can alternatively be regarded as the dual variables of problem w and T and γ1>0 and γ2>0 are penalty parameters which are constant that should be chosen properly. We convert it into the maximum minimization problem and simplify(56)minu,w,Tmaxβ1,β2∫ΩT2,1dx+12θu−vL22+γ12w−∇u+β1γ1L22+γ22T−Gw+β2γ2L22.In the following, we use ADMM to solve each subproblem about (56). First, the minimum energy function Eu,w,T is as follows:(57)minu,w,TEu,w,T.The maximum energy functionEβ1,β2 is as follows:(58)maxβ1,β2Eβ1,β2.Then, we solve minimization and maximization problems by using ADMM.(1) Subminimization with respect to. T. We fixed uandw in order to determine an exact solution of T; then, we have the minimization of T with respect to uandw as follows:(59)limT∫ΩT21dx+γ22T−Gw+β2γ2L22.In order to solve T, we use the soft threshing operator defined in [35] and obtain(60)T=maxGw−β2γ21−1γ2,0Gw−β2/γ2Gw−β2/γ21.For the corresponding Lagrangian multiplierβ2, we update it by using the gradient ascent method:(61)β2=β2+γ2T−Gw.(2) Subminimization with respect tow. We fixed uandT in order to determine an exact solution of w; the minimization of T with respect to uandw is as follows:(62)minwγ12w−∇u+β1γ1L22+γ22T−Gw+β2γ2L22,and then we calculate the equation about the w and obtain(63)γ1w−∇u+β1γ1−γ2GTT−Gw+β2γ2=0.Simplifying equation (63), we can obtain(64)w=γ2GT+β2/γ2+γ1∇u−β1/γ1γ1+γ2G2.For the corresponding Lagrangian multiplierβ1, we update it by using the gradient ascent method:(65)β1=β1+γ1w−∇u.(3) Subminimization with respect tou. We fixed Tandw in order to determine an exact solution of u; we begin by considering the minimization of (55) with respect to u, and then,(66)minu12θu−vL22+γ12w−∇u+β1γ1L22.Then, using the divergence theorem in the equation about u, we obtain(67)1θu−v+γ1divw−∇u+β1γ1=0.Similarly, due to periodic boundary condition being imposed, we can solveu efficiently by the fast Fourier transform (FFT). This yields the following solution:(68)u=ℱ−1ℱv−θγ1ℱdivw+β1/γ11−θγ1ℱΔ.with ℱ being the FFT operator and ℱ−1 being the complex conjugate. In (68), ℱdivw+β1/γ1=ℱ∇1+∘ℱwx+β1x/γ1+ℱ∇2+∘ℱwy+β1y/γ1, where “∘” denotes the componentwise multiplication and the division is also componentwise.


## 3.2. Algorithm for Solving (52)

In this section, we solve model (52) about v:(69)v≔minv12θu−v22+λ∫ΩI0−c12−I0−c22v+ανvdx.The details are as follows: as introduced in (23) and since the constrained variable v=u, we calculate the minimization problem and get(70)θλI0−c12−I0−c22+v−u=0,that is,(71)vx=ux−θλI0−c12−I0−c22.Thus, the v-minimization can be achieved through the following update:(72)vx=minmaxux−θλI0−c12−I0−c22,0,1.Then, the framework of the algorithm which solves problem (52) can be illustrated as follows.In this section, we give experimental justification of our proposed model and compare it with ACWE model [8] and Bresson method [24]. These experiments show that our method is robust in segmentation of objects in a range of images that have noise, have low contrast, and are directional. The proposed model was implemented by Matlab7 on a computer with Intel Core2 Duo2.2 GHz CPU, 2GB RAM, and Windows XP operating system. In numerical implementations, we only consider to use proposed model (52) for the basic image segmentation problems. In fact, our model has been applied to synthetic and real images in this section.


## 3.3. Experimental Results of Real Images

In this section, the effectiveness of the segmentation models and algorithms is verified by some synthetic images. We first segment every original image into four images by adding 0.1, 0.15, and 0.2 random noise. The computer is Windows 7, a 64-bit operating system. The four original images with their three-level noise have been processed, and the results are illustrated in Figures1–4. Slight tuning of model parameters is necessary between images (but not between models for a single image).Figure 1
Illustration of without noisy image and three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure 2
Illustration of without noisy natural image and three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure 3
Illustration of without noisy natural image as well as three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure 4
Illustration of without noisy natural image as well as three noisy image and the performance of all the three models by 3000 iterations. First row: the result of three methods for the original image segmentation without noise. Second row: the result of three methods for the original image segmentation with 0.1 noise. Third row: the result of three methods for the original image segmentation with 0.15 noise. Fourth row: the result of three methods for the original image segmentation with 0.2 noise.Figure5 gives four original images which will be segmented. Figures 1–4 show the segmentation results of the proposed method and the related methods. The first line of Figure 1 shows the segmentation of the original pictures, and the second and third lines are the segmentation results by adding 0.1 noise and 0.15 noise to the original image, respectively; we can see that our method can be used to successfully segment for the bottom part of the image, and the other two methods all failed; from the fourth line, although we have not fully segmented, we can see that our segmentation is more full at several locations below the image.Figure 5
The given four original images.Figure2 shows that our method is successful when the image has no noise and 0.1 random noise. When the image is added 0.15 random noise and 0.2 random noise, the CV method failed in the three methods, and compared with the Bresson method, we can see l and p provide more details in the lion’s neck. In Figure 3, compared with the CV method, we have a more obvious advantage for the arc of the aircraft’s tail and the head of the aircraft. And compared with the Bresson method, for the letter A on the tail, we are more clear. From Figure 4, we can see that our method is better.


## 3.4. Experimental Results of Synthetic Images

We give three synthetic images and add three different degrees of random noise in this part, and the effectiveness of the segmentation models and algorithms proposed in the following is verified by comparison with related models. We test four typical real images and tested each image with three different degrees of random noise. Figure6 shows that the size of the synthetic images is 128×128 and the pixel values of these images are only 0.4 and 0.8.Figure 6
Sizes of the synthetic images is128×128; the pixel values of these images are only 0.4 and 0.8.
(a)(b)(c)Before doing these experiments, we make some explanations as follows:(1)
For the above three synthetic images, we have two pixel values of 0.4 and 0.8.(2)
We use the index error rate (err) to measure the similarity of the initial image and the segmentation, which are defined as follows, whereP represent the test image, I represent the initial image:(73)errP,I=P≠II.From the definition, we can obtain that the closer the value of err is to 0, the closer the imageP is to the image I (Figures 7–9).Figure 7
Segmentation results to the synthetic image by using Bresson method [24], Chan–Vese [8], and our Algorithm 2 for different levels of noises. First row: the result of three methods for the original image segmentation with 0.1 noise. Second row: the result of three methods for the original image segmentation with 0.15 noise. Third row: the result of three methods for the original image segmentation with 0.2 noise.Figure 8
Segmentation results to the synthetic image by using Bresson method [24], Chan–Vese [8], and our Algorithm 2 for different levels of noises. First row: the result of three methods for the original image segmentation with 0.1 noise. Second row: the result of three methods for the original image segmentation with 0.15 noise. Third row: the result of three methods for the original image segmentation with 0.2 noise.Figure 9
Segmentation results to the synthetic image by using Bresson method [24], Chan–Vese [8], and our Algorithm 2 for different levels of noises. First row: the result of three methods for the original image segmentation with 0.1 noise. Second row: the result of three methods for the original image segmentation with 0.15 noise. Third row: the result of three methods for the original image segmentation with 0.2 noise.Algorithm 2: Solving the model (52).
Inputθ,γ1,γ1>0, set the starting values w=T, γ1=γ2, the constants c1 and c2 are updated periodically every 10 iterations, let k=0 and start k−th iteration, which includes the following steps, until a stopping criterion maxun+1−un,vn+1−vn≤ε is satisfied:(i)
UpdateT by using Tk+1=maxGwk−β2k/γ21−1/γ2,0Gwk−β2k/γ2/Gwk−β2k/γ21,(ii)
Updateβ2 by using β2k+1=β2k+γ2Tk+1−Gwk,(iii)
Update thew, solving wk by using wk+1=γ2GTk+1+β2k+1/γ2+γ1uk−β1k/γ1/γ1+γ2G2,(iv)
Updateβ1 by using β1k+1=β1k+γ1wk−∇uk,(v)
Updateu by using u=ℱ−1ℱvk−θγ1ℱdivwk+1+β1k+1/γ1/1−θγ1ℱΔ,(vi)
Updatev by using v=minmaxuk+1−θλI0−c12−I0−c22,0,1,(vii)
Letk=k+1 go to the k+1 iteration until converge.The error rate for Figure7 by using Algorithm 2 is shown in Table 1.Table 1
The err value in Figure7 by using the Algorithm 2.
Image with 0.1 noiseImage with 0.15 noiseImage with 0.2 noiseChan−Vese method0.00850.01500.0356Bresson method8.5 449e − 040.00310.0143Our method7.3 242e − 040.00290.0120Bold indicates that the err rate is the smallest compared with other two methods.The error rate for the synthetic image is shown in Table2.Table 2
The err value in Figure8 through Algorithm 2.
Image with 0.1 noiseImage with 0.15 noiseImage with 0.2 noiseChan−Vese method0.01390.03310.0598Bresson method0.00310.01140.0199Our method0.00140.00930.0190The error rate for the synthetic image is shown in Table3.Table 3
The err value in Figure9 by Algorithm 2.
Image with 0.1 noiseImage with 0.15 noiseImage with 0.2 noiseChan−Vese method0.00900.02440.0422Bresson method0.00130.00450.0114Our method0.00100.00400.0109The bold values illustrate that our method is the most effective since it is smallest.The closer the value of err is to 0, the better the quality of the segmentation is. The err tests can be seen in Tables1 to 3 where all the three models are tested on simple synthetic images, and these results validate quantitatively the outstanding performance of the Chan–Vese model and the Bresson method in comparison with the competing models.


## 4. Discussion

In this paper, we have proposed a new variational model suitable for segmentation a range of images with a blurred edge, a certain degree of noise, and a directional texture.L2,1 norm is made for the boundary indication function, and different weights are taken for the x direction and y direction of the boundary at a certain place of the image. Experimental results on both real and synthetic images demonstrated that our method is very robust and efficient. We will further improve the proposed model to extract the colorful image in our future work.



---

*Source: 2852399-2021-12-17.xml*