# Poissonian Image Deconvolution via Sparse and Redundant Representations and Framelet Regularization

Authors: Yu Shi; Houzhang Fang; Guoyou Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Poissonian image deconvolution is a key issue in various applications, such as astronomical imaging, medical imaging, and electronic microscope imaging. A large amount of literature on this subject is analysis-based methods. These methods assign various forward measurements of the image. Meanwhile, synthesis-based methods are another well-known class of methods. These methods seek a reconstruction of the image. In this paper, we propose an approach that combines analysis with synthesis methods. The method is proposed to address Poissonian image deconvolution problem by minimizing the energy functional, which is composed of a sparse representation prior over a learned dictionary, the data fidelity term, and framelet based analysis prior constraint as the regularization term. The minimization problem can be efficiently solved by the split Bregman technique. Experiments demonstrate that our approach achieves better results than many state-of-the-art methods, in terms of both restoration accuracy and visual perception.


---

## Body

## 1. Introduction

Poissonian image deconvolution appears in various applications such as astronomical imaging [1], medical imaging [2], and electronic microscope imaging [3]. It aims to reconstruct a high quality image u from the degraded image g. Mathematically, the process of Poissonian image deconvolution can be generally modeled by
(1)
g
=
𝒫
(
H
u
)
,
where H
∈
ℝ
N
×
N denotes the point spread function (PSF), u
∈
ℝ
N denotes the unknown image to be estimated, 𝒫 denotes the Poisson noise process, g
∈
ℝ
N denotes the blurred noisy image, and N is the length of image vector which is stacked by columns.It is known that Poissonian image deconvolution is a typical ill-posed inverse problem. In general, the solution of (1) is not unique. Prior knowledge of image, including analysis-based and synthesis-based priors, can be used to address this problem. For an overview of the two classes of priors, we refer to [4]. Analysis-based priors are frequently used as regularization term in energy functional E
(
u
)
(2)
E
(
u
)
=
inf
⁡
u
∫
Ω
H
u
-
g
log
⁡
⁡
H
u
d
x
+
R
(
u
)
,
where R
(
u
) is the regularization constraint term. Two main analysis-based methods have been proposed to solve problem (2): the total variation (TV) [5, 6] based methods and wavelet frames-based methods. TV-based methods have shown good performance on blurred images for discontinuous solution and edge-preserving advantages. Many authors have combined the TV regularization term with variant methods to solve problem (2). For example, Sawatzky et al. combined expectation-maximization algorithm with TV regularization [7] in positron emission tomography. Setzer et al. considered using TV regularization term with split Bregman techniques [8]. TV-based methods can remove noise in flat region effectively. However, they lose details in texture region. Furthermore, the inherent problem of TV-based methods is that they favor a piecewise constant solution, which causes staircase effects in flat region.To overcome disadvantages of TV-based methods, several authors are interested in frame-based methods that can be regarded as analysis regularization methods as well. It is assumed that majority of nature images have sparse approximation under some redundant tight frame systems including wavelet transform, Gabor transform, curvelets, and framelets [9]. Chen and Cheng considered a hybrid variational model [10] that describes a cartoon part by TV and a detailed part by the sparse representation over the wavelet basis. Figueiredo and Bioucas-Dias proposed an alternating direction method [11] using frame-based analysis regularizer. Fang et al. presented blind Poissonian images deconvolution with framelet regularization [12]. He considered framelet-based analysis prior and combination of framelet and total variation analysis priors in [13] as well.Dictionary learning [14–20] based sparse representation has been extensively studied, and significant progress has been made in the development of modern dictionaries learning for synthesis-based sparse representation of image. Natural images contain repeated patterns such as texture region, flat region, and smooth region. Hence, they can be described as linear combinations of a few atoms from a dictionary. Dictionary learning for synthesis model in image denoising was proposed by Elad and Aharon. They presented the method via sparse and redundant representations over learned dictionaries, called K-SVD [16]. This method is performed to remove the additive noise. Huang et al. attempted to use the K-SVD method [17] for multiplicative noise removal which followed. It was proven that K-SVD method could perform well just by transforming the multiplicative noise to additive noise in logarithmic domain. Similarly, Poisson noise can be transformed to additive Gaussian noise with unitary variance by Anscombe root transformation (ART) [18] and then removed by K-SVD method [15]. In fact, the variance stabilization often has been questioned for the poor numerical results. Therefore, the denoising algorithm should be specifically designed for Poisson noise. Ma et al. proposed dictionary learning for Poisson noise removal [21].We still have work to do, on dictionary learning for Poisson noise removal, because TV regularizer and dictionary learning (refer to [21]) do not match well. TV regularizer can reduce artifacts in smooth regions caused by patch-based priors of dictionary learning; however, it smooths the fine details in texture region. Motivated by analysis and synthesis based methods, we combine analysis with synthesis method for Poissonian image deconvolution. We choose B-spline framelet as the analysis regularization term, because framelet transform has the ability of multiple-resolution analysis. Different framelet masks reflect different orders of difference operators, which can adaptively capture multiscale edge structures in images. Besides, synthesis-based method is integrated into the proposed method using a sparse representation prior over a learned dictionary. At last, split Bergman method is used to address optimization problem. Bergman variables make the proposed method stable and regularization parameters do not need to be updated in Bregman iteration process.The rest of this paper is organized as follows. The K-SVD algorithm and B-spline framelet are briefly reviewed in Section2. The Poissionan image deconvolution model is proposed in Section 3. Numerical experiments are given to demonstrate the qualities of the restored images in Section 4. Conclusions are given in Section 5.


## 2. Review of Prior Works

In this section, we review K-SVD algorithm and B-spline framelet briefly. K-SVD algorithm is a classical sparse representation method, and it will be used in our proposed model. Framelet draws many attentions by its multiscale characteristic.


### 2.1. K-SVD Algorithm

We consider the image patchx of size n
×
n and define dictionary D
∈
ℝ
n
×
k, where k
>
n. For each image patch, x can be represented sparsely over the dictionary; that is, the solution of
(3)
α
^
=
arg
min
⁡
α
∥
α
∥
0
subject
to
D
α
≈
x
is very sparse. The notation ∥
α
∥
0 stands for the count of the nonzero entries in α. Considering the image patch x contaminated by an additive zero-mean white Gaussian noise with standard deviation δ, the MAP estimator for denoising this image patch is built by solving
(4)
α
^
=
arg
min
⁡
α
∥
D
α
-
x
∥
2
2
+
μ
∥
α
∥
0
.
If each image patch x of image X can be described by (4), the MAP estimator of image X will be described by
(5)
α
^
i
j
,
X
^
=
arg
min
⁡
α
i
j
,
X
λ
∥
X
-
Y
∥
2
2
+
∑
i
j
‍
∥
D
α
i
j
-
R
i
j
X
∥
2
2
+
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
0
,
where Y denotes the degraded image and each image patch x
i
j
=
R
i
j
X, R
i
j is an n
×
N matrix that extracts the (
i
j
) block from the image X. Instead of addressing both unknowns α
^
i
j and X together, a block-coordinate minimization algorithm is used: (i)
initializationX
=
Y, D
= overcomplete DCT dictionary;
(ii)
repeatJ times:
(1)
sparse coding stage: use any pursuit algorithm to seek the optimalα
^
i
j, by approximating the solution of
(6)
∀
i
j
min
⁡
α
i
j
∥
α
i
j
∥
0
s
.
t
.
∥
D
α
i
j
-
R
i
j
X
∥
2
2
≤
ϵ
2
;
(2)
dictionary update stage: for each columnl
=
1
,
2
,
…
,
k in D, and α
i
j
≠
0, compute the representation error
(7)
e
i
j
l
=
R
i
j
X
i
j
-
∑
m
≠
l
‍
d
m
α
i
j
(
m
)
and update the dictionary column d
~
l by SVD decomposition of the representation error;
(iii)
Given allα
^
i
j, update X:
(8)
X
^
=
arg
min
⁡
X
λ
∥
X
-
Y
∥
2
2
+
∑
i
j
‍
∥
D
α
^
i
j
-
R
i
j
X
∥
2
2
;
the simple quadratic term has a closed-form solution of the form
(9)
X
^
=
(
λ
I
+
∑
i
j
R
i
j
T
R
i
j
)
-
1
(
λ
Y
+
∑
i
j
R
i
j
T
D
α
^
i
j
)
.


### 2.2. Analysis Operator-B-Spline Framelet

In this subsection, we will briefly introduce the theory of frame and B-spline framelet [22]. A countable subset of X
(
Ψ
)
=
{
ψ
1
,
ψ
2
,
…
,
ψ
n
}
⊂
L
2
(
ℝ
) will be called a frame of L
2
(
ℝ
), if it satisfies
(10)
A
∥
f
∥
2
≤
∑
ψ
n
∈
Ψ
|
〈
f
,
ψ
n
〉
|
2
≤
B
∥
f
∥
2
,
∀
f
∈
L
2
(
ℝ
)
,
where 〈
·
,
·
〉 is the inner product of L
2
(
ℝ
), f is a signal, and A and B are positive constants. If A
=
B, we call it tight frame. That is to say, f can be perfectly reconstructed by
(11)
f
=
∑
ψ
n
∈
Ψ
〈
f
,
ψ
n
〉
ψ
n
,
∀
f
∈
L
2
(
ℝ
)
.Furthermore, wavelet frame is constructed by wavelet transform using time and scale parameters sampling. It is defined as(12)
X
(
Ψ
)
=
{
ψ
u
,
s
=
1
s
ψ
(
t
-
u
s
)
}
.If the wavelet frame satisfies the tight frame condition, functionsψ
∈
Ψ are called framelets which belong to the wavelet tight frame X
(
Ψ
). To construct a set of framelets, one starts from a refinement mask h of a refinable function ϕ (also called a scaling function) satisfying
(13)
ϕ
^
(
2
·
)
=
h
ϕ
^
,
where ϕ
^ is the Fourier transform of ϕ. A wavelet tight frame can be constructed by an appropriate set of framelets Ψ
=
{
ψ
1
,
ψ
2
,
…
,
ψ
n
} that can be described in the Fourier domain as
(14)
ψ
^
i
(
2
·
)
=
h
i
ϕ
^
,
i
=
1,2
,
…
,
n
.B-spline tight framelet is constructed by piecewise linear B-spline functionϕ and refinement mask h
0
=
cos
⁡
2
(
ω
). The corresponding low-pass filter is
(15)
h
0
=
1
4
[
1,2
,
1
]
.Two frameletsψ
1
,
ψ
2 can be defined by the framelet masks h
1
=
-
(
2
i
/
2
)
sin
(
ω
) and h
2
=
sin
2
(
ω
/
2
), and the corresponding high-pass filters are
(16)
h
1
=
2
4
[
1,0
,
-
1
]
,
h
2
=
1
4
[
-
1,2
,
-
1
]
.Subsequently, the framelet decomposition off can be easily obtained. The decomposition process which also can be called wavelet transform is defined as
(17)
W
f
(
u
,
s
)
=
〈
f
,
ψ
u
,
s
〉
,
where W is the transform matrix. In our paper, W is called the analysis operator associated with 1-tight (Parseval) frame. We have W
T
W
=
I, in discrete domain, and W and W
T denote the discrete wavelet tight frame decomposition and reconstruction, respectively. More details about the properties of frame and framelet can be found in [23].


## 2.1. K-SVD Algorithm

We consider the image patchx of size n
×
n and define dictionary D
∈
ℝ
n
×
k, where k
>
n. For each image patch, x can be represented sparsely over the dictionary; that is, the solution of
(3)
α
^
=
arg
min
⁡
α
∥
α
∥
0
subject
to
D
α
≈
x
is very sparse. The notation ∥
α
∥
0 stands for the count of the nonzero entries in α. Considering the image patch x contaminated by an additive zero-mean white Gaussian noise with standard deviation δ, the MAP estimator for denoising this image patch is built by solving
(4)
α
^
=
arg
min
⁡
α
∥
D
α
-
x
∥
2
2
+
μ
∥
α
∥
0
.
If each image patch x of image X can be described by (4), the MAP estimator of image X will be described by
(5)
α
^
i
j
,
X
^
=
arg
min
⁡
α
i
j
,
X
λ
∥
X
-
Y
∥
2
2
+
∑
i
j
‍
∥
D
α
i
j
-
R
i
j
X
∥
2
2
+
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
0
,
where Y denotes the degraded image and each image patch x
i
j
=
R
i
j
X, R
i
j is an n
×
N matrix that extracts the (
i
j
) block from the image X. Instead of addressing both unknowns α
^
i
j and X together, a block-coordinate minimization algorithm is used: (i)
initializationX
=
Y, D
= overcomplete DCT dictionary;
(ii)
repeatJ times:
(1)
sparse coding stage: use any pursuit algorithm to seek the optimalα
^
i
j, by approximating the solution of
(6)
∀
i
j
min
⁡
α
i
j
∥
α
i
j
∥
0
s
.
t
.
∥
D
α
i
j
-
R
i
j
X
∥
2
2
≤
ϵ
2
;
(2)
dictionary update stage: for each columnl
=
1
,
2
,
…
,
k in D, and α
i
j
≠
0, compute the representation error
(7)
e
i
j
l
=
R
i
j
X
i
j
-
∑
m
≠
l
‍
d
m
α
i
j
(
m
)
and update the dictionary column d
~
l by SVD decomposition of the representation error;
(iii)
Given allα
^
i
j, update X:
(8)
X
^
=
arg
min
⁡
X
λ
∥
X
-
Y
∥
2
2
+
∑
i
j
‍
∥
D
α
^
i
j
-
R
i
j
X
∥
2
2
;
the simple quadratic term has a closed-form solution of the form
(9)
X
^
=
(
λ
I
+
∑
i
j
R
i
j
T
R
i
j
)
-
1
(
λ
Y
+
∑
i
j
R
i
j
T
D
α
^
i
j
)
.


## 2.2. Analysis Operator-B-Spline Framelet

In this subsection, we will briefly introduce the theory of frame and B-spline framelet [22]. A countable subset of X
(
Ψ
)
=
{
ψ
1
,
ψ
2
,
…
,
ψ
n
}
⊂
L
2
(
ℝ
) will be called a frame of L
2
(
ℝ
), if it satisfies
(10)
A
∥
f
∥
2
≤
∑
ψ
n
∈
Ψ
|
〈
f
,
ψ
n
〉
|
2
≤
B
∥
f
∥
2
,
∀
f
∈
L
2
(
ℝ
)
,
where 〈
·
,
·
〉 is the inner product of L
2
(
ℝ
), f is a signal, and A and B are positive constants. If A
=
B, we call it tight frame. That is to say, f can be perfectly reconstructed by
(11)
f
=
∑
ψ
n
∈
Ψ
〈
f
,
ψ
n
〉
ψ
n
,
∀
f
∈
L
2
(
ℝ
)
.Furthermore, wavelet frame is constructed by wavelet transform using time and scale parameters sampling. It is defined as(12)
X
(
Ψ
)
=
{
ψ
u
,
s
=
1
s
ψ
(
t
-
u
s
)
}
.If the wavelet frame satisfies the tight frame condition, functionsψ
∈
Ψ are called framelets which belong to the wavelet tight frame X
(
Ψ
). To construct a set of framelets, one starts from a refinement mask h of a refinable function ϕ (also called a scaling function) satisfying
(13)
ϕ
^
(
2
·
)
=
h
ϕ
^
,
where ϕ
^ is the Fourier transform of ϕ. A wavelet tight frame can be constructed by an appropriate set of framelets Ψ
=
{
ψ
1
,
ψ
2
,
…
,
ψ
n
} that can be described in the Fourier domain as
(14)
ψ
^
i
(
2
·
)
=
h
i
ϕ
^
,
i
=
1,2
,
…
,
n
.B-spline tight framelet is constructed by piecewise linear B-spline functionϕ and refinement mask h
0
=
cos
⁡
2
(
ω
). The corresponding low-pass filter is
(15)
h
0
=
1
4
[
1,2
,
1
]
.Two frameletsψ
1
,
ψ
2 can be defined by the framelet masks h
1
=
-
(
2
i
/
2
)
sin
(
ω
) and h
2
=
sin
2
(
ω
/
2
), and the corresponding high-pass filters are
(16)
h
1
=
2
4
[
1,0
,
-
1
]
,
h
2
=
1
4
[
-
1,2
,
-
1
]
.Subsequently, the framelet decomposition off can be easily obtained. The decomposition process which also can be called wavelet transform is defined as
(17)
W
f
(
u
,
s
)
=
〈
f
,
ψ
u
,
s
〉
,
where W is the transform matrix. In our paper, W is called the analysis operator associated with 1-tight (Parseval) frame. We have W
T
W
=
I, in discrete domain, and W and W
T denote the discrete wavelet tight frame decomposition and reconstruction, respectively. More details about the properties of frame and framelet can be found in [23].


## 3. Sparse and Redundant Representation on Dictionary Learning and Framelet

### 3.1. Proposed Model

Inspired by K-SVD for dictionary learning and the split Bregman method, we propose a minimization model based on the sparse representation of the image over a certain unknown dictionary and analysis prior. The proposed model can be described by(18)
min
⁡
α
,
u
〈
1
,
H
u
-
g
log
⁡
H
u
〉
+
λ
∥
W
u
∥
1
+
τ
∑
i
j
‍
∥
R
i
j
u
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
,
where λ
,
τ are the positive regularization parameters. Among these parameters, τ is an important parameter. It can reduce the workload of parameter setting drastically. The effect can be seen in parameter setting analysis at the end of this section. ∥
·
∥
1 is the l
1-norm that sums the absolute values of a vector. The sparse solution of (1) can be approximated by l
1-norm regularizer. The l
1-norm becomes popular in recent years. It is a more robust measurement which gives suitable penalty on flat region or edge. Besides, it fully embodies the sparsity. Most of all, it is easy to be computed. The notation ∥
W
u
∥
1 stands for the analysis-based regularizer. The framelet transform W is generated by the piecewise linear B-splines framelet filter masks for its simplicity and effectiveness in image deblurring.


### 3.2. Alternating Split Bergman Algorithm

Split Bergman algorithm [9, 24] has shown good performance on splitting the seemingly unsolved variable optimization problem into number of easy subproblems. Considering that the general unconstrained minimization problems have the form of
(19)
min
⁡
u
∈
ℝ
N
Φ
(
u
)
+
Ψ
(
C
u
)
,
where C
:
ℝ
N
→
ℝ
M is a bounded linear operator and both functions Φ
:
ℝ
N
→
ℝ
∪
{
+
∞
} and Ψ
:
ℝ
M
→
ℝ
∪
{
+
∞
} are proper, convex, and lsc, it is assumed that (19) has a solution and 0
∈
int
⁡
(
C
dom
⁡
(
Φ
)
-
dom
⁡
(
Ψ
)
). To derive the split Bregman algorithm, problem (19) is transformed into the constrained minimization problem
(20)
min
⁡
u
∈
ℝ
N
,
d
∈
ℝ
M
Φ
(
u
)
+
Ψ
(
d
)
s
.
t
.
C
u
=
d
.Using augmented Lagrangian method, the split Bregman iteration for (19) is defined as
(21)
(
u
(
k
+
1
)
,
d
(
k
+
1
)
)
=
arg
min
⁡
u
∈
ℝ
N
,
d
∈
ℝ
M
{
Φ
(
u
)
+
Ψ
(
d
)
+
1
2
γ
∥
C
u
-
d
∥
2
2
}
,
b
(
k
+
1
)
=
b
(
k
)
+
1
γ
(
C
u
(
k
+
1
)
-
d
(
k
+
1
)
)
.The alternating split Bregman algorithm is proposed to minimize the (21) alternatingly with respect to u and d:
(22)
u
(
k
+
1
)
=
arg
min
⁡
u
∈
ℝ
N
{
Φ
(
u
)
+
1
2
γ
∥
b
(
k
)
+
C
u
-
d
(
k
)
∥
2
2
}
,
d
(
k
+
1
)
=
arg
min
⁡
d
∈
ℝ
M
{
Ψ
(
d
)
+
1
2
γ
∥
b
(
k
)
+
C
u
(
k
+
1
)
-
d
∥
2
2
}
,
b
(
k
+
1
)
=
b
(
k
)
+
C
u
(
k
+
1
)
-
d
(
k
+
1
)
.It is important to note that the penalty parameter1
/
(
2
γ
) should not be too large. If the parameter becomes very large, the minimization subproblems will become increasingly ill-conditioned, and it will lead to numerical problems [25].If we use the alternating split Bergman algorithm directly, it will cause a little problem. Letingd
1
=
H
u and d
2
=
W
u, we can obtain the minimization problem about variables u
,
α
,
D
,
d
1
,
d
2
(23)
min
⁡
D
,
α
,
d
1
,
d
2
,
u
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
λ
∥
d
2
∥
1
+
τ
∑
i
j
∥
R
i
j
u
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
μ
i
j
∥
α
i
j
∥
1
+
1
2
ρ
∥
b
1
+
H
u
-
d
1
∥
2
2
+
1
2
β
∥
b
2
+
W
u
-
d
2
∥
2
2
.Supposing thatD
,
α
,
d
1
,
d
2 are known, the solution with respect to u can be transformed to
(24)
min
⁡
u
τ
∑
i
j
‍
∥
R
i
j
u
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
+
1
2
ρ
∥
b
1
+
H
u
-
d
1
∥
2
2
+
1
2
β
∥
b
2
+
W
u
-
d
2
∥
2
2
.Obviously, it is a quadratic term. We can obtain the linear equation ofu by its closed form(25)
u
^
=
(
1
ρ
H
T
(
d
1
-
b
1
)
+
1
β
W
T
(
d
2
-
b
2
)
+
τ
∑
i
j
R
i
j
T
D
α
i
j
)
×
(
1
ρ
H
T
H
+
1
β
W
T
W
+
τ
∑
i
j
R
i
j
T
R
i
j
)
-
1
.The value of term about variableD is large, and other elements are so small that can be ignored. Thus, u
^ is about the same as the following form:
(26)
u
^
≈
∑
i
j
R
i
j
T
D
α
i
j
∑
i
j
R
i
j
T
R
i
j
.In [16], it has proved that learning dictionary using K-SVD can deal with the additive noise. However, it cannot deal with the Poisson noise directly. Here, we add the special constraint d
3
=
u to address this difficult problem. The special constraint comes from the idea which combines the split Bregman algorithm with an alternating minimization algorithm; refer to [26]. It is shown that this special alternating optimization algorithm simplifies the complex problem. Moreover, it solves the actual problem that K-SVD algorithm is only performed on additive noise denoising or multiplicative noise denoising, and multiplicative noise can be transformed to the additive noise in logarithm domain [17].


### 3.3. Proposed Algorithm

The difficulties in seeking variables are nonquadratic data term and nonseparablel
1-norm term. The alternating optimization algorithm can be used to solve these difficulties by introducing auxiliary variables d
1
=
H
u
,
d
2
=
W
u
,
d
3
=
u. The minimization problem becomes
(27)
min
⁡
α
,
D
,
d
1
,
d
2
,
d
3
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
λ
∥
d
2
∥
1
+
τ
∑
i
j
‍
∥
R
i
j
d
3
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
s
.
t
.
d
1
=
H
u
,
d
2
=
W
u
,
d
3
=
u
.Here, the alternating optimization algorithm is a special split Bregman algorithm. It is to be noted thatd
3 replaced u in the sparse representation of the image over a certain unknown dictionary.The proposed model is finally defined as(28)
min
⁡
α
,
D
,
u
,
d
1
,
d
2
,
d
3
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
λ
∥
d
2
∥
1
+
τ
∑
i
j
‍
∥
R
i
j
d
3
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
+
1
2
ρ
∥
b
1
+
H
u
-
d
1
∥
2
2
+
1
2
β
∥
b
2
+
W
u
-
d
2
∥
2
2
+
1
2
γ
∥
b
3
+
u
-
d
3
∥
2
2
,
where ρ
,
β
,
γ are the positive penalty parameters. The minimization problem (28) can be decoupled into several subproblems via split Bregman iteration.(1) Givend
1
(
k
)
,
d
2
(
k
)
,
d
3
(
k
)
,
b
1
(
k
)
,
b
2
(
k
)
,
b
3
(
k
), the minimization of (28) with respect to u is transformed to
(29)
min
⁡
u
1
2
ρ
∥
b
1
(
k
)
+
H
u
-
d
1
(
k
)
∥
2
2
+
1
2
β
∥
b
2
(
k
)
+
W
u
-
d
2
(
k
)
∥
2
2
+
1
2
γ
∥
b
3
(
k
)
+
u
-
d
3
(
k
)
∥
2
2
.According to its closed form, we can obtain the linear equation ofu
(
k
+
1
)
(30)
u
(
k
+
1
)
=
(
1
ρ
H
T
(
d
1
(
k
)
-
b
1
(
k
)
)
+
1
β
W
T
(
d
2
(
k
)
-
b
2
(
k
)
)
M
+
1
γ
(
d
3
(
k
)
-
b
3
(
k
)
)
)
×
(
1
ρ
H
T
H
+
1
β
W
T
W
+
1
γ
)
-
1
.
The convolution of H can be easily performed in FFT domain.(2) Givenb
1
(
k
)
,
u
(
k
+
1
), the subproblem with respect to d
1 becomes
(31)
min
⁡
d
1
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
1
2
ρ
∥
b
1
(
k
)
+
H
u
(
k
+
1
)
-
d
1
∥
2
2
.Thend
1 can be obtained via solving the corresponding Euler-Lagrange equation of functional (31)
(32)
d
1
(
k
+
1
)
=
1
2
(
(
b
1
(
k
)
+
H
u
(
k
+
1
)
-
ρ
)
2
+
4
ρ
g
b
1
(
k
)
+
H
u
(
k
+
1
)
-
ρ
m
m
+
(
b
1
(
k
)
+
H
u
(
k
+
1
)
-
ρ
)
2
+
4
ρ
g
)
.(3) Givenb
2
(
k
)
,
u
(
k
+
1
), the third subproblem leads to a l
1-norm problem
(33)
min
⁡
d
2
λ
∥
d
2
∥
1
+
1
2
β
∥
b
2
(
k
)
+
W
u
(
k
+
1
)
-
d
2
∥
2
2
.The shrinkage operator is used to solve problem (31). Supposed that S
λ
(
t
)
=
sign
⁡
(
t
)
(
|
t
|
-
λ
)
+, d
2
(
k
+
1
)
=
S
λ
β
(
W
u
(
k
+
1
)
+
b
2
(
k
)
)
+, d
2
(
k
+
1
) is expanded,
(34)
d
2
(
k
+
1
)
=
max
⁡
{
∥
W
u
(
k
+
1
)
+
b
2
(
k
)
∥
1
-
λ
β
,
0
}
W
u
(
k
+
1
)
+
b
2
(
k
)
∥
W
u
(
k
+
1
)
+
b
2
(
k
)
∥
2
2
.(4) Givenu
(
k
+
1
)
,
d
3
(
k
)
,
D, the sparse representation coefficients {
α
i
j
} can be obtained by solving
(35)
min
⁡
α
i
j
∑
i
j
∥
R
i
j
d
3
(
k
)
-
D
α
i
j
∥
2
2
+
∑
i
j
μ
i
j
∥
α
i
j
∥
1
.Comparing (35) with (6), it can be found that (35) is the sparse coding stage of K-SVD algorithm. Any pursuit algorithm can be used to compute the representation vectors α
^
i
j.(5) Givenα
^
i
j
,
D
,
u
(
k
+
1
)
,
b
3
(
k
), the last subproblem with respect to d
3 is
(36)
min
⁡
d
3
τ
∑
i
j
‍
∥
R
i
j
d
3
-
D
α
^
i
j
∥
2
2
+
1
2
γ
∥
b
3
(
k
)
+
u
(
k
+
1
)
-
d
3
∥
2
2
.This is also a linear equation problem; similar to (33), we can obtain
(37)
d
3
(
k
+
1
)
=
(
1
γ
τ
I
+
2
∑
i
j
‍
R
i
j
T
R
i
j
)
-
1
×
(
1
γ
τ
(
b
3
(
k
)
+
u
(
k
+
1
)
)
+
2
∑
i
j
‍
R
i
j
T
D
α
^
i
j
)
.(6) Givenb
1
(
k
)
,
b
2
(
k
)
,
b
3
(
k
)
,
d
1
(
k
+
1
)
,
d
2
(
k
+
1
)
,
d
3
(
k
+
1
)
,
u
(
k
+
1
), the Bregman auxiliary variables b
1
,
b
2
,
b
3 are updated as follows:
(38)
b
1
(
k
+
1
)
=
b
1
(
k
)
+
H
u
(
k
+
1
)
-
d
1
(
k
+
1
)
,
b
2
(
k
+
1
)
=
b
2
(
k
)
+
W
u
(
k
+
1
)
-
d
2
(
k
+
1
)
,
b
3
(
k
+
1
)
=
b
3
(
k
)
+
u
(
k
+
1
)
-
d
3
(
k
+
1
)
.Besides, the dictionary learning is also an important procedure on sparse and redundant representation. We update the dictionaryD at the last step. The D is updated in an outer loop for saving computational time as suggested in [21]. Therefore, the loop times are divided into inner loop for main process and outside loop only for dictionary updating.Taking into account the above equations, the proposed Poissonian image deconvolution algorithm is summarized as shown in Algorithm1.Algorithm 1:Poissonian image deconvolution via sparse and redundant representations and framelet regularization.
Initialize: Initialization: u
0
=
g, d
1
0
=
H
g, d
2
0
=
W
g, d
3
0
=
g, D
=
DCT, and
b
1
0
=
b
2
0
=
b
3
0. Set ite
r
O
:
=
0, ite
r
I
:
=
0.
for  ite
r
O
=
0
to
N
O  do
while  stopping criterions is not satisfied  do
Step 1. Update u
(
k
+
1
) using (30);
Step 2. Update d
1
(
k
+
1
) using (32);
Step 3. Update d
2
(
k
+
1
) using (34);
Step 4. Update d
3
(
k
+
1
) using (37);
Step 5. Update sparse representation coefficients {
α
i
j
} using (6);
Step 6. Update the Bregman auxiliary variables b
1
(
k
+
1
)
,
b
2
(
k
+
1
)
,
b
3
(
k
+
1
);
ite
r
I
=
ite
r
I
+
1;
end
Step 7. Update Dictionary D;
endOne of the stopping criterions is calculated for each iteration by the “normalized step difference energy” (NSDE): NSDE
=
∥
u
(
k
)
-
u
(
k
-
1
)
∥
2
2
/
∥
u
(
k
)
∥
2
2, where u
(
k
) and u
(
k
-
1
) denote the image vector at the kth and (
k
-
1
)th iteration, respectively. The other stopping criterion is k
≤
N
I, where N
I is the maximum inner loop times. Each of the stopping criterion is satisfied; the inner loop will be terminated. In our paper, we set NSDE
≤
0.00001, N
I
=
40, and N
O
=
10. Besides, some conditions should be satisfied to obtain a sound deconvolution result. For example, the restored image u should be normalized, and it should be nonnegative in the process of iteration as we set u
=
max
⁡
⁡
(
u
,
0
).The proposed method involves many penalty parameters. Parameters choosing plays an important role in efficiency and stability of the algorithm. There are some skills to determinate these parameters. (1) In many literatures,ρ
,
β
,
γ can be equal to each other. To eliminate numerical problems, the three parameters should not be too small. (2) The value of ρ is the larger the better in practice. In our paper, we set ρ
∈
[
400,10000
]. We can see that λ
β is an important parameters group in the process of framelet regularization. The value of group λ
β belongs to empirical range [
3,10
]. The last parameters group is τ
γ which belongs to empirical range [
1,5
].


## 3.1. Proposed Model

Inspired by K-SVD for dictionary learning and the split Bregman method, we propose a minimization model based on the sparse representation of the image over a certain unknown dictionary and analysis prior. The proposed model can be described by(18)
min
⁡
α
,
u
〈
1
,
H
u
-
g
log
⁡
H
u
〉
+
λ
∥
W
u
∥
1
+
τ
∑
i
j
‍
∥
R
i
j
u
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
,
where λ
,
τ are the positive regularization parameters. Among these parameters, τ is an important parameter. It can reduce the workload of parameter setting drastically. The effect can be seen in parameter setting analysis at the end of this section. ∥
·
∥
1 is the l
1-norm that sums the absolute values of a vector. The sparse solution of (1) can be approximated by l
1-norm regularizer. The l
1-norm becomes popular in recent years. It is a more robust measurement which gives suitable penalty on flat region or edge. Besides, it fully embodies the sparsity. Most of all, it is easy to be computed. The notation ∥
W
u
∥
1 stands for the analysis-based regularizer. The framelet transform W is generated by the piecewise linear B-splines framelet filter masks for its simplicity and effectiveness in image deblurring.


## 3.2. Alternating Split Bergman Algorithm

Split Bergman algorithm [9, 24] has shown good performance on splitting the seemingly unsolved variable optimization problem into number of easy subproblems. Considering that the general unconstrained minimization problems have the form of
(19)
min
⁡
u
∈
ℝ
N
Φ
(
u
)
+
Ψ
(
C
u
)
,
where C
:
ℝ
N
→
ℝ
M is a bounded linear operator and both functions Φ
:
ℝ
N
→
ℝ
∪
{
+
∞
} and Ψ
:
ℝ
M
→
ℝ
∪
{
+
∞
} are proper, convex, and lsc, it is assumed that (19) has a solution and 0
∈
int
⁡
(
C
dom
⁡
(
Φ
)
-
dom
⁡
(
Ψ
)
). To derive the split Bregman algorithm, problem (19) is transformed into the constrained minimization problem
(20)
min
⁡
u
∈
ℝ
N
,
d
∈
ℝ
M
Φ
(
u
)
+
Ψ
(
d
)
s
.
t
.
C
u
=
d
.Using augmented Lagrangian method, the split Bregman iteration for (19) is defined as
(21)
(
u
(
k
+
1
)
,
d
(
k
+
1
)
)
=
arg
min
⁡
u
∈
ℝ
N
,
d
∈
ℝ
M
{
Φ
(
u
)
+
Ψ
(
d
)
+
1
2
γ
∥
C
u
-
d
∥
2
2
}
,
b
(
k
+
1
)
=
b
(
k
)
+
1
γ
(
C
u
(
k
+
1
)
-
d
(
k
+
1
)
)
.The alternating split Bregman algorithm is proposed to minimize the (21) alternatingly with respect to u and d:
(22)
u
(
k
+
1
)
=
arg
min
⁡
u
∈
ℝ
N
{
Φ
(
u
)
+
1
2
γ
∥
b
(
k
)
+
C
u
-
d
(
k
)
∥
2
2
}
,
d
(
k
+
1
)
=
arg
min
⁡
d
∈
ℝ
M
{
Ψ
(
d
)
+
1
2
γ
∥
b
(
k
)
+
C
u
(
k
+
1
)
-
d
∥
2
2
}
,
b
(
k
+
1
)
=
b
(
k
)
+
C
u
(
k
+
1
)
-
d
(
k
+
1
)
.It is important to note that the penalty parameter1
/
(
2
γ
) should not be too large. If the parameter becomes very large, the minimization subproblems will become increasingly ill-conditioned, and it will lead to numerical problems [25].If we use the alternating split Bergman algorithm directly, it will cause a little problem. Letingd
1
=
H
u and d
2
=
W
u, we can obtain the minimization problem about variables u
,
α
,
D
,
d
1
,
d
2
(23)
min
⁡
D
,
α
,
d
1
,
d
2
,
u
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
λ
∥
d
2
∥
1
+
τ
∑
i
j
∥
R
i
j
u
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
μ
i
j
∥
α
i
j
∥
1
+
1
2
ρ
∥
b
1
+
H
u
-
d
1
∥
2
2
+
1
2
β
∥
b
2
+
W
u
-
d
2
∥
2
2
.Supposing thatD
,
α
,
d
1
,
d
2 are known, the solution with respect to u can be transformed to
(24)
min
⁡
u
τ
∑
i
j
‍
∥
R
i
j
u
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
+
1
2
ρ
∥
b
1
+
H
u
-
d
1
∥
2
2
+
1
2
β
∥
b
2
+
W
u
-
d
2
∥
2
2
.Obviously, it is a quadratic term. We can obtain the linear equation ofu by its closed form(25)
u
^
=
(
1
ρ
H
T
(
d
1
-
b
1
)
+
1
β
W
T
(
d
2
-
b
2
)
+
τ
∑
i
j
R
i
j
T
D
α
i
j
)
×
(
1
ρ
H
T
H
+
1
β
W
T
W
+
τ
∑
i
j
R
i
j
T
R
i
j
)
-
1
.The value of term about variableD is large, and other elements are so small that can be ignored. Thus, u
^ is about the same as the following form:
(26)
u
^
≈
∑
i
j
R
i
j
T
D
α
i
j
∑
i
j
R
i
j
T
R
i
j
.In [16], it has proved that learning dictionary using K-SVD can deal with the additive noise. However, it cannot deal with the Poisson noise directly. Here, we add the special constraint d
3
=
u to address this difficult problem. The special constraint comes from the idea which combines the split Bregman algorithm with an alternating minimization algorithm; refer to [26]. It is shown that this special alternating optimization algorithm simplifies the complex problem. Moreover, it solves the actual problem that K-SVD algorithm is only performed on additive noise denoising or multiplicative noise denoising, and multiplicative noise can be transformed to the additive noise in logarithm domain [17].


## 3.3. Proposed Algorithm

The difficulties in seeking variables are nonquadratic data term and nonseparablel
1-norm term. The alternating optimization algorithm can be used to solve these difficulties by introducing auxiliary variables d
1
=
H
u
,
d
2
=
W
u
,
d
3
=
u. The minimization problem becomes
(27)
min
⁡
α
,
D
,
d
1
,
d
2
,
d
3
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
λ
∥
d
2
∥
1
+
τ
∑
i
j
‍
∥
R
i
j
d
3
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
s
.
t
.
d
1
=
H
u
,
d
2
=
W
u
,
d
3
=
u
.Here, the alternating optimization algorithm is a special split Bregman algorithm. It is to be noted thatd
3 replaced u in the sparse representation of the image over a certain unknown dictionary.The proposed model is finally defined as(28)
min
⁡
α
,
D
,
u
,
d
1
,
d
2
,
d
3
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
λ
∥
d
2
∥
1
+
τ
∑
i
j
‍
∥
R
i
j
d
3
-
D
α
i
j
∥
2
2
+
τ
∑
i
j
‍
μ
i
j
∥
α
i
j
∥
1
+
1
2
ρ
∥
b
1
+
H
u
-
d
1
∥
2
2
+
1
2
β
∥
b
2
+
W
u
-
d
2
∥
2
2
+
1
2
γ
∥
b
3
+
u
-
d
3
∥
2
2
,
where ρ
,
β
,
γ are the positive penalty parameters. The minimization problem (28) can be decoupled into several subproblems via split Bregman iteration.(1) Givend
1
(
k
)
,
d
2
(
k
)
,
d
3
(
k
)
,
b
1
(
k
)
,
b
2
(
k
)
,
b
3
(
k
), the minimization of (28) with respect to u is transformed to
(29)
min
⁡
u
1
2
ρ
∥
b
1
(
k
)
+
H
u
-
d
1
(
k
)
∥
2
2
+
1
2
β
∥
b
2
(
k
)
+
W
u
-
d
2
(
k
)
∥
2
2
+
1
2
γ
∥
b
3
(
k
)
+
u
-
d
3
(
k
)
∥
2
2
.According to its closed form, we can obtain the linear equation ofu
(
k
+
1
)
(30)
u
(
k
+
1
)
=
(
1
ρ
H
T
(
d
1
(
k
)
-
b
1
(
k
)
)
+
1
β
W
T
(
d
2
(
k
)
-
b
2
(
k
)
)
M
+
1
γ
(
d
3
(
k
)
-
b
3
(
k
)
)
)
×
(
1
ρ
H
T
H
+
1
β
W
T
W
+
1
γ
)
-
1
.
The convolution of H can be easily performed in FFT domain.(2) Givenb
1
(
k
)
,
u
(
k
+
1
), the subproblem with respect to d
1 becomes
(31)
min
⁡
d
1
〈
1
,
d
1
-
g
log
⁡
d
1
〉
+
1
2
ρ
∥
b
1
(
k
)
+
H
u
(
k
+
1
)
-
d
1
∥
2
2
.Thend
1 can be obtained via solving the corresponding Euler-Lagrange equation of functional (31)
(32)
d
1
(
k
+
1
)
=
1
2
(
(
b
1
(
k
)
+
H
u
(
k
+
1
)
-
ρ
)
2
+
4
ρ
g
b
1
(
k
)
+
H
u
(
k
+
1
)
-
ρ
m
m
+
(
b
1
(
k
)
+
H
u
(
k
+
1
)
-
ρ
)
2
+
4
ρ
g
)
.(3) Givenb
2
(
k
)
,
u
(
k
+
1
), the third subproblem leads to a l
1-norm problem
(33)
min
⁡
d
2
λ
∥
d
2
∥
1
+
1
2
β
∥
b
2
(
k
)
+
W
u
(
k
+
1
)
-
d
2
∥
2
2
.The shrinkage operator is used to solve problem (31). Supposed that S
λ
(
t
)
=
sign
⁡
(
t
)
(
|
t
|
-
λ
)
+, d
2
(
k
+
1
)
=
S
λ
β
(
W
u
(
k
+
1
)
+
b
2
(
k
)
)
+, d
2
(
k
+
1
) is expanded,
(34)
d
2
(
k
+
1
)
=
max
⁡
{
∥
W
u
(
k
+
1
)
+
b
2
(
k
)
∥
1
-
λ
β
,
0
}
W
u
(
k
+
1
)
+
b
2
(
k
)
∥
W
u
(
k
+
1
)
+
b
2
(
k
)
∥
2
2
.(4) Givenu
(
k
+
1
)
,
d
3
(
k
)
,
D, the sparse representation coefficients {
α
i
j
} can be obtained by solving
(35)
min
⁡
α
i
j
∑
i
j
∥
R
i
j
d
3
(
k
)
-
D
α
i
j
∥
2
2
+
∑
i
j
μ
i
j
∥
α
i
j
∥
1
.Comparing (35) with (6), it can be found that (35) is the sparse coding stage of K-SVD algorithm. Any pursuit algorithm can be used to compute the representation vectors α
^
i
j.(5) Givenα
^
i
j
,
D
,
u
(
k
+
1
)
,
b
3
(
k
), the last subproblem with respect to d
3 is
(36)
min
⁡
d
3
τ
∑
i
j
‍
∥
R
i
j
d
3
-
D
α
^
i
j
∥
2
2
+
1
2
γ
∥
b
3
(
k
)
+
u
(
k
+
1
)
-
d
3
∥
2
2
.This is also a linear equation problem; similar to (33), we can obtain
(37)
d
3
(
k
+
1
)
=
(
1
γ
τ
I
+
2
∑
i
j
‍
R
i
j
T
R
i
j
)
-
1
×
(
1
γ
τ
(
b
3
(
k
)
+
u
(
k
+
1
)
)
+
2
∑
i
j
‍
R
i
j
T
D
α
^
i
j
)
.(6) Givenb
1
(
k
)
,
b
2
(
k
)
,
b
3
(
k
)
,
d
1
(
k
+
1
)
,
d
2
(
k
+
1
)
,
d
3
(
k
+
1
)
,
u
(
k
+
1
), the Bregman auxiliary variables b
1
,
b
2
,
b
3 are updated as follows:
(38)
b
1
(
k
+
1
)
=
b
1
(
k
)
+
H
u
(
k
+
1
)
-
d
1
(
k
+
1
)
,
b
2
(
k
+
1
)
=
b
2
(
k
)
+
W
u
(
k
+
1
)
-
d
2
(
k
+
1
)
,
b
3
(
k
+
1
)
=
b
3
(
k
)
+
u
(
k
+
1
)
-
d
3
(
k
+
1
)
.Besides, the dictionary learning is also an important procedure on sparse and redundant representation. We update the dictionaryD at the last step. The D is updated in an outer loop for saving computational time as suggested in [21]. Therefore, the loop times are divided into inner loop for main process and outside loop only for dictionary updating.Taking into account the above equations, the proposed Poissonian image deconvolution algorithm is summarized as shown in Algorithm1.Algorithm 1:Poissonian image deconvolution via sparse and redundant representations and framelet regularization.
Initialize: Initialization: u
0
=
g, d
1
0
=
H
g, d
2
0
=
W
g, d
3
0
=
g, D
=
DCT, and
b
1
0
=
b
2
0
=
b
3
0. Set ite
r
O
:
=
0, ite
r
I
:
=
0.
for  ite
r
O
=
0
to
N
O  do
while  stopping criterions is not satisfied  do
Step 1. Update u
(
k
+
1
) using (30);
Step 2. Update d
1
(
k
+
1
) using (32);
Step 3. Update d
2
(
k
+
1
) using (34);
Step 4. Update d
3
(
k
+
1
) using (37);
Step 5. Update sparse representation coefficients {
α
i
j
} using (6);
Step 6. Update the Bregman auxiliary variables b
1
(
k
+
1
)
,
b
2
(
k
+
1
)
,
b
3
(
k
+
1
);
ite
r
I
=
ite
r
I
+
1;
end
Step 7. Update Dictionary D;
endOne of the stopping criterions is calculated for each iteration by the “normalized step difference energy” (NSDE): NSDE
=
∥
u
(
k
)
-
u
(
k
-
1
)
∥
2
2
/
∥
u
(
k
)
∥
2
2, where u
(
k
) and u
(
k
-
1
) denote the image vector at the kth and (
k
-
1
)th iteration, respectively. The other stopping criterion is k
≤
N
I, where N
I is the maximum inner loop times. Each of the stopping criterion is satisfied; the inner loop will be terminated. In our paper, we set NSDE
≤
0.00001, N
I
=
40, and N
O
=
10. Besides, some conditions should be satisfied to obtain a sound deconvolution result. For example, the restored image u should be normalized, and it should be nonnegative in the process of iteration as we set u
=
max
⁡
⁡
(
u
,
0
).The proposed method involves many penalty parameters. Parameters choosing plays an important role in efficiency and stability of the algorithm. There are some skills to determinate these parameters. (1) In many literatures,ρ
,
β
,
γ can be equal to each other. To eliminate numerical problems, the three parameters should not be too small. (2) The value of ρ is the larger the better in practice. In our paper, we set ρ
∈
[
400,10000
]. We can see that λ
β is an important parameters group in the process of framelet regularization. The value of group λ
β belongs to empirical range [
3,10
]. The last parameters group is τ
γ which belongs to empirical range [
1,5
].


## 4. Experimental Results

In this section, in order to evaluate the performance of the proposed Poissonian image deconvolution method, we test and compare it with framelet regularization based method for Poissonian image deconvolution using Bregman iteration (PIDSB-FA) [13] method and dictionary learning approach and total variation regularization (TV-LD) [21] method. The ratio of peak signal and noise (PSNR) is frequently used to measure the restored image quality:
(39)
PSNR
=
10
log
⁡
10
(
L
2
NMSE
)
,
where L denotes the maximum value of pixel in an image, NMSE
=
∥
u
-
u
^
∥
2
2
/
∥
u
∥
2
2, and u and u
^ denote the original image and restored image, respectively.The numerical experiments are implemented with Matlab 7.14.0 on a dual core personal computer, 3.10 GHz processor and 2 G physical memory.Example 1.
To test the performance of the proposed method, numerical experiments are conducted on two kinds of test images. The first kind of test images is medical images [21], brain image with size of 380
×
366 and neck image with size of 300
×
300. The other kind of test images is nature images, house image with size of 300
×
300, Lena image with size of 521
×
521, and Bara image with size of 512
×
512 (Figure 1). All images are normalized in the interval [
0,1
].Original images. (a) brain, (b) neck, (c) Lena, (d) Bara, (e) house.
(a)
(b)
(c)
(d)
(e)These images are clear and original images. To simulate the Poissonian image deconvolution, these images should be degraded by the blur kernel and contaminated by Poisson noise. Gaussian blur kernel and uniform blur kernel are used to test the performance of deblurring. We choose9
×
9 Gaussian blur kernel with standard deviation 1 and 5
×
5 uniform blur kernel. Three Poisson noise levels are considered, I
max
⁡
=
600
,
1000
,
2000, where I
max
⁡ denotes the peak intensity of Poisson noise. The lower peak intensity corresponds to the stronger Poisson noise. In general, the simulated process has three steps. At first, the original image is scaled by I
max
⁡
/
U
max
⁡, where U
max
⁡ denotes the maximal value of the image. Subsequently, blurring process is added by convolving the normalized image with blur kernel. At last, Poisson noise is added to the burred image.The range of the parameters is given at the end of Section3. In this section, we offer the detailed parameters setting. The value of parameters pair λ
β is 10 for all blur kernels and peak intensities of Poisson noise. We set τ
γ
=
3 for Gaussian blur kernel, ρ
=
3000
,
ρ
=
8000
,
ρ
=
10000 for Gaussian blur kernel, and Poisson noise with I
max
⁡
=
600
,
1000
,
2000, respectively. We take ρ
=
10000
,
τ
γ
=
1, for uniform blur kernel and all peak intensities of Poisson noise.In all experiments, we set the parameters of K-SVD [16] as follows: iteration number is 30, the average error passing the threshold in OMP is ϵ
=
C
δ, C
=
1.15, δ
=
5, the size of dictionaries is 16
×
256, and the size of image patches is 4
×
4. Large patch size will produce redundant information and small patch size will lead to useful information missing. Besides, the larger image patch size is, the longer time will be taken [27]. The final dictionary learned at the last iteration of the proposed method for Example 2 is shown in (Figure 3).To illustrate the influence of the parameterρ on the performance of the proposed method, we note down the PSNR value of restored image when each parameter ρ is changing in (Figure 2). The plot tests and verifies the fact that parameter 1
/
ρ is not too high. Besides, the parameter receives good performance in a wide range.Figure 2
PSNR versus the value of parameter of the proposed method on Brain image contaminated by the Gaussian blur and Poisson noise withI
max
⁡
=
600.Figure 3
Dictionary learned by the proposed method for Example2.The parameters setting for PIDSB-FA are offered as follows:γ
λ
=
10, λ
=
0.025, and NSDE
=
1
0
-
4, and we set parameters for LDTV as suggested in [21].Restoration of brain is shown in Figure4. The brain image is blurred by uniform kernel of size 5
×
5 and contaminated by Poisson noise with I
max
⁡
=
600. Restoration of house is shown in Figure 5. The house image is blurred by Gaussian kernel of size 9
×
9 and contaminated by Poisson noise with I
max
⁡
=
2000. The image restored by PIDSB-FA method contains much details, no matter useful or useless. LD-TV reduces the useless details; however, it also smooths the useful details. The proposed method gives the sound restored results. It overcomes the problems which are caused by PIDSB-FA and LD-TV methods.Restoration of the brain image. (a) Degraded image is blurred by uniform kernel of size5
×
5 and contaminated by Poisson noise with I
max
⁡
=
600 (PSNR = 23.20 dB). (b) Restored image by the PIDSB-FA method (PSNR = 28.82 dB). (c) Restored image by the LD-TV method (PSNR = 28.08 dB). (d) Restored image by the proposed method (PSNR = 29.51 dB).
(a)
(b)
(c)
(d)Restoration of the Bara image. (a) Degraded house image is blurred by Gaussian kernel of size9
×
9 and contaminated by Poisson noise with I
max
⁡
=
2000 (PSNR = 25.14 dB). (b) Restored image by the PIDSB-FA method (PSNR = 25.68 dB). (c) Restored image by the LD-TV method (PSNR = 25.55 dB). (d) Restored image by the proposed method (PSNR = 26.29 dB).
(a)
(b)
(c)
(d)The PSNR values of degraded images and restored images by the three methods are compared in Table1. The improvement in terms of PSNR values obtained by PIDSB-FA is about the same as LD-TV. Experiment results describe that both analysis operator using piecewise linear B-splines framelet filter and synthesis operator using dictionary learning are good operators for image restoration. The proposed method combines the two operators. It can be found from Table 1 that the proposed method has achieved the highest PSNR on medical and nature images. Obviously, it remains the advantages of two operators.Table 1
PSNR (dB) values of three methods.
Image
I
max
⁡
Type
Degraded image
PIDSB-FA [13]
LD-TV [21]
Proposed method
Brain
600
Gaussian
27.06
32.19
32.73
32.80
Uniform
23.20
28.82
28.98
29.51
1000
Gaussian
27.48
33.33
33.67
33.94
Uniform
23.38
29.27
29.54
29.87
2000
Gaussian
27.84
34.33
34.35
35.07
Uniform
23.52
29.67
29.83
31.36
Neck
600
Gaussian
26.94
30.19
30.34
30.37
Uniform
23.62
27.21
27.35
27.50
1000
Gaussian
27.43
33.33
33.67
33.94
Uniform
23.82
27.43
27.48
28.12
2000
Gaussian
27.81
31.26
31.26
32.19
Uniform
24.00
27.57
27.66
29.19
Lena
600
Gaussian
28.84
32.89
32.81
33.45
Uniform
27.05
31.30
31.07
32.04
1000
Gaussian
30.02
33.76
33.63
34.08
Uniform
27.80
31.64
31.47
31.77
2000
Gaussian
31.18
34.34
34.30
34.88
Uniform
28.46
31.85
31.77
32.55
Bara
600
Gaussian
24.46
25.26
25.30
25.55
Uniform
22.65
24.50
24.41
24.54
1000
Gaussian
24.84
25.37
25.30
25.66
Uniform
22.89
24.56
24.49
24.57
2000
Gaussian
25.15
25.68
25.55
26.29
Uniform
23.09
24.59
24.57
25.04
House
600
Gaussian
27.98
32.10
32.19
32.92
Uniform
25.88
31.08
31.13
32.04
1000
Gaussian
29.02
32.73
33.00
33.66
Uniform
26.49
31.35
31.63
32.35
2000
Gaussian
30.02
33.34
33.68
34.53
Uniform
27.03
31.53
31.99
33.09Example 2.
The restoration of the Lena image is blurred by customized blur ([0 1 1 1 1 0 0; 1 1 0 0 0 1 0; 0 0 0 0 1 1 0; 0 0 1 1 0 0 0; 0 1 1 0 0 1 1; 0 1 0 0 1 1 0; 0 1 1 1 1 0 0]/22) and then contaminated by Poisson noise withI
max
⁡
=
2000. The parameters setting for proposed method are ρ
=
10000
,
λ
β
=
3, and τ
γ
=
1. Restoration of the Lena image by three methods is described in Figure 6. Lena image is a widely used image for its rich repetition pattern like texture region and strong sparse representation like large flat region and smooth region. The restored image by PIDSB-FA method preserves the rich detains in texture region such as decorations of the hat and also draws the artifacts into flat region such as Lena’s cheek. The restored image by LD-TV is relatively more smooth in flat region; however, it also misses some fine details in texture region. Meanwhile, the proposed method performs well on highlighting the texture region and smoothing the flat region. On the other hand, the proposed method receives better PSNR values than the other two methods.Restoration of the Lena image. (a) Blurred by using the customized kernel and then corrupted by Poisson noise (PSNR = 26.8674 dB). (b) Restored image by the PIDSB-FA method (PSNR = 30.62 dB). (c) Restored image by the LD-TV method (PSNR = 30.66 dB). (d) Restored image by the proposed method (31.74 dB).
(a)
(b)
(c)
(d)The above experiments describe that the proposed method gives better restoration of various images under different kinds of blurring kernels and Poisson noise with multilevel peak intensities. Besides, the stability is also the chief concerned matter in the proposed method. Figure7 shows that the NSDE and PSNR value of proposed method is changing with iterations on the brain image with Gaussian blur kernel and I
max
⁡
=
2000 of Poisson noise. The plot of NSDE describes the convergence, and PSNR value of restored image increases progressively until being stable.Figure 7
NSDE and PSNR values versus the number of iterations of the proposed method.


## 5. Conclusions

In this study, we proposed and implemented a Poissonian deconvolution method via sparse and redundant representations and framelet. The proposed method includes three terms: the Poisson data term, the synthesis term, and the analysis term. Dictionary learning for synthesis model and B-spline framelet-based analysis prior are integrated perfectly. The proposed method keeps the satisfied result of dictionary learning method; meanwhile, it takes full advantage of the multiscale characteristic of framelet-based methods. The seemingly unsolved variable optimization problem is efficiently addressed by using the specific split Bregman method. Furthermore, parameters are divided into groups which simplify the parameter setting greatly. Numerical comparisons on a various of simulated images show that the proposed method outperforms the other two latest methods. However, the dictionary learning is a time-consuming process. Therefore, improving the proposed method on dictionary learning in terms of computation time is the future research.



---

*Source: 917040-2014-01-16.xml*