Hybrid regularizers-based adaptive anisotropic diffusion for image denoising

To eliminate the staircasing effect for total variation filter and synchronously avoid the edges blurring for fourth-order PDE filter, a hybrid regularizers-based adaptive anisotropic diffusion is proposed for image denoising. In the proposed model, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^{-1}$$\end{document}H-1-norm is considered as the fidelity term and the regularization term is composed of a total variation regularization and a fourth-order filter. The two filters can be adaptively selected according to the diffusion function. When the pixels locate at the edges, the total variation filter is selected to filter the image, which can preserve the edges. When the pixels belong to the flat regions, the fourth-order filter is adopted to smooth the image, which can eliminate the staircase artifacts. In addition, the split Bregman and relaxation approach are employed in our numerical algorithm to speed up the computation. Experimental results demonstrate that our proposed model outperforms the state-of-the-art models cited in the paper in both the qualitative and quantitative evaluations.


Introduction
With the popularity of image sensor, digital images play a key role in people's daily life. Unfortunately, images are ineluctably contaminated by noise during acquisition, transmission, and storage. Therefore, image denoising is still an open and complex problem in image processing and computer vision (Chatterjee and Milanfar 2010). Image denoising aims to recovering the original image u from the observed noisy image u 0 , where u 0 = u + n, and n is the zero-mean Gaussian white noise with standard deviation σ.
During the past three decades, lots of approaches for removing noise have been developed from linear models to nonlinear models. Linear models perform well in the smooth area. However, they don't preserve edges and corners. To overcome the disadvantages of the linear denoising models, nonlinear denoising models have been developed which have a good balance between noise removal and edge-preserving. Nonlinear models based on variation (Rudin et al. 1992) and partial differential equation (PDE) (Perona and Malik 1990) have been widely used for image denoising. The best known variational denoising model is the total variation (TV) model proposed by Rudin et al. (1992), which minimizes the following equation, where ⊆ R 2 is a bounded open domain with Lipschitz boundary, ∇ denotes the gradient operator, |∇u| is the TV regularization term, (u − u 0 ) 2 is the fidelity term, > 0 is the regularization parameter, which measures the trade off between the regularization term and the fidelity term. The classical TV model is efficient for removing noise and preserving the edges. However, it possesses some undesirable properties in the recovered image under some circumstances, such as the staircasing effect. To overcome the deficiency of the original TV model (Strong 1997), developed the adaptive TV regularization based variational model as, where g(x) is an adaptive edge-stopping function, which is defined in Strong (1997) as follow, where K > 0 is a threshold parameter for balancing the noise removal and edge preservation, and G ρ (x) is the Gaussian filter with standard deviation ρ. Seen from (3), g(x) is smaller near the edges and larger away from the boundaries, so the model of (2) has the capability of preserving the edges while removing noise because the diffusion is stopped across edges.
In addition, Nikolova replaced the ℓ 2 -norm with the ℓ 1 -norm in the fidelity term of TV model in Nikolova (2002). Osher et al. (2005) proposed an iterative regularization method of TV model. Chen et al. (2010) presented an adaptive total variation method based on the difference curvature. Wang et al. (2011) put forward a modified TV model.
Numerical experiments demonstrate that the models mentioned above have good performance in the terms of the trade-off between removing noise and preserving the edges. Unfortunately, the staircasing effects appear in the recovered image owing to using the TV-norm as the regularization term. To overcome this shortcoming, high-order PDE filters have been proposed and applied for image denoising successfully (Lysaker et al. 2003;Liu et al. 2011). One of the most classical fourth-order PDEs (LLT) is introduced by Lysaker et al. (2003) where ∇ 2 denotes the Laplacian operator. However, a major challenge is that higherorder PDEs blur the edges during image denoising.
To make use of the advantages of both TV filter and high-order PDE filters, some hybrid regularization models are recently proposed, which combined the second-order partial differential equations and the fourth-order partial differential equations (Oh et al. 2013). Li et al. (2007) proposed the adaptive image denoising model based on hybrid regularizers combining the advantages of TV model and LLT model as follows, where g(x) also denotes the edge-stopping function defined as in (3). The results of experiments indicate that the model of (5) performs better than the pure second-order or hight-order models.
In recent years, efficient computational algorithms for solving the denoising models have emerged in large numbers, for instance, fixed point iteration, gradient descent methods, primal-dual methods, relaxation methods, Bregman iteration and split Bregman method, and so on. These methods are efficient for image denoising while preserving the edges.
Inspired by Li et al. (2007) and Liu (2015), we propose a novel adaptive anisotropic diffusion model, incorporating the advantages of the total variation filter and the fourthorder filter, and develop an efficient computational algorithm. The main contributions of our paper can be generalized as follows. First of all, the hybrid regularization term of the novel model is composed of total variation regularization and a fourth-order filter. The fidelity term uses the H −1 -norm as opposed to the more commonly used ℓ 1 -norm or ℓ 2 -norm. The two above-mentioned filters can be adaptively selected according to the diffusion function. When the pixels locate at the edges, the total variation filter is selected to filter the image, which can preserve the edges. When the pixels belong to the flat regions, the fourth-order filter is adopted to smooth the image, which can eliminate the staircase artifacts. Another main contribution is that the split Bregman and relaxation approach are successively employed in our numerical algorithm to speed up the computation. Experimental results demonstrate that our proposed model achieves higher quality in both the qualitative and quantitative aspects than that of the state-of-the-art models cited in the paper.
The remainder of this paper is organized as follows. In "Preliminaries" section, we give some definitions. In "The new model and algorithms" section, we give the proposed model and numerical implementation in detail. The experimental results are given in "Experiments" section. Finally, this paper is concluded in the fifth section.

Preliminaries
In this section, we give a brief overview of some necessary notations and definitions for the proposed model, which will be used in the subsequent sections.
Definition 1 (Chen and Wunderli 2002). Let be an open bounded subset of R n (n ≥ 2) with Lipschitz boundary. Given a function u ∈ L 1 (�). Then the total variation of u in is defined as, where div is the divergence operator, C 1 c (�, R n ) is the subset of continuously differentiable vector functions of compact support contained in , and L ∞ (�) is the essential supremum norm.

The new model and algorithms
The proposed model Meyer analyzed that there exists no oscillation function in the space L 2 (�) and a weaker H −1 -norm is appropriate to represent textured or oscillatory patterns (Meyer 2001), so that we replace ℓ 2 -norm of the fidelity term (u 0 − u) with H −1 -norm. Therefore, a novel adaptive image denoising model is proposed, where u and u 0 are the recovered image and the noisy image, respectively. Seen from Eq. (12), the H −1 -norm is considered as the fidelity term and the regularization term is composed of a total variation regularization and a fourth-order filter in the proposed model.
where Gaussian filter G ρ (x) pre-smooths the noisy image. The larger standard deviation σ of the noise is, the larger standard deviation ρ of Gaussian filter is. We set ρ = Cσ, where C lies between 0 to 1. When g(x) → 0, it means that the pixels locate at the edges. Then total variation filter is selected to filter the image, which can preserve the edges. When g(x) → 1, it means that the pixels belong to the flat regions. Then the fourthorder filter is adopted to smooth the image, which can eliminate the staircase artifacts. We replace g(x) with g in the next part of this article. Figure 1 shows the results of image denoising by our proposed model and model from Li et al. (2007), which demonstrates that the model whose fidelity term uses the H −1 -norm yields better results in image denoising since H −1 -norm is appropriate to represent textured or oscillatory patterns.

The numerical algorithm for the proposed model
We apply a split Bregman method (Cai et al. 2009) to solve Eq. (12). The idea of split Bregman method is to use splitting operator and Bregman iteration to solve various inverse problems (Goldstein and Osher 2009).
We turn Eq. (12) into the following constrained minimization problem by introducing an auxiliary variable z, The method of solving the constrained minimization problem is that it may be transformed into the unconstrained minimization problem, so the constrained problem (14) can be turned into the following unconstrained problem by introducing an auxiliary variable b, By taking advantage of split Bregman method, Eq. (15) can be solved iteratively according to the following equations, where k is the number of iterations.

Solve the first subproblem in Eq. (16)
At present, the Euler-Lagrange equation method is usually used to solve the problem similarity to the first subproblem in Eq. (16). However, it works slowly. To accelerate the computation speed, the split Bregman algorithm and relaxation algorithm are adopted to solve the first subproblem in Eq. (16). First, we define |∇u| = �∇ x u� 1 + �∇ y u� 1 , and |∇ 2 u| = � x u� 1 + � y u� 1 , and then the first subproblem in Eq. (16) can be rewritten as follows,  Li et al. (2007), f noise, detecting by model from Li et al. (2007) where ∇ x , ∇ y , x and y are the first-order difference operators and the second-order difference operators, respectively. All the difference operators are approximated using following formulas: where M × N represents the image size. Second, we introduce four auxiliary variables υ x , υ y , ω x , and ω y , and then Eq. (17) can be transformed into the following constrained optimization problem, with υ x = ∇ x u, υ y = ∇ y u, ω x = � x u, and ω y = � y u.
The above constrained problem (22) are turned into the unconstrained minimization problem, where the parameters α > 0 and , and apply the split Bregman method, Eq. (23) can be solved by following equations, with the update equations, f k+1 According to the relaxation algorithm (Jia et al. 2011), we may define, So, we have where ∇ T x , ∇ T y , T x and T y are respectively the adjoint operators of ∇ x , ∇ y , x and y . ∇ T x and ∇ T y have the following discrete forms, Definitely, T x = x and T y = y .

Solve the second subproblem in Eq. (16)
For the second subproblem in Eq. (16), we derive the Euler-Lagrange equation with respect to z, which is as follows, This is a linear equation, so additional operator split (AOS) iteration and Gauss-Seidel (GS) iteration can be used to solve Eq. (30). We use AOS iteration to solve this equation.
In summary, the proposed algorithm for image denoising can be described as follows,

Experiments
In this section, we experimentally compare our proposed model with the state-of-the-art models. All experiments are performed under Matlab R2009a on a PC with an Intel CPU of 1.7 GHz and 4 GB memory. Six grayscale images viz. "Manmade", "Lena", "Peppers", "Barbara", "Cameraman", and "House" are selected as testing examples for both qualitative and quantitative evaluations. The original test images are shown in Fig. 2. The performances of all methods are compared quantitatively by using the peak signal to noise ratio (PSNR), structural similarity index measure (SSIM) (Wang et al. 2004), multi-scale structural similarity index (MS-SSIM) (Wang et al. 2003), and feature-similarity index (FSIM) (Zhang et al. 2011). In addition, we also compare the computing time and iterations of six models. PSNR is defined as follows, (31) PSNR = 10 × log 10 255 2 MSE (db), where u and ū are respectively the recovered image and the original image. Generally, the larger the value of the PSNR, the better the performance. However, PSNR is inconsistent with human visual judgments. SSIM, MS-SSIM, and FSIM are close to the human vision system, so we also use them to assess the noise removal quality. SSIM is defined by, where µ u and σ 2 u are the mean and variance of u, respectively, σ uū is the covariance of u and ū, and c 1 and c 2 are two constants to avoid instability. MS-SSIM is defined by, where the luminance distortion l i (u,ū), the contrast distortion c i (u,ū) and the structure distortion s i (u,ū) at scale i between images u and ū are defined as follows, where µ u and µū represent the mean intensity of u and ū at scale i; σ u (resp. σū) is the standard deviation of u (and ū) at scale i, and σ uū is the covariance between u and ū at scale i. c 1 , c 2 and c 3 are three small constants to avoid instability. In this paper, the values of the exponents α M , β i and γ i are set as the same as those in Wang et al. (2003). FSIM is defined by, where S L (x) at each location x is the similarity measure, which is defined as product of the similarity function S PC (x) on Phase Congruency (PC) and similarity function S G (x) on Gradient Magnitude (GM). S PC (x) and S G (x) are defined as follows, where PC u and PCū denote the PC maps extracted from u and ū, respectively, and G u and Gū denote the GM maps extracted from u and ū, respectively; T 1 and T 2 are two small positive constants to avoid instability. et al. SpringerPlus (2016) 5:404 The termination condition for all experiments is defined as follows, where u n and u n+1 are respectively denoising results at nth and (n + 1)th iteration, and ε is a given positive number. We set ε = 10 −3 in the experiments. Figure 3 shows the results for the 8-bit gray-scale synthetic image with size 320 × 320 pixels, which is corrupted by zero-mean Gaussian white noise with σ = 30. In our experiments, we use the trial-and-error method for determining the optimal parameters. We set µ = 0.3, α = β = 0.15, t = 0.2, K = 0.005, and ρ = 1 in our algorithm, while all parameter values in TV model Rudin et al. (1992), LLT model Lysaker et al. (2003), non-local means (NLM) model Buades et al. (2005), BLS-GSM Portilla et al.   2003), and hybrid model Liu (2015) are chosen manually by trial-and-error method to ensure the best results. Figure 3a, b are the original and noisy images, respectively. Figure 3c-h show that the denoising results of TV model, LLT model, NLM model, BLS-GSM, hybrid model, and our proposed model, respectively. Figure 3c indicates that there exists the staircase effects in the recovered images by TV model. Although LLT model has the advantage on relieving the staircase effect, the edges are blurred and there are serious speckles in the recovered image. The computational efficiency of NLM model is very low. There exits the edges blurring in the result of BLS-GSM from Fig. 3f. We also find that there still exits a little stair-case effect by hybrid model from Fig. 3g. Our proposed model relieves the staircasing effects and avoids the edges blurring. Table 1 shows the PSNR, SSIM, MS-SSIM, and FSIM values corresponding to Fig. 3. From Table 1 and Fig. 3, it is shown that our proposed model produces the best result. At the same time, the proposed method takes less computational time than LLT model, NLM model and hybrid model, but the computational time of our proposed method is slightly inferior to that of TV model and BLS-GSM.
In order to demonstrate that our model can also work well for natural images, the next experiments are conducted for different images corrupted by zero-mean Gaussian white noise with σ = 30. The experimental results are shown in Figs. 4, 5, 6, 7 and 8, where the experimental results of TV model, LLT model, NLM model, BLS-GSM, hybrid model, and our proposed model are illustrated, respectively. The figures show that our model produces the visually most appealing results among the six models. The quantitative PSNR, SSIM, MS-SSIM, and FSIM values are presented in Fig. 9, which depicts that the performance of our proposed model are better than those of TV model, LLT model, NLM model, BLS-GSM, and hybrid model. In order to verify the better performance of our proposed method, Fig. 10 shows the enlarged regions cropped from Fig. 4.
We also use six images corrupted with different levels of Gaussian noise to examine the performance of the proposed model and the alternative models. Tables 2, 3, 4 and 5 give the PSNR, SSIM, MS-SSIM, and FSIM values obtained by the proposed model and the alternative models, respectively. From Tables 2, 3, 4 and 5, it can be observed that our model is greater than or equal to the other five models in PSNR, SSIM, MS-SSIM, and FSIM for the same standard deviation, which demonstrate that our model provides the best noise removal performance at different noise level.

Conclusions
To eliminate the so-called staircase effect in total variation filter and avoid the edges blurring for fourth-order PDE filter, we propose an adaptive anisotropic diffusion model for image denoising, which is composed of a hybrid regularization term combining a total variation filter and a fourth-order filter and the fidelity term using the H −1 -norm. We also develop an efficient algorithm to solve our proposed model. Numerical