A Weberized Total Variation Regularization-Based Image Multiplicative Noise Removal Algorithm

Multiplicative noise removal is of momentous signiﬁcance in coherent imaging systems and various image processing applications. This paper proposes a new nonconvex variational model for multiplicative noise removal under the Weberized total variation (TV) regularization framework. Then, we propose and investigate another surrogate strictly convex objective function for Weberized TV regularization-based multiplicative noise removal model. Finally, we propose and design a novel way of fast alternating optimizing algorithm which contains three subminimizing parts and each of them permits a closed-form solution. Our experimental results show that our algorithm is e ﬀ ective and e ﬃ cient to ﬁlter out multiplicative noise while well preserving the feature details.


Introduction
Image denoising is one of the fundamental problems in image processing and computer vision.Multiplicative noise appears in various image processing applications, for example, in synthetic aperture radar (SAR), ultrasound imaging, single particle emission computed tomography (SPECT), and positron emission tomography (PET) [1].Hence, multiplicative noise removal is of momentous significance in coherent imaging systems and various image processing applications.
The essential idea for image denoising is to filter out noise in an image without losing significant features such as edges and textures.However, one of the challenges in image denoising is its ill-posed nature.To cope with this problem, a large number of approaches have been proposed, most of them under the regularization or the Bayesian frameworks [2,3].These approaches are supported on some form of a priori knowledge or regularization about the original image to be estimated.Some of these methods, including Markov random field priors, wavelet-based priors or regularization [4,5], curvelet-based diffusion [6], and total variation (TV) regularization [7][8][9] are considered the state-of-the-art.Among these approaches, variational functional regularization and partial differential equations-(PDE-) based models have become international popular issues.These models provide a good theoretical foundation to the image denoising task and other inverse problems such as image segmentation and image inpainting, and so forth.And they also simulate the human vision well because of its close relationship to multi-scale analysis theory [10].
1.1.Background.Although other types of noise (e.g., impulse or Poisson noise) have also been studied in the literature of image processing, the term "image denoising" is usually devoted to the problem associated with additive Gaussian noise.The goal of image denoising is to recover ideal u from the observed noisy data u 0 = u + v, where u 0 is an observed image defined on Ω, and Ω ⊂ R 2 denotes an open bounded set with Lipschitz boundary, and v denotes the additive Gaussian noise.To cope with the ill-posed nature of denoising, the regularization techniques are often utilized to solve such a problem.Specifically, the regularization functional-based denoising is given by u = arg min u {E smooth (u) + λE data (u, u 0 )}, (1) where E data (u, u 0 ) is the image fidelity term depending upon the additive noise model, which penalizes the inconsistency between the underestimated recovery image and the acquired noisy image, E smooth (u) is the regularization term which imposes some priori constraints on the original image and EURASIP Journal on Advances in Signal Processing to a great degree determines the quality of the recovery image, and λ is the regularization parameter which controls the trade-off between the image fidelity term E data (u, u 0 ) and the regularization term E smooth (u) [2,3,9,11].The classical model is the minimizing total variational functional [7] E smooth (u) = TV(u) = Ω |∇u|dx with least square data fidelity E data (u, u 0 ) = Ω (u − u 0 ) 2 dx.We refer the readers to [2,3] and references herein for an overview of the subject [9].
In this paper, we focus on the issue of multiplicative noise removal.Specifically, we are interested in the denoising of SAR images.According to [12] and other references, the noise in the observed SAR image is a type of multiplicative noise which is called speckle.And the image formation model is where u 0 is the observed image, u is the original SAR image, and v is the noise which follows a Gamma Law with mean one with its probability density function given by where L is the number of looks (in general, an integer coefficient) and Γ(•) is a Gamma function.Speckle is one of the more complex image noise models.It is signal independent, non-Gaussian, and spatially dependent.Hence, speckle denoising is a very challenging problem compared with additive Gaussian noise.

Prior Works on Variational Approaches.
Multiplicative noise removal methods have been discussed in many reports.Popular methods include the Lee method [13], the Kuan method [14], and various anisotropic diffusion-based methods [15][16][17], which will not be addressed in this paper.We will focus on the variational approaches-based multiplicative noise removal, especially that our researches will emphasis on the TV model-based methods.
To the best of our knowledge, there exist only several variational approaches devoted to multiplicative noise removal problem.The first total variation-based multiplicative noise removal model was presented by Rudin et al. [18], which used a constrained optimization approach with two Lagrange multipliers.The authors derived a denoising model (RLO) as follows: and designed a gradient projection-based algorithm.Following the maximum a posteriori (MAP) estimator for multiplicative Gamma noise, Aubert and Aujol [19] introduced a nonconvex model (AA): where E data (u, u 0 ) = Ω (log u + u 0 /u)dx is the data fitting term and E smooth (u) = TV(u) = Ω |∇u|dx is the total variation (TV) regularization term.They also used gradient descend method to solve AA model.Recently, Shi and Osher [20] adopted a noisy observation log u 0 = log u + log v and the data term of the AA model then derived the TV minimization model for multiplicative noise removal problems.They applied a corresponding relaxed inverse scale space flow as denoising technique.Huang et al. [21] proposed a strictly convex model by letting w = log u and replacing TV(w) by TV(u) based on the prior work in [19], and they proposed a simpler alternating minimization algorithm by adding a quadratic term to the model.A variational model involving curvelet coefficients for cleaning multiplicative Gamma noise was considered in [22].Steidl and Teuber [1] introduced a variational restoration model consisting of the I-divergence as data fitting term and the total variation seminorm as regularizer.They applied Douglas-Rachford splitting techniques, respectively, alternating split Bregman methods to denoising.

Contributions.
The main aim of this paper is to propose and investigate another nonconvex variational model for multiplicative noise removal which inspiring from the Weberized TV regularization method [23,24].We also incorporate another way of surrogate functional model to recover image edges.We develop an alternating minimization algorithm to find the minimizer of such an objective function efficiently.Our experimental results show that our proposed method has good performance for multiplicative noise removal.The outline of this paper is as follows.In Section 2, we introduce our Weberized total variation-based multiplicative noise removal model.In Section 3, we propose a surrogate functional model which can be viewed as an approximating way of the new nonconvex model.In Section 4, we develop a novel fast algorithm for minimizing such surrogate functional model.In Section 5, we show experimental results to demonstrate the quality of the denoised images and the efficiency of our proposed algorithm.Finally, concluding remarks are given in Section 6.

Weberized TV Regularization Model for Multiplicative Noise Removal
All images are eventually perceived and interpreted by the Human Visual System (HVS).As a result, many researchers have found that human vision psychology and psychophysics play an important role in image processing.The classical example is the using of the Just Noticeable Difference Model (JND) in image coding and watermarking techniques [25,26].In these fields, the JND model is used to control the visual perceptual distortion during the coding procedure and watermark embedding.Weber's law was first described in 1834 by German physiologist Weber [27].The law reveals the universal influence of the background stimulus u on human's sensitivity to the intensity increment |∇u|, or so-called JND, in the perception of both sound and light: According to Weber's law, when the mean intensity of the background is increasing with a higher value, the intensity increment |∇u| also has higher value.In literature [23], the author gave a complete analysis and report for the Gaussian denoising problem: He proposed a nonconvex variational model: The essential idea of the above model is the use of Weber's law.For the optimizing model ( 8), the author had given complete existence and uniqueness theorem in the following natural admissible space S(Ω) defined as and presented a fast gradient descent algorithm.Inspired from the Weberized TV regularization method, we propose a nonconvex Weberized TV regularization-based multiplicative noise removal model: where the first term E w (u) := TV(log u) = Ω (|∇u|/u)dx is the well-known Weberized TV regularization term, while the second one is the nonconvex data fidelity term.Comparing the model in (10) with the model in (5), it is very interesting to see that the model in ( 10) is essentially similar to Weberized TV Regularization for Gaussian noise removal model (8), and the difference is that the data fidelity term is adjusted by AA model's data fidelity term.Furthermore, the proposed model ( 10) can be understood and interpreted from the statistical perspective using Bayesian formulation, and this is straightforward application of MAP theory as illustrated in [19,21].
In addition, one may wonder whether the minimizing problem of (10) has the unique solution or not.It is very important to point that existence and uniqueness of this problem can also be valid and can be proved.The proof can be made very similar to the discussion for the model ( 8), which can be found in [13].

The Proposed Surrogate Functional Model and Mathematical Analysis
3.1.The Proposed Surrogate Functional Model.In this subsection, we will propose another extension strictly convex objective function which can be viewed as a surrogate functional model for minimizing Weberized TV regularization problems (10).To do this, we first give the formal equilibrium Euler-Lagrange equation of problem (10).For the sake of convenience, let us denote that Φ(u) = 1/u.
Proof.Note that in the first integral of the second line, by the standard computation of the operator div(Φ(u)∇u/|∇u|), it is easy to prove that div Φ(u) Then we use the standard computation of Calculus of Variation E → E + δE: where ds denotes the arc-length element of the boundary.This completes the proof.
Since u > 0, Φ(u) = 1/u, then the Euler-Lagrange equation ( 12) can be rewritten equivalently as EURASIP Journal on Advances in Signal Processing Suppose that there exists a function M(u) which satisfies M (u) = (u − u 0 )/u, then it is easy to find the function M(u) = u − u 0 log u, and then ( 15) is also the corresponding Euler-Lagrange equation of another objective functional or a new reference energy E R (u) which is defined as follows: Note that optimization problems given in (10) and ( 16) have the same equilibrium Euler-Lagrange equation (15).In Section 3.2, we can prove that the solution of optimizing model ( 16) admits a unique solution.In other words, on one hand, the unique solution of optimization problems given in ( 10) is also a solution of optimizing model (16).
On the other hand, optimizing model ( 16) admits a unique solution which satisfies (15).Taking into consideration two facts mentioned above, it is clear to see that the equivalence of optimization problems given in ( 10) and ( 16) can be guaranteed.Hence, model ( 16) can be viewed as a surrogate functional for minimizing the Weberized TV regularization problem (10).

Existence and Uniqueness Results
. As discussed in Section 3.1, to show the equivalence of optimization problems given in ( 10) and ( 16), we need to prove the existence and uniqueness of solution in optimizing model (16).To do this, we briefly recall some notations and basic preliminaries of BV(Ω) (see [3,24]).Let C p 0 (Ω) denote the space of real-valued functions, p is continuously differentiable with compact support, let L p (Ω) denote the space of Lebesgue measurable functions u such that Ω |u| p dx < ∞ and L ∞ (Ω) the space of Lebesgue measurable functions u such that there exists a constant c with |u(x)| ≤ c a.e.x ∈ Ω. Definition 1.Let BV(Ω) be a space of function of u ∈ L 1 (Ω) such that the following quantity: We summarize below the lower semicontinuity and compactness properties of BV(Ω) [3,28] that we will use in the proof.
n=1 is a sequence in BV(Ω) satisfying sup n u n BV(Ω) < ∞, then there exists a subsequence {u nj } ∞ j=1 and a function u ∈ BV(Ω) such that u nj → u in L 1 (Ω) as j → ∞.
Using the lower semicontinuity and compactness of BV(Ω), we can conclude the following existence and uniqueness theorem.
Theorem 1.Let u 0 ∈ L ∞ (Ω) be a positive, bounded function with inf Ω u 0 > 0, then the minimizing problem of energy functional in (16) Proof.See the appendix.

Proposed Fast Algorithm
4.1.Motivation.In this subsection, we will develop a fast multiplicative noise removal algorithm for the optimizing energy functional of ( 16).Our algorithm is designed using the well-known variable-splitting and penalty techniques in optimization [20][21][22].An auxiliary variable w is firstly introduced to transfer u out of the nondifferentiable term TV(u), and the difference between w and u is penalized, yielding the following approximation model of model ( 16): with a sufficiently large penalty parameter α 2 > 0. It is well known that the solution of ( 19) converges to that of ( 16) as α 2 → ∞.While either one of the two variables u and w is fixed, problem (19) can be solved by alternative iterative subproblem minimizing: (ii) w = arg min The main advantage of this procedure is that the proposed model takes advantage of using the total variation minimization scheme to remove the multiplicative noise.It can be solved by many TV denoising fast methods such as the Chambolle projection algorithm [29], primal-dual method, and the lagged diffusivity fixed point method [30][31][32][33].In this paper, we will use another fast algorithm proposed in [34,35] to deal with TV denoising problem (21).Similar to the previous discussion, we introduce another vector variable v = (v 1 , v 2 ) to approximate ∇w, and then subproblem (21) is approximated by where α 1 > 0 is the penalty parameter, and the solution of ( 22) converges to that of (21) as α 1 → ∞.

Discrete Scheme.
Let us consider the discrete scheme of the problem.Let matrix u 0 ∈ R N×N represent a two dimensional gray-scale digital image, where each component (u 0 ) i, j is the intensity value of pixel (i, j) for i, j ∈ {1, 2, . . ., N}.
Similarly, we let matrix u represent the unknown image to be restored.Let us define the forward finite difference operator with periodic boundary adjustments for i = N and j = N, and D is a block-circulant-circulant-block (BCCB) matrix.Under the above definitions, the corresponding discrete form of ( 20) and ( 22) is (i) u-subproblem: (ii) w-subproblem: 4.2.1.u-Subproblem.Since this is strictly convex problem, it is equivalent to solve or (let γ = α 3 /(2α 2 )) find the meaningful positive root for this equation, that is,

w-Subproblem.
It can be solved by alternately minimizing the objective function with respect to v while fixing w, and vice versa: It is not difficult to verify that the v-subproblem permits a closed-form solution [34]: For a fixed v, let μ = α 1 /α 2 , then w-subproblem (29) is quadratic in w and the minimizer w is given by the normal equation [34,35] Under the periodic boundary condition for w, D 1 , and D 1 are block circulant.So D T 1 D 1 and D T 2 D 2 are all block circulant.Therefore, the Hessian matrix on the left-hand side of ( 31) can be diagonalized by two-dimensional discrete Fourier transform F. Using the convolution theorem of Fourier transforms, we can write where the symbol " * " denotes complex conjugacy, "•" denotes component-wise multiplication, and the division is component-wise as well.With a slight abuse of notation, we have used F(D 1 ) for the Fourier transform of the function represented by D 1 in the convolution D 1 u (and similarly for D 2 ).Since all quantities but v 1 and v 2 are constant, computing w from (32) involves two FFTs and one inverse FFT, once the constant quantities are computed.

Parameters Choice.
There are three parameters α 1 , α 2 , and α 3 or equivalently α 1 , γ = α 3 /(2α 2 ), and μ = α 1 /α 2 involved in the iterative procedure.Firstly, from the connection between model ( 16) and model (19), how well the solution of (19) approximates that of (16) or its constrained equivalent depends on the magnitude of α 1 , α 2 , which determines the amount of penalty applied to the discrepancy between u and w and also between ∇w and v in the squared L2-distance.Hence, the magnitude of α 1 , α 2 must have larger enough value.
Secondly, α 3 is the regularization parameter which controls the trade-off between the image fidelity term and the regularization term.We dynamically compute the value of α 3 according to the variance of the recovered noise which matches that of our prior knowledge [7].The Gamma distributed noise has the mean and variance as follows: Notice that the Euler-Lagrange equation of the optimizing problem of ( 16) is (15).The solution procedure uses a parabolic equation with time as an evolution parameter.This means that we solve for t > 0. We merely multiply the first equation of ( 34) by ((u 0 − u)/u) and integrate by parts over Ω.If steady state has been reached, the left side of the first equation of (34) vanishes, then we have 4.4.Proposed Fast Algorithm.In summary, according to the iterative process analysis and parameters choice which were discussed in Sections 4.1 and 4.2, respectively, we can propose a fast algorithm framework as follows.Two loops of iterations are contained in this algorithm framework.In the outer loop, we use the continued method by increasing the parameters α 1 and α 2 to achieve the good convergence.At the same time, in the inner loop, we update the parameter α 3 in order to match the variance of the recovered noise.
According to the optimality conditions of ( 20), (28), and ( 29), the stopping criterion of the proposed algorithm is proposed in the following.Let The inner loop is terminated once where Res measures the total residual and ζ > 0 is a prescribed tolerance.
The complete resulting algorithm is summarized in Algorithm 1.

Convergence Analysis.
The convergence of the quadratic penalty method as the penalty parameter goes to infinity is well known (e.g., see [36,Theorem 17.1]).That is, as α 1 , α 2 → ∞ the solution of (19) converges to that of (16), and the solution of ( 22) converges to that of (21), respectively.In this subsection, we present the convergence results of Algorithm 1 for fixed α 1 and α 2 without proofs since these results are rather straightforward application of [34,Theorem 3.4].For the sake of completeness, we first present some necessary definitions and then give the main convergence results.
Algorithm 1 where 0 • (0/0) = 0 is followed.An then we define an Ndimensional vector shrinkage operator S(u; v) : In addition, we define a linear operator h : Using the definition of S and h, we can rewrite the three iterative steps in Algorithm 1 into then the convergence results of Algorithm 1 can be proved as discussed in [34].Firstly, as proved by Proposition 3.1 in [34], the nonexpansiveness of the shrinkage operator S can be ensured.Then, it is easy to check that M is nonsingularity, thus another symmetric positive definite matrix T = DM −1 D T is well defined and the spectral radius ρ(T) < 1.Thus the operator h is also nonexpansive.Secondly, the objective function in ( 25) is convex, bounded below, and coercive.Hence, for any fixed α 1 > 0, the sequence {(v (m) , w (m) )} generated by (40), (41), and (42) from any star point w (0) converges to a solution of (25) (see [34,Theorem 3.4]).Furthermore, the convergence of {u (m) } to some u * follows directly from (40), hence we can conclude that for fixed α 1 and α 2 , the sequence {(u (m) , w (m) )} generated by Algorithm 1 from any start point (u (0) , w (0) ) converges to a fixed point (u * , w * ), which is the solution of ( 24) and ( 25).

Some Complexity Notes.
It is clear that the complexity of the proposed algorithm mainly includes three parts.
The calculation in ( 27) and ( 30) has a linear-time complexity of order O(N 2 ) for an N-by-N image.Hence, the usubproblem and v-subproblem can be solved quickly.The solution of the w-subproblem (32) requires three fast Fourier transforms and one inverse transform and a total complexity in the order of O(N 2 log(N 2 )) = O(N 2 log(N)).

Numerical Results and Performance Analysis
In this section, we demonstrate the effectiveness of our proposed Algorithm 1 in image denoising.The numerical results are compared with those obtained by the "HMW" method proposed by Huang et al. [21], "AA" method proposed by Aubert and Aujol [19], and the "RLO" method proposed by Rudin et al.
where u, u * , and u 0 are the original, the restored, and the observed images, respectively.The solution of model ( 4) by the "RLO" method is obtained by using the following gradient projection iterative scheme [18]: Here, the two Lagrange multipliers λ and μ are dynamically updated to satisfy the constraints (as explained in [18]).
The solution of the "AA" model ( 5) is obtained by using the following explicit iterative scheme [19]: Here, the regularization parameter λ is dynamically updated (as explained in [7]).In the "RLO" and "AA" methods, ε is set to be 10 −4 , and Δt is set to be a small positive number to ensure the convergence of the iterative scheme.The solution of the "HMW" model [21] min is obtained by using the following alternating minimization algorithm: The corresponding nonlinear Euler-Lagrange equation of subproblem (47) was solved using the Newton method.The Chambolle projection algorithm was employed in the denoising subproblem (47), and then the restored image was computed by exp(w).The stopping criterion of the "HMW" method, the "AA" method, and the "RLO" method is that the PSNR value of the restored image reaches its maximum.
Our numerical experiments use several natural and manmade images such as "Lena" (Figure 1  show the restored "Lena," synthetic image and SAR images, respectively.In these experiments, the initial guess is set to be u (0) = u 0 .It is clear that the restoration results by the proposed method are visually much better than those by the "HMW" method, the "AA" method and the "RLO" method, especially when the noise variance is large, that is, when L is small.The effectiveness of our

2 ,ISNR = 10 log 10 u 0 − u 2 2 u
[18].In the tests, each pixel of an original image is degraded by a noise which follows a Gamma distribution with density function in (3) and v being specified to have mean 1 and standard deviation 1/ √ L. The noise level is controlled by the value of L in the experiments.To compare the performance of the algorithms mentioned above, we compute the quality of restored images by the peak signal-to-noise ratio (PSNR), the improved SNR (ISNR), and the relative error (ReErr) of the restored image defined by PSNR = 10 log 10 MN max {u} 2 u * − u 2

Figure 9 :
Figure 9: The 100th line of the original, noisy, and restored images for denoising Figure 2(b).(a) The noisy image; (b) the slice restored by the proposed method; (c) the slice restored by the "HMW" method; (d) the slice restored by the "AA" method; (e) the slice restored by the "RLO" method.Here the blue line is the original image, and the red line is the restored image.

Figure 10 :
Figure 10: The 128th line of the original, noisy, and restored images for denoising Figure 3(a).(a) The noisy image; (b) the slice restored by the proposed method; (c) the slice restored by the "HMW" method; (d) the slice restored by the "AA" method; (e) the slice restored by the "RLO" method.Here the blue line is the original image, and the red line is the restored image.

Figure 3 (Figure 11 :
Figure 11: The 128th line of the original, noisy, and restored images for denoising Figure 3(b).(a) The noisy image; (b) the slice restored by the proposed method; (c) the slice restored by the "HMW" method; (d) the slice restored by the "AA" method; (e) the slice restored by the "RLO" method.Here the blue line is the original image, and the red line is the restored image.