# Variable Splitting Based Method for Image Restoration with Impulse Plus Gaussian Noise

Authors: Tingting Wu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Image denoising is a fundamental problem in realm of image processing. A large amount of literature is dedicated to restoring an image corrupted by a certain type of noise. However, little literature is concentrated on the scenario of mixed noise removal. In this paper, based on the model of two-phase method for image denoising proposed by Cai et al. (2008) and the idea of variable splitting, we are capable of decomposing the image denoising problem into subproblems with closed form. Numerical results illustrate the validity and robustness of the proposed algorithms, especially for restoring the images contaminated by impulse plus Gaussian noise.


---

## Body

## 1. Introduction

When recording a real-world scene by camera, we desire to acquire an ideal image which is a faithful representation of the scene. However, most observed images are more or less involved in blurry and noise. Hence, deblurring and denoising to the observed image are fundamental aspects in image processing. Letx∈Rn1×n2 be an ideal image scaled in [0,1]. Hereafter, we represent an n1×n2 image x as x∈Rn (n=n1n2) by stacking its columns as a vector. H∈Rn×n is matrix representation of spatially invariant point spread function (PSF) which characterizes the blurry effects on imaging system and n∈Rn is additive noise with some statistical distributions (e.g., Gaussian noise follows normal distribution, impulse noise follows binomial distribution or uniform distribution, and uniform white noise follows uniform distribution). Accordingly, the observed image x0∈Rn involved in the blurring matrix H and the additive noise n can be boiled down to(1)x0=Hx+n.The main task of image restoration is to recover the ideal imagex from the observed image x0. There are two difficulties in finishing this task. The first one is that the singular values of H decay asymptotically to zero which results in a large condition number, and as a consequence (1) is an ill-posed inverse problem in mathematical sense. The second difficulty is that the inevitably additive noise n in the observed image makes (1) difficult to handle. These two difficulties make the linear least squares solvers generally fail in deriving ideal image x and regularization strategy is therefore introduced to improve the numerical performances, for example, Tikhonov regularization [1], Mumford-Shah regularization [2], and total variation (TV) regularization [3]. The main superiority of TV regularization is that it can regularize images without blurring their edges. The common models utilizing TV regularization for image restoration are the constrained model (2)min∇x1s.t.Hx-x0p≤σ,p=1,2,∞,and it is equivalent to (or rather equivalent under a proper τ, see [4]) the unconstrained model(3)minxτ∇x1+Hx-x0pp,p=1,2,∞,where τ, σ are parameters depending on noise variance; ∇ is first-order derivative operator and ∇·1 is the TV-norm (details in Section 2); ·p denotes the p-norm in Euclidean space (if p=∞, ·pp herein should be comprehended as ·∞). Commonly, p is hinged on the statistical distribution of the additive noise n; for example, p=1 for impulse noise removal, p=2 for Gaussian noise removal, and p=∞ for uniform white noise removal. It is therefore pivotal to figure out statistical information of n before opting for mathematical model.Recently, Huang et al. [5] proposed a modified TV regularization model for Gaussian noise removal(4)minx,yτ∇x1+μ2x-y22+12Hy-x022,where τ, μ are positive parameters. Intuitively, a new fitting term x-y2 is added to the unconstrained model (3). Essentially, by smoothing of the nonsmooth function ∇x1 via the Moreau approximation, the objective function of the following problem (5)minyτ∇y1+Hy-x022can be replaced by(6)minx,yτ∇x1+μ2x-y22+12Hy-x022,which is identical with (4). They employed the alternating minimization (AM) algorithm to handle model (4) in variables x and y, respectively. The numerical results therein indicated that the model was attractive. Lately, they extended model (4) to the mixed noise (impulse plus Gaussian noise) removal in [6]. The extended model reads(7)minx,yτ∇x1+μ2x-y22+1pPAHy-x0pp,p=1,2,where A⊂Ω≔{1,2,…,n} contains the noise-free pixels (pixels are not contaminated by impulse noise); PA is the characteristic function of set A. Specifically, for any x∈Rn, j∈{1,2,…,n}, PA(x) can be defined as (8)PAxj=0,ifj∈Ω∖A,1,ifj∈A,which is convex and nonsmooth. Then, by using the incomplete data set in A, TV minimization method is constructed for image restoration and suitable values can be filled in the detected noisy pixel locations.As interposition, we expand in this paragraph about mixed noise removal. Gaussian noise [3] generally results from the analog-to-digital conversion of measured voltages and it corrupts an image by exerting subtle perturbations on the gray scale of all pixels. Impulse noise [7] typically arises from the malfunctioning arrays in camera sensors or faulty memory locations in hardware and it affects an image by seriously devastating the gray scales of some pixels (but keeping the other pixels noise-free). Salt-and-pepper noise [8] and random-valued noise [9] are two categories of impulse noise. Let dmin,dmax be the dynamic interval of image x. For given noise level (also the so-called noise intensity) s∈(0,1), salt-and-pepper noise and random-valued noise contaminate image x, respectively, by the following:(i)
Salt-and-pepper noise (see [10]) (9)Nimpxi=dmin,with  probabilitys2,dmax,with  probabilitys2,xi,with  probability1-s.(ii)
Random-valued noise (see [11]) (10)Nimpxi=di,with  probabilitys,xi,with  probability1-s,whereNimp:Rn→Rn represents the procedure of contaminating image x with impulse noise and di satisfies uniform distribution in interval [dmin,dmax].Figure 1 shows distinctions of image contaminated by salt-and-pepper noise and random-valued noise. Visually, the salt-and-pepper noise corrupted image possesses black-white speckles while the random-valued noise corrupted image has speckles with various gray scales.Figure 1
Impulse noise. (a)50×50 white-black ideal image. (b) Contaminated by salt-and-pepper noise with s=10%. (c) Contaminated by random-valued noise with s=10%.
(a) (b) (c)The observed imagex0 is contaminated by mixed noise via(11)x0=NimpHx+ng, where ng is additive Gaussian noise. It is tricky to retrieve the ideal image x in (11) from its observation x0 because the statistical distribution of noise in x0 is superimposed.Cai et al. [12] proposed a two-phase method for mixed noise removal. They first detect the noise-free set A by median-type filter [10, 13] and then execute the procedure of image deblurring on set A. The model therein is(12)minxτ∫Ω∖Γ∇x2dx+μ∫Γdσ+∫AHx-x0pdx,p=1,2,where Γ⊂Ω is the set of image edges and the option for p is determined experimentally. Model (12) is attractive to render restoration with fine structures in image, such as neat edges and textures. However, the nonconvexity of objective function in model (12) makes it difficult to be minimized. A convex approximation of model (12) is considered therein. However, the computational effort on minimizing that convex approximation is time-consuming. Some literature is devoted to the case of mixed noise removal recently because of its realistic applications, for example, [6, 12, 14–16].As aforementioned, Huang et al. [6] proposed a fast method for mixed noise removal by solving model (7). The two-phase pattern is exploited as in [12]. Concretely, one has the following.(i)
Phase  1: Noise Detection. Median-type filters are favorable for impulse noise detection. Adaptive median filter (AMF) [10] is adopted for salt-and-pepper noise detection and adaptive center-weighted median filter (ACWMF) [13] is chosen for random-valued noise detection. Assuming that x0 is the observed image and y0 is the restored image by median-type filter, the set of pixels corrupted by impulse noise, denoted by N, is detected via
(a)
salt-and-pepper noise:N=i∣yi0≠xi0,xi0∈{dmin,dmax},i=1,2,…,n.(b)
random-valued noise:N=i∣yi0≠xi0,i=1,2,…,n.(ii)
Phase 2: Minimization on model (7). Model (7) is solved by AM algorithm. Specifically, the x-related subproblem is solved by Chambolle’s projection algorithm [17] and the y-related subproblem is solved by preconditioned conjugate gradient (PCG) method. The numerical results in [6] illustrated that their algorithm is efficient and promising. It is about 8–10 times (resp., 15–20 times) faster than the algorithm proposed in [12] for impulse noise (resp., mixed noise removal). However, the properties of model (7) can be further exploited by variable splitting strategy and the resulting subproblems can be tractably tackled with closed-form solutions. This is also the main goal of our paper. Since model (4) is an exceptional case of model (7) with A=Ω and p=2, we concentrate on handling model (7) in this paper by variable splitting strategy.The contributions of this paper are as follows. Firstly, we extend models (4) and (7) to deal with different kinds of image denoising, that is, Gaussian noise, salt-and-pepper noise plus Gaussian noise, and random-valued noise plus Gaussian noise in various ways of blur kernels. Secondly, the variable splitting strategies with efficient closed-form solution are investigated. Concretely, we propose the recursions of optimizations (4) and (7) by the classical alternating direction method (ADM). Furthermore, we contemplate coping with (7) with three-block ADM scheme. To guarantee the convergence of algorithm, the ADM with Gaussian-back-substitution is considered for handling model (7). Overall, compared with some standard methods, our proposed models and algorithms exhibit the attractive performance of restoring the images contaminated by impulse plus Gaussian noise.The rest of the paper is organized as follows. In Section2, we summarize some basic notations and properties that will be useful in the subsequent sections. In Sections 3 and 4, algorithms based on variable splitting are presented and followed by some numerical experiments in Section 5. The numerical comparisons with the algorithm for Gaussian noise removal in [5] and the algorithm for mixed noise removal in [6] are implemented therein. Finally, we end the paper with some conclusions in Section 6.


## 2. Preliminaries

In this section, we summarize some basic concepts and properties that will be used in subsequent analysis. Letx∈Rn be a vector and p≥1 be an integer. By xp≔∑i=1nxip1/p, we denote the p-norm of x and x denotes the Euclidean norm for brevity. The inner product in Euclidean space is denoted ·,·. For any positive definite matrix M∈Rn×n, the matrix norm of x is defined by xM≔x,Mx. We denote by I the identity operator. For any u=(u1,u2)∈Rn×Rn, we define a vector u∈Rn with the ith component ui=(u1)i2+(u2)i2. With the definition of p-norm, it follows that u1=∑i=1nui is the TV-norm of u. Let f:Rn→[-∞,+∞] be a function. The domain (resp., epigraph) of f is denoted by dom⁡f≔x∈Rn∣fx<+∞ (resp., epi⁡f≔(x,y)∈Rn×R∣f(x)≤y). If dom⁡f is nonempty, the function f is said to be proper. The function f is called lower semicontinuous (l.s.c.) if epi⁡f is a closed set in Rn×R, and f is a convex function if epi⁡f is convex set in Rn×R. The subdifferential of f, denoted by ∂f:Rn→2Rn, is defined (13)∂fx≔ξ∈Rn∣fy≥fx+ξ,y-x,∀y∈Rn.For any β>0, the proximity function of f(u)=βu1 is the well-known soft-thresholding (also the so-called shrinkage). With special denotation S:Rn×Rn→Rn×Rn, it is defined as(14)Sβu=u-min⁡u,βuu,where all operations are calculated in componentwise and (u/u)i is taken as zero if ui=0. Let X⊆Rn be nonempty set. The indicator function of a nonempty set X⊂Rn is denoted by (15)ιXx≔0,ifx∈X,+∞,ifx∉X.The proximity operator of f, denoted by proxf:Rn→Rn, is defined as [18, 19] (16)proxfy≔argminxfx+12x-y2,∀y∈Rn.Specially, the projection operator on a nonempty closed convex set V⊆Rn, denoted by ΠV, is the proximity operator of ιV.In TV based image processing, the boundary conditions are indispensable to be taken into account for avoiding unwanted artifacts near the boundary. A common technique is assuming a certain behavior of the sharp image outside the boundary. The conventional boundary conditions are zero, periodic, and reflexive boundary conditions. We consider the reflexive boundary condition herein. By denoting∇≔∇1∇2 as the first-order derivative operator, the horizontal direction first-order derivative ∇1 and the vertical direction first-order derivative ∇2 under reflexive boundary condition are defined as (for notational convenience, x∈Rn1×n2 is used for this definition)(17)∇1xi,j=xi+1,j-xi,j,ifi<n1,0,ifi=n1,∇2xi,j=xi,j+1-xi,j,ifj<n2,0,ifj=n2.By denoting C as the discrete cosine transform (DCT) and C-1 as its inverse, ∇T∇ (also the so-called Laplacian operator) under reflective boundary condition can be diagonalized by DCT (see [20]) as(18)∇T∇=C-1ΛDC, where ΛD∈Rn×n is a diagonal matrix. Similarly, the blurring matrix H (under the circumstance that its PSF h is double symmetry) can be diagonalized by DCT (see [21, Chapter 4]) as(19)H=C-1ΛHC,where ΛH∈Rn×n is a diagonal matrix.


## 3. Alternating Direction Method for Model (7)

Consider the following optimization problem with special structure:(20)minf1x+f3zs.t.Ax+Cz=b,x∈X,z∈Z,where f1:Rm1→(-∞,+∞] and f3:Rm3→(-∞,+∞] are l.s.c. proper convex functions; A∈Rl×m1 and C∈Rl×m3 are given matrices; b∈Rl is a known vector; X⊆Rm1 and Z⊆Rm3 are nonempty closed convex sets (we denote f3(z) and C instead of f2(y) and B in (20) as being merely for convenience of analysis in subsequent sections). Suppose that (x∗,z∗) is an optimal solution of (20). For some ξi∗∈∂fi(x∗)(i=1,3), it follows from the first-order optimality condition in optimization that (21)x′-x∗,ξ1∗-ATλ∗≥0,∀x′∈X,z′-z∗,ξ3∗-CTλ∗≥0,∀z′∈Z,λ′-λ∗,Ax∗+Cz∗-b≥0,∀λ′∈Rl,where λ∈Rl is the Lagrange multiplier of (20) [4, 22]. Owing to the special structure of (20), that is, both objective function and linear constraint are separable in variables x and z, the methods based on variable splitting are recommended, for example, the classical ADM which was proposed originally in [23] and was studied intensively in literature (e.g., [24–28]). Numerous practical applications in diversified realms have illustrated its robustness and efficiency (see [29] for an overview of ADM). By introducing the augmented Lagrangian function of (20) on X×Z×Rl as (22)Lx,z,λ=f1x+f3z-λ,Ax+Cz-b+β2Ax+Cz-b2,where β>0 is the penalty parameter of the linear constraint, the recursions of the well-known ADM can be summarized as Algorithm 1.Algorithm 1:Alternating direction method.
Choose arbitraryβ>0, z0∈Rm3 and λ0∈Rl for  k=1,2,…   do (1)xk+1=argminx⁡Lx,zk,λk (2)zk+1=argminz⁡Lxk+1,z,λk (3)λk+1=λk-βAxk+1+Czk+1-bWe now concentrate on dealing with model (7) by ADM, that is, the second phase in [6]. The case p=2 in model (7) is considered in this section. By introducing auxiliary variable v∈Rn, we rewrite model (7) as(23)minτ∇x1+μ2v2+12PAHy-x02s.t.v=x-y.By denoting z≔(y,v), (24)A≔I,C≔-I,-I,b≔0,the linear constraint in (23) can be reformulated as Ax+Cz=b. Furthermore, by setting X≔Rn, Z≔Rn×Rn, and(25)f1x≔τ∇x1,f3z≔μ2v2+12PAHy-x02,the optimization (23) complies with the pattern of optimization (20) and hence can be handled by ADM.We elaborate on the recursions for optimization (23) by utilizing ADM as follows. Firstly, the augmented Lagrangian function for (23) can be specified as (26)Lx,y,v,λ=τ∇x1+μ2v2+12PAHy-x02-λ,x-y-v+β2x-y-v2.(i)
Thex-related subproblem for (23) exploiting ADM amounts to(27)xk+1=argminxτβ∇x1+12x-yk-vk-λkβ2, which can be approximated efficiently by Chambolle’s projection algorithm [17]. (ii)
Thez-related subproblem (involved in the duplets (y,v)) for (23) utilizing ADM is (28)yk+1,vk+1=argminy,v⁡μv2+PAHy-x02+βy+v-xk+1+λkβ2.According to the first-order optimality condition and the fact that PATPA=PA, it follows that (yk+1,vk+1) satisfies (29)μ+βv+βy=βxk+1-λk,βv+βI+HTPAHy=βxk+1-λk+HTPAx0.Consequently, yk+1 and vk+1 can be solved orderly by(30)μ+βHTPAH+μβIyk+1=μ+βHTPAx0+μβxk+1-λk,vk+1=1μHTPAHyk+1-x0.We hence only need to concentrate on tackling the first linear system in (30). It can be solved by preconditioned conjugate gradient (PCG) method as suggested in [6] because its coefficient matrix is positive definite.Remark 1.
Specially, ifA=Ω, the model (7) with p=2 reduces to the model (4) proposed in [5]. Because of PA=I under this circumstance, the x-related subproblem is identical to (27) while the z-related subproblem decays to(31)μ+βHTH+μβIyk+1=μ+βHTx0+μβxk+1-λk,vk+1=1μHTHyk+1-x0.Likewise, PCG method is a better option for the first linear system in (31). Specially, if the PSF of blurring matrix H is doubly symmetric, it can be diagonalized by DCT (see [21]). Hence, linear system (31) can be easily solved.Note that whenp=1 is adopted in model (7), the resulted recursion for x-related subproblem is identical to the case of p=2 in (27). However, the z-related subproblem is quite different with the case of p=2. The z-related subproblem of model (7) with p=1 is (32)yk+1,vk+1=argminy,v⁡μ2v2+PAHy-x01+β2y+v-xk+1+λkβ2,which is tricky to be solved because of the nondifferential term PA(Hy-x0)1 (although it can be solved by iterative methods, the computational effort is intensive). It hence inspires our further investigation in following section.


## 4. Alternating Direction Method with Three Variables for Model (7)

Although the algorithm of ADM in Section3 is valid for model (7) with p=2 and it is favorable to handle model (4) for Gaussian noise removal (i.e., model (7) with p=2 and PA=I), we have deadlock over the model (7) with p=1. Consequently, we contemplate coping with model (7) by ADM with three separable variables.A general extension of optimization (20) is(33)minf1x+f2y+f3zs.t.Ax+By+Cz=b,x∈X,y∈Y,z∈Z,where fi:Rmi→(-∞,+∞], i=1,2,3, are lower semicontinuous proper convex functions; A∈Rl×m1, B∈Rl×m2, and C∈Rl×m3 are given matrices; b∈Rl is a known vector; X∈Rm1, Y∈Rm2, and Z∈Rm3 are nonempty closed convex sets. Suppose that (x∗,y∗,z∗) is an optimal solution of optimization (33). Accordingly, based on the first-order optimality condition in optimization and for some ξi∗∈∂fi(x∗)(i=1,2,3), it follows that (34)x′-x∗,ξ1∗-ATλ∗≥0,∀x′∈X,y′-y∗,ξ2∗-BTλ∗≥0,∀y′∈Y,z′-z∗,ξ3∗-CTλ∗≥0,∀z′∈Z,λ′-λ∗,Ax∗+By∗+Cz∗-b≥0,∀λ′∈Rl,where λ∈Rl is the Lagrange multiplier of (33). By defining the augmented Lagrangian function of (33) as (35)Lx,y,z,λ=f1x+f2y+f3z-λ,Ax+By+Cz-b+β2Ax+By+Cz-b2with β>0 being the penalty parameter of the linear constraint, the recursions of ADM with three separable variables for (33) are summarized as in Algorithm 2.Algorithm 2:Alternating direction method with three separable variables.
Choose arbitraryβ>0, y0∈Rm2, z0∈Rm3 and λ0∈Rl for  k=1,2,…  do (1)xk+1=argminx⁡L(x,yk,zk,λk) (2)yk+1=argminy⁡L(xk+1,y,zk,λk) (3)zk+1=argminz⁡L(xk+1,yk+1,z,λk) (4)λk+1=λk-βAxk+1+Byk+1+Czk+1-bWe now consider utilizing Algorithm2 to solve the model (7). It will be shown that all the resulting subproblems for model (7) exploiting Algorithm 2 possess closed-form solutions.By introducing auxiliary variablesu∈Rn×Rn, v∈Rn, and w∈Rn, the model (7) can be reformulated as(36)minτu1+μ2v2+1pPAwpp,p=1,2,s.t.u=∇xv=x-yw=Hy-x0.By analogical manipulations as in Section 3, we define z≔(u,v,w),(37)A=∇I0,B=0-IH,C=-I000-I000-I,b=00x0.Thereby, the linear constraint in (36) can be realigned as in the form of Ax+By+Cz=b. Furthermore, by denoting X≔Rn, Y≔Rn, and Z≔(Rn×Rn)×Rn×Rn, (38)f1x≔0,f2y≔0,f3z≔τu1+μ2v2+1pPAwpp,the objective function in (36) can be split into three-block functions in variables x, y, and z, respectively. Overall, optimization (36) favors the pattern of (33). Similarly, the augmented Lagrangian function of (36) can be specified as(39)Lx,y,z,λ=τu1+μ2v2+1pPAwpp-λ1,∇x-u-λ2,x-y-v-λ3,Hy-w-x0+β2∇x-u2+β2x-y-v2+β2Hy-w-x02,where λ=(λ1,λ2,λ3)∈(Rn×Rn)×Rn×Rn is Lagrange multiplier.We now elucidate the procedures of handling optimization (36) by Algorithm 2, and it will be shown that all the resulting subproblems possess closed-form solutions (either linear least squares problem or problem involved in soft-thresholding).(i) Thex-related subproblem for (36) utilizing Algorithm 2 can be formulated as (40)xk+1=argminx⁡∇x-uk-λ1kβ2+x-yk-vk-λ2kβ2,which essentially is a linear least squares problem. By first-order optimality condition in optimization, it follows that xk+1 satisfies the linear system (41)β∇T∇+Ixk+1=λ2k+βyk+vk+∇Tλ1k+βuk.As aforementioned, when the first-order partial derivatives operator ∇ is defined as in (17), the above linear system is tractable to be solved by DCT as in (18).(ii) They-related subproblem for (36) exploiting Algorithm 2 can be rewritten as (42)yk+1=argminyxk+1-y-vk-λ2kβ2+Hy-wk-x0-λ3kβ2,which is also a linear least squares problem. Likewise, by the first-order optimality condition, it follows that (43)βI+HTHyk+1=HTλ3k+βwk+x0-λ2k+βxk+1-vk,which is homogeneous to the first linear system in (31). Accordingly, it can be handled as promoted in Remark 1.(iii) Thez-related subproblem for (36) applying Algorithm 2 is involved in the triplets (u,v,w) as follows:(44)zk+1=arg⁡minz⁡τu1+μ2v2+1pPAwpp+β2∇xk+1-u-λ1kβ2+β2xk+1-yk+1-v-λ2kβ2+β2Hyk+1-w-x0-λ3kβ2.It is obvious that the unconstrained optimization (44) is separable in variables u, v, and w. Hence, zk+1=(uk+1,vk+1,wk+1) can be acquired simultaneously (in parallel manner if possible) as follows: (a)
Theu-related subproblem is (45)uk+1=argminu⁡τu1+β2u-∇xk+1+λ1kβ2,and its closed-form solution is obtained by the soft-thresholding defined in (14). Concretely, we have uk+1=Sτ/β∇xk+1-λ1k/β.(b)
Thev-related subproblem is (46)vk+1=argminv⁡μ2v2+β2xk+1-yk+1-v-λ2kβ2,which can be facilely solved byvk+1=(1/(μ+β))[β(xk+1-yk+1)-λ2k].(c)
Thew-related subproblem is (47)wk+1=argminw⁡1pPAwpp+β2Hyk+1-w-x0-λ3kβ2.Combining the definition ofPA, the following follows. If p=2, (48)wk+1i=11+ββHyk+1-x0-λ3ki,ifi∈A,Hyk+1-x0-λ3kβi,ifi∈Ω∖A.Ifp=1, (49)wk+1i=S1/βHyk+1-x0-λ3kβi,ifi∈A,Hyk+1-x0-λ3kβi,ifi∈Ω∖A.Remark 2.
Literally, we take a fixed penaltyβ for the foregoing inferences. In practice, it is attractive to choose β dynamically, for example, the continuation strategy (e.g., [30]) and self-adaptive strategy (e.g., [28]). Moreover, it is reasonable to choose β differently in the last three terms of (39) so as to penalize those three different linear constraints in (36); that is, the augmented Lagrangian function in (39) can be rewritten as (50)Lx,y,z,λ=τu1+μ2v2+1pPAwpp-λ1,∇x-u-λ2,x-y-v-λ3,Hy-w-x0+β12∇x-u2+β22x-y-v2+β32Hy-w-x02,where βi>0, i=1,2,3, are penalties corresponding to different linear constraints. Hence, the recursions for x-, y-, and z-related subproblems can be regained accordingly. Generally, those skills on penalties can numerically render impressive performance (see, e.g., [7, 31, 32] and references therein).Recently, A series variant algorithms for the ADM with more than two variables [33–35] are proposed, which are proved to be convergent and illustrate attractive numerical performances in tackling optimization problems with separable structure. Generally, the three-block ADM may not converge, which was demonstrated by Chen et al. [36] using counterexamples. Nevertheless, if all the objective functions are strongly convex, Han and Yuan [37] proved the global convergence of the three-block ADM scheme. Then, this condition was relaxed and only two or more functions in the objective are required to be strongly convex to ensure the convergence. More recently, it was studied in [38, 39] that the global convergence of the three-block ADM scheme can be guaranteed only when one of the objective functions is strongly convex. However, neither of the three-block functions in the objective function in (36) is strongly convex and then the convergence of Algorithm 2 is still ambiguous. Under some special circumstances (e.g., the matrix involved in the linear constraint is the identity matrix or it is sparse), the performances of some variant algorithms can be matchable with that of Algorithm 2, even fairly close to that of Algorithm 2 [36]. We recommend herein the ADM with Gaussian-back-substitution algorithm proposed in [40] for handling model (7).By denoting variableP≔(y,z,λ) and defining the sparse matrix M as (51)M=I0-BTB-1BTB-1HT00I00000I00000I00000I,where H is blurring matrix and B is defined in (37), ADM with Gaussian-back-substitution algorithm for the concrete optimization (36) is summarized in Algorithm 3.Algorithm 3:Alternating direction method with Gaussian-back-substitution for (36).
Choose arbitraryα∈[0.5,1), β>0, y0∈Rm2, z0∈Rm3 andλ0∈Rl for  k=1,2,…  do (1)xk+1=argminx⁡L(x,yk,zk,λk) (2)y~k=argminy⁡L(xk+1,y,zk,λk) (3)z~k=argminz⁡L(xk+1,y~k,z,λk) (4)λ~k=λk-βAxk+1+By~k+Cz~k-b (5)Pk+1=Pk+αM(P~k-Pk)The augmented Lagrangian functionL(x,y,z,λ) in Algorithm 3 is as defined in (39). The detailed recursions for solving (36) utilizing Algorithm 3 can be inferred analogically as the scenario of using Algorithm 2. The convergence analysis of Algorithm 3 was expanded in [40]. The option on parameter α≈1 is recommended for practical numerical experiments.


## 5. Numerical Experiments

We test the aforementioned algorithms on the fundamental image processing problem, image denoising (both Gaussian noise removal and mixed noise removal). We first concentrate on the scenario of Gaussian noise removal in Section5.1, with numerical comparisons to the algorithm in [5] (abbr, HNW08) and then turn to the case of mixed noise removal in Section 5.2, with numerical comparisons to the algorithm in [6] (abbr, HNW09). We abbreviate hereafter Algorithms 1, 2, and 3 as “ADM2,” “ADM3,” and “ADMGB,” respectively. All the codes are written inmatlab 7.9 and are run on a desktop computer with Intel(R) Core (TM) 2.66 HZ and 2 G memory. All the testing images for numerical experiments in this section are illustrated in Figure 2. To measure the quality of restored image, we define the signal-to-noise ratio (SNR) as (52)SNR=20log10⁡xx¯-x,where x¯ is the restoration of the ideal image x. The stopping criterion for testing algorithms is set as(53)xk+1-xkxk+1<Tol.As mentioned in [12], it is an old and open problem for choosing the optimal parameters in algorithms; we therefore tone those parameters (e.g., τ, μ, and β) by trial and error. As suggested in Remark 2, we opt for the penalty β differently in Algorithms 2 and 3 for various linear constraints.Figure 2
Ideal images:128×128 shape.jpg, 256×256 cameraman.png, 512×512 Barbara.png, 512×512 man.png, 512×512 boat.png, and 512×512 hill.png.


### 5.1. Gaussian Noise Removal

In this subsection, we consider the case of Gaussian noise removal, that is, the model (4). The ideal images are degraded by various blurs and are further contaminated by additive Gaussian noise with noise variance σ2. Those procedures can be realized via fspecial and imnoise inmatlab Image Processing Toolbox (IPT). The concrete parameters for both commands in experiments are summarized in Table 1.Table 1
Problems for Gaussian noise removal.
Blur Uniform Gaussian Out-of-focus Shape H size = 7 hsize = 7, sigma = 2 radius = 3 Cameraman Hsize = 9 hsize = 9, sigma = 2 radius = 7 Barbara Hsize = 11 hsize = 11, sigma = 5 radius = 9We test HWN08 (codes are available athttp://www.math.hkbu.edu.hk/~mng/) and ADM2 in this numerical experiments. We choose the weights τ=10-4 and μ=1 in model (4) for both testing algorithms. The penalty β=2×10-3 is adopted for ADM2. Both HNW08 and ADM2 require Chambolle’s projection algorithm [17] to approximate the solution of x-related subproblem. We take 10 steps for this approximation at each iteration. The initial point x0 is chosen as the observed degraded image and z0=0 and λ0=0. The stopping criterion in (53) is set as Tol=10-4. The numerical results are reported in Table 2. The datum therein, such as “62/1.94/23.91,” represents the iteration number, computing time in seconds and SNR, respectively. We can see that ADM is faster than AM algorithm. The reason is that the dual variable λ is exploited and is updated at each iteration. Figure 3 illustrated some devastated images and their restorations by ADM2.Table 2
Computational results on Gaussian noise removal by solving model (4).
Blur σ2 Algorithm Shape Cameraman Barbara Uniform 0.005 HNW08 62/1.94/23.91 61/7.33/21.77 51/31.63/18.23 ADM2 24/0.75/24.00 26/3.59/21.68 26/17.38/18.19 0.01 HNW08 53/1.67/20.87 51/6.72/19.91 52/32.61/17.70 ADM2 30/0.92/21.69 29/4.05/20.12 31/20.41/17.68  Gaussian 0.005 HNW08 82/2.56/20.89 80/9.70/20.28 58/31.22/18.02 ADM2 35/1.05/21.22 34/4.70/20.39 30/19.88/18.01 0.01 HNW08 58/1.80/17.93 54/6.48/19.53 57/35.38/17.64 ADM2 35/1.08/19.12 36/4.92/19.74 34/22.55/17.61  Out-of-focus 0.005 HNW08 58/1.83/23.37 66/7.91/20.43 61/32.91/18.03 ADM2 20/0.77/23.43 29/3.98/20.43 33/22.11/18.01 0.01 HNW08 48/1.23/19.34 55/6.59/18.81 57/31.56/17.39 ADM2 23/0.80/20.38 34/4.61/19.09 34/22.25/17.46Figure 3
(a) Image devastated by out-of-focus blur and Gaussian noise (σ2=0.01). (b) Restored images by ADM2.
(a) (b)


### 5.2. Mixed Noise Removal

We now test the case of mixed noise removal, that is, the model (7), and report numerical comparisons between HNW09, ADM2, ADM3, and ADMGB. The ideal image of cameraman is degraded by Gaussian blur withhsize = 7 andsigma = 2 and out-of-focus blur with radius = 3, respectively. The blurred images are further corrupted by additive Gaussian noise with noise variance as σ2=0.005 and impulse noise with various noise level s. The median-type filters are attractive in removing impulse noise. We choose AMF for salt-and-pepper noise detection while ACWMF for random-valued noise detection as in [6].We first experiment on the scenario of salt-and-pepper noise plus Gaussian noise removal. The parameters in model (7) are taken as τ=10-4, μ=1, and p=2 for all testing algorithms. Moreover, we choose β=10-3 for ADM2; β1=10-3, β2=10-3, and β3=10-1 for ADM3 and ADMGB. The window size for median-type filter AMF is set as 19 as in [12]. Because the x- and y-related subproblems in both HNW09 and ADM2 should be approximated by Chambolle’s projection method [17] and PCG method, respectively, we set 5-step Chambolle’s approximation for x-related subproblem while 10-4 as tolerance for pcg inmatlab for y-related subproblem. The initial points in experiments are taken y0=0, z0=0, and λ0=0 for all testing algorithms and the stopping criterion in (53) is set as Tol=5×10-4. The numerical results are listed in Table 3. We plot the evolutions of SNR with respect to iteration numbers and computing time in Figure 4. It manifests the efficiency of ADM3 and ADMGB as the resulting subproblems possess closed-form solutions. Figure 5 illustrates some restored images by ADM3 with different impulse noise level. We also test algorithms on other images for salt-and-pepper noise plus Gaussian noise removal. The ideal 512×512 images are degraded by out-of-focus blur with radius = 7 and contaminated by salt-and-pepper noise with noise level 70% plus Gaussian noise with noise variance σ2=0.005. Figure 6 shows the restored images by ADM3.Table 3
Computational results on salt-and-pepper noise plus Gaussian noise removal by solving model (7).
Blur s HNW09 ADM2 ADM3 ADMGB Gaussian 30% 23/39.08/20.65 14/30.72/20.68 45/6.53/20.95 46/10.08/20.92 50% 29/49.66/20.41 15/32.61/20.42 50/7.66/20.57 51/11.28/20.59 70% 36/61.31/20.09 17/36.58/20.10 66/9.95/20.12 67/13.55/20.13 90% 56/89.73/18.74 30/64.38/18.71 143/21.48/18.77 154/28.22/18.62  Out-of-focus 30% 19/32.78/22.66 13/27.95/22.30 46/6.67/22.59 47/10.13/22.48 50% 24/40.27/21.95 14/30.08/21.72 53/7.91/21.96 54/11.20/21.92 70% 35/59.88/20.89 17/37.61/20.74 85/11.89/20.82 94/19.78/20.89 90% 57/97.95/18.43 33/71.48/18.57 144/21.27/18.58 146/30.77/18.42Figure 4
Evolutions of SNR with respect to iteration numbers and computing time (cameraman image degraded by Gaussian blur and contaminated by salt-and-pepper noise (s=70%) plus Gaussian noise (σ2=0.005)).Figure 5
(a) Images degraded by Gaussian blur and contaminated by salt-and-pepper noise (s=30%, 50%, 70%, and 90%) plus Gaussian noise (σ2=0.005). (b) Filtered images by AMF. (c) Restored images by ADM3.
(a) (b) (c)Figure 6
(a) Devastated images by out-of-focus blur (radius = 7) and salt-and-pepper noise (s=70%) plus Gaussian noise (σ2=0.005). (b) Restored images by ADM3. From left to right: Barbara (96/80.39/17.77), man (95/80.13/20.98), boat (95/79.80/21.04), and hill (93/78.55/21.26).
(a) (b)We now experiment on the scenario of random-valued noise plus Gaussian noise removal. Generally, Gaussian noise mixed with random-valued noise is more intractable to be removed than that of Gaussian noise mixed with salt-and-pepper noise. The reason is that it is difficult to distinguish pixels contaminated by Gaussian noise with those contaminated by small value random-valued noise. We tackle model (7) by choosing τ=3×10-3, μ=1, and p=2 for all testing algorithms. We take β=5×10-2 for ADM2 and β1=5×10-2, β2=5×10-2, and β3=1 for ADM3 and ADMGB. The median-type filter ACWMF is run 4 times in phase of noise detection as in [12]. The x- and y-related subproblems in HNW09 and ADM2 are still approximated by Chambolle’s projection method and PCG method as the scenario of random-valued noise plus Gaussian noise removal. Likewise, the initial points are also taken as y0=0, z0=0, and λ0=0 and the stopping criterion in (53) is set as Tol=5×10-4. The numerical results are listed in Table 4. Figure 7 plots the evolutions of SNR with respect to iteration numbers and computing time. Some restored images by ADM3 are illustrated in Figure 8. Those 512×512 images are also tested for random-valued noise plus Gaussian noise removal. All images are degraded by out-of-focus blur with radius = 7 and contaminated by random-valued noise with noise level 40% and Gaussian noise with noise variance σ2=0.005. Figure 9 shows the restored images by ADM3.Table 4
Computational results on random-valued noise plus Gaussian noise removal by solving model (7).
Blur s HNW09 ADM2 ADM3 ADMGB Gaussian 20% 30/40.95/19.08 24/43.78/19.29 42/6.28/19.28 43/8.77/19.25 30% 31/41.80/18.91 24/41.80/19.01 43/6.28/19.14 45/9.64/19.12 40% 31/39.66/18.50 25/42.42/18.61 43/6.33/18.54 45/9.17/18.77 50% 33/41.25/17.40 27/47.44/17.22 47/6.89/17.28 49/10.03/17.26  Out-of-focus 20% 27/35.69/19.46 21/38.78/19.77 38/5.78/19.73 40/8.33/19.73 30% 27/36.95/19.20 22/40.64/19.45 40/5.91/19.40 42/8.80/19.50 40% 30/40.58/18.48 23/42.22/18.95 41/6.19/18.93 45/9.30/18.66 50% 39/48.56/17.05 25/44.56/17.15 45/6.63/17.67 49/10.19/17.20Figure 7
Evolutions of SNR with respect to iteration numbers and computing time (cameraman image degraded by out-of-focus blur and contaminated by random-valued noise (s=40%) plus Gaussian noise (σ2=0.005)).Figure 8
(a) Devastated images by out-of-focus blur (radius = 7) and random-valued noise (s=20%, 30%, 40%, and 50%) plus Gaussian noise (σ2=0.005). (c) Restored images by ADM3.
(a) (b) (c)Figure 9
(a) Devastated images by out-of-focus blur (radius = 7) and random-valued noise (s=40%) plus Gaussian noise (σ2=0.005). (b) Restored images by ADM3. From left to right: Barbara (42/35.17/17.17), man (42/34.89/19.19), boat (41/33.92/18.90), and hill (42/35.13/19.66).
(a) (b)


## 5.1. Gaussian Noise Removal

In this subsection, we consider the case of Gaussian noise removal, that is, the model (4). The ideal images are degraded by various blurs and are further contaminated by additive Gaussian noise with noise variance σ2. Those procedures can be realized via fspecial and imnoise inmatlab Image Processing Toolbox (IPT). The concrete parameters for both commands in experiments are summarized in Table 1.Table 1
Problems for Gaussian noise removal.
Blur Uniform Gaussian Out-of-focus Shape H size = 7 hsize = 7, sigma = 2 radius = 3 Cameraman Hsize = 9 hsize = 9, sigma = 2 radius = 7 Barbara Hsize = 11 hsize = 11, sigma = 5 radius = 9We test HWN08 (codes are available athttp://www.math.hkbu.edu.hk/~mng/) and ADM2 in this numerical experiments. We choose the weights τ=10-4 and μ=1 in model (4) for both testing algorithms. The penalty β=2×10-3 is adopted for ADM2. Both HNW08 and ADM2 require Chambolle’s projection algorithm [17] to approximate the solution of x-related subproblem. We take 10 steps for this approximation at each iteration. The initial point x0 is chosen as the observed degraded image and z0=0 and λ0=0. The stopping criterion in (53) is set as Tol=10-4. The numerical results are reported in Table 2. The datum therein, such as “62/1.94/23.91,” represents the iteration number, computing time in seconds and SNR, respectively. We can see that ADM is faster than AM algorithm. The reason is that the dual variable λ is exploited and is updated at each iteration. Figure 3 illustrated some devastated images and their restorations by ADM2.Table 2
Computational results on Gaussian noise removal by solving model (4).
Blur σ2 Algorithm Shape Cameraman Barbara Uniform 0.005 HNW08 62/1.94/23.91 61/7.33/21.77 51/31.63/18.23 ADM2 24/0.75/24.00 26/3.59/21.68 26/17.38/18.19 0.01 HNW08 53/1.67/20.87 51/6.72/19.91 52/32.61/17.70 ADM2 30/0.92/21.69 29/4.05/20.12 31/20.41/17.68  Gaussian 0.005 HNW08 82/2.56/20.89 80/9.70/20.28 58/31.22/18.02 ADM2 35/1.05/21.22 34/4.70/20.39 30/19.88/18.01 0.01 HNW08 58/1.80/17.93 54/6.48/19.53 57/35.38/17.64 ADM2 35/1.08/19.12 36/4.92/19.74 34/22.55/17.61  Out-of-focus 0.005 HNW08 58/1.83/23.37 66/7.91/20.43 61/32.91/18.03 ADM2 20/0.77/23.43 29/3.98/20.43 33/22.11/18.01 0.01 HNW08 48/1.23/19.34 55/6.59/18.81 57/31.56/17.39 ADM2 23/0.80/20.38 34/4.61/19.09 34/22.25/17.46Figure 3
(a) Image devastated by out-of-focus blur and Gaussian noise (σ2=0.01). (b) Restored images by ADM2.
(a) (b)


## 5.2. Mixed Noise Removal

We now test the case of mixed noise removal, that is, the model (7), and report numerical comparisons between HNW09, ADM2, ADM3, and ADMGB. The ideal image of cameraman is degraded by Gaussian blur withhsize = 7 andsigma = 2 and out-of-focus blur with radius = 3, respectively. The blurred images are further corrupted by additive Gaussian noise with noise variance as σ2=0.005 and impulse noise with various noise level s. The median-type filters are attractive in removing impulse noise. We choose AMF for salt-and-pepper noise detection while ACWMF for random-valued noise detection as in [6].We first experiment on the scenario of salt-and-pepper noise plus Gaussian noise removal. The parameters in model (7) are taken as τ=10-4, μ=1, and p=2 for all testing algorithms. Moreover, we choose β=10-3 for ADM2; β1=10-3, β2=10-3, and β3=10-1 for ADM3 and ADMGB. The window size for median-type filter AMF is set as 19 as in [12]. Because the x- and y-related subproblems in both HNW09 and ADM2 should be approximated by Chambolle’s projection method [17] and PCG method, respectively, we set 5-step Chambolle’s approximation for x-related subproblem while 10-4 as tolerance for pcg inmatlab for y-related subproblem. The initial points in experiments are taken y0=0, z0=0, and λ0=0 for all testing algorithms and the stopping criterion in (53) is set as Tol=5×10-4. The numerical results are listed in Table 3. We plot the evolutions of SNR with respect to iteration numbers and computing time in Figure 4. It manifests the efficiency of ADM3 and ADMGB as the resulting subproblems possess closed-form solutions. Figure 5 illustrates some restored images by ADM3 with different impulse noise level. We also test algorithms on other images for salt-and-pepper noise plus Gaussian noise removal. The ideal 512×512 images are degraded by out-of-focus blur with radius = 7 and contaminated by salt-and-pepper noise with noise level 70% plus Gaussian noise with noise variance σ2=0.005. Figure 6 shows the restored images by ADM3.Table 3
Computational results on salt-and-pepper noise plus Gaussian noise removal by solving model (7).
Blur s HNW09 ADM2 ADM3 ADMGB Gaussian 30% 23/39.08/20.65 14/30.72/20.68 45/6.53/20.95 46/10.08/20.92 50% 29/49.66/20.41 15/32.61/20.42 50/7.66/20.57 51/11.28/20.59 70% 36/61.31/20.09 17/36.58/20.10 66/9.95/20.12 67/13.55/20.13 90% 56/89.73/18.74 30/64.38/18.71 143/21.48/18.77 154/28.22/18.62  Out-of-focus 30% 19/32.78/22.66 13/27.95/22.30 46/6.67/22.59 47/10.13/22.48 50% 24/40.27/21.95 14/30.08/21.72 53/7.91/21.96 54/11.20/21.92 70% 35/59.88/20.89 17/37.61/20.74 85/11.89/20.82 94/19.78/20.89 90% 57/97.95/18.43 33/71.48/18.57 144/21.27/18.58 146/30.77/18.42Figure 4
Evolutions of SNR with respect to iteration numbers and computing time (cameraman image degraded by Gaussian blur and contaminated by salt-and-pepper noise (s=70%) plus Gaussian noise (σ2=0.005)).Figure 5
(a) Images degraded by Gaussian blur and contaminated by salt-and-pepper noise (s=30%, 50%, 70%, and 90%) plus Gaussian noise (σ2=0.005). (b) Filtered images by AMF. (c) Restored images by ADM3.
(a) (b) (c)Figure 6
(a) Devastated images by out-of-focus blur (radius = 7) and salt-and-pepper noise (s=70%) plus Gaussian noise (σ2=0.005). (b) Restored images by ADM3. From left to right: Barbara (96/80.39/17.77), man (95/80.13/20.98), boat (95/79.80/21.04), and hill (93/78.55/21.26).
(a) (b)We now experiment on the scenario of random-valued noise plus Gaussian noise removal. Generally, Gaussian noise mixed with random-valued noise is more intractable to be removed than that of Gaussian noise mixed with salt-and-pepper noise. The reason is that it is difficult to distinguish pixels contaminated by Gaussian noise with those contaminated by small value random-valued noise. We tackle model (7) by choosing τ=3×10-3, μ=1, and p=2 for all testing algorithms. We take β=5×10-2 for ADM2 and β1=5×10-2, β2=5×10-2, and β3=1 for ADM3 and ADMGB. The median-type filter ACWMF is run 4 times in phase of noise detection as in [12]. The x- and y-related subproblems in HNW09 and ADM2 are still approximated by Chambolle’s projection method and PCG method as the scenario of random-valued noise plus Gaussian noise removal. Likewise, the initial points are also taken as y0=0, z0=0, and λ0=0 and the stopping criterion in (53) is set as Tol=5×10-4. The numerical results are listed in Table 4. Figure 7 plots the evolutions of SNR with respect to iteration numbers and computing time. Some restored images by ADM3 are illustrated in Figure 8. Those 512×512 images are also tested for random-valued noise plus Gaussian noise removal. All images are degraded by out-of-focus blur with radius = 7 and contaminated by random-valued noise with noise level 40% and Gaussian noise with noise variance σ2=0.005. Figure 9 shows the restored images by ADM3.Table 4
Computational results on random-valued noise plus Gaussian noise removal by solving model (7).
Blur s HNW09 ADM2 ADM3 ADMGB Gaussian 20% 30/40.95/19.08 24/43.78/19.29 42/6.28/19.28 43/8.77/19.25 30% 31/41.80/18.91 24/41.80/19.01 43/6.28/19.14 45/9.64/19.12 40% 31/39.66/18.50 25/42.42/18.61 43/6.33/18.54 45/9.17/18.77 50% 33/41.25/17.40 27/47.44/17.22 47/6.89/17.28 49/10.03/17.26  Out-of-focus 20% 27/35.69/19.46 21/38.78/19.77 38/5.78/19.73 40/8.33/19.73 30% 27/36.95/19.20 22/40.64/19.45 40/5.91/19.40 42/8.80/19.50 40% 30/40.58/18.48 23/42.22/18.95 41/6.19/18.93 45/9.30/18.66 50% 39/48.56/17.05 25/44.56/17.15 45/6.63/17.67 49/10.19/17.20Figure 7
Evolutions of SNR with respect to iteration numbers and computing time (cameraman image degraded by out-of-focus blur and contaminated by random-valued noise (s=40%) plus Gaussian noise (σ2=0.005)).Figure 8
(a) Devastated images by out-of-focus blur (radius = 7) and random-valued noise (s=20%, 30%, 40%, and 50%) plus Gaussian noise (σ2=0.005). (c) Restored images by ADM3.
(a) (b) (c)Figure 9
(a) Devastated images by out-of-focus blur (radius = 7) and random-valued noise (s=40%) plus Gaussian noise (σ2=0.005). (b) Restored images by ADM3. From left to right: Barbara (42/35.17/17.17), man (42/34.89/19.19), boat (41/33.92/18.90), and hill (42/35.13/19.66).
(a) (b)


## 6. Conclusions

The variable splitting is a hot-investigated strategy in both variational inequalities and optimization. Numerous practical applications illustrate the attractive and efficient performance of this strategy. We illustrate in this paper that the variable splitting can be efficient and promising in image restoration. Note that the numerical effects for mixed noise removal utilizing model (7) with p=1 are similar as that of p=2 (see also [6]). As stated in this paper, we only consider the image denoising with blurry. However, the mixed noise may be occurring in image inpainting, superresolution, and decomposition, even in the area of video processing, which may be difficult to be handled. They are thus some of our future research topics.



---

*Source: 3151303-2016-11-24.xml*