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

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. Let [mathematical expression not reproducible] be an ideal image scaled in [0, 1]. Hereafter, we represent an [n.sub.1] x [n.sub.2] image x as x [member of] [R.sup.n] (n = [n.sub.1][n.sub.2]) by stacking its columns as a vector. H [member of] [R.sup.nxn] is matrix representation of spatially invariant point spread function (PSF) which characterizes the blurry effects on imaging system and n [member of] [R.sup.n] 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 [x.sup.0] e [R.sup.n] involved in the blurring matrix H and the additive noise n can be boiled down to

[x.sup.0] = Hx + n. (1)

The main task of image restoration is to recover the ideal image x from the observed image [x.sup.0]. 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

[mathematical expression not reproducible], (2)

and it is equivalent to (or rather equivalent under a proper [tau], see [4]) the unconstrained model

[mathematical expression not reproducible], (3)

where [tau], [sigma] are parameters depending on noise variance; [nabla] is first-order derivative operator and [parallel][absolute value of [nabla] x ][parallel][sub.1] is the TV-norm (details in Section 2); [parallel] x [parallel][sub.p] denotes the p-norm in Euclidean space (if p = [infinity], [parallel] x [parallel][sup.p.sub.p] herein should be comprehended as [parallel x [parallel][sub.[infinity]]). 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 = [infinity] 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

[mathematical expression not reproducible], (4)

where [tau], [mu] are positive parameters. Intuitively, a new fitting term [parallel] x - y [parallel][sup.2] is added to the unconstrained model (3). Essentially, by smoothing of the nonsmooth function [parallel][absolute value of [nabla]x][parallel][sub.1] via the Moreau approximation, the objective function of the following problem

[mathematical expression not reproducible] (5)

can be replaced by

[mathematical expression not reproducible], (6)

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

[mathematical expression not reproducible], (7)

where A [subset] [OMEGA] := {1, 2, ..., n} contains the noise-free pixels (pixels are not contaminated by impulse noise); [P.sub.A] is the characteristic function of set A. Specifically, for any x [member of] [R.sup.n], j [member of] {1, 2,..., n}, [P.sub.A](x) can be defined as

[mathematical expression not reproducible], (8)

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 [[d.sub.min], [d.sub.max]] be the dynamic interval of image x. For given noise level (also the so-called noise intensity) s [member of] (0,1), salt-and-pepper noise and random-valued noise contaminate image x, respectively, by the following:

(i) Salt-and-pepper noise (see [10])

[mathematical expression not reproducible]. (9)

(ii) Random-valued noise (see [11])

[mathematical expression not reproducible], (10)

where [N.sub.imp] : [R.sup.n] [right arrow] [R.sup.n] represents the procedure of contaminating image x with impulse noise and [d.sub.i] satisfies uniform distribution in interval [[d.sub.min], [d.sub.max]].

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.

The observed image [x.sup.0] is contaminated by mixed noise

[x.sup.0] = [N.sub.imp] (Hx + [n.sub.g]), (11)

where [n.sub.g] is additive Gaussian noise. It is tricky to retrieve the ideal image x in (11) from its observation [x.sup.0] because the statistical distribution of noise in [x.sup.0] 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

[mathematical expression not reproducible], (12)

where [GAMMA] [subset] [OMEGA] 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 [x.sup.0] is the observed image and [y.sup.0] 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 | [y.sup.0.sub.i] [not equal to] [x.sup.0.sub.i], [x.sup.0.sub.i] [member of] {[d.sub.min], [d.sub.max]}, i = 1,2, ..., n}.

(b) random-valued noise: N = {i | [y.sup.0.sub.i] = [x.sup.0.sub.i], 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 = [OMEGA] 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 Section 2, 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. Let x [member of] [R.sup.n] be a vector and p [greater than or equal to] 1 be an integer. By [parallel]x[parallel][sub.p] := [([[summation].sup.n.sub.i=1] [absolute value of [x.sub.i]][sup.p]).sup.1/p], we denote the p-norm of x and [parallel]x[parallel] denotes the Euclidean norm for brevity. The inner product in Euclidean space is denoted (x, x). For any positive definite matrix M [member of] [R.sup.nxn], the matrix norm of x is defined by [parallel]x[parallel][sub.M] := [square root of (x,Mx)]. We denote by I the identity operator. For u = ([u.sub.1], [u.sub.2]) [member of] [R.sup.n] x [R.sup.n], we define a vector [absolute value of u] [member of] [R.sup.n] with the ith component [absolute value of u][sub.i] = [square root of [([u.sub.1]).sup.2.sub.i] + [([u.sub.2]).sup.2.sub.i]. With the definition of p-norm, it follows that [parallel]| u [parallel][sub.1] = [[summation].sup.n.sub.i=1] [absolute value of u][sub.i] is the TV-norm of u. Let f : [R.sup.n] [right arrow] [-[infinity], +[infinity]] be a function. The domain (resp., epigraph) of f is denoted by dom f := {x [member of] [R.sup.n] | f(x) < +[infinity] (resp., epi f := {(x,y) [member of] [R.sup.n] x R | f(x) [less than or equal to] 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 [R.sup.n] x R, and f is a convex function if epi f is convex set in [R.sup.n] x R. The subdifferential of f, denoted by [mathematical expression not reproducible], is defined

[mathematical expression not reproducible]. (13)

For any [beta] > 0, the proximity function of f(u) = [beta][parallel[absolute value of u][parallel][sub.1] is the well-known soft-thresholding (also the so-called shrinkage). With special denotation S : [R.sup.n] x [R.sup.n] [right arrow] [R.sup.n] x [R.sup.n], it is defined as

[S.sub.[beta]] (u) = u- min {[absolute value of u], [beta]} u/[absolute value of u], (14)

where all operations are calculated in componentwise and [(u/[absolute value of u]).sub.i] is taken as zero if [absolute value of u][sub.i] = 0. Let X [subset or equal to] [R.sup.n] be nonempty set. The indicator function of a nonempty set X [subset] [R.sup.n] is denoted by

[mathematical expression not reproducible]. (15)

The proximity operator of f, denoted by [prox.sub.f] : [R.sup.n] [right arrow] [R.sup.n], is defined as [18, 19]

[mathematical expression not reproducible]. (16)

Specially, the projection operator on a nonempty closed convex set V [subset or equal to] [R.sup.n], denoted by [[PI].sub.V], is the proximity operator of [i.sub.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 [mathematical expression not reproducible] as the first-order derivative operator, the horizontal direction first-order derivative [[nabla].sub.1] and the vertical direction first-order derivative [[nabla].sub.2] under reflexive boundary condition are defined as (for notational convenience, [mathematical expression not reproducible] is used for this definition)

[mathematical expression not reproducible]. (17)

By denoting C as the discrete cosine transform (DCT) and [C.sup.-1] as its inverse, [[nabla].sup.T][nabla] (also the so-called Laplacian operator) under reflective boundary condition can be diagonalized by DCT (see [20]) as

[[nabla].sup.T][nabla] = [C.sup.-1][[LAMBDA].sub.D]C, (18)

where [[LAMBDA].sub.D] [member of] [R.sup.nxn] 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

H = [C.sup.-1] [[LAMBDA].sub.H] C, (19)

where [[LAMBDA].sub.H] [member of] [R.sup.nxn] is a diagonal matrix.

3. Alternating Direction Method for Model (7)

Consider the following optimization problem with special structure:

[mathematical expression not reproducible], (20)

where [mathematical expression not reproducible] are l.s.c. proper convex functions; [mathematical expression not reproducible] are given matrices; b [member of] [R.sup.l] is a known vector; [mathematical expression not reproducible] are nonempty closed convex sets (we denote [f.sub.3](z) and C instead of [f.sub.2](y) and B in (20) as being merely for convenience of analysis in subsequent sections). Suppose that ([x.sup.*], [z.sup.*]) is an optimal solution of (20). For some [[xi].sup.*.sub.i] [member of] [partial derivative][f.sub.i]([x.sup.*]) (i = 1,3), it follows from the first- order optimality condition in optimization that

[mathematical expression not reproducible], (21)

where [lambda] [member of] [R.sup.l] 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 x Z x [R.sup.l] as

[mathematical expression not reproducible], (22)

where [beta] > 0 is the penalty parameter of the linear constraint, the recursions of the well-known ADM can be summarized as Algorithm 1.

We 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 [member of] [R.sup.n], we rewrite model (7) as

[mathematical expression not reproducible]. (23)

By denoting z := (y, v),

A := I,

C := [-I,-I],

b := 0, (24)

the linear constraint in (23) can be reformulated as Ax + Cz = b. Furthermore, by setting X = [R.sup.n], Z := [R.sup.n] x [R.sup.n], and

[mathematical expression not reproducible], (25)

the optimization (23) complies with the pattern of optimization (20) and hence can be handled by ADM.
```Algorithm 1: Alternating direction method.

Choose arbitrary [mathematical expression not reproducible]
for k = 1,2, ... do
(1) [mathematical expression not reproducible]
(2) [mathematical expression not reproducible]
(3) [[lambda].sup.k+1] = [[lambda].sup.k] - [beta]
(A[x.sup.k+1] + C[z.sup.k+1] - b)
```

We elaborate on the recursions for optimization (23) by utilizing ADM as follows. Firstly, the augmented Lagrangian function for (23) can be specified as

[mathematical expression not reproducible]. (26)

(i) The x-related subproblem for (23) exploiting ADM amounts to

[mathematical expression not reproducible], (27)

which can be approximated efficiently by Chambolle's projection algorithm [17].

(ii) The z-related subproblem (involved in the duplets (y, v)) for (23) utilizing ADM is

[mathematical expression not reproducible]. (28)

According to the first-order optimality condition and the fact that [P.sup.T.sub.A] [P.sub.A] = [P.sub.A], it follows that ([y.sup.k+1], [v.sup.k+1]) satisfies

[mathematical expression not reproducible]. (29)

Consequently, [y.sup.k+1] and [v.sup.k+1] can be solved orderly by

[mathematical expression not reproducible]. (30)

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, if A = [OMEGA], the model (7) with p = 2 reduces to the model (4) proposed in [5]. Because of [P.sub.A] = I under this circumstance, the x-related subproblem is identical to (27) while the z-related subproblem decays to

[mathematical expression not reproducible]. (31)

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 when p =1is 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

[mathematical expression not reproducible], (32)

which is tricky to be solved because of the nondifferential term [parallel][P.sub.A](Hy - [x.sup.0])[parallel][sub.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 Section 3 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 [P.sub.A] = 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

[mathematical expression not reproducible], (33)

where [mathematical expression not reproducible] are lower semicontinuous proper convex functions; [mathematical expression not reproducible] are nonempty closed convex sets. Suppose that ([x.sup.*], [y.sup.*], [z.sup.*]) is an optimal solution of optimization (33). Accordingly, based on the first-order optimality condition in optimization and for some [[xi].sup.*.sub.i] [member of] [partial derivative][f.sub.i]([x.sup.*]) (i = 1,2,3), it follows that

[mathematical expression not reproducible], (34)

where [lambda] [member of] [R.sup.l] is the Lagrange multiplier of (33). By defining the augmented Lagrangian function of (33) as

[mathematical expression not reproducible] (35)

with [beta] > being the penalty parameter of the linear constraint, the recursions of ADM with three separable variables for (33) are summarized as in Algorithm 2.

We now consider utilizing Algorithm 2 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 variables u [member of] [R.sup.n] x [R.sup.n], v [member of] [R.sup.n], and w [member of] [R.sup.n], the model (7) can be reformulated as

[mathematical expression not reproducible]. (36)
```Algorithm 2: Alternating direction method with three separable
variables.

Choose arbitrary [mathematical expression not reproducible]
for k = 1,2, ... do
(1) [mathematical expression not reproducible]
(2) [mathematical expression not reproducible]
(3) [mathematical expression not reproducible]
(4) [[lambda].sup.k+1] = [[lambda].sup.k] - [beta](A[x.sup.k+1]
+ B[y.sup.k+1] + C[z.sup.k+1] - b)
```

By analogical manipulations as in Section 3, we define z := (u, v, w),

[mathematical expression not reproducible]. (37)

Thereby, the linear constraint in (36) can be realigned as in the form of Ax + By + Cz = b. Furthermore, by denoting X := [R.sup.n], y = [R.sup.n], and Z := ([R.sup.n] x [R.sup.n]) x [R.sup.n] x [R.sup.n],

[mathematical expression not reproducible], (38)

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

[mathematical expression not reproducible]. (39)

where [lambda] = ([[lambda].sub.1], [[lambda].sub.2], [[lambda].sub.3]) [member of] ([R.sup.n] x [R.sup.n]) x [R.sup.n] x [R.sup.n] 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) The x-related subproblem for (36) utilizing Algorithm 2 can be formulated as

[mathematical expression not reproducible], (40)

which essentially is a linear least squares problem. By first-order optimality condition in optimization, it follows that [x.sup.k+1] satisfies the linear system

[mathematical expression not reproducible]. (41)

As aforementioned, when the first-order partial derivatives operator [nabla] is defined as in (17), the above linear system is tractable to be solved by DCT as in (18).

(ii) The y-related subproblem for (36) exploiting Algorithm 2 can be rewritten as

[mathematical expression not reproducible], (42)

which is also a linear least squares problem. Likewise, by the first-order optimality condition, it follows that

[mathematical expression not reproducible], (43)

which is homogeneous to the first linear system in (31). Accordingly, it can be handled as promoted in Remark 1.

(iii) The z-related subproblem for (36) applying Algorithm 2 is involved in the triplets (w, v, w) as follows:

[mathematical expression not reproducible]. (44)

It is obvious that the unconstrained optimization (44) is separable in variables m, v, and w. Hence, [z.sup.k+1] = ([u.sup.k+1], [v.sup.k+1], [w.sup.k+1]) can be acquired simultaneously (in parallel manner if possible) as follows:

(a) The u-related subproblem is

[mathematical expression not reproducible], (45)

and its closed-form solution is obtained by the soft-thresholding defined in (14). Concretely, we have [u.sp.k+1] = [S.sub.[tau]/[beta]]/([nabla][x.sup.k+1] - [[lambda].sup.k.sub.1]/[beta]).

(b) The v-related subproblem is

[mathematical expression not reproducible], (46)

which can be facilely solved by = (1/([mu] + [beta])) [[beta]([x.sup.k+1] - [y.sup.k+1]) - [[lambda].sup.k.sub.2]].

(c) The w-related subproblem is

[mathematical expression not reproducible]. (47)

Combining the definition of [P.sub.A], the following follows. If p = 2,

[mathematical expression not reproducible]. (48)

If p = 1,

[mathematical expression not reproducible]. (49)

Remark 2. Literally, we take a fixed penalty [beta] for the foregoing inferences. In practice, it is attractive to choose [beta] dynamically, for example, the continuation strategy (e.g., [30]) and self-adaptive strategy (e.g., [28]). Moreover, it is reasonable to choose [beta] 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

[mathematical expression not reproducible], (50)

where [[beta].sub.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 variable P := (y, z, X) and defining the sparse matrix M as

[mathematical expression not reproducible], (51)

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 [mathematical expression not reproducible]
and [[lamda].sup.0] [member of] [R.sup.l]
for k = 1,2, ... do
(1) [mathematical expression not reproducible]
(2) [mathematical expression not reproducible]
(3) [mathematical expression not reproducible]
(4) [mathematical expression not reproducible]
(5) [P.sup.k+1] = [P.sup.k] + [alpha]M([[??].sup.k]
- [P.sup.k])
```

The augmented Lagrangian function L(x, y, z, X) 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 [alpha] [approximately equal to] 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 Section 5.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 in matlab 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

SNR = 20 [log.sub.10] [parallel]x[parallel]/ [parallel][bar.x] - x[parallel], (52)

where [bar.x] is the restoration of the ideal image x. The stopping criterion for testing algorithms is set as

[parallel] [x.sup.k+1] - [x.sup.k] [parallel]/ [parallel][x.sup.k+1] [parallel] < Tol. (53)

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., [tau], [mu], and [beta]) by trial and error. As suggested in Remark 2, we opt for the penalty [beta] differently in Algorithms 2 and 3 for various linear constraints.

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 [[sigma].sub.2]. Those procedures can be realized via fspecial and imnoise in matlab Image Processing Toolbox (IPT). The concrete parameters for both commands in experiments are summarized in Table 1.

We test HWN08 (codes are available at http://www.math .hkbu.edu.hk/~mng/) and ADM2 in this numerical experiments. We choose the weights [tau] = [10.sup.-4] and [mu] = 1 in model (4) for both testing algorithms. The penalty [beta] = 2 x [10.sup.-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 [x.sup.0] is chosen as the observed degraded image and [z.sup.0] = 0 and [[lambda].sup.0] = 0. The stopping criterion in (53) is set as Tol = [10.sup.-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 [lambda] is exploited and is updated at each iteration. Figure 3 illustrated some devastated images and their restorations by ADM2.

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 with hsize = 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 [[sigma].sup.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 [tau] = [10.sup.-4], [mu] = 1, and p = 2 for all testing algorithms. Moreover, we choose [beta] = [10.sup.-3] for ADM2; [[beta].sub.1] = [10.sup.-3], [beta] = [10.sup.-3], and [beta] = [10.sup.-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.sup.-4] as tolerance for pcg in MATLAB for y-related subproblem. The initial points in experiments are taken [y.sup.0] = 0, = [z.sup.0] and [[lambda].sup.0] = 0 for all testing algorithms and the stopping criterion in (53) is set as Tol = 5 x [10.sup.-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 x 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 [[sigma].sup.2] = 0.005. Figure 6 shows the restored images by ADM3.

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 [tau] = 3 x [10.sup.-3], [mu] = 1, and p = 2 for all testing algorithms. We take [beta] = 5 x [10.sup.-2] for ADM2 and [[beta].sub.1] = 5 x [10.sup.-2], [[beta].sub.2] = 5 x [10.sup.-2], and [[beta].sub.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 [y.sup.0] = 0, [z.sup.0] = 0, and [[lambda].sup.0] = 0 and the stopping criterion in (53) is set as Tol =5 x [10.sup.-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 x 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 [[sigma].sup.2] = 0.005. Figure 9 shows the restored images by ADM3.

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.

http://dx.doi.org/10.1155/2016/3151303

Competing Interests

The author declares that there are no competing interests.

Acknowledgments

The author was supported by the Natural Science Foundation of China (Grants nos. 11501301 and 11401319), Jiangsu Planned Projects for Postdoctoral Research Funds (Grant no. 1501071B), and the Foundation of Jiangsu Key Lab for NSLSCS (Grant no. 201601).

References

[1] A. Tikhonov and V. Arsenin, Solution of Ill-Posed Problems, John Wiley & Sons, New York, NY, USA, 1977.

[2] D. Mumford and J. Shah, "Optimal approximations by piecewise smooth functions and associated variational problems," Communications on Pure and Applied Mathematics, vol. 42, no. 5, pp. 577-685, 1989.

[3] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms," Physica D, vol. 60, no. 1-4, pp. 259-268, 1992.

[4] R. T. Rockafellar, Convex Analysis, Princeton University Press, Princeton, NJ, USA, 1970.

[5] Y. M. Huang, M. K. Ng, and Y.-W. Wen, "A fast total variation minimization method for image restoration," Multiscale Modeling & Simulation, vol. 7, no. 2, pp. 774-795, 2008.

[6] Y.-M. Huang, M. K. Ng, and Y.-W. Wen, "Fast image restoration methods for impulse and Gaussian noises removal," IEEE Signal Processing Letters, vol. 16, no. 6, pp. 457-460, 2009.

[7] D. R. Han, X. M. Yuan, and W. X. Zhang, "An augmented Lagrangian based parallel splitting method for separable convex minimization with applications to image processing," Mathematics of Computation, vol. 83, no. 289, pp. 2263-2291, 2014.

[8] R. H. Chan, C.-W. Ho, and M. Nikolova, "Salt-and-pepper noise removal by median-type noise detectors and detail-preserving regularization," IEEE Transactions on Image Processing, vol. 14, no. 10, pp. 1479-1485, 2005.

[9] Y. Q. Dong, R. H. Chan, and S. F. Xu, "A detection statistic for random-valued impulse noise," IEEE Transactions on Image Processing, vol. 16, no. 4, pp. 1112-1120, 2007.

[10] H. Hwang and R. A. Haddad, "Adaptive median filters: new algorithms and results," IEEE Transactions on Image Processing, vol. 4, no. 4, pp. 499-502, 1995.

[11] T. Chen and H. R. Wu, "Space variant median filters for the restoration of impulse noise corrupted images," IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, vol. 48, no. 8, pp. 784-789, 2001.

[12] J.-F. Cai, R. H. Chan, and M. Nikolova, "Two-phase approach for deblurring images corrupted by impulse plus Gaussian noise," Inverse Problems and Imaging, vol. 2, no. 2, pp. 187-204, 2008.

[13] S.-J. Ko and Y. H. Lee, "Center weighted median filters and their applications to image enhancement," IEEE Transactions on Circuits and Systems, vol. 38, no. 9, pp. 984-993, 1991.

[14] C. Cocianu and A. Stan, "Neural architectures for correlated noise removal in image processing," Mathematical Problems in Engineering, vol. 2016, Article ID 6153749, 14 pages, 2016.

[15] S. Morillas, V. Gregori, and A. Hervas, "Fuzzy peer groups for reducing mixed Gaussian-impulse noise from color images," IEEE Transactions on Image Processing, vol. 18, no. 7, pp. 1452-1466, 2009.

[16] M.-G. Shama, T.-Z. Huang, J. Liu, and S. Wang, "A convex total generalized variation regularized model for multiplicative noise and blur removal," Applied Mathematics and Computation, vol. 276, pp. 109-121, 2016.

[17] A. Chambolle, "An algorithm for total variation minimization and applications," Journal of Mathematical Imaging and Vision, vol. 20, no. 1-2, pp. 89-97, 2004.

[18] P. L. Combettes and V. R. Wajs, "Signal recovery by proximal forward-backward splitting," Multiscale Modeling & Simulation, vol. 4, no. 4, pp. 1168-1200, 2005.

[19] J.-J. Moreau, "Fonctions convexes duales et points proximaux dans un espace hilbertien," Comptes Rendus de l'Academie des Sciences, vol. 255, pp. 2897-2899, 1962.

[20] M. J. Buckley, "Fast computation of a discretized thin-plate smoothing spline for image data," Biometrika, vol. 81, no. 2, pp. 247-258, 1994.

[21] P. C. Hansen, J. G. Nagy, and D. P. O'Leary, Deblurring Images: Matrices, Spectra, and Filtering, vol. 3 of Fundamentals of Algorithms, SIAM, Philadelphia, Pa, USA, 2006.

[22] D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods, Academic Press, New York, NY, USA, 1982.

[23] D. Gabay and B. Mercier, "A dual algorithm for the solution of nonlinear variational problems via finite element approximation," Computers & Mathematics with Applications, vol. 2, no. 1, pp. 17-40, 1976.

[24] M. Fukushima, "Application of the alternating direction method of multipliers to separable convex programming problems," Computational Optimization and Applications, vol. 1, no. 1, pp. 93-111, 1992.

[25] R. Glowinski and P. Le Tallec, Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics, vol. 9 of SIAM Studies in Applied Mathematics, SIAM, Philadelphia, Pa, USA, 1989.

[26] D. R. Han, H. J. He, H. Yang, and X. M. Yuan, "A customized Douglas-Rachford splitting algorithm for separable convex minimization with linear constraints," Numerische Mathematik, vol. 127, no. 1, pp. 167-200, 2014.

[27] D. R. Han and X. M. Yuan, "Local linear convergence of the alternating direction method of multipliers for quadratic programs," SIAM Journal on Numerical Analysis, vol. 51, no. 6, pp. 3446-3457, 2013.

[28] B. S. He, L.-Z. Liao, D. R. Han, and H. Yang, "A new inexact alternating directions method for monotone variational inequalities," Mathematical Programming, vol. 92, no. 1, pp. 103-118, 2002.

[29] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," Foundations and Trends[R] in Machine Learning, vol. 3, no. 1, pp. 1-122, 2011.

[30] E. T. Hale, W. Yin, and Y. Zhang, "Fixed-point continuation for [l.sub.1]-minimization: methodology and convergence," SIAM Journal on Optimization, vol. 19, no. 3, pp. 1107-1130, 2008.

[31] M. K. Ng, P. Weiss, and X. Yuan, "Solving constrained total-variation image restoration and reconstruction problems via alternating direction methods," SIAM Journal on Scientific Computing, vol. 32, no. 5, pp. 2710-2736, 2010.

[32] Z. N. Zhu, G. Ch. Cai, and Y.-W. Wen, "Adaptive box-constrained total variation image restoration using iterative regularization parameter adjustment method," International Journal of Pattern Recognition and Artificial Intelligence, vol. 29, no. 7, Article ID 1554003, 20 pages, 2015.

[33] C. Chen, R. H. Chan, S. Ma, and J. Yang, "Inertial proximal ADMM for linearly constrained separable convex optimization," SIAM Journal on Imaging Sciences, vol. 8, no. 4, pp. 2239-2267, 2015.

[34] C. H. Chen, M. Li, X. Liu, and Y. Y. Ye, "On the convergence of multi-block alternating direction method of multipliers and block coordinate descent method," https://arxiv.org/abs/1508 .00193.

[35] W. H. Yang and D. R. Han, "Linear convergence of the alternating direction method of multipliers for a class of convex optimization problems," SIAM Journal on Numerical Analysis, vol. 54, no. 2, pp. 625-640, 2016.

[36] C. H. Chen, B. S. He, Y. Y. Ye, and X. M. Yuan, "The direct extension of ADMM for multi-block convex minimization problems is not necessarily convergent," Mathematical Programming, vol. 155, no. 1-2, pp. 57-79, 2016.

[37] D. R. Han and X. M. Yuan, "A note on the alternating direction method of multipliers," Journal of Optimization Theory and Applications, vol. 155, no. 1, pp. 227-238, 2012.

[38] X. J. Cai, D. R. Han, and X. M. Yuan, "On the convergence of the direct extension of ADMM for three-block separable convex minimization models with one strongly convex function," Computational Optimization and Applications, 2016.

[39] M. Li, D. F. Sun, and K.-C. Toh, "A convergent 3-block semi-proximal ADMM for convex minimization problems with one strongly convex block," Asia-Pacific Journal of Operational Research, vol. 32, no. 4, Article ID 1550024, 19 pages, 2015.

[40] B. He, M. Tao, and X. Yuan, "Alternating direction method with Gaussian back substitution for separable convex programming," SIAM Journal on Optimization, vol. 22, no. 2, pp. 313-340, 2012.

Tingting Wu

College of Science, Nanjing University of Posts and Telecommunications, Nanjing 210023, China

Correspondence should be addressed to Tingting Wu; wutt@njupt.edu.cn

Received 9 June 2016; Revised 7 October 2016; Accepted 23 October 2016

Caption: Figure 1: Impulse noise. (a) 50 x 50 white-black ideal image. (b) Contaminated by salt-and-pepper noise with s = 10%. (c) Contaminated by random-valued noise with s = 10%.

Caption: Figure 2: Ideal images: 128 x 128 shape, 256 x 256 cameraman.png, 512 x 512 Barbara.png, 512 x 512 man.png, 512 x 512 boat.png, and 512x512 hill.png.

Caption: Figure 3: (a) Image devastated by out-of-focus blur and Gaussian noise ([[sigma].sup.2] = 0.01). (b) Restored images by ADM2.

Caption: Figure 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 ([[sigma].sup.2] = 70%) plus Gaussian noise

Caption: Figure 5: (a) Images degraded by Gaussian blur and contaminated by salt-and-pepper noise (s = 30%, 50%, 70%, and 90%) plus Gaussian noise ([[sigma].sup.2] = 0.005). (b) Filtered images by AMF. (c) Restored images by ADM3.

Caption: Figure 6: (a) Devastated images by out-of-focus blur (radius = 7) and salt-and-pepper noise (s = 70%) plus Gaussian noise ([[sigma].sup.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).

Caption: Figure 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 ([[sigma].sup.2] = 0.005)).

Caption: 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 ([[sigma].sup.2] = 0.005). (c) Restored images by ADM3.

Caption: Figure 9: (a) Devastated images by out-of-focus blur (radius = 7) and random-valued noise (s = 40%) plus Gaussian noise ([[sigma].sup.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).
```Table 1: Problems for Gaussian noise removal.

Blur         Uniform           Gaussian          Out-of-focus

Shape       Hsize = 7    hsize = 7, sigma = 2     radius = 3
Cameraman   Hsize = 9    hsize = 9, sigma = 2     radius = 7
Barbara     Hsize = 11   hsize = 11, sigma = 5    radius = 9

Table 2: Computational results on salt-and-pepper noise plus
Gaussian noise removal by solving model (7).

Blur           [[sigma].sup.2]   Algorithm       Shape

0.005          HNW08     62/1.94/23.91
0.01           HNW08     53/1.67/20.87

0.005          HNW08     82/2.56/20.89
0.01           HNW08     58/1.80/17.93

HNW08     58/1.83/23.37
0.01           HNW08     48/1.23/19.34

Blur           [[sigma].sup.2]     Cameraman        Barbara

0.005        61/7.33/21.77   51/31.63/18.23
Uniform                          26/3.59/21.68   26/17.38/18.19
0.01         51/6.72/19.91   52/32.61/17.70
29/4.05/20.12   31/20.41/17.68

0.005        80/9.70/20.28   58/31.22/18.02
Gaussian                         34/4.70/20.39   30/19.88/18.01
0.01         54/6.48/19.53   57/35.38/17.64
36/4.92/19.74   34/22.55/17.61

66/7.91/20.43   61/32.91/18.03
Out-of-focus        0.005        29/3.98/20.43   33/22.11/18.01
0.01         55/6.59/18.81   57/31.56/17.39
34/4.61/19.09   34/22.25/17.46

Table 3: Computational results on salt-and-pepper noise plus Gaussian
noise removal by solving model (7).

30%   23/39.08/20.65   14/30.72/20.68    45/6.53/20.95
Gaussian       50%   29/49.66/20.41   15/32.61/20.42    50/7.66/20.57
70%   36/61.31/20.09   17/36.58/20.10    66/9.95/20.12
90%   56/89.73/18.74   30/64.38/18.71   143/21.48/18.77

30%   19/32.78/22.66   13/27.95/22.30    46/6.67/22.59
Out-of-focus   50%   24/40.27/21.95   14/30.08/21.72    53/7.91/21.96
70%   35/59.88/20.89   17/37.61/20.74   85/11.89/20.82
90%   57/97.95/18.43   33/71.48/18.57   144/21.27/18.58

30%   46/10.08/20.92
Gaussian       50%   51/11.28/20.59
70%   67/13.55/20.13
90%   154/28.22/18.62

30%   47/10.13/22.48
Out-of-focus   50%   54/11.20/21.92
70%   94/19.78/20.89
90%   146/30.77/18.42

Table 4: Computational results on salt-and-pepper noise plus Gaussian
noise removal by solving model (7).

20%   30/40.95/19.08   24/43.78/19.29   42/6.28/19.28
Gaussian       30%   31/41.80/18.91   24/41.80/19.01   43/6.28/19.14
40%   31/39.66/18.50   25/42.42/18.61   43/6.33/18.54
50%   33/41.25/17.40   27/47.44/17.22   47/6.89/17.28

20%   27/35.69/19.46   21/38.78/19.77   38/5.78/19.73
Out-of-focus   30%   27/36.95/19.20   22/40.64/19.45   40/5.91/19.40
40%   30/40.58/18.48   23/42.22/18.95   41/6.19/18.93
50%   39/48.56/17.05   25/44.56/17.15   45/6.63/17.67

20%   43/8.77/19.25
Gaussian       30%   45/9.64/19.12
40%   45/9.17/18.77
50%   49/10.03/17.26

20%   40/8.33/19.73
Out-of-focus   30%   42/8.80/19.50
40%   45/9.30/18.66
50%   49/10.19/17.20
```