Multiplicative noise removal using primal–dual and reweighted alternating minimization

Multiplicative noise removal is an important research topic in image processing field. An algorithm using reweighted alternating minimization to remove this kind of noise is proposed in our preliminary work. While achieving good results, a small parameter is needed to avoid the denominator vanishing. We find that the parameter has important influence on numerical results and has to be chosen carefully. In this paper a primal–dual algorithm is designed without the artificial parameter. Numerical experiments show that the new algorithm can get a good visual quality, overcome staircase effects and preserve the edges, while maintaining high signal-to-noise ratio.

improvement of the iteratively reweighted algorithm is introduced without the artificial parameter, The rest of the paper is organized as follows: "Iteratively reweighted model with TV" section briefly introduces the IR-TV model as well as the classical iteratively reweighted algorithm. "Solution to the model" section presents the new algorithm of IR-TV model, which is based on the primal-dual optimization and without the artificial parameter. In "Numerical experiment" section the effectiveness of the proposed algorithm is verified through numerical experiments. Finally we conclude in "Conclusion" section.

Iteratively reweighted model with TV
Suppose the degraded image f (x) = u(x)v(x), x ∈ �, where the original image u(x) is a real piecewise smooth function defined on a bounded domain ⊂ R 2 , and the multiplicative noise v(x) is assumed to obey Gamma distribution with mean 1 In Eq. (1), Γ (·) is a Gamma function with variance 1/L. Iteratively reweighted l 1 regularization minimization problem attempts to find a local minimum of concave penalty functions that more closely resembles the l 0 regularization problem (Simon and Lai 2009;Candes et al. 2008). In our previous work Wang et al. (2012), we put forward an iteratively reweighted model where z(x) = log u(x) and φ(z) = |∇z| ,regularizer parameter μ is a constant connected with the intensity of noise, g(x) is a nonnegative weight function which controls the strength of smoothing. According to the classical iteratively reweighted algorithm, we choose where n is the number of outer iteration. It is obvious that the larger |∇z|, the weaker smoothing strength is, thus the noise is removed while the edges are preserved.
The classical algorithm to Eq. (2) attempts to find a local minimum of a concave function, whereas in each iteration the algorithm simply requires to solve a convex optimization problem, In order to prevent the zero denominator, Eq. (3) usually be revised to The parameter ε (n) provides the stability for iterations. The choice of ε (n) has a significant effect on the result of the denoising. Therefore it needs to carefully adjusted. It will lead to poor denoising results with a inappropriate ε (n) (Wang et al. 2012;Simon and Lai 2009).
In next section, we propose a novel algorithm to solve Eq. (2). First the splitting method is used to transform the original equation into two corresponding equations. Then the primal-dual algorithm and the Euler-Lagrange method are applied to solve these two subproblems respectively. As Huang et al. (2009) has mentioned, let us consider the splitting form of Eq. (2) where w is an auxiliary function, The parameter γ measures the amount of regularization to a denoising image, which is large enough to make w be close to z. In our experiment, γ = 19 is chosen. The main advantage of the proposed method is that the TV norm can be used in the noise removal process in an efficient manner. And Eq. (5) can be splitted into two equations This is an alternating minimization algorithm. The first step of the method is to apply a weighted TV denoising scheme to the image generated by the previous multiplicative noise removal step. The second step of the method is to solve a part of the optimization problem.

Solution to the model
In this paper, a primal-dual algorithm (Bertsekas et al. 2006;Bertsekas 2011) is applied to iteratively reweighted model 6a. Convex close set K is defined by, where {·} denotes convex close set of {·}.
Let X, Y be two finite dimention real vector spaces, the corresponding norm defined as �·� = �·, ·� 1/2 , where �·, ·� is the inner product. Gradient operator ∇ : X → Y is continuous linear operator, the corresponding norm defined as We introduce a divergence operator div : X → Y , the adjoint of divergence operator is defined by ∇ * = −div. Then we introduce dual variable p = (p 1 , p 2 ), which divergence is divp = ∂p 1 ∂x 1 + ∂p 2 ∂x 2 , we have The regularizator of 6a is and 6a can be transformed into for every w ∈ X and > 0, J ( w) = J (w) holds, so J is one-homogeneous. By the Legendre-Fenchel transform, we can obtain with J * (v) is the "characteristic function" of a closed convex set K: Since J * * = J, we recover The Euler equation for (7) is where ∂J is the "sub-differential" of J. Writing this as we get that q = 2γ z (n−1) − w /µ is the minimizer of q − 2γ z (n−1) /µ 2 + 2γ µ J * (q). Since J * is given by (3), the solution of problem (6) is simply given by Therefore the problem to compute w (n) become a problem to compute the nonlinear projection q = π µK / 2γ z (n−1) . Consider the following problem: Following the standard arguments in convex analysis (Chambolle 2004;Chambolle and Pock 2011), the Karush-Kuhn-Tucker conditions yield the existence of a Lagrange multiplier α i,j (x) ≥ 0, such that constraint problem (9) become to, Notice constraint problem |p| ≤ g in Eq. (10). For any x, α(x) ≥ 0, if |p| 2 < g 2 ,then α(x) = 0; If |p| 2 = g 2 , we see that in any case Then Substituting (11) into (10) gives, We thus propose the following semi-implicit gradient descent (or fixed point) algorithm. We choose τ > 0, let p 0 = 0 and for any n ≥ 0, The denominator of Eq. (13) is greater than zero, which avoids the appearance of the rectified parameters, and of course does not need to be adjusted. The method can be seen a new method to solve nonconvex problem. We need to calculate the boundary of the norm div .
In the following, we will prove the convergence of µ 2γ divp m . Let s = lim m→∞ divp m − 2γ z (n−1) µ , and p be the limit of a converging subsequence p m k of {p m }. Letting p ′ be the limit of p m k +1 , we have and repeating the previous calculations we see It holds η i,j = p ′ i,j −p i,j δt = 0, for any i, j,i.e., p ′ =p So we can deduce which is the Euler equation for a solution of (8). One can deduce that p solves (8) and that µ/2γ divp is the projection of π µK /2γ z (n−1) . Since this projection is unique, we deduce that all the sequence µ/2γ divp n converges to π µK / 2γ z (n−1) as κ 2 ≤ 8. 6b is equivalence to solve the nonlinear system

Numerical experiment
We compare our algorithm on eliminating staircase effect and preserving the detail to SO model, HNW model and classic iteratively reweighted total variation (CWTV). Signal to Noise Ratio (SNR) of the denoising image to the corresponding true image is defined as where X is the denoised image and X is the true image. We stop algorithm while attaining maximum SNR. The test images are, "Shape1", "Shape2", "Barbara", "Lena256", "Cameraman", "Phantom". The multiplicative noise with standard variance (NSV) of 1/30 and 1/10 are considered in our experiments. Table 1 shows the effect of artificial parameter ε n to denoising results of classic iteratively reweighted isotropous total variation method. Table 2 is the comparison of denoising results on SNR. From Table 1, we can explicitly see that suitable artificial parameter ε n can obtain better denoising results than some other models (such as SO model, HNW model), while unsuitable artificial parameter ε n obtain lower SNR than other models. New algorithm can obtain the highest SNR than SO model, HNW model and classic iteratively reweighted method. Moreover the new algorithm is not affected by this parameter.

Experiment 1: Comparison on eliminating staircase effect
"Shape1" is used as a test image in this experiment, the multiplicative noise intensity is standard variance 1/10. In our algorithm, µ = 0.013 and the number of inner iteration is set 30, the denoising SNR result can achieve 12.3856 dB. Figure 1 is the denoising results. Comparing Fig. 1c-f, we can see, staircase effect is restrained in the alternative splitting minimizating algorithm (HNW model and our algorithm), and the transition of smooth region in the new model has a good visual effect. Moreover, we can clearly find new model can preserve edge and detail better than SO model, HNW model. The edge and details of the restored images are preserved because of the action of the weighted function. In Fig. 1 short widthways lines in our methods can be restored more number than SO model and HNW model.  Experiment 2: Detail preserving "Shape 2" and "Lena256" images are contaminated by multiplicative noise with standard variance 1/10. Figures 2 and 3 are the denoising results. In our algorithm to "Shape 2", µ = 0.015 and the number of inner iteration is set 30, the denoising SNR result can achieve 16.0540 dB. We can see the denoising results is better than the SO model and HNW model. In our algorithm to "Lena256", µ = 0.0025 and the number of inner iteration is same as the experiment 1, and the denoising SNR result can achieve 13.9022 dB.
The preserved detail of our algorithm is better than the SO model and HNW model, especially the feather on the cap. On the edge of the image, the derivative of image edges is bigger, then weight function value becomes little and the degree of polishing is weakened to the edges. thus the edges are preserved; On the other hand, The derivative of the smooth regions is much small, weighted function is large, which strengthen the smoothing to relatively smooth regions, thus the noise is removed. Compare to Figs. 2 and 3c-f, it is obvious that the denoising results of proposed algorithm can keep details better.

Conclusion
We study a new algorithm on iteratively reweighted to remove multiplicative noise model. An alternating minimization method is employed to solve the proposed model. And a Chambolle projection algorithm to iteratively reweighted model is proposed. Our experimental results have shown that the quality of images restored by the proposed method is quite good, especially on preserving the detail and restraining the staircase effect. Moreover the proposed algorithm provides an approach to solve the non-convex problem.