# A New Variational Model for Segmenting Objects of Interest from Color Images

Authors: Yanli Zhai; Boying Wu; Dazhi Zhang; Jiebao Sun
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

We propose a new variational model for segmenting objects of interest from color images. This model is inspired by the geodesic active contour model, the region-scalable fitting model, the weighted bounded variation model and the active contour models based on the Mumford-Shah model. In order to segment desired objects in color images, the energy functional in our model includes a discrimination function that determines whether an image pixel belongs to the desired objects or not. Compared with other active contour models, our new model cannot only avoid the usual drawback in the level set approach but also detect the objects of interest accurately. Moreover, we investigate the new model mathematically and establish the existence of the minimum to the new energy functional. Finally, numerical results show the effectiveness of our proposed model.


---

## Body

## 1. Introduction

Image segmentation is one of the fundamental problems in image processing and computer vision. Recently, variational methods have been extensively studied for image segmentation because of their flexibility in modeling and advantages in numerical implementation. In generally these methods can be categorized into two classes: edge-based models [1–3] and region-based models [4–6].Edge-based models utilize the image gradient to stop evolving contours on object boundaries. The geodesic active contour (GAC) model is one of the most well-known edge-based models. It was proposed in [2] and widely used in practice [7]. However, this model has one major disadvantage: given an initial curve, during the evolution, the energy may evolve to its bad local minimizers. Region-based models have some advantages over the edge-based models. Firstly, region-based models do not utilize the image gradient, so they have better performance for images with weak object boundaries. Secondly, these models are significantly less sensitive to the location of initial contours. The Chan-Vese (CV) model [5] is one of the most popular region-based models. It is successful for an image with two regions, each of which has a distinct mean of pixel intensities. In order to handle images with multiple regions, Vese and Chan proposed the piecewise constant (PC) models [8], in which multiple regions can be represented by multiple level set functions. And yet, these PC models are not very successful for images with intensity inhomogeneity. In order to segment images with intensity inhomogeneity efficiently, local intensity information was incorporated into the active contours models [9–11]. For example, Li et al. [9] proposed the region-scalable fitting model which draws upon local intensity means. This model is able to segment object boundaries accurately. In papers [12, 13], both local intensity means and variances are used to characterize the local intensity distribution in their proposed active contour models. However, these local intensity means and variances must be defined empirically. In addition, similar forms of local intensity means and variances were also introduced in [14] for a statistical interpretation of Mumford-Shah functional [15]. In a word, these methods have a certain capability of handling intensity inhomogeneity. Moreover, there are some models that make full use of advantages of the edge-based and the region-based models, such as the model [7] that unifies the GAC model and the CV model. This model can obtain good results when the contrast between meaningful objects and the background is low.When the aim is to detect all the objects, the above models are very efficient. However, if only the objects of interest (OOI) are concerned, the segmentation of OOI using active contour methods becomes difficult. In order to detect OOI, the properties of OOI must be known mathematically. However, it is not an easy task to find out a suitable mathematical description of these characters. Many segmentation methods have used shape or texture characters of OOI. For example, the active shape models (ASM) [16], the geodesic active contour method incorporated with shape information [17], and the supervised texture segmentation method [18].In this paper, inspired by the weighted bounded variation model [19], the geodesic active contour model, the region-scalable fitting model and the active contour models based on the Mumford-Shah model [7], we propose a new model which can be applied to segmenting objects of interest from color images. In order to segment desired objects in color images, the energy functional in our model includes a discrimination function [20] which determines whether an image pixel belongs to the desired objects or not. Moreover, our new model uses not only the edge detector which contains information concerning boundaries of desired objects but also the spatially varying fitting functions which are used to approximate the image intensities. Compared with other active contour models, our method cannot only have the advantages of both edge-based and region-based models, but also avoid the usual drawback in the level set approach (e.g., initiation and reinitiation). In particular, our new model can detect the desired objects accurately. Finally, we investigate this new model mathematically: we establish the existence of the minimum to the energy functional, analyze the property of it and implement the numerical algorithm efficiently.The remainder of the paper is organized as follows: in Section2, we show some background. In Section 3, our new model is proposed. Theoretical results, iterating schemes and experimental results are also given in this section. Finally, we conclude our paper in Section 4.


## 2. Background

### 2.1. The Region-Scalable Fitting Model

LetΩ⊂ℝ2 be the image domain and f:Ω→ℝ be a given gray-scale image. The region-scalable fitting model is defined by minimizing the following energy functional:(2.1)λ1∫Ω∫ΩKσ(x-y)|f(y)-f1(x)|2H(ϕ(y))dydx+λ2∫Ω12(|∇ϕ(x)|2-1)dx+λ3∫Ω∫ΩKσ(x-y)|f(y)-f2(x)|2(1-H(ϕ(y)))dydx+∫ΩH′(ϕ)|∇ϕ(x)|dx,
where Kσ is a Gaussian kernel and f1,f2 are two functions that fit image intensities near the point x. Moreover, ϕ is the level set function embedding the evolving active contour c={x:ϕ(x)=0} and H(ϕ) is the Heaviside function.This model does not need to reinitializeϕ periodically during the evolution because of the second term of (2.1). If λ2=0, (2.1) is equivalent to
(2.2)∂ϕ∂t=H′[div(∇ϕ|∇ϕ|)-λ1∫ΩKσ(x-y)|f(x)-f1(y)|2dy+λ3∫ΩKσ(x-y)|f(x)-f2(y)|2dy].
The steady state solution of this gradient flow is(2.3)∂ϕ∂t=div(∇ϕ|∇ϕ|)-λ1∫ΩKσ(x-y)|f(x)-f1(y)|2dy+λ3∫ΩKσ(x-y)|f(x)-f2(y)|2dy.
And it is known that (2.3) is the gradient descent flow of the following energy:(2.4)∫Ω|∇ϕ|dx+∫Ω(λ1∫ΩKσ(x-y)|f(x)-f1(y)|2dy-λ3∫ΩKσ(x-y)|f(x)-f2(y)|2dy)ϕ(x)dx.


### 2.2. Some Related Models

In [2], the geodesic active contour model (GAC) is defined by the following minimization problem:(2.5)minc{EGAC(c)=∫0L(c)gds},
where L(c) is the length of the curve c and the function g is an edge indicator function that vanishes at object boundaries. For a color image f=(f1,f2,f3), a stopping function g(x) was proposed in [21]
(2.6)g=11+β∧2,
where ∧ is the largest eigenvalue of the structure tensor metric gij in the spatial-spectral space, and(2.7)gij=(1+Rx12+Bx12+Gx12Rx1Rx2+Bx1Bx2+Gx1Gx2Rx1Rx2+Bx1Bx2+Gx1Gx21+Rx22+Bx22+Gx22),
where x=(x1,x2) and R, G, B represent the pixel values of Red, Green, and Blue after Gaussian convolution, respectively, that is, R=Gσ1*f1, G=Gσ1*f2 and B=Gσ1*f3.In order to segment the desired objects in color images, a novel stopping function depending on both the discrimination function of OOI and the image gradient was proposed in [20]:(2.8)g=11+(Gσ1*α)∧2,
where Gσ1 is a Gaussian kernel and the discrimination function α(x) can be defined by the following steps in [20].(1)
m sample pixels are chosen from the OOI. The color information of the ith sample pixel is denoted by (Ji1,Ji2,Ji3). Therefore,  J=(Jij) is a m×3 matrix.(2)
Compute the correlation coefficient matrix ofJ as ρ=(ρi,j)3×3.(3)
Compute the eigenvaluesλ̃1,λ̃2,λ̃3 of ρ and the corresponding eigenvectors w1,w2,w3. And these eigenvalues satisfy λ̃1≥λ̃2≥λ̃3.(4)
ComputeP=(P1,P2,P3)=J·(w1,w2,w3).(5)
The confidence interval is constructed as[p¯1-(sp1/m)t1-ϱ/2,p¯1+(sp1/m)t1-ϱ/2] with the degree of confidence 1-ϱ(0<ϱ<1), where p¯1=(1/m)∑i=1mPi1, sp1=(1/(m-1))∑i=1m(Pi1-p¯1)2 and t1-ϱ/2 is the t-value of t-distribution with m-1 degrees of freedom.(6)
Computea=(a1,a2,a3)=f·(w1,w2,w3). If a1(x) belongs to the above interval, α(x)=1. If a1(x) does not belong to the above interval, α(x)=0.


## 2.1. The Region-Scalable Fitting Model

LetΩ⊂ℝ2 be the image domain and f:Ω→ℝ be a given gray-scale image. The region-scalable fitting model is defined by minimizing the following energy functional:(2.1)λ1∫Ω∫ΩKσ(x-y)|f(y)-f1(x)|2H(ϕ(y))dydx+λ2∫Ω12(|∇ϕ(x)|2-1)dx+λ3∫Ω∫ΩKσ(x-y)|f(y)-f2(x)|2(1-H(ϕ(y)))dydx+∫ΩH′(ϕ)|∇ϕ(x)|dx,
where Kσ is a Gaussian kernel and f1,f2 are two functions that fit image intensities near the point x. Moreover, ϕ is the level set function embedding the evolving active contour c={x:ϕ(x)=0} and H(ϕ) is the Heaviside function.This model does not need to reinitializeϕ periodically during the evolution because of the second term of (2.1). If λ2=0, (2.1) is equivalent to
(2.2)∂ϕ∂t=H′[div(∇ϕ|∇ϕ|)-λ1∫ΩKσ(x-y)|f(x)-f1(y)|2dy+λ3∫ΩKσ(x-y)|f(x)-f2(y)|2dy].
The steady state solution of this gradient flow is(2.3)∂ϕ∂t=div(∇ϕ|∇ϕ|)-λ1∫ΩKσ(x-y)|f(x)-f1(y)|2dy+λ3∫ΩKσ(x-y)|f(x)-f2(y)|2dy.
And it is known that (2.3) is the gradient descent flow of the following energy:(2.4)∫Ω|∇ϕ|dx+∫Ω(λ1∫ΩKσ(x-y)|f(x)-f1(y)|2dy-λ3∫ΩKσ(x-y)|f(x)-f2(y)|2dy)ϕ(x)dx.


## 2.2. Some Related Models

In [2], the geodesic active contour model (GAC) is defined by the following minimization problem:(2.5)minc{EGAC(c)=∫0L(c)gds},
where L(c) is the length of the curve c and the function g is an edge indicator function that vanishes at object boundaries. For a color image f=(f1,f2,f3), a stopping function g(x) was proposed in [21]
(2.6)g=11+β∧2,
where ∧ is the largest eigenvalue of the structure tensor metric gij in the spatial-spectral space, and(2.7)gij=(1+Rx12+Bx12+Gx12Rx1Rx2+Bx1Bx2+Gx1Gx2Rx1Rx2+Bx1Bx2+Gx1Gx21+Rx22+Bx22+Gx22),
where x=(x1,x2) and R, G, B represent the pixel values of Red, Green, and Blue after Gaussian convolution, respectively, that is, R=Gσ1*f1, G=Gσ1*f2 and B=Gσ1*f3.In order to segment the desired objects in color images, a novel stopping function depending on both the discrimination function of OOI and the image gradient was proposed in [20]:(2.8)g=11+(Gσ1*α)∧2,
where Gσ1 is a Gaussian kernel and the discrimination function α(x) can be defined by the following steps in [20].(1)
m sample pixels are chosen from the OOI. The color information of the ith sample pixel is denoted by (Ji1,Ji2,Ji3). Therefore,  J=(Jij) is a m×3 matrix.(2)
Compute the correlation coefficient matrix ofJ as ρ=(ρi,j)3×3.(3)
Compute the eigenvaluesλ̃1,λ̃2,λ̃3 of ρ and the corresponding eigenvectors w1,w2,w3. And these eigenvalues satisfy λ̃1≥λ̃2≥λ̃3.(4)
ComputeP=(P1,P2,P3)=J·(w1,w2,w3).(5)
The confidence interval is constructed as[p¯1-(sp1/m)t1-ϱ/2,p¯1+(sp1/m)t1-ϱ/2] with the degree of confidence 1-ϱ(0<ϱ<1), where p¯1=(1/m)∑i=1mPi1, sp1=(1/(m-1))∑i=1m(Pi1-p¯1)2 and t1-ϱ/2 is the t-value of t-distribution with m-1 degrees of freedom.(6)
Computea=(a1,a2,a3)=f·(w1,w2,w3). If a1(x) belongs to the above interval, α(x)=1. If a1(x) does not belong to the above interval, α(x)=0.


## 3. Our Proposed Model

Letf=(f1,f2,f3):Ω→ℝ3 be a given color image where Ω⊂ℝ2 is the image domain (a bounded open set). Inspired by the GAC model, the energy functional (2.4), and the active contour models based on the Mumford-Shah model [7], our new model is constructed. This model is to minimize the following energy functional w.r.t. u∈g-BV[0,1](Ω),fi1,fi2(i=1,2,3):(3.1)E(u,fi1,fi2)=|u|g-BV(Ω)+13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)u(x)dx+13∫Ω(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)(1-u(x))dx,
where λi1,λi2(i=1,2,3) are positive constants, |u|g-BV(Ω)=sup{∫Ωu(x)divφ(x)dx:φ∈C01(Ω),|φ|≤g}, g is a diffusion coefficient defined as the formula (2.8), Kσ,Gσ1 are the Gaussian kernels, and α(x) is a discrimination function. The purpose of constructing such a discrimination function is to derive the characteristics of desired objects so that these characteristics can be shown in the energy functional. This is done by analyzing m sample pixels chosen from the desired objects. The Principal Components Analysis (PCA) and interval estimation [22] are used in the analysis. By using PCA, a new set of variables which is called principal components, is obtained. By using the interval estimation, we can define an interval that covers many samples for each principal component. In our problem, at most the first two principal components are used. Without loss of generality, we only assume the first principal component is used. An interval [a¯,b¯] for this principal component can be constructed by the interval estimation. When every pixel x of the color image is projected from its RGB values to the first principal component axis, we get a new value ũ(x). If the value is within the interval, the pixel is probabilistically regarded as a pixel in the desired objects. Thus the discrimination function α(x) based on the color information is constructed as(3.2)α(x)={1,a¯≤ũ(x)≤b¯,0,others.In order to understand more details of constructing this discrimination function, please see [23].In the following, we analyze the above model from two aspects. Firstly, for any givenu, according to the necessary condition of the minimization problem, it is known that the functions fi1(y),fi2(y)(i=1,2,3) must satisfy the following equations:(3.3)fi1(y)∫Ω(Gσ1*α(x))Kσ(x-y)u(x)dx=∫Ω(Gσ1*α(x))Kσ(x-y)fi(x)u(x)dx,fi2(y)∫Ω(Gσ1*α(x))Kσ(x-y)(1-u(x))dx=∫Ω(Gσ1*α(x))Kσ(x-y)fi(x)(1-u(x))dx,
where Kσ takes larger values at the points near the center point y, and decreases to 0 as x goes away from y. Therefore, fi1(y),fi2(y)(i=1,2,3) are allowed to vary in space.Furthermore, for any givenfi1,fi2(i=1,2,3), the model (3.1) can be converted into a simpler form. That is(3.4)minu∈g-BV[0,1](Ω)E1(u)=|u|g-BV(Ω)+13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)u(x)dx-13∫Ω(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)u(x)dx.Ifu is limited to a characteristic function 1Ωc, the energy functional (3.4) can be changed into the following form:(3.5)E1(u)=∫cgds+13∫Ωc(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)dx+13∫Ω∖Ωc(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)dx+C,
where C is a constant.In this case, the model (3.4) is equal to the following constrained minimization problem: (3.6)minc∫cgds,
when approximating fi of desired objects with spatially varying fitting functions fi1,fi2(i=1,2,3).The above analysis shows that our new model uses both the edge detector which contains information concerning boundaries of desired objects, and the spatially varying fitting functionsfi1(y),fi2(y) (i=1,2,3) which are used to approximate the image intensities of desired objects. Just because of this, our model can segment the desired objects very accurately.


### 3.1. Mathematical Results

In [7, 24], the proof of the existence of models has not been given. In the following, we will state the existence of the minimizer to the energy functional (3.4) and analyze the property of it.Theorem 3.1.
For any givenfi1,fi2∈L∞(Ω) (i=1,2,3), there exists a function u∈g-BV[0,1](Ω) minimizing the energy functional E1 in (3.4).Proof.
Let(3.7)r(x)=13(Gσ1*α(x))(∑i=13(λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy-λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)).
By fi∈L∞(Ω)(i=1,2,3),Gσ1*α∈L∞(Ω), we get r(x)∈L∞(Ω). Since u∈g-BV[0,1](Ω) and |∫Ωr(x)udx|≤M1<∞, we have
(3.8)E1(u)=|u|g-BV+∫Ωr(x)udx≥-M1.
Assume infE1=β̃≥-M1 and {un} is the minimizing sequence of (3.4) in g-BV[0,1](Ω), that is, limn→∞E1(un)=β̃. So there is a positive constant M2 such that
(3.9)|E1(un)|=||un|g-BV+∫Ωr(x)undx|<M2.
The structure tensor metricgij is symmetric positive and ∧ is the largest eigenvalue of gij. Thus(3.10)0≤∧≤trace(gij)≤2+c1∥f1∥L∞(Ω)2+c2∥f2∥L∞(Ω)2+c3∥f3∥L∞(Ω)2,
where ci≥0(i=1,2,3).
Sincefi∈L∞(Ω)(i=1,2,3), Gσ1*α≥0 and Gσ1*α∈L∞(Ω), g≥1/(1+C2) where C2 is a positive constant. That is, g≥δ(δ=1/(1+C2)). Then(3.11)δ|un|BV≤M2+|∫Ωr(x)undx|≤M2+M1=M3.
So |un|BV is bounded.
Because{un}∈g-BV[0,1](Ω)={u:u(x)∈[0,1]foreveryx∈Ω}⋂g-BV(Ω), the BV-norm of un is bounded. Thus, there is a subsequence, also denoted by {un}, and u*∈BV(Ω) such that un→u* strongly in L1(Ω).
Moreover, according to the formulaun→u* in L1(Ω), we know that there is a subsequence, also denoted by {un}, satisfying limn→∞un(x)=u*(x) almost everywhere for x∈Ω. Since un(x)∈[0,1] for any x∈Ω, u*(x)∈[0,1] almost everywhere for x∈Ω.
AssumeΩ0⊂Ω, s.t. m(Ω0)=0 and u*(x)∈[0,1], limn→∞un(x)=u*(x) for any x∈Ω∖Ω0⊂Ω. Let(3.12)u¯*(x)={u*(x),ifx∈Ω∖Ω0,1,ifx∈Ω0.
Then
(3.13)∫Ω|un-u¯*|dx=∫Ω∖Ω0|un-u¯*|dx=∫Ω∖Ω0|un-u*|dx=∫Ω|un-u*|dx.
That is, un→u¯* strongly in L1(Ω) and u¯*(x)∈[0,1] for every x∈Ω.
Furthermore, by the lower semicontinuity for theg-BV space, we get(3.14)|u¯*|g-BV≤limn→∞inf|un|g-BV<∞.
So u¯*∈g-BV[0,1](Ω). According to the dominated convergence theorem, we know
(3.15)∫Ωr(x)u¯*dx=limn→∞∫Ωr(x)undx.
The formulas (3.14)-(3.15) imply the weak lower semicontinuity of the energy functional E1(3.16)E1(u¯*)≤limn→∞infE1(un)=β̃.
Therefore, the infimum is attained by u¯*∈g-BV[0,1](Ω) and it is a minimum of the energy functional E1.Similar to [7, 25], we can obtain a property of minimizers to the energy functional (3.4).Theorem 3.2.
Letfi(x)(i=1,2,3),g(x)∈[0,1]. For any given fi1,fi2 (i=1,2,3), if u(x) is any minimum of the energy functional E1 defined in (3.4), then for almost every μ∈[0,1], the characteristic function 1Ωc=1{x:u(x)>μ} is also a global minimum of the functional E1 where c is the boundary of the set Ωc.In addition, according to the above theorem, we knowΩc={x:u(x)>μ} is a minimizer of the following minimization problem
(3.17)minc∫cgds+13∫Ωc(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)dx+13∫Ω∖Ωc(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)dx.


### 3.2. Numerical Implementation

In the numerical algorithm, we do not deal with the new variational model (3.1) directly since too many equations need to be computed. In order to improve computational efficiency, we use the algorithm framework of the paper [24] to deal with our model. That is, we minimize the energy functional E by alternating the following steps.(1)
Consideringu fixed, compute fi1 and fi2(i=1,2,3) by using the formula (3.3).(2)
Consideringfi1 and fi2(i=1,2,3) fixed, update u by using the iterative schemes of the minimization problem (3.4). When a steady state is found, the final segmentation is obtained by thresholding u at any level in [0,1] (in our experiments, we choose μ=0.5).In the following, we will give the iterative schemes of the minimization problem (3.4). To solve this minimization problem, we firstly change it into the following unconstrained minimization problem:(3.18)minuE2(u)=|u|g-BV(Ω)+ζ∫Ων(u)dx+13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)u(x)dx-13∫Ω(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)u(x)dx,
where ν(u)=max{0,2|u-0.5|-1} is an exact penalty function provided that the constant ζ is chosen large enough. According to the paper [7], it is known that this unconstrained minimization problem has the same set of minimizers as the minimization problem (3.4). The energy functional E2 is convex, so it doesn't possess local minimizers. Hence, any minimizer of E2 is global.Based on [26–28], a convex regularization is used (3.19)minu,v13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)vdx+|u|g-BV(Ω)+12θ∥u-v∥2+∫Ω-13(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)v+ζν(v)dx,
where θ is chosen to be small enough so that we almost have u=v. Thus, we consider the iterative schemes of the above weak approximation as the iterative schemes of step (2). Since this functional is convex, its minimum can be computed by minimizing this functional with respect to u and v separately. That is(1)
v being fixed, we search for u as a solution of
(3.20)minu|u|g-BV(Ω)+12θ∥u-v∥2,(2)
u being fixed, we search for v as a solution of
(3.21)minv13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)vdx+12θ∥u-v∥2+∫Ω-13(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)v+ζν(v)dx.According to [29], we know the solution of (3.20) can be given by(3.22)u=v-θdivp⃗,
where p⃗=(p1,p2) is given by (3.23)g(x)∇(θdivp⃗-v)-|∇(θdivp⃗-v)|p⃗=0.
The previous equation can be solved by the fixed point method(3.24)p⃗n+1=p⃗n+δt∇(divp⃗n-v/θ)1+δt/g(x)|∇(divp⃗n-v/θ)|,p⃗0=(0,0).
Moreover, the solution of (3.21) is given by(3.25)v=min{max{u(x)-θr(x),0},1},
where
(3.26)r(x)=13(Gσ1*α(x))(∑i=13(λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy-λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)).In this paper, our model is solved numerically by alternative minimization schemes and not by using the classical Euler-Lagrange equations method. When using Euler-Lagrange equations, one must consider∫Ωg(x)|∇u|2+ɛ2dx instead of the weighted TV-norm, where the small parameter ɛ is necessary to prevent numerical instabilities. The direct result of this regularization parameter is the obligation to use a small temporal step to ensure a correct minimization process. Thus, a large number of iterations to reach the steady state solution is necessary. In other words, although it is correct, the segmentation process remains slow. In the algorithm of our paper, the iteration process (3.22)–(3.24) is very similar to the standard Chambolle's projection algorithm. The only difference is the appearance of the factor g(x). This iteration process has two important advantages: firstly, it can be implemented very fast since it uses many useful convex optimization tools. Secondly, it is more faithful to the continuous formulation of the energy since it does not use the additional artificial parameter ɛ. Therefore, our algorithm is fast and doesn't need so many iterations to reach the steady state. Moreover, the algorithm of our paper is numerically more efficient than classical level set methods since the level set methods need to initialize the active contour in a distance function and reinitialize it periodically during the evolution. In fact, the algorithm of our paper belongs to the algorithm framework mentioned in [24]. And there are many papers [7, 28] that have analyzed the complexity of this kind of algorithms. So we just analyze it simply in this paper. If more details are wanted, please see [7, 24, 28].


### 3.3. Experimental Results

Our method can be implemented on many synthetic and natural images. First, we consider two different cases Figures1(a) and 2(a). Figure 1(a) is an image containing three different objects. The blue rectangle is our desired object. Figure 2(a) is not a simple image where the bird is our interested object. Figures 1(b), 2(b), 1(c), and 2(c) display the final active contours and u got by using our new model, respectively. From these experimental results, we find our model can detect OOI, regardless of other objects.An example of an image segmentation. (a) The original image (our desired object is the blue rectangle). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1).
(a)(b)(c)An example of an image segmentation. (a) The original image (our desired object is the bird). (b–c) The final active contour andu obtained by using our model (m=6, δt=1/8,  λi1=λi2=0.001(i=1,2,3), θ=0.1).
(a)(b)(c)In the following, we compare our new model with the model [23] in Figures 3, 4, and 5. Figures 3(a), 4(a), and 5(a) are the original images. Our purpose is to segment the plane in Figure 3(a), the yellow leaves in Figure 4(a) and the green leaves of oranges in Figure 5(a). Figures 3(b), 4(b), 5(b), 3(c), 4(c), and 5(c) denote the final active contours and u obtained by using our new model, respectively. According to these results, we see that our model can detect the boundaries of the desired objects very accurately. And these results are hard to achieve by using the model in [23].An example of an image segmentation. (a) The original image (our desired object is the plane). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1). (d) The final active contour obtained by using the model in [23].
(a)(b)(c)(d)An example of an image segmentation. (a) The original image (our desired objects are the yellow leaves). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1). (d) The final active contour obtained by using the model in [23].
(a)(b)(c)(d)An example of an image segmentation. (a) The original image (our desired objects are the leaves). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1). (d) The final active contour obtained by using the model in [23].
(a)(b)(c)(d)


## 3.1. Mathematical Results

In [7, 24], the proof of the existence of models has not been given. In the following, we will state the existence of the minimizer to the energy functional (3.4) and analyze the property of it.Theorem 3.1.
For any givenfi1,fi2∈L∞(Ω) (i=1,2,3), there exists a function u∈g-BV[0,1](Ω) minimizing the energy functional E1 in (3.4).Proof.
Let(3.7)r(x)=13(Gσ1*α(x))(∑i=13(λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy-λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)).
By fi∈L∞(Ω)(i=1,2,3),Gσ1*α∈L∞(Ω), we get r(x)∈L∞(Ω). Since u∈g-BV[0,1](Ω) and |∫Ωr(x)udx|≤M1<∞, we have
(3.8)E1(u)=|u|g-BV+∫Ωr(x)udx≥-M1.
Assume infE1=β̃≥-M1 and {un} is the minimizing sequence of (3.4) in g-BV[0,1](Ω), that is, limn→∞E1(un)=β̃. So there is a positive constant M2 such that
(3.9)|E1(un)|=||un|g-BV+∫Ωr(x)undx|<M2.
The structure tensor metricgij is symmetric positive and ∧ is the largest eigenvalue of gij. Thus(3.10)0≤∧≤trace(gij)≤2+c1∥f1∥L∞(Ω)2+c2∥f2∥L∞(Ω)2+c3∥f3∥L∞(Ω)2,
where ci≥0(i=1,2,3).
Sincefi∈L∞(Ω)(i=1,2,3), Gσ1*α≥0 and Gσ1*α∈L∞(Ω), g≥1/(1+C2) where C2 is a positive constant. That is, g≥δ(δ=1/(1+C2)). Then(3.11)δ|un|BV≤M2+|∫Ωr(x)undx|≤M2+M1=M3.
So |un|BV is bounded.
Because{un}∈g-BV[0,1](Ω)={u:u(x)∈[0,1]foreveryx∈Ω}⋂g-BV(Ω), the BV-norm of un is bounded. Thus, there is a subsequence, also denoted by {un}, and u*∈BV(Ω) such that un→u* strongly in L1(Ω).
Moreover, according to the formulaun→u* in L1(Ω), we know that there is a subsequence, also denoted by {un}, satisfying limn→∞un(x)=u*(x) almost everywhere for x∈Ω. Since un(x)∈[0,1] for any x∈Ω, u*(x)∈[0,1] almost everywhere for x∈Ω.
AssumeΩ0⊂Ω, s.t. m(Ω0)=0 and u*(x)∈[0,1], limn→∞un(x)=u*(x) for any x∈Ω∖Ω0⊂Ω. Let(3.12)u¯*(x)={u*(x),ifx∈Ω∖Ω0,1,ifx∈Ω0.
Then
(3.13)∫Ω|un-u¯*|dx=∫Ω∖Ω0|un-u¯*|dx=∫Ω∖Ω0|un-u*|dx=∫Ω|un-u*|dx.
That is, un→u¯* strongly in L1(Ω) and u¯*(x)∈[0,1] for every x∈Ω.
Furthermore, by the lower semicontinuity for theg-BV space, we get(3.14)|u¯*|g-BV≤limn→∞inf|un|g-BV<∞.
So u¯*∈g-BV[0,1](Ω). According to the dominated convergence theorem, we know
(3.15)∫Ωr(x)u¯*dx=limn→∞∫Ωr(x)undx.
The formulas (3.14)-(3.15) imply the weak lower semicontinuity of the energy functional E1(3.16)E1(u¯*)≤limn→∞infE1(un)=β̃.
Therefore, the infimum is attained by u¯*∈g-BV[0,1](Ω) and it is a minimum of the energy functional E1.Similar to [7, 25], we can obtain a property of minimizers to the energy functional (3.4).Theorem 3.2.
Letfi(x)(i=1,2,3),g(x)∈[0,1]. For any given fi1,fi2 (i=1,2,3), if u(x) is any minimum of the energy functional E1 defined in (3.4), then for almost every μ∈[0,1], the characteristic function 1Ωc=1{x:u(x)>μ} is also a global minimum of the functional E1 where c is the boundary of the set Ωc.In addition, according to the above theorem, we knowΩc={x:u(x)>μ} is a minimizer of the following minimization problem
(3.17)minc∫cgds+13∫Ωc(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)dx+13∫Ω∖Ωc(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)dx.


## 3.2. Numerical Implementation

In the numerical algorithm, we do not deal with the new variational model (3.1) directly since too many equations need to be computed. In order to improve computational efficiency, we use the algorithm framework of the paper [24] to deal with our model. That is, we minimize the energy functional E by alternating the following steps.(1)
Consideringu fixed, compute fi1 and fi2(i=1,2,3) by using the formula (3.3).(2)
Consideringfi1 and fi2(i=1,2,3) fixed, update u by using the iterative schemes of the minimization problem (3.4). When a steady state is found, the final segmentation is obtained by thresholding u at any level in [0,1] (in our experiments, we choose μ=0.5).In the following, we will give the iterative schemes of the minimization problem (3.4). To solve this minimization problem, we firstly change it into the following unconstrained minimization problem:(3.18)minuE2(u)=|u|g-BV(Ω)+ζ∫Ων(u)dx+13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)u(x)dx-13∫Ω(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)u(x)dx,
where ν(u)=max{0,2|u-0.5|-1} is an exact penalty function provided that the constant ζ is chosen large enough. According to the paper [7], it is known that this unconstrained minimization problem has the same set of minimizers as the minimization problem (3.4). The energy functional E2 is convex, so it doesn't possess local minimizers. Hence, any minimizer of E2 is global.Based on [26–28], a convex regularization is used (3.19)minu,v13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)vdx+|u|g-BV(Ω)+12θ∥u-v∥2+∫Ω-13(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)v+ζν(v)dx,
where θ is chosen to be small enough so that we almost have u=v. Thus, we consider the iterative schemes of the above weak approximation as the iterative schemes of step (2). Since this functional is convex, its minimum can be computed by minimizing this functional with respect to u and v separately. That is(1)
v being fixed, we search for u as a solution of
(3.20)minu|u|g-BV(Ω)+12θ∥u-v∥2,(2)
u being fixed, we search for v as a solution of
(3.21)minv13∫Ω(Gσ1*α(x))(∑i=13λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy)vdx+12θ∥u-v∥2+∫Ω-13(Gσ1*α(x))(∑i=13λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)v+ζν(v)dx.According to [29], we know the solution of (3.20) can be given by(3.22)u=v-θdivp⃗,
where p⃗=(p1,p2) is given by (3.23)g(x)∇(θdivp⃗-v)-|∇(θdivp⃗-v)|p⃗=0.
The previous equation can be solved by the fixed point method(3.24)p⃗n+1=p⃗n+δt∇(divp⃗n-v/θ)1+δt/g(x)|∇(divp⃗n-v/θ)|,p⃗0=(0,0).
Moreover, the solution of (3.21) is given by(3.25)v=min{max{u(x)-θr(x),0},1},
where
(3.26)r(x)=13(Gσ1*α(x))(∑i=13(λi1∫ΩKσ(x-y)(fi(x)-fi1(y))2dy-λi2∫ΩKσ(x-y)(fi(x)-fi2(y))2dy)).In this paper, our model is solved numerically by alternative minimization schemes and not by using the classical Euler-Lagrange equations method. When using Euler-Lagrange equations, one must consider∫Ωg(x)|∇u|2+ɛ2dx instead of the weighted TV-norm, where the small parameter ɛ is necessary to prevent numerical instabilities. The direct result of this regularization parameter is the obligation to use a small temporal step to ensure a correct minimization process. Thus, a large number of iterations to reach the steady state solution is necessary. In other words, although it is correct, the segmentation process remains slow. In the algorithm of our paper, the iteration process (3.22)–(3.24) is very similar to the standard Chambolle's projection algorithm. The only difference is the appearance of the factor g(x). This iteration process has two important advantages: firstly, it can be implemented very fast since it uses many useful convex optimization tools. Secondly, it is more faithful to the continuous formulation of the energy since it does not use the additional artificial parameter ɛ. Therefore, our algorithm is fast and doesn't need so many iterations to reach the steady state. Moreover, the algorithm of our paper is numerically more efficient than classical level set methods since the level set methods need to initialize the active contour in a distance function and reinitialize it periodically during the evolution. In fact, the algorithm of our paper belongs to the algorithm framework mentioned in [24]. And there are many papers [7, 28] that have analyzed the complexity of this kind of algorithms. So we just analyze it simply in this paper. If more details are wanted, please see [7, 24, 28].


## 3.3. Experimental Results

Our method can be implemented on many synthetic and natural images. First, we consider two different cases Figures1(a) and 2(a). Figure 1(a) is an image containing three different objects. The blue rectangle is our desired object. Figure 2(a) is not a simple image where the bird is our interested object. Figures 1(b), 2(b), 1(c), and 2(c) display the final active contours and u got by using our new model, respectively. From these experimental results, we find our model can detect OOI, regardless of other objects.An example of an image segmentation. (a) The original image (our desired object is the blue rectangle). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1).
(a)(b)(c)An example of an image segmentation. (a) The original image (our desired object is the bird). (b–c) The final active contour andu obtained by using our model (m=6, δt=1/8,  λi1=λi2=0.001(i=1,2,3), θ=0.1).
(a)(b)(c)In the following, we compare our new model with the model [23] in Figures 3, 4, and 5. Figures 3(a), 4(a), and 5(a) are the original images. Our purpose is to segment the plane in Figure 3(a), the yellow leaves in Figure 4(a) and the green leaves of oranges in Figure 5(a). Figures 3(b), 4(b), 5(b), 3(c), 4(c), and 5(c) denote the final active contours and u obtained by using our new model, respectively. According to these results, we see that our model can detect the boundaries of the desired objects very accurately. And these results are hard to achieve by using the model in [23].An example of an image segmentation. (a) The original image (our desired object is the plane). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1). (d) The final active contour obtained by using the model in [23].
(a)(b)(c)(d)An example of an image segmentation. (a) The original image (our desired objects are the yellow leaves). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1). (d) The final active contour obtained by using the model in [23].
(a)(b)(c)(d)An example of an image segmentation. (a) The original image (our desired objects are the leaves). (b–c) The final active contour andu obtained by using our model (m=7, δt=1/8, λi1=λi2=0.001(i=1,2,3), θ=0.1). (d) The final active contour obtained by using the model in [23].
(a)(b)(c)(d)


## 4. Conclusion

This paper describes a new variational model for segmenting desired objects in color images. This model is inspired by the GAC model, the region-scalable fitting model, the weighted bounded variation model and the active contour models based on the Mumford-Shah model. In order to segment objects of interest from color images, the energy functional in our model includes a discrimination function which determines whether an image pixel belongs to the desired objects or not. Compared with other active contour models, our new model cannot only avoid the usual drawback in the level set approach but also detect the desired objects accurately. Our numerical results confirm the effectiveness of our algorithm. Moreover, we establish the existence of the minimum to the new energy functional and analyze the property of it.



---

*Source: 950405-2010-06-28.xml*