# Implicit Active Contour Model with Local and Global Intensity Fitting Energies

Authors: Xiaozeng Xu; Chuanjiang He
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We propose a new active contour model which integrates a local intensity fitting (LIF) energy with an auxiliary global intensity fitting (GIF) energy. The LIF energy is responsible for attracting the contour toward object boundaries and is dominant near object boundaries, while the GIF energy incorporates global image information to improve the robustness to initialization of the contours. The proposed model not only can provide desirable segmentation results in the presence of intensity inhomogeneity but also allows for more flexible initialization of the contour compared to the RSF and LIF models, and we give a theoretical proof to compute a unique steady state regardless of the initialization; that is, the convergence of the zero-level line is irrespective of the initial function. This means that we can obtain the same zero-level line in the steady state, if we choose the initial function as a bounded function. In particular, our proposed model has the capability of detecting multiple objects or objects with interior holes or blurred edges.


---

## Body

## 1. Introduction

Implicit active contour models have been extensively studied and successfully used in image segmentation [1–3]. The basic idea is to evolve a contour under some constraints to extract the desired object. According to the nature of constraints, the existing active contour models can be categorised into two classes: edge-based models [4–7], and region-based models [8–17]. Each of them has its own pros and cons; the choice of them in applications depends on different characteristics of images. In this study, we focus on region-based models and consider images with intensity inhomogeneity.Unlike edge-based models that utilize typically an edge indicator depending on image gradient to perform contour extraction, region-based models usually use global and/or local statistics inside regions rather than gradient on edges to find a partition of image domain. They generally have better performances in the presence of weak or discontinuous boundaries than edge-based models. Early popular region-based models tend to rely on intensity homogeneous (roughly constant or smooth) statistics in each of the regions to be segmented. Therefore, they either lack the ability to deal with intensity inhomogeneity like the PC (piecewise constant) model [8] or are excessively expensive in computation like the PS (piecewise smooth) model [9].To handle intensity inhomogeneity efficiently, some localized region-based models [11–16] have been proposed recently. For example, Li et al. [12] recently proposed a region-scalable fitting (RSF) active contour model (originally termed as local binary fitting (LBF) model [11]). The RSF model draws upon intensity information in spatially varying local regions depending on a scale parameter, so it is able to deal with intensity inhomogeneity efficiently. Very recently, Zhang et al. [15] proposed a novel active contour model driven by local image fitting energy, which also can handle intensity inhomogeneity efficiently. The experimental results in [15] show that this model is more efficient than the RSF model, while yielding similar results. However, these two models easily get stuck in local minimums for the most of contour initializations. This makes it need user intervention to define the initial contours professionally.In this study, based on the PC model [8] and RSF model [12], we propose a new active contour model, which integrates a local intensity fitting (LIF) energy with an auxiliary global intensity fitting (GIF) energy. The LIF energy is responsible for attracting the contour toward object boundaries and is dominant near object boundaries, while the GIF energy incorporates global image information to improve the robustness to initialization of the contours. The proposed model can efficiently handle intensity inhomogeneity, while allowing for more flexible initialization and maintaining the subpixel accuracy.The remainder of this paper is organized as follows. Section2 briefly reviews the PC model [8] and RSF model [12]. Section 3 introduces the proposed model. Section 4 presents experimental results using a set of synthetic and real images. This paper is summarized in Section 5.


## 2. Related Works

### 2.1. Piecewise Constant Model by Chan and Vese

In [8], Chan and Vese (CV) proposed a region-based active contour model that relies on intensity homogeneous (roughly constant) statistics in each of the regions to be segmented. In the CV model, a contour C is represented implicitly by the zero-level set of a Lipschitz function ϕ:Ω→R, which is called a level set function. In what follows, we let the level set function ϕ take positive and negative values inside and outside the contour C, respectively.LetI:Ω⊂R2→R be an input image, and let Hε be the regularized Heaviside function; the energy functional of the CV model is defined as
(1)ECV(c1,c2,ϕ)=λ1∫Ω|I-c1|2Hε(ϕ)dx+λ1∫Ω|I-c2|2(1-Hε(ϕ))dx+ν∫Ω|∇Hε(ϕ)|dx,
where λ1,λ2>0, ν>0 are constants. The regularized version of H(z) is chosen as
(2)Hε(z)=12(1+2πarctan(zε)).c1 and c2 are the global averages of the image intensities in the region {x:ϕ(x)>0} and {x:ϕ(x)<0}, respectively; that is,
(3)c1(ϕ)=∫ΩI(x)Hε(ϕ(x))dx∫ΩHε(ϕ(x))dx,c2(ϕ)=∫ΩI(x)(1-Hε(ϕ(x)))dx∫Ω(1-Hε(ϕ(x)))dx.The solution of the CV model in fact leads to a piecewise constant segmentation of the original imageI(x):
(4)IG(x)=c1Hε(ϕ(x))+c2(1-Hε(ϕ(x))),
where c1 and c2 are the averages of the image intensities in the region {x:ϕ(x)>0} and {x:ϕ(x)<0}, respectively. Such constants may be far away from the original image data, if the intensities outside or inside the contour C={x:ϕ(x)=0} are not homogeneous. As a result, the CV model generally fails to segment images with intensity inhomogeneity.


### 2.2. Region-Scalable Fitting Model

In order to improve the performance of the global CV model [8] and the PS model [9] on images with inhomogeneity, Li et al. [11, 12] recently proposed a novel region-based active contour model. They introduced a kernel function and defined the following energy functional:
(5)ERSF(f1,f2,ϕ)=λ1∫[∫Kσ(x-y)|I(y)-f1(x)|2×Hε(ϕ(y))dy∫Kσ(x-y)|I(y)-f1(x)|2]dx+λ2∫[∫Kσ(x-y)|I(y)-f1(x)|2×(1-Hε(ϕ(y)))dy∫Kσ(x-y)|I(y)-f1(x)|2]dx,+ν∫Ω|∇Hε(ϕ(x))|dx+μ∫Ω12(|∇ϕ(x)|-1)2dx,
where Kσ is a Gaussian kernel with standard deviation σ, and f1(x) and f2(x) are two smooth functions that approximate the local image intensities inside and outside the contour, respectively. They are computed by
(6)f1(x)=∫ΩKσ(x-y)I(y)Hε(ϕ(y))dy∫ΩKσ(x-y)Hε(ϕ(y))dy,f2(x)=∫ΩKσ(x-y)I(y)(1-Hε(ϕ(y)))dy∫ΩKσ(x-y)(1-Hε(ϕ(y)))dy.The solution of the RSF model leads to a piecewise smooth approximation of the original imageI(x):
(7)IL(x)=f1(x)Hε(ϕ(x))+f2(x)(1-Hε(ϕ(x))),
where the smooth functions f1(x) and f2(x) are computed by (6). f1(x) and f2(x) are the averages of local intensities on the two sides of the contour, which are different from the constants c1 and c2 in the CV model, the averages of the image intensities on the two sides of the contour. Therefore, the RSF model can deal with intensity inhomogeneity efficiently. However, it is sensitive to contour initialization (initial locations, sizes, and shapes).


## 2.1. Piecewise Constant Model by Chan and Vese

In [8], Chan and Vese (CV) proposed a region-based active contour model that relies on intensity homogeneous (roughly constant) statistics in each of the regions to be segmented. In the CV model, a contour C is represented implicitly by the zero-level set of a Lipschitz function ϕ:Ω→R, which is called a level set function. In what follows, we let the level set function ϕ take positive and negative values inside and outside the contour C, respectively.LetI:Ω⊂R2→R be an input image, and let Hε be the regularized Heaviside function; the energy functional of the CV model is defined as
(1)ECV(c1,c2,ϕ)=λ1∫Ω|I-c1|2Hε(ϕ)dx+λ1∫Ω|I-c2|2(1-Hε(ϕ))dx+ν∫Ω|∇Hε(ϕ)|dx,
where λ1,λ2>0, ν>0 are constants. The regularized version of H(z) is chosen as
(2)Hε(z)=12(1+2πarctan(zε)).c1 and c2 are the global averages of the image intensities in the region {x:ϕ(x)>0} and {x:ϕ(x)<0}, respectively; that is,
(3)c1(ϕ)=∫ΩI(x)Hε(ϕ(x))dx∫ΩHε(ϕ(x))dx,c2(ϕ)=∫ΩI(x)(1-Hε(ϕ(x)))dx∫Ω(1-Hε(ϕ(x)))dx.The solution of the CV model in fact leads to a piecewise constant segmentation of the original imageI(x):
(4)IG(x)=c1Hε(ϕ(x))+c2(1-Hε(ϕ(x))),
where c1 and c2 are the averages of the image intensities in the region {x:ϕ(x)>0} and {x:ϕ(x)<0}, respectively. Such constants may be far away from the original image data, if the intensities outside or inside the contour C={x:ϕ(x)=0} are not homogeneous. As a result, the CV model generally fails to segment images with intensity inhomogeneity.


## 2.2. Region-Scalable Fitting Model

In order to improve the performance of the global CV model [8] and the PS model [9] on images with inhomogeneity, Li et al. [11, 12] recently proposed a novel region-based active contour model. They introduced a kernel function and defined the following energy functional:
(5)ERSF(f1,f2,ϕ)=λ1∫[∫Kσ(x-y)|I(y)-f1(x)|2×Hε(ϕ(y))dy∫Kσ(x-y)|I(y)-f1(x)|2]dx+λ2∫[∫Kσ(x-y)|I(y)-f1(x)|2×(1-Hε(ϕ(y)))dy∫Kσ(x-y)|I(y)-f1(x)|2]dx,+ν∫Ω|∇Hε(ϕ(x))|dx+μ∫Ω12(|∇ϕ(x)|-1)2dx,
where Kσ is a Gaussian kernel with standard deviation σ, and f1(x) and f2(x) are two smooth functions that approximate the local image intensities inside and outside the contour, respectively. They are computed by
(6)f1(x)=∫ΩKσ(x-y)I(y)Hε(ϕ(y))dy∫ΩKσ(x-y)Hε(ϕ(y))dy,f2(x)=∫ΩKσ(x-y)I(y)(1-Hε(ϕ(y)))dy∫ΩKσ(x-y)(1-Hε(ϕ(y)))dy.The solution of the RSF model leads to a piecewise smooth approximation of the original imageI(x):
(7)IL(x)=f1(x)Hε(ϕ(x))+f2(x)(1-Hε(ϕ(x))),
where the smooth functions f1(x) and f2(x) are computed by (6). f1(x) and f2(x) are the averages of local intensities on the two sides of the contour, which are different from the constants c1 and c2 in the CV model, the averages of the image intensities on the two sides of the contour. Therefore, the RSF model can deal with intensity inhomogeneity efficiently. However, it is sensitive to contour initialization (initial locations, sizes, and shapes).


## 3. The Proposed Model

### 3.1. Description of Our Model

Given a positive constantω (0≤ω≤1), from (4) and (7), we define the following energy functional:
(8)E(ϕ,c1,c2,f1,f2)=ωEG(ϕ)+(1-ω)EL(ϕ)=ω(12∫Ω|I(x)-IG(x)|2dx)+(1-ω)(12∫Ω|I(x)-IL(x)|2dx),
where
(9)IG(x)=c1Hε(ϕ(x))+c2(1-Hε(ϕ(x))),IL(x)=f1(x)Hε(ϕ(x))+f2(x)(1-Hε(ϕ(x))).Keepingϕ fixed and minimizing the functional E(ϕ,c1,c2,f1,f2) with respect to c1, c2, f1, and f2, we have
(10)c1=∫ΩI(x)Hε(ϕ(x))dx∫ΩHε(ϕ(x))dx,c2=∫ΩI(x)(1-Hε(ϕ(x)))dx∫Ω(1-Hε(ϕ(x)))dx,f1(x)=∫ΩKσ(x-y)I(y)Hε(ϕ(y))dy∫ΩKσ(x-y)Hε(ϕ(y))dy,f2(x)=∫ΩKσ(x-y)I(y)(1-Hε(ϕ(y)))dy∫ΩKσ(x-y)(1-Hε(ϕ(y)))dy.Keepingc1, c2, f1, and f2 fixed and minimizing the functional E(ϕ,c1,c2,f1,f2) with respect to  ϕ, we obtain the corresponding gradient descent flow equation:
(11)∂ϕ∂t=δε(ϕ)[ω((I-IG)(c1-c2))+(1-ω)((I-IL)(f1-f2))]=δε(ϕ)[ωFG(ϕ)+(1-ω)FL(ϕ)]=F(ϕ),
where
(12)δε(z)=Hε′(z)=1πεε2+z2.Like the CV and RSF models, our model is also implemented using an alternative procedure: for each iteration and the corresponding level set functionϕn, we first compute the fitting values  ci(ϕn) and fi(ϕn) and then obtain ϕn+1 by minimizing E(ϕ,c1(ϕn),c2(ϕn),f1(ϕn),f2(ϕn)) with respect to ϕ. This process is repeated until the zero-level set of ϕn+1 is exactly on the object boundary.In the following, we first discuss the properties ofFG(ϕ) and FL(ϕ) and then analyze the behavior of (11).


### 3.2. Properties ofFG(ϕ) and FL(ϕ)

For the sake of simplicity, we state and prove the properties ofFG(ϕ) and FL(ϕ) only for an image consisting of only two distinct gray levels:
(13)I(x)={g1,x∈ω,g2,x∈Ω∖ω,
where g1,g2≥0 with g1≠g2, ω and Ω∖ω represent the objects of interest and the background, respectively.Theorem 1.
LetI(x) be an image by (13). Then one has
(14)FG(ϕ)={(Mn-Nm)(g1-g2)2×[(M-N)(N-n)Hε+N((M+n)-(m+N))(1-Hε)]×(N2(M-N)2)-1,inω-(Mn-Nm)(g1-g2)2×[n(M-N)Hε+N(m-n)(1-Hε)]×(N2(M-N)2)-1,inΩ∖ω,
and so,
(15)sign(FG(ϕ))={+sign(Mn-Nm),inω,-sign(Mn-Nm),inΩ∖ω,
where
(16)M=|Ω|,m=|ω|,N=|{ϕ>0}|,n=|ω∩{ϕ>0}|,
in which |Ω| is the area of the region Ω and similarly for others.Remarks.(i) Due to the discrete nature of image, |Ω| is in fact the number of pixels in the image I(x), and similarly to others.(ii) The cases ofn=m, n=N, and n=0 correspond to the zero-level lines of ϕ(x) which are encircling the object (ω), inside the object and within the background (Ω∖ω), respectively. The cases of N≠n≠m and N=n=m correspond to the zero-level line of ϕ(x) that are partially inside the object and exactly on the object edge, respectively.(iii) The significance of (15) is that the function FG(ϕ) has the opposite sign in ω (object) and Ω∖ω (background), respectively.The proof of Theorem1 is provided in Appendix A. The following result will be used in the proof of Theorem 7, which guarantees that the evolution by (10) converges to a unique stable value after finite time.Corollary 2.
LetI(x) be an image by (13). Then one has the following.(i)
IfMn-Nm≥0, then
(17)FG(ϕ)≥(Mn-Nm)(g1-g2)2M4/4,inω,FG(ϕ)≤-(Mn-Nm)(g1-g2)2M4/4,inΩ∖ω.(ii)
IfMn-Nm≤0, then
(18)FG(ϕ)≤(Mn-Nm)(g1-g2)2M4/4,inω,FG(ϕ)≥-(Mn-Nm)(g1-g2)2M4/4,inΩ∖ω.The proof of Corollary2 is given in Appendix A.We call the property in Theorem1 an adaptive sign-changing property of FG(ϕ). Such property also holds for FL(ϕ), which can be expressed by the following theorem.Theorem 3.
LetI(x) be an image given by (13). Then one has
(19)FL(ϕ)={(Pq-Qp)(g1-g2)2×[(P-Q)(Q-q)Hε(ϕ)+Q((P+q)-(p+Q))(1-Hε(ϕ))]×(Q2(P-Q)2)-1,inω-(Pq-Qp)(g1-g2)2×[q(P-Q)Hε(ϕ)+Q(p-q)(1-Hε(ϕ))]×(Q2(P-Q)2)-1,inΩ∖ω,
and so,
(20)sign(FL(ϕ))={+sign(Pq-Qp),inω,-sign(Pq-Qp),inΩ∖ω,
where
(21)P=∫ΩKσ(x-ξ)dξ,p=∫ωKσ(x-ξ)dξ,Q=∫{ϕ>0}Kσ(x-ξ)dξ,q=∫ω∩{ϕ>0}Kσ(x-ξ)dξ.
This theorem shows that the local forceFL(ϕ) has exactly opposite sign in ω (object) and in Ω∖ω (background).The following result together with Corollary2 will be used in the proof of Theorem 7 mentioned later.Corollary 4.
Under the assumption of Theorem3, one has the following.(i)
IfPq-Qp≥0, then
(22)FL(ϕ)≥S(Pq-Qp)(g1-g2)2P4/4,inω,FL(ϕ)≤-S(Pq-Qp)(g1-g2)2P4/4,inΩ∖ω.(ii)
IfPq-Qp≤0, then
(23)FL(ϕ)≤S(Pq-Qp)(g1-g2)2P4/4,inω,FL(ϕ)≥-S(Pq-Qp)(g1-g2)2P4/4,inΩ∖ω,
where
(24)S=min{Kσ(x-ξ):ξ∈Ω)}.The proofs of Theorem3 and Corollary 4 are similar to Theorem 1 and Corollary 2, respectively; see Appendix B for details.


### 3.3. Behavior of Our Model

In this section, we analyze the behavior of our model (10) for image segmentation. We will show that the zero-level line of the evolution function starting with a bounded function finally comes to a unique steady state, which separates the object from the background.Due to the discrete nature of image, the continuous equation (10) can be discretized in space with space step 1 (implied pixel spacing), but also in time with a time step Δt as follows:
(25)ϕi,jk+1=ϕi,jk+Δt·F(ϕi,jk),0≤i≤p-1,0≤j≤q-1,
where p×q corresponds to the image size, Ii,j=I(i,j) and ϕi,jk=ϕ(kΔt,i,j) with k≥0 and ϕi,j0=ϕ0(i,j).Corollary 5.
LetMn0-N0m=|Ω||ω∩{ϕ0>0}|-|{ϕ0>0}||ω|≠0.(i)
IfMn0-N0m>0, then
(26)FG(ϕk)>(Mn0-N0m)×(g1-g2)2M4/4=A>0,inω(k≥1,k∈z+),FG(ϕk)<-(Mn0-N0m)×(g1-g2)2M4/4=-A<0,inΩ∖ω(k≥1,k∈z+).(ii)
IfMn0-N0m<0, then
(27)FG(ϕk)<(Mn0-N0m)×(g1-g2)2M4/4=A<0,inω(k≥1,k∈z+),FG(ϕk)>-(Mn0-N0m)×(g1-g2)2M4/4=-A>0,inΩ∖ω(k≥1,k∈z+),
where
(28)A=(Mn0-N0m)(g1-g2)2M4/8.We provide the proof of Corollary5 in Appendix C. Similar analysis can prove the following corollary, Corollary 6.Corollary 6.
LetPq0-Q0p≠0.(i)
IfPq0-Q0p>0, then
(29)FL(ϕk)>(Pq0-Q0p)×(g1-g2)2P4/4=B>0,inω(k≥1,k∈z+),FL(ϕk)<-(Pq0-Q0p)×(g1-g2)2P4/4=-B<0,inΩ∖ω(k≥1,k∈z+).(ii)
IfPq0-Q0p<0, then
(30)FL(ϕk)<(Pq0-Q0p)×(g1-g2)2P4/4=B<0,inω(k≥1,k∈z+),FL(ϕk)>-(Pq0-Q0p)×(g1-g2)2P4/4=-B>0,inΩ∖ω(k≥1,k∈z+),
where
(31)B=(Pq0-Q0p)(g1-g2)2P4/8.Theorem 7.
Ifϕ0(x) is a bounded function with (Mn0-N0m)(Pq0-Q0p)>0, then there exists a positive integer K such that, for k≥K,
(32)ϕi,jk+1>0,inω,ϕi,jk+1<0,inΩ∖ω,ifMn0-N0m>0,Pq0-Q0p>0,(33)ϕi,jk+1<0,inω,ϕi,jk+1>0,inΩ∖ω,ifMn0-N0m<0,Pq0-Q0p<0.We provide the proof of Theorem7 in Appendix D.The significance of Theorem7 is that if we choose ϕ0(x) such that (Mn0-N0m)(Pq0-Q0p)>0, then the zero-level line of ϕ(x) starting with such initial function ϕ0(x) can finally come to a unique steady state, which separates object from the background.Remark. Mathematically there does exist ϕ0(x) such that (Mn0-N0m)(Pq0-Q0p)≤0; in practice, however, we always can guarantee that
(34)(Mn0-N0m)(Pq0-Q0p)>0.


### 3.4. Discussion of Initial Function

Theorem7 guarantees that the proposed model computes a unique steady state regardless of the initialization; that is, the convergence of the zero-level line {ϕ=0} is irrespective of the initial function. This means that we can obtain the same zero-level line in the steady state if we choose the initial function as a bounded function ϕ0 with (Mn0-N0m)(Pq0-Q0p)>0.In applications, the initial functionϕ0 can be defined via a simple curve (closed curve or line segment) in image domain. For example, we can choose the initial function as a signed distance to a circle, which is widely used in most of image segmentation models with level set methods. For the proposed model, however, we prefer to define the initial function ϕ0(x) as a piecewise constant or constant function as follows.(1)
If the curveC is a closed curve (e.g., circle or square), then ϕ0(x) is defined by
(35)ϕ0(x)={+ρ,(x)∈inside(C),-ρ,(x)∈outside(C),
where ρ≠0 is a constant.(2)
If the curveC is a line segment that partitions the image domain Ω into two disjoint regions Ω1 and Ω2 (e.g., the left and right half regions of image domain Ω, respectively) with Ω=Ω1∪Ω2 and Ω1∩Ω2=Φ, then ϕ0(x,y) is defined by
(36)ϕ0(x)={+ρ,(x)∈Ω1,-ρ,(x)∈Ω2,
where ρ≠0 is a constant.(3)
We can also define a zero function as follows:(37)ϕ0(x)=0.Next, we prove the fact that, withϕ0(x)=0, ϕ(x) becomes a sign-changing function and satisfy the condition of (Mn1-N1m)(Pq1-Q1p)>0 after the first iteration.Theorem 8.
If the initial functionϕ0(x)=0 in Ω and
(38)c10(x)=f10(x)=2·meanx∈ΩI(x),c20(x)=f20(x)=0,
then after the first iteration, one has
(39)sign(F(0))={-1,x∈Ω1,1,x∈Ω2,
where Ω1∩Ω2=Φ, Ω1≠Φ, Ω2≠Φ.We provide the proof of Theorem8 in Appendix E.By (25) and Theorem 8, we have
(40)ϕ1(x)=δε(0)meanx∈ΩI(x)×(I-meanx∈ΩI(x)){<0,x∈Ω1,>0,x∈Ω2.
Therefore, ϕ1(x) became a sign-changing function. Then, using two distinct gray levels of (13) and the above demonstration, we have
(41)(Mn1-N1m)(Pq1-Q1p)=(m(M-N1))(p(P-Q1))=mp(M-N1)(P-Q1)>0.


## 3.1. Description of Our Model

Given a positive constantω (0≤ω≤1), from (4) and (7), we define the following energy functional:
(8)E(ϕ,c1,c2,f1,f2)=ωEG(ϕ)+(1-ω)EL(ϕ)=ω(12∫Ω|I(x)-IG(x)|2dx)+(1-ω)(12∫Ω|I(x)-IL(x)|2dx),
where
(9)IG(x)=c1Hε(ϕ(x))+c2(1-Hε(ϕ(x))),IL(x)=f1(x)Hε(ϕ(x))+f2(x)(1-Hε(ϕ(x))).Keepingϕ fixed and minimizing the functional E(ϕ,c1,c2,f1,f2) with respect to c1, c2, f1, and f2, we have
(10)c1=∫ΩI(x)Hε(ϕ(x))dx∫ΩHε(ϕ(x))dx,c2=∫ΩI(x)(1-Hε(ϕ(x)))dx∫Ω(1-Hε(ϕ(x)))dx,f1(x)=∫ΩKσ(x-y)I(y)Hε(ϕ(y))dy∫ΩKσ(x-y)Hε(ϕ(y))dy,f2(x)=∫ΩKσ(x-y)I(y)(1-Hε(ϕ(y)))dy∫ΩKσ(x-y)(1-Hε(ϕ(y)))dy.Keepingc1, c2, f1, and f2 fixed and minimizing the functional E(ϕ,c1,c2,f1,f2) with respect to  ϕ, we obtain the corresponding gradient descent flow equation:
(11)∂ϕ∂t=δε(ϕ)[ω((I-IG)(c1-c2))+(1-ω)((I-IL)(f1-f2))]=δε(ϕ)[ωFG(ϕ)+(1-ω)FL(ϕ)]=F(ϕ),
where
(12)δε(z)=Hε′(z)=1πεε2+z2.Like the CV and RSF models, our model is also implemented using an alternative procedure: for each iteration and the corresponding level set functionϕn, we first compute the fitting values  ci(ϕn) and fi(ϕn) and then obtain ϕn+1 by minimizing E(ϕ,c1(ϕn),c2(ϕn),f1(ϕn),f2(ϕn)) with respect to ϕ. This process is repeated until the zero-level set of ϕn+1 is exactly on the object boundary.In the following, we first discuss the properties ofFG(ϕ) and FL(ϕ) and then analyze the behavior of (11).


## 3.2. Properties ofFG(ϕ) and FL(ϕ)

For the sake of simplicity, we state and prove the properties ofFG(ϕ) and FL(ϕ) only for an image consisting of only two distinct gray levels:
(13)I(x)={g1,x∈ω,g2,x∈Ω∖ω,
where g1,g2≥0 with g1≠g2, ω and Ω∖ω represent the objects of interest and the background, respectively.Theorem 1.
LetI(x) be an image by (13). Then one has
(14)FG(ϕ)={(Mn-Nm)(g1-g2)2×[(M-N)(N-n)Hε+N((M+n)-(m+N))(1-Hε)]×(N2(M-N)2)-1,inω-(Mn-Nm)(g1-g2)2×[n(M-N)Hε+N(m-n)(1-Hε)]×(N2(M-N)2)-1,inΩ∖ω,
and so,
(15)sign(FG(ϕ))={+sign(Mn-Nm),inω,-sign(Mn-Nm),inΩ∖ω,
where
(16)M=|Ω|,m=|ω|,N=|{ϕ>0}|,n=|ω∩{ϕ>0}|,
in which |Ω| is the area of the region Ω and similarly for others.Remarks.(i) Due to the discrete nature of image, |Ω| is in fact the number of pixels in the image I(x), and similarly to others.(ii) The cases ofn=m, n=N, and n=0 correspond to the zero-level lines of ϕ(x) which are encircling the object (ω), inside the object and within the background (Ω∖ω), respectively. The cases of N≠n≠m and N=n=m correspond to the zero-level line of ϕ(x) that are partially inside the object and exactly on the object edge, respectively.(iii) The significance of (15) is that the function FG(ϕ) has the opposite sign in ω (object) and Ω∖ω (background), respectively.The proof of Theorem1 is provided in Appendix A. The following result will be used in the proof of Theorem 7, which guarantees that the evolution by (10) converges to a unique stable value after finite time.Corollary 2.
LetI(x) be an image by (13). Then one has the following.(i)
IfMn-Nm≥0, then
(17)FG(ϕ)≥(Mn-Nm)(g1-g2)2M4/4,inω,FG(ϕ)≤-(Mn-Nm)(g1-g2)2M4/4,inΩ∖ω.(ii)
IfMn-Nm≤0, then
(18)FG(ϕ)≤(Mn-Nm)(g1-g2)2M4/4,inω,FG(ϕ)≥-(Mn-Nm)(g1-g2)2M4/4,inΩ∖ω.The proof of Corollary2 is given in Appendix A.We call the property in Theorem1 an adaptive sign-changing property of FG(ϕ). Such property also holds for FL(ϕ), which can be expressed by the following theorem.Theorem 3.
LetI(x) be an image given by (13). Then one has
(19)FL(ϕ)={(Pq-Qp)(g1-g2)2×[(P-Q)(Q-q)Hε(ϕ)+Q((P+q)-(p+Q))(1-Hε(ϕ))]×(Q2(P-Q)2)-1,inω-(Pq-Qp)(g1-g2)2×[q(P-Q)Hε(ϕ)+Q(p-q)(1-Hε(ϕ))]×(Q2(P-Q)2)-1,inΩ∖ω,
and so,
(20)sign(FL(ϕ))={+sign(Pq-Qp),inω,-sign(Pq-Qp),inΩ∖ω,
where
(21)P=∫ΩKσ(x-ξ)dξ,p=∫ωKσ(x-ξ)dξ,Q=∫{ϕ>0}Kσ(x-ξ)dξ,q=∫ω∩{ϕ>0}Kσ(x-ξ)dξ.
This theorem shows that the local forceFL(ϕ) has exactly opposite sign in ω (object) and in Ω∖ω (background).The following result together with Corollary2 will be used in the proof of Theorem 7 mentioned later.Corollary 4.
Under the assumption of Theorem3, one has the following.(i)
IfPq-Qp≥0, then
(22)FL(ϕ)≥S(Pq-Qp)(g1-g2)2P4/4,inω,FL(ϕ)≤-S(Pq-Qp)(g1-g2)2P4/4,inΩ∖ω.(ii)
IfPq-Qp≤0, then
(23)FL(ϕ)≤S(Pq-Qp)(g1-g2)2P4/4,inω,FL(ϕ)≥-S(Pq-Qp)(g1-g2)2P4/4,inΩ∖ω,
where
(24)S=min{Kσ(x-ξ):ξ∈Ω)}.The proofs of Theorem3 and Corollary 4 are similar to Theorem 1 and Corollary 2, respectively; see Appendix B for details.


## 3.3. Behavior of Our Model

In this section, we analyze the behavior of our model (10) for image segmentation. We will show that the zero-level line of the evolution function starting with a bounded function finally comes to a unique steady state, which separates the object from the background.Due to the discrete nature of image, the continuous equation (10) can be discretized in space with space step 1 (implied pixel spacing), but also in time with a time step Δt as follows:
(25)ϕi,jk+1=ϕi,jk+Δt·F(ϕi,jk),0≤i≤p-1,0≤j≤q-1,
where p×q corresponds to the image size, Ii,j=I(i,j) and ϕi,jk=ϕ(kΔt,i,j) with k≥0 and ϕi,j0=ϕ0(i,j).Corollary 5.
LetMn0-N0m=|Ω||ω∩{ϕ0>0}|-|{ϕ0>0}||ω|≠0.(i)
IfMn0-N0m>0, then
(26)FG(ϕk)>(Mn0-N0m)×(g1-g2)2M4/4=A>0,inω(k≥1,k∈z+),FG(ϕk)<-(Mn0-N0m)×(g1-g2)2M4/4=-A<0,inΩ∖ω(k≥1,k∈z+).(ii)
IfMn0-N0m<0, then
(27)FG(ϕk)<(Mn0-N0m)×(g1-g2)2M4/4=A<0,inω(k≥1,k∈z+),FG(ϕk)>-(Mn0-N0m)×(g1-g2)2M4/4=-A>0,inΩ∖ω(k≥1,k∈z+),
where
(28)A=(Mn0-N0m)(g1-g2)2M4/8.We provide the proof of Corollary5 in Appendix C. Similar analysis can prove the following corollary, Corollary 6.Corollary 6.
LetPq0-Q0p≠0.(i)
IfPq0-Q0p>0, then
(29)FL(ϕk)>(Pq0-Q0p)×(g1-g2)2P4/4=B>0,inω(k≥1,k∈z+),FL(ϕk)<-(Pq0-Q0p)×(g1-g2)2P4/4=-B<0,inΩ∖ω(k≥1,k∈z+).(ii)
IfPq0-Q0p<0, then
(30)FL(ϕk)<(Pq0-Q0p)×(g1-g2)2P4/4=B<0,inω(k≥1,k∈z+),FL(ϕk)>-(Pq0-Q0p)×(g1-g2)2P4/4=-B>0,inΩ∖ω(k≥1,k∈z+),
where
(31)B=(Pq0-Q0p)(g1-g2)2P4/8.Theorem 7.
Ifϕ0(x) is a bounded function with (Mn0-N0m)(Pq0-Q0p)>0, then there exists a positive integer K such that, for k≥K,
(32)ϕi,jk+1>0,inω,ϕi,jk+1<0,inΩ∖ω,ifMn0-N0m>0,Pq0-Q0p>0,(33)ϕi,jk+1<0,inω,ϕi,jk+1>0,inΩ∖ω,ifMn0-N0m<0,Pq0-Q0p<0.We provide the proof of Theorem7 in Appendix D.The significance of Theorem7 is that if we choose ϕ0(x) such that (Mn0-N0m)(Pq0-Q0p)>0, then the zero-level line of ϕ(x) starting with such initial function ϕ0(x) can finally come to a unique steady state, which separates object from the background.Remark. Mathematically there does exist ϕ0(x) such that (Mn0-N0m)(Pq0-Q0p)≤0; in practice, however, we always can guarantee that
(34)(Mn0-N0m)(Pq0-Q0p)>0.


## 3.4. Discussion of Initial Function

Theorem7 guarantees that the proposed model computes a unique steady state regardless of the initialization; that is, the convergence of the zero-level line {ϕ=0} is irrespective of the initial function. This means that we can obtain the same zero-level line in the steady state if we choose the initial function as a bounded function ϕ0 with (Mn0-N0m)(Pq0-Q0p)>0.In applications, the initial functionϕ0 can be defined via a simple curve (closed curve or line segment) in image domain. For example, we can choose the initial function as a signed distance to a circle, which is widely used in most of image segmentation models with level set methods. For the proposed model, however, we prefer to define the initial function ϕ0(x) as a piecewise constant or constant function as follows.(1)
If the curveC is a closed curve (e.g., circle or square), then ϕ0(x) is defined by
(35)ϕ0(x)={+ρ,(x)∈inside(C),-ρ,(x)∈outside(C),
where ρ≠0 is a constant.(2)
If the curveC is a line segment that partitions the image domain Ω into two disjoint regions Ω1 and Ω2 (e.g., the left and right half regions of image domain Ω, respectively) with Ω=Ω1∪Ω2 and Ω1∩Ω2=Φ, then ϕ0(x,y) is defined by
(36)ϕ0(x)={+ρ,(x)∈Ω1,-ρ,(x)∈Ω2,
where ρ≠0 is a constant.(3)
We can also define a zero function as follows:(37)ϕ0(x)=0.Next, we prove the fact that, withϕ0(x)=0, ϕ(x) becomes a sign-changing function and satisfy the condition of (Mn1-N1m)(Pq1-Q1p)>0 after the first iteration.Theorem 8.
If the initial functionϕ0(x)=0 in Ω and
(38)c10(x)=f10(x)=2·meanx∈ΩI(x),c20(x)=f20(x)=0,
then after the first iteration, one has
(39)sign(F(0))={-1,x∈Ω1,1,x∈Ω2,
where Ω1∩Ω2=Φ, Ω1≠Φ, Ω2≠Φ.We provide the proof of Theorem8 in Appendix E.By (25) and Theorem 8, we have
(40)ϕ1(x)=δε(0)meanx∈ΩI(x)×(I-meanx∈ΩI(x)){<0,x∈Ω1,>0,x∈Ω2.
Therefore, ϕ1(x) became a sign-changing function. Then, using two distinct gray levels of (13) and the above demonstration, we have
(41)(Mn1-N1m)(Pq1-Q1p)=(m(M-N1))(p(P-Q1))=mp(M-N1)(P-Q1)>0.


## 4. Implementation and Experimental Results

### 4.1. Implementation

In tradition PDE-based methods, a certain diffusion term is usually included in the evolution equation to regularize the evolving function, but which increases the computational cost. Recently, [16] proposed a novel scheme to regularize the evolving function, that is, Gaussian filtering the evolving function after each of iterations. We adopt this scheme to regularize the evolving function u at each of iteration; that is, ϕn=Gρ*ϕn, where p controls the regularization strength. Such regularization procedure has some advantages over the traditional regularization term, such as computational efficiency and better smoothing effect; see [16] for more explanations.The main procedures of the proposed algorithm can be summarized as follows.(1)
Initialize the evolving functionϕ0(x).(2)
Computeci(ϕ) and fi(ϕ), (i=1,2) using (10).(3)
Evolve the functionϕ according to (11).(4)
Regularize the functionϕ with a Gaussian filter, that is; ϕn=Gρ*ϕn.(5)
Check whether the evolution is finished. If not, return to step 2.


### 4.2. Experimental Results

In this section, we show experimentally that the proposed model not only can provide desirable segmentation results in the presence of intensity inhomogeneity but also allows for more flexible initialization of the contour compared to the RSF and LIF models.In our numerical experiments, for our model, we choose the parameters as follows:ε=1.5 for the regularized Dirac function, Δt=0.025 (time step), h=1 (space step), ω=0.1, and n=5. The MATLAB source code of the LIF model algorithm is downloaded at http://www4.comp.polyu.edu.hk/~cslzhang/code/LIF.zip, and the RSF model algorithm is downloaded at http://www.engr.uconn.edu/~cmli/code/RSF_v0_v0.1.rar. All experiments were run under MATLAB 2010b on a PC with Dual 2.7 GHz processor.We will utilize two region overlap metrics to evaluate the performances of the three models quantitatively. The two metrics are the ratio of segmentation error (RSE) [18] and the dice similarity coefficient (DSC) [14, 18, 19]. If S1 and S2 represent a given baseline foreground region (e.g., true object) and the foreground region found by the model, respectively, then the two metrics are defined as follows:
(42)RSE=N(S1∖S2)+N(S2∖S1)N(Ω),DSC=2N(S1∩S2)N(S1)+N(S2),
where N(·) indicates the number of pixels in the enclosed region, and Ω is image domain. The closer the DSC value to 1, the better the segmentation. Since N(S1∖S2)+N(S2∖S1) is the number of the pixels mistakenly classified by the model, lower RSE means that there are fewer pixels misclassified; that is, the image can be segmented more accurately. Thus, a perfect segmentation will give DSC=1, RSE=0.In the first example (Figures1–3), we mainly verify the computation of a unique steady state of the zero-level line of ϕ, starting with three types of representative initial functions, that is, signed distance function, piecewise constant functions by (26)-(27), and zero function, respectively. The top row of Figure 1 shows the evolution of an active contour (i.e., zero-level line {ϕ=0}), with the function ϕ initialized to a signed distance function, piecewise constant functions by (26)-(27) with ρ=1, and zero function, respectively, while the bottom row shows the corresponding evolution of ϕ. With such initializations, the zero-level line of ϕ converges to the same steady state.Figure 1
Segmentations of our model for two real images withϕ initialized as different functions. Columns 1 to 4: signed distance function, piecewise constant functions by (26)-(27) with ρ=1, and zero function.Figure 2
Applications of our model to four images (slug, cell, ventriculus sinister MR, and real plane images). The curve evolution process from the initial contour (in the first column) to the final contour (in the fourth column) is shown in every row for the corresponding image.Figure 3
Comparisons of three models. Rows 1–3: RSF model, LIF model, and our model withϕ0=0. For the RSF and LIF models, initial and final contours are shown in green and red color, respectively. Columns 1–4: the 4 pictures indicate (a), (b), (c), (d), respectively.In Figure2, we show that our model has the capability of detecting multiple objects or objects with interior holes or blurred edges, only starting with a zero function. The contours (zero-level lines) evolution processes are shown in the second column to the forth column.We choose the segmentation results of the RSF and LIF models as baseline foreground regions and then compute DSC values for the corresponding images. The RSE and DSC values for the four real images are given in Table1, from which we can see that the proposed algorithm works as well as the RSF and LIF models for images with intensity inhomogeneity.Table 1
RSE and DSC for RSF and LIF models and our model.
RSF model
LIF model
RSE (%)
DSC (%)
RSE (%)
DSC (%)
Figure3(a)
0.27
97.17
2.99
98.91
Figure3(b)
2.71
96.37
0.12
99.13
Figure3(c)
0.49
87.95
0.45
88.06
Figure3(d)
0.16
97.66
0.76
98.91The experimental results shown in Figure4 validate that our method can also achieve subpixel segmentation accuracy as in [15]. As can be seen from Figures 4(b) and 4(d), both models achieve subpixel segmentation accuracy of the finger boundaries. The final contour accurately reflects the true hand shape.Figure 4
Segmentation results of both models for a hand phantom. Upper row: LIF model. Lower row: our model withϕ0=0. Column 2: zoomed view of the narrow parts in hand fingers.


## 4.1. Implementation

In tradition PDE-based methods, a certain diffusion term is usually included in the evolution equation to regularize the evolving function, but which increases the computational cost. Recently, [16] proposed a novel scheme to regularize the evolving function, that is, Gaussian filtering the evolving function after each of iterations. We adopt this scheme to regularize the evolving function u at each of iteration; that is, ϕn=Gρ*ϕn, where p controls the regularization strength. Such regularization procedure has some advantages over the traditional regularization term, such as computational efficiency and better smoothing effect; see [16] for more explanations.The main procedures of the proposed algorithm can be summarized as follows.(1)
Initialize the evolving functionϕ0(x).(2)
Computeci(ϕ) and fi(ϕ), (i=1,2) using (10).(3)
Evolve the functionϕ according to (11).(4)
Regularize the functionϕ with a Gaussian filter, that is; ϕn=Gρ*ϕn.(5)
Check whether the evolution is finished. If not, return to step 2.


## 4.2. Experimental Results

In this section, we show experimentally that the proposed model not only can provide desirable segmentation results in the presence of intensity inhomogeneity but also allows for more flexible initialization of the contour compared to the RSF and LIF models.In our numerical experiments, for our model, we choose the parameters as follows:ε=1.5 for the regularized Dirac function, Δt=0.025 (time step), h=1 (space step), ω=0.1, and n=5. The MATLAB source code of the LIF model algorithm is downloaded at http://www4.comp.polyu.edu.hk/~cslzhang/code/LIF.zip, and the RSF model algorithm is downloaded at http://www.engr.uconn.edu/~cmli/code/RSF_v0_v0.1.rar. All experiments were run under MATLAB 2010b on a PC with Dual 2.7 GHz processor.We will utilize two region overlap metrics to evaluate the performances of the three models quantitatively. The two metrics are the ratio of segmentation error (RSE) [18] and the dice similarity coefficient (DSC) [14, 18, 19]. If S1 and S2 represent a given baseline foreground region (e.g., true object) and the foreground region found by the model, respectively, then the two metrics are defined as follows:
(42)RSE=N(S1∖S2)+N(S2∖S1)N(Ω),DSC=2N(S1∩S2)N(S1)+N(S2),
where N(·) indicates the number of pixels in the enclosed region, and Ω is image domain. The closer the DSC value to 1, the better the segmentation. Since N(S1∖S2)+N(S2∖S1) is the number of the pixels mistakenly classified by the model, lower RSE means that there are fewer pixels misclassified; that is, the image can be segmented more accurately. Thus, a perfect segmentation will give DSC=1, RSE=0.In the first example (Figures1–3), we mainly verify the computation of a unique steady state of the zero-level line of ϕ, starting with three types of representative initial functions, that is, signed distance function, piecewise constant functions by (26)-(27), and zero function, respectively. The top row of Figure 1 shows the evolution of an active contour (i.e., zero-level line {ϕ=0}), with the function ϕ initialized to a signed distance function, piecewise constant functions by (26)-(27) with ρ=1, and zero function, respectively, while the bottom row shows the corresponding evolution of ϕ. With such initializations, the zero-level line of ϕ converges to the same steady state.Figure 1
Segmentations of our model for two real images withϕ initialized as different functions. Columns 1 to 4: signed distance function, piecewise constant functions by (26)-(27) with ρ=1, and zero function.Figure 2
Applications of our model to four images (slug, cell, ventriculus sinister MR, and real plane images). The curve evolution process from the initial contour (in the first column) to the final contour (in the fourth column) is shown in every row for the corresponding image.Figure 3
Comparisons of three models. Rows 1–3: RSF model, LIF model, and our model withϕ0=0. For the RSF and LIF models, initial and final contours are shown in green and red color, respectively. Columns 1–4: the 4 pictures indicate (a), (b), (c), (d), respectively.In Figure2, we show that our model has the capability of detecting multiple objects or objects with interior holes or blurred edges, only starting with a zero function. The contours (zero-level lines) evolution processes are shown in the second column to the forth column.We choose the segmentation results of the RSF and LIF models as baseline foreground regions and then compute DSC values for the corresponding images. The RSE and DSC values for the four real images are given in Table1, from which we can see that the proposed algorithm works as well as the RSF and LIF models for images with intensity inhomogeneity.Table 1
RSE and DSC for RSF and LIF models and our model.
RSF model
LIF model
RSE (%)
DSC (%)
RSE (%)
DSC (%)
Figure3(a)
0.27
97.17
2.99
98.91
Figure3(b)
2.71
96.37
0.12
99.13
Figure3(c)
0.49
87.95
0.45
88.06
Figure3(d)
0.16
97.66
0.76
98.91The experimental results shown in Figure4 validate that our method can also achieve subpixel segmentation accuracy as in [15]. As can be seen from Figures 4(b) and 4(d), both models achieve subpixel segmentation accuracy of the finger boundaries. The final contour accurately reflects the true hand shape.Figure 4
Segmentation results of both models for a hand phantom. Upper row: LIF model. Lower row: our model withϕ0=0. Column 2: zoomed view of the narrow parts in hand fingers.


## 5. Conclusion

In this study, we propose a new active contour model integrating a local intensity fitting (LIF) energy with an auxiliary global intensity fitting (GIF) energy. The LIF energy is responsible for attracting the contour toward object boundaries and is dominant near object boundaries, while the GIF energy incorporates global image information to improve the robustness to initialization of the contours. The proposed model can efficiently handle intensity inhomogeneity, while allowing for more flexible initialization and maintaining the subpixel accuracy. The utility model has the advantages of allowing for more flexible initialization of the contour and the capability of detecting multiple objects or objects with interior holes or blurred edges. But [14] has not been implemented.



---

*Source: 367086-2013-05-14.xml*