Sage Journals: Discover world-class research

Abstract

Regularization is a key element of variational models in image processing. To overcome the weakness of models based on total variation, various high order (typically second order) regularization models have been proposed and studied recently. Among these, Euler's elastica energy based regularizer is perhaps the most interesting in terms of both mathematical and physical justifications. More importantly its success has been proven in applications; however it has been a major challenge to develop fast and effective algorithms.

In this paper we propose a new idea for deriving a primal dual algorithm, based on Legendre–Fenchel transformations, for representing the elastica regularizer. Combined with an augmented Lagrangian for-mulation, we are able to derive an equivalent unconstrained optimization that has fewer variables to work with than previous works based on splitting methods. We shall present our algorithms for both the image restoration problem and the image segmentation model. The idea applies to other models where the elastica regularizer is required. Numerical experiments show that the proposed method can produce highly competitive results with better efficiency.

Keywords

Euler’s elastica Legendre–Fenchel transformation primal-dual technique augmented Lagrangian method image denoising image segmentation

Introduction

We are primarily concerned with the challenging problem of developing effective algorithms for the Euler elastica regularizer. The latter is often used in imaging applications including image restoration, which is briefly reviewed below to motivate the new work.

Image restoration is about acquiring and recovering unknown u = u(x) (without any restrictions) from an observed image by z = z(x), x ∈ Ω ⊂ ℝ ^d where Ω is the bounded domain and has a Lipschitz boundary. Here we consider planar images, i.e. d = 2 and x = (x₁, x₂). All results and discussions will be applicable to d > 2 and other models. An additive image restoration model assumes $z = A u + η 0$ with η₀ representing some unknown Gaussian noise of mean zero and deviation σ, and $A$ a blurring operator.^1–6 The additive image restoration model minimises the fidelity to z and leads to the least-square problem

min_{u} \int Ω | A u - z | 2 d x

(1)

according to the maximum likelihood principle.⁷ Here

A = I

for image denoising. Problem (1) is in general ill-posed due to non-uniqueness, therefore, how to effectively solve it becomes a fundamental task.

The classical regularization technique by Tikhonov et. al⁸ is to add a smoothing regularization term into the energy functional to derive the following minimization problem. The resulting well-posed model admits an unique solution. This classical model cannot preserve image edges, though it is simple to use. Based on the total variation (TV) regularizer, the Rudin-Osher-Fatemi (ROF)⁴ model preserves the image edges by seeking solutions of piecewise constant functions in the space of bounded variation functions (BV) and, hence, is widely used. A variety of methods based on TV regularization have been developed to deal with imaging problems such as image restoration,^9–12 image registration,^13–15 image decomposition,^3,16,17 image inpainting^18–21 and image segmentation.^22, ²³

One disadvantage of the ROF model is that it yields so-called blocky (staircase) effects, in restoring smooth images in applications where edges are not the main features.^24–27 Another disadvantage of the model is the loss of image contrast.²⁸ It should be remarked that the recently popular method by the iterative regularization technique²⁹ can reduce the staircasing effect and improve on the image contrast to some extent; it also provides fast implementation.

The elastica was discovered by Euler in 1744.³⁰ Euler’s elastica energy function defined by

TeeV (u):= \int Ω [a + | \nabla \cdot n | 2] | \nabla u | d x

(2)

was firstly used for computer vision by Mumford,³¹ where

n = \nabla u | \nabla u |

. In some applications, often, a valid a priori assumption is that level lines are mostly smooth in addition to small curvature, which suggests the use of curvature-penalizing priors, such as in image reconstruction, inpainting,³² denoising³³ or segmentation.³⁴ Hence the curvature-dependent functionals have gained importance during the past decade. However, few effective methods exist for such variational problems.

The rest of this paper is organized as follows. In the next section, we briefly review two recent models usingl Euler’s elastica for image processing and their associated numerical algorithms. In the section “The proposed algorithm for Eulers elastica regularization” we present our new algorithm for the models reviewed in the previous section. These are: image restoration and segmentation, respectively. The next section provides some numerical results to illustrate the effectiveness of our new algorithm. The final section draws some conclusions.

Previous works based on Euler’s elastica

As a regularizer, Euler’s elastica can be found in and potentially applied to many fields of imaging models. Below we review two models that interest us most.

Euler’s elastica denoising model³³

Shen et al.³² proposed to interpolate a gray-valued image by extending its lines of constant intensity (isophotes) in the inpainting domain. Their model for image denoising takes the form

min u (α \int Ω [a + b | κ (u) | 2] | \nabla u | d x + 12 \int Ω | u - z | 2 d x)

(3)

where the mean curvature

κ (u) = \nabla \cdot \nabla u | \nabla u |

is the key quantity for Euler’s elastica energy. This approach is in the spirit of the early works of Mumford³¹ to model the ability of a visual system to complete edges in occlusion. Equation (3) leads to a highly nonlinear fourth order Euler-Lagrange equation which does not have known efficient solvers. In Shen et al.³² a slow time-marching solver was used.

Instead of solving it directly, Tai et al.³³ developed an augmented Lagrangian method that avoids the complicated higher-order derivatives. First rewrite (3) as

min v, u, m, p, n (α \int Ω [a + b | \nabla \cdot n | 2] | p | d x + 12 \int Ω (v - z) 2 d x) s . t . v = u, p = \nabla u, n = m, | p | = m \cdot p, | m | \in [- 1, 1]

(4)

Second, the augmented Lagrangian method solves the following

L (v, u, m, p, n, β 1, β 2, β 3, β 4) = 12 \int Ω (v - z) d x + α \int Ω [a + b | \nabla \cdot n | 2] | p | d x + r 1 \int Ω (| p | - p \cdot m) 2 d x + \int Ω β 1 (| p | - p \cdot m) d x + r 2 \int Ω | p - \nabla u | 2 d x + \int Ω β 2 \cdot (p - \nabla u) d x + r 3 \int Ω (v - u) 2 d x + \int Ω β 3 (v - u) d x + r 4 \int Ω | n - m | 2 d x + \int Ω β 3 \cdot (n - m) d x + δ D (m)

(5)

Further by alternating minimization, a series of subproblems (each much simpler to solve) are formulated. Note that, instead of the single variable u in (3), the total number of variables in (5) is 14—a substantial increase in complexity.

A Euler’s elastica-based segmentation model

The widely used Chan-Vese segmentation model³⁵ minimizes the length (via TV) of detected objects. Minimizing Euler’s elastica³⁴ instead of TV leads to a new segmentation model capable of integrating missing or broken parts to form complete meaningful objects and capturing objects with tiny but elongated structures. Their Euler’s elastica-based segmentation model takes the following form

min φ (α \int Ω [a + b | \nabla \cdot \nabla H (φ) | \nabla H (φ) | | 2] | \nabla H (φ) | d x + 12 \int Ω [H (φ) | f - c 1 | 2 + (1 - H (φ)) | f - c 2 | 2] d x)

(6)

where λ, a, b are positive parameters, and φ is the level set function whose zero level curve defines the object.

This regularization was originally proposed and used in the famous work of segmentation with depth by Nitzberg et al.³⁶

Letting u = H(φ) ∈{0, 1} and relaxing it to u ∈ [0, 1], the model (6) can be written as

min u (α \int Ω [a + b | \nabla \cdot \nabla u | \nabla u | | 2] | \nabla u | d x + 12 \int Ω u | f - c 1 | 2 + (1 - u) | f - c 2 | 2 d x)

(7)

which is amenable to simplification to the following constrained minimization problem

min u, p, n, c 1, c 2 (α \int Ω [a + b | \nabla \cdot n | 2] | p | d x + 12 \int Ω u | f - c 1 | 2 + (1 - u) | f - c 2 | 2 d x) s . t . p = \nabla u, n = \nabla u | \nabla u |, u \in [0, 1] .

(8)

Then the augmented Lagrangian functional

L (v, u, p, n, m, c 1, c 2, β 1, β 2, β 3, β 4) = 12 \int Ω v | f - c 1 | 2 + (1 - v) | f - c 2 | 2 d x) + α \int Ω [a + b | \nabla \cdot n | 2] | p | d x + r 1 \int Ω (| p | - p \cdot m) 2 d x + \int Ω β 1 (| p | - p \cdot m) d x + r 2 \int Ω | p - \nabla u | 2 d x + \int Ω β 2 \cdot (p - \nabla u) d x + r 3 \int Ω (v - u) 2 d x + \int Ω β 3 (v - u) d x + δ D (v) + r 4 \int Ω | n - m | 2 d x + \int Ω β 4 \cdot (n - m) d x + δ D (m)

(9)

is proposed, where δ_D(v) and δ_D( m ) are the characteristic functions on the sets D and R = { m ∈ L²(Ω): | m | ≤ 1 a.e. in Ω}. Further each subproblem can be made simple to solve.³⁴ However one notes that the total number of variables (9) is 16, which is again quite substantial.

The proposed algorithm for Euler’s elastica regularization

We shall take a different approach to develop algorithms for the above models. Our idea is to reformulate Euler’s elastica regularization term TeeV(u) first. This is achieved by a transformation.

Definition 1 (Legendre–Fenchel transformation)

Let $S$ be a normed vector space and $F : S \to ℝ \cup {+ \infty}$ be any functional. Then

F * (φ) = sup u \in S {〈 φ, u 〉 - F (u)}

(10)

is called the Legendre–Fenchel transformation of F.

The Legendre–Fenchel transformation $F * : S * \to ℝ \cup {+ \infty}$ is convex, regardless whether or not F is convex. For vectors $y \in ℝ d$ , the well-known equality used in dual formualtions $| y | 1 = max_{| ω | 1 \leq 1} ω \cdot y$ is related to and quite different from a Legendre–Fenchel transformation for $f (x) = | x | 1 f * (s) = max_{x} {s T x - | x | 1} = δ R (s) .$ Similarly for $g (x) = | x - x 0 | 1$ , $g * (s) = max_{x} {s T x - | x - x 0 | 1} = δ R (s) + s T x 0$ . For simple functionals, we can also compute the transform explicitly, e.g. if $f (η) = 12 \int Ω η 2 (x) d x η \in L 2 (Ω)$

f * (φ) = sup η \in L 2 (Ω) 〈 φ, η 〉 - 12 \int Ω η 2 d x = 12 \int Ω φ 2 d x

(11)

as the minimizer η^*= φ is found directly.

In this work, we shall use the Legendre–Fenchel transformation to study TeeV(u) from (2).

Lemma 1

Based the Legendre–Fenchel transformation of $F (φ) = 12 \int Ω φ 2 (x) d x$ from (11), it holds that

12 \int Ω (\nabla \cdot n) 2 | \nabla u | d x = sup φ \in L 2 (Ω) \int Ω (\nabla \cdot n) φ (x) | \nabla u | d x - 12 \int Ω | φ | 2 | \nabla u | d x

Proof

Take $φ = (\nabla \cdot n) \sqrt{| \nabla u |}$ . Then (11) becomes

12 \int Ω (\nabla \cdot n) 2 | \nabla u | d x = sup φ \in L 2 (Ω) \int Ω (\nabla \cdot n) \sqrt{| \nabla u |} φ (x) d x - 12 \int Ω | φ | 2 d x = sup φ \sqrt{| \nabla u |} \in L 2 (Ω) \int Ω (\nabla \cdot n) \sqrt{| \nabla u |} φ (x) \sqrt{| \nabla u |} d x - 12 \int Ω | φ | 2 | \nabla u | d x = sup φ \in L 2 (Ω) \int Ω (\nabla \cdot n) φ (x) | \nabla u | d x - 12 \int Ω | φ | 2 | \nabla u | d x

Consequently from Lemma, we obtain a useful reformulation for (2)

which will be used below. We give two more simple lemmas to aid further discussion.

Lemma 2

For any vector $x \in ℝ d$ , the unit vector $y = x | x |$ ( y = 0 when x = 0 ) is equivalent to $P y x = x$ , where $P y = yy T$ is a projection matrix.

Lemma 3

The formula for the normal vector $n = \nabla u | \nabla u |$ is equivalent to $(nn T - I) \nabla u = 0$ .

Finally, the Euler elastica functional TeeV(u) can be split into the following scheme

{TeeV (u) = P 1 (u, n):= sup ϕ \in L 2 (Ω) \int Ω [a + b (2 (\nabla \cdot n) ϕ - ϕ 2)] (\nabla u \cdot n) dx, s . t . c (u, n):= (nn T - I) \nabla u = 0

(13)

in which all terms are differentiable.

A new algorithm for image denoising

Now we are ready to present our augmented Lagrangian-based primal dual approach based on the Legendre–Fenchel transformation, which serves as an alternative algorithm to (5) for the elastica based model (3).

The approach considered here differs from existing augmented Lagrangian approaches for the solution of the same problem; indeed, the augmented Lagrangian functional we use here contains only one Lagrange multiplier (instead of three), and three associated augmentation (dual) terms (instead of four).

Following equation (13), the denoising problem (3) takes the equivalent form We next reformulate the constrained minimization in (14) by the method of augmented Lagrangian multipliers. This leads to the unconstrained optimization problem with the multiplier update

q k + 1 = q k + γc (u, n)

(15)

For the sake of effective computation, we use the Legendre–Fenchel transformation (L²-transformation) as it was done in

P 1 (u, n)

-energy in the same way to derive

min (u, n) \in W {α P 1 (u, n) + γ sup β \in L 2 (Ω) (\int Ω 〈 c (u, n) + q k γ, β 〉 - 12 | β | 2 d x) + 12 \int Ω | u - z | 2 d x}

(16)

where the notation 〈·〉 is for vector operations.

Our goal is to develop a fast algorithm to seek the saddle-point of the energy functional

L (u, n, ϕ, β):= α \int Ω [a + b (2 〈 \nabla \cdot n, ϕ 〉 - | ϕ | 2)] 〈 \nabla u, n 〉 d x + γ (\int Ω 〈 c (u, n) + q k γ, β 〉 - 12 | β | 2 d x) + 12 \int Ω | u - z | 2 d x

with multiplier update (15). Since a saddle point of the the above augmented Lagrangian functional relates to a solution of optimization problem (14), we propose to use a majoration minimization/maximization (MM) approximation to define the relaxed functional associated with (15) as follows

L k (u, n, ϕ, β) : = α \int Ω [a + b (2 〈 \nabla \cdot n, ϕ 〉 - | ϕ | 2)] 〈 \nabla u, n 〉 d x + 12 \int Ω | u - z | 2 d x + γ (\int Ω 〈 c (u, n) + q k γ, β 〉 - 12 | β | 2 d x) + τ 2 \int Ω | u - u k | 2 d x + τ 2 \int Ω | n - n k | 2 d x - τ 2 \int Ω | ϕ - ϕ k | 2 d x - τ 2 \int Ω | β - β k | 2 d x

We minimize/maximize the corresponding functional for each of the variable functions u, n , ϕ, β by fixing the other ones, and after the advance of all the variable functions, we update the Lagrange multipliers accordingly based on the optimization theory. The procedure will be repeated until the variable functions are all convergent, which means a saddle point of the augmented Lagrangian functional is obtained.

The corresponding subproblems are given by

The remainder of this subsection describes in detail each step of the process.

u-subproblem

In u-subproblem, the energy functional

α \int Ω h a, b (n, ϕ) 〈 \nabla u, n 〉 d x + γ \int Ω 〈 (nn T - I) \nabla u, β 〉 d x

becomes

\int Ω 〈 \nabla u, F α, γ, a, b (n, ϕ, β) 〉 d x : = \int Ω 〈 \nabla u, α h a, b (n, ϕ) n + γ (〈 n, β 〉 n - β) 〉 d x

using the integral formula by part, the left hand term derives

\int Ω 〈 \nabla u, F α, γ, a, b (n, ϕ, β) 〉 d x = \int \partial Ω (uF α, γ, a, b (n, ϕ, β)) \cdot ν d x - \int Ω u div (F α, γ, a, b (n, ϕ, β)) d x

where ν is the outer normal direction along the boundary ∂Ω. Choosing

n \cdot ν | \partial Ω = 0

β \cdot ν | \partial Ω = 0

and using the first variation, we get the solution of u-subproblem

u = τ 1 + τ u k + 1 1 + τ (z + α \nabla \cdot (h a, b (n, ϕ) n) + γ \nabla \cdot (〈 n, β 〉 n - β))

n-subproblem

We first consider the first-order variation of the functional $φ 1 (n):= \int Ω 〈 \nabla \cdot n, ϕ 〉〈 \nabla u, n 〉 d x = \int Ω 〈 ϕ \nabla u, \nabla \cdot nn 〉 d x$ :

lim t \to 0 φ 1 (n + t m) - φ 1 (n) t = \int Ω 〈 \nabla \cdot n + t \nabla \cdot m, ϕ 〉 〈 \nabla u, n + t m 〉 dx - \int Ω 〈 \nabla \cdot n, ϕ 〉 〈 \nabla u, n 〉 dx t = \int Ω 〈 \nabla \cdot m, ϕ 〉 〈 \nabla u, n 〉 dx + \int Ω 〈 \nabla \cdot n, ϕ 〉 〈 \nabla u, m 〉 dx = \int Ω 〈 - \nabla (〈 \nabla u, n 〉 ϕ) + 〈 \nabla \cdot n, ϕ 〉 \nabla u, m 〉 dx + \int \partial Ω (〈 \nabla u, n 〉 ϕ) (m \cdot ν) dx

Further the first-order variation of functional

φ 2 (n):= \int Ω 〈 nn T \nabla u, β 〉 d x

can be derived as

lim t \to 0 φ 2 (n + t m) - φ 2 (n) t = \int Ω 〈 (n + t m) (n + t m) T \nabla u, β 〉 - 〈 nn T \nabla u, β 〉 d x t = \int Ω 〈 mn T \nabla u, β 〉 + 〈 nm T \nabla u, β 〉 d x = \int Ω 〈 (n T \nabla u) β + (n T β) \nabla u, m 〉 d x

Since we have equipped

n \cdot ν | \partial Ω = 0

on boundary ∂Ω, then the test function m should also satisfy

m \cdot ν | \partial Ω = 0

on boundary ∂Ω. Hence we obtain the Euler-Lagrangian equation

α (a - b | ϕ | 2) \nabla u - 2 α b (\nabla (〈 \nabla u, n 〉 ϕ) - 〈 \nabla \cdot n, ϕ 〉 \nabla u) + γ ((n T \nabla u) β + (n T β) \nabla u) + τ (n - n k) = 0

(18)

The solutions of ϕ-subproblem and β -subproblem can be easily computed by using the first variation in same way. The minimization problems for u, ϕ, and β can be done by some simple arithmetic calculations at each grid point. There is no need to solve any equations. The minimization problem for n needs to solve a linear equation by iterative scheme with very low cost, however such linear equation has the closed-form solution at each grid point. In general the following special Lemma can be used.

Lemma 4

Assume that $1 + μ T ω \neq 0 and (1 + μ T ω) 2 - | μ | 2 | ω | 2 \neq 0$ , then

(I + ω μ T + μ ω T) - 1 x = x - μ T x 1 + μ T ω ω - ((1 + μ T ω) ω T x - | ω | 2 μ T x) μ + (| μ | 2 | ω | 2 μ T x 1 + μ T ω - | μ | 2 ω t x) ω (1 + μ T ω) 2 - | μ | 2 | ω | 2

(19)

Proof

Firstly, let A be $I + ω μ T$ , from the Sherman-Morrison formula $(A + μ ω T) - 1 = A - 1 - A - 1 μ ω T A - 1 1 + ω T A - 1 μ (1 + ω T A - 1 μ \neq 0)$ , one has $(A + μ ω T) - 1 x = A - 1 x - A - 1 μ ω T A - 1 x 1 + ω T A - 1 μ$ , further $A - 1 = (I + ω μ T) - 1 = I - ω μ T 1 + ω T μ$ .

We address the optimization problem in (16) using an alternating iterative scheme with respect to u, ϕ, n , β and multiplier q . The complete resulting image denoising algorithm is summarized in Algorithm. Note that some steps will be reused when performing more complicated tasks such as image segmentation or inpainting.

Algoritham 1 (Augmented Lagrangian Primal-Dual method (ALPD)-elastica denoising)

Step 1. Input an observed image z. Set the initial multiplier q ⁰= 0. Given initial solutions $u 0 = z, n 0 = 0, ϕ 0 = 0, β 0 = 0$ . Set $θ 1 = τ / (1 + τ)$ , where α, γ are regularizing parameters, τ is a Lipschitz parameter;

Step 2: for k ≥ 0

u - subproblem : u k + 1 / 2 = z + α \nabla \cdot (h a, b (n k, ϕ k) n k) + γ \nabla \cdot (〈 n k, β k 〉 n k - β k) u k + 1 = θ 0 u k + (1 - θ 0) u k + 1 / 2 n - subproblem : Semi - implicititerativesolution (Lemma 4) n k + j / m = (I + ω k (μ k) T + μ k (ω k) T) - 1 x k + (j - 1) / m with ω k = γ / τ \nabla u k + 1, μ = β k, 1 \leq j \leq m \leq 5 where x k + (j - 1) / m = n k - α / τ [(a - b | ϕ k | 2) \nabla u k + 1 - 2 b (\nabla (〈 \nabla u k + 1, n k + (j - 1) / m 〉 ϕ k) - 〈 \nabla \cdot n k + (j - 1) / m, ϕ k 〉 \nabla u k + 1)] On x \in Ω, if 〈 n k + 1 (x), \nabla u k + 1 (x) 〉 < 0, then n k + 1 (x) = - n k + 1 (x); ϕ - subproblem : θ k = τ / (τ + 2 α b 〈 \nabla u k + 1 n k + 1 〉), ϕ k + 1 = θ k ϕ k + (1 - θ k) \nabla \cdot n k + 1 β - subproblem : β k + 1 / 2 = γ (c (u k + 1 n k + 1) + q k γ) β k + 1 = θ 0 β k + (1 - θ 0) β k + 1 / 2 multiplierupdate : q k + 1 = q k + γc (u k + 1, n k + 1) ⇆

Step 3. Stop if a given criterion is valid, otherwise go to Step 2.

Image segmentation

Euler’s elastica as a new regularization of segmentation contour was originally proposed and used in the famous in-depth work by Nitzberg et al³⁶ Recently, Zhu et al.³⁴ have developed a numerical method to minimize Euler’s elastica dependent functionals by using the augmented Lagrangian method. In their work, they also considered Chan-Vese’s model with the substitution of Euler’s elastica for the length term. As discussed in the introduction, the purpose of our present work is to present some new properties and implementation of this modified Chan-Vese’s model through numerical study with the aid of the Legendre–Fenchel transformation.³⁴

We know that the almost theoretical results for the segmentation case can be obtained in the same way as the denoising case of (16). Further, Euler’s elastica-dependent segmentation problem (8) can be rewritten as the following constraint optimization minimization in the equivalent form

min u, n, c 1, c 2 (α P 1 (u, n) + \int Ω g (u, f, c 1, c 2) d x) s . t . c (u, n) = 0

(20)

where

g (u, f, c 1, c 2) = u (f - c 1) 2 + (1 - u) (f - c 2) 2

. In order to efficiently solve the above constrained optimization problem, we employ the augmented Lagrangian functional and the Legendre–Fenchel transformation to obtain the following problem

min u, n, c 1, c 2 {α P 1 (u, n) + \int Ω g (u, f, c 1, c 2) d x + γ sup β \in L 2 (Ω) (\int Ω 〈 c (u, n) + q k γ, β 〉 - 12 | β | 2 d x)}

which can be solved in a similar way by adjusting the notation accordingly. In fact, the above optimization problem becomes a maximization-minimization problem, the subproblems with respect to ϕ, n , β and multiplier q can be computed in the same way as described above. Moreover for u-subproblems, we first compute the first variation of the above functional with respect to variant u and then obtain the Euler-Lagrangian equation with the boundary conditions

n \cdot ν | \partial Ω = 0

and

β \cdot ν | \partial Ω = 0

as follows

((f - c 1) 2 - (f - c 2) 2) - α \nabla \cdot (h a, b (n, ϕ) n) - γ \nabla \cdot (〈 n, β 〉 n - β) + τ (u - u k) = 0

which simplifies to an explicit solution

u = u k + 1 / τ (- ((f - c 1) 2 - (f - c 2) 2) + α \nabla \cdot (h a, b (n, ϕ) n) + γ \nabla \cdot (〈 n, β 〉 n - β))

With all the prior and likelihood information, the final segmentation decision is drawn by some simple thresholding ς, i.e., Ω _in = {x: u(x) ≥ ς} and Ω _out = {x: u(x) < ς}. The overall segmentation algorithm is summarized in Algorithm 2.

Algorithm 2 (ALPD-elastica segmentation)

Step 1. Given an initial image f, the initial multiplier q ⁰= 0 . Input the initial two-phase function u⁰, set initial solutions $n 0 = 0, ϕ 0 = 0, β 0 = 0 . Set θ 0 = τ / (1 + τ)$ θ₀ =τ/(1 + τ), where α, γ are regularizing parameters, τ is a Lipschitz parameter;

Step 2. for k ≥ 0

c 1, c 2 - meanvalue : c_{1}^{k} = \int Ω u k fdx \int Ω u k dx, c_{2}^{k} = \int Ω (1 - u k) fdx \int Ω (1 - u k) dx; u - subproblem : u k + 1 / 2 = u k + 1 / τ (- ((f - c_{1}^{k}) 2 - (f - c_{2}^{k}) 2) + α \nabla \cdot (h a, b (n k, ϕ k) n k) + γ \nabla \cdot (〈 n k, β k 〉 n k - β k)), u k + 1 = Proj [0, 1] (u k + 1 / 2); n - subproblem : Semi - implicititerativesolution (Lemma 4) n k + j / m = (I + ω k (μ k) T + μ k (ω k) T) - 1 x k + (j - 1) / m with ω k = γ / τ \nabla u k + 1, μ = β k, 1 \leq j \leq m \leq 5, wherex k + (j - 1) / m = n k - α / τ [(a - b | ϕ k | 2) \nabla u k + 1 - 2 b (\nabla (〈 \nabla u k + 1, n k + (j - 1) / m 〉 ϕ k) - 〈 \nabla \cdot n k + (j - 1) / m ϕ k 〉 \nabla u k + 1)], Onx \in Ω, if 〈 n k + 1 (x), \nabla u k + 1 (x) 〉 < 0, then n k + 1 (x) = - n k + 1 (x) ϕ - subproblem : θ k = τ / (τ + 2 α b 〈 \nabla u k + 1 n k + 1 〉) ϕ k + 1 = θ k ϕ k + (1 - θ k) \nabla \cdot n k + 1 β - subproblem : β k + 1 / 2 = γ (c (u k + 1 n k + 1) + q k γ) β k + 1 = θ 0 β k + (1 - θ 0) β k + 1 / 2 multiplierupdate : q k + 1 = q k + γc (u k + 1, n k + 1) L \to

Step 3. Stop if a given criterion is valid, otherwise go to Step 2.

Numerical results

In this section, we present some numerical results from applying our proposed algorithms. The experimental results include two parts: surface denoising and image segmentation. Let us first introduce the discrete setting which we will use in the rest of this section.

Discrete setting

We consider a regular Cartesian grid of size $m \times n : {(ih, jh) : 1 \leq i \leq m, 1 \leq j \leq n}$ , where h denotes the size of the spacing and (i, j) denote the indices of the discrete locations (ih, jh) in the image domain Ω: = [0, 1]. Let us define the dicrete gradient operator applied to a function u at the point (i, j)

(\nabla u) ij = ((\nabla u) ij 1, (\nabla u) ij 2)

where

Let us define the discrete divergence operator applied to a vector

v = (v 1, v 2)

at the grid point (i, j)

and define the vector inner product

〈 s i, j, t i, j 〉 = s_{ij}^{1} t_{ij}^{1} + s_{ij}^{2} t_{ij}^{2}

for any n × m × 2 matrices s and t .

Image denoising

We now use numerical simulations to illustrate the effectiveness of our denoising Algorithm 1 just developed. The numerical techniques are based on augmented Lagrangian and Legendre–Fenchel transformation. We follow the previous subsection and discretize the gradient and diffusion operators.

In order to demonstrate the functionality of the variant of our proposing algorithm, we applied it against a large variety of test problems. These problems include synthetic and natural images. For all these test problems (here all images are grayscale, and the original images are scaled to the range [0, 1]), the values of the noise function are uniformly distributed with mean variation σ > 0. The main compared approach is naturally the Tai-Hahn-Chung (THC) method proposed by Tai et al.³³ because our method in this work is inspired by it. Their approach is different from ours since it implies three extra augmentation functionals and the related Lagrange multipliers. The solution results restoring the synthetic images (“Smooth” and “Anthony”) natural image (“Pepper”) with the noise level $σ = 10255$ are shown in Figure 1. One can clearly see that the two algorithms only differ in small magnitudes around some parts on the edges of the shapes. Our approach is more practically efficient of eliminating noise for the very smooth image, and obtaining a comparable denoising quality when a natural image is tested.

Figure 1.

The denoising results for a synthetic image and a real image. The true and noisy images for three examples are listed in the first row; denoised and residual images by the THC method proposed by Tai et al.³³ are shown in the middle row; our results are shown in the third row. Test examples from left to right: “Smooth,” “Anthony” and “Pepper,” respectively. (a) True. (b) Noise. (c) THC. (d) Error. (e) ALPD. (f) Error. (g) True. (h) Noise. (i) THC. (j) Error. (k) ALPD. (l) Error. (m) True. (n) Noise. (o) THC. (p) Error. (q) ALPD. (r) Error.

In Table 1, we first compare the restoration quality (via psnr) and efficiency (via Iters) of two approaches by testing the above examples; in each approach four different image sizes are used to perform the major advantage of our numerical scheme. Experimental results show our algorithm can usually save the outer iteration to achieve the same relative error. Besides, the subproblems for the variable n in THC and ALPD require inner iterations for the convergence, while the approximations at every node (i, j) given by our algorithm are independent from each other instead of treating it with global-dependent update by THC algorithm. Further the inner iteration is only once for u-subproblem, ϕ-subproblem, β -subproblem and q -subproblem in our proposed scheme. Therefore, our algorithm can save computational cost in both inner iteration and outer iteration, which makes it superior to the previous algorithm.

Table 1.

Performance comparison of different examples and different image sizes by using our algorithm and THC algorithm.

	Smooth				Anthony				Pepper
	ALPD		THC³³		ALPD		THC³³		ALPD		THC³³
	psnr	Iters	psnr	Iters	psnr	Iters	psnr	Iters	psnr	Iters	psnr	Iters
64 × 64	39.61	78	39.67	128	31.96	112	31.78	145	31.75	69	31.71	113
128 × 128	42.52	69	42.34	104	33.17	109	33.25	138	32.51	86	32.47	129
256 × 256	43.15	65	42.95	98	34.42	101	34.23	143	33.75	82	33.56	142
512 × 512	43.23	77	42.91	119	36.29	85	36.13	127	34.26	75	34.15	121

THC: Tai-Hahn-Chung; ALPD: Augmented Lagrangian Primal-Dual method.

In Figure 2, plots of objective function versus iteration numbers (Iters), CPU times versus image sizes and noise variations versus psnr are shown for the “Smooth” example with the size 256 × 256 in Table 1. From the third plot (see Figure 2(c)), we can see that both algorithms produce the restored results of good quality with similar psnr, but the objective function in (3) of the proposed algorithm reduce faster and more stably compared to the THC algorithm (see Figure 2(a)). Our computational cost (see Figure 2(b)) is lower compared to the THC algorithm, which demonstrates that the proposed algorithm is more reliable and efficient in practice.

Figure 2.

The first plot (a) lists the plots of objective function values versus iterations for the example “Smooth” using our algorithm and the THC algorithm; the second plot (b) presents the corresponding plots of CPU times versus image sizes. The plot (c) shows the small difference of psnr between the two algorithms restoring different variation images.

It is well known that noise exists in computed tomography (CT) images. Figure 3 shows the results obtained while applying the implementation against CT images with noise. We conduct two experiments of Euler’s elastica model on real liver CT data. Three slices of liver and bone CT images are selected from two data sets, which are displayed at the top of Figure 3. The restored images are shown at the bottom of Figure 3. As illustrated by these two experiments, our proposed scheme can successfully remove the noise in CT images.

Figure 3.

The performance of our method removing the noise contained in CT medical images. (a) CT with noise. (b) Bone I with noise. (c) Bone II with noise. (d) CT denoising. (e) Bone I denoising. (f) Bone II denoising.

Image segmentation

In this subsection, we describe the application of the proposed algorithm solving Euler’s elastica regularized Mumford–Shah model with two-phase clustering on real camera and medical images.

The first experiment consists of a two-phase segmentation in a natural image shown in Figure 4(a). The objective is to distinguish the cameraman and the background. Figure 4(b) and 4(d) show the segmentation contour and the sampled regions by our proposed algorithm. Figure 4(c) and 4(e) show the segmentation contour and the sampled regions by Zhu-Tai-Chan (ZTC’s) algorithm.³⁴ In Figure 5, a liver CT image is used to compare both algorithms. The result demonstrates that the tissue target in Figure 5(a) is detected without false alarm by using the proposed algorithm and ZTC’s algorithm (see Figure 5(b)–5(e)). We see that, with the empirically optimal parameters, the differences in solution’s visual quality between each algorithm is very small, though our new and primal-dual algorithm is slightly better (more details available). However, in Table 2, where the iteration number and CPU times of the both algorithms are tested, our new algorithm performs more efficiently (with fewer iterations and less computation time).

Figure 4.

Euler’s elastica based image segmentation: camera-man natural image. (a) Original image. (b) The contour {x : u(x) = 0.5} by our algorithm. (c) The contour {x : u(x) = 0.5} by ZTC algorithm.³⁴ (d) Phase function u by our algorithm. (e) Phase function u by ZTC algorithm.³⁴

Figure 5.

Euler’s elastica-based image segmentation: liver CT image. (a) Original image. (b) The contour {x : u(x) = 0.5} by our algorithm. (c) The contour {x : u(x) = 0.5} by ZTC algorithm.³⁴ (d) Phase function u by our algorithm. (e) Phase function u by ZTC algorithm.³⁴

Table 2.

Iterations (Iters) and CPU times (CPU(s)) of different examples and different image sizes by using our Algorithm 2 and the ZTC algorithm.

size	Camera (see Figure 4(a))				Liver-CT (see Figure 5(a))				Bone-CT (see (a))
	ALPD		ZTC ³⁴		ALPD		ZTC ³⁴		ALPD		ZTC ³⁴
	Iters	CPU(s)	Iters	CPU(s)	Iters	CPU(s)	Iters	CPU(s)	Iters	CPU(s)	Iters	CPU(s)
256 × 256	58	8.47	105	25.12	45	7.78	88	22.31	66	9.06	123	27.94
512 × 512	61	29.33	112	90.51	43	28.63	102	101.81	72	33.15	135	128.72

ZTC: Zhu-Tai-Chan..

Figure 6.

Euler’s elastica-based image segmentation: bone micro CT image. (a) Original image. (b) The contour {x : u(x) = 0.5} by our algorithm. (c) The contour {x : u(x) = 0.5} by ZTC algorithm.³⁴

In some cases, one might want to do simulations on a sample containing trabecular bone. However, it is notable that some of the connectivity in the trabeculae bone is lost, which is quite important for some applications and affects the mechanical stiffness of this part. The problem is that preserving the small feature of connected complete trabeculae bone (the small feature being the thin trabeculae) is not possible when using a popular CV model based on total variation regularization. To illustrate the advantage of Euler’s elastica in integrating missing or broken parts,³⁴ we apply our algorithm to an image with incomplete trabeculae, as shown in Figure 6(a). Our method gives more satisfactory results in Figure 6 since our method works well at connecting the missing information and preserving more details in the bone image. It clearly shows that the proposed method is suitable for segmentation of trabeculae bone.

Conclusions

A novel implementation formula solving Euler elastica regularizer for image denoising and segmentation is proposed in this paper. There are two novelties. One is that a new constraint $(nn T - I) \nabla u = 0$ is introduced to represent the normal vector $n = \nabla u | \nabla u |$ , which can reduce the singularity of variant u and computational cost. The other is the Legendre–Fenchel transformations and augmented Lagrangian technique are incorporated into the Euler elastica energy functional. Furthermore, the proposed method is faster and easier to implement than the other surface evolution based and augmented Lagrangian based methods. According to our experimental results, we find that our approach is competitive with other state-of-the-art algorithms.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The research of the authors was supported by the UK EPSRC grant (EP/K036939/1) and the National Natural Science Foundation of China (NSFC project 11301447).

References

Liu

Lai

Pericleous

. On a modified diffusion model for noise removal. J Algorithm Comput Technol 2012; 6: 35–58.

Chen

Piccolomini

Zama

. Iterative constrained minimization for vectorial tv image deblurring. J Math Imaging Vision 2016; 54(2): 240–255.

Osher

Sole

Vese

. Image decomposition and restoration using total variation minimization and the H-1 norm. Multiscale Model Simul 2003; 1(3): 349–370.

Rudin

Osher

Fatemi

. Nonlinear total variation based noise removal algorithms. Physica D 1992; 60(1): 259–268.

Savage

Chen

. An improved and accelerated non-linear multigrid method for total-variation denoising. Int J Comput Math 2005; 82(8): 1001–1015.

Vogel

Oman

. Fast, robust total variation-based reconstruction of noisy, blurred images. IEEE Trans Image Process 1998; 7(6): 813–824.

Geman

. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Trans Pattern Anal Mach Intell 1984; 6(6): 721–741.

Tikhonov

Arsenin

. Solutions of ill-posed problems, New York: John Wiley & Sons, 1977.

Acar

Vogel

. Analysis of bounded variation penalty methods for ill-posed problems. Inverse Prob 1994; 10(6): 1217–1229.

10.

Agarwal

Gribok

Abidi

. Image restoration using l-1 norm penalty function. Inverse Prob Sci Eng 2007; 15(8): 785–809.

11.

Aujol

. Some first-order algorithms for total variation based image restoration. J Math Imaging Vision 2009; 34(3): 307–327.

12.

Zhang

Chen

. An iterative Lagrange multiplier method for constrained total-variation-based image denoising. SIAM J Numer Anal 2012; 50(3): 983–1003.

13.

Hömke

Frohn-Schauf

Henn

Total variation based image registration. In: Tai

Lie

Chan

(eds). Image Processing Based on Partial Differential Equations, Berlin: Springer, 2007, pp. 343–361.

14.

Frohn-Schauf

Henn

Witsch

. Multigrid based total variation image registration. Comput Vis Sci 2008; 11: 101–113.

15.

Pock

Urschler

Zach

. A duality based algorithm for TV-L1-optical-flow image registration. In: Medical Image Computing and Computer-Assisted Intervention, Springer, 2007, pp. 511–518.

16.

Garnett

Meyer

. Image decompositions using bounded variation and generalized homogeneous besov spaces. Appl Comput Harmon Anal 2007; 23(1): 25–56.

17.

Duval

Aujol

Vese

. Projected gradient based color image decomposition. Lect Notes Comput Sci 2009; 5567: 295–306.

18.

Guo W and Qiao LH. Inpainting based on total variation. In: International Conference on Wavelet Analysis and Pattern Recognition, Beijing, China, 2–4 November 2007, pp.939–943. IEEE: ICWAPR.

19.

Getreuer

. Total variation inpainting using split Bregman. IPOL 2012; 2: 147–157.

20.

Ghate

Achaliya

Raveendran

. An algorithm of total variation for image inpainting. IJCER 2012; 1(3): 124–130.

21.

Chan

Yip

Park

. Simultaneous total variation image inpainting and blind deconvolution. Int J Imaging Syst Technol 2005; 15(1): 92–102.

22.

Bresson

Esedoglu

Vandergheynst

. Fast global minimization of the active contour/snake model. J Math Imaging Vision 2007; 28(2): 151–167.

23.

Unger M, Pock T, Trobin W, et al. TVSeg-interactive total variation based image segmentation. In: BMVC. 2008 British machine vision conference, Leeds, pp.1–10.

24.

Chan

Marquina

Mulet

. High-order total variation-based image restoration. SIAM J Sci Comput 2000; 22(2): 503–516.

25.

Lysaker

Osher

Tai

. Noise removal using smoothed normals and surface fitting. IEEE Trans Image Process 2004; 13(10): 1345–1357.

26.

Chang

Tai

Xing

. A compound algorithm of denoising using second-order and fourth-order partial differential equations. NM-TMA 2009; 2(4): 353–376.

27.

Brito-Loeza

Chen

. Multigrid algorithm for high order denoising. SIAM J Imaging Sci 2010; 3: 363–389.

28.

Lysaker

Lundervold

Tai

. Noise removal using fourth-order partial differential equation with application to medical magnetic resonance images in space and time. IEEE Trans Image Process 2003; 12(12): 1579–1590.

29.

Osher

Burger

Goldfarb

. An iterative regularization method for total variation-based image restoration. Multiscale Model Simul 2005; 4(2): 460–489.

30.

Levien R. The elastica: a mathematical history, https://www2.eecs.berkeley.edu/Pubs/TechRpts/2008/EECS-2008-103.pdf (2008, accessed 6 September 2016).

31.

Mumford

Elastica and computer vision. In: Bajaj

(ed). Algebraic Geometry and its Applications, New York: Springer, 1994, pp. 491–506.

32.

Shen

Kang

Chan

. Euler’s elastica and curvature-based inpainting. SIAM J Appl Math 2003; 63(2): 564–592.

33.

Tai

Hahn

Chung

. A fast algorithm for Euler’s elastica model using augmented Lagrangian method. SIAM J Imag Sci 2011; 4(1): 313–344.

34.

Zhu

Tai

Chan

. Image segmentation using Euler’s elastica as the regularization. SIAM J Sci Comput 2013; 57(2): 414–438.

35.

Chan

Vese

. Active contours without edges. IEEE Trans Image Process 2001; 10(2): 266–277.

36.

Nitzberg

Mumford

Shiota

. Filtering, segmentation and depth (Lecture Notes in Computer Science, 662), Berlin, Heidelberg: Springer, 1993.

A new augmented Lagrangian primal dual algorithm for elastica regularization

Abstract

Keywords

Introduction

Previous works based on Euler’s elastica

Euler’s elastica denoising model33

A Euler’s elastica-based segmentation model

The proposed algorithm for Euler’s elastica regularization

Definition 1 (Legendre–Fenchel transformation)

Lemma 1

Proof

Lemma 2

Lemma 3

A new algorithm for image denoising

u-subproblem

n-subproblem

Lemma 4

Proof

Algoritham 1 (Augmented Lagrangian Primal-Dual method (ALPD)-elastica denoising)

Image segmentation

Algorithm 2 (ALPD-elastica segmentation)

Numerical results

Discrete setting

Image denoising

Image segmentation

Conclusions

Footnotes

Declaration of conflicting interests

Funding

References

Euler’s elastica denoising model³³