# Fast Total-Variation Image Deconvolution with Adaptive Parameter Estimation via Split Bregman Method

Authors: Chuan He; Changhua Hu; Wei Zhang; Biao Shi; Xiaoxiang Hu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The total-variation (TV) regularization has been widely used in image restoration domain, due to its attractive edge preservation ability. However, the estimation of the regularization parameter, which balances the TV regularization term and the data-fidelity term, is a difficult problem. In this paper, based on the classical split Bregman method, a new fast algorithm is derived to simultaneously estimate the regularization parameter and to restore the blurred image. In each iteration, the regularization parameter is refreshed conveniently in a closed form according to Morozov’s discrepancy principle. Numerical experiments in image deconvolution show that the proposed algorithm outperforms some state-of-the-art methods both in accuracy and in speed.


---

## Body

## 1. Introduction

Digital image restoration, which aims at recovering an estimate of the original scene from the degraded observation, is a recurrent task with many real-world applications, for example, remote sensing, astronomy, and medical imaging. During acquisition, the observed images are often degraded by relative motion between the camera and the original scene, defocusing of the lens system, atmospheric turbulence, and so forth. In most cases, the degradation can be modeled as a spatially linear shift invariant system, where the original image is convolved by a spatially invariant point spread function (PSF) and contaminated with Gaussian white noise [1].Without loss of generality, we assume that the digital gray-scale images used throughout this paper have anm
×
n domain and are represented by m
n vectors formed by stacking up the image matrix rows. So the (
i
,
j
)th pixel becomes the (
(
i
-
1
)
n
+
j
)th entry of the vector. Then, in general, the degradation process can be modeled as the following discrete linear inverse problem:
(1)
f
=
H
u
clean
+
n
,
where f and u
clean are the observed image and the original image, respectively, both expressed in vectorial form, H is the convolution operator in accordance with the spatially invariant PSF, which is assumed to be known, and n is a vector of zero mean Gaussian white noise of variance σ
2. In most cases, H is ill-conditioned so that directly estimating u
clean from f is of no possibility. The solution of (1) is highly sensitive to noise in the observed image and it becomes a well-known ill-posed linear inverse problem (IPLIP). The inverse filtering in a least square form, which tries to solve this problem directly, usually results in an estimation of no usability.If we get some prior knowledge such as prior distribution or sparse quality about the original image, we can incorporate such information into the restoration process via some sort of regularization [2]. This makes the solution of IPLIP possible. A large class of regularization approaches leads to the following minimization problem:
(2)
min
⁡
u
{
Φ
(
u
)
+
λ
2
∥
H
u
-
f
∥
2
2
}
,
where u is the estimate of u
clean and λ is the so-called regularization parameter. The first term of (2) represents the regularization term, whereas the second represents the data-fidelity term. The regularization has the quality of numerical stabilizing and encourages the result to have some desirable properties. The positive regularization parameter λ plays the role of balancing the relative weight of the two terms.Among the various regularization methods, the total-variation (TV) regularization is famed for its attractive edge preservation ability. It was introduced into image restoration by Rudin et al. [3] in 1992. From then on, the TV regularization has been arousing significant attention [4–7], and, so far, it has resulted in several variants [8–10]. The objective functional of the TV restoration problem is given by
(3)
min
⁡
u
{
∑
i
∥
D
i
u
∥
2
+
λ
2
∥
H
u
-
f
∥
2
2
}
,
where the first term is the so-called TV seminorm of u and D
i
u (its detailed definition is in Section 2) is the discrete gradient of u at pixel i. In minimization functional (3), the TV is either isotropic if ∥
·
∥ is 2-norm or anisotropic if it is 1-norm. We emphasize here that our method is applicable to both isotropic and anisotropic cases. However, we will only treat the isotropic one for simplicity, since the treatment for the other one is completely analogous. Despite the advantage of edge preservation, the minimization of functional (3) is troublesome and it has no closed form solution at all. Various methods have been proposed to minimize (3), including time-marching schemes [3], primal-dual based methods [11–13], fixed point iteration approaches [14], and variable splitting algorithms [15–17]. In particular, the split Bregman method adopted in this paper is an instance of the variable splitting based algorithms.Another critical issue in TV regularization is the selection of the regularization parameterλ, since it plays a very important role. If λ is too large, the regularized solution will be undersmoothed, and, on the contrary, if λ is too small, the regularized solution will not fit the observation properly. Most works in the literature only consider a fixed λ and, when applying these methods to image restoration problems, one should adjust λ manually to get a satisfying solution. So far, a few strategies are proposed for the adaptive estimation of parameter λ, for example, the L-curve method [18], the variational Bayesian approach [19], the generalized cross-validation (GCV) method [20], and Morozov’s discrepancy principle [21].If the noise level is available or can be estimated first, Morozov’s discrepancy principle is a good choice for the selection ofλ. According to this rule, the TV image restoration problem can be described as
(4)
min
⁡
u
∑
i
∥
D
i
u
∥
2
s
.
t
.
u
∈
S
,
where S
:
=
{
u
:
∥
H
u
-
f
∥
2
2
≤
c
} with c
=
τ
m
n
σ
2 is the feasible set in accordance with the discrepancy principle. Although it is much easier to solve the unconstrained problem (3) than the constrained problem (4), formulation (4) has a clear physics meaning (c is proportional to the noise variance) and this makes the estimation of λ easier. In fact, referring to the theory of Lagrangian methods, if u is a solution of constrained problem (4), it will also be a solution of (3) for a particular choice of λ
≥
0, which is the Lagrangian multiplier corresponding to the constraint in (4). To minimize (4), we have either u
∈
S for λ
=
0 or
(5)
∥
H
u
-
f
∥
2
2
=
c
for λ
>
0. In fact, if λ
=
0, minimizing (3) is equivalent to minimizing ∑
i
∥
D
i
u
∥
2, which means that the solution is a constant image. Obviously, this will not happen to a nature image. Therefore, only λ
>
0 will happen in practical applications.There exists no closed form solution of functional (3) or (4), and, up to now, several papers pay attention to the numerical solving of problem (4). In [22], the authors provided a modular solver to update λ for making use of existing methods for the unconstrained problems. Afonso et al. [17] proposed an alternating direction method of multipliers (ADMM) based approach and suggested using Chambolle’s dual method [23] to adaptively restore the degraded image. In [13], Wen and Chan proposed a primal-dual based method to solve the constrained problem (4). The minimization problem was transformed into a saddle point problem of the primal-dual model of (4), and then the proximal point method [24] was applied to find the saddle point. When dealing with the updating of λ, they resorted to a Newton’s inner iteration. All these methods mentioned above have the same limitation: in order to adaptively update λ, an inner iteration is introduced, and this results in extra computing cost.In this paper, based on the split Bregman scheme, we propose a fast algorithm to solve the constrained TV restoration problem (4). When referring to the variable splitting technique, we introduce two auxiliary variables to represent D
u and the TV norm, respectively, and therefore the constrained problem (4) can be solved efficiently with a separable structure without any inner iteration. Differing from the previous works focusing on the adaptive regularization parameter estimation in TV restoration problems, our method involves no inner iteration and adjusts the regularization parameter in a closed form in each iteration. Thus, fast computation speed is achieved. The simulation results in TV restoration problems indicate that our method outperforms some famous methods in accuracy and especially in speed. According to the equivalence of split Bregman method, ADMM, and Douglas-Rachford splitting algorithm under the assumption of linear constraints [25–27], our algorithm can also be seen as an instance of ADMM or Douglas-Rachford splitting algorithm.In the rest of this paper, the basic notation is presented in Section2. Section 3 gives the derivation leading to the proposed algorithm and some practical parameter setting strategies. In Section 4, several experiments are reported to demonstrate the effectiveness of our algorithm. Finally, Section 5 draws a short conclusion of this paper.


## 2. Basic Notation

Let us describe the notation that we will be using throughout this paper. Euclidean spaceR
m
n is denoted as P, whereas Euclidean space R
m
n
×
m
n is denoted as T
:
=
P
×
P. The ith components of x
∈
P and y
∈
T are denoted as x
i
∈
R and y
i
=
(
y
i
(
1
)
,
y
i
(
2
)
)
T
∈
R
2, respectively. Define inner products 〈
x
,
x
〉
P
=
∑
i
x
i
x
i, 〈
y
,
y
〉
T
=
∑
i
∑
k
=
1
2
y
i
(
k
)
y
i
(
k
), and norm ∥
x
∥
2
=
〈
x
,
x
〉
P, ∥
y
∥
2
=
〈
y
,
y
〉
T. For each u
∈
P, we define D
i
u
:
=
[
(
D
(
1
)
u
)
i
,
(
D
(
2
)
u
)
i
]
T, with
(6)
(
D
(
1
)
u
)
i
:
=
{
u
i
+
n
-
u
i
,
if
1
≤
i
≤
n
(
m
-
1
)
,
u
mod
⁡
(
i
,
n
)
-
u
i
.
otherwise
,
(7)
(
D
(
2
)
u
)
i
:
=
{
u
i
+
1
-
u
i
,
if
mod
⁡
(
i
,
n
)
≠
0
,
u
i
-
n
+
1
-
u
i
.
otherwise
,
where D
(
1
)
,
D
(
2
)
∈
R
m
n
×
m
n are m
n
×
m
n matrices in the vertical and horizontal directions, and obviously it holds that D
(
1
)
u
∈
P
and
D
(
2
)
u
∈
P. D
i
∈
R
2
×
m
n is a tow-row matrix formed by stacking the ith rows of D
(
1
) and D
(
2
) together. Define the global first-order finite difference operator as D
:
=
[
(
D
(
1
)
)
T
,
(
D
(
2
)
)
T
]
T
∈
R
2
m
n
×
m
n and we consider D
u
∈
T. From (6) and (7), we see that the periodic boundary condition is assumed here.Given a convex functionalJ
(
z
), the subdifferential ∂
J
(
z
1
) of J
(
z
) at z
1 is defined as
(8)
∂
J
(
z
1
)
:
=
{
q
∈
P
:
〈
q
,
z
-
z
1
〉
≤
J
(
z
)
-
J
(
z
1
)
,
∀
z
∈
P
}
.
And the Bregman distance between z and z
1 is defined as
(9)
D
J
(
z
1
)
=
J
(
z
)
-
J
(
z
1
)
-
〈
q
,
z
-
z
1
〉
.
From the definition of Bregman distance, we learn that it is positive all the time.


## 3. Methodology

### 3.1. Deduction of the Proposed Algorithm

We refer to the variable splitting technique [28] for liberating the discrete operator D
i
u out from nondifferentiability and simplifying the regularization parameter’s updating. An auxiliary variable x
∈
P is introduced for H
u, and another auxiliary variable y
∈
T is introduced to represent D
u (or y
i
∈
R
2 for D
i
u, resp.). Therefore, functional (3) is equivalent to
(10)
min
⁡
u
,
x
,
y
{
∑
i
∥
y
i
∥
2
+
λ
2
∥
x
-
f
∥
2
2
}
subject
to
H
u
=
x
,
y
i
=
D
i
u
,
i
=
1,2
,
…
,
m
n
.
Define Bregman functional
(11)
J
(
u
,
x
,
y
)
=
{
∑
i
∥
y
i
∥
2
+
λ
2
∥
x
-
f
∥
2
2
}
.
Then the Bregman distance of J
(
u
,
x
,
y
) is
(12)
D
J
(
p
u
k
,
p
x
k
,
p
y
k
)
(
u
,
x
,
y
;
u
k
,
x
k
,
y
k
)
=
J
(
u
,
x
,
y
)
-
J
(
u
k
,
x
k
,
y
k
)
-
〈
p
u
k
,
u
-
u
k
〉

-

〈
p
x
k
,
x
-
x
k
〉
-
〈
p
y
k
,
y
-
y
k
〉
.
According to the split Bregman method [16, 29], we obtain the following iterative scheme:
(13)
(
u
k
+
1
,
x
k
+
1
,
y
k
+
1
)
=
argmin
u
,
x
,
y
{
D
J
(
p
u
k
,
p
x
k
,
p
y
k
)
(
u
,
x
,
y
;
u
k
,
x
k
,
y
k
)
+
β
1
2
∥
x
-
H
u
∥
2
2
+
β
2
2
∥
y
-
D
u
∥
2
2
}
,
(14)
p
u
k
+
1
=
p
u
k
+
β
1
H
T
(
x
k
+
1
-
H
u
k
+
1
)
+
β
2
D
T
(
y
k
+
1
-
D
u
k
+
1
)
,
(15)
p
x
k
+
1
=
p
x
k
+
β
1
(
H
u
k
+
1
-
x
k
+
1
)
,
(16)
p
y
k
+
1
=
p
y
k
+
β
2
(
D
u
k
+
1
-
y
k
+
1
)
,
if we define that
(17)
p
u
0
:
=
-
β
1
H
T
b
0
-
β
2
D
T
d
0
p
x
0
:
=
β
1
b
0
p
y
0
:
=
β
2
d
0
,
for any elements b
0
∈
P and  d
0
∈
T, and then, according to (14)–(16), it holds that
(18)
p
u
k
=
-
β
1
H
T
b
k
-
β
2
D
T
d
k
p
x
k
=
β
1
b
k
p
y
k
=
β
2
d
k
k
=
0,1
,
…
and we obtain the following iterative scheme:
(19)
(
u
k
+
1
,
x
k
+
1
,
y
k
+
1
)
=
arg
min
⁡
u
,
x
,
y
{
+
∑
i
∥
y
i
∥
2
+
β
2
2
∥
y
-
D
u
-
d
k
∥
2
2
λ
2
∥
x
-
f
∥
2
2
+
β
1
2
∥
x
-
H
u
-
b
k
∥
2
2
+
∑
i
∥
y
i
∥
2
+
β
2
2
∥
y
-
D
u
-
d
k
∥
2
2
}
,
b
k
+
1
=
b
k
+
H
u
k
+
1
-
x
k
+
1
,
d
k
+
1
=
d
k
+
D
u
k
+
1
-
y
k
+
1
.
In iterative scheme (19), the problem yielding (
u
k
+
1
,
x
k
+
1
,
y
k
+
1
) exactly is difficult, since it needs an inner iterative scheme. Here, we adopt the alternating direction method (ADM) to approximately calculate u
k
+
1, x
k
+
1, and y
k
+
1 in each iteration and this leads to the following iterative framework:
(20)
u
k
+
1
=
arg
min
⁡
u
{
β
1
2
∥
x
k
-
H
u
-
b
k
∥
2
2
+
β
2
2
∥
y
k
-
D
u
-
d
k
∥
2
2
}
,
(21)
y
k
+
1
=
arg
min
⁡
y
{
∑
i
∥
y
i
∥
2
+
β
2
2
∥
y
-
D
u
-
d
k
∥
2
2
}
,
(22)
x
k
+
1
=
arg
min
⁡
x
{
λ
k
+
1
2
∥
x
-
f
∥
2
2
+
β
1
2
∥
x
-
H
u
k
+
1
-
b
k
∥
2
2
}
,
(23)
b
k
+
1
=
b
k
+
H
u
k
+
1
-
x
k
+
1
,
(24)
d
k
+
1
=
d
k
+
D
u
k
+
1
-
y
k
+
1
.
In the following, we will discuss how to solve problems (20)–(22) efficiently.The minimization subproblem with respect tou is in the form of least square. From functional (20), we obtain
(25)
(
β
1
β
2
H
T
H
+
D
T
D
)
u
=
β
1
β
2
H
T
(
x
k
-
b
k
)
+
D
T
(
y
k
-
d
k
)
.Under the periodic boundary condition, matricesH, D
(
1
), and D
(
2
) are block-circulant, so they can be diagonalized by a Discrete Fourier Transforms (DFTs) matrix. Using the convolution theorem of Fourier Transforms, we obtain
(26)
u
k
+
1
=
ℱ
-
1
(
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
+
ℱ
*
(
D
(
1
)
)
ℱ
(
(
y
k
)
(
1
)
-
(
d
k
)
(
1
)
)
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
)
∘
(
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
H
)
+
ℱ
*
(
D
(
1
)
)
∘
ℱ
(
D
(
1
)
)
+
ℱ
*
(
D
(
2
)
)
∘
ℱ
(
D
(
2
)
)
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
)
-
1
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
/
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
)
,
where ℱ denotes the DFT, “*” denotes complex conjugate, and “∘” represents componentwise multiplication. The reciprocal notation is also componentwise here. Therefore, problem (20) can be solved by two Fast Fourier Transforms (FFTs) and one inverse FFT in O
(
m
n
log
⁡
(
m
n
)
) operations.Functional (21) is a proximal minimization problem and it can be solved componentwise by a two-dimension shrinkage as follows:
(27)
y
i
k
+
1
=
max
⁡
{
∥
D
i
u
k
+
1
+
d
i
k
∥
2
-
1
β
2
,
0
}
D
i
u
k
+
1
+
d
i
k
∥
D
i
u
k
+
1
+
d
i
k
∥
2
.
During the calculation, we employ the convention 0 × (0/0) = 0 to avoid getting results of no meaning.When dealing with problem (22), we assume that w
k
+
1
=
H
u
k
+
1
+
b
k first. It is obvious that x is λ related and it plays the role of H
u. Therefore, in each iteration, we should examine whether ∥
x
-
f
∥
2
2
≤
c holds true, that is, whether x meets the discrepancy principle.The solutions ofλ and x fall into two cases according to the range of w
k
+
1.(1)
If(28)
∥
w
k
+
1
-
f
∥
2
2
≤
c
holds true, we set λ
k
+
1
=
0 and x
k
+
1
=
w
k
+
1. Obviously this x
k
+
1 satisfies the discrepancy principle.
(2)
If∥
w
k
+
1
-
f
∥
2
2
>
c, according to the discrepancy principle, we should solve equation
(29)
∥
w
k
+
1
-
f
∥
2
2
=
c
.Since the minimization problem (22) with respect to x is quadratic, it has a closed form solution
(30)
x
k
+
1
=
(
λ
k
+
1
f
+
β
1
w
k
+
1
)
(
λ
k
+
1
+
β
1
)
.
Substituting x
k
+
1 in (29) with (30), we obtain
(31)
λ
k
+
1
=
β
1
∥
f
-
w
k
+
1
∥
2
c
-
β
1
.The above discussion can be summed up by Algorithm1.Algorithm 1:APE-SBA: Adaptive Parameter Estimation Split Bregman Algorithm.
Input: f, H, c.
(1)Initialize u
0, x
0, b
0, d
0. Set k
=
0 and β
1
>
0 and β
2
>
0
(2)
while stopping criterion is not satisfied, do
(3)Compute u
k
+
1 according to (26);
(4)Compute y
k
+
1 according to (27);
(5)
if (28) holds, then
(6)
λ
k
+
1
=
0, and x
k
+
1
=
w
k
+
1;
(7)
else
(8)Update λ
k
+
1 and x
k
+
1 according to (31) and (30);
(9)
end if
(10)
Update b
k
+
1 and d
k
+
1 according to (23) and (24);
(11)
k
=
k
+
1;
(12)
end while
(13)
return  λ
k
+
1 and u
k
+
1.In algorithm APE-SBA, by introducing the auxiliary variablex, H
u is liberated out from the constraint of the discrepancy principle, and consequently a closed form to update λ is obtained without any inner iteration. This is the major difference between APE-SBA and the methods in [13] and [17]. Since the procedure of solving (26) corresponding to the u subproblem consumes the most, the calculation cost of our algorithm is about O
(
m
n
log
⁡
(
m
n
)
) FFT operations. In fact, our algorithm is an instance of the classical split Bregman method, so the convergence of it is guaranteed by the theorem proposed by Eckstein and Bertsekas [30]. We summarize the convergence of our algorithm as follows.Theorem 1.
Forβ
1
,
β
2
>
0, the sequence {
u
k
,
x
k
,
y
k
,
b
k
,
d
k
,
λ
k
} generated by Algorithm APE-SBA from any initial point (
u
0
,
x
0
,
b
0
,
d
0
) converges to (
u
*
,
x
*
,
y
*
,
b
*
,
d
*
,
λ
*
), where (
u
*
,
x
*
,
y
*
) is a solution of the functional (10). In particular, u
* is the minimizer of functional (4), and λ
* is the Lagrange multiplier corresponding to constraint u
∈
S according to the unconstrained problem (3).


### 3.2. Parameter Setting

In this paper, the noise level is denoted by the following defined blurred signal-to-noise ratio (BSNR)(32)
BSNR
=
10
lo
g
10
(
∥
f
-
f
-
∥
2
2
m
n
σ
2
)
,
where f
- denotes the mean of f.In minimization problem (4), the noise dependent upper bound c is very important, since a good choice of it can constrain the error between the restored image and the original image to a reasonable level. To our best knowledge, the choice of this parameter is an open problem which has not been solved theoretically. One approach to choose c is referring to the equivalent degrees of freedom (DF), but the calculation of DF is a difficult problem and we can only get an estimate of it. A simple strategy of choosing c is to employ a curve approximating the relation between the noise level and τ. By fitting experimental data with a straight line, in this paper, we suggest setting
(33)
τ
=
-
0.006
×
BSNR
+
1.09
.
Besides the setting of τ, the choice of β
1 and β
2 is essential to our algorithm. We suggest setting β
1 = 10(BSNR/10−1) × β
2, where β
2 = 1. This parameter setting is obtained by large numbers of experiments. Actually, β
1, β
2 > 0 is sufficient for the convergence of the proposed algorithm, but why β
1 and β
2 play different important role when the BSNR varies? The reason is that, when the BSNR becomes higher, the distance between H
u and f is nearer. From minimization problem (10), we learn that auxiliary variable x plays the role of H
u and a higher BSNR means a larger β
1.


## 3.1. Deduction of the Proposed Algorithm

We refer to the variable splitting technique [28] for liberating the discrete operator D
i
u out from nondifferentiability and simplifying the regularization parameter’s updating. An auxiliary variable x
∈
P is introduced for H
u, and another auxiliary variable y
∈
T is introduced to represent D
u (or y
i
∈
R
2 for D
i
u, resp.). Therefore, functional (3) is equivalent to
(10)
min
⁡
u
,
x
,
y
{
∑
i
∥
y
i
∥
2
+
λ
2
∥
x
-
f
∥
2
2
}
subject
to
H
u
=
x
,
y
i
=
D
i
u
,
i
=
1,2
,
…
,
m
n
.
Define Bregman functional
(11)
J
(
u
,
x
,
y
)
=
{
∑
i
∥
y
i
∥
2
+
λ
2
∥
x
-
f
∥
2
2
}
.
Then the Bregman distance of J
(
u
,
x
,
y
) is
(12)
D
J
(
p
u
k
,
p
x
k
,
p
y
k
)
(
u
,
x
,
y
;
u
k
,
x
k
,
y
k
)
=
J
(
u
,
x
,
y
)
-
J
(
u
k
,
x
k
,
y
k
)
-
〈
p
u
k
,
u
-
u
k
〉

-

〈
p
x
k
,
x
-
x
k
〉
-
〈
p
y
k
,
y
-
y
k
〉
.
According to the split Bregman method [16, 29], we obtain the following iterative scheme:
(13)
(
u
k
+
1
,
x
k
+
1
,
y
k
+
1
)
=
argmin
u
,
x
,
y
{
D
J
(
p
u
k
,
p
x
k
,
p
y
k
)
(
u
,
x
,
y
;
u
k
,
x
k
,
y
k
)
+
β
1
2
∥
x
-
H
u
∥
2
2
+
β
2
2
∥
y
-
D
u
∥
2
2
}
,
(14)
p
u
k
+
1
=
p
u
k
+
β
1
H
T
(
x
k
+
1
-
H
u
k
+
1
)
+
β
2
D
T
(
y
k
+
1
-
D
u
k
+
1
)
,
(15)
p
x
k
+
1
=
p
x
k
+
β
1
(
H
u
k
+
1
-
x
k
+
1
)
,
(16)
p
y
k
+
1
=
p
y
k
+
β
2
(
D
u
k
+
1
-
y
k
+
1
)
,
if we define that
(17)
p
u
0
:
=
-
β
1
H
T
b
0
-
β
2
D
T
d
0
p
x
0
:
=
β
1
b
0
p
y
0
:
=
β
2
d
0
,
for any elements b
0
∈
P and  d
0
∈
T, and then, according to (14)–(16), it holds that
(18)
p
u
k
=
-
β
1
H
T
b
k
-
β
2
D
T
d
k
p
x
k
=
β
1
b
k
p
y
k
=
β
2
d
k
k
=
0,1
,
…
and we obtain the following iterative scheme:
(19)
(
u
k
+
1
,
x
k
+
1
,
y
k
+
1
)
=
arg
min
⁡
u
,
x
,
y
{
+
∑
i
∥
y
i
∥
2
+
β
2
2
∥
y
-
D
u
-
d
k
∥
2
2
λ
2
∥
x
-
f
∥
2
2
+
β
1
2
∥
x
-
H
u
-
b
k
∥
2
2
+
∑
i
∥
y
i
∥
2
+
β
2
2
∥
y
-
D
u
-
d
k
∥
2
2
}
,
b
k
+
1
=
b
k
+
H
u
k
+
1
-
x
k
+
1
,
d
k
+
1
=
d
k
+
D
u
k
+
1
-
y
k
+
1
.
In iterative scheme (19), the problem yielding (
u
k
+
1
,
x
k
+
1
,
y
k
+
1
) exactly is difficult, since it needs an inner iterative scheme. Here, we adopt the alternating direction method (ADM) to approximately calculate u
k
+
1, x
k
+
1, and y
k
+
1 in each iteration and this leads to the following iterative framework:
(20)
u
k
+
1
=
arg
min
⁡
u
{
β
1
2
∥
x
k
-
H
u
-
b
k
∥
2
2
+
β
2
2
∥
y
k
-
D
u
-
d
k
∥
2
2
}
,
(21)
y
k
+
1
=
arg
min
⁡
y
{
∑
i
∥
y
i
∥
2
+
β
2
2
∥
y
-
D
u
-
d
k
∥
2
2
}
,
(22)
x
k
+
1
=
arg
min
⁡
x
{
λ
k
+
1
2
∥
x
-
f
∥
2
2
+
β
1
2
∥
x
-
H
u
k
+
1
-
b
k
∥
2
2
}
,
(23)
b
k
+
1
=
b
k
+
H
u
k
+
1
-
x
k
+
1
,
(24)
d
k
+
1
=
d
k
+
D
u
k
+
1
-
y
k
+
1
.
In the following, we will discuss how to solve problems (20)–(22) efficiently.The minimization subproblem with respect tou is in the form of least square. From functional (20), we obtain
(25)
(
β
1
β
2
H
T
H
+
D
T
D
)
u
=
β
1
β
2
H
T
(
x
k
-
b
k
)
+
D
T
(
y
k
-
d
k
)
.Under the periodic boundary condition, matricesH, D
(
1
), and D
(
2
) are block-circulant, so they can be diagonalized by a Discrete Fourier Transforms (DFTs) matrix. Using the convolution theorem of Fourier Transforms, we obtain
(26)
u
k
+
1
=
ℱ
-
1
(
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
+
ℱ
*
(
D
(
1
)
)
ℱ
(
(
y
k
)
(
1
)
-
(
d
k
)
(
1
)
)
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
)
∘
(
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
H
)
+
ℱ
*
(
D
(
1
)
)
∘
ℱ
(
D
(
1
)
)
+
ℱ
*
(
D
(
2
)
)
∘
ℱ
(
D
(
2
)
)
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
)
-
1
(
+
ℱ
*
(
D
(
2
)
)
ℱ
(
(
y
k
)
(
2
)
-
(
d
k
)
(
2
)
)
(
β
1
/
β
2
)
ℱ
*
(
H
)
∘
ℱ
(
x
k
-
b
k
)
)
,
where ℱ denotes the DFT, “*” denotes complex conjugate, and “∘” represents componentwise multiplication. The reciprocal notation is also componentwise here. Therefore, problem (20) can be solved by two Fast Fourier Transforms (FFTs) and one inverse FFT in O
(
m
n
log
⁡
(
m
n
)
) operations.Functional (21) is a proximal minimization problem and it can be solved componentwise by a two-dimension shrinkage as follows:
(27)
y
i
k
+
1
=
max
⁡
{
∥
D
i
u
k
+
1
+
d
i
k
∥
2
-
1
β
2
,
0
}
D
i
u
k
+
1
+
d
i
k
∥
D
i
u
k
+
1
+
d
i
k
∥
2
.
During the calculation, we employ the convention 0 × (0/0) = 0 to avoid getting results of no meaning.When dealing with problem (22), we assume that w
k
+
1
=
H
u
k
+
1
+
b
k first. It is obvious that x is λ related and it plays the role of H
u. Therefore, in each iteration, we should examine whether ∥
x
-
f
∥
2
2
≤
c holds true, that is, whether x meets the discrepancy principle.The solutions ofλ and x fall into two cases according to the range of w
k
+
1.(1)
If(28)
∥
w
k
+
1
-
f
∥
2
2
≤
c
holds true, we set λ
k
+
1
=
0 and x
k
+
1
=
w
k
+
1. Obviously this x
k
+
1 satisfies the discrepancy principle.
(2)
If∥
w
k
+
1
-
f
∥
2
2
>
c, according to the discrepancy principle, we should solve equation
(29)
∥
w
k
+
1
-
f
∥
2
2
=
c
.Since the minimization problem (22) with respect to x is quadratic, it has a closed form solution
(30)
x
k
+
1
=
(
λ
k
+
1
f
+
β
1
w
k
+
1
)
(
λ
k
+
1
+
β
1
)
.
Substituting x
k
+
1 in (29) with (30), we obtain
(31)
λ
k
+
1
=
β
1
∥
f
-
w
k
+
1
∥
2
c
-
β
1
.The above discussion can be summed up by Algorithm1.Algorithm 1:APE-SBA: Adaptive Parameter Estimation Split Bregman Algorithm.
Input: f, H, c.
(1)Initialize u
0, x
0, b
0, d
0. Set k
=
0 and β
1
>
0 and β
2
>
0
(2)
while stopping criterion is not satisfied, do
(3)Compute u
k
+
1 according to (26);
(4)Compute y
k
+
1 according to (27);
(5)
if (28) holds, then
(6)
λ
k
+
1
=
0, and x
k
+
1
=
w
k
+
1;
(7)
else
(8)Update λ
k
+
1 and x
k
+
1 according to (31) and (30);
(9)
end if
(10)
Update b
k
+
1 and d
k
+
1 according to (23) and (24);
(11)
k
=
k
+
1;
(12)
end while
(13)
return  λ
k
+
1 and u
k
+
1.In algorithm APE-SBA, by introducing the auxiliary variablex, H
u is liberated out from the constraint of the discrepancy principle, and consequently a closed form to update λ is obtained without any inner iteration. This is the major difference between APE-SBA and the methods in [13] and [17]. Since the procedure of solving (26) corresponding to the u subproblem consumes the most, the calculation cost of our algorithm is about O
(
m
n
log
⁡
(
m
n
)
) FFT operations. In fact, our algorithm is an instance of the classical split Bregman method, so the convergence of it is guaranteed by the theorem proposed by Eckstein and Bertsekas [30]. We summarize the convergence of our algorithm as follows.Theorem 1.
Forβ
1
,
β
2
>
0, the sequence {
u
k
,
x
k
,
y
k
,
b
k
,
d
k
,
λ
k
} generated by Algorithm APE-SBA from any initial point (
u
0
,
x
0
,
b
0
,
d
0
) converges to (
u
*
,
x
*
,
y
*
,
b
*
,
d
*
,
λ
*
), where (
u
*
,
x
*
,
y
*
) is a solution of the functional (10). In particular, u
* is the minimizer of functional (4), and λ
* is the Lagrange multiplier corresponding to constraint u
∈
S according to the unconstrained problem (3).


## 3.2. Parameter Setting

In this paper, the noise level is denoted by the following defined blurred signal-to-noise ratio (BSNR)(32)
BSNR
=
10
lo
g
10
(
∥
f
-
f
-
∥
2
2
m
n
σ
2
)
,
where f
- denotes the mean of f.In minimization problem (4), the noise dependent upper bound c is very important, since a good choice of it can constrain the error between the restored image and the original image to a reasonable level. To our best knowledge, the choice of this parameter is an open problem which has not been solved theoretically. One approach to choose c is referring to the equivalent degrees of freedom (DF), but the calculation of DF is a difficult problem and we can only get an estimate of it. A simple strategy of choosing c is to employ a curve approximating the relation between the noise level and τ. By fitting experimental data with a straight line, in this paper, we suggest setting
(33)
τ
=
-
0.006
×
BSNR
+
1.09
.
Besides the setting of τ, the choice of β
1 and β
2 is essential to our algorithm. We suggest setting β
1 = 10(BSNR/10−1) × β
2, where β
2 = 1. This parameter setting is obtained by large numbers of experiments. Actually, β
1, β
2 > 0 is sufficient for the convergence of the proposed algorithm, but why β
1 and β
2 play different important role when the BSNR varies? The reason is that, when the BSNR becomes higher, the distance between H
u and f is nearer. From minimization problem (10), we learn that auxiliary variable x plays the role of H
u and a higher BSNR means a larger β
1.


## 4. Numerical Results

In this section, two experiments are presented to demonstrate the effectiveness of the proposed method. They were performed under MATLAB v7.8.0 and Windows 7 on a PC with Intel Core (TM) i5 CUP (3.20 GHz) and 8 GB of RAM. The improved signal-to-noise ratio (ISNR) is used to measure the quality of the restoration results. It is defined as(34)
ISNR
=
10
lo
g
10
(
∥
f
-
u
clean
∥
2
2
∥
u
-
u
clean
∥
2
2
)
.
During the experiments, the four images shown in Figure 1 were used; they are named Cameraman, Lena, Shepp-Logan phantom,and Abdomen all of size 256 × 256.Figure 1
Test images:Cameraman, Lena, Shepp-Logan phantom, and Abdomen of size 256 × 256.


### 4.1. Experiment 1

In this experiment, we examine whether the regularization parameter is well estimated by the prosed algorithm. We compare APE-SBA with some famous TV-based methods in the literature and they are denoted by BFO [5], BMK [19], and LLN [20]. We make use of MATLAB commands “fspecial (“average”, 9)” and “fspecial (“Gaussian”, [9 9], 3)” to blur the Lena, Cameraman,andShepp-Logan phantom images first, and then the images are contaminated with Gaussian noises such that the BSNRs of the observed images are 20 dB, 30 dB, and 40 dB. We adopt ∥
u
k
+
1
-
u
k
∥
2
2
/
∥
u
k
∥
2
2
≤
1
0
-
6 as the stopping criteria for our algorithm, where u
k is the restored image in the kth iteration.Table1 presents the ISNRs of the restoration results of different methods. Symbol “—” means that the results are not presented in the original reference, and bold type numbers represent the best results among the four methods. From Table 1, we see that our algorithm is more competitive than the other three and only in one case our result is worse than but close to the best. This also indicates that the regularization parameter obtained by our method is good.Table 1
ISNRs obtained by different methods.
BSNR
Method
Lena
Cameraman
Shepp-Logan
ISNR (dB)
ISNR (dB)
ISNR (dB)
9 × 9 uniform blur
20
BFO [5]
4.05
3.27
6.25
BMK [19]
3.72
2.42
3.01
LLN [20]
3.15
2.88
—
APE-SBA
4.09
3.88
7.60
30
BFO [5]
5.43
5.69
10.49
BMK [19]
5.89
5.41
7.77
LLN [20]
4.43
5.57
—
APE-SBA
5.97
5.87
11.56
40
BFO [5]
6.22
8.46
16.39
BMK [19]
8.42
8.57
13.69
LLN [20]
6.92
7.86
—
APE-SBA
8.11
8.60
17.80
9 × 9 Gaussian blur
20
BFO [5]
2.99
2.21
4.24
BMK [19]
2.87
1.72
1.85
LLN [20]
2.57
1.82
—
APE-SBA
3.10
2.61
5.92
30
BFO [5]
3.82
3.59
7.21
BMK [19]
3.87
2.63
4.31
LLN [20]
4.17
3.43
—
APE-SBA
4.20
4.17
8.87
40
BFO [5]
4.41
5.78
10.27
BMK [19]
4.78
3.39
6.69
LLN [20]
5.44
5.02
—
APE-SBA
5.97
6.38
11.08


### 4.2. Experiment 2

In this subsection, we compare our algorithm with the other two state-of-the-art algorithms: the primal-dual based method in [13], named AutoRegSel, and the ADMM based method in [17], named C-SALSA. The stopping criterion of all methods is ∥
u
k
+
1
-
u
k
∥
2
2
/
∥
u
k
∥
2
2
≤
1
0
-
6 or the number of iterations is larger than 1000. We consider the three image restoration problems adopted in [17]. In the first problem, the PSF is a 9 × 9 uniform blur with noise variance 0.562 (Prob. 1); in the second problem, the PSF is a 9 × 9 Gaussian blur with noise variance 2 (Prob. 2); in the third problem, the PSF is given by h
i
,
j
=
1
/
(
1
+
i
2
+
j
2
) with noise variance 2 (Prob. 3), where i
,
j
=
-
7
,
…
,
7.The plots of ISNR (in dB) versus runtime (in second) are shown in Figure2. Table 2 presents the ISNR values, the number of iterations, and the total runtime to reach convergence. We again use the bold type numbers to represent the best results. From the results, we see that APE-SBA produces the best ISNRs compared with the other methods within the least runtime. Besides, in most cases, APE-SBA obtains the best ISNR within the least iterations. Only when dealing with the Abdomen image under Prob. 2, APE-SBA takes more iterations but less runtime to reach convergence than C-SALSA, and the total iteration number for these two is close to each other. For achieving the adaptive image restoration, both C-SALSA and AutoRegSel introduce in an inner iterative scheme, whereas APE-SBA contains no inner iteration. Obviously, the superiority in speed of our method will be enlarged when the image size becomes larger. Figure 3 shows the blurred image and the restored results by different methods in Prob. 2 of the Abdomen image. Our algorithm results in the best ISNR, and, for other problems in Experiment 2, we obtain the similar results.Table 2
Comparison between different methods in terms of ISNR, iterations, and runtime.
Problem
Method
Abdomen 256 × 256
Lena 256 × 256
ISNR (dB)
Iterations
Runtime (s)
ISNR (dB)
Iterations
Runtime (s)
Prob. 1
APE-SBA
9.63
261
2.31
7.44
201
1.78
AutoRegSel [13]
9.41
435
6.53
7.39
392
5.94
C-SALSA [17]
9.14
773
20.44
6.98
658
17.42
Prob. 2
APE-SBA
5.54
551
4.90
4.08
421
3.75
AutoRegSel [13]
5.24
855
12.92
3.96
1000
15.52
C-SALSA [17]
5.00
533
14.03
3.66
492
12.89
Prob. 3
APE-SBA
8.87
263
2.34
6.78
197
1.75
AutoRegSel [13]
8.54
414
6.26
6.65
404
6.11
C-SALSA [17]
8.20
422
11.08
6.20
507
13.29ISNR versus runtime for the (left)Abdomen image and (right) Lena image, which are blurred by a 9 × 9 uniform blur with noise variance 0.562 (first row), by a 9 × 9 Gaussian blur with noise variance 2 (second row), and by PSF given by h
i
j
=
1
/
(
1
+
i
2
+
j
2
)  (
i
,
j
=
-
7
,
…
,
7
.
) with noise variance 2 (third row).
(a)
(b)
(c)
(d)
(e)
(f)The observed image (a) which is degraded by a 9 × 9 Gaussian blur with noise variance 2, and the restored images by APE-SBA (b), by AutoRegSel (c), and by C-SALSA (d) of theAbdomen image under Prob. 2.
(a)
(b)
(c)
(d)


## 4.1. Experiment 1

In this experiment, we examine whether the regularization parameter is well estimated by the prosed algorithm. We compare APE-SBA with some famous TV-based methods in the literature and they are denoted by BFO [5], BMK [19], and LLN [20]. We make use of MATLAB commands “fspecial (“average”, 9)” and “fspecial (“Gaussian”, [9 9], 3)” to blur the Lena, Cameraman,andShepp-Logan phantom images first, and then the images are contaminated with Gaussian noises such that the BSNRs of the observed images are 20 dB, 30 dB, and 40 dB. We adopt ∥
u
k
+
1
-
u
k
∥
2
2
/
∥
u
k
∥
2
2
≤
1
0
-
6 as the stopping criteria for our algorithm, where u
k is the restored image in the kth iteration.Table1 presents the ISNRs of the restoration results of different methods. Symbol “—” means that the results are not presented in the original reference, and bold type numbers represent the best results among the four methods. From Table 1, we see that our algorithm is more competitive than the other three and only in one case our result is worse than but close to the best. This also indicates that the regularization parameter obtained by our method is good.Table 1
ISNRs obtained by different methods.
BSNR
Method
Lena
Cameraman
Shepp-Logan
ISNR (dB)
ISNR (dB)
ISNR (dB)
9 × 9 uniform blur
20
BFO [5]
4.05
3.27
6.25
BMK [19]
3.72
2.42
3.01
LLN [20]
3.15
2.88
—
APE-SBA
4.09
3.88
7.60
30
BFO [5]
5.43
5.69
10.49
BMK [19]
5.89
5.41
7.77
LLN [20]
4.43
5.57
—
APE-SBA
5.97
5.87
11.56
40
BFO [5]
6.22
8.46
16.39
BMK [19]
8.42
8.57
13.69
LLN [20]
6.92
7.86
—
APE-SBA
8.11
8.60
17.80
9 × 9 Gaussian blur
20
BFO [5]
2.99
2.21
4.24
BMK [19]
2.87
1.72
1.85
LLN [20]
2.57
1.82
—
APE-SBA
3.10
2.61
5.92
30
BFO [5]
3.82
3.59
7.21
BMK [19]
3.87
2.63
4.31
LLN [20]
4.17
3.43
—
APE-SBA
4.20
4.17
8.87
40
BFO [5]
4.41
5.78
10.27
BMK [19]
4.78
3.39
6.69
LLN [20]
5.44
5.02
—
APE-SBA
5.97
6.38
11.08


## 4.2. Experiment 2

In this subsection, we compare our algorithm with the other two state-of-the-art algorithms: the primal-dual based method in [13], named AutoRegSel, and the ADMM based method in [17], named C-SALSA. The stopping criterion of all methods is ∥
u
k
+
1
-
u
k
∥
2
2
/
∥
u
k
∥
2
2
≤
1
0
-
6 or the number of iterations is larger than 1000. We consider the three image restoration problems adopted in [17]. In the first problem, the PSF is a 9 × 9 uniform blur with noise variance 0.562 (Prob. 1); in the second problem, the PSF is a 9 × 9 Gaussian blur with noise variance 2 (Prob. 2); in the third problem, the PSF is given by h
i
,
j
=
1
/
(
1
+
i
2
+
j
2
) with noise variance 2 (Prob. 3), where i
,
j
=
-
7
,
…
,
7.The plots of ISNR (in dB) versus runtime (in second) are shown in Figure2. Table 2 presents the ISNR values, the number of iterations, and the total runtime to reach convergence. We again use the bold type numbers to represent the best results. From the results, we see that APE-SBA produces the best ISNRs compared with the other methods within the least runtime. Besides, in most cases, APE-SBA obtains the best ISNR within the least iterations. Only when dealing with the Abdomen image under Prob. 2, APE-SBA takes more iterations but less runtime to reach convergence than C-SALSA, and the total iteration number for these two is close to each other. For achieving the adaptive image restoration, both C-SALSA and AutoRegSel introduce in an inner iterative scheme, whereas APE-SBA contains no inner iteration. Obviously, the superiority in speed of our method will be enlarged when the image size becomes larger. Figure 3 shows the blurred image and the restored results by different methods in Prob. 2 of the Abdomen image. Our algorithm results in the best ISNR, and, for other problems in Experiment 2, we obtain the similar results.Table 2
Comparison between different methods in terms of ISNR, iterations, and runtime.
Problem
Method
Abdomen 256 × 256
Lena 256 × 256
ISNR (dB)
Iterations
Runtime (s)
ISNR (dB)
Iterations
Runtime (s)
Prob. 1
APE-SBA
9.63
261
2.31
7.44
201
1.78
AutoRegSel [13]
9.41
435
6.53
7.39
392
5.94
C-SALSA [17]
9.14
773
20.44
6.98
658
17.42
Prob. 2
APE-SBA
5.54
551
4.90
4.08
421
3.75
AutoRegSel [13]
5.24
855
12.92
3.96
1000
15.52
C-SALSA [17]
5.00
533
14.03
3.66
492
12.89
Prob. 3
APE-SBA
8.87
263
2.34
6.78
197
1.75
AutoRegSel [13]
8.54
414
6.26
6.65
404
6.11
C-SALSA [17]
8.20
422
11.08
6.20
507
13.29ISNR versus runtime for the (left)Abdomen image and (right) Lena image, which are blurred by a 9 × 9 uniform blur with noise variance 0.562 (first row), by a 9 × 9 Gaussian blur with noise variance 2 (second row), and by PSF given by h
i
j
=
1
/
(
1
+
i
2
+
j
2
)  (
i
,
j
=
-
7
,
…
,
7
.
) with noise variance 2 (third row).
(a)
(b)
(c)
(d)
(e)
(f)The observed image (a) which is degraded by a 9 × 9 Gaussian blur with noise variance 2, and the restored images by APE-SBA (b), by AutoRegSel (c), and by C-SALSA (d) of theAbdomen image under Prob. 2.
(a)
(b)
(c)
(d)


## 5. Conclusions

We developed a split Bregman based algorithm to solve the TV image restoration/deconvolution problem. Unlike some other methods in the literature, without any inner iteration, our method achieves the updating of the regularization parameter and the restoration of the blurred image simultaneously, by referring to the operator splitting technique and introducing two auxiliary variables for both the data-fidelity term and the TV regularization term. Therefore, the algorithm can run without any manual interference. The numerical results have indicated that the proposed algorithm outperforms some state-of-the-art methods in both speed and accuracy.



---

*Source: 617026-2014-02-17.xml*