A new augmented Lagrangian primal dual algorithm for elastica regularization

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 ∫ Ω | A u - z | 2 d x

(1)

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, 23}

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 ):= ∫ Ω [ a + | ∇ · n | 2 ] | ∇ u | d x

(2)

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 ( α ∫ Ω [ a + b | κ ( u ) | 2 ] | ∇ u | d x + 1 2 ∫ Ω | u - z | 2 d x )

(3)

where the mean curvature κ ( u ) = ∇ · ∇ u | ∇ 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 ( α ∫ Ω [ a + b | ∇ · n | 2 ] | p | d x + 1 2 ∫ Ω ( v - z ) 2 d x ) s . t . v = u , p = ∇ u , n = m , | p | = m · p , | m | ∈ [ - 1 , 1 ]

(4)

Second, the augmented Lagrangian method solves the following

L ( v , u , m , p , n , β 1 , β 2 , β 3 , β 4 ) = 1 2 ∫ Ω ( v - z ) d x + α ∫ Ω [ a + b | ∇ · n | 2 ] | p | d x + r 1 ∫ Ω ( | p | - p · m ) 2 d x + ∫ Ω β 1 ( | p | - p · m ) d x + r 2 ∫ Ω | p - ∇ u | 2 d x + ∫ Ω β 2 · ( p - ∇ u ) d x + r 3 ∫ Ω ( v - u ) 2 d x + ∫ Ω β 3 ( v - u ) d x + r 4 ∫ Ω | n - m | 2 d x + ∫ Ω β 3 · ( 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 φ ( α ∫ Ω [ a + b | ∇ · ∇ H ( φ ) | ∇ H ( φ ) | | 2 ] | ∇ H ( φ ) | d x + 1 2 ∫ Ω [ 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 ( α ∫ Ω [ a + b | ∇ · ∇ u | ∇ u | | 2 ] | ∇ u | d x + 1 2 ∫ Ω 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 ( α ∫ Ω [ a + b | ∇ · n | 2 ] | p | d x + 1 2 ∫ Ω u | f - c 1 | 2 + ( 1 - u ) | f - c 2 | 2 d x ) s . t . p = ∇ u , n = ∇ u | ∇ u | , u ∈ [ 0 , 1 ] .

(8)

Then the augmented Lagrangian functional

L ( v , u , p , n , m , c 1 , c 2 , β 1 , β 2 , β 3 , β 4 ) = 1 2 ∫ Ω v | f - c 1 | 2 + ( 1 - v ) | f - c 2 | 2 d x ) + α ∫ Ω [ a + b | ∇ · n | 2 ] | p | d x + r 1 ∫ Ω ( | p | - p · m ) 2 d x + ∫ Ω β 1 ( | p | - p · m ) d x + r 2 ∫ Ω | p - ∇ u | 2 d x + ∫ Ω β 2 · ( p - ∇ u ) d x + r 3 ∫ Ω ( v - u ) 2 d x + ∫ Ω β 3 ( v - u ) d x + δ D ( v ) + r 4 ∫ Ω | n - m | 2 d x + ∫ Ω β 4 · ( 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 → ℝ ∪ { + ∞ } be any functional. Then

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

(10)

The Legendre–Fenchel transformation F * : S * → ℝ ∪ { + ∞ } is convex, regardless whether or not F is convex. For vectors y ∈ ℝ d , the well-known equality used in dual formualtions | y | 1 = max | ω | 1 ≤ 1 ω · 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 ( η ) = 1 2 ∫ Ω η 2 ( x ) d x η ∈ L 2 ( Ω )

f * ( φ ) = sup η ∈ L 2 ( Ω ) 〈 φ , η 〉 - 1 2 ∫ Ω η 2 d x = 1 2 ∫ Ω φ 2 d x

(11)

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 ( φ ) = 1 2 ∫ Ω φ 2 ( x ) d x from (11), it holds that

1 2 ∫ Ω ( ∇ · n ) 2 | ∇ u | d x = sup φ ∈ L 2 ( Ω ) ∫ Ω ( ∇ · n ) φ ( x ) | ∇ u | d x - 1 2 ∫ Ω | φ | 2 | ∇ u | d x

Proof

Take φ = ( ∇ · n ) | ∇ u | . Then (11) becomes

1 2 ∫ Ω ( ∇ · n ) 2 | ∇ u | d x = sup φ ∈ L 2 ( Ω ) ∫ Ω ( ∇ · n ) | ∇ u | φ ( x ) d x - 1 2 ∫ Ω | φ | 2 d x = sup φ | ∇ u | ∈ L 2 ( Ω ) ∫ Ω ( ∇ · n ) | ∇ u | φ ( x ) | ∇ u | d x - 1 2 ∫ Ω | φ | 2 | ∇ u | d x = sup φ ∈ L 2 ( Ω ) ∫ Ω ( ∇ · n ) φ ( x ) | ∇ u | d x - 1 2 ∫ Ω | φ | 2 | ∇ u | d x

Consequently from Lemma, we obtain a useful reformulation for (2) OPEN IN VIEWER which will be used below. We give two more simple lemmas to aid further discussion.

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

{ TeeV ( u ) = P 1 ( u , n ):= sup ϕ ∈ L 2 ( Ω ) ∫ Ω [ a + b ( 2 ( ∇ · n ) ϕ - ϕ 2 ) ] ( ∇ u · n ) dx , s . t . c ( u , n ):= ( nn T - I ) ∇ 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 OPEN IN VIEWER We next reformulate the constrained minimization in (14) by the method of augmented Lagrangian multipliers. This leads to the unconstrained optimization problem OPEN IN VIEWER 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 ) ∈ W { α P 1 ( u , n ) + γ sup β ∈ L 2 ( Ω ) ( ∫ Ω 〈 c ( u , n ) + q k γ , β 〉 - 1 2 | β | 2 d x ) + 1 2 ∫ Ω | 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 , ϕ , β ):= α ∫ Ω [ a + b ( 2 〈 ∇ · n , ϕ 〉 - | ϕ | 2 ) ] 〈 ∇ u , n 〉 d x + γ ( ∫ Ω 〈 c ( u , n ) + q k γ , β 〉 - 1 2 | β | 2 d x ) + 1 2 ∫ Ω | 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 , ϕ , β ) : = α ∫ Ω [ a + b ( 2 〈 ∇ · n , ϕ 〉 - | ϕ | 2 ) ] 〈 ∇ u , n 〉 d x + 1 2 ∫ Ω | u - z | 2 d x + γ ( ∫ Ω 〈 c ( u , n ) + q k γ , β 〉 - 1 2 | β | 2 d x ) + τ 2 ∫ Ω | u - u k | 2 d x + τ 2 ∫ Ω | n - n k | 2 d x - τ 2 ∫ Ω | ϕ - ϕ k | 2 d x - τ 2 ∫ Ω | β - β 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 OPEN IN VIEWER

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

u-subproblem

In u-subproblem, the energy functional

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

becomes

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

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

∫ Ω 〈 ∇ u , F α , γ , a , b ( n , ϕ , β ) 〉 d x = ∫ ∂ Ω ( uF α , γ , a , b ( n , ϕ , β ) ) · ν d x - ∫ Ω u div ( F α , γ , a , b ( n , ϕ , β ) ) d x

where ν is the outer normal direction along the boundary ∂Ω. Choosing n · ν | ∂ Ω = 0 , β · ν | ∂ Ω = 0 and using the first variation, we get the solution of u-subproblem

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

n-subproblem

We first consider the first-order variation of the functional φ 1 ( n ):= ∫ Ω 〈 ∇ · n , ϕ 〉 〈 ∇ u , n 〉 d x = ∫ Ω 〈 ϕ ∇ u , ∇ · nn 〉 d x :

lim t → 0 φ 1 ( n + t m ) - φ 1 ( n ) t = ∫ Ω 〈 ∇ · n + t ∇ · m , ϕ 〉 〈 ∇ u , n + t m 〉 dx - ∫ Ω 〈 ∇ · n , ϕ 〉 〈 ∇ u , n 〉 dx t = ∫ Ω 〈 ∇ · m , ϕ 〉 〈 ∇ u , n 〉 dx + ∫ Ω 〈 ∇ · n , ϕ 〉 〈 ∇ u , m 〉 dx = ∫ Ω 〈 - ∇ ( 〈 ∇ u , n 〉 ϕ ) + 〈 ∇ · n , ϕ 〉 ∇ u , m 〉 dx + ∫ ∂ Ω ( 〈 ∇ u , n 〉 ϕ ) ( m · ν ) dx

Further the first-order variation of functional φ 2 ( n ):= ∫ Ω 〈 nn T ∇ u , β 〉 d x can be derived as

lim t → 0 φ 2 ( n + t m ) - φ 2 ( n ) t = ∫ Ω 〈 ( n + t m ) ( n + t m ) T ∇ u , β 〉 - 〈 nn T ∇ u , β 〉 d x t = ∫ Ω 〈 mn T ∇ u , β 〉 + 〈 nm T ∇ u , β 〉 d x = ∫ Ω 〈 ( n T ∇ u ) β + ( n T β ) ∇ u , m 〉 d x

Since we have equipped n · ν | ∂ Ω = 0 on boundary ∂Ω, then the test function m should also satisfy m · ν | ∂ Ω = 0 on boundary ∂Ω. Hence we obtain the Euler-Lagrangian equation

α ( a - b | ϕ | 2 ) ∇ u - 2 α b ( ∇ ( 〈 ∇ u , n 〉 ϕ ) - 〈 ∇ · n , ϕ 〉 ∇ u ) + γ ( ( n T ∇ u ) β + ( n T β ) ∇ 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 ω ≠ 0 and ( 1 + μ T ω ) 2 - | μ | 2 | ω | 2 ≠ 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 μ ≠ 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 + α ∇ · ( h a , b ( n k , ϕ k ) n k ) + γ ∇ · ( 〈 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 = γ / τ ∇ u k + 1 , μ = β k , 1 ≤ j ≤ m ≤ 5 where x k + ( j - 1 ) / m = n k - α / τ [ ( a - b | ϕ k | 2 ) ∇ u k + 1 - 2 b ( ∇ ( 〈 ∇ u k + 1 , n k + ( j - 1 ) / m 〉 ϕ k ) - 〈 ∇ · n k + ( j - 1 ) / m , ϕ k 〉 ∇ u k + 1 ) ] On x ∈ Ω , if 〈 n k + 1 ( x ) , ∇ u k + 1 ( x ) 〉 < 0 , then n k + 1 ( x ) = - n k + 1 ( x ) ; ϕ - subproblem : θ k = τ / ( τ + 2 α b 〈 ∇ u k + 1 n k + 1 〉 ) , ϕ k + 1 = θ k ϕ k + ( 1 - θ k ) ∇ · 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 ) + ∫ Ω 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 ) + ∫ Ω g ( u , f , c 1 , c 2 ) d x + γ sup β ∈ L 2 ( Ω ) ( ∫ Ω 〈 c ( u , n ) + q k γ , β 〉 - 1 2 | β | 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 · ν | ∂ Ω = 0 and β · ν | ∂ Ω = 0 as follows

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

which simplifies to an explicit solution

u = u k + 1 / τ ( - ( ( f - c 1 ) 2 - ( f - c 2 ) 2 ) + α ∇ · ( h a , b ( n , ϕ ) n ) + γ ∇ · ( 〈 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 = ∫ Ω u k fdx ∫ Ω u k dx , c 2 k = ∫ Ω ( 1 - u k ) fdx ∫ Ω ( 1 - u k ) dx ; u - subproblem : u k + 1 / 2 = u k + 1 / τ ( - ( ( f - c 1 k ) 2 - ( f - c 2 k ) 2 ) + α ∇ · ( h a , b ( n k , ϕ k ) n k ) + γ ∇ · ( 〈 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 = γ / τ ∇ u k + 1 , μ = β k , 1 ≤ j ≤ m ≤ 5 , wherex k + ( j - 1 ) / m = n k - α / τ [ ( a - b | ϕ k | 2 ) ∇ u k + 1 - 2 b ( ∇ ( 〈 ∇ u k + 1 , n k + ( j - 1 ) / m 〉 ϕ k ) - 〈 ∇ · n k + ( j - 1 ) / m ϕ k 〉 ∇ u k + 1 ) ] , Onx ∈ Ω , if 〈 n k + 1 ( x ) , ∇ u k + 1 ( x ) 〉 < 0 , then n k + 1 ( x ) = - n k + 1 ( x ) ϕ - subproblem : θ k = τ / ( τ + 2 α b 〈 ∇ u k + 1 n k + 1 〉 ) ϕ k + 1 = θ k ϕ k + ( 1 - θ k ) ∇ · 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 →

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 × n : { ( ih , jh ) : 1 ≤ i ≤ m , 1 ≤ j ≤ 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)

( ∇ u ) ij = ( ( ∇ u ) ij 1 , ( ∇ u ) ij 2 )

where OPEN IN VIEWER Let us define the discrete divergence operator applied to a vector v = ( v 1 , v 2 ) at the grid point (i, j) OPEN IN VIEWER 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 σ = 10 255 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. OPEN IN VIEWER

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.

OPEN IN VIEWER

Table 1.

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

	Smooth				Anthony				Pepper
	ALPD		THC 33		ALPD		THC 33		ALPD		THC 33
	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. OPEN IN VIEWER

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. OPEN IN VIEWER

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). OPEN IN VIEWER

OPEN IN VIEWER

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. ³⁴

OPEN IN VIEWER

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 34		ALPD		ZTC 34		ALPD		ZTC 34
	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..

OPEN IN VIEWER

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 ) ∇ u = 0 is introduced to represent the normal vector n = ∇ u | ∇ 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.