# Variational Level Set Method for Two-Stage Image Segmentation Based on Morphological Gradients

Authors: Zemin Ren
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We use variational level set method and transition region extraction techniques to achieve image segmentation task. The proposed scheme is done by two steps. We first develop a novel algorithm to extract transition region based on the morphological gradient. After this, we integrate the transition region into a variational level set framework and develop a novel geometric active contour model, which include an external energy based on transition region and fractional order edge indicator function. The external energy is used to drive the zero level set toward the desired image features, such as object boundaries. Due to this external energy, the proposed model allows for more flexible initialization. The fractional order edge indicator function is incorporated into the length regularization term to diminish the influence of noise. Moreover, internal energy is added into the proposed model to penalize the deviation of the level set function from a signed distance function. The results evolution of the level set function is the gradient flow that minimizes the overall energy functional. The proposed model has been applied to both synthetic and real images with promising results.


---

## Body

## 1. Introduction

Image segmentation is a functional problem and complex task in image processing and computer vision. Its goal is to change the representation of a given image into something that is more meaningful and easier to analyze. However, image segmentation is yet a difficult task since it requires a semantic understanding of the image. To perform the image segmentation task, many successful techniques including geometric active contour models using the level set method [1] have been presented.Active contour model, proposed by Kass et al. [2], has been proved to be an efficient framework for image segmentation. However, these models suffer from the sensitivity to initial conditions and the difficulties associated with topological changes like the merging and splitting of the evolving curve. In order to overcome these problems, implicit active contour models, that is, active contour models in a level set formulation, have been proposed for image segmentation. The basic idea of the implicit active model is that an active contour is implicitly represented by the zero level set of a function in higher dimension (called level set function), and then the level set function is deformed according to an evolution partial differential equation. In the view of mathematics, implicit active contour models can be categorized into two categories: one is pure PDE model [1, 3–6] whose evolution equation is directly constructed; another is the variational level set model [7–10] whose evolution equation is derived from the minimization problem for the energy functional defined on the level set function. In this study, we focus on the variational level set methods.The variational level set methods have been well established and widely used in many applications of image processing. In field of image segmentation, it can be formulated by minimizing energy functional defined on the level function by gradient decent method or other methods. The level set functions evolve in keeping with a partial differential equation, which is derived from the minimization of the energy functional. The popular piecewise constant (PC) model [7] is a typical variational-based level set method, which aims to minimize the Mumford-Shah functional [8]. Chan and Vese [7] utilize the global image statistics inside and outside the evolving curve rather than the gradients on the boundaries to formulate and achieve good performance in image segmentation task. However, the PC model usually fails to segment the images with intensity inhomogeneity since it assumes that the intensities in each region always remain constant. To overcome the limitation of the PC model, many efficient schemes have been proposed in related literatures [9–11]. Wang et al. formulate a local PC model using local statistical function and extend this model for texture segmentation by an extended structure tensor [11]. Li et al. proposed an implicit active contours model based on local binary fitting energy, which used the local information as constraint and well work on the image with intensity inhomogeneities [12, 13]. He et al. propose an improved region-scalable fitting model based on the “mollifying” kernel and local entropy [14]. There is a problem appearing in these implicit models: contour initialization. This means that the segmentation results generally depend on the selection of initial contours. Recently, Li et al. design a new external energy to make the proposed model robust to initialization or free of manual initialization, which is denoted as MCY model in this paper [15].To solve the limitation of variational level set method, many researchers try to propose efficient schemes combined with other methods [16]. In this paper, we incorporate the transition region into the variational framework, which is different from other methods. Transition region is geometrically located between object and background and composed of pixels having intermediate gray levels between those of object and background [17]. Due to this characteristic, transition region-based thresholding has been developed quickly in image segmentation area. Gerbrands demonstrates the existence of transition region in an image for the first time. After this, Zhang and Gerbrands introduce the transition region into thresholding segmentation method [18]. Yan et al. present a local entropy based transition region extraction and thresholding method [19]. Li and Liu propose a modified local entropy based transition region extraction method by integrating local complexity and local variance [20]. These references show that transition region has been used in thresholding segmentation well. However, thresholding techniques generally introduce discontinuous contour. In this work, we incorporate the transition region into the variational framework and develop a novel image segmentation method.In this paper, we propose a two-stage variational level set model for image segmentation. The proposed segmentation scheme is done in two steps. We first develop a novel algorithm to extract transition region based on the morphological gradient. After this, we integrate the transition region into a variational level set framework and develop a novel geometric active contour model. The energy functional for the proposed model consists of three parts: external term, regularization term, and internal term. The external energy based on transition region is used to drive the zero level set toward the desired image features, such as object boundaries. Due to this external energy, the proposed model allows for more flexible initialization. The fractional order edge indicator function is incorporated into the length regularization term to diminish the influence of noise. Moreover, internal energy is added into the proposed model to penalize the deviation of the level set function from a signed distance function. The results evolution of the level set function is the gradient flow that minimizes the overall energy functional. The proposed model has been applied to both synthetic and real images with promising results.The remainder of the paper is organized as follows. In Section2, we simply review the related knowledge of morphological gradients and fractional order differentiation. The proposed model and analysis are introduced in Section 3. We present the numerical scheme and experimental results in Section 4. The conclusion is presented in Section 5.


## 2. Related Works

### 2.1. Morphological Gradients

Mathematical morphology is a theory and technique for the analysis and processing of geometrical structures, based on set theory, lattice theory, topology, and random functions, which can be dated back to 1964. It is a typical nonlinear method in image processing, which has been widely applied in image denoising, edge detection, and image restoration. In mathematical morphology, dilation and erosion are two fundamental operators. Other operators can be defined by these operators, like morphological gradient. LetI
(
x
)
:
Ω
→
R be a grayscale image and let S be a grayscale structuring element which is used to detect the entire image. Thus, the dilation operator of image I
(
x
) can be computed by
(1)
(
I
⊕
S
)
(
x
)
=
sup
⁡
y
∈
E
[
I
(
y
)
+
S
(
x
-
y
)
]
;
and the erosion operator of image I
(
x
) can be computed by
(2)
(
I
⊙
S
)
(
x
)
=
inf
⁡
y
∈
E
[
I
(
y
)
+
S
(
y
-
x
)
]
.
According to the definition of dilation, there are two results for image performed by this operator: (
1
) if the value of structuring element is positive, the output image is brighter than the input one; (
2
) it will reduce the black detail part of the input image. Conversely, there are two results after the erosion: (
1
) if the value of structuring element is positive, the output image is blacker than the input one; (
2
) it will reduce the bright detail part of the input image.Morphological gradients can be defined by the difference between the dilation and the erosion of a given image. In other words, the morphological gradients are of enhancing variation of pixel intensity in a given neighborhood. Thus, it is useful for edge detection and image segmentation. For a grayscale imageI, the morphological gradient has the following formulation:
(3)
Grad
(
I
)
≡
(
I
⊕
S
)
-
(
I
⊙
S
)
.
We also denote G
=
Grad
(
I
).


### 2.2. Fractional Order Differentiation Formulation

As a generalization of integer order derivative, the fractional order derivative can date back to correspondence between Leibniz and Hospital in 1695. The fractional order derivative has been recognized as a powerful model mythology and applied in many fields such as image processing. Up to now, the reader is aware that more than one fractional order derivative exists in literatures. In this paper, we give the following formulation about the fractional order differentiation.For a given continuous function of a single variablef
(
x
)
∈
L
2
(
R
), the Fourier transform of f
(
x
) is defined as
(4)
f
^
(
w
)
=
∫
R
f
(
t
)
exp
⁡
⁡
(
-
j
w
t
)
d
t
,
where j
=
-
1. In this paper, we use f
n
(
t
) instead of f
(
n
)
(
t
) to represent the nth derivative of f
(
t
). According to the differentiation property of Fourier transform, the equivalent formulation of the nth derivative in the frequency domain is
(5)
f
n
(
t
)
⟷
F
(
f
n
(
t
)
)
=
(
j
w
)
n
f
^
(
w
)
,
where ↔ denotes the Fourier transform pair. Clearly, the right-hand side of the above expression is meaningful for any number n. Thus, we can define the fractional order derivative of f with the order ν whose real part is greater than zero in the Fourier domain by
(6)
D
ν
f
(
t
)
=
F
-
1
(
(
j
w
)
ν
f
^
(
w
)
)
,
Re
ν
>
0
,
where F
-
1 is the inverse Fourier transform operator. Similarly, fractional order partial derivatives of f
(
x
,
y
)
∈
L
2
(
R
2
)
∩
C
(
R
2
) are defined as
(7)
D
x
ν
f
(
x
,
y
)
=
F
-
1
(
(
j
w
1
)
ν
f
^
(
w
1
,
w
2
)
)
,
D
y
ν
f
(
x
,
y
)
=
F
-
1
(
(
j
w
2
)
ν
f
^
(
w
1
,
w
2
)
)
.Naturally, the fractional gradient operator can be defined as(8)
∇
ν
f
=
(
D
x
ν
f
,
D
y
ν
f
)
with the associated norm |
∇
ν
f
|
=
(
D
x
ν
f
)
2
+
(
D
y
ν
f
)
2. According to this norm, the fractional order edge indicator function is defined by the following formulation generally
(9)
g
(
|
∇
ν
I
|
)
=
1
1
+
c
2
|
∇
ν
I
|
2
.It can be used for edge detection and also enhance selectivity and robustness the presence of noise. The fractional order edge indicator function is proposed in image denoising [21] and image upsampling [22].


## 2.1. Morphological Gradients

Mathematical morphology is a theory and technique for the analysis and processing of geometrical structures, based on set theory, lattice theory, topology, and random functions, which can be dated back to 1964. It is a typical nonlinear method in image processing, which has been widely applied in image denoising, edge detection, and image restoration. In mathematical morphology, dilation and erosion are two fundamental operators. Other operators can be defined by these operators, like morphological gradient. LetI
(
x
)
:
Ω
→
R be a grayscale image and let S be a grayscale structuring element which is used to detect the entire image. Thus, the dilation operator of image I
(
x
) can be computed by
(1)
(
I
⊕
S
)
(
x
)
=
sup
⁡
y
∈
E
[
I
(
y
)
+
S
(
x
-
y
)
]
;
and the erosion operator of image I
(
x
) can be computed by
(2)
(
I
⊙
S
)
(
x
)
=
inf
⁡
y
∈
E
[
I
(
y
)
+
S
(
y
-
x
)
]
.
According to the definition of dilation, there are two results for image performed by this operator: (
1
) if the value of structuring element is positive, the output image is brighter than the input one; (
2
) it will reduce the black detail part of the input image. Conversely, there are two results after the erosion: (
1
) if the value of structuring element is positive, the output image is blacker than the input one; (
2
) it will reduce the bright detail part of the input image.Morphological gradients can be defined by the difference between the dilation and the erosion of a given image. In other words, the morphological gradients are of enhancing variation of pixel intensity in a given neighborhood. Thus, it is useful for edge detection and image segmentation. For a grayscale imageI, the morphological gradient has the following formulation:
(3)
Grad
(
I
)
≡
(
I
⊕
S
)
-
(
I
⊙
S
)
.
We also denote G
=
Grad
(
I
).


## 2.2. Fractional Order Differentiation Formulation

As a generalization of integer order derivative, the fractional order derivative can date back to correspondence between Leibniz and Hospital in 1695. The fractional order derivative has been recognized as a powerful model mythology and applied in many fields such as image processing. Up to now, the reader is aware that more than one fractional order derivative exists in literatures. In this paper, we give the following formulation about the fractional order differentiation.For a given continuous function of a single variablef
(
x
)
∈
L
2
(
R
), the Fourier transform of f
(
x
) is defined as
(4)
f
^
(
w
)
=
∫
R
f
(
t
)
exp
⁡
⁡
(
-
j
w
t
)
d
t
,
where j
=
-
1. In this paper, we use f
n
(
t
) instead of f
(
n
)
(
t
) to represent the nth derivative of f
(
t
). According to the differentiation property of Fourier transform, the equivalent formulation of the nth derivative in the frequency domain is
(5)
f
n
(
t
)
⟷
F
(
f
n
(
t
)
)
=
(
j
w
)
n
f
^
(
w
)
,
where ↔ denotes the Fourier transform pair. Clearly, the right-hand side of the above expression is meaningful for any number n. Thus, we can define the fractional order derivative of f with the order ν whose real part is greater than zero in the Fourier domain by
(6)
D
ν
f
(
t
)
=
F
-
1
(
(
j
w
)
ν
f
^
(
w
)
)
,
Re
ν
>
0
,
where F
-
1 is the inverse Fourier transform operator. Similarly, fractional order partial derivatives of f
(
x
,
y
)
∈
L
2
(
R
2
)
∩
C
(
R
2
) are defined as
(7)
D
x
ν
f
(
x
,
y
)
=
F
-
1
(
(
j
w
1
)
ν
f
^
(
w
1
,
w
2
)
)
,
D
y
ν
f
(
x
,
y
)
=
F
-
1
(
(
j
w
2
)
ν
f
^
(
w
1
,
w
2
)
)
.Naturally, the fractional gradient operator can be defined as(8)
∇
ν
f
=
(
D
x
ν
f
,
D
y
ν
f
)
with the associated norm |
∇
ν
f
|
=
(
D
x
ν
f
)
2
+
(
D
y
ν
f
)
2. According to this norm, the fractional order edge indicator function is defined by the following formulation generally
(9)
g
(
|
∇
ν
I
|
)
=
1
1
+
c
2
|
∇
ν
I
|
2
.It can be used for edge detection and also enhance selectivity and robustness the presence of noise. The fractional order edge indicator function is proposed in image denoising [21] and image upsampling [22].


## 3. The Proposed Two-Stage Model

In this section, we will present two steps for the proposed segmentation scheme: transition region extraction method and variational level set segmentation model.


### 3.1. Transition Region Extraction Method Based on Morphological Gradient

Recently, transition region based method has been developed in related references such as local entropy based transition region extraction method and gray level difference based transition region extraction method. In this work, we present the morphological gradients based transition region extraction method and then use the mean of transition region to compute the weight of the external energy in the variational level set model. Detailed process of the proposed transition extraction algorithm is shown in Algorithm1.Algorithm 1:Image transition region extraction.
Input: Data I and grayscale structuring element S.
Processing:  (1) Compute the morphological gradients G of the data I according to (3);
(2) Design the following threshold G
T for the transition region extraction;
G
T
=
κ
×
G
max
⁡ with G
max
⁡
=
max
∀
(
i
,
j
)
⁡
G
(
i
,
j
)
(3) Extract the transition region TR via the following way:
TR
(
i
,
j
)
=
{
1
if
G
(
i
,
j
)
≥
G
T
0
otherwie;
Output: the mean of transition region T
* for TR.The process of extraction is different from other methods, which is established on domain of the morphological gradients. Meanwhile, we will deal with transition region in a different way and incorporate it into the variational framework.


### 3.2. The Proposed Variational Model Based on Transition Region

In level set methods, a closed curveC
⊂
Ω is presented implicitly by the zero level set of a Lipschitz function ϕ
:
Ω
→
R, with the following properties:
(10)
ϕ
(
x
,
y
)
>
0
,
(
x
,
y
)
∈
inside
(
C
)
ϕ
(
x
,
y
)
=
0
,
(
x
,
y
)
∈
C
ϕ
(
x
,
y
)
<
0
,
(
x
,
y
)
∈
outside
(
C
)
.
The proposed model is established on level set function, and its evolution is controlled by three forces, an external force, adaptive regularizing force, and internal force.


#### 3.2.1. External Energy Term Based on Morphological Gradient

In image segmentation, an external energy can be used to move the zero level curve toward the objects boundaries, which depend on the image information under the level set method. In this subsection, we define an external energyA
(
ϕ
) based on morphological gradients to drive the level set function to have opposite sign besides the edges.For the given imageI
:
Ω
→
R, we can define the following function:
(11)
w
(
I
,
T
*
)
=
G
σ
*
(
T
*
-
I
)
,
where T
* is the mean of transition region, obtained by Algorithm 1 as described in Section 3.1, and G
σ is the Gaussian kernel function with the standard deviation σ. In the following, we would show that the weight function w
(
I
,
T
*
) has the opposite sign beside edges in the transition region, which can be named as adaptive sign property. Without loss of generality, we assume that the gray-level of the object is brighter than that in the background. In the transition region, if the point (
x
,
y
) is located at the region of the object, the gray value of the given image I for the point (
x
,
y
) is lower than the mean of transition region T
*. Thus, the weight function of weighted area functional w
(
I
,
T
*
) for the point (
x
,
y
) is negative. In contrast, the value of the weight function is positive in the region of the background. By the above analysis, we can obtain that the value of the weight function is of adaptive sign property.Next, we propose a novel variational formulation based on this weighted function. This will serve as an external energy for the proposed model, which can be defined as(12)
A
(
ϕ
)
=
∫
Ω
w
(
I
,
T
*
)
H
(
-
ϕ
)
d
x
d
y
,
where w
(
I
,
T
*
) is the weight function defined by (11) and H
(
·
) is the Heaviside function. The gradient flow of A
(
ϕ
) is given by the Gateaux derivative:
(13)
∇
A
(
ϕ
)
=
-
w
(
I
,
T
*
)
δ
(
ϕ
)
.According to the adaptive sign property ofw
(
I
,
T
*
), the level set function ϕ moves with speed of |
w
(
I
,
T
*
)
δ
(
ϕ
)
| in the opposite direction (up or down) automatically, which coincides with the property of level set function.


#### 3.2.2. Adaptive Length Regularization Term Based on Fractional Order Edge Indicator Function

As the existing level set evolution models, for our model it is also necessary to add a regularizing term into the evolution to smooth the level set functionϕ. In the spirit of the Mumford-Shah (MS) functional, most of models focus on penalizing the length of contours as the regularization. For our model, we pursue the idea of MS functional and define the following adaptive regularizing force based on the fractional order edge indicator function for the level set function ϕ:
(14)
L
(
ϕ
)
=
∫
Ω
1
1
+
c
2
|
∇
ν
I
|
2
|
∇
H
(
ϕ
)
|
d
x
d
y
=
∫
Ω
1
1
+
c
2
|
∇
ν
I
|
2
δ
(
ϕ
)
|
∇
ϕ
|
d
x
d
y
,
where δ
(
·
) is Dirac delta function and I
:
Ω
→
R is the observed image. For the minimum point ϕ, the gradient flow of the objective energy functional L
(
ϕ
) is given by
(15)
∇
L
(
ϕ
)
=
-
δ
(
ϕ
)
div
⁡
(
g
(
|
∇
ν
I
|
)
∇
ϕ
|
∇
ϕ
|
)
,
where g
(
|
∇
ν
I
|
)
=
1
/
(
1
+
c
2
|
∇
ν
I
|
2
). In order to understand the geometric meaning of the objective energy functional L
(
ϕ
), we suppose that zero level set of ϕ can be represented by a differentiable parameterized curve C
(
p
),  p
∈
[
0,1
]. Then, functional L
(
ϕ
) in (14) shows the length of the zero level curve of ϕ in the conformal metric d
s
=
g
(
C
(
p
)
)
|
C
′
(
p
)
|
d
p. Thus, the zero level curves move with speed g
(
∇
ν
I
) and stop on the desired boundary. In this term, the fractional order edge indicator function can release the influence of noise for computing the length of the zero level curves, which is different from the generally edge indicator function as used in PM model [23]. The effect of this edge indicator function is from the fractional order differentiation. There are two basic ideals behind the selection of fractional order differentiation. First, the intensity inhomogeneity is slowly varying in the image domain. Its spectrum in frequency domain will be concentrated in the low-frequency area. Fractional order differentiation is able to preserve and enhance the low frequency information, which is superior to the first order differentiation. Second, fractional order differentiation can improve immunity to noise, which can be interpreted in terms of robustness to noise in general. This property is very important for medical images and infrared images that suffer from noise. According to these reasons, we incorporate this edge indicator function into the proposed model.


#### 3.2.3. Internal Energy Functional

In many situations, the level set function will develop shocks, very sharp and/or flat shape during the evolution, which in turn makes further computation highly inaccurate in numerical approximations [9]. To avoid these problems, it is necessary to keep the level set function as an approximate signed distance function in the processing of the evolution, especially in the neighborhood around the zero level set [24]. An alternative approach is to solve the eikonal equation |
∇
ϕ
|
=
1. Some fast algorithms have been proposed for solving this eikonal equation. However, these algorithms are generally time-consuming. In this paper, we add an extra internal energy to replace solving the eikonal equation, which can be written as
(16)
P
(
ϕ
)
=
1
2
∫
Ω
(
|
∇
ϕ
(
x
)
|
-
1
)
2
d
x
with the gradient flow
(17)
∇
P
(
ϕ
)
=
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
-
Δ
ϕ
.This extra internal energy avoids directly using the reinitialization step to keep the level set function as a signed distance function. Actually, this term is more like a metric to characterize how close a functionϕ is to a signed distance function, which plays an important role in our model.


#### 3.2.4. The Proposed Energy Functional

With the above defined forces, that is, adaptive length regularization termL
(
ϕ
), weighted external energy A
(
ϕ
), and P
(
ϕ
), the overall objective energy functional for the level set function ϕ is expressed as
(18)
E
(
ϕ
)
=
λ
L
(
ϕ
)
+
γ
A
(
ϕ
)
+
μ
P
(
ϕ
)
=
λ
∫
Ω
g
δ
(
ϕ
)
|
∇
ϕ
|
d
x
d
y
+
γ
∫
Ω
w
(
I
,
T
*
)
H
(
-
ϕ
)
d
x
d
y
+
μ
2
∫
Ω
(
|
∇
ϕ
|
-
1
)
2
d
x
d
y
.By the linearity of gradient operator, the gradient flow of the functionalE
(
ϕ
) is given by
(19)
∇
E
(
ϕ
)
=
-
λ
δ
(
ϕ
)
div
⁡
(
g
∇
ϕ
|
∇
ϕ
|
)
-
γ
w
(
I
,
T
*
)
δ
(
ϕ
)
-
μ
[
Δ
ϕ
-
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
]
.We obtain the gradient decent flow equation for the functionalE
(
ϕ
):
(20)
∂
ϕ
∂
t
=
λ
δ
(
ϕ
)
div
⁡
(
g
∇
ϕ
|
∇
ϕ
|
)
+
γ
w
(
I
,
T
*
)
δ
(
ϕ
)
+
μ
[
Δ
ϕ
-
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
]
.For practical and feasible implementation,H
ε
(
x
) is chosen as a noncompactly supported, smooth, and strictly monotone approximation of H
(
x
), which can be formulated as
(21)
H
ε
(
x
)
=
1
2
|
1
+
2
π
arctan
|
x
ε
|
|
,
ε
⟶
0
.The regularized versionδ
ε
(
x
) of the Dirac function δ
(
x
) is correspondingly written as follows:
(22)
δ
ε
(
x
)
=
1
π
·
ε
ε
2
+
x
2
.


#### 3.2.5. Further Analysis of the Proposed Energy Functional

In this subsection, we would verify the effectiveness of the two new terms incorporated in the new level set method. The external energy term is a force which moves the zero level set curve toward the object boundaries. If without this term, the proposed model does not have any other forces to finish this work, and its experimental results will be like the first column of Figure1. The adaptive regularization term based on fractional order differentiation is used to diminish the influence of noise and smooth the level set function. If we reduce this term, the proposed model will not finish segmentation task for images with strong noise. We give the experimental results on two images: original image and degraded image with strong noise, as shown in Figure 1. As we can see, the proposed model cannot segment image with strong noise when it loses the help of the adaptive regularization term.Segmentations results on the regularization term. First column: original image and noisy images; second column: the proposed model withoutL
(
ϕ
); third column: the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.1. Transition Region Extraction Method Based on Morphological Gradient

Recently, transition region based method has been developed in related references such as local entropy based transition region extraction method and gray level difference based transition region extraction method. In this work, we present the morphological gradients based transition region extraction method and then use the mean of transition region to compute the weight of the external energy in the variational level set model. Detailed process of the proposed transition extraction algorithm is shown in Algorithm1.Algorithm 1:Image transition region extraction.
Input: Data I and grayscale structuring element S.
Processing:  (1) Compute the morphological gradients G of the data I according to (3);
(2) Design the following threshold G
T for the transition region extraction;
G
T
=
κ
×
G
max
⁡ with G
max
⁡
=
max
∀
(
i
,
j
)
⁡
G
(
i
,
j
)
(3) Extract the transition region TR via the following way:
TR
(
i
,
j
)
=
{
1
if
G
(
i
,
j
)
≥
G
T
0
otherwie;
Output: the mean of transition region T
* for TR.The process of extraction is different from other methods, which is established on domain of the morphological gradients. Meanwhile, we will deal with transition region in a different way and incorporate it into the variational framework.


## 3.2. The Proposed Variational Model Based on Transition Region

In level set methods, a closed curveC
⊂
Ω is presented implicitly by the zero level set of a Lipschitz function ϕ
:
Ω
→
R, with the following properties:
(10)
ϕ
(
x
,
y
)
>
0
,
(
x
,
y
)
∈
inside
(
C
)
ϕ
(
x
,
y
)
=
0
,
(
x
,
y
)
∈
C
ϕ
(
x
,
y
)
<
0
,
(
x
,
y
)
∈
outside
(
C
)
.
The proposed model is established on level set function, and its evolution is controlled by three forces, an external force, adaptive regularizing force, and internal force.


### 3.2.1. External Energy Term Based on Morphological Gradient

In image segmentation, an external energy can be used to move the zero level curve toward the objects boundaries, which depend on the image information under the level set method. In this subsection, we define an external energyA
(
ϕ
) based on morphological gradients to drive the level set function to have opposite sign besides the edges.For the given imageI
:
Ω
→
R, we can define the following function:
(11)
w
(
I
,
T
*
)
=
G
σ
*
(
T
*
-
I
)
,
where T
* is the mean of transition region, obtained by Algorithm 1 as described in Section 3.1, and G
σ is the Gaussian kernel function with the standard deviation σ. In the following, we would show that the weight function w
(
I
,
T
*
) has the opposite sign beside edges in the transition region, which can be named as adaptive sign property. Without loss of generality, we assume that the gray-level of the object is brighter than that in the background. In the transition region, if the point (
x
,
y
) is located at the region of the object, the gray value of the given image I for the point (
x
,
y
) is lower than the mean of transition region T
*. Thus, the weight function of weighted area functional w
(
I
,
T
*
) for the point (
x
,
y
) is negative. In contrast, the value of the weight function is positive in the region of the background. By the above analysis, we can obtain that the value of the weight function is of adaptive sign property.Next, we propose a novel variational formulation based on this weighted function. This will serve as an external energy for the proposed model, which can be defined as(12)
A
(
ϕ
)
=
∫
Ω
w
(
I
,
T
*
)
H
(
-
ϕ
)
d
x
d
y
,
where w
(
I
,
T
*
) is the weight function defined by (11) and H
(
·
) is the Heaviside function. The gradient flow of A
(
ϕ
) is given by the Gateaux derivative:
(13)
∇
A
(
ϕ
)
=
-
w
(
I
,
T
*
)
δ
(
ϕ
)
.According to the adaptive sign property ofw
(
I
,
T
*
), the level set function ϕ moves with speed of |
w
(
I
,
T
*
)
δ
(
ϕ
)
| in the opposite direction (up or down) automatically, which coincides with the property of level set function.


### 3.2.2. Adaptive Length Regularization Term Based on Fractional Order Edge Indicator Function

As the existing level set evolution models, for our model it is also necessary to add a regularizing term into the evolution to smooth the level set functionϕ. In the spirit of the Mumford-Shah (MS) functional, most of models focus on penalizing the length of contours as the regularization. For our model, we pursue the idea of MS functional and define the following adaptive regularizing force based on the fractional order edge indicator function for the level set function ϕ:
(14)
L
(
ϕ
)
=
∫
Ω
1
1
+
c
2
|
∇
ν
I
|
2
|
∇
H
(
ϕ
)
|
d
x
d
y
=
∫
Ω
1
1
+
c
2
|
∇
ν
I
|
2
δ
(
ϕ
)
|
∇
ϕ
|
d
x
d
y
,
where δ
(
·
) is Dirac delta function and I
:
Ω
→
R is the observed image. For the minimum point ϕ, the gradient flow of the objective energy functional L
(
ϕ
) is given by
(15)
∇
L
(
ϕ
)
=
-
δ
(
ϕ
)
div
⁡
(
g
(
|
∇
ν
I
|
)
∇
ϕ
|
∇
ϕ
|
)
,
where g
(
|
∇
ν
I
|
)
=
1
/
(
1
+
c
2
|
∇
ν
I
|
2
). In order to understand the geometric meaning of the objective energy functional L
(
ϕ
), we suppose that zero level set of ϕ can be represented by a differentiable parameterized curve C
(
p
),  p
∈
[
0,1
]. Then, functional L
(
ϕ
) in (14) shows the length of the zero level curve of ϕ in the conformal metric d
s
=
g
(
C
(
p
)
)
|
C
′
(
p
)
|
d
p. Thus, the zero level curves move with speed g
(
∇
ν
I
) and stop on the desired boundary. In this term, the fractional order edge indicator function can release the influence of noise for computing the length of the zero level curves, which is different from the generally edge indicator function as used in PM model [23]. The effect of this edge indicator function is from the fractional order differentiation. There are two basic ideals behind the selection of fractional order differentiation. First, the intensity inhomogeneity is slowly varying in the image domain. Its spectrum in frequency domain will be concentrated in the low-frequency area. Fractional order differentiation is able to preserve and enhance the low frequency information, which is superior to the first order differentiation. Second, fractional order differentiation can improve immunity to noise, which can be interpreted in terms of robustness to noise in general. This property is very important for medical images and infrared images that suffer from noise. According to these reasons, we incorporate this edge indicator function into the proposed model.


### 3.2.3. Internal Energy Functional

In many situations, the level set function will develop shocks, very sharp and/or flat shape during the evolution, which in turn makes further computation highly inaccurate in numerical approximations [9]. To avoid these problems, it is necessary to keep the level set function as an approximate signed distance function in the processing of the evolution, especially in the neighborhood around the zero level set [24]. An alternative approach is to solve the eikonal equation |
∇
ϕ
|
=
1. Some fast algorithms have been proposed for solving this eikonal equation. However, these algorithms are generally time-consuming. In this paper, we add an extra internal energy to replace solving the eikonal equation, which can be written as
(16)
P
(
ϕ
)
=
1
2
∫
Ω
(
|
∇
ϕ
(
x
)
|
-
1
)
2
d
x
with the gradient flow
(17)
∇
P
(
ϕ
)
=
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
-
Δ
ϕ
.This extra internal energy avoids directly using the reinitialization step to keep the level set function as a signed distance function. Actually, this term is more like a metric to characterize how close a functionϕ is to a signed distance function, which plays an important role in our model.


### 3.2.4. The Proposed Energy Functional

With the above defined forces, that is, adaptive length regularization termL
(
ϕ
), weighted external energy A
(
ϕ
), and P
(
ϕ
), the overall objective energy functional for the level set function ϕ is expressed as
(18)
E
(
ϕ
)
=
λ
L
(
ϕ
)
+
γ
A
(
ϕ
)
+
μ
P
(
ϕ
)
=
λ
∫
Ω
g
δ
(
ϕ
)
|
∇
ϕ
|
d
x
d
y
+
γ
∫
Ω
w
(
I
,
T
*
)
H
(
-
ϕ
)
d
x
d
y
+
μ
2
∫
Ω
(
|
∇
ϕ
|
-
1
)
2
d
x
d
y
.By the linearity of gradient operator, the gradient flow of the functionalE
(
ϕ
) is given by
(19)
∇
E
(
ϕ
)
=
-
λ
δ
(
ϕ
)
div
⁡
(
g
∇
ϕ
|
∇
ϕ
|
)
-
γ
w
(
I
,
T
*
)
δ
(
ϕ
)
-
μ
[
Δ
ϕ
-
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
]
.We obtain the gradient decent flow equation for the functionalE
(
ϕ
):
(20)
∂
ϕ
∂
t
=
λ
δ
(
ϕ
)
div
⁡
(
g
∇
ϕ
|
∇
ϕ
|
)
+
γ
w
(
I
,
T
*
)
δ
(
ϕ
)
+
μ
[
Δ
ϕ
-
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
]
.For practical and feasible implementation,H
ε
(
x
) is chosen as a noncompactly supported, smooth, and strictly monotone approximation of H
(
x
), which can be formulated as
(21)
H
ε
(
x
)
=
1
2
|
1
+
2
π
arctan
|
x
ε
|
|
,
ε
⟶
0
.The regularized versionδ
ε
(
x
) of the Dirac function δ
(
x
) is correspondingly written as follows:
(22)
δ
ε
(
x
)
=
1
π
·
ε
ε
2
+
x
2
.


### 3.2.5. Further Analysis of the Proposed Energy Functional

In this subsection, we would verify the effectiveness of the two new terms incorporated in the new level set method. The external energy term is a force which moves the zero level set curve toward the object boundaries. If without this term, the proposed model does not have any other forces to finish this work, and its experimental results will be like the first column of Figure1. The adaptive regularization term based on fractional order differentiation is used to diminish the influence of noise and smooth the level set function. If we reduce this term, the proposed model will not finish segmentation task for images with strong noise. We give the experimental results on two images: original image and degraded image with strong noise, as shown in Figure 1. As we can see, the proposed model cannot segment image with strong noise when it loses the help of the adaptive regularization term.Segmentations results on the regularization term. First column: original image and noisy images; second column: the proposed model withoutL
(
ϕ
); third column: the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.2.1. External Energy Term Based on Morphological Gradient

In image segmentation, an external energy can be used to move the zero level curve toward the objects boundaries, which depend on the image information under the level set method. In this subsection, we define an external energyA
(
ϕ
) based on morphological gradients to drive the level set function to have opposite sign besides the edges.For the given imageI
:
Ω
→
R, we can define the following function:
(11)
w
(
I
,
T
*
)
=
G
σ
*
(
T
*
-
I
)
,
where T
* is the mean of transition region, obtained by Algorithm 1 as described in Section 3.1, and G
σ is the Gaussian kernel function with the standard deviation σ. In the following, we would show that the weight function w
(
I
,
T
*
) has the opposite sign beside edges in the transition region, which can be named as adaptive sign property. Without loss of generality, we assume that the gray-level of the object is brighter than that in the background. In the transition region, if the point (
x
,
y
) is located at the region of the object, the gray value of the given image I for the point (
x
,
y
) is lower than the mean of transition region T
*. Thus, the weight function of weighted area functional w
(
I
,
T
*
) for the point (
x
,
y
) is negative. In contrast, the value of the weight function is positive in the region of the background. By the above analysis, we can obtain that the value of the weight function is of adaptive sign property.Next, we propose a novel variational formulation based on this weighted function. This will serve as an external energy for the proposed model, which can be defined as(12)
A
(
ϕ
)
=
∫
Ω
w
(
I
,
T
*
)
H
(
-
ϕ
)
d
x
d
y
,
where w
(
I
,
T
*
) is the weight function defined by (11) and H
(
·
) is the Heaviside function. The gradient flow of A
(
ϕ
) is given by the Gateaux derivative:
(13)
∇
A
(
ϕ
)
=
-
w
(
I
,
T
*
)
δ
(
ϕ
)
.According to the adaptive sign property ofw
(
I
,
T
*
), the level set function ϕ moves with speed of |
w
(
I
,
T
*
)
δ
(
ϕ
)
| in the opposite direction (up or down) automatically, which coincides with the property of level set function.


## 3.2.2. Adaptive Length Regularization Term Based on Fractional Order Edge Indicator Function

As the existing level set evolution models, for our model it is also necessary to add a regularizing term into the evolution to smooth the level set functionϕ. In the spirit of the Mumford-Shah (MS) functional, most of models focus on penalizing the length of contours as the regularization. For our model, we pursue the idea of MS functional and define the following adaptive regularizing force based on the fractional order edge indicator function for the level set function ϕ:
(14)
L
(
ϕ
)
=
∫
Ω
1
1
+
c
2
|
∇
ν
I
|
2
|
∇
H
(
ϕ
)
|
d
x
d
y
=
∫
Ω
1
1
+
c
2
|
∇
ν
I
|
2
δ
(
ϕ
)
|
∇
ϕ
|
d
x
d
y
,
where δ
(
·
) is Dirac delta function and I
:
Ω
→
R is the observed image. For the minimum point ϕ, the gradient flow of the objective energy functional L
(
ϕ
) is given by
(15)
∇
L
(
ϕ
)
=
-
δ
(
ϕ
)
div
⁡
(
g
(
|
∇
ν
I
|
)
∇
ϕ
|
∇
ϕ
|
)
,
where g
(
|
∇
ν
I
|
)
=
1
/
(
1
+
c
2
|
∇
ν
I
|
2
). In order to understand the geometric meaning of the objective energy functional L
(
ϕ
), we suppose that zero level set of ϕ can be represented by a differentiable parameterized curve C
(
p
),  p
∈
[
0,1
]. Then, functional L
(
ϕ
) in (14) shows the length of the zero level curve of ϕ in the conformal metric d
s
=
g
(
C
(
p
)
)
|
C
′
(
p
)
|
d
p. Thus, the zero level curves move with speed g
(
∇
ν
I
) and stop on the desired boundary. In this term, the fractional order edge indicator function can release the influence of noise for computing the length of the zero level curves, which is different from the generally edge indicator function as used in PM model [23]. The effect of this edge indicator function is from the fractional order differentiation. There are two basic ideals behind the selection of fractional order differentiation. First, the intensity inhomogeneity is slowly varying in the image domain. Its spectrum in frequency domain will be concentrated in the low-frequency area. Fractional order differentiation is able to preserve and enhance the low frequency information, which is superior to the first order differentiation. Second, fractional order differentiation can improve immunity to noise, which can be interpreted in terms of robustness to noise in general. This property is very important for medical images and infrared images that suffer from noise. According to these reasons, we incorporate this edge indicator function into the proposed model.


## 3.2.3. Internal Energy Functional

In many situations, the level set function will develop shocks, very sharp and/or flat shape during the evolution, which in turn makes further computation highly inaccurate in numerical approximations [9]. To avoid these problems, it is necessary to keep the level set function as an approximate signed distance function in the processing of the evolution, especially in the neighborhood around the zero level set [24]. An alternative approach is to solve the eikonal equation |
∇
ϕ
|
=
1. Some fast algorithms have been proposed for solving this eikonal equation. However, these algorithms are generally time-consuming. In this paper, we add an extra internal energy to replace solving the eikonal equation, which can be written as
(16)
P
(
ϕ
)
=
1
2
∫
Ω
(
|
∇
ϕ
(
x
)
|
-
1
)
2
d
x
with the gradient flow
(17)
∇
P
(
ϕ
)
=
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
-
Δ
ϕ
.This extra internal energy avoids directly using the reinitialization step to keep the level set function as a signed distance function. Actually, this term is more like a metric to characterize how close a functionϕ is to a signed distance function, which plays an important role in our model.


## 3.2.4. The Proposed Energy Functional

With the above defined forces, that is, adaptive length regularization termL
(
ϕ
), weighted external energy A
(
ϕ
), and P
(
ϕ
), the overall objective energy functional for the level set function ϕ is expressed as
(18)
E
(
ϕ
)
=
λ
L
(
ϕ
)
+
γ
A
(
ϕ
)
+
μ
P
(
ϕ
)
=
λ
∫
Ω
g
δ
(
ϕ
)
|
∇
ϕ
|
d
x
d
y
+
γ
∫
Ω
w
(
I
,
T
*
)
H
(
-
ϕ
)
d
x
d
y
+
μ
2
∫
Ω
(
|
∇
ϕ
|
-
1
)
2
d
x
d
y
.By the linearity of gradient operator, the gradient flow of the functionalE
(
ϕ
) is given by
(19)
∇
E
(
ϕ
)
=
-
λ
δ
(
ϕ
)
div
⁡
(
g
∇
ϕ
|
∇
ϕ
|
)
-
γ
w
(
I
,
T
*
)
δ
(
ϕ
)
-
μ
[
Δ
ϕ
-
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
]
.We obtain the gradient decent flow equation for the functionalE
(
ϕ
):
(20)
∂
ϕ
∂
t
=
λ
δ
(
ϕ
)
div
⁡
(
g
∇
ϕ
|
∇
ϕ
|
)
+
γ
w
(
I
,
T
*
)
δ
(
ϕ
)
+
μ
[
Δ
ϕ
-
div
⁡
(
∇
ϕ
|
∇
ϕ
|
)
]
.For practical and feasible implementation,H
ε
(
x
) is chosen as a noncompactly supported, smooth, and strictly monotone approximation of H
(
x
), which can be formulated as
(21)
H
ε
(
x
)
=
1
2
|
1
+
2
π
arctan
|
x
ε
|
|
,
ε
⟶
0
.The regularized versionδ
ε
(
x
) of the Dirac function δ
(
x
) is correspondingly written as follows:
(22)
δ
ε
(
x
)
=
1
π
·
ε
ε
2
+
x
2
.


## 3.2.5. Further Analysis of the Proposed Energy Functional

In this subsection, we would verify the effectiveness of the two new terms incorporated in the new level set method. The external energy term is a force which moves the zero level set curve toward the object boundaries. If without this term, the proposed model does not have any other forces to finish this work, and its experimental results will be like the first column of Figure1. The adaptive regularization term based on fractional order differentiation is used to diminish the influence of noise and smooth the level set function. If we reduce this term, the proposed model will not finish segmentation task for images with strong noise. We give the experimental results on two images: original image and degraded image with strong noise, as shown in Figure 1. As we can see, the proposed model cannot segment image with strong noise when it loses the help of the adaptive regularization term.Segmentations results on the regularization term. First column: original image and noisy images; second column: the proposed model withoutL
(
ϕ
); third column: the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)


## 4. Numerical Algorithm and Results

### 4.1. Numerical Implementation of the Model

The partial differential equation in the continuous domain presented in (20) can be implemented using a simple finite difference method in numerical scheme. All the spatial derivatives ∂
ϕ
/
∂
x and ∂
ϕ
/
∂
y are approximated by the central difference and the temporal partial derivatives ∂
ϕ
/
∂
t are approximated by the forward difference. To approximate (20) numerically, we recall first the notations in the finite differences schemes. Let Δ
t be the time step, let h be the space step, and let (
x
i
,
y
j
)
=
(
i
h
,
j
h
) be the grid points. We also denote ϕ
i
,
j
n
=
ϕ
(
x
i
,
y
j
,
n
Δ
t
) as an approximation of the level set function ϕ
(
x
,
y
,
t
). The central differences of spatial partial derivatives are written as the following notations:
(23)
Δ
x
ϕ
i
,
j
=
ϕ
i
+
1
,
j
-
ϕ
i
-
1
,
j
2
h
,
Δ
y
ϕ
i
,
j
=
ϕ
i
,
j
+
1
-
ϕ
i
,
j
-
1
2
h
.Then, the evolution equation of our model in (20) can be discredited using the forward difference as the following formulation:
(24)
ϕ
i
,
j
n
+
1
-
ϕ
i
,
j
n
Δ
t
=
Q
(
ϕ
i
,
j
n
)
,
where Q
(
ϕ
i
,
j
n
) is the numerical approximation of the right hand side of (20). The corresponding curvature κ
=
div
⁡
(
∇
ϕ
/
|
∇
ϕ
|
) in the Q
(
ϕ
i
,
j
n
) can be discredited as
(25)
k
i
,
j
n
=
Δ
x
[
Δ
x
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
+
Δ
y
[
Δ
y
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
.Then, (20) is implemented as follows:
(26)
ϕ
i
,
j
n
+
1
-
ϕ
i
,
j
n
Δ
t
=
λ
δ
(
ϕ
i
,
j
n
)
l
i
,
j
n
+
γ
w
(
I
,
T
*
)
i
,
j
δ
(
ϕ
i
,
j
n
)
+
μ
(
Δ
x
(
Δ
x
ϕ
i
,
j
n
)
+
Δ
y
(
Δ
y
ϕ
i
,
j
n
)
-
k
i
,
j
n
)
,
where
(27)
l
i
,
j
n
=
Δ
x
[
1
1
+
c
2
|
∇
ν
I
i
,
j
|
·
Δ
x
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
+
Δ
y
[
1
1
+
c
2
|
∇
ν
I
i
,
j
|
·
Δ
y
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
.


### 4.2. Experimental Results

In this subsection, we present the experimental results of our model on a variety of synthetic and real images and compare them with the MCY model [15] and RSF model [13]. Unless otherwise specified, the level set function ϕ
(
x
,
y
,
t
) is simply initialized to a constant function ϕ
0
(
x
,
y
)
=
0 for all the experiments. We use the same parameters of time step Δ
t
=
5, λ
=
15, and γ
=
30 for all the images in this paper.The first experiment is to perform the proposed mode on synthetic images with increased strength of noise shown in Figure2. This group contains four images: original image and three noisy images generated by adding Gaussian white noise with variances 0.02, 0.1, and 0.2. We use κ
=
0.8 for these four images. The results of our method are shown in the second row of Figure 2. From this, we can see that the proposed model can finish the task of segment well for not only original image but also degraded images with strong noise. In addition, we compare our model with the MCY model to test the sensitivity to initialization. The experimental results for ϕ
0
=
-
1,0
,
1 are shown in Figure 3. We observe that the proposed model successfully extracts the object for ϕ
0
=
-
1,0
,
1, while the MCY model does not finish the task of segment for the case of ϕ
0
=
-
1,0. Thus, the proposed model is of more flexible initialization. From Figure 3, we also see that the robustness to noise for the proposed model is better than that for the MCY model from the third and fourth rows of Figure 3.Segmentation results by both models: the proposed model and RSF model on the synthetic image with different noise. First row: original image and noisy images (Gaussian white noise with zero mean with different variances: 0.02, 0.1, and 0.2 for (b)–(d)); second row: the proposed model; third row: RSF model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)The comparisons of the MCY model (the first and the third row) and the proposed model (the second and the fourth row) on segmenting the synthetic images in Figures2(a) and 2(b), starting with different constant functions (ϕ
0
=
-
1,0
,
1 from left column to right column).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)The second experiment is to perform the proposed model on four infrared ship images shown in Figure4. This type of image has common characteristic, that is, burry images with high strength of noise. This brings very great difficulties for image segmentation. Here, we advise that the value of κ belong to interval [
0.1,0.5
] and use κ
=
0.15,0.25,0.25,0.5 for the given image. From the first experiment, we can find that the proposed model is robust to the noise, which is used to extract the object for this kind image. The results of the MCY model and the RSF model are also shown in Figure 4. We observed that the proposed model achieves the segmentation task.Segmentation results on four ship images by three models. First row: original images; second row: the MCY model; third row: the RSF model; fourth row: the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)The third experiment is to perform the proposed model on another kind of image, skin image. To analyze skin lesions, it is necessary to accurately locate and isolate the lesions. However, it is a challenging task since the skin image contains many things, such as cutaneous melanoma, grid texture, and hair on the human skin surface. In this experiment, we test the proposed model on six skin images. The segmentation results by the proposed model and the other two models are shown in Figure5. We can see that the MCY model and the RSF model cannot deal with this kind of image, while the proposed model achieves satisfactory segmentation results for these six images.Segmentation results on skin images by three models. First column: original image; second column: the MCY model; third column: the RSF model; fourth column: the proposed model.
(a)
(b)
(c)
(d)The final experiment is to perform the proposed model on the other three real images from different modalities: a breast cyst image with complex tissue background, an aerial image with cognitive contours, and a MR image of ventriculus cordis with complex structure and weak boundary. We list the results by the proposed model, the MCY model, and the RSF model in Figure6. As we can see, the proposed model extracts the objects for these three images. They illustrate the advantages of the proposed model: the ability to deal with weak boundary and complex background.Segmentation results on other types of image by three models. First row: an aerial image, a breast cyst image, and a MR image of ventriculus cordis; second row: the MCY model; third row: the RSF model; fourth row: the proposed model.
(a)
(b)
(c)
(d)


## 4.1. Numerical Implementation of the Model

The partial differential equation in the continuous domain presented in (20) can be implemented using a simple finite difference method in numerical scheme. All the spatial derivatives ∂
ϕ
/
∂
x and ∂
ϕ
/
∂
y are approximated by the central difference and the temporal partial derivatives ∂
ϕ
/
∂
t are approximated by the forward difference. To approximate (20) numerically, we recall first the notations in the finite differences schemes. Let Δ
t be the time step, let h be the space step, and let (
x
i
,
y
j
)
=
(
i
h
,
j
h
) be the grid points. We also denote ϕ
i
,
j
n
=
ϕ
(
x
i
,
y
j
,
n
Δ
t
) as an approximation of the level set function ϕ
(
x
,
y
,
t
). The central differences of spatial partial derivatives are written as the following notations:
(23)
Δ
x
ϕ
i
,
j
=
ϕ
i
+
1
,
j
-
ϕ
i
-
1
,
j
2
h
,
Δ
y
ϕ
i
,
j
=
ϕ
i
,
j
+
1
-
ϕ
i
,
j
-
1
2
h
.Then, the evolution equation of our model in (20) can be discredited using the forward difference as the following formulation:
(24)
ϕ
i
,
j
n
+
1
-
ϕ
i
,
j
n
Δ
t
=
Q
(
ϕ
i
,
j
n
)
,
where Q
(
ϕ
i
,
j
n
) is the numerical approximation of the right hand side of (20). The corresponding curvature κ
=
div
⁡
(
∇
ϕ
/
|
∇
ϕ
|
) in the Q
(
ϕ
i
,
j
n
) can be discredited as
(25)
k
i
,
j
n
=
Δ
x
[
Δ
x
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
+
Δ
y
[
Δ
y
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
.Then, (20) is implemented as follows:
(26)
ϕ
i
,
j
n
+
1
-
ϕ
i
,
j
n
Δ
t
=
λ
δ
(
ϕ
i
,
j
n
)
l
i
,
j
n
+
γ
w
(
I
,
T
*
)
i
,
j
δ
(
ϕ
i
,
j
n
)
+
μ
(
Δ
x
(
Δ
x
ϕ
i
,
j
n
)
+
Δ
y
(
Δ
y
ϕ
i
,
j
n
)
-
k
i
,
j
n
)
,
where
(27)
l
i
,
j
n
=
Δ
x
[
1
1
+
c
2
|
∇
ν
I
i
,
j
|
·
Δ
x
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
+
Δ
y
[
1
1
+
c
2
|
∇
ν
I
i
,
j
|
·
Δ
y
ϕ
i
,
j
n
(
Δ
x
ϕ
i
,
j
n
)
2
+
(
Δ
y
ϕ
i
,
j
n
)
2
]
.


## 4.2. Experimental Results

In this subsection, we present the experimental results of our model on a variety of synthetic and real images and compare them with the MCY model [15] and RSF model [13]. Unless otherwise specified, the level set function ϕ
(
x
,
y
,
t
) is simply initialized to a constant function ϕ
0
(
x
,
y
)
=
0 for all the experiments. We use the same parameters of time step Δ
t
=
5, λ
=
15, and γ
=
30 for all the images in this paper.The first experiment is to perform the proposed mode on synthetic images with increased strength of noise shown in Figure2. This group contains four images: original image and three noisy images generated by adding Gaussian white noise with variances 0.02, 0.1, and 0.2. We use κ
=
0.8 for these four images. The results of our method are shown in the second row of Figure 2. From this, we can see that the proposed model can finish the task of segment well for not only original image but also degraded images with strong noise. In addition, we compare our model with the MCY model to test the sensitivity to initialization. The experimental results for ϕ
0
=
-
1,0
,
1 are shown in Figure 3. We observe that the proposed model successfully extracts the object for ϕ
0
=
-
1,0
,
1, while the MCY model does not finish the task of segment for the case of ϕ
0
=
-
1,0. Thus, the proposed model is of more flexible initialization. From Figure 3, we also see that the robustness to noise for the proposed model is better than that for the MCY model from the third and fourth rows of Figure 3.Segmentation results by both models: the proposed model and RSF model on the synthetic image with different noise. First row: original image and noisy images (Gaussian white noise with zero mean with different variances: 0.02, 0.1, and 0.2 for (b)–(d)); second row: the proposed model; third row: RSF model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)The comparisons of the MCY model (the first and the third row) and the proposed model (the second and the fourth row) on segmenting the synthetic images in Figures2(a) and 2(b), starting with different constant functions (ϕ
0
=
-
1,0
,
1 from left column to right column).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)The second experiment is to perform the proposed model on four infrared ship images shown in Figure4. This type of image has common characteristic, that is, burry images with high strength of noise. This brings very great difficulties for image segmentation. Here, we advise that the value of κ belong to interval [
0.1,0.5
] and use κ
=
0.15,0.25,0.25,0.5 for the given image. From the first experiment, we can find that the proposed model is robust to the noise, which is used to extract the object for this kind image. The results of the MCY model and the RSF model are also shown in Figure 4. We observed that the proposed model achieves the segmentation task.Segmentation results on four ship images by three models. First row: original images; second row: the MCY model; third row: the RSF model; fourth row: the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)The third experiment is to perform the proposed model on another kind of image, skin image. To analyze skin lesions, it is necessary to accurately locate and isolate the lesions. However, it is a challenging task since the skin image contains many things, such as cutaneous melanoma, grid texture, and hair on the human skin surface. In this experiment, we test the proposed model on six skin images. The segmentation results by the proposed model and the other two models are shown in Figure5. We can see that the MCY model and the RSF model cannot deal with this kind of image, while the proposed model achieves satisfactory segmentation results for these six images.Segmentation results on skin images by three models. First column: original image; second column: the MCY model; third column: the RSF model; fourth column: the proposed model.
(a)
(b)
(c)
(d)The final experiment is to perform the proposed model on the other three real images from different modalities: a breast cyst image with complex tissue background, an aerial image with cognitive contours, and a MR image of ventriculus cordis with complex structure and weak boundary. We list the results by the proposed model, the MCY model, and the RSF model in Figure6. As we can see, the proposed model extracts the objects for these three images. They illustrate the advantages of the proposed model: the ability to deal with weak boundary and complex background.Segmentation results on other types of image by three models. First row: an aerial image, a breast cyst image, and a MR image of ventriculus cordis; second row: the MCY model; third row: the RSF model; fourth row: the proposed model.
(a)
(b)
(c)
(d)


## 5. Conclusion

A novel image segmentation model is presented in this paper, which is established by two steps. Firstly, we present a novel transition region extraction method based on morphological gradient. After this, we integrate the transition region into a variational level set framework and develop a novel geometric active contour model. Since the transition region and fraction order indicator function are used, the proposed algorithm allows for more flexible initialization and diminishes the influence of noise. The proposed model has been applied to both synthetic and real images with promising results. The experimental results confirm the effectiveness of the proposed model for image segmentation.



---

*Source: 145343-2014-10-21.xml*