# Total Variation Image Restoration Method Based on Subspace Optimization

Authors: XiaoGuang Liu; XingBao Gao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The alternating direction method is widely applied in total variation image restoration. However, the search directions of the method are not accurate enough. In this paper, one method based on the subspace optimization is proposed to improve its optimization performance. This method corrects the search directions of primal alternating direction method by using the energy function and a linear combination of the previous search directions. In addition, the convergence of the primal alternating direction method is proven under some weaker conditions. Thus the convergence of the corrected method could be easily obtained since it has same convergence with the primal alternating direction method. Numerical examples are given to show the performance of proposed method finally.


---

## Body

## 1. Introduction

Digital image restoration has a wide application in various areas including Navigation, Aerospace, and Biomedicine (see [1–3] and the references therein). In general, the relationship between the original image f^∈Rp and the observed image g∈Rq is(1)g=Af^+b,where b∈Rq is the additive noise and the spatial-invariant matrix A∈Rq×p represents the degradation system caused by problems such as motion blur, distortion radiation, and distortion wavelets in seismic imaging, p=p1p2 with p1 and p2 being the number of rows and columns, respectively, when images are expressed as a matrix. In particular, matrix A represents Block Toeplitz-plus-Hankel matrix (with Toeplitz-plus-Hankel blocks) or Block Toeplitz matrix (with Toeplitz blocks) when Neumann boundary condition and the zero boundary condition are used, respectively [4].The objective of the image restoration is to estimate the original imagef^ according to some a priori knowledge about the degradation system A, the additive noise b, and the observed image g. However, it often tends to be very ill-conditioned when the reverse process of model (1) is only used to get an ideal estimation f∗. Thus one of the effective ways to solve these problems is to combine some a priori information of the original image and define the regularization solution; that is, f∗ is a minimizer of the following cost (energy) function:(2)Jf=θAf-g+αΦf,where θ:Rq→R is the measure of the difference between Af and g. The regularization term Φ embodies a priori information and a regularization parameter α>0 is used to control the tradeoff between the terms θ and Φ.In general,θ(x)=x22=∑i=1nxi2 (the square of l2-norm of vector x∈Rn) when the additive noise meets Gaussian distribution. θ(x)=x1=∑i=1n|xi| (l1-norm of vector x∈Rn) when the additive noise meets non-Gaussian distributions such as uniform and speckle.The regularization termΦ(f)=∑i∈IϕDif2, where I=1,2,…,p, ϕ:R→R+=t∈R:t≥0 is the potential function, Di:Rp→Rd is the difference operator which can be seen as a d×p matrix and used to create the difference vector between ith pixel and its d neighboring pixels. Here Di=Di1;Di2 with (3)Di1j=0∀j∈I,ifi∈I1=p2,2p2,…,p1p2,Di1i=-1,Di1i+1=1,Di1j=0∀j∉i,i+1,ifi∉I1,(4)Di2j=0∀j∈I,ifi∈I2=p1-1p2+1,p1-1p2+2,…,p1p2,Di2i=-1,Di2i+p2=1,Di2j=0∀j∉i,i+p2,ifi∉I2.In image restoration, the potential functionϕ plays a key role so that it is intensively studied in recent decades [5–20]. Two classes of regularization terms are well known. One is the Tikhonov class [5–7], for example, Φ(f)=∑i∈IDif22. In these cases, the minimum point of (2) can be easily solved due to differentiability. However, these methods tend to make restored images overly smooth and could not protect the clear boundaries.The other class is based on total variation (TV) regularization [8]:(5)Φf=fTV=∑i∈IDif2=∑i∈IDi1f2+Di2f2.In comparison to the first class, an obvious benefit is that the clear boundaries can be recovered well. Thus many effective methods are proposed to solve the TV deblurring problem (see [8–23]). However, they suffer from some numerical difficulties, since the TV regularization is nondifferentiable. In order to overcome the drawback, in [16], the TV regularization is modified as (6)Φf=fTV,ε=∑i∈IDi1f2+Di2f2+ε, where 0<ε≪1. Obviously, the regularization is differentiable and many classical optimization methods can be applied [24, 25]. Unfortunately, experimental results show that ε must be small enough to keep the restoration quality, but the efficiency of algorithms will be reduced with smaller ε [16, 20].By using the following modified TV regularization,(7)Φf=∑i∈IDi1f+Di2f, a primal-dual active set method is proposed in [17], but it has no rotational invariant [8].Recently, the alternating direction method [26] is used to solve the l2-TV and l1-TV image restoration models, and better performance has been obtained [18–23, 27]. It divides original problem into some simple subproblems and solves them alternatively. The downside is that the search directions are not accurate enough which could reduce the performance of algorithms (see the analysis in Section 3).Based on the considerations above, we will take the algorithms in [18, 20] as examples and present a method to improve the performance of the alternating direction method used in l2-TV and l1-TV models. By means of the subspace optimization [28], it corrects the search direction of the primal alternating direction method by using the energy function and a linear combination of the previous search directions. In addition, the convergence of the primal alternating direction method is proven under some weaker conditions. Thus the convergence of the corrected method could be easily obtained by the equivalence between them. Numerical examples are given to shown the performance of proposed method.The outline of this paper is as follows. Some preliminaries are stated in Section2. A l1-TV image restoration algorithm based on subspace optimization is proposed in Section 3. The convergence analysis and numerical examples are given in Sections 4 and 5, respectively. Finally, some concluding remarks are given in Section 6.


## 2. Preliminaries

In this section, we first propose some basic definitions and properties; then the alternating direction method to solvel1-TV model in [20] will be introduced. The alternating direction method in [18] (l2-TV model) will be omitted here since it is more simpler. But the comparison experimental results of the method and its corrected method still will be listed in Section 5.2. In the following, let Dh and Dv be the one-side difference matrix on the horizontal direction and vertical direction, respectively, D=Dh,Dv⊤, Null(A) is the null space of matric A, xi is the ith iteration, and x(i) is the ith element of vector x.Definition 1 (see [20]).
An operatorP is called nonexpansive if for any x1, x2∈X⊂Rn, we have (8)Px1-Px22≤x1-x22.Specially, P is called β-averaged nonexpansive if there exists some nonexpansive operator A and β∈(0,1) such that P=(1-β)I+βA, where I is the identity operator.Lemma 2 (see [20]).
Ifρ>0, then the minimizer of ψs=x-s2+ρ|s| is given by the following:(9)νx=x-ρ2,x>ρ2,0,x≤ρ2,x+ρ2,x<-ρ2.We know from [20] that the operator ν(x) is nonexpansive.Lemma 3 (see [29]).
Letφ be convex and semicontinuous, β>0 and (10)x^=arg⁡minx⁡y-x22+βφx.Define S such that x^=S(y) for each y; then S is (1/2)-averaged nonexpansive.Definition 4 (see [25]).
A functionφ:Rn→R (i) is said to be proper over a set X⊂Rn if φ(x)<+∞ for at least one x∈X and φ(x)>-∞ for all x∈X and (ii) is said to be coercive over a set X⊂Rn if for every sequence {xk}⊂X such that xk2→∞ we have (11)limk→∞⁡φxk⁡=∞.When X=Rn, we say that φ is coercive on Rn.Lemma 5 (see [30]).
Letφ:Rn→R be a closed, proper, and coercive function. Then the set of minima of φ over Rn is nonempty and compact.
In [20], the l1-TV image restoration model is(12)arg⁡minf⁡Af-g1+αfTV.To avoid the numerical difficulties caused by the nondifferentiability, two auxiliary variables ω∈Rq and η∈Rp are used to cope with the nonsmooth terms Af-g1 and fTV. Then (12) has been transformed as(13)arg⁡minf,ω,η⁡Ff,ω,η=Af-ω22+αηTV+α1ω-g1+α2η-f22,where the parameters α1>0 and α2>0 are used to ensure the closeness of g and ω, f and η, respectively. When applied to (13), the alternating direction method is(14)HAωi,ηi≔fi+1=arg⁡minf⁡Af-ωi22+α2ηi-f22,Hl1fi+1≔ωi+1=arg⁡minω⁡Afi+1-ω22+α1ω-g1,HTVfi+1≔ηi+1=arg⁡minη⁡α2η-fi+122+αηTV.
The first step in (14) is equivalent to solving the nonsingular linear system: (15)A⊤A+α2If=A⊤ωi+α2ηi.It can be solved directly; that is, (16)f=A⊤A+α2I-1A⊤ωi+α2ηi.To let computational cost be low, many classical optimization and numerical methods such as CG method and preconditioned iteration methods can be used to solve it [24, 25, 31]. As in [4, 18], we suppose that the Neumann boundary conditions are used and the blurring function is symmetric in this paper; then A is a (block) Toeplitz-plus-Hankel matrix and it could be diagonalized by a discrete cosine transform matrix. Thus we could solve the nonsingular linear system utilizing lower cost, since the inverse of blurring matrix A can be computed by using fast cosine transforms; please see [4, 18, 27] for more details.Since(17)Afi+1-ω22+α1ω-g1=∑k=1qAfi+1k-ωk2+α1ωk-gk,the second step in (14) is equivalent to solving the minimizers of q functions in the form of ν(s)=x-s2+α1|s|. We can know from Lemma 2 that the kth element of its solution ωi+1 is (18)ωki+1=Afi+1k-α12,Afi+1-gk>α12,gk,Afi+1-gk≤α12,Afi+1k+α12,Afi+1-gk<-α12.The third step in (14) is a l2-TV model with the blurring matrix A is the identity matrix, that is, image denoising problem. Many methods can solve it effectively [13–15], and Chambolle’ projection algorithm [13] is used in [20].


## 3. The Method Based on Subspace Optimization

Based on the principle of alternating direction method above, it is easy to know that the solutionli+1=[fi+1;ωi+1;ηi+1] of (14) is different from the solution of (13), so the performance of solving the minima of (13) may be lowered. Thus a method based on subspace optimization [28] will be proposed to improve the optimization performance of (14) in this section.Firstly, we introduce subspace optimization method. Letl¯i+1=[f¯i+1;ω¯i+1;η¯i+1] denote the solution corrected by subspace optimization method, and (19)Si+1=t0i+1,t1i+1,…,tmi+1,where t0i+1 = li+1 − l¯i and tki+1 = l¯i+1-k − l¯i-k  (1≤k≤m≤i), then(20)l¯i+1=l¯i+Si+1πi+1,(21)πi+1=arg⁡minπ⁡Fl¯i+Si+1π.Supposing πi+1∈Rm+1 is selected as (1,0,…,0) in (20), then the equation l¯i+1=li+1 is valid. In fact, we can know from (21) that πi+1 is a minimum point of the cost function: F(πi+1) = F(l¯i+Si+1π); thus F(l¯i+1)≤F(li+1) must be true. It is helpful to improve the performance of original alternating direction method on every iteration, and the key problem is how to solve (21) now. To reduce the computational cost, let Si+1=[t0i+1] in this paper.Next, the correction method of (14) will be proposed according to three cases.Case 1 (ω(i)≠g(i) and Di1η≠0 or Di2η≠0 for ∀i∈I).
Obviously, the functionF(f,ω,η) is smooth; we have (22)∇F=2A⊤Af-ω+2α2f-η2ω-Af+α1∑i∈Iωi-giωi-giei2α2η-f+α∑i∈Ixi-1/2yizi⊤,∇2F=2A⊤A+2α2I-2A⊤-2α2I-2A2I0-2α2I02α2I+∂2ηTV∂η2,where ei denotes the unit column vector with its ith element being 1 and the others being zero, xi=(Di1η)2+(Di2η)2, yi=Di1η+Di2η, and zi=Di1+Di2. It is expensive to compute ∂2ηTV/∂η2, so we let ∂2ηTV/∂η2=0 to improve the efficiency of the algorithm because ∇2F is independent of variable l now.
On the other hand, sinceπi+1 is an optimum of the convex function F, we can know from (21) and Taylor series expansion that(23)Si+1⊤∇Fl¯i+Si+1πi+1=0,∇Fl¯i+Si+1πi+1=∇Fl¯i+∇2FSi+1πi+1.Equation (23) implies that(24)Si+1⊤∇Fl¯i+∇2FSi+1πi+1=0.Obviously, vector Si+1≠0; then (25)12Si+1⊤∇2FSi+1=Afi+1-f¯i2+ωi+1-ω¯i2-2ωi+1-ω¯i⊤Afi+1-f¯i+α2fi+1-f¯i2+ηi+1-η¯i2-2ηi+1-η¯i⊤fi+1-f¯i=Afi+1-f¯i-ωi+1+ω¯i22+α2fi+1-f¯i-ηi+1+η¯i22≥0.Thus (Si+1)⊤∇2FSi+1 = 0⇔Afi+1-f¯i-ωi+1+ω¯i22 + α2fi+1-f¯i-ηi+1+η¯i22 = 0⇔Afi+1-f¯i = ωi+1-ω¯i and fi+1-f¯i = ηi+1-η¯i. The two equalities do not hold in image restoration since the dimension of vector f is very huge. Therefore (Si+1)⊤∇2FSi+1>0, and(26)πi+1=-Si+1⊤∇2FSi+1-1Si+1⊤∇Fl¯ifrom (24).Case 2 (I3=i∣Di1η=Di2η=0,i∈I is nonempty).
Now,ηTV is nonsmooth. Let (27)ηTV+ε=∑i∈I3Di1η2+Di2η2+ε1/2+∑i∈I/I3Di1η2+Di2η21/2,then limε→0⁡ηTV+ε = ηTV and ∂ηTV+ε/∂ηi=(xi+ε)-1/2yizi⊤=0  (yi=0) when i∈I3. Thus ∂ηTV/∂η(i)=0 is set when i∈I3.Case 3 (∃i∈I such that ω(i)=g(i)).
The termα1ω(i)-g(i)1 is not differential at this case; the Huber function is often used to replace it [25]. For simplicity, ∂ω-g1/∂ω(i)=0 is selected when ω(i) = g(i).
In summary, the proposed algorithm is as follows.
AlgorithmStep 1.
Initializel¯0.Step 2.
Solve(i)
fi+1=arg⁡minf⁡Af-ωi22+α2ηi-f22,(ii)
ωi+1=arg⁡minω⁡Afi+1-ω22+α1ω-g1,(iii)
ηi+1=arg⁡minη⁡α2η-fi+122+αηTV.Step 3.
Computeπi+1 using (26), and let l¯i+1=l¯i+Si+1πi+1.Step 4.
Iffi+1-fi2/fi+12<10-4, then stop; otherwise go to Step 2.


## 4. Convergence Analysis

As known from (20) and (26), when limi→∞t0i2=0 and one of the sequences li and l¯i is convergent, then the other sequence is also convergent, where t0i = li − l¯i-1. Thus a simple proof of the convergence of algorithm (14) will be given below without the condition that Hl1(HA(·,η)) and HTV(HA(ω,·)) are asymptotically regular which is needed in [20], where(28)ωi+1=Hl1fi+1=Hl1HAωi,ηi=T1ωi,ηi,ηi+1=HTVfi+1=HTVHAωi,ηi=T2ωi,ηi,(29)T1·=Hl1HA·,η,T2·=HTVHAω,·.Lemmas 6 and 7 below appeared in [18, 20], but the proofs given by us are simpler.Lemma 6.
SupposeNull(A)∩Null(D)={0}; then the sets of fixed points of operators T1 and T2 are nonempty, respectively.Proof.
It is easy to see that the objection functionF in (13) satisfies the conditions of Lemma 5 [18, 20]. Then F has at least one minimizer (f¯,ω¯,η¯) that cannot be decreased by the alternating scheme (14). Thus (30)f¯=HAω¯,η¯=arg⁡min⁡F·,ω¯,η¯,ω¯=Hl1f¯=arg⁡min⁡Ff¯,·,η¯,η¯=HTVf¯=arg⁡min⁡Ff¯,ω¯,·,(31)ω-=Hl1HAω-,η-=T1ω-,η-=HTVHAω-,η-=T2η-.Therefore ω-, η- are the fixed points of T1 and T2, respectively.Lemma 7.
The operatorsT1 and T2 are nonexpansive.Proof.
From Lemmas2 and 3, we have(32)T1ω1-T1ω22=Hl1HAω1,η-Hl1HAω2,η2=ψAHAω1,η-ψAHAω2,η2≤AHAω1,η-AHAω2,η2=AA⊤A+α2I-1A⊤ω1+α2η-AA⊤A+α2I-1A⊤ω2+α2η2=AA⊤A+α2I-1A⊤ω1-ω22≤ω1-ω22,T2η1-T2η22=HTVHAω,η1-HTVHAω,η22≤HAω,η1-HAω,η22=AA⊤A+α2I-1A⊤ω+α2η1-AA⊤A+α2I-1A⊤ω+α2η22=AA⊤A+α2I-1A⊤η1-η22≤η1-η22.
This completes the proof.Lemma 8.
Let{ωi} and {ηi} be generated by (28); then ∑k=1∞ωik+1-ωik22 and ∑k=1∞‍ηik+1-ηik22 are bounded, where ωik and ηik are the subsequence of {ωi} and {ηi}, respectively.Proof.
ForF(f,ωik,ηik) in (13), we have (33)Ffik,ωik,ηik=Ffik+1,ωik,ηik+fik-fik+1⊤∂F∂ffik+1,ωik,ηik+12fik-fik+1⊤A⊤A+α2Ifik+1,ωik,ηikfik-fik+1.Since fik+1 is the minimizer of F(f,ωik,ηik), so (34)∂F∂ffik+1,ωik,ηik=0.Then (35)Ffik,ωik,ηik-Ffik+1,ωik,ηik≥α2fik+1-fik22.Notice that F(fik+1,ωik+1,ηik+1) ≤ F(fik+1,ωik,ηik), we get (36)Ffik,ωik,ηik-Ffik+1,ωik+1,ηik+1≥α2fik+1-fik22.From Lemma 2, Hl1 is nonexpansive; that is, (37)fik+1-fik22≥Hl1fik+1-Hl1fik22=ωik+1-ωik22,and thus (38)Ffik,ωik,ηik-Ffik+1,ωik+1,ηik+1≥α2ωik+1-ωik22.On the other hand, (39)Ffik+1,ωik+1,ηik+1≤Ffik+1,ωik+1,ηik+1,and we have (40)∑k=1∞Ffik,ωik,ηik-Ffik+1,ωik+1,ηik+1=Ffi1,ωi1,ηi1+∑k→∞Ffik+1,ωik+1,ηik+1-Ffik+1,ωik+1,ηik+1-Ffik+1+1,ωik+1+1,ηik+1+1<Ffi1,ωi1,ηi1.Namely, ∑k=1∞‍ωik+1-ωik22 is bounded. Similarly, one can verify that ∑k=1∞‍ηik+1-ηik22 is also bounded.
Now we state and prove the convergence result of (14).Theorem 9.
SupposeNull(A)∩Null(D)={0}; the sequence {(ωi,ηi)} generated by (14) with any initial point (ω0,η0) converges to a solution (ω∗,η∗) of (13).Proof.
From Lemma6, let ω¯ be any fixed point of T1; then(41)ωi-ω¯2=T1ωi-1-T1ω¯2≤ωi-1-ω¯2≤⋯≤ω0-ω¯2from Lemma 7. Thus the sequence {ωi} is bounded, and there exists a subsequence {ωik} with limk→∞ωik=ω∗. On the other hand, ∑k=1∞‍ωik+1-ωik22 is bounded from Lemma 8; thus (42)limk→∞ωik+1-ωik22=0.Let k→∞, then T1ω∗-ω∗2=0; that is, ω∗ is also a fixed point of T1. Similarly, to prove (41), one can also verify that the sequence ωi-ω∗2 is nonincreasing. Thus limk→∞ωi-ω∗2=0 since limk→∞ωik=ω∗.
Following the argument above, one can prove that the sequence{ηi} generated by (14) with any initial point η0 converges to a solution η∗ of (13).Remark 10.
Ifx∈Null(D), then x(i)=c for all i∈I, where c is a nonzero constant. Since matrix A≥0 in image restoration, so Ax is a nonzero vector. It follows that the assumption Null(A)∩Null(D)={0} holds in general.


## 5. Experimental Results

In this section, some numerical results will be provided to show the performance of our method. As in [18, 20], the tested images are selected as Cameraman of 256 × 256 and Lena of 512 × 512. The algorithm in [20] and ours will be compared in Section 5.1 (l1-TV model). The algorithm in [18] and its corrected algorithm will be compared in Section 5.2 (l2-TV model). All the computational tasks are performed using MATLAB 2016a with Core(TM)2CPU with 2.83 GHz and 3.87 GB of RAM. The average value of ten tests would be selected.The tested blurring function is chosen to be truncated 2D Gaussian function:(43)hs,t=exp⁡-s2-t22σ2,-3≤s,t≤3.Here three sets of parameters are chosen: (1) the support being equal to 5×5(σ=1); (2) the support being equal to 7×7(σ=1.5); (3) the support being equal to 9×9(σ=2).In all runs, CPU time is used to compare the efficiency, signal-to-noise ratio (SNR), and peak signal-to-noise ratio (PSNR):(44)SNR=20log10⁡f^2fi+1-f^2,PSNR=-20log10⁡fi+1-f^2p1p2,which are used to measure the quality of the restored images. The stopping criterion of the algorithms should satisfies (45)fi+1-fi2fi+12≤10-4.


### 5.1. Comparison Experiment for thel1-TV Model

In this subsection, a comparison of the algorithm in [20] with the proposed method is made under different non-Gaussian additive noises. Firstly, for the uniform noise, we let α = 0.008, α1 = 0.05, and α2 = 0.2 as in [20]. The uniform random numbers appear in the intervals (0, 0.05), (0, 0.1), and (0, 0.2), respectively. Specially, [20] points out that the noise added to the blur image is independent and identically distributed, so the same set of parameters could be selected for different images when the same kind of noise is considered. This could reduce the computational cost of searching regularization parameters. For the speckle noise, as in [20, 23], let α ≡ 1, α1 = 2.5 × 104, the continuation scheme about α2 will be used, easier subsequences with smaller α2 can be solved quickly, and the later subproblems can also be solved relatively quickly with warm starts from previous solutions. The noise variance will be selected as 0, 0.01, and 0.05 respectively. To be fair, the same methods will be used to solve the primal alternating direction method in this paper.Now, we select two cases that have significant improvement of the visual sense to illustrate our method’s restoration quality. Under different conditions, more image restoration results with uniform noise and speckle noise will be summarized in Tables1 and 2, respectively. Figure 1(a) shows the original Camera image. Figure 1(b) shows the observed Camera image, where the noise is the uniform noise, the random numbers are in the interval (0, 0.2), and the support is equal to 9 × 9. Figure 1(c) shows the restored Camera image by the algorithm in [20]. Figure 1(d) shows the restored Camera image by our method. Figure 2(a) shows the original Lena image. Figure 2(b) shows the observed Lena image, where the noise is the speckle noise, the noise variance is 0.01, the support is equal to 7×7. Figure 2(c) shows the restored Lena image by the algorithm in [20]. Figure 2(d) shows the restored Lena image by our method.Table 1
Summary results for the uniform noise.
Image
Noise
Blur
SNR
PSNR
CPU
[20]
Ours
[20]
Ours
[20]
Ours
Camera
(0, 0.05)
5 × 5
21.25
21.46
26.83
26.84
5.49
4.84
7 × 7
19.35
19.37
24.93
24.92
6.82
5.99
9 × 9
18.28
18.30
23.86
23.88
8.47
7.99
(0, 0.1)
5 × 5
18.12
18.13
23.70
23.73
6.12
5.93
7 × 7
17.20
17.21
22.78
22.77
8.14
7.09
9 × 9
16.52
16.52
22.11
22.16
8.95
8.17
(0, 0.2)
5 × 5
12.86
12.86
18.45
18.57
7.57
7.07
7 × 7
12.90
12.91
18.48
18.48
7.82
7.00
9 × 9
12.81
12.88
18.39
18.46
9.14
8.03
Lena
(0, 0.05)
5 × 5
24.29
24.31
29.94
29.93
29.84
27.60
7 × 7
23.42
23.46
29.07
29.11
39.13
35.62
9 × 9
22.60
22.67
28.26
28.26
44.05
37.16
(0, 0.1)
5 × 5
19.37
19.37
25.03
25.07
30.50
24.35
7 × 7
19.18
19.15
24.84
24.84
38.24
33.65
9 × 9
18.90
18.92
24.55
24.57
45.47
40.36
(0, 0.2)
5 × 5
13.65
13.74
19.31
19.36
44.06
39.56
7 × 7
13.58
13.63
19.24
19.33
41.28
36.33
9 × 9
13.25
13.25
18.99
18.97
45.49
40.68Table 2
Summary results for the speckle noise.
Image
Noise
Blur
SNR
PSNR
CPU
[20]
Ours
[20]
Ours
[20]
Ours
Camera
0
5 × 5
22.06
22.12
27.64
27.64
5.09
4.36
7 × 7
19.94
19.93
25.52
25.57
7.55
6.05
9 × 9
18.72
18.74
23.89
23.90
8.13
6.88
0.01
5 × 5
20.32
20.33
25.90
25.93
5.88
5.02
7 × 7
19.22
19.27
24.83
24.84
8.30
6.91
9 × 9
18.28
18.28
24.30
24.35
8.13
7.04
0.05
5 × 5
15.77
15.79
21.35
21.35
9.48
8.25
7 × 7
14.95
14.96
20.54
20.55
9.16
7.25
9 × 9
12.63
12.62
18.21
18.24
9.87
8.23
Lena
0
5 × 5
28.33
28.34
33.99
34.01
29.94
27.61
7 × 7
26.07
26.10
31.72
31.72
37.44
32.56
9 × 9
24.59
24.62
30.25
30.61
45.29
38.65
0.01
5 × 5
23.64
23.69
29.30
29.32
36.86
32.58
7 × 7
23.86
23.94
28.92
29.00
42.32
36.32
9 × 9
23.24
23.24
28.89
28.92
50.97
43.65
0.05
5 × 5
18.23
18.25
23.89
23.90
61.20
55.63
7 × 7
16.58
16.57
22.24
22.31
61.31
56.86
9 × 9
16.48
16.52
20.14
20.14
71.48
63.53Figure 1
The restored Camera images by different methods with the uniform noise. The random numbers are in the interval (0, 0.2), the support being equal to9×9. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [20]. (d) Image restored by our method.
(a)
(b)
(c)
(d)Figure 2
The restored Lena images by different methods with the speckle noise. The noise variance is 0.01, the support being equal to7×7. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [20]. (d) Image restored by our method.
(a)
(b)
(c)
(d)From Figures1 and 2, we see that the quality of images restored by our method is improved obviously compared to those restored by the method of [20] in the two cases. For the uniform noise and speckle noise, Tables 1 and 2 show that the SNR and PSNR of the images restored by our method both are higher than the results of [20] in most cases. For the same requirement of relative error, our proposed algorithm is clearly faster than the competing algorithm in [20]. Referring to Sections 2 and 3, the proposed algorithm uses subspace optimization method to correct the search direction, which could ensure F(l¯i+1)≤F(li+1) is valid on the ith iteration. So it is clear that the proposed method is more efficient. The analysis is verified by the experiments results above.


### 5.2. Comparison Experiment for thel2-TV Model

In the subsection, the algorithm in [18] and ours will be compared. Being different from [19, 20], [18] involves a fitting of the auxiliary variable to Df which has superior performance compared to fitting f [18, 23, 32]. The additive noise is selected as Gaussian noise. Here four sets of the standard deviation of noise are chosen: 0.01, 0.05, 0.1, and 0.2. Same methods appearing in [18] will be used to solve the primal alternating direction method.Similarly, we will first select two cases that have significant improvement of the visual sense to illustrate our method’s restoration quality. Figures3(a) and 4(a) show the original images. Figures 3(b) and 4(b) show the observed images. Figure 3(c) shows the restored Camera image by algorithm in [18]; Figure 3(d) shows the restored Camera image by our method, where the standard deviation of noise is 0.1 and the support is equal to 9 × 9. Figure 4(c) shows the restored Lena image by algorithm in [18]; Figure 4(d) shows the restored Lena image by our method, where the standard deviation of noise is 0.2 and the support is equal to 9×9. Table 3 shows more image restoration results with the different support and standard deviation of noise.Table 3
Summary results for Gaussian noise.
Image
Noise
Blur
SNR
PSNR
CPU
[18]
Ours
[18]
Ours
[18]
Ours
Camera
0.01
5 × 5
21.99
21.98
27.57
27.61
2.45
2.31
7 × 7
19.90
19.93
25.49
25.49
2.67
2.32
9 × 9
18.69
18.69
24.27
24.28
2.73
2.63
0.05
5 × 5
20.64
20.69
26.22
26.24
1.83
1.62
7 × 7
19.04
19.12
24.62
24.66
7.60
7.02
9 × 9
18.01
18.07
23.59
23.59
2.57
2.26
0.1
5 × 5
17.97
17.99
23.55
23.58
1.84
1.56
7 × 7
17.06
17.11
22.64
22.63
2.20
2.12
9 × 9
16.41
16.40
22.00
22.61
2.22
1.98
0.2
5 × 5
13.37
17.38
18.95
18.94
1.78
1.51
7 × 7
13.11
13.10
18.70
18.73
1.68
1.53
9 × 9
12.91
12.93
18.49
18.61
1.86
1.68
Lena
0.01
5 × 5
28.05
28.10
33.71
33.73
11.89
10.32
7 × 7
25.90
25.91
31.56
31.62
12.93
10.33
9 × 9
24.47
24.49
30.13
30.15
13.32
10.92
0.05
5 × 5
24.09
24.06
29.74
29.77
10.39
9.66
7 × 7
23.12
23.12
28.78
28.79
11.26
9.65
9 × 9
22.33
22.42
27.99
28.12
13.18
11.56
0.1
5 × 5
19.50
19.55
25.16
25.17
9.14
8.68
7 × 7
19.19
19.23
24.85
24.84
10.26
9.33
9 × 9
18.86
18.88
24.52
24.51
11.79
9.91
0.2
5 × 5
13.81
13.81
19.47
19.44
7.14
6.56
7 × 7
13.83
13.86
19.49
19.83
8.47
7.82
9 × 9
13.78
14.01
19.44
19.87
9.78
7.97Figure 3
The restored Camera images by different methods with Gaussian noise. The standard deviation of noise is 0.1; the support is9×9. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [18]. (d) Image restored by our method.
(a)
(b)
(c)
(d)Figure 4
The restored Lena images by different methods with Gaussian noise. The standard deviation of noise is 0.2; the support is 9 × 9. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [18]. (d) Image restored by our method.
(a)
(b)
(c)
(d)As in Section5.1, Figures 3 and 4 show that the quality of restoration images by using our method is better than those restored by the algorithm of [18] in the two cases. We also can know from Table 3 that the SNR and PSNR of the images restored by our method both are higher than the results of [18] in most cases. With the same reasons mentioned in l1-TV model, all computational time required by our method is significantly less than that required by the algorithm in [18]. So the performance of the algorithm in [18] also has been improved through the correction.


## 5.1. Comparison Experiment for thel1-TV Model

In this subsection, a comparison of the algorithm in [20] with the proposed method is made under different non-Gaussian additive noises. Firstly, for the uniform noise, we let α = 0.008, α1 = 0.05, and α2 = 0.2 as in [20]. The uniform random numbers appear in the intervals (0, 0.05), (0, 0.1), and (0, 0.2), respectively. Specially, [20] points out that the noise added to the blur image is independent and identically distributed, so the same set of parameters could be selected for different images when the same kind of noise is considered. This could reduce the computational cost of searching regularization parameters. For the speckle noise, as in [20, 23], let α ≡ 1, α1 = 2.5 × 104, the continuation scheme about α2 will be used, easier subsequences with smaller α2 can be solved quickly, and the later subproblems can also be solved relatively quickly with warm starts from previous solutions. The noise variance will be selected as 0, 0.01, and 0.05 respectively. To be fair, the same methods will be used to solve the primal alternating direction method in this paper.Now, we select two cases that have significant improvement of the visual sense to illustrate our method’s restoration quality. Under different conditions, more image restoration results with uniform noise and speckle noise will be summarized in Tables1 and 2, respectively. Figure 1(a) shows the original Camera image. Figure 1(b) shows the observed Camera image, where the noise is the uniform noise, the random numbers are in the interval (0, 0.2), and the support is equal to 9 × 9. Figure 1(c) shows the restored Camera image by the algorithm in [20]. Figure 1(d) shows the restored Camera image by our method. Figure 2(a) shows the original Lena image. Figure 2(b) shows the observed Lena image, where the noise is the speckle noise, the noise variance is 0.01, the support is equal to 7×7. Figure 2(c) shows the restored Lena image by the algorithm in [20]. Figure 2(d) shows the restored Lena image by our method.Table 1
Summary results for the uniform noise.
Image
Noise
Blur
SNR
PSNR
CPU
[20]
Ours
[20]
Ours
[20]
Ours
Camera
(0, 0.05)
5 × 5
21.25
21.46
26.83
26.84
5.49
4.84
7 × 7
19.35
19.37
24.93
24.92
6.82
5.99
9 × 9
18.28
18.30
23.86
23.88
8.47
7.99
(0, 0.1)
5 × 5
18.12
18.13
23.70
23.73
6.12
5.93
7 × 7
17.20
17.21
22.78
22.77
8.14
7.09
9 × 9
16.52
16.52
22.11
22.16
8.95
8.17
(0, 0.2)
5 × 5
12.86
12.86
18.45
18.57
7.57
7.07
7 × 7
12.90
12.91
18.48
18.48
7.82
7.00
9 × 9
12.81
12.88
18.39
18.46
9.14
8.03
Lena
(0, 0.05)
5 × 5
24.29
24.31
29.94
29.93
29.84
27.60
7 × 7
23.42
23.46
29.07
29.11
39.13
35.62
9 × 9
22.60
22.67
28.26
28.26
44.05
37.16
(0, 0.1)
5 × 5
19.37
19.37
25.03
25.07
30.50
24.35
7 × 7
19.18
19.15
24.84
24.84
38.24
33.65
9 × 9
18.90
18.92
24.55
24.57
45.47
40.36
(0, 0.2)
5 × 5
13.65
13.74
19.31
19.36
44.06
39.56
7 × 7
13.58
13.63
19.24
19.33
41.28
36.33
9 × 9
13.25
13.25
18.99
18.97
45.49
40.68Table 2
Summary results for the speckle noise.
Image
Noise
Blur
SNR
PSNR
CPU
[20]
Ours
[20]
Ours
[20]
Ours
Camera
0
5 × 5
22.06
22.12
27.64
27.64
5.09
4.36
7 × 7
19.94
19.93
25.52
25.57
7.55
6.05
9 × 9
18.72
18.74
23.89
23.90
8.13
6.88
0.01
5 × 5
20.32
20.33
25.90
25.93
5.88
5.02
7 × 7
19.22
19.27
24.83
24.84
8.30
6.91
9 × 9
18.28
18.28
24.30
24.35
8.13
7.04
0.05
5 × 5
15.77
15.79
21.35
21.35
9.48
8.25
7 × 7
14.95
14.96
20.54
20.55
9.16
7.25
9 × 9
12.63
12.62
18.21
18.24
9.87
8.23
Lena
0
5 × 5
28.33
28.34
33.99
34.01
29.94
27.61
7 × 7
26.07
26.10
31.72
31.72
37.44
32.56
9 × 9
24.59
24.62
30.25
30.61
45.29
38.65
0.01
5 × 5
23.64
23.69
29.30
29.32
36.86
32.58
7 × 7
23.86
23.94
28.92
29.00
42.32
36.32
9 × 9
23.24
23.24
28.89
28.92
50.97
43.65
0.05
5 × 5
18.23
18.25
23.89
23.90
61.20
55.63
7 × 7
16.58
16.57
22.24
22.31
61.31
56.86
9 × 9
16.48
16.52
20.14
20.14
71.48
63.53Figure 1
The restored Camera images by different methods with the uniform noise. The random numbers are in the interval (0, 0.2), the support being equal to9×9. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [20]. (d) Image restored by our method.
(a)
(b)
(c)
(d)Figure 2
The restored Lena images by different methods with the speckle noise. The noise variance is 0.01, the support being equal to7×7. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [20]. (d) Image restored by our method.
(a)
(b)
(c)
(d)From Figures1 and 2, we see that the quality of images restored by our method is improved obviously compared to those restored by the method of [20] in the two cases. For the uniform noise and speckle noise, Tables 1 and 2 show that the SNR and PSNR of the images restored by our method both are higher than the results of [20] in most cases. For the same requirement of relative error, our proposed algorithm is clearly faster than the competing algorithm in [20]. Referring to Sections 2 and 3, the proposed algorithm uses subspace optimization method to correct the search direction, which could ensure F(l¯i+1)≤F(li+1) is valid on the ith iteration. So it is clear that the proposed method is more efficient. The analysis is verified by the experiments results above.


## 5.2. Comparison Experiment for thel2-TV Model

In the subsection, the algorithm in [18] and ours will be compared. Being different from [19, 20], [18] involves a fitting of the auxiliary variable to Df which has superior performance compared to fitting f [18, 23, 32]. The additive noise is selected as Gaussian noise. Here four sets of the standard deviation of noise are chosen: 0.01, 0.05, 0.1, and 0.2. Same methods appearing in [18] will be used to solve the primal alternating direction method.Similarly, we will first select two cases that have significant improvement of the visual sense to illustrate our method’s restoration quality. Figures3(a) and 4(a) show the original images. Figures 3(b) and 4(b) show the observed images. Figure 3(c) shows the restored Camera image by algorithm in [18]; Figure 3(d) shows the restored Camera image by our method, where the standard deviation of noise is 0.1 and the support is equal to 9 × 9. Figure 4(c) shows the restored Lena image by algorithm in [18]; Figure 4(d) shows the restored Lena image by our method, where the standard deviation of noise is 0.2 and the support is equal to 9×9. Table 3 shows more image restoration results with the different support and standard deviation of noise.Table 3
Summary results for Gaussian noise.
Image
Noise
Blur
SNR
PSNR
CPU
[18]
Ours
[18]
Ours
[18]
Ours
Camera
0.01
5 × 5
21.99
21.98
27.57
27.61
2.45
2.31
7 × 7
19.90
19.93
25.49
25.49
2.67
2.32
9 × 9
18.69
18.69
24.27
24.28
2.73
2.63
0.05
5 × 5
20.64
20.69
26.22
26.24
1.83
1.62
7 × 7
19.04
19.12
24.62
24.66
7.60
7.02
9 × 9
18.01
18.07
23.59
23.59
2.57
2.26
0.1
5 × 5
17.97
17.99
23.55
23.58
1.84
1.56
7 × 7
17.06
17.11
22.64
22.63
2.20
2.12
9 × 9
16.41
16.40
22.00
22.61
2.22
1.98
0.2
5 × 5
13.37
17.38
18.95
18.94
1.78
1.51
7 × 7
13.11
13.10
18.70
18.73
1.68
1.53
9 × 9
12.91
12.93
18.49
18.61
1.86
1.68
Lena
0.01
5 × 5
28.05
28.10
33.71
33.73
11.89
10.32
7 × 7
25.90
25.91
31.56
31.62
12.93
10.33
9 × 9
24.47
24.49
30.13
30.15
13.32
10.92
0.05
5 × 5
24.09
24.06
29.74
29.77
10.39
9.66
7 × 7
23.12
23.12
28.78
28.79
11.26
9.65
9 × 9
22.33
22.42
27.99
28.12
13.18
11.56
0.1
5 × 5
19.50
19.55
25.16
25.17
9.14
8.68
7 × 7
19.19
19.23
24.85
24.84
10.26
9.33
9 × 9
18.86
18.88
24.52
24.51
11.79
9.91
0.2
5 × 5
13.81
13.81
19.47
19.44
7.14
6.56
7 × 7
13.83
13.86
19.49
19.83
8.47
7.82
9 × 9
13.78
14.01
19.44
19.87
9.78
7.97Figure 3
The restored Camera images by different methods with Gaussian noise. The standard deviation of noise is 0.1; the support is9×9. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [18]. (d) Image restored by our method.
(a)
(b)
(c)
(d)Figure 4
The restored Lena images by different methods with Gaussian noise. The standard deviation of noise is 0.2; the support is 9 × 9. (a) Original image. (b) Observed image. (c) Image restored by algorithm in [18]. (d) Image restored by our method.
(a)
(b)
(c)
(d)As in Section5.1, Figures 3 and 4 show that the quality of restoration images by using our method is better than those restored by the algorithm of [18] in the two cases. We also can know from Table 3 that the SNR and PSNR of the images restored by our method both are higher than the results of [18] in most cases. With the same reasons mentioned in l1-TV model, all computational time required by our method is significantly less than that required by the algorithm in [18]. So the performance of the algorithm in [18] also has been improved through the correction.


## 6. Concluding Remarks

Taking thel1-TV model in image restoration as an example, a modified method is proposed to overcome drawbacks of the primal alternating direction method in this paper. We have illustrated that our method about how to correct the search direction improves the optimization performance. In addition, the convergence of the primal alternating direction method has been proven under some weaker conditions, and thus the convergence of proposed method is easily obtained by the equivalence between them. The experimental results based on two models in [18, 20] show that the proposed method could enhance the quality of restored images in most cases and the efficiency of algorithms has been significantly improved. In fact, our method can be applied to many other cases optimized alternatively.



---

*Source: 6921742-2018-01-21.xml*