New Regularization Models for Image Denoising with a Spatially Dependent Regularization Parameter

Abstract

We consider simultaneously estimating the restored image and the spatially dependent regularization parameter which mutually benefit from each other. Based on this idea, we refresh two well-known image denoising models: the LLT model proposed by Lysaker et al. (2003) and the hybrid model proposed by Li et al. (2007). The resulting models have the advantage of better preserving image regions containing textures and fine details while still sufficiently smoothing homogeneous features. To efficiently solve the proposed models, we consider an alternating minimization scheme to resolve the original nonconvex problem into two strictly convex ones. Preliminary convergence properties are also presented. Numerical experiments are reported to demonstrate the effectiveness of the proposed models and the efficiency of our numerical scheme.

1. Introduction

Image denoising is a fundamental problem in image processing and computer vision. In many real-world applications, it forms a significant preliminary step for subsequent image processing operations, such as object recognition, medical image analysis, surveillance, and many more. During acquisition and transmission, images are often corrupted by Gaussian noise. In this problem, the degradation process is modeled as 𝑓=𝑢_exact+𝜂, (1) where 𝑓, 𝑢_exact, and 𝜂 represent the observed image, the original image, and the additive white Gaussian noise, respectively.

The objective of image denoising is to compute a good estimate of 𝑢_exact from 𝑓. To obtain a reasonable approximated solution from (1), the regularization method, generally used as a numerical technique for stabilizing inverse problems, has been increasingly applied to image denoising over the past decades. A large class of regularization based denoising methods unifies under the following framework: min_𝑢𝐸(𝑢)=𝐹(𝑢)+𝛼𝑅(𝑢), (2) where 𝑢 denotes the restored image to be estimated. In the right-hand side of (2), 𝐹(𝑢) represents the fidelity term, which measures the closeness of the estimate to the data. The functional 𝑅(𝑢) is called the regularization term pushing 𝑢 to exhibit some a priori expected features. Parameter 𝛼 in (2) is known as the regularization parameter, which balances the tradeoff between the two terms.

One key topic in regularization methods is the choice of the regularizer. Among many regularization based denoising models, the total variation regularization, proposed by Rudin, Osher, and Fatemi (ROF) [1], has won tremendous success due to its edge-preserving property. In the ROF model, the restoration result 𝑢 is generated by solving the following minimization problem: min_𝑢𝐸(𝑢)=_Ω𝑢−𝑓²𝑑𝑥+𝛼𝑅_ROF(𝑢), (3) where the TV-term 𝑅_ROF(𝑢)=∫_Ω|𝐷𝑢| denotes the total variation of 𝑢 (see [2] for more details) and Ω denotes the image domain. A remarkable aspect of the ROF model is that the TV-term does not penalize the discontinuities in 𝑢; see, for example, [3]. This property allows us to restore the edges of the original image. However, the main disadvantage of the ROF model is the so-called staircase effect (smooth regions are transformed into piecewise constant regions), a phenomenon long observed in the literature [4–6]. As a consequence, the restoration result is unsatisfactory to the eye due to the loss of texture details and the generation of artificial edges that do not exist in the true image. To overcome this effect, models based on high-order PDEs have been proposed in the literature; see, for example, [7–11]. For instance, Lysaker, Lundervold, and Tai [8] proposed a fourth-order PDE model (termed as the LLT model) with the form min_𝑢𝐸(𝑢)=_Ω𝑢−𝑓²𝑑𝑥+𝛼𝑅_LLT(𝑢), (4) where 𝑅_LLT(𝑢)=∫_Ω|∇²𝑢|𝑑𝑥 is a fourth-order filter. From a theoretical point of view [12], it has been shown that fourth-order PDEs are superior to second-order PDEs in some aspects, including avoiding the staircase effect. However, this type of filters usually blurs the edges of the original image and suffers from the speckle effect in homogeneous regions. Since the ROF model and the LLT model are of both merits and drawbacks, it may be desirable to promote solutions that simultaneously exhibit properties that are enforced by both regularizers. This is the basic idea of [13] proposed by Lysaker and Tai, in which they studied an iterative algorithm based on combining the results generated by (3) and (4). But their method needs to solve two separate PDEs and their combination is not quite intuitive. Li et al. [14] proposed a hybrid model in the following form: min_𝑢𝐸(𝑢)=_Ω𝑢−𝑓²𝑑𝑥+𝛼_Ω1−𝑔|𝐷𝑢|+_Ω𝑔∇²𝑢𝑑𝑥, (5) where the function 𝑔 is chosen as an edge detector 1/(1+𝛾+𝑘|∇𝐺_𝜎∗𝑓|²) to control the relative weights of the two regularizers (see [14] for details). By their selection of 𝑔, the second-order filter plays the dominant role where |∇𝐺_𝜎∗𝑓| is large (regions including sharp features), whereas the fourth-order filter plays the dominant role where |∇𝐺_𝜎∗𝑓| is small (homogeneous regions). Therefore, (5) reaps the benefits of both regularizers. More related works on the combination of the ROF model and the LLT model can be found in [15–17].

Another crucial issue in the regularization process is the suitable selection of the regularization parameter. In regularization models, the regularization parameter controls the relative weights of the data fidelity and regularization terms. However, due to the inhomogeneous distribution of cartoon, textures, and small details in an image (in terms of variance), a global constant parameter may not be suitable for these features of different scales. According to the cartoon pyramid model studied in [18], as the regularization parameter goes from small to large values, the corresponding solutions generated by the ROF model range from undersmoothed (textures are preserved, while noise remains almost unchanged in homogeneous regions) to oversmoothed (noise is reduced well, but significant details are lost). This suggests that a spatially dependent regularization parameter including, instead of a single constraint, a group of constraints adapted to different regions of the image is desired to obtain higher quality results. To this end, denoising methods with a spatially dependent regularization parameter have been extensively studied; see, for example, [18–22]. In these works, the automated selection strategies of the regularization parameter are based on local variance measures [18, 20–22] and the local statistical characteristics of the noise [21, 22].

Note that, in models (3), (4), and (5), the selection of the regularization parameter never accommodates the information of the current restored image. This triggers us to seek a new approach for regularization parameter estimation. Our main motivation is to simultaneously estimate both the restored image and the regularization parameter which mutually benefit from each other during the denoising procedure. We propose a general model with the following form: min_𝑢,𝑔𝐸𝑢,𝑔=𝐹(𝑢)+𝐹𝑔+𝑅𝑢,𝑔, (6) where 𝑔 denotes a spatially dependent regularization parameter, 𝐹(𝑢) and 𝐹(𝑔) are, respectively, fidelity terms for the two variables, and the binary function 𝑅(𝑢,𝑔) represents a regularizer. The advantages of (6) are presented as follows. First, instead of a global constant regularization parameter, a spatially dependent regularization parameter is more flexible to image features of different scales. Second, since the estimation of the two variables is derived simultaneously, the regularization parameter is changing more reasonablly by exploiting the more accuracy restored image instead of the observed images corrupted by noise. Here we remark that (6) is a general-purpose model which can incorporate various classical methods. Specially we focus on the fourth-order filter in (4) and the hybrid regularizer in (5). The refresh models are named as Models 1 and 2, respectively. Model 1 can suppress the speckle effect caused by the LLT model while less overregularizing textures. Model 2 has the advantage of better restoring textures and homogeneous regions while preserving edges. To overcome the nonconvexity of our models, we utilize an alternating minimization scheme to resolve the original nonconvex problem into two strictly convex ones. Thus, our models can be asymptotically solved.

The outline of the rest of the paper is as follows. In the next section, we give notations and discretizations. In Section 3, we introduce our models with some discussions. The following Section 4 presents the numerical scheme for solving the proposed models. In Section 5, numerical experiments are given to demonstrate the performance of the proposed models. Finally, we conclude the paper in Section 6.

2. Notations and Discretizations

From now on, we will restrict our attention to the discrete setting. We first introduce some notations. Without loss of generality, we assume that all the images in this paper are grayscale and have a square domain. Then we represent an image 𝑢 as an 𝑛×𝑛 matrix, where 𝑢_𝑖,𝑗 represents the intensity value of 𝑢 at pixel (𝑖,𝑗), for 𝑖,𝑗=1,…,𝑛. For the sake of simplicity, we assume that the image is periodically extended, and then the FFT can be adopted in our algorithm. It should be pointed out that adaption to other boundary conditions is not difficult in principle. In the rest of this paper, we let ‖⋅‖₂, ‖⋅‖_F, and ∘ denote the 2-norm, the Frobenius norm, and the Hadamard product, respectively. Let 𝑋 denote the Euclidean space ℝ^𝑛×𝑛. The usual inner product and Euclidean norm of 𝑋 are denoted as ⟨⋅,⋅⟩_𝑋 and ‖⋅‖_𝑋, respectively. Denote by 𝑌 the space 𝑋×𝑋 equipped with the inner product ⟨⋅,⋅⟩_𝑌 leading canonically to the norm ‖⋅‖_𝑌, that is; for 𝑝=(𝑝¹,𝑝²)∈𝑌 and 𝑞=(𝑞¹,𝑞²)∈𝑌, 𝑝,𝑞_𝑌=𝑝¹,𝑞¹_𝑋+𝑝²,𝑞²_𝑋,𝑝_𝑌=𝑝,𝑝_𝑌. (7) Moreover, for 𝑦=(𝑦¹,𝑦²)∈𝑌, |𝑦| denotes the 𝑛×𝑛 matrix whose element |𝑦|_𝑖,𝑗 is equal to ‖𝑦_𝑖,𝑗‖₂ with 𝑦_𝑖,𝑗=(𝑦¹_𝑖,𝑗,𝑦²_𝑖,𝑗). We denote the space 𝑌×𝑌 as 𝑍. The definitions of the inner product ⟨⋅,⋅⟩_𝑍, the norm ‖⋅‖_𝑍, and |𝑧| are analogous to those of 𝑌.

Now we introduce a discretized version of some necessary operators. For 𝑢∈𝑋, we introduce forward and backward difference operators as follows: 𝐷^∘+_𝑥𝑢_𝑖,𝑗=𝑢_𝑖,𝑗+1−𝑢_𝑖,𝑗𝑢,1≤𝑗≤𝑛−1,_𝑖,1−𝑢_𝑖,𝑛,𝑗=𝑛,_∘𝐷⁺_𝑦𝑢_𝑖,𝑗=𝑢_𝑖+1,𝑗−𝑢_𝑖,𝑗𝑢,1≤𝑖≤𝑛−1,_1,𝑗−𝑢_𝑛,𝑗,𝑖=𝑛,_∘𝐷⁻_𝑥𝑢_𝑖,𝑗=𝑢_𝑖,1−𝑢_𝑖,𝑛𝑢,𝑗=1,_𝑖,𝑗−𝑢_{𝑖,𝑗−1},2≤𝑗≤𝑛,_∘𝐷⁻_𝑦𝑢_𝑖,𝑗=𝑢_1,𝑗−𝑢_𝑛,𝑗𝑢,𝑖=1,_𝑖,𝑗−𝑢_{𝑖−1,𝑗},2≤𝑖≤𝑛. (8) Second-order difference operators can be expressed by using a recursive application of first-order difference operators; that is, the operator 𝐷^∘−+_𝑥𝑥 is defined by 𝐷^∘−+_𝑥𝑥𝑢_𝑖,𝑗∶=𝐷^∘𝐷⁻_𝑥^∘+_𝑥𝑢_𝑖,𝑗. (9) Other second order difference operators used in this paper such as 𝐷^∘+−_𝑥𝑥, 𝐷^∘++_𝑥𝑦, 𝐷^∘++_𝑦𝑥, 𝐷^∘−−_𝑥𝑦, 𝐷^∘−−_𝑦𝑥, 𝐷^∘−+_𝑦𝑦, and 𝐷^∘+−_𝑦𝑦 can be similarly defined. Based on the definitions above, we can introduce the discrete gradient operator. For 𝑢∈𝑋, ∇𝑢 is a vector of 𝑌, which is given by ∇𝑢=𝐷^∘+_𝑥𝑢,𝐷^∘+_𝑦𝑢. (10) The discrete total variation of 𝑢 is then given by 𝑅_ROF(𝑢)=_{1≤𝑖,𝑗≤𝑛}|∇𝑢|_𝑖,𝑗. (11) To discretize the fourth-order filter in the LLT model (4), we have to introduce the discrete Hessian operator. Here we adopt the definition in [23]. The discrete Hessian operator is a mapping 𝐻∶𝑋→𝑍, and for 𝑢∈𝑋, 𝐻𝑢 is defined by 𝐻𝑢=^∘𝐷⁻⁺_𝑥𝑥𝑢𝐷^∘++_𝑥𝑦𝑢_∘𝐷⁺⁺_𝑦𝑥𝑢𝐷^∘−+_𝑦𝑦𝑢. (12) Then, 𝑅_LLT(𝑢) can be discretized as 𝑅_LLT(𝑢)=_{1≤𝑖,𝑗≤𝑛}|𝐻𝑢|_𝑖,𝑗. (13) We also introduce two important operators: div∶𝑌→𝑋 and 𝐻^∗∶𝑍→𝑋, that is, the adjoint operator of −∇ and 𝐻, respectively. By analogy with the continuous setting, for 𝑢∈𝑋, 𝑦∈𝑌, and 𝑧∈𝑍, we want them to satisfy div𝑦,𝑢_𝑋=𝑦,−∇𝑢_𝑌,𝐻^∗𝑧,𝑢_𝑋=⟨𝑧,𝐻𝑢⟩_𝑍. (14) Then they are formulated as follows: div𝑦=𝐷^∘−_𝑥𝑦¹+𝐷^∘−_𝑦𝑦²,𝐻^∗𝑧=𝐷^∘+−_𝑥𝑥𝑧¹¹+𝐷^∘−−_𝑦𝑥𝑧¹²+𝐷^∘−−_𝑥𝑦𝑧²¹+𝐷^∘+−_𝑦𝑦𝑧²². (15) Finally, the composite operators Δ=div⋅∇ and 𝐻^∗𝐻 are also used.

3. The Proposed Models and Discussion

3.1. The Proposed Models

Arising from model (6) and exploiting the benefits of the LLT model (4) and the hybrid model (5), we consider and study the following models.

Model 1. One has min_{𝑢,𝑔∈𝑋}𝐸₁=1𝑢,𝑔2𝑢−𝑓²_𝑋𝑀+𝛼𝑔−𝛼⋅𝟏²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝑔²_𝑖,𝑗(𝐾(|𝐻𝑢|))_𝑖,𝑗. (16)

Model 2. One has min_{𝑢,𝑔1,𝑔2∈𝑋}𝐸₂𝑢,𝑔₁,𝑔₂=12𝑢−𝑓²_𝑋+𝛼𝑔₁−𝑀𝛼⋅𝟏²_𝑋+𝑔₂−𝑀𝛼⋅𝟏²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝑔₁²_𝑖,𝑗(𝐾(|∇𝑢|))_𝑖,𝑗+𝑔₂²_𝑖,𝑗(𝐾(|𝐻𝑢|))_𝑖,𝑗. (17) In the above two models, 𝛼 and 𝑀 are positive parameters, 𝟏 represents the matrix whose elements are equal to 1, and 𝐾∶𝑋→𝑋 denotes a discrete mean filter. More precisely, we choose an odd integer 𝑟 and define an 𝑟-by-𝑟 window centered at pixel (𝑖,𝑗) (with a periodic extension at the boundary); that is, Ω^𝑟_𝑖,𝑗=(𝑠,𝑡)∶min(|𝑠−𝑖|,𝑛−|𝑠−𝑖|)≤𝑟−12,min𝑡−𝑗,𝑛−𝑡−𝑗≤𝑟−12, (18) where 𝑠,𝑡=1,…,𝑛. Then for 𝑢∈𝑋, 𝐾(𝑢) is given by (𝐾(𝑢))_𝑖,𝑗=1𝑟²_{()𝑠,𝑡∈Ω𝑟𝑖,𝑗}𝑢_𝑠,𝑡. (19) It is immediately clear that if 𝑟=1, 𝐾 coincides with the identity operator. We start our investigations of the proposed models by the following lemma.

Lemma 1. For any 𝑢,𝑣∈𝑋, one has _{1≤𝑖,𝑗≤𝑛}𝑢_𝑖,𝑗(𝐾(𝑣))_𝑖,𝑗=_{1≤𝑖,𝑗≤𝑛}(𝐾(𝑢))_𝑖,𝑗𝑣_𝑖,𝑗. (20)

Proof. We define the following function 𝑘∶ℝ²→ℝ: 𝑘=1𝑥,𝑦𝑟²ifmin(|𝑥|,𝑛−|𝑥|)≤𝑟−12,min𝑦,𝑛−𝑦≤𝑟−12,0else. (21) From the definition of 𝑘, it is immediate to see that 𝑘(𝑥,𝑦)=𝑘(−𝑥,−𝑦). Thus, we have _{1≤𝑖,𝑗≤𝑛}𝑢_𝑖,𝑗(𝐾(𝑣))_𝑖,𝑗=_{1≤𝑖,𝑗≤𝑛}𝑢_{𝑖,𝑗1≤𝑠,𝑡≤𝑛}𝑘𝑖−𝑠,𝑗−𝑡𝑣_𝑠,𝑡=_{1≤𝑠,𝑡≤𝑛}𝑣_{s,𝑡1≤𝑖,𝑗≤𝑛}𝑘𝑠−𝑖,𝑡−𝑗𝑢_𝑖,𝑗=_{1≤𝑠,𝑡≤𝑛}(𝐾(𝑢))_𝑠,𝑡𝑣_𝑠,𝑡=_{1≤𝑖,𝑗≤𝑛}(𝐾(𝑢))_𝑖,𝑗𝑣_𝑖,𝑗. (22)

According to Lemma 1, we can rewrite the regularizers in our models in another equivalent form. In (16), for example, we rewrite _{1≤𝑖,𝑗≤𝑛}𝑔²_𝑖,𝑗(𝐾(|𝐻𝑢|))_𝑖,𝑗=_{1≤𝑖,𝑗≤𝑛}𝐾𝑔∘𝑔_𝑖,𝑗|𝐻𝑢|_𝑖,𝑗. (23)

Next we study the existence result for the solution of our models.

Theorem 2. Assume that 𝛼>0, 𝑀>0, and 𝑓∈𝑋. Then both Models 1 and 2 have a global minimizer.

Proof. First, under our assumption, both models are bounded from below, which implies that the minimization problems are the correct setting.
It is immediate to see that the functions 𝐸₁∶𝑋×𝑋→ℝ and 𝐸₂∶𝑋×𝑋×𝑋→ℝ are proper, coercive, and continuous. Then, the existence of a global minimizer is deduced by Weierstrass’ theorem [24].

We end this subsection with the following proposition which describes the relationship between the variables in a solution of our models.

Proposition 3. (i) Assume that (𝑢^∗,𝑔^∗)∈𝑋×𝑋 is a (local) minimizer of Model 1. Then (𝑢^∗,𝑔^∗) satisfies 𝑢^∗=argmin_𝑢𝐸₁𝑢,𝑔^∗,𝑔^∗=argmin_𝑔𝐸₁𝑢^∗,𝑔. (24)
(ii) Assume that (𝑢^∗,𝑔^∗₁,𝑔^∗₂)∈𝑋×𝑋×𝑋 is a (local) minimizer of Model 2. Then (𝑢^∗,𝑔^∗₁,𝑔^∗₂) satisfies 𝑢^∗=argmin_𝑢𝐸₂𝑢,𝑔^∗₁,𝑔^∗₂,𝑔^∗₁,𝑔^∗₂=argmin_𝑔1,𝑔2𝐸₂𝑢^∗,𝑔₁,𝑔₂. (25)

Proof. (i) Under our assumption, there exists some 𝜖>0 such that 𝐸₁𝑢^∗,𝑔^∗≤𝐸₁𝑢,𝑔∀𝑢,𝑔∈𝑋with𝑢−𝑢^∗_𝑋≤𝜖,𝑔−𝑔^∗_𝑋≤𝜖. (26) Then, we have 𝐸₁𝑢^∗,𝑔^∗≤𝐸₁𝑢,𝑔^∗∀𝑢∈𝑋with𝑢−𝑢^∗_𝑋≤𝜖, (27) which implies that 𝑢^∗ is a local minimizer of the following minimization problem: min_𝑢∈𝑋𝐸₁𝑢,𝑔^∗=12𝑢−𝑓²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝐾𝑔^∗∘𝑔^∗_𝑖,𝑗|𝐻𝑢|_𝑖,𝑗. (28) According to the nonnegativity of (𝐾(𝑔^∗∘𝑔^∗))_𝑖,𝑗, we can deduce that the objective function of (28) is strictly convex. Then, 𝑢^∗ is the unique global minimizer of (28). Similarly, we can obtain that 𝑔^∗ is the unique global minimizer of the following strictly convex minimization problem: min_𝑔∈𝑋𝐸₁𝑢^∗𝑀,𝑔=𝛼𝑔−𝛼⋅𝟏²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝑔²_𝑖,𝑗𝐾𝐻𝑢^∗_𝑖,𝑗. (29) Hence we complete the proof of (i).
(ii) Using the same technique, (ii) can be similarly proved, and we omit the details.

3.2. Discussion on the Proposed Models

In this subsection we discuss the strengths of the proposed models by investigating the numerical behavior of the solution. Proposition 3 implies that there is an interaction between the restored image and the regularization parameter in a solution of our models, which can achieve joint optimality. In other words, both variables benefit from each other. Assume that (𝑢^∗,𝑔^∗) and (𝑢^∗,𝑔^∗₁,𝑔^∗₂) are solutions of Models 1 and 2, respectively. According to Proposition 3, 𝑔^∗, 𝑔^∗₁, and 𝑔^∗₂ can be calculated directly from 𝑢^∗. Since the minimization subproblems with respect to 𝑔^∗ and 𝑔^∗₂ are actually the same, we just consider 𝑔^∗₁ and 𝑔^∗₂ for the sake of brevity. They are given by 𝑔^∗_1𝑖,𝑗=𝑀𝛼+(𝐾(|∇𝑢^∗|))_𝑖,𝑗,𝑔^∗_2𝑖,𝑗=𝑀𝛼+(𝐾(|𝐻𝑢^∗|))_𝑖,𝑗. (30) We are interested in the numerical behavior of 𝑔^∗₁ and 𝑔^∗₂ corresponding to regions of 𝑢^∗ with different scales (e.g., texture, flat, and ramp regions). Figure 1(a) shows a 256×256 synthetic image. The corresponding 𝑔^∗₁ and 𝑔^∗₂ are shown in Figures 1(b) and 1(c), respectively, where 𝑢^∗ is set to be the original image. Form Figure 1, numerical behavior of 𝑔^∗₁ and 𝑔^∗₂ is summarized as follows.(i)For texture regions of 𝑢^∗, 𝑔^∗₁ and 𝑔^∗₂ are both small (darker regions indicate smaller value).(ii)For flat regions of 𝑢^∗, 𝑔^∗₁ and 𝑔^∗₂ are both large.(iii)For ramp regions of 𝑢^∗, 𝑔^∗₁ is inversely proportional to the gradient of 𝑢^∗, whereas 𝑔^∗₂ is large.

(a)

(b)

(c)

Figure 1 

The numerical behavior of 𝑔∗1 and 𝑔∗2; here 𝛼=0.05, 𝑀=0.05, and 𝑟=7: (a) the original image, (b) 𝑔∗1, and (c) 𝑔∗2.

Below we pay attention to the restored image 𝑢^∗. As we have mentioned in the Introduction, values of the regularization parameter control the relative weights of the fidelity and regularization terms. More precisely, small values lead to little regularization, whereas large values result in overregularization. Based on the behavior analysis above, in Model 1, textures of 𝑢^∗ are not compromised due to the corresponding small values of 𝑔^∗, and the speckle effect caused by the LLT model is removed from flat and ramp regions of 𝑢^∗ due to the corresponding large values of 𝑔^∗. Therefore, Model 1 can produce higher quality restoration results than the LLT model. Model 2 inherits the advantages of Model 1 and exhibits some new ones. Edges of 𝑢^∗ are well preserved by total variation regularization. At the same time, the staircase effect is suppressed in ramp regions due to the corresponding small values of 𝑔^∗₁ (the fourth-order filter plays the dominate role). Therefore, Model 2 incorporates the strengths of both regularizers while avoiding their drawbacks. Finally, we remark that, compared with the hybrid model (5), Model 2 is more reasonable. Observe that (5) utilizes a fixed weighting function to combine the two regularizers. However, the weighting function is computed from the noisy image, which leads to inaccuracy. On the other hand, the regularization parameter in Model 2 depends on the restored image, which is less affected by noise.

4. Algorithms

In this section, we formulate the numerical scheme for solving our models. Our basic idea is to use the alternating minimization scheme which is described in Algorithms 1 and 2.

Input:  𝑓, 𝛼, 𝑀, and 𝑟. Output:  𝑢𝑘. Initialization:  𝑔0=0.   while  not converged do    Step  1. Given 𝑔𝑘−1, computing 𝑢𝑘 by solving   min𝑢∈𝑋𝐸1𝑢,𝑔𝑘−1=12𝑢−𝑓2𝑋+∑1≤𝑖,𝑗≤𝑛𝐾𝑔𝑘−1∘𝑔𝑘−1𝑖,𝑗|𝐻𝑢|𝑖,𝑗.    (a)    Step  2. Given 𝑢𝑘, computing 𝑔𝑘 by solving   min𝑔∈𝑋𝐸1𝑢𝑘𝑀,𝑔=𝛼𝑔−𝛼⋅𝟏2𝑋+∑1≤𝑖,𝑗≤𝑛𝑔2𝑖,𝑗𝐾𝐻𝑢𝑘𝑖,𝑗.      (b)   end while

Algorithm 1 

The alternating minimization method for solving Model 1.

Input:  𝑓, 𝛼, 𝑀, and 𝑟. Output:  𝑢𝑘. Initialization:  𝑔01=0 and 𝑔02=0.   while not converged do    Step  1. Given 𝑔1𝑘−1 and 𝑔2𝑘−1, computing 𝑢𝑘 by solving      min𝑢∈𝑋𝐸2𝑢,𝑔1𝑘−1,𝑔2𝑘−1=12𝑢−𝑓2𝑋+∑1≤𝑖,𝑗≤𝑛𝐾𝑔1𝑘−1∘𝑔1𝑘−1𝑖,𝑗|∇𝑢|𝑖,𝑗             +∑1≤𝑖,𝑗≤𝑛𝐾𝑔2𝑘−1∘𝑔2𝑘−1𝑖,𝑗|𝐻𝑢|𝑖,𝑗.       (a)    Step  2. Given 𝑢𝑘, computing 𝑔𝑘1 and 𝑔𝑘2 by solving     min𝑔1,𝑔2∈𝑋𝐸2𝑢𝑘,𝑔1,𝑔2=𝛼𝑔1−𝑀𝛼⋅𝟏2𝑋+𝑔2−𝑀𝛼⋅𝟏2𝑋             +∑1≤𝑖,𝑗≤𝑛𝑔12𝑖,𝑗𝐾∇𝑢𝑘𝑖,𝑗+𝑔22𝑖,𝑗𝐾𝐻𝑢𝑘𝑖,𝑗.  (b)   end while

Algorithm 2 

The alternating minimization method for solving Model 2.

Next we would like to show the convergence of our algorithms.

Theorem 4. (i) Let {(𝑢^𝑘,𝑔^𝑘)} be the sequence derived from Algorithm 1. Then {(𝑢^𝑘,𝑔^𝑘)} converges to (𝑢,𝑔)∈𝑋×𝑋 (up to a subsequence), and for any (𝑢,𝑔)∈𝑋×𝑋, one has 𝐸₁𝑢,𝑔≤𝐸₁𝑢,𝑔,𝐸₁𝑢,𝑔≤𝐸₁𝑢,𝑔. (31)
(ii) Let {(𝑢^𝑘,𝑔^𝑘₁,𝑔^𝑘₂)} be the sequence derived from Algorithm 2. Then {(𝑢^𝑘,𝑔^𝑘₁,𝑔^𝑘₂)} converges to (𝑢,𝑔₁,𝑔₂)∈𝑋×𝑋×𝑋 (up to a subsequence), and for any (𝑢,𝑔₁,𝑔₂)∈𝑋×𝑋×𝑋, one has 𝐸₂𝑢,𝑔₁,𝑔₂≤𝐸₂𝑢,𝑔₁,𝑔₂,𝐸₂𝑢,𝑔₁,𝑔₂≤𝐸₂𝑢,𝑔₁,𝑔₂. (32)

Proof. (i) First, we can easily deduce the following inequality from Algorithm 1: 𝐸₁𝑢^𝑘+1,𝑔^𝑘+1≤𝐸₁𝑢^𝑘+1,𝑔^𝑘≤𝐸₁𝑢^𝑘,𝑔^𝑘. (33) Then the sequence {𝐸₁(𝑢^𝑘,𝑔^𝑘)} is bounded, and there exists a constant 𝑁≥0 such that 𝐸₁𝑢^𝑘,𝑔^𝑘≤𝑁,∀𝑘∈ℕ. (34) The above inequality reads as 12𝑢^𝑘−𝑓²_𝑋+𝛼𝑔^𝑘−𝑀𝛼⋅𝟏²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝑔^𝑘2_𝑖,𝑗𝐾𝐻𝑢^𝑘_𝑖,𝑗≤𝑁, (35) which implies that {(𝑢^𝑘,𝑔^𝑘)} is uniformly bounded in 𝑋×𝑋. Then we can find a subsequence {(𝑢^𝑛𝑘,𝑔^𝑛𝑘)}⊂{(𝑢^𝑘,𝑔^𝑘)} and (𝑢,𝑔)∈𝑋×𝑋 such that they satisfy 𝑢^𝑛𝑘,𝑔−−−−→^{𝑛𝑘𝑋×𝑋}𝑢,𝑔. (36)
On the other hand, for any 𝑢∈𝑋, we have 𝐸₁𝑢^𝑛𝑘+1,𝑔^𝑛𝑘+1≤𝐸₁𝑢^𝑛𝑘+1,𝑔^𝑛𝑘+1≤𝐸₁𝑢^𝑛𝑘+1,𝑔^𝑛𝑘≤𝐸₁𝑢,𝑔^𝑛𝑘. (37) Recall that the function 𝐸₁∶𝑋×𝑋→ℝ is continuous; we obtain 𝐸₁𝑢,𝑔≤𝐸₁𝑢,𝑔, (38) by letting 𝑘 tend to infinity. Similarly, for any 𝑔∈𝑋, we have 𝐸₁𝑢,𝑔≤𝐸₁𝑢,𝑔. (39) Hence we complete the proof of (i).
(ii) Using the same technique, (ii) can be similarly proofed, and we omit the details.

It is obvious that (b) in Algorithm 1 and (b) in Algorithm 2 can be easily solved. To solve (a) in Algorithm 1 and (a) in Algorithm 2, we use the alternating direction method (ADM) [25], which is closely related to the augmented Lagrangian method [23], the Douglas-Rachford splitting algorithm [26], and the split Bregman method [27]. For the sake of brevity, we only present the ADM procedure for solving (a) in Algorithm 2. The ADM procedure for solving (a) in Algorithm 1 can be analogously derived. We rewrite (a) in Algorithm 2 by introducing auxiliary variables 𝑝 and 𝑞 as follows: min_{𝑢∈𝑋,𝑝∈𝑌,𝑞∈𝑍}12𝑢−𝑓²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝐾𝑔₁^𝑘−1∘𝑔₁^𝑘−1_𝑖,𝑗𝑝_𝑖,𝑗+_{1≤𝑖,𝑗≤𝑛}𝐾𝑔₂^𝑘−1∘𝑔₂^𝑘−1_𝑖,𝑗𝑞_𝑖,𝑗subjectto𝑝=∇𝑢,𝑞=𝐻𝑢. (40) Then the problem fits the framework of the alternating direction method. By using the Lagrangian multipliers 𝜆₁∈𝑌 and 𝜆₂∈𝑍, the augmented Lagrangian function of (40) is given by 𝐿𝑢,𝑝,𝑞,𝜆₁,𝜆₂=12𝑢−𝑓²_𝑋+_{1≤𝑖,𝑗≤𝑛}𝐾𝑔₁^𝑘−1∘𝑔₁^𝑘−1_𝑖,𝑗𝑝_𝑖,𝑗+_{1≤𝑖,𝑗≤𝑛}𝐾𝑔₂^𝑘−1∘𝑔₂^𝑘−1_𝑖,𝑗𝑞_𝑖,𝑗+𝜆₁,𝑝−∇𝑢_𝑌+𝜆₂,𝑞−𝐻𝑢_𝑍+𝛽2𝑝−∇𝑢²_𝑌+𝑞−𝐻𝑢²_𝑍, (41) where 𝛽>0 is the penalty parameter for the linear constraints. Then the minimization of (40) is achieved by an iterative process: in each iteration, we minimize the augmented Lagrangian function (41) with respect to 𝑢, 𝑝, and 𝑞, given the other two updated in previous iteration, and then update 𝜆₁ and 𝜆₂. The computational procedure is presented in Algorithm 3.

Input:  𝑓, 𝑔1𝑘−1, 𝑔2𝑘−1, 𝑟, 𝛽, and 𝐿. Output:  𝑢𝑘, 𝑝𝑘, 𝑞𝑘, 𝜆𝑘1, and 𝜆𝑘2. Initialization:  𝑝𝑘,0=𝑝𝑘−1, 𝑞𝑘,0=𝑞𝑘−1, 𝜆1𝑘,0=𝜆1𝑘−1, and 𝜆2𝑘,0=𝜆2𝑘−1.   while not converged and 𝑙<𝐿  do    Step  1. Given 𝑝𝑘,𝑙−1, 𝑞𝑘,𝑙−1, 𝜆1𝑘,𝑙−1, and 𝜆2𝑘,𝑙−1, computing 𝑢𝑘,𝑙 by       𝑢𝑘,𝑙=ℱ−1ℱ𝑓−div𝜆1𝑘,𝑙−1+𝛽𝑝𝑘,𝑙−1+𝐻∗𝜆2𝑘,𝑙−1+𝛽𝑞𝑘,𝑙−11−𝛽ℱ(Δ)+𝛽ℱ(𝐻∗𝐻),        (a)    where ℱ and ℱ−1 denote the discrete Fourier transform and the inverse discrete Fourier    transform, respectively, and Fourier transforms of operators Δ and 𝐻∗𝐻 are regarded as the    transforms of their corresponding convolution kernels.    Step  2. Given 𝑢𝑘,𝑙, 𝜆1𝑘,𝑙−1, and 𝜆2𝑘,𝑙−1, update 𝑝𝑘,𝑙 and 𝑞𝑘,𝑙 by the two-dimensional shrinkage       𝑝𝑘,𝑙𝑖,𝑗=max∇𝑢𝑘,𝑙−𝜆1𝑘,𝑙−1𝛽𝑖,𝑗−𝐾𝑔1𝑘−1∘𝑔1𝑘−1𝑖,𝑗𝛽,0∇𝑢𝑘,𝑙−𝜆1𝑘,𝑙−1/𝛽𝑖,𝑗∇𝑢𝑘,𝑙−𝜆1𝑘,𝑙−1/𝛽𝑖,𝑗,      (b)       𝑞𝑘,𝑙𝑖,𝑗=max𝐻𝑢𝑘,𝑙−𝜆2𝑘,𝑙−1𝛽𝑖,𝑗−𝐾𝑔2𝑘−1∘𝑔2𝑘−1𝑖,𝑗𝛽,0𝐻𝑢𝑘,𝑙−𝜆2𝑘,𝑙−1/𝛽𝑖,𝑗𝐻𝑢𝑘,𝑙−𝜆2𝑘,𝑙−1/𝛽𝑖,𝑗,    (c)    respectively, where 0⋅(0/0)=0 is assumed.    Step  3. Given 𝑢𝑘,𝑙, 𝑝𝑘,𝑙, and 𝑞𝑘,𝑙, update 𝜆1𝑘,𝑙 and 𝜆2𝑘,𝑙 by          𝜆1𝑘,𝑙=𝜆1𝑘,𝑙−1+𝛽𝑝𝑘,𝑙−∇𝑢𝑘,𝑙,               (d)          𝜆2𝑘,𝑙=𝜆2𝑘,𝑙−1+𝛽𝑞𝑘,𝑙−𝐻𝑢𝑘,𝑙,               (e)    respectively. end while

Algorithm 3 

The alternating minimization method for solving (41).

Now we make some remarks on the ADM procedure. First, we observe that every step in Algorithm 3 has a closed-form solution, and then the alternating direction method can be efficiently implemented. Moreover, for a convex objective function with linear constraints, the convergence of the alternating direction method is guaranteed; see, for example, [28, 29]. Second, to obtain an exact solution of (a) in Algorithm 2, one needs to let 𝐿 (the maximum iteration number of the ADM procedure) tend to infinity; that is, 𝑢^𝑘=𝑢^𝑘,∞. Since the regularization parameters are not optimal, in practice, it is not necessary to compute 𝑢^𝑘 exactly. For the sake of computational efficiency, only several ADM iterations are implemented. Third, aiming at faster convergence of the ADM procedure, we initialize each iteration with the auxiliary variables and the Lagrangian multipliers updated in the previous iteration. Numerical examples will be given in Section 5.2 to demonstrate the efficiency of our numerical scheme.

5. Numerical Experiments

In this section, we provide numerical results to illustrate the effectiveness of the proposed models. In Section 5.1, we compare the performance of the proposed models with the ROF model (3), the LLT model (4), and the hybrid model (5). In Section 5.2, we give some criterions on choosing the parameters in our algorithms. All the experiments are performed under Windows 7 and MATLAB R2010a (Version 7.10.0.499) running on a desktop with an Intel Core i3-2130 CPU at 3.40 GHz and 4 GB memory.

5.1. Comparative Results

Four test images shown in Figure 1 are considered for the comparative experiment. The intensity range of the original images is scaled to [0,1]. The degraded images are corrupted by white Gaussian noise with the noise level 0.1, which are shown in Figures 3(a), 4(a), 5(a), and 6(a). The quality of the restored images is assessed using the relative error (ReErr) and the signal-to-noise ratio (SNR). They are defined as ReErr=𝑢−𝑢_𝑐F‖𝑢‖_F,SNR=10log₁₀‖𝑢−𝑢‖_F𝑢−𝑢_𝑐F, (42) where 𝑢, 𝑢, and 𝑢_𝑐 are the original image, the mean intensity value of 𝑢, and the restored image, respectively.

Now we present implementation details of the comparative experiment. In our algorithms, parameters 𝑀 and 𝐿 are fixed and set to be 0.07 and 10, respectively. For parameters 𝛼 and 𝑟, their values are reported in the corresponding figures. For solving models (3), (4), and (5), we use the alternating direction method (ADM) [22]. For a fair comparison, all of the parameters in these three models are optimized to achieve the best restoration with respect to the ReErr and the SNR values. In all the ADM procedures, we set the penalty parameter 𝛽 to be 1 as a default value for most of digital images with intensity range in [0,1]. We terminate all the algorithms by the following stopping criterion: 𝑢_𝑖−1−𝑢_𝑖F𝑢_𝑖−1F≤10⁻⁴. (43)

In Figures 3, 4, 5, and 6, we display the resulting images restored by different models. The corresponding relative errors, SNR values (dB), and computational time (in seconds) are listed in Table 1. From these results, we find that our models restore better images than the ROF model (3), the LLT model (4), and the hybrid model (5), both visually and quantitatively. In the images restored by the ROF model, we observe that sharp edges are preserved, but the staircase effect is clearly present (e.g., the homogeneous regions in Figures 5 and 6). For the restored results of the LLT model, the staircase effect is suppressed, but some speckle effect is introduced in homogeneous regions (e.g., the background of Figures 3, 4, 5, and 6). We also note that both models (3) and (4) compromise details and textures due to using the global constant regularization parameter. For the images restored by the hybrid model, visual artifacts are almost suppressed, but details are not recovered well (e.g., the ramp region in Figure 3). This is due to the inaccurate estimation of their weighting function from the observed image. Our methods, on the other hand, restore the homogeneous regions significantly better while preserving more details. To better understand the behavior of our methods, some zoomed-in local results and contour plots are shown in Figures 7, 8, and 9. We can find that the homogeneous regions are almost smooth (e.g., the background of Figure 3, the flat and ramp regions in Figure 7, and the face of Barbara in Figure 9). The contour plots given in Figure 8 also visually illustrate the above observation. At the same time, details and textures are preserved clearly without being overregularized (e.g., the ramp region in Figure 3, the textures in Figure 7, and the textures on the scarf in Figure 9). Furthermore, we find that Model 2 performs better than Model 1 with respect to preserving edges (e.g., the contours around the lips and nose of Lena in Figure 8). In Figure 10, we show the final values of the regularization parameters in our models compared with the weighting functions in the hybrid model. It is clear that Figures 10(a) and 10(e) are noisy and inaccurate, whereas Figures 10(b)–10(d) and 10(f)–10(h) are less affected by noise. Moreover, we observe that the values of the regularization parameters are small in detail and texture regions, and they are large in homogeneous regions. Therefore, our methods are superior with respect to removing noise while preserving details.

Table 1 

Summarized results of the comparative experiment.

5.2. Parameter Study

In this section, we present some criterions on choosing the parameters which are necessary to start up our algorithms. In our experiments, the parameter settings are the same as those in Section 5.1, except the testing parameters.

The window size 𝑟 is used to control the smoothness of the spatially dependent regularization parameter. We test our algorithms for different values of 𝑟 varying from 1 to 19. Figure 11 shows the plots of the SNR values of the restored images. We see from the plots that a small value of 𝑟 (leads to a sharp regularization parameter) is suitable for images with sharp features (e.g., 𝑟=1 for Figure 2(a) and 𝑟=3 for Figure 2(c)), whereas a moderate value of 𝑟 (leads to a relatively smooth regularization parameter) is needed for images containing textures (e.g., 𝑟=7 for Figures 2(b) and 2(d)). For large values of 𝑟 (𝑟≥9), we find that the changes of the SNR values are not significant. However, if 𝑟 becomes too large, then the regularization parameter is oversmooth, which compromises image details.

(a)

(b)

(c)

(d)

Figure 2 

Original images: (a) synthetic image 1 (128×128), (b) synthetic image 2 (256×256), (c) “Lena” (256×256), and (d) “Barbara” (256×256).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Results of different models for the synthetic image 1: (a) the noisy image, the restored images by (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1 (𝛼=0.1 and 𝑟=1), and (f) our Model 2 (𝛼=0.18 and 𝑟=1).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Results of different models for the synthetic image 2: (a) the noisy image, the restored images by (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1 (𝛼=0.21 and 𝑟=7), and (f) our Model 2 (𝛼=0.35 and 𝑟=7).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5 

Results of different models for the image “Lena”: (a) the noisy image, the restored images by (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1 (𝛼=0.26 and 𝑟=3), and (f) our Model 2 (𝛼=0.36 and 𝑟=3).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Results of different models for the image “Barbara”: (a) the noisy image, the restored images by (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1 (𝛼=0.28 and 𝑟=7), and (f) Model 2 (𝛼=0.39 and 𝑟=7).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Zoomed-in local results corresponding to Figures 4(b)–4(f): (a) the original image, results of (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1, and (f) our Model 2.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8 

Zoomed-in local results corresponding to Figures 6(b)–6(f): (a) the original image, results of (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1, and (f) our Model 2.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 

Contour plots corresponding to Figures 5(b)–5(f): (a) the ideal contour plot, results of (b) the ROF model, (c) the LLT model, (d) the hybrid model, (e) our Model 1, and (f) our Model 2.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 10 

Final values of the regularization parameters in the proposed models and weighting functions in the hybrid model corresponding to Figures 4(d)–4(f) and 6(d)–6(f). (a) 𝑔 in the hybrid model, (b) 𝑔∗ in Model 1, (c) 𝑔∗1 in Model 2, (d) 𝑔∗2 in Model 2, (e) 𝑔 in the hybrid model, (f) 𝑔∗ in Model 1, (g) 𝑔∗1 in Model 2, and (h) 𝑔∗2 in Model 2.

(a)

(b)

Figure 11 

SNR for results restored by our models with different window size 𝑟: (a) Model 1 and (b) Model 2.

The parameters 𝑀 and 𝛼 jointly measure the values of the spatially dependent regularization parameter. The difficulty of tuning these parameters is that they depend not only on the noise level but also on the images. For this reason, trial by error for the two parameters is used. More precisely, we simplify the tuning process by fixing 𝑀 and searching for the best 𝛼 of each image.

We also study the behavior of the proposed numerical scheme with respect to the setting of 𝐿 and our initialization. For the setting of Figure 2(c), Table 2 shows the summarized results of our models for different values of 𝐿 with our initialization and zero initialization (the auxiliary variables and the Lagrangian multipliers are initialized to zero). We see from Table 2 that when using zero initialization, only sufficiently large 𝐿 (e.g., 𝐿=30 and 40) can provide good restoration result with increasing computational time. Our initialization, on the other hand, produces results as effectively as zero initialization with slightly less computation time. In addition, there is no considerable difference between the relative errors and the SNR values for different values of 𝐿. Since a too small or too large 𝐿 yields comparatively high computational cost, a moderate 𝐿 (e.g., 𝐿=10) is recommended in our algorithms.

Table 2 

Summarized results of the proposed models for different values of 𝐿 with our initialization and zero initialization.

6. Concluding Remarks

We have proposed two regularization models for image denoising, in which both the restored image and the spatially dependent regularization parameter are estimated simultaneously. By the construction of our models, both variables mutually benefit from each other during the denoising process, which can achieve joint optimality. The numerical behavior of the regularization parameter and the strengths of our models have been discussed in detail. The proposed models, which are nonconvex, can be asymptotically solved by the proposed alternating minimization scheme. Finally, the numerical results indicated that our methods outperform several popular methods with respect to both noise removal and detail preservation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research is supported by 973 Program (2013CB329404), NSFC (61170311), Chinese Universities Specialized Research Fund for the Doctoral Program (20110185110020), and Sichuan Province Sci. and Tech. Research Project (2012GZX0080).