A Fast High-Order Total Variation Minimization Method for Multiplicative Noise Removal

Abstract

Multiplicative noise removal problem has received considerable attention in recent years. The total variation regularization method for the solution of the noise removal problem can preserve edges well but has the sometimes undesirable staircase effect. In this paper, we propose a fast high-order total variation minimization method to restore multiplicative noisy images. The proposed method is able to preserve edges and at the same time avoid the staircase effect in the smooth regions. An alternating minimization algorithm is employed to solve the proposed high-order total variation minimization problem. We discuss the convergence of the alternating minimization algorithm. Some numerical results show that the proposed method gives restored images of higher quality than some existing multiplicative noise removal methods.

1. Introduction

Image denoising problem has been widely studied in the areas of image processing. The problem includes additive noise removal and multiplicative noise removal. Most of the literature deals with the additive noise model [1–4]. Given a noisy image 𝑔=𝑢+𝑣, where 𝑢 is the original image and 𝑣 is the noise, the denoising problem is to recover 𝑢 from the observed image 𝑔. The additive noise model has been extensively studied. We refer the reader to [5–11] and references therein for a review of various methods.

In this paper, we consider the problem of seeking the original image from a multiplicative noisy image. In the multiplicative noise model, the recorded image 𝑔 is the multiplication of the original image 𝑢 and the noise 𝑣: 𝑔=𝑢𝑣. (1) This problem arises in ultrasound imaging, synthetic aperture radar (SAR) and sonar (SAS), laser imaging, and magnetic field inhomogeneity in MRI [12, 13]. In this work, we concentrate on the assumption that 𝑣 follows a Gamma distribution, which commonly occurs in SAR. With respect to a given resolution cell of the imaging device, the complex SAR image is computed after the SAR system receives the coherent sum of reflected monochromatic microwaves. The complex SAR image at one point is usually modeled by a complex zero-mean circular Gaussian density (i.e., the real and imaginary parts are independent Gaussian variable with a common variance). The intensity image is defined as the square of the magnitude of the corresponding complex SAR image; its variables follow exponentially distribution independently [14, 15]. Since the complex observation is acquired in highly coherent system, the intensity image occurs to have peculiar granular appearance (known as speckle). The conventional SAR system reduces the speckle effect of the intensity image by the multilooking process (𝐿-look, in the case of 𝐿-looks), that is, averaging the observed intensity images from slightly different angles of the same resolution cell [16].

Due to the coherent nature of these image acquisition processes, nearly all the information of the original image may vanish when it is distorted by the multiplicative noise. In general, the standard additive denoising method, so prevalent in image processing, is inadequate. Therefore, it is necessary to devise efficient and reliable algorithms for recovering the true image from the observed multiplicative noisy image. Until the past decade, a few variational approaches devoted to multiplicative noise removal have been proposed. The early variational approach for multiplicative noise removal is the one by Rudin et al. [17]. According to the statistical properties of the multiplicative noise 𝑣, the recovery of the image 𝑢 was based on solving the following constrained optimization problem [17]: min_𝑢Ω|𝐷𝑢|𝑑𝑥𝑑𝑦subjectto:_Ω𝑔𝑢𝑑𝑥𝑑𝑦=1,_Ω𝑔𝑢−1²𝑑𝑥𝑑𝑦=𝜎², (2) where ∫_Ω|𝐷𝑢|𝑑𝑥𝑑𝑦 is denoted by the seminorm on boundary variation (BV) space that coincides with ∫_Ω𝑢²_𝑥+𝑢²_𝑦𝑑𝑥𝑑𝑦 when 𝑢 is smooth. The two constraints state that the mean of the noise is equal to 1, and the variance is equal to 𝜎². In (2), only basic statistical properties, the mean, and the variance of the noise 𝑣 are considered, which somehow limits the restored results.

By using a maximum a posteriori (MAP) estimator, Aubert and Aujol [18] proposed a function whose minimizer corresponds to the denoised image to be recovered. This function is min_𝑢Ω𝑔log𝑢+𝑢𝑑𝑥𝑑𝑦+𝛼_Ω|𝐷𝑢|𝑑𝑥𝑑𝑦, (3) where the total variation of 𝑢 is utilized as the regularization term, and 𝛼 is the regularization parameter which controls the trade-off between a good fit of 𝑔 and a smoothness requirement due to the total variation regularization. Although the functional they proposed is not convex, they still proved the existence of a minimizer and showed the capability of their model on some numerical examples.

As a result of the drawback of the objective function (3) that is not convex for all 𝑢, the obtained solution is likely not the global optimal solution of (3). To overcome this problem, Huang et al. [19] proposed and studied a strictly convex objective function for multiplicative noise removal problem. In [19], the authors introduced an auxiliary variable 𝑧=log𝑢 in (3). With the new variable, the proposed unconstrained total variation denoising problem is described as follows: min_𝑧,𝑤Ω𝑧+𝑔𝑒^−𝑧𝑑𝑥𝑑𝑦+𝛼₁‖𝑧−𝑤‖²_𝐿2+𝛼_2Ω|𝐷𝑤|𝑑𝑥𝑑𝑦, (4) where 𝛼₁ and 𝛼₂ are positive regularization parameters. The parameter 𝛼₁ measures the trade-off between an image obtained by a maximum likelihood estimation from the first term and a total variation denoised image 𝑤. The parameter 𝛼₂ measures the amount of regularization to a denoised image 𝑤. The main advantage of the proposed method is that the total variation norm can be used in the noise removal process in an efficient manner. Therefore the method proposed in [19] has the ability to preserve edges very well in the denoised image.

Recently, Bioucas-Dias and Figueiredo [20] proposed an efficient multiplicative noise removal method by using variable splitting and constrained optimization. After converting the multiplicative noise model into an additive one by taking logarithms, they used variable splitting to obtain an equivalent constrained problem and then solved this optimization problem by using the augmented Lagrangian method. A set of experiments has shown that the proposed method, which they named MIDAL (multiplicative image denoising by augmented Lagrangian), yields state-of-the-art results both in terms of speed and denoising performance. In [21], Steidl and Teuber considered a variational restoration model consisting of the 𝐼-divergence as data fitting term and the total variation seminorm or nonlocal means as regularizer for removing multiplicative noise. They proposed to compute the minimizers of their restoration functionals by using Douglas-Rachford splitting techniques. For a particular splitting, they presented a semi-implicit scheme to solve the involved nonlinear systems of equations and proved its 𝑄-linear convergence.

It is known that the total variation denoising method is a PDE-based technique that preserves edges well but has the sometimes undesirable staircase effect, namely, the transformation of smooth regions (ramps) into piecewise constant regions (stairs) [22–26]. Therefore, the approaches involving the traditional total variation often cause staircase effect since they favor solutions that are piecewise constant. Staircase solutions fail to satisfy the evaluation of visual quality and they can develop false edges that do not exist in the true image. To attenuate the staircase effect, there is a growing interest in the literature for replacing the total variation norm by the high-order total variation norm. The motivation behind such attempt is to restore potentially a wider of images, which comprise more than merely piecewise constant regions. The majority of the high-order norms involve second-order differential operators because piecewise-vanishing second-order derivatives lead to piecewise-linear solutions that better fit smooth intensity changes; see [27] for more details. There are two classes of second-order regularization methods for image restoration problems. The first class employs a second-order regularizer in a standalone way. For example, in [28], the authors considered a fourth-order partial differential equations (PDE) model for noise removal and employed the dual algorithm of Chambolle for solving the high-order problems. In [29], Hu and Jacob applied higher degree total variation (HDTV) regularization for image recovery. The second class is to combine the TV norm with a second-order regularizer. For example, a technique in [30, 31] combining the TV filter with a fourth-order PDE filter was proposed to preserve edges and at the same time avoid the staircase effect in smooth regions for noise removal. In [32], Papafitsoros and Schönlieb considered a high-order model involving convex functions of the TV and the TV of the first derivatives for image restoration problems and used the split Bregman method to solve numerically the corresponding discretized problem.

In this paper, we present a fast high-order total variation minimization method for multiplicative noise removal problem. This technique substantially reduces the staircase effect, while preserving sharp jump discontinuities (edges). An alternating minimization algorithm is employed to solve the proposed high-order total variation minimization model. We discuss the convergence of the proposed alternating minimization algorithm. Some numerical results show that the proposed method gives restored images of higher quality than some existing multiplicative noise removal methods.

The organization of this paper is outlined as follows. In the next section, we present the high-order total variation minimization method for multiplicative noise removal problem. In Section 3, an alternating minimization algorithm is employed to find the minimizer of the proposed minimization problem. We analyze the convergence of the proposed alternating minimization algorithm in Section 4. Some numerical experiments are given to illustrate the performance of the proposed algorithm in Section 5. Concluding remarks are given in the last section.

2. A High-Order Total Variation Multiplicative Denoising Model

In this section, our aim is to propose a strictly convex objective function for restoring images distorted by the multiplicative noise. The most important is that we apply a high-order total variation to the objective function to achieve sharp edges and staircase reduction efficiently. We introduce our multiplicative denoising model from the statistical perspective using the Bayesian formula.

Let 𝑔, 𝑢, and 𝑣 denote samples of instances of some random variables 𝐹,  𝑈, and 𝑉. Moreover, we assume that the random variable 𝑣 on each pixel is mutually independent and identically distributed (i.i.d). Moreover, the random variable 𝑣 in each pixel follows Gamma distribution; that is, its probability density function is 𝑃_𝑉()=𝑀𝑣;𝐿,𝑀^𝐿𝑣^𝐿−1Γ(𝐿)𝑒^−𝑀𝑣for𝑣>0, (5) where Γ(⋅) is the usual Gamma-function, and 𝐿 and 𝑀 denote the shape scale and inverse parameters in the Gamma distribution, respectively. We note that the mean of a gamma-distributed variable is 𝐿/𝑀 and the variance is 𝐿/𝑀². As the multiplicative noise, in this work we assume that the mean of 𝑣 equals 1, and then we have 𝐿=𝑀.

According to the posteriori estimation, the restored image 𝑢 can be determined by 𝑢=argmax_𝑢𝑃_𝑈∣𝐺𝑢∣𝑔. (6) From the Bayes rule, we have 𝑃_𝑈∣𝐺(𝑢∣𝑔)=𝑃_𝐺∣𝑈(𝑔∣𝑢)𝑃_𝑈(𝑢)/𝑃_𝐺(𝑔). It then follows that 𝑢=argmax_𝑢𝑃_𝐺∣𝑈𝑔∣𝑢𝑃_𝑈(𝑢)𝑃_𝐺𝑔. (7) Using the proposition in [18]: 𝑃_𝑉(𝑔/𝑢)(1/𝑢)=𝑃_𝐺∣𝑈(𝑔∣𝑢), we obtain 𝑃_𝐺∣𝑈=𝐿𝑔∣𝑢^𝐿𝑔^𝐿−1𝑢^𝐿Γ(𝐿)𝑒^{−𝐿𝑔/𝑢}. (8) Taking the logarithm transformation into account, we assume that the image prior 𝑃_𝑈(𝑢) is given by the following: 𝑃_𝑈(𝑢)=𝑃_𝑈∣𝑊(𝑢∣𝑤)𝑃_𝑊(𝑤) (9) with 𝑃_𝑈∣𝑊(𝑢∣𝑤)∝𝑒^{−𝛼1‖log𝑢−𝑤‖22},𝑃_𝑊(𝑤)=𝑒^{−𝛼2‖𝑤‖HTV}, (10)(11) where 𝛼₁ and 𝛼₂ are two positive scalars and ‖𝑤‖_HTV is a high-order total variation of 𝑤. Here we assume that the difference between log𝑢 and 𝑤 follows a Gaussian distribution, and 𝑤 follows a high-order total variation prior. Thus, maximizing 𝑃_𝑈∣𝐺(𝑢∣𝑔) amounts to minimizing the log-likelihood: −log𝑃_𝑈∣𝐺𝑢∣𝑔=−log𝑃_𝐺∣𝑈𝑃𝑔∣𝑢−log_𝑈∣𝑊(𝑢∣𝑤)−log𝑃_𝑊(𝑤)+log𝑃_𝐺𝑔. (12) Since 𝑃_𝐺(𝑔) is constant, (6) can be rewritten as the following optimization problem: 𝑢=argmin_𝑢−log𝑃_𝐺∣𝑈𝑔∣𝑢−log𝑃_𝑈∣𝑊(𝑃𝑢∣𝑤)−log_𝑊(𝑤). (13) The previous computation leads us to propose the following high-order total variation model for restoring images corrupted by the Gamma noise by introducing a new variable 𝑧=log𝑢: min_𝑧,𝑤𝒥(𝑧,𝑤)≐min_𝑧,𝑤Ω𝑧+𝑔𝑒^−𝑧𝑑𝑥𝑑𝑦+𝛼₁‖𝑧−𝑤‖²_𝐿2+𝛼_2Ω𝐷²𝑤𝑑𝑥𝑑𝑦, (14) where 𝛼₁ and 𝛼₂ are positive regularization parameters. In the above model, 𝐷²𝑤 denotes the Hessian of 𝑤, and |𝐷²𝑤|=𝑤²_𝑥𝑥+𝑤²_𝑥𝑦+𝑤²_𝑦𝑥+𝑤²_𝑦𝑦 is just the Frobenius norm of the Hessian 𝐷²𝑤. The parameter 𝛼₁ measures the trade-off between an image obtained by a maximum likelihood estimation from the first term and a high-order total variation denoised image 𝑤. The parameter 𝛼₂ measures the amount of regularization to a denoised image 𝑤. The main advantage of the proposed method is that the high-order total variation norm can be used in the noise removal process in an efficient manner. Therefore, the proposed method has the ability to preserve edges very well and reduce staircase effect in the denoised image.

Moreover, it is obvious that when 𝑢 includes an edge, 𝑧 also contains an edge at the same point; that is, the logarithm transformation preserves image edges; see [19] for more details. Based on this idea, we can view 𝑧 as an image in the logarithm domain. We can apply a high-order total variation regularization to 𝑧 to recover its edges and reduce staircase effect. It is interesting to note that 𝑧 in the proposed model can be positive, zero, or negative, and the corresponding 𝑢=𝑒^𝑧 is still positive, although 𝑢 in the objective function (3) should be positive. Especially, one can find that when 𝑧 in the proposed model is a large negative value, 𝑢 is just close to zero.

3. The Alternating Minimization Algorithm

In this section, an alternating minimization algorithm is employed to solve the proposed high-order total variation minimization problem efficiently. To solve (14), we need to consider the following two minimization subproblems. (i)For 𝑤 fixed, determine the solution of min_𝑧Ω𝑧+𝑔𝑒^−𝑧𝑑𝑥𝑑𝑦+𝛼₁‖𝑧−𝑤‖²_𝐿2. (15) (ii)For 𝑧 fixed, determine the solution of min_𝑤𝛼₁‖𝑧−𝑤‖²_𝐿2+𝛼_2Ω𝐷²𝑤𝑑𝑥𝑑𝑦. (16)In the discrete setting, minimizing the first problem amounts to solve the following optimization problem: min_𝑧^𝑛2₁𝑧_𝑖+𝑔𝑒^−𝑧𝑖+𝛼₁‖𝑧−𝑤‖²₂. (17) Straightforward computations can show that the minimizer of the above problem is the solution of the following nonlinear systems: 1−𝑔𝑒^−𝑧𝑖+2𝛼₁𝑧_𝑖−𝑤_𝑖=0for𝑖=1,2,…,𝑛². (18) There are 𝑛² decoupled nonlinear equations to be solved. The second derivative with respect to 𝑧 of the first term in (17) is equal to 𝑔𝑒^−𝑧, which is always greater than zero. Here we assume that 𝑔>0 in the multiplicative noise model (1). Therefore, the first term of the objective function is strictly convex for all 𝑧. Hence, we know that the objective function (17) is also strictly convex. Therefore, the corresponding nonlinear equation (18) has a unique solution and the solution can be determined by using the Newton method very efficiently.

Subproblem (16) is a high-order total variation regularization process for image denoising. Some effective algorithms for the minimizer of the high-order total variation problem have been proposed in the literature. For example, a lot of PDE-based methods are widely used to preserve edges and reduce the staircase effect for the high-order total variation problem; see [6, 22, 23, 30] for more details. In particular, inspired by the work from Chambolle [5], a dual algorithm and its convergence analysis for minimization of the high-order model are presented in [24, 28]. The algorithm overcomes the numerical difficulties related to the nondifferentiability of the high-order model. Some numerical results show that the convergence of the dual algorithm is faster than the gradient descent time-marching algorithm.

In the present paper, we employ the dual algorithm proposed by [24, 28] for the high-order total variation denoising problem (16). The discrete version of problem (16) is represented as min_𝑤𝛼₁𝛼₂‖𝑧−𝑤‖²₂+‖𝑤‖_HTV. (19) In the dual algorithm, we need to solve the following minimization problem: min_𝑝𝑧−𝛽div²𝑝,subjectto𝑝_𝑗,𝑘≤1,∀1≤𝑗,𝑘≤𝑛, (20) where 𝛽=𝛼₂/2𝛼₁ and 𝑝=^𝑝11𝑝12_𝑝21𝑝22 with 𝑝^ℎ𝑘∈𝑅^𝑛×𝑛 for ℎ,𝑘=1,2. This problem can be solved by a semi-implicit fixed point iteration: 𝑝^(𝑖+1)_𝑗,𝑘=𝑝^(𝑖)_𝑗,𝑘−𝜏𝐴^(𝑖)_𝑗,𝑘1+𝜏𝐴^(𝑖)_𝑗,𝑘, (21) where 𝑝⁰=0, 𝐴^(𝑖)_𝑗,𝑘=(∇²(div²𝑝^(𝑖)−𝑧/𝛽))_𝑗,𝑘, and 𝜏≤1/64 is the step size. The operator ∇² is defined as ∇²𝑤=^{𝑤𝑥𝑥𝑤𝑥𝑦}_{𝑤𝑦𝑥𝑤𝑦𝑦}. Here, div²𝑝^(𝑖)_𝑗,𝑘=𝐷_𝑥𝑥((𝑝¹¹)^(𝑖)_𝑗,𝑘)+𝐷⁻_𝑦𝑥((𝑝¹²)^(𝑖)_𝑗,𝑘)+𝐷⁺_𝑥𝑦((𝑝²¹)^(𝑖)_𝑗,𝑘)+𝐷_𝑦𝑦((𝑝²²)^(𝑖)_𝑗,𝑘); see [24, 28] for more details. After the minimizer 𝑝^∗ of the constrained optimization problem in (20) is determined, the denoised image 𝑤 can be computed by 𝑤=𝑧−𝛽div²𝑝^∗. In addition, it is shown that if 𝜏≤1/64, then div²𝑝^(𝑖) is convergent and 𝑝^(𝑖) with the dual method is the solution of the minimization problem (20) as 𝑖→∞.

Starting from an initial guess 𝑤⁽⁰⁾, the alternating minimization method computes a sequence of iterates 𝑧⁽¹⁾,𝑤⁽¹⁾,𝑧⁽²⁾,𝑤⁽²⁾,…,𝑧^(𝑘),𝑤^(𝑘),… (22) such that 𝑧^(𝑘)≐𝒮_ℎ𝑤^(𝑘−1)=argmin_𝑧^𝑛2₁𝑧_𝑖+𝑔𝑒^−𝑧𝑖+𝛼₁𝑧−𝑤^(𝑘−1)2₂,𝑤^(𝑘)≐𝒮_𝑝𝑧^(𝑘)=argmin_𝑤𝛼₁𝛼₂𝑧^(𝑘)−𝑤²₂+‖𝑤‖_HTV (23) for 𝑘=1,2,….

We are now in a position to describe the alternating minimization algorithm for the multiplicative noise removal problem.

Algorithm 1. The alternating minimization algorithm for multiplicative noise removal:
Choose an initial guess 𝑤⁽⁰⁾, and 𝑘=1.(1) For 𝑤^(𝑘−1) fixed, employ the Newton iterative method to compute 𝑧^(𝑘)=argmin_𝑧^𝑛2₁𝑧_𝑖+𝑔𝑒^−𝑧𝑖+𝛼₁𝑧−𝑤^(𝑘−1)2₂. (24)(2) For 𝑧^(𝑘) fixed, use a fixed point iteration to compute 𝑝^∗=argmin_𝑝𝑧^(𝑘)−𝛽div²𝑝 (25) and obtain 𝑤^(𝑘)=𝑧^(𝑘)−𝛽div²𝑝^∗. (26)(3) Check the stopping criteria. If a stopping criterion is satisfied, then exit with an approximate restored image 𝑢^(𝑘)=𝑒^𝑧(𝑘); otherwise, let 𝑘=𝑘+1 and go to step (1).

4. Convergence Analysis

The main aim of this section is to analyze the convergence of the proposed alternating minimization algorithm. We first remark that 𝒥(𝑧,𝑤) in (14) is strictly convex, as the sum of a convex function and of two strictly convex functions; see [33] for more details. Hence, it suffices to show that there exists a unique minimizer of the objective function 𝒥(𝑧,𝑤). We denote by 𝒥^∗ the minimum value of 𝒥(𝑧,𝑤). In the following we show that our alternating minimization algorithm gives asymptotically the solution of the discrete problem associated with (14). We have the following theorem.

Theorem 2. For any initial guess 𝑤⁽⁰⁾∈ℝ^𝑛2, (𝑧^(𝑘),𝑤^(𝑘)) generated by Algorithm 1 converges to a minimizer of the objective function 𝒥(𝑧,𝑤) as 𝑘→∞.

Proof. It follows from the alternating iterative process in Algorithm 1 that 𝒥𝑧^(𝑘),𝑤^(𝑘)≥𝒥𝑧^(𝑘+1),𝑤^(𝑘)≥𝒥𝑧^(𝑘+1),𝑤^(𝑘+1). (27)
It is obvious that the sequence 𝒥(𝑧^(𝑘),𝑤^(𝑘)) is nonincreasing. We note that it is bounded from below by the minimum value 𝒥^∗. Let 𝜃₁(𝑧)=∑^𝑛2₁(𝑧_𝑖+𝑔𝑒^−𝑧𝑖), 𝜃₂(𝑤)=‖𝑤‖_HTV and 𝐷_Φ(𝑧,𝑤)=𝛼₁‖𝑧−𝑤‖²₂. Since 𝜃₁ and 𝜃₂ are convex, Φ is convex and differentiable and 𝐷_Φ is the associated Bregman distance, we know from [34, 35] that 𝒥(𝑧,𝑤)=𝜃₁(𝑧)+𝜃₂(𝑤)+𝐷_Φ(𝑧,𝑤) has the five-point property. Hence, 𝒥(𝑧^(𝑘),𝑤^(𝑘)) thus converges in ℝ to 𝒥^∗. So we have lim_𝑘→∞𝒥𝑧^(𝑘),𝑤^(𝑘)=𝒥^∗. (28) We note that if the objective function 𝒥(𝑧,𝑤) is coercive, the sequence (𝑧^(𝑘),𝑤^(𝑘)) is bounded since the sequence 𝒥(𝑧^(𝑘),𝑤^(𝑘)) converges.
Now we need to show that the objective function 𝒥(𝑓,𝑢,𝑣) is coercive. Let 𝐷_𝑥𝑥, 𝐷_𝑥𝑦, 𝐷_𝑦𝑥, and 𝐷_𝑦𝑦 be the second order difference matrices in the 𝑥𝑥, 𝑥𝑦, 𝑦𝑥, and 𝑦𝑦 direction, respectively, and 𝐷𝑆=^𝑥𝑥𝐷_𝑥𝑦𝐷_𝑦𝑥𝐷_𝑦𝑦. (29) It is not difficult to obtain the lower bound of the discrete high-order total variation as follows: ‖𝑤‖_HTV=_{1≤𝑗,𝑘≤𝑛}(𝑤_𝑥𝑥)_𝑗,𝑘²+(𝑤_𝑥𝑦)_𝑗,𝑘²+(𝑤_𝑦𝑥)_𝑗,𝑘²+(𝑤_𝑦𝑦)_𝑗,𝑘²≥12_{1≤𝑗,𝑘≤𝑛}𝑤_{𝑥𝑥𝑗,𝑘}+𝑤_{𝑥𝑦𝑗,𝑘}+𝑤_{𝑦𝑥𝑗,𝑘}+𝑤_{𝑦𝑦𝑗,𝑘}=12‖𝑆𝑤‖₁. (30) Let 𝒥₁(𝑧,𝑤)=𝛼₁‖𝑧−𝑤‖²₂ and 𝒥₂(𝑧,𝑤)=(𝛼₂/2)‖𝑆𝑤‖₁. Denote Ω=(𝑧,𝑤)𝒥₁(𝑧,𝑤)+𝒥₂(𝑧,𝑤)=0. (31) If (𝑧,𝑤)∈Ω, we obtain 𝑧=𝑤. As ‖(𝑧,𝑤)‖₂→∞, it is easy to show that 𝒥(𝑧,𝑤)=^𝑛2_𝑖𝑧_𝑖+𝑔𝑒^−𝑧𝑖⟶∞. (32) On the other hand, we obtain ^𝑛2_𝑖𝑧_𝑖+𝑔𝑒^−𝑧𝑖≥^𝑛2_𝑖log𝑔_𝑖+1 (33) since it is a strictly convex function. By using the above inequality, we have 𝒥(𝑧,𝑤)≥^𝑛2_𝑖log𝑔_𝑖+1+𝒥₁(𝑧,𝑤)+𝒥₂(𝑧,𝑤). (34) Hence, if (𝑧,𝑤)∉Ω with ‖(𝑧,𝑤)‖₂→∞, it is not difficult to obtain that 𝒥₁(𝑧,𝑤)+𝒥₂(𝑧,𝑤)→∞. Thus, we get that 𝒥(𝑧,𝑤) also tends to infinity. It follows from the definition of coercive that 𝒥(𝑧,𝑤) is coercive.
Since the sequence (𝑧^(𝑘),𝑤^(𝑘)) is bounded, we can thus extract a convergent subsequence {(𝑧^(𝑘𝑖),𝑤^(𝑘𝑖))} from {(𝑧^(𝑘),𝑤^(𝑘))} such that lim_𝑖→∞{(𝑧^(𝑘𝑖),𝑤^(𝑘𝑖))}=(𝑧,𝑤). Then we obtain 𝑚=𝒥(𝑧,𝑤). Moreover, we have, for all 𝑘_𝑖∈𝑁 and all 𝑧, 𝐽𝑧,𝑤^(𝑘𝑖)≥𝒥𝑧^(𝑘𝑖+1),𝑧^(𝑘𝑖), (35) and for all 𝑘_𝑖∈ℕ and all 𝑤, 𝒥𝑧^(𝑘𝑖+1),𝑤≥𝒥𝑧^(𝑘𝑖+1),𝑤^(𝑘𝑖+1). (36)
Let us denote by 𝑧^∼ a cluster point of {𝑧^(𝑘𝑖+1)}. We may immediately obtain from (27) that 𝒥^∗=𝒥(𝑧,𝑤)=𝒥𝑧^∼,𝑤. (37) Thus, we have that ∑^𝑛2_𝑖(𝑧_𝑖+𝑔𝑒⁻𝑧_𝑖)+𝛼₁‖𝑧−𝑤‖²₂=∑^𝑛2_𝑖(𝑧^∼_𝑖+𝑔𝑒^{−∼𝑧𝑖})+𝛼₁‖𝑧^∼−𝑤‖²₂. We conclude that 𝑧^∼=𝑧 since 𝑧^∼ and 𝑧 are the minimizers of the strictly convex function 𝒮_ℎ(𝑤) defined in (23). Hence, lim_{𝑘𝑖→∞}𝑧^(𝑘𝑖+1)=𝑧. Following a similar analysis, we can show that lim_{𝑘𝑖→∞}𝑤^(𝑘𝑖+1)=𝑤. By passing to the limit in (35) and (36), we have 𝒥(𝑤)(𝑤)𝒥()(𝑤)𝑧,≥𝒥𝑧,,∀𝑧,𝑧,𝑤≥𝒥𝑧,,∀𝑤. (38) We know that (38) can be rewritten as follows: 𝒥(𝑧,𝑤)=min_𝑧𝒥(𝑤),𝒥(𝑤)𝑧,𝑧,=min_𝑤𝒥(𝑧,𝑤). (39) From the definition of 𝒥(𝑧,𝑤), (39) is equivalent to the following: (𝑤),(𝑤),0∈𝜓𝑧,0∈𝜙𝑧, (40) where 𝜓(𝑧,𝑤)=1−𝑔𝑒^−𝑧+2𝛼₁(𝑧−𝑤) and 𝜙(𝑧,𝑤)=2𝛼₁(𝑧−𝑤)+𝛼₂𝜕(‖𝑤‖_HTV). The subdifferential of 𝒥(𝑧,𝑤) at (𝑧,𝑤) is given by (𝑤)=𝜓(𝑤)𝜙(𝑤)𝜕𝒥𝑧,𝑧,𝑧,. (41) Therefore, we have 00∈𝜕𝒥(𝑧,𝑤), (42) which implies that (𝑧,𝑤) is the minimizer of 𝒥(𝑧,𝑤). Hence the whole sequence 𝒥(𝑧^(𝑘),𝑤^(𝑘)) converges towards 𝒥^∗, the unique minimum of 𝒥(𝑧,𝑤). We deduce that the sequence (𝑧^(𝑘),𝑤^(𝑘)) converges to (𝑧,𝑤), the minimizer of 𝒥(𝑧,𝑤) as 𝑘 tends to infinity.

5. Numerical Experiments

In this section we provide numerical results from multiplicative noise removal problem to demonstrate the performance of the proposed regularization method. All computations of the present paper were carried out in Matlab 7.10. The results were obtained by running the Matlab codes on an Intel Core i3 CPU (2.27 GHz, 2.27 GHz) computer with RAM of 2048 M. The initial guess is chosen to be the noisy image in all tests.

The quality of the restoration results with different methods is compared quantitatively by using the Peak-Signal-to-Noise Ratio (PSNR), the relative error (ReErr) and Structural SIMilarity index (SSIM). They are defined as follows: PSNR=20log₁₀255𝑛𝑢^(𝑘)−𝑢_true2,𝑢ReErr=^(𝑘)−𝑢_true2𝑢_true2, (43) where 𝑢_true and 𝑢^(𝑘) are the ideal image and the restored image of order 𝑛 respectively, and SSIM=2𝜇_𝑢𝜇^∼_𝑢+𝐶₁2𝜎_𝑢^∼_𝑢+𝐶₂𝜇²_𝑢+𝜇²_∼𝑢+𝐶₁𝜎²_𝑢+𝜎²_∼𝑢+𝐶₂, (44) where 𝜇_𝑢 and 𝜇^∼_𝑢 are averages of 𝑢 and 𝑢^∼, respectively. 𝜎_𝑢 and 𝜎^∼_𝑢 are the variance of 𝑢 and 𝑢^∼, respectively. 𝜎_𝑢^∼_𝑢 is the covariance of 𝑢 and 𝑢^∼. The positive constants 𝐶₁ and 𝐶₂ can be thought of as stabilizing constants for near-zero denominator values. In general, a high PSNR-value or a small ReErr-value indicates that the restoration is more accurate. The SSIM is a well-known quality metric used to measure the similarity between two images. This method developed by Wang et al. [36] is based on three specific statistical measures that are much closer to how the human eye perceives differences between images. The whiter SSIM map, the closer two images are. Therefore, we use the SSIM map to reveal areas of high or low similarity between two images in this work.

We note that a straightforward high-order total variation regularization scheme for multiplicative noise removal is to solve the following optimization problem: min_𝑢Ω𝑔log𝑢+𝑢𝑑𝑥𝑑𝑦+𝛼_Ω𝐷²𝑢𝑑𝑥𝑑𝑦, (45) where 𝛼 is a positive regularization parameter which measures the trade-off between a good fit and a regularized solution. The Euler-Lagrange equation for minimization problem (45) amounts to solve (omitting the boundary conditions): 𝑢−𝑔𝛼𝑢²+𝑢_𝑥𝑥𝐷²𝑢_𝑥𝑥+𝑢_𝑥𝑦𝐷²𝑢_𝑦𝑥+𝑢_𝑦𝑥𝐷²𝑢_𝑥𝑦+𝑢_𝑦𝑦𝐷²𝑢_𝑦𝑦=0. (46) Hence, the time-marching (TM) algorithm solving the high-order total variation model is as follows: 𝑢_𝑡=𝑔−𝑢𝛼𝑢²−𝑢_𝑥𝑥𝐷²𝑢+𝜖_𝑥𝑥−𝑢_𝑥𝑦𝐷²𝑢+𝜖_𝑦𝑥−𝑢_𝑦𝑥𝐷²𝑢+𝜖_𝑥𝑦−𝑢_𝑦𝑦𝐷²𝑢+𝜖_𝑦𝑦, (47) where the parameter 𝜖=10⁻⁵ is introduced to avoid divisions by zero.

In the first test, we compare our method with the straightforward high-order total variation regularization scheme. In the time-marching algorithm for solving the straightforward high-order total variation regularization problem, we set the step size as 0.1 in order to obtain a stable iterative procedure. The algorithm is stopped when the maximum number of iterations is reached. In addition, we carry out many experiments with different 𝛼-values in the model (45), the ones with the best results are presented in this work. For the proposed high-order total variation method (HighTV method), we terminate the iterations for the method and accept 𝑢^(𝑘) as the computed approximation of the ideal image as soon as the maximum number of allowed outer iterations has been carried out or the relative differences between consecutive iterates 𝑢⁽¹⁾,𝑢⁽²⁾,𝑢⁽³⁾,… satisfy 𝑢^(𝑘+1)−𝑢^(𝑘)₂𝑢^(𝑘+1)₂<10⁻⁴. (48) In the proposed HighTV regularization method, there are two regularization parameters. We know that the quality of the restored image is highly depended on the regularization parameters. Similar to [19], we fix 𝛼₁=19 in the proposed method in order to reduce the computational time in search for good regularization parameters. We determine the best value of 𝛼₂ such that the relative error of such a restored image with respect to such an ideal image is the smallest.

In the first experiment, we consider the “Cameraman” image of size 256×256. The original image in Figure 1(a) is distorted by the Gamma noise with 𝐿=10 and 𝐿=20. The noisy images are shown in Figures 1(b) and 1(c). It is seen from Figures 1(b) and 1(c) that the smaller 𝐿 is, the more noisy the observed images are. When 𝐿=10, the relative error (ReErr) between the noisy image and the original image is 0.3150. When 𝐿=20, the ReErr between the noisy image and the original image is 0.2234.

(a)

(b)

(c)

Figure 1 

The related images in Test 1. (a) Original image. (b) Noisy image with 𝐿=10. (c) Noisy image with 𝐿=20.

In Figures 2(a)–2(d), we show the restoration results of the two different methods for 𝐿=10 and 𝐿=20. According to the figures, we see that the proposed algorithm can provide better visual quality of restored images than those by the time-marching algorithm. In Figures 3(a) and 3(b), we compare the speed of convergence of the two algorithms for 𝐿=10 and 𝐿=20. In the figures, the 𝑥-axis is the number of iterations and the 𝑦-axis is the the relative error between the restored image and the original image. We see that the speed of convergence of the proposed algorithm is faster and the relative error of the proposed algorithm is smaller. In Table 1, we compare the results of the computational time and the the relative error when the best restored images are obtained. From Table 1, we observe that the computational time required by our method is lesser than that of the time-marching algorithm.

Table 1 

Restoration results with TM and High TV methods for Test 1.

(a)

(b)

(c)

(d)

Figure 2 

Restored results for Test 1. (a) Restored image by the TM algorithm for 𝐿=10. (b) Restored image by the high-order method for 𝐿=10. (c) Restored image by the TM algorithm for 𝐿=20. (d) Restored image by the high-order method for 𝐿=20.

(a)

(b)

Figure 3 

Comparison of the TM algorithm and the proposed method for Test 1. (a) 𝐿=10. (b) 𝐿=20.

In the following, we compare the proposed HighTV method with the one proposed in [18] (AA method), the one proposed in [19] (HNW method) and the one proposed in [20] (MIDAL method). For the AA method, we use the time-marching algorithm with 𝜖=10⁻⁴ to solve the minimization models as proposed in [18]. Similarly, we set the step size as 0.1 in order to obtain a stable iterative procedure. The algorithm is stopped when the maximum number of iterations is reached. In addition, we carry out many experiments with different 𝛼-values in the model (3), the ones with the best results are presented in this work. For the HNW and MIDAL methods, we use the same stopping rule as the HighTV method.

For the MIDAL method, the authors have verified experimentally that setting two parameter be equal yields good results. Therefore, we only select 𝛼 by searching for the value leading to the smallest relative error of such a restored image with respect to the ideal image. In the HNW method, there are two regularization parameters. In [19], the authors fix 𝛼₁=19 in all tests in order to reduce the computational time in search for good regularization parameters. We use the best value of 𝛼₂ such that the relative error of such a restored image with respect to such an ideal image is the smallest. In this way, we have a fair comparison since we compare about the best restorations for four different methods.

The aim of the second test is to give evidence of the effectiveness of employing the propose algorithm over the other three algorithms for multiplicative noise removal problems. Especially, we show that the proposed method can reduce the staircase effect. In this test, we consider the well-known “Lena” image of size 256×256. The original image in Figure 4(a) is distorted by the Gamma noise with 𝐿=3 and 𝐿=13. The noisy images are shown in Figures 4(b) and 4(c). When 𝐿=3, the PSNR and the relative error (ReErr) between the noisy image and the original image are 11.95 dB and 0.5768. When 𝐿=13, the PSNR and the ReErr between the noisy image and the original image are 18.37 dB and 0.2765.

(a)

(b)

(c)

Figure 4 

The related images in Test 2. (a) Original image. (b) Noisy image with 𝐿=3. (c) Noisy image with 𝐿=13.

The restored images by the AA method, the HNW method, the MIDAL method and the proposed method are shown in Figures 5(a)–5(h). From these figures, compared with the AA method, the HNW method and the MIDAL method, the proposed approach yields better results in image restoration since it avoids the staircase effect of the general total variation methods while at same time preserving edges as well as the total variation methods. In Figure 6, we have enlarged some details of the eight restored images. It is clear that the models involving the general total variation transform smooth regions into piecewise constant regions. As it is seen in the zoomed parts, the proposed method outperforms another three methods since the proposed model process smooth regions better than another three models. For the comparison of the performance quantitatively and the computational efficiency, in Table 2, we give their restoration results in SNRs and ReErrs. We observe from Table 2 that both the SNR and ReErr values of the restored images by the proposed method are better than those by the AA method, HNW method and the MIDAL method.

Table 2 

Restoration results with four different methods for Test 2 and Test 3.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 5 

Restored results by different methods for Test 2. (a) Restored image by the AA method for 𝐿=3. (b) Restored image by the HNW method for 𝐿=3. (c) Restored image by the MIDAL method for 𝐿=3. (d) Restored image by our method for 𝐿=3. (e) Restored image by the AA method for 𝐿=13. (f) Restored image by the HNW method for 𝐿=13. (g) Restored image by the MIDAL method for 𝐿=13. (h) Restored image by our method for 𝐿=13.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 6 

A small portion of Lena image is emphasized to show the difference of the four methods for Test 2. (a) Original image. (b) Noisy image with 𝐿=3. (c) Noisy image with 𝐿=13. (d) Restored image by the AA method for 𝐿=3. (e) Restored image by the HNW method for 𝐿=3. (f) Restored image by the MIDAL method for L=3. (g) Restored image by our method for 𝐿=3. (h) Restored image by the AA method for 𝐿=13. (i) Restored image by the HNW method for 𝐿=13. (j) Restored image by the MIDAL method for 𝐿=13. (k) Restored image by our method for 𝐿=13.

In the third test, the “barche” image is considered. The original image in Figure 7(a) is distorted by the Gamma noise with 𝐿=6 and 𝐿=16, see Figures 7(b) and 7(c). When 𝐿=6, the PSNR and the ReErr between the noisy image and the original image are 13.20 dB and 0.4090. When 𝐿=16, the PSNR and the ReErr between the noisy image and the original image are 17.44 dB and 0.2514.

(a)

(b)

(c)

Figure 7 

The related images in Test 3. (a) Original image. (b) Noisy image with 𝐿=6. (c) Noisy image with 𝐿=16.

The information of restored images by the AA method, the HNW method, the MIDAL method and the proposed method are displayed in Figures 8(a)–8(h). From the visual quality of restored images, the proposed regularization method is quite competitive with the other three algorithms. In addition, in Figures 9(a)–9(h) we show the SSIM maps of the restored images. It is not difficult to observe that the SSIM maps by the proposed method are whiter than the AA method, the HNW method and the MIDAL method, which means that the proposed method behaves slightly better. The restoration results of these four different methods in SNRs, ReErrs are also presented in Table 2. It is shown from Table 2 that both the SNR and ReErr values of the restored images by the proposed method are better than those by the AA method, the HNW method and the MIDAL method.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8 

Restored results by different methods for Test 3. (a) Restored image by the AA method for 𝐿=6. (b) Restored image by the HNW method for 𝐿=6. (c) Restored image by the MIDAL method for 𝐿=6. (d) Restored image by our method for 𝐿=6. (e) Restored image by the AA method for 𝐿=16. (f) Restored image by the HNW method for 𝐿=16. (g) Restored image by the MIDAL method for 𝐿=16. (h) Restored image by our method for 𝐿=16.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 9 

Comparison of the SSIM maps for Test 3. (a) The SSIM map by the AA method for 𝐿=6. (b) The SSIM map by the HNW method for 𝐿=6. (c) The SSIM map by the MIDAL method for 𝐿=6. (d) The SSIM map by our method for 𝐿=6. (e) The SSIM map by the AA method for 𝐿=16. (f) The SSIM map by the HNW method for 𝐿=16. (g) The SSIM map by the MIDAL method for 𝐿=16. (h) The SSIM map by our method for 𝐿=16.

6. Conclusions

In this paper, we investigate a fast and efficient way to restore blurred and noisy images with a high-order total variation minimization technique. An alternating minimization algorithm is employed to solve the high-order total variation model. The edges in the restored image can be preserved quite well and the staircase effect is reduced effectively in the proposed algorithm. Our experimental results show that the performance of the proposed method is superior to that of some existing total variation restoration methods.

Acknowledgments

This research is supported by NSFC (10871034, 61170311), Sichuan Province Sci. & Tech. Research Project (2011JY0002, 12ZC1802), Jiangsu Province Research Project (12KJD110001) and Research Grant KQ12013 from HuaiHai Institute of Technology.