Image Restoration Based on the Hybrid Total-Variation-Type Model

Abstract

We propose a hybrid total-variation-type model for the image restoration problem based on combining advantages of the ROF model with the LLT model. Since two 𝐿¹-norm terms in the proposed model make it difficultly solved by using some classically numerical methods directly, we first employ the alternating direction method of multipliers (ADMM) to solve a general form of the proposed model. Then, based on the ADMM and the Moreau-Yosida decomposition theory, a more efficient method called the proximal point method (PPM) is proposed and the convergence of the proposed method is proved. Some numerical results demonstrate the viability and efficiency of the proposed model and methods.

1. Introduction

Image restoration is of momentous significance in coherent imaging systems and various image processing applications. The goal is to recover the real image from the deteriorated image, for example, image denoising, image deblurring, image inpainting, and so forth; see [1–3] for details.

For the additive noisy image, many denoising models have been proposed based on PDEs or variational methods over the last decades [1–3]. The essential idea for this class of models is to filter out the noise in an image while preserving significant features such as edges and textures. However, due to the ill-posedness of the restoration problem, we have to employ some regularization methods [4] to overcome it. The general form of regularization methods consists in minimizing an energy functional of the following form: 𝐹(𝑢)=𝜆2𝑢−𝑢₀²_𝒳+ℛ(𝑢) (1.1) in Banach space 𝒳, where 𝜆 is the regularization parameter, ℛ, is a regularization term, 𝑢₀ is the observed image, and 𝑢 is the image to be restored. The development of the energy functional (1.1) actually profits from the ROF model [5] which is of the following form: min_𝑢𝜆₁2𝑢−𝑢₀²_𝐿2+|𝑢|_𝐵𝑉, (1.2) where |𝑢|_𝐵𝑉∶=∫_Ω|∇𝑢|d𝑥. In the case without confusion, for simplification we omit the open set Ω with Lipschitz boundary for ‖⋅‖_𝐿²_(Ω) and |⋅|_𝐿¹_(Ω). Due to edge-preserving property of the term |𝑢|_𝐵𝑉, this model has been extended to other sorts of image processing problems such as min_𝑢𝜆₂2𝐾𝑢−𝑢₀²_𝐿2+|𝑢|_𝐵𝑉 (1.3) for image deblurring [6] and min_𝑢𝜆₃2𝑢−𝑢₀²_{𝐿2(Ω⧵𝐷)}+|𝑢|_𝐵𝑉 (1.4) for image inpainting [2, 7], where 𝐾 is a blurring operator and 𝐷 is the inpainting domain. Furthermore, this model was applied to restore the multiplicative noisy image which usually appears in various image processing applications such as in laser images, ultrasound images [8], synthetic aperture radar (SAR) [9], and medical ultrasonic images [10]. One of the models based on 𝐵𝑉 was proposed by Huang et al. (HNW) [11] with the form min_𝑧𝛼_Ω𝑧+𝑓𝑒^−𝑧d𝑥+|𝑧|_𝐵𝑉, (1.5) where 𝑧 is the logarithmic transformation of 𝑢 and |𝑧|_𝐵𝑉 keeps the total variation property related to |𝑒^𝑧|_𝐵𝑉 and 𝛼>0. Here 𝑢 satisfies 𝑓=𝑢𝜂 for the noise 𝜂.

However, as we all know the |⋅|_𝐵𝑉 term usually reduces the computational solution of the above models to be piecewise constant, which is also called the staircasing effect in smooth regions of the image. The staircase effect implies to produce new edges that do not exist in the true image so that the restored image is unsatisfactory to the eye. To overcome this drawback, some high-order models [12–15] have been proposed such as a model proposed by Lysaker, Lundervold, and Tai (the LLT model) with the following form: min_𝑢𝛽₁2𝑢−𝑢₀²_𝐿2+|𝑢|_𝐵𝑉², (1.6) where |𝑢|_𝐵𝑉²∶=∫_Ω|∇²𝑢|d𝑥 for the Hessian operator ∇². However, these classes of models can blur the edges of the image in the course of restoration. Therefore, it is a natural choice to combine advantages of the ROF model and the LLT model if we want to preserve edges while avoiding the staircase effect in smooth regions. One of convex combinations between the 𝐵𝑉 and 𝐵𝑉² was proposed by Lysaker et al. [13] to restore the image with additive noise. But their model is not quite intuitive due to lack of gradient information in the weighting function. Since the edge detector function 𝑔∶=𝑔(𝑢₀)∶=1/(1+𝜚|∇(𝑢₀∗𝐺_𝜎(𝑥))|) with 𝐺_𝜎(𝑥)=(1/2𝜋𝜎²)exp(−(𝑥²/2𝜎²)) can depict the information of edges, we can employ it as a balance function; that is, we can apply the following model: min_𝑢𝛽2𝑢−𝑢₀²_𝐿2+1−𝑔𝑢_𝐵𝑉+𝑔𝑢_𝐵𝑉2 (1.7) to restore the noisy image [16]. Obviously, |(1−𝑔)𝑢|_𝐵𝑉 tends to be predominant where edges most likely appear and |𝑔𝑢|_𝐵𝑉² tends to be predominant at the locations with smooth signals. Based on the advantages of the hybrid model (1.7), we also extend it to the image restoration models (1.3)–(1.5) in this paper.

Another topic for image restoration is to find some efficient methods to solve the above proposed models. In fact, there are many different methods based on PDE or convex optimization to solve the minimization problem (1.1) by means of the specific models (1.2)–(1.6). For example, for the purpose of solving the ROF model (1.2) or the LLT model (1.6), we can use the gradient descent method [5, 13], the Chambolle’s dual method [13, 17, 18], the primal and dual method [19–21], the second order cone programming method [22], the multigrid method [23], operator splitting method [24–26], the inverse scale method [27], and so forth. However, different from the ROF model (1.2) and the LLT model (1.6), the model (1.7) includes two 𝐿¹-norm terms which make it solved more difficultly. More generally, the model (1.7) can be fell into the following framework: min_𝑢𝜃2𝑢−𝑓²_𝐿2+ℎ₁Λ₁𝑢+ℎ₂Λ₂𝑢, (1.8) where ℎ₁,ℎ₂∶𝒳→R are proper, convex, and lower semicontinuous (l.s.c.), Λ₁ and Λ₂ are bounded linear operators, and 𝜃 is a parameter. A specific form of (1.8) was considered by Afonso and Bioucas-Dias in [28] where a Bregman iterative method was proposed to solve the model with the combination of the 𝐿¹-norm term and the total variation (TV) term. Actually this splitting Bregman method is formally equivalent to the alternating direction method of multipliers (ADMM) [24, 29–34]. However, the ADMM ineluctably tends to solve some subproblems which correspond to the related modified problems. Furthermore, these make us obtain the numerical results by requiring much more computational cost. In order to obtain an efficient numerical method, it is a natural choice to avoid solving the related subproblems. In this paper, we propose a proximal point method (PPM) which can be deduced from the ADMM. This deduction is based on the connection that the sum of the projection operator and the shrinkage operator is equal to the identity operator; it is known as the Moreau-Yosida decomposition Theorem  31.5 in [35]. Then the PPM not only keeps the advantages of the ADMM but also requires much less computational cost. This implies that the PPM is much more fast and efficient, especially for the larger scale images. Furthermore, using the monotone operator theory, we give the convergence analysis of the proposed method. Moreover, we extend the PPM to solve the model (1.7) to image deblurring, image inpainting, and the multiplicative noisy image restoration. Experimental results show that the restored images generated by the proposed models and methods are desirable.

The paper is organized as follows. In Section 2, we recall some knowledge related to convex analysis. In Section 3, we first propose the ADMM to solve the problem (1.8) and then give the PPM to improve this method. In Section 4, we give some applications by using the proposed algorithms and also compare the related models and the proposed methods. Some concluding remarks are given in Section 5.

2. Notations and Definitions

Let us describe some notations and definitions used in this paper. For simplifying, we use 𝑋=𝑅^𝑛. Usually, we can set 𝑛=2 for the gray-scale images. The related contents can be referred to [1, 27, 36–43].

Definition 2.1. The operator 𝐴∶𝑋→𝑅 is monotone if it satisfies 𝑦₁−𝑦₂,𝑥₁−𝑥₂≥0 (2.1) for all 𝑦₁∈𝐴(𝑥₁) and 𝑦₂∈𝐴(𝑥₂), and 𝐴 is maximal if it is not strictly contained in any other monotone operator on 𝑋.

Definition 2.2. Let 𝐺∶𝑋→𝑅∪{+∞} be a convex function. The subdifferential of 𝐺 at a point 𝑦∈𝑋 is defined by 𝜕𝐺𝑦=𝜂∶𝐺𝑦−𝐺(𝑥)≤𝜂,𝑦−𝑥,∀𝑥∈𝑋, (2.2) where 𝜂∈𝜕𝐺(𝑦) is called a subgradient. It reduces to the classical gradient if 𝐺(𝑥) is differentiable.

Definition 2.3. Assume that 𝐴 is a maximal monotone operator. Denote by 𝐻_𝑐 its Yosida approximation: 𝐻_𝑐𝐴=1𝑐𝐼−𝐽_𝑐𝐴, (2.3) where 𝐼 denotes the identity operator [38, 39] and 𝐽_𝑐𝐴 is the resolvent of 𝐴 with the form as 𝐽_𝑐𝐴=(𝐼+𝑐𝐴)⁻¹. (2.4)

Definition 2.4. Let 𝐺∶𝑋→𝑅 be a convex proper function. The conjugate function of 𝐺 is defined by 𝐺^∗𝑥^∗=sup_𝑥∈𝑋𝑥,𝑥^∗−𝐺(𝑥) (2.5) for all 𝑥^∗∈𝑋.

Definition 2.5. Let 𝐺∶𝑋→𝑅∪{∞} be convex and 𝑡>0. The proximal mapping to 𝐺∶𝑋→𝑋 is defined by Prox_𝑡𝐺(𝑥)∶=argmin_𝑦𝐺𝑦+12𝑡𝑦−𝑥²_𝐿2 (2.6) for 𝑦∈𝑋.

Definition 2.6. The projection operator 𝒫_ℬ𝜏(⋅)∶𝑋→𝑋 onto the closed disc ℬ_𝜏∶={𝑥∈𝑋∶|𝑥|_𝐿¹≤1/𝜏} is defined by 𝒫_ℬ𝜏(𝑥)=𝑥|𝑥|_𝐿¹min|𝑥|_𝐿¹,1𝜏, (2.7) where 𝑥∈𝑋 and 𝜏>0.

Definition 2.7. The shrinkage operator 𝑆_𝜏(⋅)∶𝑋→𝑋 is defined by 𝒮_𝜏(𝑥)=𝑥|𝑥|_𝐿¹max|𝑥|_𝐿¹−1𝜏,0, (2.8) where we use the convention 0/0=0.

Remark 2.8. It is obvious that the function 𝐺∶𝑋→𝑅 and its conjugate function 𝐺^∗ satisfy the following relationship: Prox_𝑡𝐺(𝑥)+Prox_1/𝑡𝐺^∗(𝑥)=𝑥, (2.9) for 𝑡>0. Especially, the projection operator 𝒫_ℬ𝜏(⋅) and the shrinkage operator 𝑆_𝜏(⋅) satisfy 𝒫_ℬ𝜏(𝑥)+𝑆_𝜏(𝑥)=𝑥 (2.10) for any 𝑥∈𝑋. In fact, this corresponds to the classic Moreau-Yosida decomposition Theorem  31.5 in [35].

3. The Alternating Direction Method of Multipliers (ADMM) and the Proximal Point Method (PPM)

Variable splitting methods such as the ADMM [29–32, 44] and the operator splitting methods [42, 45, 46] have been recently used in the image, signal, and data processing community. The key of this class of methods is to transform the original problem into some subproblems so that we can easily solve these subproblems by employing some numerical methods. In this section we first consider to use the ADMM to solve the general minimization problem (1.8). However, the computational cost of the ADMM is tediously increased due to its looser form. In order to overcome this drawback, we thus change this method to a compacter form called the proximal method based on the relationship (2.9) in Remark 2.8.

3.1. The Alternating Direction Method of Multipliers (ADMM)

We now consider the following constrained problem: min_{𝑢,𝑣,𝑤}𝜃2𝑢−𝑓²_𝐿2+ℎ₁(𝑣)+ℎ₂(𝑧),s.t.𝑣=Λ₁𝑢,𝑧=Λ₂𝑢, (3.1) which is clearly equivalent to the constrained problem (1.8) in the feasible set {(𝑢,𝑣,𝑧)∶𝑣=Λ₁𝑢,and𝑧=Λ₂𝑢}. Throughout the following subsections, we always assume that Λ₁ and Λ₂ are a surjective map. It seems that the problem (3.1) including three variables looks more complex than the original unconstrained problem (1.8). In fact, this problem can be solved more easily under the condition that ℎ₁ and ℎ₂ are nondifferentiable. In the augmented Lagrangian framework [33, 34, 47, 48], the problem (3.1) equivalently solves the following minimization Lagrangian function: min_{𝑢,𝑣,𝑧,𝜁1,𝜁2,𝜇1,𝜇2}𝐿𝑢,𝑣,𝑧,𝜁₁,𝜁₂,𝜇₁,𝜇₂=𝜃2𝑢−𝑓²_𝐿2+ℎ₁(𝑣)+𝜁₁,Λ₁+𝜇𝑢−𝑣¹2Λ₁𝑢−𝑣²_𝐿2+ℎ₂(𝑧)+𝜁₂,Λ₂+𝜇𝑢−𝑧²2Λ₂𝑢−𝑧²_𝐿2, (3.2) where 𝜁_𝑖 is the Lagrangian multiplier and 𝜇_𝑖 is the penalty parameter for 𝑖=1, 2. Then we can use the following augmented Lagrangian method (ALM): 𝑢^𝑛,𝑣^𝑛,𝑧^𝑛=argmin_{𝑢,𝑣,𝑧}𝐿𝑢,𝑣,𝑧,𝜁₁^𝑛−1,𝜁₂^𝑛−1,𝜇₁,𝜇₂,𝜁^𝑛₁=𝜁₁^𝑛−1+𝜇₁Λ₁,𝜁𝑢−𝑣_𝑛2=𝜁₂^𝑛−1+𝜇₂Λ₂𝑢−𝑧 (3.3a)(3.3b)(3.3c)with choosing the original values 𝜁⁰₁ and 𝜁⁰₂ to solve (3.2). If we set 𝑑^𝑛_𝑖=𝜁^𝑛_𝑖/𝜇_𝑖 for 𝑖=1, 2 and omit the terms which are independent of (𝑢^𝑛,𝑣^𝑛,𝑧^𝑛) in (3.3a), the above strategy (3.3a)–(3.3c) can be written as 𝑢^𝑛,𝑣^𝑛,𝑧^𝑛=argmin_{𝑢,𝑣,𝑧}𝜃2𝑢−𝑓²_𝐿2+ℎ₁(𝑣)+𝜇₁2Λ₁𝑢−𝑣+𝑑₁^𝑛−12_𝐿2+ℎ₂(𝑧),+𝜇₂2Λ₂𝑢−𝑧+𝑑₂^𝑛−12_𝐿2,𝑑^𝑛₁=𝑑₁^𝑛−1+Λ₁𝑢^𝑛−𝑣^𝑛,𝑑^𝑛₂=𝑑₂^𝑛−1+Λ₂𝑢^𝑛−𝑧^𝑛 (3.4a)(3.4b)(3.4c)for the original values 𝑑⁰₁ and 𝑑⁰₂. Then we can use the following ADMM to solve (3.4a)–(3.4c).

Algorithm 3.1 (ADMM for solving (3.4a)–(3.4c)). (1) Choose the original values: 𝑣⁰, 𝑧⁰, 𝑑⁰₁, and 𝑑⁰₂. Set 𝜃, 𝜇₁, 𝜇₂>0 and 𝑛=1.
(2) Compute (𝑢^𝑛,𝑣^𝑛,𝑧^𝑛,𝑑^𝑛₁,𝑑^𝑛₂) by

𝑢^𝑛=argmin_𝑢𝜃2𝑢−𝑓²_𝐿2+𝜇₁2Λ₁𝑢−𝑣^𝑛−1+𝑑₁^𝑛−12_𝐿2+𝜇₂2Λ₂𝑢−𝑧^𝑛−1+𝑑₂^𝑛−12_{𝐿2=ℎ(𝑢)},𝑣^𝑛=argmin_𝑣ℎ₁(𝑣)+𝜇₁2Λ₁𝑢^𝑛−𝑣+𝑑₁^𝑛−12_𝐿2,𝑧^𝑛=argmin_𝑧ℎ₂(𝑧)+𝜇₂2Λ₂𝑢^𝑛−𝑧+𝑑₂^𝑛−12_𝐿2,𝑑^𝑛₁=𝑑₁^𝑛−1+Λ₁𝑢^𝑛−𝑣^𝑛,𝑑^𝑛₂=𝑑₂^𝑛−1+Λ₂𝑢^𝑛−𝑧^𝑛. (3.5a)(3.5b)(3.5c)(3.5d)(3.5e)

(3) If the stop criterion is not satisfied, set 𝑛∶=𝑛+1 and go to step (2).

Since ℎ(𝑢) is differentiable and strictly convex, we can get the unique solution of (3.5a) which satisfies 𝜃𝑢^𝑛−𝑓+𝜇₁Λ^∗₁Λ₁𝑢^𝑛−𝑣^𝑛−1+𝑑₁^𝑛−1+𝜇₂Λ^∗₂Λ₂𝑢^𝑛−𝑧^𝑛−1+𝑑₂^𝑛−1=0,0=𝑤^𝑛₁+𝜇₁𝑣^𝑛−𝑑₁^𝑛−1−Λ₁𝑢^𝑛,0=𝑤^𝑛₂+𝜇₂𝑧^𝑛−𝑑₂^𝑛−1−Λ₂𝑢^𝑛𝑑_𝑛1,=𝑑₁^𝑛−1+Λ₁𝑢^𝑛−𝑣^𝑛,𝑑^𝑛₂=𝑑₂^𝑛−1+Λ₂𝑢^𝑛−𝑧^𝑛, (3.6a)(3.6b)(3.6c)(3.6d)(3.6e)where Λ^∗₁ and Λ^∗₂ are the adjoint operators of Λ₁ and Λ₂, 𝑤^𝑛₁∈𝜕ℎ₁(𝑣^𝑛) and 𝑤^𝑛₂∈𝜕ℎ₂(𝑧^𝑛), respectively. It follows that the solution 𝑢^𝑛 in (3.6a) can be directly obtained by the following explicit formulation: 𝑢^𝑛=𝜃𝐼+𝜇₁Λ^∗₁Λ₁+𝜇₂Λ^∗₂Λ_2𝑀⁻¹𝜃𝑓+𝜇₁Λ^∗₁𝑣^𝑛−1−𝑑₁^𝑛−1+𝜇₂Λ^∗₂𝑧^𝑛−1−𝑑₂^𝑛−1. (3.7) However, when the operator 𝑀 is ill-posed, the solution is unsuitable or unavailable. Hence we have to go back to (3.6a) and to employ some iteration strategy such as the Gauss-Seidel method to solve this equation. On the other hand, it is obvious that (3.5b) and (3.5c) can be looked at as the proximal mapping, so the solutions of the minimization problems (3.5b) and (3.5c) can be obviously written as 𝑣^𝑛=Prox_1/𝜇1ℎ₁Λ₁𝑢^𝑛+𝑑₁^𝑛−1,𝑧^𝑛=Prox_1/𝜇2ℎ₂Λ₂𝑢^𝑛+𝑑₂^𝑛−1. (3.8)

Theorem 3.2. Assume that (𝑢^∗,𝑣^∗,𝑧^∗,𝑤^∗₁,𝑤^∗₂) is the saddle point of the Lagrange function ℒ𝑢,𝑣,𝑧,𝑤₁,𝑤₂=𝜃2𝑢−𝑓²_𝐿2+ℎ₁(𝑣)+𝑤₁,Λ₁𝑢−𝑣+ℎ₂(𝑧)+𝑤₂,Λ₂𝑢−𝑧. (3.9) Then 𝑢^∗ is the solution of the minimization problem (1.8). Furthermore, the sequence {(𝑢^𝑛,𝑣^𝑛,𝑧^𝑛)} generated by Algorithm 3.1 converges to (𝑢^∗,𝑣^∗,𝑧^∗).

Notice that Algorithm 3.1 can be actually looked at as the split Bregman method [25]. The based idea of this method is to introduce some intermediate variables so as to transform the original problem into some subproblems which are easily solved. The connection between the split Bregman method and the ADMM has been shown in [29, 49]. However, our algorithm considers the sum of three convex functions, which is more general than the related algorithms in [25, 49]. Furthermore, it must be noted that 𝑣 and 𝑧 are completely separated in (3.4a), so the two subproblems (3.5b) and (3.5c) are parallel. Therefore the convergence results of the ADMM can be applied here.

3.2. The Proximal Point Method

Though the ADMM in Algorithm 3.1 can effectively solve the original problem (3.1), we have to solve five subproblems. This actually makes this method suffer from a looser form as in [25, 45] so that it can badly affect its numerical computation efficiency. In this subsection, we propose a compacter form comparing to the ADMM. This formation called the PPM by using the relationship (2.9) in Remark 2.8 can reduce the original five subproblems of the ADMM in Algorithm 3.1 to solve three subproblems, thus it can improve computation cost of the ADMM. Now we have rewritten (3.5a)–(3.5e) with a little variation as the following form: 𝑣^𝑛=argmin_𝑣ℎ₁(𝑣)+𝜇₁2𝑣−𝑑₁^𝑛−1−Λ₁𝑢^𝑛−12_𝐿2,𝑧^𝑛=argmin_𝑧ℎ₂(𝑧)+𝜇₂2𝑧−𝑑₂^𝑛−1−Λ₂𝑢^𝑛−12_𝐿2,𝑑^𝑛₁=𝑑₁^𝑛−1+Λ₁𝑢^𝑛−1−𝑣^𝑛,𝑑^𝑛₂=𝑑₂^𝑛−1+Λ₂𝑢^𝑛−1−𝑧^𝑛,𝑢^𝑛=argmin_𝑢𝜃2𝑢−𝑓²_𝐿2+𝜇₁2Λ₁𝑢−𝑣^𝑛+𝑑^𝑛₁²_𝐿2+𝜇₂2Λ₂𝑢−𝑧^𝑛+𝑑^𝑛₂²_𝐿2 (3.10a)(3.10b)(3.10c)(3.10d)(3.10e)with the first order optimality conditions given by 𝑣^𝑛=prox_1/𝜇1ℎ₁Λ₁𝑢^𝑛−1+𝑑₁^𝑛−1,𝑧^𝑛=prox_1/𝜇1ℎ₂Λ₂𝑢^𝑛−1+𝑑₂^𝑛−1,𝑑^𝑛₁=𝑑₁^𝑛−1+Λ₁𝑢^𝑛−1−𝑣^𝑛,𝑑^𝑛₂=𝑑₂^𝑛−1+Λ₂𝑢^𝑛−1−𝑧^𝑛,𝜃𝑢^𝑛−𝑓+𝜇₁Λ^∗₁Λ₁𝑢^𝑛−𝑣^𝑛+𝑑^𝑛₁+𝜇₂Λ^∗₂Λ₂𝑢^𝑛−𝑧^𝑛+𝑑^𝑛₂=0. (3.11a)(3.11b)(3.11c)(3.11d)(3.11e)If (3.11e) is replaced by 𝜃𝑢^𝑛−𝑓+𝜇₁Λ^∗₁Λ₁𝑢^𝑛−1−𝑣^𝑛+𝑑₁^𝑛−1+𝜇₂Λ^∗₂Λ₂𝑢^𝑛−1−𝑧^𝑛+𝑑₂^𝑛−1=0, (3.12) it follows from (3.11a)–(3.11e) and Moreau-Yosida decomposition Theorem  31.5 in [35] that 𝑑^𝑛₁=prox_𝜇1ℎ^∗₁Λ₁𝑢^𝑛−1+𝑑₁^𝑛−1,𝑑^𝑛₂=prox_𝜇2ℎ^∗₂Λ₂𝑢^𝑛−1+𝑑₂^𝑛−1,𝑢^𝑛𝜇=𝑓−¹𝜃Λ^∗₁𝑑^𝑛₁+𝜇₂𝜃Λ^∗₂𝑑^𝑛₂. (3.13) So we propose the following algorithm to solve (3.1).

Algorithm 3.3 (PPM for solving (3.1)). (1) Choose the original values: 𝑑⁰₁, 𝑑⁰₂, and 𝑢⁰. Set 𝜃, 𝜇₁, 𝜇₂>0 and 𝑛=1.
(2) Compute (𝑑^𝑛₁,𝑑^𝑛₂,𝑢^𝑛) by (3.13).
(3) If the stop criterion is not satisfied, set 𝑛∶=𝑛+1 and go to step (2).

Lemma 3.4. Set 𝑥₁, 𝑥₂∈𝑋 and 𝐴 is a maximal monotone operator, then the operators 𝐻_𝑐𝐴 and 𝐽_𝑐𝐴 satisfy 1𝑐²𝐽_𝑐𝐴(𝑥₁)−𝐽_𝑐𝐴(𝑥₂)²_𝐿2+𝐻_𝑐𝐴(𝑥₁)−𝐻_𝑐𝐴(𝑥₂)²_𝐿2≤1𝑐²𝑥₁−𝑥₂²_𝐿2. (3.14)

Theorem 3.5. Assume that ℎ₁(𝑥) and ℎ₂(𝑢) are convex and proper. If 𝜇₁∈(0,1/𝜃‖Λ₁‖²) and 𝜇₂∈(0,1/𝜃‖Λ₂‖²), here ‖⋅‖∶=max{‖𝐾𝑥‖_𝐿²∶𝑥∈𝑋 with ‖𝑥‖_𝐿²≤1} for a continuous linear operator 𝐾, then the sequence {(𝑢^𝑛,𝑑^𝑛₁,𝑑^𝑛₂)} generated by Algorithm 3.3 converges to the limit point (𝑑₁,𝑑₂,𝑢). Furthermore, the limit point 𝑢 is the solution of (1.8).

4. Some Applications in Image Restoration

In Section 4.1, we apply the ADMM and the above PPM to the image denoising problem. Here we also compare the proposed hybrid model with the ROF model and the LLT model. Then, based on the proposed hybrid model, we set the PPM as a basic method to solve image deblurring, image inpainting, and image denoising for the multiplicative noise in the last three subsections. For simplicity, we assume that the image region Ω is squared with the size 𝑀×𝑀 and set 𝑆=𝑅^𝑀×𝑀, 𝑇=𝑆×𝑆, and 𝑍=𝑇×𝑇 as in [17]. The usual scalar product can be denoted as ⟨𝑝¹,𝑝²⟩_𝑇∶=∑^𝑀_𝑖=1∑^𝑀_𝑗=1𝑝¹_𝑖,𝑗𝑝²_𝑖,𝑗 for 𝑝¹,𝑝²∈𝑇 and (𝐩,𝐪)_𝑍=∑^𝑀_𝑖,𝑗=1(𝑝¹¹_𝑖,𝑗𝑞¹¹_𝑖,𝑗+𝑝¹²_𝑖,𝑗𝑞¹²_𝑖,𝑗+𝑝²¹_𝑖,𝑗𝑞²¹_𝑖,𝑗+𝑝²²_𝑖,𝑗𝑞²²_𝑖,𝑗) for 𝐩,𝐪∈𝑍. The 𝐿¹ norm of 𝐩=(𝑝₁,𝑝₂)∈𝑇 is defined by |𝐩|_ℓ¹=𝑝²₁+𝑝²₂ and the ℓ¹ norm of 𝐪=(𝑞₁,𝑞₂;𝑞₃,𝑞₄)∈𝑍 is defined by |𝐪|_ℓ¹=𝑞²₁+𝑞²₂+𝑞²₃+𝑞²₄. If 𝑢∈𝑆, we use ∇⁺=(∇⁺_𝑥,∇⁺_𝑦)∈𝑇 to denote the first order forward difference operator with ∇⁺_𝑥𝑢_𝑖,𝑗=𝑢_𝑖+1,𝑗−𝑢_𝑖,𝑗∇for1≤𝑖<𝑁,0for𝑖=𝑁,^+𝑦𝑢_𝑖,𝑗=𝑢_𝑖,𝑗+1−𝑢_𝑖,𝑗for1≤𝑗<𝑁,0for𝑗=𝑁, (4.1) and use ∇⁻=(∇⁻_𝑥,∇⁻_𝑦) to denote the first order backward difference operator with ∇⁻_𝑥𝑢_𝑖,𝑗=−𝑢_𝑖,𝑗𝑢,for𝑖=1,_𝑖,𝑗−𝑢_{𝑖−1,𝑗}𝑢for1<𝑖<𝑁,∇𝑖,𝑗for𝑖=𝑁,^−𝑦𝑢_𝑖,𝑗=−𝑢_𝑖,𝑗𝑢for𝑗=1,_𝑖,𝑗−𝑢_{𝑖,𝑗−1}𝑢for1<𝑗<𝑁,_𝑖,𝑗for𝑗=𝑁, (4.2) for 𝑖,𝑗=1,…,𝑁. Based on the first order difference operators, we can give the second order difference operator as follows: ∇²𝑢_𝑖,𝑗=∇⁻_𝑥∇⁺_𝑥𝑢_𝑖,𝑗∇⁺_𝑥∇⁺_𝑦𝑢_𝑖,𝑗∇⁻_𝑦∇⁻_𝑥𝑢_𝑖,𝑗∇⁻_𝑦∇⁺_𝑦𝑢_𝑖,𝑗∈𝑍. (4.3) Using the same approach, we can define some other second order operators such as ∇⁺_𝑥(∇⁻_𝑥𝑢_𝑖,𝑗), ∇⁻_𝑥(∇⁻_𝑦𝑢_𝑖,𝑗), ∇⁺_𝑦(∇⁺_𝑥𝑢_𝑖,𝑗), and ∇⁺_𝑦(∇⁻_𝑦𝑢_𝑖,𝑗). Then we can give the first order and the second divergence operators as div𝑢_𝑖,𝑗=∇⁻_𝑥𝑢_𝑖,𝑗+∇⁻_𝑦𝑢_𝑖,𝑗,div²𝑢_𝑖,𝑗=∇⁺_𝑥∇⁻_𝑥𝑢_𝑖,𝑗+∇⁻_𝑥∇⁻_𝑦𝑢_𝑖,𝑗+∇⁺_𝑦∇⁺_𝑥𝑢_𝑖,𝑗+∇⁺_𝑦∇⁻_𝑦𝑢_𝑖,𝑗. (4.4) Furthermore, if we set 𝐩∈𝑇 and 𝐪∈𝑍, it is easy to deduce that div𝑝²₂≤8‖𝐩‖²_𝑇,div²𝑞²₂≤64‖𝐪‖²_𝑍. (4.5)

Remark 4.1. The related examples in the following subsections are performed using Windows 7 and Matlab 2009(a) on a desktop with Intel Core i5 processor at 2.4 GHz and 4 GB memory. All of the parameters for related models are chosen by trial and empirically which can yield better restored images. On the other hand, we should notice that it is not very expensive when we use the ADMM and the PPM to get 𝑢^𝑛, but the total computational effort of one outer iteration requiring many inner steps can be very huge. In order to reduce the computational effort and keep fair comparison of these two methods, we so simplify the inner-outer iterative framework by performing only one-step in inner iteration. It is obvious that these sets are very efficient from the following numerical experiences.

4.1. Image Denoising for the Additive Noise

In this subsection, we consider to use the ADMM and the PPM to solve (1.7) for restoring the additive noisy image. If we set Λ₁=∇ and Λ₂=∇², then the algorithms are proposed as follows.

Algorithm 4.2 (ADMM to solve (1.7)). (1) Choose the original 𝑑⁰₁=𝑣⁰=𝟎∈𝑇 and 𝑑⁰₂=𝑧⁰=𝟎∈𝑍. Set 𝜃, 𝜇₁, 𝜇₂>0 and 𝑛=1.
(2) Compute (𝑢^𝑛,𝑣^𝑛,𝑧^𝑛,𝑑^𝑛₁,𝑑^𝑛₂) by 𝜃𝐼−𝜇₁Δ+𝜇₂Δ²𝑢^𝑛=𝜃𝑓+𝜇₁div𝑣^𝑛−1−𝑑₁^𝑛−1−𝜇₂div²𝑧^𝑛−1−𝑑₂^𝑛−1𝑣^𝑛=∇𝑢^𝑛+𝑑₁^𝑛−1∇𝑢^𝑛+𝑑₁^𝑛−1_ℓ1⋅max∇𝑢^𝑛+𝑑₁^𝑛−1_ℓ1−1−𝑔(𝑥)𝜇₁,𝑧,0_𝑛=∇²𝑢^𝑛+𝑑₂^𝑛−1∇²𝑢^𝑛+𝑑₂^𝑛−1_ℓ1⋅max∇²𝑢^𝑛+𝑑₂^𝑛−1_ℓ1−𝑔(𝑥)𝜇₂,𝑑,0_𝑛1=𝑑₁^𝑛−1+∇𝑢^𝑛−𝑣^𝑛,𝑑^𝑛₂=𝑑₂^𝑛−1+∇²𝑢^𝑛−𝑧^𝑛, (4.6a)(4.6b)(4.6c)(4.6d)(4.6e)where Δ=div∘∇ and Δ²=div²∘∇².
(3) If the stop criterion is not satisfied, set 𝑛∶=𝑛+1 and go to step (2).

Remark 4.3. For the first subproblem (4.6a), we can use the Gauss-Seidel method as shown in [25] to get the solution. However, in this paper, we use the following strategy: 𝑢^𝑛_𝑖,𝑗=𝜃𝑓_𝑖,𝑗+𝜇₁div𝑣^𝑛−1_𝑖,𝑗−𝑑^𝑛−1_1,𝑖,𝑗+4𝜇₁𝑢^𝑛−1_𝑖,𝑗−𝜇₂div²𝑧^𝑛−1_𝑖,𝑗−𝑑^𝑛−1_2,𝑖,𝑗+28𝜇₂𝑢^𝑛−1_𝑖,𝑗𝜃+4𝜇₁+28𝜇₂, (4.7) where some information of operators Δ and Δ² related to 𝑢 is used. The formulas (4.6b) and (4.6c) of Algorithm 4.2 can be easily deduced from min_𝑣1−𝑔𝑣_ℓ1+𝜇₁2𝑣−∇𝑢^𝑛−𝑑₁^𝑛−12₂,min_𝑧𝑔𝑣_ℓ1+𝜇₂2𝑧−∇²𝑢^𝑛−𝑑₂^𝑛−12₂. (4.8) Furthermore, following form Theorem 3.2, we can also deduce that the sequence {𝑢^𝑛} generated by Algorithm 4.2 converges to the solution of (1.7).

Algorithm 4.4 (PPM to solve (1.7)). (1) Choose the original 𝑑⁰₁=𝟎∈𝑆, 𝑑⁰₂=𝟎∈𝑍 and 𝑢⁰=𝑓. Set 𝜃, 𝜇₁, 𝜇₂>0 and 𝑛=1.
(2) Compute (𝑑^𝑛₁,𝑑^𝑛₂,𝑢^𝑛) by 𝑑^𝑛₁=∇𝑢^𝑛−1+𝑑₁^𝑛−1∇𝑢^𝑛−1+𝑑₁^𝑛−1_ℓ1⋅min∇𝑢^𝑛−1+𝑑₁^𝑛−1_ℓ1,𝑑,1_𝑛2=∇²𝑢^𝑛−1+𝑑₂^𝑛−1∇²𝑢^𝑛−1+𝑑₂^𝑛−1_ℓ1⋅min∇²𝑢^𝑛−1+𝑑₂^𝑛−1_ℓ1,𝑢,1_𝑛𝜇=𝑓−¹𝜃1−𝑔div𝑑^𝑛₁+𝜇₂𝜃𝑔div²𝑑^𝑛₂. (4.9)
(3) If the stop criterion is not satisfied, set 𝑛∶=𝑛+1 and go to step (2).

For Algorithm 4.4, based on the relations in (4.5) and Theorem 3.5, we have the following result.

Theorem 4.5. If 𝜇₁∈(0,1/8𝜃) and 𝜇₂∈(0,1/64𝜃), then the sequence {𝑢^𝑛} generated by Algorithm 4.4 converges to the solution of (1.7).

Remark 4.6. For the above two algorithms, we can also set 𝑔 as a constant between 0 and 1. In fact, it is easy to find that the algorithms correspond to solving the ROF model or the LLT model, respectively, when 𝑔=0 or 𝑔=1. At this time, the iteration strategy can be simplified as 𝑢^𝑛_𝑖,𝑗=𝜃𝑓_𝑖,𝑗+𝜇₁div𝑣^𝑛−1_𝑖,𝑗−𝑑^𝑛−1_1,𝑖,𝑗+4𝜇₁𝑢^𝑛−1_𝑖,𝑗𝜃+4𝜇₁,𝑢^𝑛_𝑖,𝑗=𝜃𝑓_𝑖,𝑗−𝜇₂div²𝑧^𝑛−1_𝑖,𝑗−𝑑^𝑛−1_2,𝑖,𝑗+28𝜇₂𝑢^𝑛−1_𝑖,𝑗𝜃+28𝜇₂ (4.10) for these two models, respectively. Furthermore, when 𝑔∈(0,1), these two algorithms correspond to solve the model which is the convex combination of the ROF model and the LLT model.

Example 4.7. In this example, we compare the ADMM with the PPM for solving the ROF model (1.4), the LLT model (1.6), and the hybrid model (1.7). The original images with three different sizes shown in Figure 1 are added to the Gaussian white noise with the standard deviation 𝜎=15.3. Iterations were terminated when the stop conditions ‖𝑢^𝑛+1−𝑢^𝑛‖₂/‖𝑢^𝑛‖₂≤𝜀 are met. It is easy to find the related results from Table 1 that the PPM is faster than the ADMM. Especially, the average CPU time of the PPM compared with that of the ADMM can save about 50% for the ROF model and the LLT model. It saves about 40% for the hybrid model.

Table 1 

The related results in Example 4.7.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1 

The original images and the noisy images with three different sizes in Example 4.7.

Example 4.8. In this example, the noisy image is added to the Gaussian white noisy with the standard deviation 𝜎=12. The algorithms will be stopped after 100 iterations. We compare the results generated by the ROF model, the LLT model, the convex combination of the ROF model and the hybrid model. As we can see from Table 2, the hybrid model get the lowest MSE and the highest SNR; these imply that the hybrid model can give the best restored image. On the other hand, it is easy to find that the ROF model makes staircasing effect appear and the LLT model leads to edge blurring. In fact, they are based on the fact that the restored model by the ROF model is piecewise constant on large areas and the LLT model as a higher model damps oscillations much faster in the region of edges. For the convex combined model and the hybrid model, they can efficiently suppress these two drawbacks. Furthermore, the hybrid model is more efficient than the convex combined model, because the hybrid model uses the edge detector function which can efficiently coordinate edge information. It should be noticed that here we use the Chambolle’s strategy [17] to solve the convex combined model so that it is slower than the strategy by using the PPM to solve the hybrid model. To focus on these facts, we present some zoomed-in local results and select a slice of the images which meets contours and the smooth regions shown in Figures 2 and 3.

Table 2 

The related data in Example 4.8 (here R.P. is the regularization parameter).

(a)

(b)

Figure 2 

The original and the noisy image in Example 4.8.

(a)

(b)

(c)

(d)

Figure 3 

The related restored images in Example 4.8.

4.2. Other Applications

In this subsection, we extend the hybrid model to other classes of restoration problems. As we can see in Section 4.1, the hybrid model has some advantages compared with the ROF model and the LLT model. Since the PPM is faster and more efficient than the ADMM, we only employ the PPM to solve the related image restoration model. Now we first consider the following iteration strategy: 𝑧^𝑛=ℱ𝑢^𝑛=∶argmin_𝑧𝛽2𝑢^𝑛−𝑧²_𝐿2(𝑧),𝑢+ℎ_𝑛+1=𝒢𝑧^𝑛∶=argmin_𝑢𝛽2𝑢−𝑧^𝑛2_𝐿2+_Ω1−𝑔∇𝑢d𝑥+_Ω𝑔∇²𝑢d𝑥_𝒟(𝑢), (4.11a)(4.11b)where ℎ(𝑧)∈𝐶¹(Ω). It is easy to find that the minimization problem 𝒟(𝑢) in (4.11b) is coercive and strictly convex, so the subproblem (4.11b) has a unique solution. Based on [50, Lemma 2.4], we also deduce that the operator 𝒢 is (1/2)-averaged nonexpansive.

Theorem 4.9. Assume that the functional ℳ()=𝛽𝑢,𝑧2‖𝑢−𝑧‖²_𝐿2+ℎ(𝑧)+_Ω1−𝑔∇𝑢d𝑥+_Ω𝑔∇²𝑢d𝑥 (4.12) is convex, coercive, and bounded below, then the sequence {(𝑢^𝑛,𝑧^𝑛)} generated by (4.11a)-(4.11b) converges to a point (𝑢^∗,𝑧^∗). Furthermore, when 𝑢^∗=𝑧^∗, 𝑢^∗ is the solution of min_𝑢ℎ(𝑢)+_Ω1−𝑔∇𝑢d𝑥+_Ω𝑔∇²𝑢d𝑥. (4.13)

Remark 4.10. It should be noticed that the following three models satisfy the conditions of Theorem 4.9, so we can use the strategy (4.11a)-(4.11b) to solve these models.

4.2.1. Image Deblurring

Now we apply the hybrid model to the image deblurring problems with the following formula: min_𝑢𝜆2𝐾𝑢−𝑢₀²_𝐿2+|𝑢|_{𝐵𝑉(1−𝑔)}+|𝑢|_𝐵𝑉²_𝑔. (4.14) Because of the compactness, the blurring operator 𝐾 is not continuously invertible, which implies that we cannot directly use the proximal method to solve minimization problem (4.14). However, based on the penalty method [48], we can transform the problem (4.14) to solve the following problem: 𝑧^𝑛=𝜇𝜆𝐾^𝑇𝐾+𝐼⁻¹𝜇𝜆𝐾^𝑇𝑓+𝑢^𝑛,𝑢^𝑛+1∶=argmin_𝑢12𝜇𝑢−𝑧^𝑛2_𝐿2+_Ω1−𝑔∇𝑢dx+_Ω𝑔∇²𝑢d𝑥. (4.15a)(4.15b)The key of (4.15a)-(4.15b) is to separate the blurring operator 𝐾 from the original problem (4.14) so that we can avoid the ill-posed operator 𝐾. Furthermore, it is obvious that the problem (4.15b) satisfie Theorem 4.9, so we can use the proximal method to solve it.

Example 4.11. In this experiment, we use the image Lena, which is blurred with a Gaussian kernel of “hsize = 3” and in which is added the Gaussian white noise with the standard deviation 𝜎=0.02. The related images and SNRs are arranged in Figure 4. As we can see in Example 4.7, the proximal method tends to be stable before 100 iterations, so here we fix the algorithm to be stopped when the iteration attains 100. It is easy to find that the hybrid model can give a better image quality than the other two models. Especially, we also observe some staircasing effect in Figure 4(b) and edge blurring in Figure 4(c). However, all of the drawbacks can be efficiently suppressed in Figure 1(d).

(a) Original image

(b) ROF model

(c) LLT model

(d) Hybrid model

Figure 4 

The related SNRs in Example 4.11. (a) 11.2166; (b) 14.3378; (c) 14.4778; (d) 14.5118.

4.2.2. Image Inpainting

Here we consider to use the hybrid model for the image inpainting problem with the following form: min_𝑢𝜆2𝑢−𝑢₀²_{𝐿2(Ω⧵𝐷)}+1−𝑔𝑢_𝐵𝑉+𝑔𝑢_𝐵𝑉2, (4.16) where 𝐷 is the inpainting domain and 𝜆 satisfies 𝜆=𝜆,if𝑢∈Ω⧵𝐷,0,if𝑢∈𝐷. (4.17) If we introduce an auxiliary variable 𝑧, based on the penalty method [48] again, then the solution of the problem (4.15a)-(4.15b) can be approximated by solving the following problem: 𝑢^𝑛𝑧∶=^𝑛+𝜆𝜇𝑢₀1+𝑧𝜆𝜇_𝑛+1,∶=argmin_𝑧12𝜇𝑧−𝑢^𝑛2_𝐿2+_Ω1−𝑔∇𝑧d𝑥+_Ω𝑔∇²𝑧d𝑥. (4.18a)(4.18b)

Example 4.12. In this example, we show the results of real image inpainting in Figures 6 and 7. We consider to use these three models to restore images shown in Figures 5(c) and 5(d), which are, respectively, contaminated by the mask image and the noise with the standard deviation 𝜎=12. The parameters of these three models are chosen for getting better restored images and the algorithms will be stopped after 500 iterations (see also Table 3). It is obvious that these three models can efficiently restore the original deteriorated images. However, there are much more inappropriate information restored by the ROF model than the LLT model and the hybrid model. These advantages is based on the fact that the fourth-order linear diffusion (the LLT model and the hybrid model) damps oscillations much faster than second-order diffusion (the ROF model). Especially, these unsuitable effects can be easily seen from the local zooming images in Figures 6 and 7. Furthermore, the restored image by the hybrid model looks more natural than the LLT model.

Table 3 

The related results in Example 4.12.

(a) Original image

(b) Mask image

(c) Painting image

(d) Noisy image

Figure 5 

The related images in Example 4.12.

Figure 6 

The related restored images corresponding to the painting image in Example 4.12.

Figure 7 

The related restored images corresponding to the noisy image in Example 4.12.

4.2.3. Image Denoising for the Multiplicative Noise

Based on the hybrid model (1.7), the multiplicative noise model (1.5) can be naturally extended to the following minimization: min_𝑧𝛼_Ω𝑧+𝑓𝑒^−𝑧d𝑥+_Ω1−𝑔∇𝑧d𝑥+_Ω𝑔∇²𝑧d𝑥. (4.19) It is obvious that the problem (4.19) can be approximated by𝑢^𝑛∶=argmin_𝑢1𝑧2𝜇^𝑛−𝑢²_𝐿2+𝛼_Ω𝑢+𝑓𝑒^−𝑢𝑧d𝑥,_𝑛+1∶=argmin_𝑧12𝜇𝑧−𝑢^𝑛2_𝐿2+_Ω1−𝑔∇𝑧d𝑥+_Ω𝑔∇²𝑧d𝑥. (4.20a)(4.20b)For the first subproblem (4.20a), its solution 𝑢^𝑛 can be approximately obtained by using the Newton method to solve the following nonlinear equation: 𝜇𝑢−𝑧^𝑛+𝛼1−𝑓𝑒^−𝑢=0. (4.21)

Example 4.13. In this example, we restore the multiplicative noisy image. The noisy Lena image shown in Figure 7 is contaminated by the Gamma noise (𝐿=33) with mean one, in which probability density function 𝑝(𝑢) is given by 𝑝(𝑢)=𝐿^𝐿𝑢^𝐿−1Γ(𝐿)𝑒^−𝐿𝑢,if𝑢>0,0,if𝑢≤0, (4.22) where 𝐿 is an integer and Γ(⋅) is a Gamma function. Based on the iteration formula (4.20a)-(4.20b), we should notice that there are two interior iterations. For employing the Newton method to solve the problem (4.21), we set the stepsize 𝑡=1 and the Newton method will be stopped when ‖𝑢^𝑛+1−𝑢^𝑛‖_ℓ²/‖𝑢^𝑛‖_ℓ²≤1.0×10⁻⁵. For solving the second subproblem (4.20b), we set the fixed iterations with 40. In fact for using the PPM to solve the ROF model, the LLT model and the hybrid model, the solutions of these models will tend to the steady. The outer iteration will be stopped after 200 iterations. The related restored images are shown in Figure 8, and it is easy to find that the hybrid model gives a better restored image than the two other models.

Figure 8 

The related images and the local zooming images in Example 4.13. SNRs: (a) 9.2017; (b) 15.6059; (c) 15.1754; (d) 16.0546.

5. Concluding Remarks

In this paper, based on the edge detector function, we proposed a hybrid model to overcome some drawbacks of the ROF model and the LLT model. Following the augmented Lagrangian method, we can employ the ADMM of multipliers to solve this hybrid model. In this paper, we mainly proposed the PPM to solve this model due to the fact that the PPM unnecessarily solves a PDE compared with the ADMM so that it is more effective than the ADMM. The convergence of the proposed method was also given. However, the convergence rate of the proposed method is only 𝑂(1/𝑘), so our future work is to improve our method with convergence rate 𝑂(1/𝑘²) based on the same strategies in [51, 52].

Appendices

A. Proof of Theorem 3.2

Proof. Following the assumption, we can find that the saddle point (𝑢^∗,𝑣^∗,𝑧^∗,𝑤^∗₁,𝑤^∗₂) of (3.9) satisfies 0=𝜃𝑢^∗−𝑓+Λ^∗₁𝑤^∗₁+Λ^∗₂𝑤^∗₂,0=𝑝^∗₁−𝑤^∗₁,0=𝑝^∗₂−𝑤^∗₂,0=Λ₁𝑢^∗−𝑣^∗,0=Λ₂𝑢^∗−𝑧^∗ (A.1) with 𝑝^∗₁∈𝜕ℎ₁(𝑣^∗) and 𝑝^∗₂∈𝜕ℎ₂(𝑧^∗). If we set 𝑑^∗₁=𝑤^∗₁/𝜇₁ and 𝑑^∗₂=𝑤^∗₂/𝜇₂, then (A.1) can be written as 0=𝜃𝑢^∗−𝑓+𝜇₁Λ^∗₁Λ₁𝑢^∗−𝑣^∗+𝑑^∗₁+𝜇₂Λ^∗₂Λ₂𝑢^∗−𝑧^∗+𝑑^∗₂,0=𝑤^∗₁+𝜇₁𝑣^∗−𝑑^∗₁−Λ₁𝑢^∗,0=𝑤^∗₂+𝜇₂𝑧^∗−𝑑^∗₂−Λ₂𝑢^∗,𝑑^∗₁=𝑑^∗₁+Λ₁𝑢^∗−𝑣^∗,𝑑^∗₂=𝑑^∗₂+Λ₂𝑢^∗−𝑧^∗. (A.2a)(A.2b)(A.2c)(A.2d)(A.2e)Then the above equations (A.2a)–(A.2e) can be rewritten as a compact form 0=𝜃𝑢^∗−𝑓+Λ^∗₁𝜕ℎ₁Λ₁𝑢^∗+Λ^∗₂𝜕ℎ₂Λ₂𝑢^∗, (A.3) which implies that 𝑢^∗ is the solution of the minimization problem (1.8).
Set 𝑢^𝑛_𝑒=𝑢^𝑛−𝑢^∗,𝑣^𝑛_𝑒=𝑣^𝑛−𝑣^∗,𝑧^𝑛_𝑒=𝑧^𝑛−𝑧^∗,𝑑^𝑛_1𝑒=𝑑^𝑛₁−𝑑^∗₁,𝑑^𝑛_2𝑒=𝑑^𝑛₂−𝑑^∗₂ (A.4) denote the related errors. Then subtracting (3.6a)–(3.6e) by (A.2a)–(A.2e) and using the similar strategy as in [24], we successively obtain 0=𝜃𝑢^𝑛_𝑒²_𝐿2+𝜇₁Λ₁𝑢^𝑛_𝑒²_𝐿2+𝜇₂Λ₂𝑢^𝑛_𝑒²_𝐿2−𝜇₁𝑣_𝑒^𝑛−1−𝑑^𝑛−1_1𝑒,Λ₁𝑢^𝑛_𝑒−𝜇₂𝑧_𝑒^𝑛−1−𝑑^𝑛−1_2𝑒,Λ₂𝑢^𝑛_𝑒,0=𝜇₁𝑣^𝑛_𝑒²_𝐿2+𝜇₁𝑣^𝑛_𝑒,𝑑^𝑛_1𝑒−𝑑^𝑛−1_1𝑒−Λ₁𝑢^𝑛_𝑒,0=𝜇₂𝑧^𝑛_𝑒²_𝐿2+𝜇₂𝑧^𝑛_𝑒,𝑑^𝑛_2𝑒−𝑑^𝑛−1_2𝑒−Λ₂𝑢^𝑛_𝑒,0=𝜇₁𝑑^𝑛_1𝑒²_𝐿2−𝜇₁𝑑^𝑛−1_1𝑒²_𝐿2+2𝜇₁𝑑^𝑛−1_1𝑒,𝑣^𝑛_𝑒−Λ₁𝑢^𝑛_𝑒−𝜇₁𝑣^𝑛_𝑒−Λ₁𝑢^𝑛_𝑒²_𝐿2,0=𝜇₂𝑑^𝑛_2𝑒²_𝐿2−𝜇₂𝑑^𝑛−1_2𝑒²_𝐿2+2𝜇₂𝑑^𝑛−1_2𝑒,𝑧^𝑛_𝑒−Λ₂𝑢^𝑛_𝑒−𝜇₂𝑧^𝑛_𝑒−Λ₂𝑢^𝑛_𝑒²_𝐿2, (A.5) where 𝑑^𝑛_1𝑒∈(1/𝜇₁)𝜕ℎ₁(𝑣^𝑛_1𝑒) and 𝑑^𝑛_2𝑒∈(1/𝜇₂)𝜕ℎ₂(𝑧^𝑛_1𝑒). Summing the above five equations, we get 𝜇₁2𝑑^𝑛−1_1𝑒²_𝐿2−𝑑^𝑛_1𝑒²_𝐿2+𝜇₂2𝑑^𝑛−1_2𝑒²_𝐿2−𝑑^𝑛_2𝑒²_𝐿2=𝜃𝑢^𝑛_𝑒²_𝐿2+𝜇₁𝑑^𝑛_1𝑒,𝑣^𝑛_𝑒+12Λ₁𝑢^𝑛_𝑒−𝑣_𝑒^𝑛−12_𝐿2+12𝑣^𝑛_𝑒²_𝐿2−12𝑣_𝑒^𝑛−12_𝐿2+𝜇₂𝑑^𝑛_2𝑒,𝑧^𝑛_𝑒+12Λ₂𝑢^𝑛_𝑒−𝑧_𝑒^𝑛−12_𝐿2+12𝑧^𝑛_𝑒²_𝐿2−12𝑧_𝑒^𝑛−12_𝐿2. (A.6) So it is easy to deduce that 𝜇₁2𝑑¹_1𝑒²_𝐿2−𝑑^𝑁_1𝑒²_𝐿2+𝜇₂2𝑑¹_2𝑒²_𝐿2−𝑑^𝑁_2𝑒²_𝐿2=^𝑁_𝑛=1𝜃𝑢^𝑛_𝑒²_𝐿2+𝜇₁𝑑^𝑛_1𝑒,𝑣^𝑛_𝑒+𝜇₁2Λ₁𝑢^𝑛_𝑒−𝑣_𝑒^𝑛−12_𝐿2+𝜇₂𝑑^𝑛_2𝑒,𝑧^𝑛_𝑒+𝜇₂2Λ₂𝑢^𝑛_𝑒−𝑧_𝑒^𝑛−12_𝐿2+𝜇₁2𝑣^𝑁_𝑒²_𝐿2−𝜇₁2𝑣¹_𝑒²_𝐿2+𝜇₂2𝑧^𝑁_𝑒²_𝐿2−𝜇₂2𝑧¹_𝑒²_𝐿2. (A.7) Following the fact that ℎ(𝑥) is strictly convex, 𝑑^𝑛_1𝑒∈(1/𝜇)𝜕ℎ₁(𝑣^𝑛_𝑒) and 𝑑^𝑛_2𝑒∈(1/𝜇)𝜕ℎ₂(𝑧^𝑛_𝑒), we can get 𝑑^𝑛_1𝑒,𝑣^𝑛_𝑒≥0,𝑑^𝑛_2𝑒,𝑧^𝑛_𝑒≥0. (A.8) Combined with (A.7), it follows that ^𝑁_𝑛=1𝑢^𝑛_𝑒²_𝐿2<∞,^𝑁_𝑛=1Λ₁𝑢^𝑛_𝑒−𝑣_𝑒^𝑛−12_𝐿2<∞,^𝑁_𝑛=1Λ₂𝑢^𝑛_𝑒−𝑧_𝑒^𝑛−12_𝐿2<∞. (A.9) Then we deduce that lim_𝑛→∞𝑢^𝑛=𝑢^∗,lim_𝑛→∞Λ₁𝑢^𝑛_𝑒−𝑣_𝑒^𝑛−12_𝐿2=0,lim_𝑛→∞Λ₂𝑢^𝑛_𝑒−𝑧_𝑒^𝑛−12_𝐿2=0. (A.10) So we have lim_𝑛→∞𝑣^𝑛−𝑣^∗2_𝐿2=0,lim_𝑛→∞𝑧^𝑛−𝑧^∗2_𝐿2=0. (A.11) Hence the proof is complete.

B. Proof of Lemma 3.4

Proof. Let 𝑥_𝑖(𝑖=1,2) be solutions of 𝑥_𝑖∈𝑦_𝑖+𝑐𝐴𝑦_𝑖, (B.1) then we can obtain 𝑥₁−𝑥₂²_𝐿2=𝑦₁−𝑦₂+𝑐𝐴𝑦₁−𝐴𝑦₂²_𝐿2=𝑦₁−𝑦₂²_𝐿2+𝑐²𝐴𝑦₁−𝐴𝑦₂²_𝐿2+2𝑐𝑦₁−𝑦₂,𝐴𝑦₁−𝐴𝑦₂≥𝑦₁−𝑦₂²_𝐿2+𝑐²𝐴𝑦₁−𝐴𝑦₂²_𝐿2=𝐽_𝑐𝐴𝑥₁−𝐽_𝑐𝐴𝑥₂²_𝐿2+𝑐²𝐻_𝑐𝐴𝑥₁−𝐻_𝑐𝐴𝑥₂²_𝐿2, (B.2) which implies that the assertion holds.

C. Proof of Theorem 3.5

Proof. Setting 𝐴₁(𝑥)∶=𝜕ℎ₁(𝑥) and 𝐴₂(𝑥)∶=𝜕ℎ₂(𝑥), then operators 𝐴₁ and 𝐴₂ are maximal monotone. Following from (3.11a)–(3.11e) and the definition of the Yosida approximation, (3.13) can be equivalently written as 𝜇₁𝑑^𝑛₁=𝐻_𝐴1/𝜇1Λ₁𝑢^𝑛−1+𝑑₁^𝑛−1,𝜇₂𝑑^𝑛₂=𝐻_𝐴2/𝜇2Λ₂𝑢^𝑛−1+𝑑₂^𝑛−1,𝑢^𝑛𝜇=𝑓−¹𝜃Λ^∗₁𝑑^𝑛₁+𝜇₂𝜃Λ^∗₂𝑑^𝑛₂. (C.1a)(C.1b)(C.1c)Then, by (C.1a)–(C.1c) and Lemma 3.4, we can deduce that 𝐽_𝐴1/𝜇1Λ₁𝑢+𝑑₁−𝐽_𝐴1/𝜇1Λ₁𝑢^𝑛−1+𝑑₁^𝑛−12_𝐿2+𝜇₁𝑑₁−𝑑^𝑛₁²_𝐿2≤Λ₁𝑢−𝑢^𝑛−1+𝑑₁−𝑑₁^𝑛−12_𝐿2,𝐽_𝐴2/𝜇2Λ₂𝑢+𝑑₂−𝐽_𝐴2/𝜇2Λ₂𝑢^𝑛−1+𝑑₂^𝑛−12_𝐿2+𝜇₂𝑑₂−𝑑^𝑛₂²_𝐿2≤Λ₂𝑢−𝑢^𝑛−1+𝑑₂−𝑑₂^𝑛−12_𝐿2. (C.2) Combine (C.2) and (C.3) with Λ₁𝑢−𝑢^𝑛−1,𝜇₁𝑑₁−𝑑₁^𝑛−1+Λ₂𝑢−𝑢^𝑛−1,𝜇₂𝑑₂−𝑑₂^𝑛−11=−𝜃𝑢−𝑢^𝑛−12_𝐿21≤−2𝜃Λ₁²Λ₁𝑢−𝑢^𝑛−12_𝐿2−12𝜃Λ₂²Λ₂𝑢−𝑢^𝑛−12_𝐿2, (C.3) which implies that 𝜇₁𝐽_𝐴1/𝜇1Λ₁𝑢+𝑑₁−𝐽_𝐴1/𝜇1Λ₁𝑢^𝑛−1+𝑑₁^𝑛−12_𝐿2+𝜇₁𝜇₁𝑑₁−𝑑^𝑛₁²_𝐿2+𝜇₂𝐽_𝐴2/𝜇2Λ₂𝑢+𝑑₂−𝐽_𝐴2/𝜇2Λ₂𝑢^𝑛−1+𝑑₂^𝑛−12_𝐿2+𝜇₂𝜇₂𝑑₂−𝑑^𝑛₂²_𝐿2≤𝜇₁−1𝜃Λ₁²Λ₁𝑢−𝑢^𝑛−12_𝐿2+𝜇₁𝑑₁−𝑑₁^𝑛−12_𝐿2+𝜇₂−1𝜃Λ₂²Λ₂𝑢−𝑢^𝑛−12_𝐿2+𝜇₂𝑑₂−𝑑₂^𝑛−12_𝐿2. (C.4) Since 𝜇₁≤1/𝜃‖Λ₁‖² and 𝜇₂≤1/𝜃‖Λ₂‖², by (C.3) we eventually get 𝑑₁−𝑑^𝑛_1𝐿2≤𝑑₁−𝑑₁^𝑛−1_𝐿2,𝑑₂−𝑑^𝑛_2𝐿2≤𝑑₂−𝑑₂^𝑛−1_𝐿2 (C.5) as long as 𝑢^𝑛≠𝑢 and lim_𝑛→∞𝑢^𝑛−1−𝑢_𝐿2=0. (C.6) Furthermore, by (C.4) we also get lim_𝑛→∞𝐽_𝐴1/𝜇1Λ₁𝑢+𝑑₁−𝐽_𝐴1/𝜇1Λ₁𝑢^𝑛−1+𝑑₁^𝑛−1_𝐿2=0,lim_𝑛→∞𝐽_𝐴2/𝜇2Λ₁𝑢+𝑑₂−𝐽_𝐴2/𝜇2Λ₂𝑢^𝑛−1+𝑑₂^𝑛−1_𝐿2=0. (C.7) Based on (C.1a) and (C.1b), using the definitions of resolvent and the Yosida approximation, it is easy to find that lim_𝑛→∞𝑑^𝑛₁−𝑑₁^𝑛−1_𝐿2=0,lim_𝑛→∞𝑑^𝑛₂−𝑑₂^𝑛−1_𝐿2=0. (C.8) Then, using Corollary 4 in page 199 of [53], we can deduce that lim_𝑛→∞𝑑^𝑛₁=𝑑₁ and lim_𝑛→∞𝑑^𝑛₂=𝑑₂.
On the other hand, if (𝑑₁,𝑑₂,𝑢) is the limit of the sequence {(𝑑^𝑛₁,𝑑^𝑛₂,𝑢^𝑛)} generated by Algorithm 3.3, it is easy to find that (3.13) can be rewritten as 𝑑₁=prox_𝜇1ℎ^∗₁Λ₁𝑢+𝑑₁,𝑑₂=prox_𝜇2ℎ^∗₂Λ₂𝑢+𝑑₂,𝜇𝑢=𝑓−¹𝜃Λ^∗₁𝑑₁+𝜇₂𝜃Λ^∗₂𝑑₂ (C.9) as 𝑛→+∞. Using the relationship (2.9) in Remark 2.8, (C.1a)–(C.1c) can be then rewritten as Λ₁𝑢=prox_1/𝜇1ℎ₁Λ₁𝑢+𝑑₁,Λ₂𝑢=prox_1/𝜇2ℎ₂Λ₂𝑢+𝑑₂,𝜇𝑢=𝑓−¹𝜃Λ^∗₁𝑑₁+𝜇₂𝜃Λ^∗₂𝑑₂. (C.10a)(C.10b)(C.10c) Following from (C.10a) and (C.10b), we can find that 𝑑₁=1𝜇₁𝜕ℎ₁Λ₁(𝑢),𝑑₂=1𝜇₂𝜕ℎ₂Λ₂(𝑢). (C.11) Submitting (C.10a)–(C.10c) into (C.11), we have 𝜃𝑢−𝑓+Λ^∗₁𝜕ℎ₁Λ₁(𝑢)+Λ^∗₂𝜕ℎ₂Λ₂(𝑢)=0, (C.12) which implies that 𝑢 is the solution of (1.8).

D. Proof of Theorem 4.9

Proof. From the assumption, the functional ℳ(𝑢,𝑧) exists in at least a minimized point denoted by (𝑢,𝑧), which implies that (𝑢,𝑧) satisfies 𝑧=ℱ(𝑢),𝑢=𝒢(𝑧). (D.1) That is to say that 𝑢 is a fixed point of the operator 𝒢∘ℱ.
Now we show that the sequence {𝑢^𝑛} generated by (4.11a)-(4.11b) converges to this fixed point 𝑢. In fact, based on the nonexpansive operators ℱ and 𝒢, for the sequence {𝑢^𝑛} we can get 𝑢^𝑛+1−𝑢_𝐿2=ℱ∘𝒢𝑢^𝑛−ℱ∘𝒢(𝑢)_𝐿2≤𝑢^𝑛−𝑢_𝐿2. (D.2) This implies that the positive sequence {‖𝑢^𝑛−𝑢‖_𝐿²} is monotonically decreasing. So we can find a subsequence {𝑢^𝑛𝑘} from {𝑢^𝑛} converging to a limit point 𝑢^∗ as 𝑛→∞. Form the continuity of ℱ and 𝒢, we then have 𝑢^∗−𝑢_𝐿2=ℱ∘𝒢𝑢^∗−𝑢_𝐿2. (D.3) Furthermore, we can get 𝑢=𝑢^∗ based on the uniqueness of the fixed point. Thus, (𝑢^∗,𝑧^∗) is the solution of (4.11a)-(4.11b), where 𝑧^∗ corresponds to 𝑢^∗. Setting 𝐽(𝑢)=∫_Ω|(1−𝑔)∇𝑢|d𝑥+∫_Ω|𝑔∇²𝑢|d𝑥, we can get 0=𝛽𝑧^∗−𝑢^∗+ℎ^′𝑧^∗,0∈𝛽𝑢^∗−𝑧^∗+𝜕𝐽𝑢^∗, (D.4) which implies that ℎ^′𝑧^∗+𝜕𝐽𝑢^∗=0. (D.5) Then 𝑢^∗ is the solution of (4.13) when 𝑢^∗=𝑧^∗.

Acknowledgments

The second author would like to thank the Mathematical Imaging and Vision Group in Division of Mathematical Sciences at Nanyang Technological University in Singapore for the hospitality during the visit. This research is supported by Singapore MOE Grant T207B2202, Singapore NRF2007IDM-IDM002-010, the NNSF of China (no. 60835004, 60872129), and the University Research Fund of Henan University (no. 2011YBZR003).