Restricted Linearized Augmented Lagrangian Method for Euler's Elastica Model

TL;DR

This paper introduces a cutting-off strategy for augmented Lagrangian algorithms applied to Euler's elastica model, improving parameter tuning, convergence speed, and image restoration quality by addressing nonlinearity and nonconvexity challenges.

Contribution

A novel cutting-off strategy is proposed to enhance augmented Lagrangian algorithms for Euler's elastica, simplifying parameter tuning and accelerating convergence.

Findings

Faster convergence demonstrated in numerical experiments.

Easier parameter tuning compared to existing methods.

Higher quality image restorations achieved.

Abstract

Euler's elastica model has been extensively studied and applied to image processing tasks. However, due to the high nonlinearity and nonconvexity of the involved curvature term, conventional algorithms suffer from slow convergence and high computational cost. Various fast algorithms have been proposed, among which, the augmented Lagrangian based ones are very popular in the community. However, parameter tuning might be very challenging for these methods. In this paper, a simple cutting-off strategy is introduced into the augmented Lagrangian based algorithms for minimizing the Euler's elastica energy, which leads to easy parameter tuning and fast convergence. The cutting-off strategy is based on an observation of inconsistency inside the augmented Lagrangian based algorithms. When the weighting parameter of the curvature term goes to zero, the energy functional boils down to the ROF model. So, a natural requirement is that its augmented Lagrangian based algorithms should also approach the augmented Lagrangian based algorithms formulated directly for solving the ROF model from the very beginning. Unfortunately, this is not the case for certain existing augmented Lagrangian based algorithms. The proposed cutting-off strategy helps to decouple the tricky dependence between the auxiliary splitting variables, so as to remove the observed inconsistency. Numerical experiments suggest that the proposed algorithm enjoys easier parameter-tuning, faster convergence and even higher quality of image restorations.

Figures

Tables

Table 1. TABLE I: The parameters used by LALM and RALM

Image	Methods	$a$	$b$	$\lambda$	$r_1$	$r_2$	$r_3$	$\gamma$	${\delta }_1$	${\delta }_2$	$tol$
lena	LALM	1	0.01	15	200	5	1.2	1e-5	1e-2	2e-2	9e-5
	RALM	1	0.01	12.5	50	1	2	1e-5	5e-2	1e-2	9e-5
peppers	LALM	2	0.01	32	200	2	1.5	1e-5	1e-2	1e-2	2e-4
	RALM	1	0.01	13	50	1	2	1e-5	5e-2	1e-2	2e-4
cameraman	LALM	1	0.01	15	150	1	1.5	1e-5	5e-2	3e-3	5e-5
	RALM	1	0.01	11.6	50	1	2	1e-5	5e-2	1e-2	5e-5

Table 2. TABLE II: PSNR, NMAD and NRMSE for images in our experiments

Image	Index	Noisy Data	LALM	RALM
lena	NRMSE	0.2797	0.0329	0.0305
$512\times 512$	NMAD	0.1632	0.0488	0.0463
	PSNR	20.0658	29.3617	29.6942
	Iter.		590	115
peppers	NRMSE	0.2473	0.0243	0.0199
$512\times 512$	NMAD	0.1725	0.0468	0.0433
	PSNR	20.1624	30.2318	31.1161
	Iter.		215	69
cameraman	NRMSE	0.1546	0.0211	0.0191
$512\times 512$	NMAD	0.1630	0.0506	0.0487
	PSNR	20.3904	29.0483	29.4845
	Iter.		712	192

Equations

$$f=u+n,$$

$$\underset{u}{\min }{\int }_{\mathrm{\Omega }}|\nabla u|+\frac{\lambda }{2}{\int }_{\mathrm{\Omega }}(f-u{)}^2,$$

$$\underset{u}{\min }{\int }_{\mathrm{\Omega }}\left(a+b{\kappa }^2|\nabla u|\right)d\mathrm{\Omega }+\frac{\lambda }{2}{\int }_{\mathrm{\Omega }}(u-f{)}^{2}d\mathrm{\Omega },$$

$$\underset{u}{\min }{\int }_{\mathrm{\Omega }}a|\nabla u|d\mathrm{\Omega }+\frac{\lambda }{2}{\int }_{\mathrm{\Omega }}(u-f{)}^{2}d\mathrm{\Omega }.$$

$$\underset{v,\mathbf{p},\mathbf{m},\mathbf{n}}{\min }{\int }_{\mathrm{\Omega }}\left(a+b(\nabla \cdot \mathbf{n}{)}^2\right)|\mathbf{p}|+\frac{\lambda }{2}{\int }_{\mathrm{\Omega }}(v-u_0{)}^2$$

$$s.t.\mathbf{p}=\nabla u,\mathbf{n}=\mathbf{m},v=u,|\mathbf{p}|=\mathbf{m}\cdot \mathbf{p},|\mathbf{m}|\le 1.$$

$$\underset{u,\mathbf{p},\mathbf{n},h}{\min }{\int }_{\mathrm{\Omega }}\left(a+b{h}^2\right)|\mathbf{p}|+\frac{\lambda }{2}{\int }_{\mathrm{\Omega }}(u-f{)}^2$$

$$s.t.\mathbf{p}=\nabla u,\mathbf{p}=|\mathbf{p}|\mathbf{n},h=\nabla \cdot \mathbf{n}.$$

$$\mathcal{L}(u,\mathbf{p},\mathbf{n},h$$

$$+{\int }_{\mathrm{\Omega }}{\mathbf{\lambda }}_1\cdot (\mathbf{p}-|\mathbf{p}|\mathbf{n})+\frac{r_1}{2}{\int }_{\mathrm{\Omega }}|\mathbf{p}-|\mathbf{p}|\mathbf{n}{|}^2$$

$$+{\int }_{\mathrm{\Omega }}{\mathbf{\lambda }}_2\cdot (\mathbf{p}-\nabla u)+\frac{r_2}{2}{\int }_{\mathrm{\Omega }}|\mathbf{p}-\nabla u{|}^2$$

$$+{\int }_{\mathrm{\Omega }}{\lambda }_3(h-\nabla \cdot \mathbf{n})+\frac{r_3}{2}{\int }_{\mathrm{\Omega }}(h-\nabla \cdot \mathbf{n}{)}^2,$$

$$u^{k+1}$$

$$+\frac{r_2}{2}{\int }_{\mathrm{\Omega }}|{\mathbf{p}}^k-\nabla u{|}^2,$$

$${\mathbf{p}}^{k+1}$$

$$+{\int }_{\mathrm{\Omega }}{\mathbf{\lambda }}_2^k\cdot (\mathbf{p}-\nabla u^{k+1})+\frac{r_2}{2}{\int }_{\mathrm{\Omega }}|\mathbf{p}-\nabla u^{k+1}{|}^2$$

$$+{\int }_{\mathrm{\Omega }}{\mathbf{\lambda }}_1^k\cdot \left(\mathbf{p}-{\mathbf{n}}^k|\mathbf{p}|\right)+\frac{r_1}{2}{\int }_{\mathrm{\Omega }}{\left|\mathbf{p}-{\mathbf{n}}^k|\mathbf{p}|\right|}^2,$$

$${\mathbf{n}}^{k+1}$$

$$+\frac{r_1}{2}{\int }_{\mathrm{\Omega }}{\left|{\mathbf{p}}^{k+1}-|{\mathbf{p}}^{k+1}|\mathbf{n}\right|}^2$$

$$+{\int }_{\mathrm{\Omega }}{\lambda }_3^k(h^k-\nabla \cdot \mathbf{n})+\frac{r_3}{2}{\int }_{\mathrm{\Omega }}(h^k-\nabla \cdot \mathbf{n}{)}^2,$$

$$h^{k+1}$$

$$+{\int }_{\mathrm{\Omega }}{\lambda }_3^k(h-\nabla \cdot {\mathbf{n}}^{k+1})+\frac{r_3}{2}{\int }_{\mathrm{\Omega }}(h-\nabla \cdot {\mathbf{n}}^{k+1}{)}^2,$$

$${\mathbf{\lambda }}_1^{k+1}={\mathbf{\lambda }}_1^k+r_1({\mathbf{p}}^{k+1}-|{\mathbf{p}}^{k+1}|{\mathbf{n}}^{k+1}),$$

$${\mathbf{\lambda }}_2^{k+1}={\mathbf{\lambda }}_2^k+r_2({\mathbf{p}}^{k+1}-\nabla u^{k+1}),$$

$${\lambda }_3^{k+1}={\lambda }_3^k+r_3(h^{k+1}-\nabla \cdot {\mathbf{n}}^{k+1}).$$

$$\mathbf{n}=\frac{\mathbf{p}}{|\mathbf{p}|}\iff \mathbf{p}=|\mathbf{p}|\mathbf{n}.$$

$${\int }_{\mathrm{\Omega }}{\left|\mathbf{n}-\frac{\mathbf{p}}{|\mathbf{p}|}\right|}^2,{\int }_{\mathrm{\Omega }}|\mathbf{p}-|\mathbf{p}|\mathbf{n}{|}^2$$

$$\underset{u,\mathbf{p},\mathbf{n},h}{\min }{\int }_{\mathrm{\Omega }}\left(a+b{h}^2\right)|\mathbf{p}|+\frac{\lambda }{2}{\int }_{\mathrm{\Omega }}(u-f{)}^2$$

$$s.t.\mathbf{p}=\nabla u,\mathbf{n}=\frac{\mathbf{p}}{|\mathbf{p}|},h=\nabla \cdot \mathbf{n},$$

$$\mathcal{L}(u,\mathbf{p},\mathbf{n},h$$

$$+{\int }_{\mathrm{\Omega }}{\mathbf{\lambda }}_1\cdot (\mathbf{n}-\frac{\mathbf{p}}{|\mathbf{p}|})+\frac{r_1}{2}{\int }_{\mathrm{\Omega }}|\mathbf{n}-\frac{\mathbf{p}}{|\mathbf{p}|}{|}^2$$

$$+{\int }_{\mathrm{\Omega }}{\mathbf{\lambda }}_2\cdot (\mathbf{p}-\nabla u)+\frac{r_2}{2}{\int }_{\mathrm{\Omega }}|\mathbf{p}-\nabla u{|}^2$$

$$+{\int }_{\mathrm{\Omega }}{\lambda }_3(h-\nabla \cdot \mathbf{n})+\frac{r_3}{2}{\int }_{\mathrm{\Omega }}(h-\nabla \cdot \mathbf{n}{)}^2,$$

$$u^{k+1}=$$

$$+\frac{r_2}{2}{\int }_{\mathrm{\Omega }}|{\mathbf{p}}^k-\nabla u{|}^2,$$

$${\mathbf{p}}^{k+1}=$$

$$+\frac{r_2}{2}{\int }_{\mathrm{\Omega }}|\mathbf{p}-\nabla u^{k+1}{|}^2,$$

$${\mathbf{n}}^{k+1}=$$

$$+{\int }_{\mathrm{\Omega }}{\lambda }_3^k(h^k-\nabla \cdot \mathbf{n})+\frac{r_3}{2}{\int }_{\mathrm{\Omega }}(h^k-\nabla \cdot \mathbf{n}{)}^2,$$

Full text

I Introduction

The classical image denoising problem is that, given a noisy image $f:\Omega\rightarrow\mathbb{R}$, where $\Omega$ is a bounded open subset of $\mathbb{R}^{2}$, a clean image $u$ is to be estimated from

[TABLE]

where $n$ is the noise. One of the most popular and successful variational methods for solving this denoising problem is developed by Rudin, Osher, and Fatemi(ROF) [1]. It is defined as the following minimization problem

[TABLE]

where $\lambda>0$ is a scalar weighting parameter. The ROF model is effective on preserving edges while removing noise thanks to the total variation term in the energy functional. Several fast algorithms have been proposed for solving the ROF model [2, 3, 4, 5].

Despite its appealing properties, the ROF model yields the staircase effect and fails to preserve image contrasts. To overcome these drawbacks, some high order variational models have been proposed [6, 7, 8, 9, 10, 11, 12]. Among them, Euler’s elastica energy functional has a wide application in image processing, such as image inpainting [13, 14], segmentation [15, 16], restoration [12, 17] and reconstruction [18, 19]. For image restoration, the Euler’s elastica model is defined as

[TABLE]

where $a,b,\kappa=\nabla\cdot\left(\frac{\nabla u}{|\nabla u|}\right)$ are the weighting parameters for the total variational term and the curvature term respectively. The Euler’s elastica functional was first introduced by Mumford, Nitzberg and Shiota [20] for segmenting images with occlusions. The conventional algorithms [21, 22, 13] for solving the Euler’s elastica model seek effective discretization schemes for solving the model’s Euler-Lagrange equation. However, due to the nonlinear and nonconvex curvature term, the involved computations are of high complexity and time-consuming.

Recently, several fast numerical algorithms [23, 24, 25, 26] have been proposed based on the augmented Lagrangian and other techniques, such as primal-dual, operator splitting or linearization techniques to convert the original problem into a few simple and easily-solved subproblems. Tai, Hahn and Chung [27] reformulated (2) as a constrained optimization problem by using the idea of augmented Lagrangian. Then, an alternating minimization method reduces the original problem into a series of subproblems that either have closed-form solutions or can be dealt with by the Fast Fourier Transform(FFT). Another related work can be found in [28], where Jianping Zhang and Ke Chen proposed an augmented Lagrangian-based primal dual approach, such that the number of auxiliary variables was reduced. To propose a simpler method with fewer parameters, Maryam Yashtini and Sung Ha Kang [14] tried to relax the curvature term and solve the associated subproblems by FFT and a shrinkage operator. Besides, they also tried to consider the Euler’s elastica energy minimization problem as a weighted ROF problem, so as to develop fast algorithms. However, the application of FFT is limited to periodic boundary conditions, which is not desirable for many image processing models. To relieve this concern, Tai [29] proposed another fast algorithm for Euler’s elastica model by replacing the FFT with an inexpensive arithmetic operation. In the same spirit, Duan et al. [30] introduced one new variable to decouple the mean curvature term. By using one sweep of the Gauss-Seidel iteration to replace the FFT operation, Duan’s method further reduced the needed computation cost. Inspired by this, Zhang et al. [31] adopted the linearization technique to approximate the quadratic terms in the Euler-Lagrange equation, such that all subproblems have closed-form solutions.

Despite the achieved promising results, the fast solvers developed so far for the Euler’s elastica model (2) might suffer from tricky parameter tuning due to the mutual dependence between the auxiliary variables. In this paper, we develop a cutting-off strategy for the augmented Lagrangian based algorithms for solving (2), which enjoys easier parameter tuning, faster convergence, and even higher quality of restorations. The proposed algorithm is based on a key observation of inconsistency inside the conventional augmented Lagrangian approach for solving (2). When $b\rightarrow 0$, the Euler’s elastica model (2) approaches the following ROF model (3)

[TABLE]

Nevertheless, the augmented Lagrangian method for solving (2) doesn’t reduce to the augmented Lagrangian method for solving (3). This will be explained and verified in the following sections. By introducing a cutting-off strategy, the restricted augmented Lagrangian method developed in this paper aims to remove such an inconsistency. Besides, we will show that a technique often used by conventional augmented Lagrangian based algorithms for simplifying the computations involving the auxiliary variable $\boldsymbol{n}$, is actually problematic. The proposed method also deals with this issue.

The remainder of this paper is organized as follows. In Section 2, the proposed method is described in detail. Numerical experiments to illustrate the effectiveness and efficiency of the proposed method are presented in Section 3. Conclusions and remarks are given in Section 4.

II The Restricted Augmented Lagrangian method for minimizing the Euler’s elastica energy

II-A Conventional Augmented Lagrangian Approaches

Tai et al. [27] applied the augmented Lagrangian method [4] to solve the Euler’s elastica model. This was done by introducing four auxiliary variables $\boldsymbol{p}$, $\boldsymbol{m}$, $v$ and $\boldsymbol{n}$ to relax the energy functional (2) as

[TABLE]

Instead of introducing four auxiliary variables, Duan et al. [30] developed an algorithm that needs just three auxiliary variables, and avoids the use of FFT, which looses the periodic boundary condition restriction. This algorithm starts by reformulating (2) as

[TABLE]

Then, the augmented Lagrangian formulation for (LABEL:DWH) reads

[TABLE]

where, $\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{2}\in\mathbb{R}^{2}$, and ${\lambda}_{3}\in\mathbb{R}$ are Lagrangian multipliers corresponding to the constraints $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$, $\boldsymbol{p}=\nabla u$ and $h=\nabla\cdot\boldsymbol{n}$, respectively, and for $i=1,2,3$, $r_{i}$ and $\lambda$ are positive scalar parameters. The above energy functional is then minimized by the alternating minimization method, i.e. given the previous iterate ($u^{k},\boldsymbol{p}^{k},\boldsymbol{n}^{k},h^{k};\boldsymbol{\lambda}_{1}^{k},\boldsymbol{\lambda}_{2}^{k},{\lambda}_{3}^{k}$), the next iterate is computed by solving the following subproblems (7)-(10) and updating the Lagrangian multipliers through (11)-(13).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The $\boldsymbol{p}$ subproblem can be solved by applying a shrinkage operator. The $u$ and $\boldsymbol{n}$ subproblems don’t possess closed-form solutions, so one sweep of Gauss-Seidel iteration is employed to compute approximate ones.

To further reduce computational complexities, Zhang et al. proposed a linearized augmented Lagrangian method [31] for Euler’s elastica model. Due to the use of the linearization technique, all the associated subproblems have closed-form solutions. In terms of image quality, Zhang’s algorithm achieves similar results as Tai’s method if parameters are chosen properly.

For the conventional algorithms mentioned above, we have the following two key observations.

Observation 1: The constraint associated with the variable $\boldsymbol{n}$ is $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$, which is a relaxation to $\boldsymbol{n}=\boldsymbol{p}/|\boldsymbol{p}|$. This is an essential technique to reduce the computational complexity for the optimization algorithms, which has been extensively utilized in the literature [27, 30, 14, 31]. This technique seems quite reasonable since the two equations are mathematically equivalent to each other when $\boldsymbol{p}\neq 0$, i.e.

[TABLE]

When involved into the optimization process, however, they could demonstrate different behaviors, i.e. the two terms

[TABLE]

are not equivalent to each other when they perform as penalty terms for $\boldsymbol{n}$ in the optimization procedure. For instance, when $\boldsymbol{p}\to\boldsymbol{0},\boldsymbol{p}\neq\boldsymbol{0}$, the second term will lose constraint over $\boldsymbol{n}$, while the first term requires that $|\boldsymbol{n}|=1,$ regardless of the value of $\boldsymbol{p}$. This observation indicates that replacing $\boldsymbol{n}=\boldsymbol{p}/|\boldsymbol{p}|$ by $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$ might be problematic. More discussions will be given later.

Observation 2: For the Euler’s elastica model, either in the form (LABEL:THC) or (LABEL:DWH), when $b\to 0$, it reduces to the ROF model, and only one auxiliary variable $\boldsymbol{p}$ and one multiplier $\boldsymbol{\lambda}$ are needed to be introduced. So, a consistent augmented Lagrangian method for the Euler’s elastica model should correspondingly reduce to the augmented Lagrangian method for solving the ROF model. However, when applying the alternating optimization technique, the subproblem for the optimization with respect to $\boldsymbol{p}$ obviously depends on $\boldsymbol{n}$, regardless of the value of parameter $b$. That’s to say, the augmented Lagrangian method for the Euler’s elastica model doesn’t immediately reduce to the augmented Lagrangian method for the ROF model when $b=0$. This is an internal inconsistency for applying the augmented Lagrangian approach directly for solving the Euler’s elastica model, which we think is the main reason responsible for the tricky coupling between the parameters and slow convergence. More discussions will be given later.

II-B The Restricted Augmented Lagrangian Method

In light of Observation 1 and Observation 2, in this subsection we will develop an efficient augmented Lagrangian based algorithm for solving the Euler’s elastica model. The proposed algorithm consists of three key ingredients. The first one accounts for Observation 1, and the original constraint $\boldsymbol{n}=\boldsymbol{p}/|\boldsymbol{p}|$ is used, instead of the relaxed version $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$. The second one accounts for Observation 2, and a simple cutting-off strategy is employed to recover the favourable property that the augmented Lagrangian methods for solving the Euler’s elastica model should reduce to the augmented Lagrangian methods for solving the ROF model when $b=0$. The third ingredient is the linearization technique utilized to reduce the computational complexity, as suggested in [31].

Let’s first reformulate the Euler’s elastica model as the following constraint minimization problem

[TABLE]

where the constraint for $\boldsymbol{n}$ takes its original form, rather than the relaxed one. Then, the augmented Lagrangian energy functional related to this minimization problem reads

[TABLE]

where, $\boldsymbol{\lambda}_{1},\boldsymbol{\lambda}_{2}\in\mathbb{R}^{2}$, and ${\lambda}_{3}\in\mathbb{R}$ are Lagrangian multipliers corresponding to the constraints $\boldsymbol{n}=\boldsymbol{p}/|\boldsymbol{p}|$, $\boldsymbol{p}=\nabla u$ and $h=\nabla\cdot\boldsymbol{n}$, respectively, and $r_{i}$, $i=1,2,3$ and $\lambda$ are positive scalars as weighting parameters. It should be noted that, to avoid the degenerate case of $\boldsymbol{p}=0$, it’s a common practice to replace $|\boldsymbol{p}|$ by $|\boldsymbol{p}|_{\epsilon}=|\boldsymbol{p}|+\epsilon$, where $\epsilon$ is a small positive constant. In this paper, we fix $\epsilon=1e-4$ in all the experiments.

As usually done in the literature, we adopt the alternating minimization method to minimize (15), just as described by (7)-(13). However, the minimization subproblem with respect to $\boldsymbol{p}$ will be modified to eliminate its connections to the variable $\boldsymbol{n}$. This is fulfilled by simply cutting-off the terms involving $\boldsymbol{n}$ in (8), which results in the following four sequential subproblems and three updating equations for the Lagrangian multipliers

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Please note that, when $b=0$, the couplings between the variable $\boldsymbol{n}$ and $\boldsymbol{p}$ will be removed, such that the subproblems for $u^{k+1}$ and $\boldsymbol{p}^{k+1}$ form a closed iteration cycle, which coincides exactly with the augmented Lagrangian method for solving the ROF model.

There is an intuitive explanation for the rationality of applying the cutting-off strategy. There is an ordered dependence between the auxiliary variables $\boldsymbol{p}$, $\boldsymbol{n}$ and $h$, which is usually ignored in the literature. The variable $\boldsymbol{p}$ can be thought of as an essential variable, while $\boldsymbol{n}$ is defined as $\boldsymbol{p}$/$|\boldsymbol{p}|$ to represent some normal field, and $h$ is introduced to denote $\nabla\cdot\boldsymbol{n}$. The role of the cutting-off strategy mentioned above is to restrict the information flowing from $\boldsymbol{n}$ to $\boldsymbol{p}$ when solving the $\boldsymbol{p}$ subproblem, which is reasonable since $\boldsymbol{p}$ should not depend on $\boldsymbol{n}$.

From now on, we will name our method as RALM (restricted augmented Lagrangian method) for brevity.

II-C The Solutions of Minimization Subproblems

In [31], Zhang et al. proposed a linearization technique to efficiently solve the subproblems derived from the augmented Lagrangian formulations. It has been shown that the linearization technique could achieve comparable results as the method proposed in [30], with much reduced computational complexity. So, the linearization technique is also adopted here to develop fast solvers for our subproblems, i.e. the $u$-subproblem and $\boldsymbol{n}$-subproblem. For easy reference, the linearization procedure [31] is described below.

For the $u$-subproblem Due to the nonlinearity of $|\boldsymbol{p}^{k}-\nabla u|^{2}$, it’s difficult to figure out an efficient solver for (16). So, we approximate this term as follows:

[TABLE]

where $\delta_{1}>0$ is a constant, which denotes the step-size for updating $u$. By applying the above approximation, the Euler-Lagrange equation of (16) can then be simplified to

[TABLE]

which admits a closed-form solution

[TABLE]

where $g_{1}$ denotes the right-hand side of (23).

For the $\boldsymbol{p}$-subproblem The $\boldsymbol{p}$-subproblem can be solved by the soft thresholding operator defined in [3], i.e. the solution of (17) can be explicitly expressed as

[TABLE]

where

[TABLE]

and the constant $c=\left(a+b\left(h^{k}\right)^{2}\right)$.

For the $\boldsymbol{n}$-subproblem Note that the constraint for $\boldsymbol{n}$ has changed back to its original form, i.e.

[TABLE]

This is the main difference from the derivation described in [31]. To handle the singularity of the operator $\nabla(\nabla\cdot)$, a quadratic penalty term shall be added. The $\boldsymbol{n}$-subproblem can then be reformulated as

[TABLE]

where $\gamma>0$ is a penalty parameter. As in solving the $u$-subproblem, we use the linearization technique to approximate the penalty term $\int_{\mathrm{\Omega}}\frac{r_{3}}{2}(h^{k}-\nabla\cdot\boldsymbol{n})^{2}$ as

[TABLE]

where $\delta_{2}>0$ is a constant. By applying the above approximation, the Euler-Lagrange equation of (18) can be simplified as

[TABLE]

which admits the following closed-form solution

[TABLE]

where $g_{2}$ denotes the right-hand side of (26).

For the $h$-subproblem The corresponding Euler-Lagrange equation for (19) can be simply derived as

[TABLE]

from which, we obtain

[TABLE]

The above solution procedures can be summarized as Algorithm 1. For easy reference, the algorithm proposed in [31] will be named as LALM (linearized augmented Lagrangian method). To validate the advantage of the constraint $\boldsymbol{n}=\boldsymbol{p}/|\boldsymbol{p}|$ over $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$, we also need adapt RALM to use the latter constraint. Such an algorithm will be named LALMn, which is summarized as Algorithm 2.

II-D Numerical discretization

For numerical implementations, all the computations in Algorithm 1 and Algorithm 2 need to be discretized on some grid. To simplify the expressions, we use same notations to denote both continuous and discrete quantities like the image $u$, gradient operator $\nabla u$, etc. The actual meaning of these operators should be clear from the context. Suppose the image we are dealing with is size of $M\times N\triangleq J$, and denote the Euclidean space $X=\mathbb{R}^{M\times N}$ and $Y=X\times X$. For given $u\in X$ and $\boldsymbol{v}\in Y$, we use the discrete forward (+) and backward (-) differential operators for the discretization of the gradient $\nabla u$ and the divergence $\nabla\cdot\boldsymbol{v}$ operators, respectively. With Neumann boundary conditions, the above mentioned operators are defined as follows

[TABLE]

with

[TABLE]

[TABLE]

for $i=1,\cdots,M$, $j=1,\cdots,N$, and

[TABLE]

where

[TABLE]

[TABLE]

for $i=1,\cdots,M$, $j=1,\cdots,N$.

III Experiments

Various experiments will be performed against LALM [31] in this section to validate the proposed algorithm. We will show that, compared with LALM, the proposed method RALM enjoys easier parameter tuning, less iteration numbers, and recovers images with higher quality. RALM will not be compared direclty with the fast algorithm proposed in [30], since it has been compared carefully with LALM in [31].

A synthesized image consisting of several circular rings will be utilized to carry out experiments to demonstrate the effect of the constraint adopted by RALM, as well as how LALM and RALM perform when the parameter $b\to 0$. Experiments on real images will be then performed to further validate the efficiency and effectiveness of the proposed method.

To compare the results quantitatively, we have computed three quality indices: PSNR, NRMSE and NMAD, which are deinfed by (31), (32) and (33), respectively. All the gray images are scaled to the range $[0,1]$ and correspondingly the display windows of all the restored images are also set to $[0,1]$. When a noisy image is needed, the Gaussian white noise with variance 0.01 is always applied. In all the following experiments, we use the relative residuals defined by (34) as the stopping criterion. To check the convergence behaviors of the competing methods, we also monitor the PSNR and other two indices: the numerical energy and and the norm of $\boldsymbol{n}$, defined by (30) and (35), respectively, to see how they evolve with the iterations.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The computing system is a desktop computer that possesses a 2.4GHz dual-core CPU and 16GB memory. The algorithms are implemented in MATLAB.

III-A Easier parameter tuning for the proposed method RALM

The Euler’s elastica model (LABEL:REE) involves three parameters: $a,b$, and $\lambda$. Due to the use of the linearization technique, six more parameters: $r_{1},r_{2},r_{3},\gamma,\delta_{1},\delta_{2}$ are introduced for both LALM and RALM. While the LALM can achieve promising results with properly chosen parameters, these parameters are tricky to tune. Thanks to the cutting-off strategy, the dependence between variables $\boldsymbol{p}$ and $\boldsymbol{n}$ is restricted to one direction, such that RALM has a low parameter sensibility. So, RALM should enjoy easier parameter tuning. To demonstrate this property, a sensitivity test is performed as follows. The two competing algorithms LALM and RALM are firstly tested on a synthesized image size of $512\times 512$ shown in the first row of Fig.1. For this image, the parameters for both the LALM and RALM are fine-tuned to achieve similar high quality restorations. As shown in Fig.1(a), both methods achieve quite similar PSNR values. Then the parameters are fixed, and the two methods are applied to three other images, namely pirate ($512\times 512$), boat ($1024\times 1024$) and montage ($256\times 256$), to see how the performances of the two methods adapt to different images.

The restored images after 200 iterations are shown in the last two columns of Fig.1, while the PSNR curves are demonstrated in Fig.2. At the first glance, both methods recover quite similar images after 200 iterations. When zooming in the resulting images and checking more closely, however, the higher quality shall be revealed for the results of RALM. For example, in the zoomed-in subregions shown at the upper right corners of the resulting images, RALM produces smoother local edges, while LALM produces small noisy structures that look like artifacts. The higher quality can also be verified by the PSNR plots shown in Fig.2. An apparent observation is that, after growing to the maximum, a sharp drop of the PSNR plots for LALM occurs for the tests with the pirate and boat images, which indicates that the optimal parameters of LALM are sensitive to different images. On the other hand, the PSNR plots of RALM suffer from much smaller drops, which indicates that the parameters of RALM are robust across different images, i.e. they are less sensitive than those of LALM. Note that for the montage image, both LALM and RALM produce rather similar results. We think this might be because the montage image shares more similar properties with the synthesized image, i.e. both images present few textures, and can be well approximated by piecewise smooth functions.

III-B Advantage Validation of the constraint $\boldsymbol{n}={\boldsymbol{p}}/{|\boldsymbol{p}|}$ over $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$

As stated previously, even the two constraints $\boldsymbol{n}=\boldsymbol{p}/|\boldsymbol{p}|$ and $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$ are mathematically equivalent to each other (when $\boldsymbol{p}\neq 0$), they might behave quite differently when being used as penalty terms. This can be illustrated by the following comparison experiment with LALMn and LALM. The former adopts the constraint $\boldsymbol{n}={\boldsymbol{p}}/{|\boldsymbol{p}|}$, just as described in Algorithm 2, while the latter adopts the constraint $\boldsymbol{p}=|\boldsymbol{p}|\boldsymbol{n}$.

The two algorithms are tested on a synthesized image consisting of several concentric circle rings, size of $512\times 512$. Fig.3(a) and (b) show the synthesized image and its noisy version, while Fig.3(c) and (d) show the restored images by LALM and LALMn, respectively. Please note that the vertical blue line shown in Fig.3(a) corresponds to the 256th column of the image, which shall be used for profile plots. From Fig.3(c) and (d), one can see that both algorithms could recover the ideal image well, and little difference could be visually identified. To examine the results more carefully, in Fig.4(a)-Fig.4(c), the line profiles of 256th columns and the norm of $\boldsymbol{n}$, which is defined by (35), are shown for the results of LALM and LALMn, respectively. All the plots are drawn up to 1000 iterations. From the line profile plots, observable differences, though still quite weak, can be identified, which shows that the LALMn plot agrees better with the ideal image than the LALM plot does. We have tried our best to fine tune the parameters of LALM, however, the difference persists and can not be diminished further.

A big difference occurs when we check the plots of the norm $\|\boldsymbol{n}\|$ produced from the two algorithms. As shown in Fig.4(c), the norm curve for LALMn converges to somewhere almost 0, while the norm curve for LALM converges to somewhere near 1. In fact, without counting the values near the edges, the norm for RALM should have converged to zero. We think that $\|\boldsymbol{n}\|\to 0$ is more reasonable. This can be justified by the following arguments. As noted in [32], instead of regarding an image as a function $f$ defined on a domain $\Omega\subset\mathbb{R}^{2}$, we can also consider the image $f$ as defining a surface in $\mathbb{R}^{3}$: $z=f(x,y)$, $(x,y)\in\Omega$. So the normal vector of the image surface takes the form of $\boldsymbol{\widetilde{n}}=(f_{x},f_{y},-1)$. For a piecewise constant region of $f$, the surface normal $\boldsymbol{\widetilde{n}}$ should be $(0,0,-1)$. In this regard, for the surface defined by the synthesized image, we should have $\boldsymbol{\widetilde{n}}=(0,0,-1)$ almost everywhere, except near the edges. Since $\boldsymbol{n}$ represents $\frac{(f_{x},f_{y})}{\sqrt{f_{x}^{2}+f_{y}^{2}+\epsilon}}$, it’s norm should be [math] almost everywhere.

III-C Inconsistency test for the conventional augmented Lagrangian formulation

As stated in the introduction section, there exists an internal inconsistency in the augmented Lagrangian formulation for the Euler’s elastica model. When the weighting parameter $b\to 0$, the Euler’s elastica model boils down to the ROF model. So, a natural requirement is that its augmented Lagrangian formulation should correspondingly boil down to the augmented Lagrangian formulation for the ROF model. This natural requirement, however, is not fulfilled by the conventional formulation, which we think is an internal inconsistency of the formulation. Since LALM is faithful to the conventional formulation, this inconsistency lives inside it. On the other hand, by introducing a cutting-off strategy, the proposed method RALM has successfully removed such an inconsistency. This can be easily seen from the pseudo-code described in Algorithm 1.

In this subsection, we will test how the inconsistency will influence the convergence behavior of LALM as well as the quality of the restored images. The synthesized image shown in Fig.3(a) is utilized again to test both LALM and RALM. With $b=0$, the Euler’s elastica model reduces to

[TABLE]

and the augmented Lagrangian functional associated with this optimization problem is

[TABLE]

which has been discussed in [4]. Since only the parameters $a,\lambda,\boldsymbol{\lambda}_{2},r_{2},\delta_{1}$ are involved in (36), it’s expected that the other left parameters $r_{1},r_{3},\delta_{2}$ should not influence the behavior of LALM and RALM. Since the connection between $\{u,p\}$ and the left auxiliary variables are established by the parameters $r_{1}$, we would like to test the methods with varied $r_{1}$. We monitor three quantitive indices: numerical energy, PSNR and relative errors, to see how they change with $r_{1}=50,500\mbox{ and }5000$, respectively. The results are shown in Fig.5. We observe that the numerical energy, PSNR and relative errors of LALM are all influenced by the values of $r_{1}$, while our algorithm RALM does not. From Fig.5(a) and Fig.5(b), one may conclude that even though the convergence behavior changes with varied $r_{1}$, it seems that both RALM and LALM converge to the same point. However, the result shown in Fig.5(c) strongly suggests that this should be just an illusion, and the final results of LALM will depend on $r_{1}$, at least to some extend.

III-D Performance tests on real images

In this subsection, we will test how the proposed RALM performs on real images against LALM. The clean images and their noisy counterparts are shown in Fig.6. The parameters used by the two competing methods, which are listed in Table I, have subjected to fine-tuning to achieve optimal performance. One can easily check that the parameters of RALM are almost immune to the changing of images, and only slight changes are needed. As comparison the parameters of LALM demonstrate strong image dependence, and significant adjustments are usually needed from image to image.

The restored results by the two methods are shown in Fig.7. From the left two columns, one can conclude that the two competing methods recover quite similar images, and little difference can be visually identified. In the zoomed-in sub-images located at the upper right corners, weak but perceptible differences can be detected. When we check the isophotes shown in the last two columns, however, the difference becomes apparent. The isophotes produced from RALM are less fragmented and visually more pleasant than LALM does. From the two zoomed-in areas of the lena image, one can easily check that for the results of LALM, there are image edges collapsing in the nose part, while for the proposed RALM, the isophotes are connected. The same phenomenon occur in the isophotes of the recovered peppers and cameraman images, shown in the last two rows of Fig.7. In the zoomed-in junction parts of the peppers, as well as in the zoomed-in parts of the cameraman’s face area, the proposed method produces smoother isophotes than LALM does.

To compare in a quantitative manner, the energy and PSNR evolutions for the two methods are monitored and shown in Fig.8. To fully demonstrate their behaviors, 1000 iterations are computed, which are beyond the need for convergence. From the plots, one can conclude that the energy and PSNR curves of RALM converge to steady state much faster than LALM does. Furthermore, compared with LALM, RALM achieves lower energy value and higher PSNR value, which suggests that the restored images by RALM possess higher quality. The quantitative indices shown in Table II further validate the above conclusion. one can check that RALM recovers higher quality images in terms of NRMSE, NMAD, and PSNR, while consumes much less iterations.

IV Conclusion

By considering the limiting case of $b\to 0$, an internal inconsistency is identified in the conventional augmented Lagrangian formulation for solving the Euler’s elastica model. Based on this key observation, a simple cutting-off strategy is introduced to the solution procedure of the $\boldsymbol{p}$ subproblem which helps to remove the inconsistency. Besides, our analysis shows that the commonly used relaxation technique for imposing the constraint for the auxiliary variable $\boldsymbol{n}$ is actually problematic, and its original form should be instead adopted. The above two observations lead to the proposed restricted augmented Lagrangian method (RALM), which enjoys easier parameter tuning, faster convergence, and better image restoration ability.

We think that the proposed **RALM ** provides a strategy for accelerating and improving other augmented Lagrangian based algorithms for solving the Euler’s elastica model or other models that need to introduce multiple auxiliary variables with ordered dependence.