# Unsupervised Joint Image Denoising and Active Contour Segmentation in Multidimensional Feature Space

Authors: Qi Ge; Xiao-Yuan Jing; Fei Wu; Jingjie Yan; Hai-Bo Li
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

We describe a new method for simultaneous image denoising and level set-based active contour segmentation using multidimensional features. We consider an image to be a surface embedded in a Riemannian manifold. By defining a metric in the embedded space, which in our case includes multidimensional image features as well as a level set-based active contour model, a minimization problem in the image space can be obtained through the Polyakov action framework. The resulting minimization problem is solved with a dual algorithm for efficiency. Benefits of this new method include the fact that it is independent of any artificial “running” parameters, and experiments using both synthetic and real images show that the method is robust with respect to noise and blurry object boundaries.


---

## Body

## 1. Introduction

Unsupervised image segmentation is an important problem with many applications in science, including medical imaging. Image segmentation is a postprocessing problem in many computer vision tasks; its aim is to divide an image into finite number of subregions. The features of different subregions are utilized as the segmentation criteria. The statistical methods, such as expectation-maximization (EM) algorithm [1] and fuzzy C-means clustering (FCM) algorithm [2], are applied in classifying the pixels based on some particular image features segmentation criteria. In general, the statistical methods achieve the classification based on only one segmentation criterion. However, there is various kinds of features in an image and the features may vary spatially. Therefore it will be not precise to use one kind of these methods. How to extract the features of an image and how to utilize these features as the segmentation criterion are significant for segmentation.Many works utilize the difference between invariable pixel intensities, as well as their spatial connectivity, in assessing whether two pixels belong to the same object. These active contour models based on the level set method [3] classify the pixels by only one image feature, that is, the image intensity based on uniform distribution [4–6]. Nevertheless, the image intensity varies spatially; thus the image intensity is not necessarily described by one kind of specific distribution. For improving the precision, the works of [7, 8] extract the multifeature to deal with more complex information content. Simultaneously, the additional artificial parameters are introduced; thus it needs the experience to set the parameters.The Polyakov action was introduced in image processing by Sochen et al. in [9]. This segmentation model is different from the other segmentation methods in two ways. First, images are represented as Riemannian manifolds embedded in a higher dimensional spatial-feature manifold. Second, the Polyakov action provides an efficient mathematical framework to embed the multifeature of images in higher-dimensional Riemannian manifolds by harmonic maps. Bresson et al. [10] propose active contour models based on the Polyakov action. These models map several kinds of features, for example, color and texture, into higher dimensional space. Because these models choose a metric with artificial parameters on the feature space, it requires careful manual parameter-tuning.In this paper, the proposed active contour model is formulated in the framework of the Polyakov action [9]. Unlike the other related works [7–9], a metric on the feature space manifold is defined by the invariant geometry of images. Consequently, the proposed method is purely based on the geometrical features of images without any artificial parameters. We implement the segmentation through two steps. First, an approximated image, removing the noise while preserving the main structures, is found in the feature space built on geometrical features of the original image. Second, the active contour is embedded into the feature space built on both the statistical and geometrical features of the approximated image. For efficiency, we solve the proposed model via the improved Chambolle dual formulation [10] of the minimization problem.The paper is organized as follows. In Section2, we introduce the mathematical framework based on the Polyakov action. In Section 3, we introduce the proposed model and the numerical algorithm of the proposed method is also summarized. In Section 4, we validate our model by some experiments on medical images. In Section 5, we end the paper by a brief conclusion.


## 2. Geometrical Framework Based on Weighted Polyakov Action

Sochen et al. introduce a general geometrical framework for low-level vision, based on the Polyakov action [9]. In this framework, images are represented as the surfaces on a Riemannian manifold. The Polyakov action is a functional that measures the weight of a mapping X
=
X
1
σ
,
…
,
X
m
σ between an n-dimensional embedded manifold (e.g., the image manifold) Σ with coordinates σ
=
σ
1
,
…
,
σ
n and the m-dimensional manifold M with the coordinates X
1
σ
,
…
,
X
m
σ, m
>
n. A Riemannian structure metric g
u
v can be introduced to measure the local distances on the embedded manifold Σ, whereas we use the metric h
i
j to measure the distance on the manifold M. To measure the weight of the mapping X
:
Σ
↦
M, the Polyakov action is used as a generalization of the L
2-norm on the embedded image to space feature manifold M:(1)
S
X
i
,
g
u
v
,
h
i
j
=
∫
g
g
μ
υ
∂
μ
X
i
∂
υ
X
j
h
i
j
d
n
σ
,where g is the determinant of the image metric tensor g
u
v and g
μ
υ is its inverse. The metric g is chosen as the induced metric, obtained by the pullback relation: g
u
v
=
h
i
j
∂
μ
X
i
∂
υ
X
j; the Polyakov energy is shortened to(2)
S
X
i
,
h
i
j
=
∫
g
d
n
σ
.In the relevant works [7, 8], the authors get the denoised image and the segmentation results by minimizing the energy functional (2) with respect to denoising and segmentation, respectively. In seminal work [9], they embed grey images in the feature x
,
y
,
I
x
,
y, where I
x
,
y is the grey intensity value for pixel x
,
y. They choose a metric h
i
j
=
d
i
a
g
1,1
,
β
2; β
>
0 is a constant. Based on this metric on feature space and the Polyakov energy, the regularization term on the intensity values is given by ∫
1
+
β
2
∇
I
2. Although it allows setting the scale of the feature dimension independently of the spatial dimensions, the accuracy of the scale is subject to the artificial parameter β.


## 3. The Active Contour Model in Multifeature Space

In this work, we utilize an improved geometrical framework based on the weighted Polyakov action without any artificial parameter. First, we get an approximated image by embedding it into the feature space constituted by the features of the original image. Second, given the approximated image, active contour is driven by embedding the level set function into the higher dimensional feature space composed of the geometrical and statistical features of the approximated image.


### 3.1. Approximating Image under an Improved Geometrical Framework

The original imageI
x
,
y is defined on the image manifold Σ with coordinates x
,
y. The approximated image u is defined on the image manifold Σ and denoted by u
x
,
y. To preserve the main edges of the original image, we extract the geometrical features of edges, ∇
S
I
=
I
x
3
+
I
x
I
y
2
,
I
y
3
+
I
y
I
x
2, derived from the anisotropic diffusion equation [2]. Considering the intensity value I
x
,
y as another feature, we build the feature space, x
,
y
,
u
,
∇
S
I
x
,
y
-
∇
S
u
x
,
y
,
I
x
,
y
-
u
x
,
y, denoted by x
,
y
,
u
,
f
1
,
f
2 for the sake of simplicity. To avoid the influence of the artificial parameter, we choose a metric tensor h
i
j on the feature space M, which is defined by the invariant geometry of the original image I. Consider h
i
j
=
d
i
a
g
1,1
,
1
,
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
,
1
/
1
+
I
x
x
+
I
y
y
+
2
I
x
y. The pullback relation yields the determinant of metric tensor g
μ
υ on manifold Σ:(3)
g
=
1
+
∇
u
2
+
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
2
∇
S
I
x
,
y
-
∇
S
u
x
,
y
2
+
1
1
+
I
x
x
+
I
y
y
+
2
I
x
y
2
I
x
,
y
-
u
x
,
y
2
.Analogizing based on the Polyakov energy (2), we get the approximated image u by minimizing the energy functional as follows:(4)
E
1
=
∫
Σ
∇
u
+
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
S
I
x
,
y
-
∇
S
u
x
,
y
d
x
d
y
+
∫
Σ
1
1
+
I
x
x
+
I
y
y
+
2
I
x
y
I
x
,
y
-
u
x
,
y
d
x
d
y
,where the weight coefficient of second term, corresponding to the third element of the metric h
i
j, denotes the coefficients of first fundamental form in differential geometry. When this weight coefficient is larger, the edge structure is enhanced in the vicinity of the edges; otherwise, smoothing the image is strengthened. The weight coefficient of the third term, corresponding to the last element of the metric h
i
j, denotes the coefficients of second fundamental form in differential geometry. Approximating the intensity I
x
,
y is strengthened when this coefficient is larger, whereas smoothing the image is strengthened when the weight is smaller.


### 3.2. Active Contour Evolution under the Improved Geometrical Framework

The active contour is represented as the zero level set functionϕ
x
,
y
=
0 on the image manifold Σ. For avoiding the effects of the intensity nonuniformity, we extract the statistical features c
1
,
c
2 on the local region of size 3
×
3, where c
1, c
2 denote the mean intensity in the local region inside and outside the zero level set. The feature space is x
,
y
,
ϕ
x
,
y
,
u
x
,
y
-
c
1
2
-
u
x
,
y
-
c
2
2, denoted by x
,
y
,
ϕ
x
,
y
,
f
3. The metric tensor defined on this feature space is h
~
i
j
=
diag
⁡
1,1
,
1
/
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
,
ϕ. The pullback relation yields the determinant of metric tensor g
~
μ
υ on manifold Σ:(5)
g
~
=
1
+
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
2
∇
ϕ
2
+
u
x
,
y
-
c
1
2
-
u
x
,
y
-
c
2
2
2
ϕ
2
.According to the Polyakov energy (2), we drive the curve evolution by minimizing the energy functional as follows:(6)
E
2
=
∫
Σ
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
ϕ
+
∫
Σ
u
x
,
y
-
c
1
2
-
u
x
,
y
-
c
2
2
ϕ
,where the weight of first term is actually an edge detector. The curve evolution tends to stop when it decreases to zero, whereas the evolution goes on.


### 3.3. Dual Algorithm

To apply the dual gradient algorithm, we introduce the dual variable,p. The total variation term in (4) and (6) can be formulated as follows:(7)
∫
Ω
∇
u
d
x
=
m
a
x
p
∈
C
1
u
,
div
⁡
p
1
,
C
1
≔
p
1
∣
p
1
∈
C
c
1
Ω
,
R
2
,
p
1
≤
1
,
∀
x
∈
Ω
,
∫
Ω
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
ϕ
d
x
=
m
a
x
p
∈
C
g
ϕ
,
div
⁡
p
2
,
C
2
≔
p
2
∣
p
2
∈
C
c
1
Ω
,
R
2
,
p
2
≤
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
,
∀
x
∈
Ω
.The approximation formulation of the energy of our model can be rewritten as(8)
m
i
n
0
≤
ϕ
≤
1
,
u
m
a
x
p
1
∈
C
1
p
2
∈
C
2
E
1
+
E
2
=
ϕ
,
div
⁡
p
2
+
u
,
div
⁡
p
1
+
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
S
I
-
∇
S
u
L
1
+
1
1
+
I
x
x
+
I
y
y
+
2
I
x
y
I
-
u
L
1
+
ϕ
,
R
u
,
c
1
,
c
2
,where R
c
1
,
c
2
,
u
=
u
-
c
1
2
-
u
-
c
2
2. We then apply the split Chambolle dual algorithm [10] to solve the optimization problem.Introducing the auxiliary variablesv
1, v
2, solving the energy functional (8) is equivalent to minimizing the problem as follows:(9)
m
i
n
u
,
ϕ
m
a
x
p
1
∈
C
1
p
2
∈
C
2
E
1
+
E
2
=
u
,
div
⁡
p
1
+
1
2
θ
I
-
u
-
v
1
L
2
2
+
α
1
∇
S
I
-
∇
S
u
L
1
-
v
1
L
2
2
+
α
2
v
1
L
1
+
ϕ
,
div
⁡
p
2
+
1
2
θ
ϕ
-
v
2
L
2
2
+
R
c
1
,
c
2
,
u
,
v
2
+
γ
η
v
2
,where the parameter θ
>
0 is chosen to be small for avoiding smearing the edges (in this paper, we choose θ
=
0.15.), v
z
=
m
a
x
0
,
2
z
-
1
/
2
-
1 is an exact penalty function provided that the constant γ is chosen large enough compared to λ such as γ
>
λ
/
2
R
L
∞
Ω, α
1
=
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y, and α
2
=
1
/
1
+
I
x
x
+
I
y
y
+
2
I
x
y. The minimization problem (9) can be divided into four subproblems as follows and can be solved alternatively.(a) Given imageI, v
1, update u. we search for u as the solution of (10)
m
i
n
u
u
,
div
⁡
p
1
+
1
2
θ
I
-
u
-
v
1
L
2
2
.The solution of (10) is given by(11)
u
=
I
-
v
1
-
θ
div
⁡
p
1
,where p
1
=
p
1
,
p
2 can be updated by fixed point method: initializing p
1
=
0 and updating (12)
p
1
k
+
1
=
p
1
k
+
δ
∇
div
⁡
p
1
k
-
I
-
v
1
/
θ
1
+
δ
∇
div
⁡
p
1
k
-
I
-
v
1
/
θ
.In this paper, we choose δ
≤
1
/
8 to ensure convergence.(b) Givenu, v
2, we search for ϕ by solving the minimization problem as follows:(13)
m
i
n
ϕ
,
v
2
ϕ
,
div
⁡
p
2
+
1
2
θ
ϕ
-
v
2
L
2
2
.The solution of (13) is given by(14)
ϕ
=
v
2
-
θ
div
⁡
p
2
.Equation (14) is solved by a fixed point method:(15)
p
2
k
+
1
=
p
2
k
+
δ
∇
div
⁡
p
2
k
-
v
2
k
/
θ
1
+
δ
∇
div
⁡
p
2
k
-
v
2
k
/
θ
,
k
=
0,1
,
2
,
…
.(c) Given the solution ofϕ, we search v
2 by solving(16)
m
i
n
ϕ
,
v
2
1
2
θ
ϕ
-
v
2
L
2
2
+
R
c
1
,
c
2
,
u
,
v
2
+
γ
υ
v
2
.The solution of (16) is given by(17)
v
2
=
min
⁡
max
⁡
ϕ
-
θ
R
c
1
,
c
2
,
u
,
0
,
1
.(d) GivenI, u, we search for v
1 as the solution of(18)
m
i
n
ϕ
,
u
1
2
θ
I
-
u
-
v
1
L
2
2
+
α
1
∇
S
I
-
∇
S
u
L
1
-
v
1
L
2
2
+
α
2
v
1
L
1
.The solution of (18) is given by(19)
v
1
=
I
-
u
+
2
α
1
θ
∇
S
I
-
∇
S
u
L
1
-
θ
α
1
,
i
f
I
-
u
≥
θ
α
1
,
I
-
u
+
2
α
1
θ
∇
S
I
-
∇
S
u
L
1
-
θ
α
1
,
i
f
I
-
u
≤
θ
α
1
,
0
,
i
f
I
-
u
≥
θ
α
1
.After v
1 is solved, it is utilized in (12) for next iteration.The algorithm of minimizing our model is described in the following.Step
1. Initialize u
0
=
ϕ
0
=
v
1
0
=
v
2
0
=
p
1
0
=
p
2
0
=
0.Step
2. Given the fixed threshold of iterations T
i
t
e
r
a
t
i
o
n
>
0, if k
=
T
i
t
e
r
a
t
i
o
n, then stop; else go to Step  3.Step
3. Do the iteration for solving subproblem.Step
3.1. While u
k
-
u
k
-
1
>
ε doStep
3.2. Update the dual variable p
1
k
→
p
1
k
+
1 by (14) and then update u
k
→
u
k
+
1 by (13); go to Step  3.3.Step
3.3. Given v
2
k, update the dual variable p
2
k
→
p
2
k
+
1 according to (15) and then update ϕ
k
→
ϕ
k
+
1 by (16); go to Step  3.4.Step
3.4. Given ϕ
k
+
1, u
k
+
1, compute R
c
1
,
c
2
,
u and update v
2
k
→
v
2
k
+
1 by (17). Go to Step  3.5.Step
3.5. Given I, u, update v
1
k
→
v
1
k
+
1 by (19); then go to Step  3.1; otherwise, go to Step  4.Step
4. End while.


## 3.1. Approximating Image under an Improved Geometrical Framework

The original imageI
x
,
y is defined on the image manifold Σ with coordinates x
,
y. The approximated image u is defined on the image manifold Σ and denoted by u
x
,
y. To preserve the main edges of the original image, we extract the geometrical features of edges, ∇
S
I
=
I
x
3
+
I
x
I
y
2
,
I
y
3
+
I
y
I
x
2, derived from the anisotropic diffusion equation [2]. Considering the intensity value I
x
,
y as another feature, we build the feature space, x
,
y
,
u
,
∇
S
I
x
,
y
-
∇
S
u
x
,
y
,
I
x
,
y
-
u
x
,
y, denoted by x
,
y
,
u
,
f
1
,
f
2 for the sake of simplicity. To avoid the influence of the artificial parameter, we choose a metric tensor h
i
j on the feature space M, which is defined by the invariant geometry of the original image I. Consider h
i
j
=
d
i
a
g
1,1
,
1
,
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
,
1
/
1
+
I
x
x
+
I
y
y
+
2
I
x
y. The pullback relation yields the determinant of metric tensor g
μ
υ on manifold Σ:(3)
g
=
1
+
∇
u
2
+
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
2
∇
S
I
x
,
y
-
∇
S
u
x
,
y
2
+
1
1
+
I
x
x
+
I
y
y
+
2
I
x
y
2
I
x
,
y
-
u
x
,
y
2
.Analogizing based on the Polyakov energy (2), we get the approximated image u by minimizing the energy functional as follows:(4)
E
1
=
∫
Σ
∇
u
+
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
S
I
x
,
y
-
∇
S
u
x
,
y
d
x
d
y
+
∫
Σ
1
1
+
I
x
x
+
I
y
y
+
2
I
x
y
I
x
,
y
-
u
x
,
y
d
x
d
y
,where the weight coefficient of second term, corresponding to the third element of the metric h
i
j, denotes the coefficients of first fundamental form in differential geometry. When this weight coefficient is larger, the edge structure is enhanced in the vicinity of the edges; otherwise, smoothing the image is strengthened. The weight coefficient of the third term, corresponding to the last element of the metric h
i
j, denotes the coefficients of second fundamental form in differential geometry. Approximating the intensity I
x
,
y is strengthened when this coefficient is larger, whereas smoothing the image is strengthened when the weight is smaller.


## 3.2. Active Contour Evolution under the Improved Geometrical Framework

The active contour is represented as the zero level set functionϕ
x
,
y
=
0 on the image manifold Σ. For avoiding the effects of the intensity nonuniformity, we extract the statistical features c
1
,
c
2 on the local region of size 3
×
3, where c
1, c
2 denote the mean intensity in the local region inside and outside the zero level set. The feature space is x
,
y
,
ϕ
x
,
y
,
u
x
,
y
-
c
1
2
-
u
x
,
y
-
c
2
2, denoted by x
,
y
,
ϕ
x
,
y
,
f
3. The metric tensor defined on this feature space is h
~
i
j
=
diag
⁡
1,1
,
1
/
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
,
ϕ. The pullback relation yields the determinant of metric tensor g
~
μ
υ on manifold Σ:(5)
g
~
=
1
+
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
2
∇
ϕ
2
+
u
x
,
y
-
c
1
2
-
u
x
,
y
-
c
2
2
2
ϕ
2
.According to the Polyakov energy (2), we drive the curve evolution by minimizing the energy functional as follows:(6)
E
2
=
∫
Σ
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
ϕ
+
∫
Σ
u
x
,
y
-
c
1
2
-
u
x
,
y
-
c
2
2
ϕ
,where the weight of first term is actually an edge detector. The curve evolution tends to stop when it decreases to zero, whereas the evolution goes on.


## 3.3. Dual Algorithm

To apply the dual gradient algorithm, we introduce the dual variable,p. The total variation term in (4) and (6) can be formulated as follows:(7)
∫
Ω
∇
u
d
x
=
m
a
x
p
∈
C
1
u
,
div
⁡
p
1
,
C
1
≔
p
1
∣
p
1
∈
C
c
1
Ω
,
R
2
,
p
1
≤
1
,
∀
x
∈
Ω
,
∫
Ω
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
ϕ
d
x
=
m
a
x
p
∈
C
g
ϕ
,
div
⁡
p
2
,
C
2
≔
p
2
∣
p
2
∈
C
c
1
Ω
,
R
2
,
p
2
≤
1
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
,
∀
x
∈
Ω
.The approximation formulation of the energy of our model can be rewritten as(8)
m
i
n
0
≤
ϕ
≤
1
,
u
m
a
x
p
1
∈
C
1
p
2
∈
C
2
E
1
+
E
2
=
ϕ
,
div
⁡
p
2
+
u
,
div
⁡
p
1
+
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y
∇
S
I
-
∇
S
u
L
1
+
1
1
+
I
x
x
+
I
y
y
+
2
I
x
y
I
-
u
L
1
+
ϕ
,
R
u
,
c
1
,
c
2
,where R
c
1
,
c
2
,
u
=
u
-
c
1
2
-
u
-
c
2
2. We then apply the split Chambolle dual algorithm [10] to solve the optimization problem.Introducing the auxiliary variablesv
1, v
2, solving the energy functional (8) is equivalent to minimizing the problem as follows:(9)
m
i
n
u
,
ϕ
m
a
x
p
1
∈
C
1
p
2
∈
C
2
E
1
+
E
2
=
u
,
div
⁡
p
1
+
1
2
θ
I
-
u
-
v
1
L
2
2
+
α
1
∇
S
I
-
∇
S
u
L
1
-
v
1
L
2
2
+
α
2
v
1
L
1
+
ϕ
,
div
⁡
p
2
+
1
2
θ
ϕ
-
v
2
L
2
2
+
R
c
1
,
c
2
,
u
,
v
2
+
γ
η
v
2
,where the parameter θ
>
0 is chosen to be small for avoiding smearing the edges (in this paper, we choose θ
=
0.15.), v
z
=
m
a
x
0
,
2
z
-
1
/
2
-
1 is an exact penalty function provided that the constant γ is chosen large enough compared to λ such as γ
>
λ
/
2
R
L
∞
Ω, α
1
=
1
+
I
x
2
+
I
y
2
+
2
I
x
I
y, and α
2
=
1
/
1
+
I
x
x
+
I
y
y
+
2
I
x
y. The minimization problem (9) can be divided into four subproblems as follows and can be solved alternatively.(a) Given imageI, v
1, update u. we search for u as the solution of (10)
m
i
n
u
u
,
div
⁡
p
1
+
1
2
θ
I
-
u
-
v
1
L
2
2
.The solution of (10) is given by(11)
u
=
I
-
v
1
-
θ
div
⁡
p
1
,where p
1
=
p
1
,
p
2 can be updated by fixed point method: initializing p
1
=
0 and updating (12)
p
1
k
+
1
=
p
1
k
+
δ
∇
div
⁡
p
1
k
-
I
-
v
1
/
θ
1
+
δ
∇
div
⁡
p
1
k
-
I
-
v
1
/
θ
.In this paper, we choose δ
≤
1
/
8 to ensure convergence.(b) Givenu, v
2, we search for ϕ by solving the minimization problem as follows:(13)
m
i
n
ϕ
,
v
2
ϕ
,
div
⁡
p
2
+
1
2
θ
ϕ
-
v
2
L
2
2
.The solution of (13) is given by(14)
ϕ
=
v
2
-
θ
div
⁡
p
2
.Equation (14) is solved by a fixed point method:(15)
p
2
k
+
1
=
p
2
k
+
δ
∇
div
⁡
p
2
k
-
v
2
k
/
θ
1
+
δ
∇
div
⁡
p
2
k
-
v
2
k
/
θ
,
k
=
0,1
,
2
,
…
.(c) Given the solution ofϕ, we search v
2 by solving(16)
m
i
n
ϕ
,
v
2
1
2
θ
ϕ
-
v
2
L
2
2
+
R
c
1
,
c
2
,
u
,
v
2
+
γ
υ
v
2
.The solution of (16) is given by(17)
v
2
=
min
⁡
max
⁡
ϕ
-
θ
R
c
1
,
c
2
,
u
,
0
,
1
.(d) GivenI, u, we search for v
1 as the solution of(18)
m
i
n
ϕ
,
u
1
2
θ
I
-
u
-
v
1
L
2
2
+
α
1
∇
S
I
-
∇
S
u
L
1
-
v
1
L
2
2
+
α
2
v
1
L
1
.The solution of (18) is given by(19)
v
1
=
I
-
u
+
2
α
1
θ
∇
S
I
-
∇
S
u
L
1
-
θ
α
1
,
i
f
I
-
u
≥
θ
α
1
,
I
-
u
+
2
α
1
θ
∇
S
I
-
∇
S
u
L
1
-
θ
α
1
,
i
f
I
-
u
≤
θ
α
1
,
0
,
i
f
I
-
u
≥
θ
α
1
.After v
1 is solved, it is utilized in (12) for next iteration.The algorithm of minimizing our model is described in the following.Step
1. Initialize u
0
=
ϕ
0
=
v
1
0
=
v
2
0
=
p
1
0
=
p
2
0
=
0.Step
2. Given the fixed threshold of iterations T
i
t
e
r
a
t
i
o
n
>
0, if k
=
T
i
t
e
r
a
t
i
o
n, then stop; else go to Step  3.Step
3. Do the iteration for solving subproblem.Step
3.1. While u
k
-
u
k
-
1
>
ε doStep
3.2. Update the dual variable p
1
k
→
p
1
k
+
1 by (14) and then update u
k
→
u
k
+
1 by (13); go to Step  3.3.Step
3.3. Given v
2
k, update the dual variable p
2
k
→
p
2
k
+
1 according to (15) and then update ϕ
k
→
ϕ
k
+
1 by (16); go to Step  3.4.Step
3.4. Given ϕ
k
+
1, u
k
+
1, compute R
c
1
,
c
2
,
u and update v
2
k
→
v
2
k
+
1 by (17). Go to Step  3.5.Step
3.5. Given I, u, update v
1
k
→
v
1
k
+
1 by (19); then go to Step  3.1; otherwise, go to Step  4.Step
4. End while.


## 4. Experimental Results

All the experiments are run with Matlab code on the PC of CPU 3.2 GHz, RAM 728 M. we show the experiments results for medical image segmentation of Chan-Vese model (CV) [11], the structure-based level set method (SLM) [12], and the region-scale fitting model (RSF) [4]. Figure 1 shows the experiments on the synthesized noisy images. This image is of size 266
×
313 with 10% white Gaussian noise.Figure 1
Segmentation for synthetic noisy image.
(a)
Original image
(b)
CV model 200th
(c)
SLM 50th
(d)
The proposed model 50th
(e)
The approximated imageAs shown in Figures1(b) and 1(c), the CV model and the SLM models are sensitive to noise. As shown in Figure 1(d), the proposed model is robust to noise and obtain the correct boundary. The CV model generates the unwanted contours because of the strong noise. Based on the edge detector function, the SLM model is more robust to the noise than the CV model. The result of the proposed model shows that it is able to extract the real object boundary even when the noise is strong.Figure2 is brain magnetic resonance image (MRI) of size 397
×
397 with 2% noise and 10% level intensity nonuniformity. The brain MRI mainly consists of three parts: the cerebrospinal fluid, the gray matter, and the white matter. The cerebrospinal fluid is the dark matter which exists in two places: the middle of the brain surrounded by the gray matter and the gap between the cranium and the brain. The task of segmenting the brain MRI is to extract the contour profile between the white matter and the gray matter. Since the contrast of boundary between the white matter and the gray matter is lower than the boundary between the cerebrospinal fluid and the gray matter, the latter is always extracted wrongly as the object boundary. As shown in Figures 2(b) and 2(c), the CV model and the SLM model cannot extract the object boundary with low-contrast. Although the RSF model can extract the boundaries with low-contrast, it also extracts some unwanted objects. The segmentation results of the proposed model show that it is robust to noise and can extract more boundaries with low-contrast.Figure 2
Segmentation for brain MRI. (a) Original image and initial level set contour. (b) Segmentation result of the CV model with 200 iterations. (c) Segmentation result of the SLM with 100 iterations. (d) Segmentation result of the RSF model with 200 iterations. (e) Segmentation result of the proposed model with 30 iterations. (f) The approximated result of the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)Figure3 is a brain MRI of size 258
×
258 with 5% noise and 40% intensity nonuniformity. As shown in Figures 3(b) and 3(c), the drawbacks of the RSF model and the SLM model still exist. Compared with the other methods, the proposed method clearly extracts more object boundaries with low-contrast between the gray matter and the white matter.Figure 3
Segmentation for brain MRI. (a) Original image and initial level set contour. (b) Segmentation result of the CV model with 200 iterations. (c) Segmentation result of the SLM with 100 iterations. (d) Segmentation result of the RSF model with 200 iterations. (e) Segmentation result of the proposed model with 30 iterations. (f) The approximated result of the proposed model.
(a)
(b)
(c)
(d)
(e)
(f)Figure4 show the segmentation results of the active contour methods based on multilayer level set functions. By the multilayer level set functions, the cerebrospinal fluid, the gray matter, and the white matter can be extracted simultaneously. In Figures 4(b) and 4(f), the results of CV model show that the cerebrospinal fluid is not extracted completely. As shown in Figures 4(c) and 4(g), the boundaries between the white and gray matters are not extracted completely by SLM model. We can observe that in Figures 4(d) and 4(h), the results of the proposed method based on multilayer level set functions show that the completed boundaries of the cerebrospinal fluid, the gray matter, and the white matter are extracted simultaneously. To show the convergence speed of the compared methods, Table 1 shows the iteration numbers and processing time in each iteration for the CV model, the SLM model based on the steepest descent method, and the proposed model with both single- and multilayer level set functions. Table 2 shows the segmentation accuracy of the compared active contour models based on multilayer level set functions.Table 1
Iteration numbers and processing time in each iteration.
CV model
SLM model
The proposed model
Single level set function
0.35 sec
0.23 sec
0.18 sec
Multilayer functions
1.47 sec
0.62 sec
0.39 sec
Iteration number for convergence
100 iterations
50 iterations
30 iterationsTable 2
Quantitative evaluation for brain MR-data in Figure4.
The CV model
The SLM model
The proposed model
Overall accuracy
65.2%
34.6%
91.5%
GM (dice metric)
78.3%
82.4%
97.2%
WM (dice metric)
69.0%
13.5%
94.7%Figure 4
Segmentation for brain MRIs by the active contour methods based on multilayer level set functions. (a) and (e) are original images and initial level set contours. (b) and (f) are segmentation results of the CV model with 100 iterations. (c) and (g) are segmentation results of the SLM with 50 iterations. (d) and (h) are the results of the proposed method with 30 iterations.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)The data in Figure5 is download from the website [13]. We also show the segmentation accuracy in Table 3 by the DICE metric [14], compared with the ground truth given in this website. We can see from Figures 5(c) and 5(d) that the CV model and the SLM model only extract the boundaries with high contrast, while the object boundaries between the gray matter and white matter are not extracted. And some CSF of the image is not extracted. The proposed method can extract the more completed CSF and can preserve the cerebral cortex more efficiently.Table 3
Quantitative evaluation for brain MR-data in Figure5.
The CV model
The SLM model
The RSF model
The proposed model
Overall accuracy
15.4%
15.7%
82.6%
88.7%
GM (dice metric)
13.5%
14.1%
79.6%
90.4%
WM (dice metric)
16.2%
16.8%
83.7%
87.9%Figure 5
Segmentation for 3D brain image: (a) original image; (b) segmentation result of the CV model with 200 iterations; (c) segmentation result of the SLM model with 100 iterations; (d) segmentation result of the RSF model with 200 iterations; (e) segmentation result of the proposed model with 30 iterations; (f) the denoising result of the proposed model; (g) the 3D original image; (h) the 3D segmentation result of the proposed method; (i) the approximated 3D image.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)


## 5. Conclusion

In this paper, we propose a new variational model for image segmentation and image denoising simultaneously. We obtain the approximated image by embedding the approximating criteria into a specific multifeature space. And then the segmentation result is obtained by embedding the active contour into another multifeature space which is composed by the segmentation criteria depending on the approximated image. The segmentation and the denoising problems are solved by the split Chambolle dual algorithm alternately. The comparisons of the other popular segmentation models demonstrate the accuracy and efficiency of the proposed model.



---

*Source: 3909645-2016-07-20.xml*