Wavelet frame based Poisson noise removal and image deblurring

Accepted Manuscript

Wavelet Frame Based Poisson Noise Removal and Image Deblurring Haimiao Zhang, Yichuan Dong, Qibin Fan PII: DOI: Reference:

S0165-1684(17)30034-8 10.1016/j.sigpro.2017.01.025 SIGPRO 6384

To appear in:

Signal Processing

Received date: Revised date: Accepted date:

30 June 2016 3 November 2016 23 January 2017

Please cite this article as: Haimiao Zhang, Yichuan Dong, Qibin Fan, Wavelet Frame Based Poisson Noise Removal and Image Deblurring, Signal Processing (2017), doi: 10.1016/j.sigpro.2017.01.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A nonconvex and non-continuous regularized variational model is proposed

CR IP T

to recover the Poisson noisy and blurred image. • The proposed nonconvex model improves the quality of restored Poisson

noisy and blurred image compared with the convex models such as total variation based model and L 1 tight wavelet frame based model.

• The proposed algorithm is efficient to implement and produces more ac-

AN US

curate solution.

• Convergence and stability of the proposed algorithms are verified from the

AC

CE

PT

ED

M

numerical perspective.

1

ACCEPTED MANUSCRIPT

Wavelet Frame Based Poisson Noise Removal and Image Deblurring✩

a School

CR IP T

Haimiao Zhanga , Yichuan Donga , Qibin Fana,∗ of Mathematics and Statistics, Wuhan University, Wuhan, Hubei, 430072, P.R. China

Abstract

AN US

The recovery of sparse signal from noisy data arises in various application fields.

One widely known approach for Gaussian noise image restoration with wavelet frame based sparse representation is the l0 norm regularized variational model. In this paper, the sparse and nonconvex noncontinuous l0 norm regularized model is proposed to recover the Poisson noise and blurred image. Then the resulted optimization problem is solved by the alternating direction method of

M

multipliers(ADMM) scheme and two different approaches are adopted to solve the ADMM scheme in the numerical experiments. Extensive simulation results

ED

verify the convergence of the proposed algorithm and indicate that the proposed l0 norm based nonconvex model is efficient and comparable with some state-ofthe-art approaches.

PT

Keywords: ADMM, Iterative hard thresholding(IHT), Conjugate Gradient(CG) method, Nonconvex regularization, Poisson noise removal

CE

2010 MSC: 65T60, 68U10, 90C90

AC

✩ This research was supported by the National Science Foundation of China under Grant 61179039, by the National Key Basic Research Development Program(973 Program) of China under Grant 2011CB707100, and by the Fundamental Research Funds for the Central Universities under Grant 2015201020202. ∗ Corresponding author Email addresses: [email protected] (Haimiao Zhang), [email protected] (Yichuan Dong), [email protected] (Qibin Fan)

Preprint submitted to Elsevier

February 3, 2017

ACCEPTED MANUSCRIPT

1. Introduction The recovery of Poisson noise corrupted and blurred images arises in various application fields such as astronomical imaging, fluorescence microscopy and

CR IP T

positron emission tomography [1]. The challenge is to design efficient numerical algorithm and to analyse the convergency property. Let u ∈ Rn×n be the

original image, A be the linear blurring operator and f be the recorded image. Then the Poisson noisy and blurred observation of an image can be expressed as f = P(Au),

AN US

(1)

where P stands for the Poisson noisy process.The likelihood of the observed image f is given by p(f |Au) =

N Y ((Au)i )fi

fi !

i=1

e−(Au)i ,

(2)

M

where (Au)i is the ith component of Au, N = n2 is the total number of pixel. Moreover, in order to make it more convenient in the following discussion, we will claim that each pixel fi is finite in the image domain and fi > 0.

ED

Instead of maximizing p(f |Au), the original clear image can be obtained by minimizing − log p(f |Au) with respect to u from the observed data f . However, this ill-posed inverse problem can only return meaningful solution by incorpo-

PT

rating some suitable priors, e.g., the sparsity representation with wavelet basis [2, 3] and the bounded total variation(TV) [4, 5]. In this paper, we propose

CE

to use the tight wavelet frame [3] to sparsely representing the image data. As a result, we can recover the Poisson noisy and blurred image by solving the

AC

following regularized variational optimization problem

5

minh1, Aui − hf, log(Au)i + λφ(u), u≥0

(3)

where h·, ·i is the inner product and λ > 0 is the regularization parameter. The

first two terms are the data fidelity term for Poisson noise, and φ(u) is the

regularization term for a priori information.

3

ACCEPTED MANUSCRIPT

There are many different choices of the regularization term φ(u) in the literature. For example, in [6, 7, 8, 9], the authors used φ(u) = T V (u), where 10

T V (u) is the total variation of u, in the Poisson noise removal model. The total

CR IP T

generalized variation (TGV) was adopted in [5] to the Poisson denoising and eliminating the staircase effect cased by the TV model. Recently, the high-order TV model with φ(u) = α||∇u||1 + β||∇2 u||1 was proposed in [10] to overcome the shortcoming of stair-case effect of TV model as well. In [11] the authors 15

chose a fractional-order and non-local composite TV regularization term to process the cartoon and texture components in the degraded Poisson noisy image

AN US

recovery. Moreover, the authors in [12] proposed a PDE-based model for outof-focus blur, Poisson noise and low axial resolution image restoration. The variance stabilizing transform based approaches[13, 14] are another routine for 20

Poisson noise removal. In this paper, we mainly consider the variational regularization models and propose a wavelet frame based model for Poisson noise image restoration.

M

Since the TV based model often causes staircase artifacts and over-smoothing at the edges in the recovered image, wavelet frame based l1 models were studied

ED

in [15, 16] to remove the Poisson noise in the degraded image. In this work, we propose a novel variational model which adopts the regularization term λφ(u) = |||λ · W u|||0 that is nonsmooth nonconvex. The l0 norm regularized model has

PT

ever been discussed in the case of additive white Gaussian noise(AWGN) image restoration and compressive sensing [17, 18, 19, 20]. The tight wavelet frame

CE

based model has the advantage of reducing the staircase effect in TV based model and recovering image with more sharp edges. Here we extend it to restore

AC

the Poisson noisy and blurred image and obtain a special form of (3) as

where |||λ · α|||0 =

minh1, Aui − hf, log(Au)i + |||λ · W u|||0 , u≥0

PL−1 P l=0

j∈I

(4)

λl,j ||αl,j ||0 . What’s more, ||αl,j ||0 is defined to

be the number of nonzero elements of tight wavelet frame coefficient αl,j and

25

called the l0 “norm”, and λl,j > 0 is the weight of αl,j at a given pixel j in band l. Numerical results will indicate the improvement of this new model to restore 4

ACCEPTED MANUSCRIPT

the Poisson noisy and blurred image when compared with other exited model that has a l1 or TV regularization term. To overcome the challenging of nonconvex and noncontinuous optimization problem with l0 term, researchers commonly turn to replace it with a convex

CR IP T

30

regularization term to approximate φ(u). In order to utilize the sparse wavelet representation of the natural image, a convex model with regularization term φ(u) = ||W u||1 is discussed in [15] for the Poisson noise image deblurring where W is the tight wavelet frame transform. Then an efficient algorithm based on 35

the alternating direction method of multipliers(ADMM) scheme is proposed to

AN US

solve the corresponding convex optimization problem.

However, the nonconvex noncontinuous l0 model is a more suitable choice for sparsity prior to model the objective image [21, 22]. The motivation of our newly proposed nonconvex l0 model for Poisson noise removal is to overcome the 40

staircase effect and oversmooth that occurred in the TV and l1 based models and at the same time better preserve the fine details in the recovered image.

M

Theoretical and numerical results in [23, 24, 25] also verify the effectiveness of nonconvex model for image restoration.

45

ED

There are a large number of literatures for Poisson noise removal. The variational model with exact data fidelity term and regularization term gone with the same routine as we discussed here can be classified into three main

PT

categories: the dictionary learning based model[26], TV regularization model[6, 7, 8, 9, 10], wavelet frame based model[15]. Comparison with other strategies

CE

for Poisson noise removal and the mixed noise image restoration model will be 50

the further research. In this work, we mainly make comparison with the several most recently proposed TV based model and wavelet frame based model to

AC

recover the Poisson noisy and blurred image. In the following sections, the l0 optimization problem (4) is discussed and

iteration scheme of the algorithm is described. There are many approaches to

55

solve the l0 minimization problem such as the fixed-point proximity algorithm in [19, 20] and the reduced dimension Frank-Wolfe algorithm in [27]. The iterative hard thresholding(IHT) algorithm [28] is simple to implement in the numerical 5

ACCEPTED MANUSCRIPT

experiment for l0 minimization problem, so we adopt it to deal with the nonconvex subproblem involved in the ADMM scheme. Two different strategies are 60

proposed to solve the ADMM based algorithm: the first one is to adopt the

CR IP T

IHT algorithm one time in each iteration which is implemented efficiently and the second one is to use the conjugate gradient(CG) algorithm [29, Chapter 5]

which produces more accurate solutions. Even though there exists convergence

results in literature, e.g.,[30, 31, 32], for ADMM scheme, they can not be directly 65

applied to the proposed nonconvex and noncontinuous regularized variational model (4). The global convergence of nonsmooth or nonconvex optimization

AN US

problem in[33, 34, 35] motivated us to analysis the novelty approaches proposed here. By the convex analysis tools, we can verify the convergence analysis results for the proposed numerical algorithm. As plenty of complex technique is 70

involved in the theoretical analysis, we only give numerical simulation about the converge of our algorithms in Section 3. More detailed investigation about the convergence analysis of the proposed algorithms can be found in [36].

M

The paper is organized as follows. We start by giving the solution of the proposed nonconvex noncontinuous l0 norm model for Poisson denoising and deblurring in Section 2. To illustrate the denoising and deblurring performance

ED

75

of the proposed algorithms, extensive numerical experiments are established in Section 3 and compared them with other recently proposed models. Finally, the

PT

conclusions are presented in Section 4.

CE

2. Solutions and Algorithms The ADMM algorithm is widely used to solve large scale linearly constrained

80

optimization problems arising from the application fields such as image process-

AC

ing [37], engineering computation and machine learning[38]. In this section, we adopt the ADMM scheme to solve the nonconvex noncontinuous Poisson image denoising and deblurring problem. For convenience, we now rewrite the tight framelet based Poisson noisy image

6

ACCEPTED MANUSCRIPT

model (4) equivalently as min µ(h1, Aui − hf, log(Au)i) + |||λ0 · Wu|||0 , u≥0

(5)

CR IP T

where µ > 0 and λ0 = λµ. Through introducing two new auxiliary variables, the above l0 minimization problem can be reformulated into a constrained optimization problem as

min{µ(h1, wi − hf, log(w)i) + |||λ0 · Wu|||0 + δRN | Au = w, u = z}, + (z) u

85

(6)

is the indicator function of RN where δRN + . The corresponding augmented La+

AN US

grangian functional of this constrained optimization problem is

L(w, u, z; p) = µ(h1, wi − hf, log wi)+ |||λ0 · Wu|||0 + ιRN (z) + γ η + ||Au − w + γ −1 p1 ||22 + ||u − z + η −1 p2 ||22 , 2 2

(7)

where p = (p1 , p2 ) is the Lagrangian multiplier, γ > 0 and η > 0 are penalty parameters. The iterative scheme of ADMM for (7) is

γ ||Auk − w + γ −1 pk1 ||22 , 2 γ η = arg min ||Au − wk+1 + γ −1 pk1 ||22 + ||u − z k + η −1 pk2 ||22 + |||λ0 · Wu|||0 , u 2 2 η = arg min δR+ (z) + ||uk+1 − z + η −1 pk2 ||22 , z 2 = pk1 + γ(Auk+1 − wk+1 ),

uk+1 z k+1 pk+1 1

ED

w>0

M

wk+1 = arg min µ(h1, wi − hf, log wi)+

pk+1 = pk2 + η(uk+1 − z k+1 ), 2

PT

(8)

where pk = (pk1 , pk2 ) is the Lagrangian multiplier.

CE

We now investigate these subproblems in (8) one by one for the Poisson noise

AC

image deblurring. The w-subproblem has a closed form solution q 1 wk+1 = [γAuk + pk1 − µ + (γAuk + pk1 − µ)2 + 4γµf ]. 2γ

(9)

The solution for z-subproblem in (8) can also be explicitly given as z k+1 = max{uk+1 + η −1 pk2 , 0}.

(10)

Note that the nonconvex and noncontinuous l0 minimization problem for u in the ADMM scheme (8) needs more effort to dispose. To simplify the 7

ACCEPTED MANUSCRIPT

presentation, we denote ξ1k = wk+1 − γ −1 pk1 and ξ2k = z k − η −1 pk2 . Furthermore, we omit the index k when clear from the context. The quadratic penalty function for u is

0

0

0

1 η ρ 1 ||Au − ξ1 ||22 + ||u − ξ2 ||22 + ||Wu − v||22 + |||λ0 · v|||0 , 2 2 2 γ

CR IP T

p(u, v) ,

0

where η = γη , ρ = γρ . In the following, we consider the solution of the nonconvex optimization problem min p(u, v).

(11)

u,v

By using the penalty decomposition(PD) algorithm [21] to (11), there are

90

AN US

mainly two alternative iteration steps

0

0

η ρ 1 ul+1 = arg min ||Au − ξ1 ||22 + ||u − ξ2 ||22 + ||Wu − vl ||22 , u 2 2 2 0 1 ρ vl+1 = arg min ||v − Wul+1 ||22 + |||λ0 · v|||0 . v 2 γ

(12) (13)

M

For (12) and (13) we have the closed form solution

(14)

vl+1 = Hλ˜ (Wul+1 ),

(15)

0

where θ = η + ρ , I is an identity matrix, and Hλ˜ (x) is the pointwise hard q q ˜ = 2λ00 = 2λ0 . thresholding operator with parameter λ ρ γρ

PT

95

0

ED

0

0

ul+1 = (A> A + θI)−1 (A> ξ1 + η ξ2 + ρ W> vl ),

Therefore, combining all the formulae (9), (10), (14) and (15), the alterna-

tive iteration scheme is described in Algorithm 1. The choosing of stopping

CE

criterion parameters, such as the maximum iteration(M axIter), the relative

error(RelErr) and the tolerance (tol) will be clarified in Section 3. Under the periodic boundary condition for u, the blurry operator A and its

transpose are all block circulant[39]. The closed form solution for u-subproblem

AC 100

can be solved by 2D discrete fast Fourier transform(FFT). So the resulted algorithm is highly efficient in the numerical implementation. In the following sections, the above Algorithm 1 is denoted by the FT-ADMM algorithm to indicate that the u-subproblem is solved by FFT.

8

ACCEPTED MANUSCRIPT

Algorithm 1 FT-ADMM for Poisson noisy image deblurring 1: Initialization: Set w 0 = 0, u0 = 0, z 0 = 0, p01 = 0, p02 = 0. Choose the

2:

Main loop:

3:

while k ≤ M axIter and RelErr > tol do

q 1 [γAuk + pk1 − µ + (γAuk + pk1 − µ)2 + 4γµf ], 2γ   0 0 = (A> A + θI)−1 A> ξ1 + η ξ2 + ρ W> Hλ˜ (Wuk ) ,

wk+1 = uk+1

CR IP T

parameters: γ > 0, µ > 0, θ > 0, η > 0, λ > 0, ρ > 0.

z k+1 = max{uk+1 + η −1 pk2 , 0},

AN US

pk+1 = pk1 + γ(Auk+1 − wk+1 ), 1 pk+1 = pk2 + η(uk+1 − z k+1 ). 2 4:

end while

5:

Output: w∗ , u∗ , z ∗ .

(16)

We also adopt another strategy to solve (12) and (13) and compare its per-

M

105

formance with the FT-ADMM algorithm. The block coordinate descent(BCD)

ED

method[21] is implemented, then the unconstrained quadratic programming problem (12) is solved by the conjugate gradient method and v-subproblem is solved by hard thresholding. This CG based approach is denoted by the CG-ADMM algorithm and described the details in Algorithm 2.

PT

110

The numerical experiments show that even though the CG-ADMM algorithm is not very efficient in terms of the CPU time, it produces more accurate

CE

solution than FT-ADMM and other alternative models. More detailed investigation of both newly proposed algorithms can be found in Section 3. To deal with the convergence properties of the proposed ADMM scheme

AC

115

based algorithm, it requires to pay some attentions to that the objective function of the augmented Lagrangian functional (7) which has a nonconvex and nonsmooth regularization term with l0 norm. Furthermore, since the u-subproblem is nonconvex, the convergence results for ADMM with convex and smooth ob-

120

jective function can not be applied directly[40, 41] to our model. While the 9

ACCEPTED MANUSCRIPT

Algorithm 2 CG-ADMM for Poisson noisy image deblurring 1: Initialization: Set w 0 = 0, u0 = 0, z 0 = 0, p01 = 0, p02 = 0. Choose the

2:

Main loop:

3:

while k ≤ M axIter and RelErr > tol do wk+1 =

4:

1 [γAuk + pk1 − µ + 2γ

CR IP T

parameters: γ > 0, µ > 0, θ > 0, η > 0, λ > 0, ρ > 0.

q (γAuk + pk1 − µ)2 + 4γµf ],

while l <= lmax do u-subproblem(12) is solved by CG algorithm.

6:

v-subproblem(13) is solved by hard thresholding.

7:

AN US

5:

end while

uk+1 = uklmax ,

z k+1 = max{uk+1 + η −1 pk2 , 0},

pk+1 = pk1 + γ(Auk+1 − wk+1 ), 1

M

pk+1 = pk2 + η(uk+1 − z k+1 ). 2 end while

9:

Output: w∗ , u∗ , z ∗ .

ED

8:

PT

convergence analysis for ADMM scheme with nonsmooth or nonconvex objective objective function can be found in [33, 42], they are again not suitable for the proposed model as well. As the objective function of the l0 minimization

CE

problem is nonconvex and noncontinuous, it is difficult to estimate the primal

125

and dual optimality gaps to further obtain the convergence results. Therefore, the existed analysis results for ADMM scheme in the literatures can not be

AC

applied in our model directly. However, with the decreasing property of the augmented Lagrangian functional, it is possible to guarantee the convergency of the proposed FT-ADMM and CG-ADMM algorithms. In this paper, through

130

the numerical experiments we will show the convergence of our algorithm from the numerical perspective.

10

ACCEPTED MANUSCRIPT

3. Numerical Simulations In this section, we conduct several experiments to illustrate the performance of the proposed algorithms: FT-ADMM and CG-ADMM for Poisson noise and blurred image restoration. Then the convergence property of the ADMM scheme

CR IP T

135

based algorithm is verified by plenty of numerical simulation results. At last, we

present some comparison results obtained from other recently proposed models

such as the tight wavelet frame based l1 model in [15] and the TV based models

in [7, 10] for Poisson noisy image deblurring. Even though there are plenty of 140

models and algorithms have been developed in recent years to solve the Poisson

AN US

noisy image restoration, here we mainly focus on the regularization based vari-

ational models. More comparisons with different models from various routines will be the future work.

The experiments are performed on a PC with an Intel Core(TM) i7 − 5500U 145

CPU(2.40GHz) and 4GB of RAM running under Windows 7.

M

3.1. Parameter Setting and Notations

The quality of restored images with different methods is measured quanti-

ED

tatively by using the Signal-to-Noise Ratio(SNR), the relative error (RelErr) and structural similarity(SSIM)[43]. The SNR, RelErr and SSIM are defined as follows:

PT

150

||u||2 ||˜ u − u||2 , RelErr(u, u ˜) , and ||˜ u − u||2 ||u||2 (2µu µu˜ + C1 )(2σu˜u + C2 ) SSIM (u, u ˜), 2 , (µu + µ2u˜ + C1 )(σu2 + σu2˜ + C2 )

CE

SN R(˜ u),20 log10

(17) (18)

where u and u ˜ are the original image and the restored image respectively, µx is

AC

the average of x, σx is the variance of x and σxy is the covariance of x and y.

The constants C1 and C2 can be thought of as stabilizing constants for near-zero denominator values.

155

In the following numerical experiments, we terminate the iterations and ob-

tain the restored images when the relative error RelErr(uk , uk+1 ) =

||uk+1 −uk ||2 ||uk ||2

tol = 3×10−3 or the maximum outer iteration(M axIter = 35) is exceeded. The 11

<

ACCEPTED MANUSCRIPT

convergency and stability of the algorithm in iteration is indicated by the relative error decreasing tendency in the numerical experiments. While the relative 160

error approaches to a limit point such as the small tolerance(tol) with respect

CR IP T

to the iteration, the algorithm is recognized as convergence. And furthermore, if the performance of restoration procedure is not seriously influenced by the experimental setting and the noise in the data, we recognized it as stable. From now on, if we use RelErr in the context, it means RelErr(u, u ˜). Otherwise, we 165

also use RelErr(uk , uk+1 ) to represent the relative difference in the iteration.

In this paper, we tested a set of images, including the well-known Camera-

AN US

man, Lena, Pepper, Baboon, Blobs and Moon. Fig. 1 shows the original image and its size. The degraded Poisson noisy and blurred images are obtained in

three steps: 1) in order to create degraded image f with different levels of 170

Poisson noise, the original clear image u is scaled by a fixed number( such as Scale = 50, 100, 200, 400); 2) the scaled image is convoluted with a given point spread function to simulate the blurring effect of the operator A; 3) the blurred

M

image is contaminated by the Poisson noise to obtain the final degraded image f . The specific image restoration task we considered here is Poisson noise removal and Gaussian deblurring, as other types of blurring kernels can be processed

ED

175

in a similar way. Since the regularization parameter ρ, η and γ in(7) and (11) are not obviously affected the performance of the recovered images for different

PT

Poisson noise level and the optimal µ in (7) can be set to 0.8 · Scale. In the following experiments we will only consider the Scale = 200 case. To make the algorithm practicable and reproducible, the parameters in the

CE

180

proposed ADMM based algorithms, FT-ADMM and CG-ADMM, are fixed in the numerical experiments. We use the piecewise linear B-spline wavelet frame

AC

in the following experiments. From the numerical experiments we observe that the coarsest decomposition level of the objective image influences the perfor-

185

mance of the proposed algorithms. Therefore, it needs to be further discussed and more detailed investigation is put in the next subsection. From the numerical simulation results of FT-ADMM and CG-ADMM with plenty of test images, we observe that the proposed algorithm is not seriously influenced by the fixed 12

ACCEPTED MANUSCRIPT

parameters γ and ρ for different input degraded data. The penalty parameters 190

in (7) and (11) are empirically set to ρ = 0.08 and γ = 0.5. Since the penalty restored images, it is set to be η = 1 × 10−6 . 3.2. Efficiency of the ADMM scheme for the l0 model

CR IP T

parameter η for positiveness of the pixels is also not sensitive to the recovery of

In this section we consider the performance of the newly proposed FT195

ADMM and CG-ADMM algorithms. Both of the proposed algorithms are based on the sparse representation with wavelet frame coefficients, so the decomposi-

AN US

tion level of the objective image will influence the recovery quality in terms of accuracy(i.e., SNR, SSIM and RelErr) and efficiency (CPU time). This can be

observed in the following experiments. The convergency property of the ADMM 200

based algorithms can be observed from the numerical perspective. Both the initialization conditions for nonconvex optimization problem and the ADMM scheme are able to influence the performance of the algorithm. So we

M

choose u0 = 0 and u0 = f in the numerical experiments to test the performance with different images. The wavelet frame decomposition level is chosen: L = 4. Then the ”Cameraman” from Fig.1 is used to report the results obtained from

ED

205

FT-ADMM and CG-ADMM. Detailed results are established in Table 1. The measure of SNR is dB and the CPU time is in seconds. It is noticeable that

PT

FT-ADMM is about three times faster than CG-ADMM in CPU time with different initialization conditions. The difference of SNR is about 1dB in both 210

cases. It is obviously that the ADMM scheme based algorithms initialized by

CE

u0 = 0 produce solutions with higher SNR and lower RelErr than those by the degraded images. The FT-ADMM is more efficient with the initialization

AC

condition u0 = 0, but the SNR is lower than the one obtained by u0 = f . While the CG-ADMM can recover more accurate (lower RelErr and higher SSIM)

215

solution with u0 = 0, it is time consuming. From the results in Table 1, we can also observe that the ADMM based

algorithm initialized by u0 = 0 needs more outer iteration in the numerical implementation to obtain a stable solution than by u0 = f . The stability and 13

ACCEPTED MANUSCRIPT

convergency of the proposed algorithm is influenced by the initialization condi220

tion. This is also the common property of the nonconvex based imaging models. Even though the convergence analysis results for nonconvex optimization advise

CR IP T

us to choose u0 = f as a initialization condition, the numerical experiments show that the proposed algorithms produce a stable local minimizer with u0 = 0. In

the following numerical experiments, we will choose u0 = f for FT-ADMM and 225

u0 = 0 for CG-ADMM in order to balance the accuracy and efficiency.

We next compare the performance of both the proposed algorithms by choosing different wavelet frame decomposition level and report the result with the

AN US

test image ”Pepper”. The recovery results are summarized in Table 2. We can notice that the recovery quality is better with larger wavelet frame decomposi230

tion level. But as a result, the computational burden is also heavy. In order to balance the computational cost and the recovery quality, we will choose L = 2 in the following experiments for our proposed FT-ADMM and CG-ADMM algorithms.

235

M

Now, we consider the convergence property of the proposed ADMM based algorithms. The relative error curve of FT-ADMM and CG-ADMM are reported

ED

in Fig. 2 for two test images: Cameraman and Lena. From the results in Fig. 2 we observe that both ADMM based algorithms converge to a stable solution in the numerical simulation. The relative difference(RelErr(uk , uk+1 )) of solution

240

PT

sequences {uk } in the iteration decreases stably in FT-ADMM for different test images, while it is not stable in CG-ADMM. This shows us that FT-ADMM has

CE

better convergency property than the CG-ADMM in numerical practice. The ADMM scheme will converge to a stable solution even though the subproblem is not accurately solved in each iteration. From the numerical experiments we

AC

observe that the ADMM based algorithm initialized by non-zero data u0 = f will

245

converge to a local minimizer more stable and quickly. The result in Fig.2 is a positive answer to verify the convergence property of the proposed algorithms. This shows the advantage of using the ADMM scheme to solve our problem at hand. In next subsection, we will adopt the newly proposed ADMM based algorithms, FT-ADMM and CG-ADMM, to compare with other state-of-the-art 14

ACCEPTED MANUSCRIPT

250

approaches for Poisson noise image deblurring. 3.3. Comparison of the l0 , l1 and TV based models for PNID

CR IP T

In this subsection, we will compare the performance of the proposed tight wavelet frame based nonconvex l0 model with the existed convex models for

Poisson noise image deblurring. The tight wavelet frame based l1 norm regu255

larized model in [15] was compared with our l0 ”norm” model and the former model was solved by the generalized inverse with linearized alternating mini-

mization(GILAM) algorithm which is depended on the ADMM scheme. Both

AN US

of the two models are proposed to penalize the sparsity of the wavelet frame coefficients of objective image. The advantages of the newly proposed nonconvex 260

one will be obvious in the following numerical results. Another representative convex models for Poisson noise removal and deblurring are the TV semi-norm regularized one proposed in [7] and the high-order TV model in [10]. All these models are solved by the ADMM scheme for better convergence properties and

M

more convenient to design efficient numerical algorithms.

We will compare the proposed methods (FT-ADMM and CG-ADMM) with

265

ED

the GILAM algorithm proposed in [15], the PNIDAL algorithm for TV model proposed in [7] and the ADMPNR for high-order TV proposed in [10]. The stoping criterion for GILAM is set to be tol = 3 × 10−3 for fairly comparison 270

PT

and other parameters are chosen manually for optimal image reconstruction results. The parameter setting for TV models, e.g. PIDAL and ADMPNR, are selected based on the original literature and then manually tuned in order to

CE

obtain optimal recovered images. Numerical simulation results are reported in Table 3, where all of the quantities are the ten times average for each test image.

AC

In Table 3, our newly proposed algorithms, FT-ADMM and CG-ADMM, are

275

denoted by l0 -FT and l0 -CG respectively. The tight wavelet frame based convex model in [15] solved with GILAM is denoted by l1 and the TV based models solved by PIDAL and ADMPNR are denoted by T V and 4th-TV, respectively.

The SNR, relative error, average iteration, SSIM and CPU time are reported in Table 3 from column three to seven. 15

ACCEPTED MANUSCRIPT

From Table 3 we can observe that the tight wavelet frame based models are

280

better performed than the TV models. More specially, the l0 -FT, l0 -CG and l1 models produce better recovery results than T V and 4th-TV’s. Moreover, for

CR IP T

the tight wavelet frame based models, the l0 -CG can recover images with higher SNR and SSIM than l0 -FT and l1 . However, the drawback of l0 -CG is that 285

it is more time consuming than other compared approaches. We should pay attention to the fact that: for the test images Lena and Pepper, even though

the SNR of l0 -FT is higher than l1 , the SSIM of l0 -FT is lower. While in other

cases, the higher SNR the corresponding model produced, the larger SSIM one

290

AN US

obtained. This indicates that SNR is not the only choice for image quality assessment and at the same time, SSIM is a better measure as recommended in [43]. From the results in Table 3, we conclude that the l0 -FT is a good choice for practice since it is stable in the numerical simulation and balances the quality of recovered images and the computational efficiency compared with all the other wavelet frame based models and the TV based models.

The GILAM, PIDAL and ADMPNR algorithms and the proposed FT-ADMM

M

295

and CG-ADMM are employed to recover the Poisson noise and Gaussian blurred

ED

images. Since the recovered images have a similar performance for tested data sets, we only establish the test images: Lena, Baboon, Blobs and Moon in Fig.3, Fig.4, Fig.5, Fig.6, respectively. To have a better visualization of the restored results, the zoomed parts of the test images from l0 -CG, l1 and high-order TV

PT

300

are also established.

CE

From the results in Fig.3 to Fig.6, we can observe that our algorithms have better performance, e.g., more visible details and sharp edges than the compared approaches. The proposed l0 model produces recovered images with better

visual effect, but l1 model in [15] has the oversmooth effect which is obvious

AC

305

at the zoomed parts in Fig.3, Fig.4 and Fig.6. While compared with the TV models, the tight wavelet frame based models produce images with more sharp edges and better visual effect. Detailed information can be found from the recovered images in Fig.3, Fig.4, Fig.5, Fig.6.

310

In Fig.3 and Fig. 4, the texture part of recovered Lena and Baboon from 16

ACCEPTED MANUSCRIPT

the wavelet frame based l1 model are obviously over-smoothed compared with the results from our l0 model. In Fig.5, the recovered image by high-order TV does not have sharp edges in the zoomed part. In Fig.6, the zoomed part of

315

CR IP T

Moon is smoothed by l1 and high-order TV models while our l0 model keeps the detailed information well.

The FT-ADMM is computationally more efficient than CG-ADMM and

GILAM, while the CG-ADMM produces more accurate solution than FT-ADMM and GILAM. Even though the TV models solved by PIDAL and ADMPNR are

more efficient in CPU time, the recovered images from wavelet frame based

models are obviously being higher quality in which case the objective images

AN US

320

contain more texture. The Poisson noise is not thoroughly removed by the TV models. This can be verified by the results in Fig.5 where plenty of geometric features still containing noise. What’s more, from Fig. 6 we can observe that the TV model solved by PIDAL always produce piecewise constant image and 325

the detailed features were lost in the recovered image. Sharp edges and more

M

detailed feature are obviously recovered by the wavelet frame based model l0 and l1 in Fig.6. The tight wavelet frame based models have the ability to re-

ED

store images that contain both cartoon and texture with better performance. To balance the computational burden and restoration quality, the FT-ADMM 330

might be a nice choice.

PT

Based on the quality measurement results in Table 3 and the visualized images in Fig.3 to Fig.6, we recommend the FT-ADMM method for Poisson

CE

noise image deblurring in practice for the fact that this approach is a good balance of the computational cost and the recovery quality in terms of CPU time and SSIM(or SNR and RelErr), respectively.

AC

335

4. Conclusions In this work we have proposed to use the l0 norm of tight framelet repre-

sentation as the sparsity prior of objective image and a nonconvex regularized variational model is adopted to restore the image which is contaminated by

17

ACCEPTED MANUSCRIPT

340

the Poisson noise and blurring kernel. Within the ADMM scheme, two different strategies are used to solve the nonconvex l0 optimization problem. The convergence property of the proposed algorithm is verified from the numerical

CR IP T

perspective. Furthermore, the numerical experiments indicate that this nonconvex model 345

produces desirable solution for the restoration task at hand and the ADMM

based algorithms converge to a stable solution within a few iterations. This

nonconvex model is comparable with the state-of-the-art approaches such as the wavelet frame based l1 model and TV based models for Poisson noise re-

350

AN US

moval and deblurring which are widely studied in many existed literatures for

Gaussian noise signal/image restoration. Future work will apply this model to deal with the large scale parallel signal recovery in the field such as X-ray computed tomography image reconstruction and perform data driven parameter estimation rule in the numerical algorithm.

355

M

References

[1] M. Bertero, P. Boccacci, G. Desider` a, G. Vicidomini, Image deblurring with

ED

poisson data: from cells to galaxies, Inverse Probl. 25 (12) (2009) 123006. [2] S. Mallat, A wavelet tour of signal processing: the sparse way, Academic

PT

press, 2008.

[3] B. Dong, Z. Shen, et al., MRA based wavelet frames and applications, IAS Lecture Notes Series, Summer Program on The Mathematics of Image

CE

360

Processing, Park City Mathematics Institute 19.

AC

[4] L. I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1) (1992) 259–268.

[5] X.-d. Wang, X.-c. Feng, W.-w. Wang, W.-j. Zhang, Iterative reweighted to-

365

tal generalized variation based poisson noise removal model, Applied Mathematics and Computation 223 (2013) 264–277.

18

ACCEPTED MANUSCRIPT

[6] T. Le, R. Chartrand, T. J. Asaki, A variational approach to reconstructing images corrupted by poisson noise, J. Math. Imaging Vis. 27 (3) (2007) 257–263. [7] M. A. T. Figueiredo, J. M. Bioucas-Dias, Restoration of poissonian im-

CR IP T

370

ages using alternating direction optimization, IEEE Trans. Image Process. 19 (12) (2010) 3133–3145.

[8] C. Brune, A. Sawatzky, F. W¨ ubbeling, T. K¨ osters, M. Burger, Forward-

backward em-tv methods for inverse problems with poisson noise, Preprint.

[9] S. Setzer, G. Steidl, T. Teuber, Deblurring poissonian images by split breg-

AN US

375

man techniques, Journal of Visual Communication and Image Representation 21 (3) (2010) 193–199.

[10] L. Jiang, J. Huang, X.-G. Lv, J. Liu, Alternating direction method for the high-order total variation-based poisson noise removal problem, Numer. Algorithms 69 (3) (2015) 495–516.

M

380

[11] Z. Zhang, J. Zhang, Z. Wei, L. Xiao, Cartoon-texture composite regulariza-

ED

tion based non-blind deblurring method for partly-textured blurred images with poisson noise, Signal Process. 116 (2015) 127–140.

PT

[12] N. Persch, A. Elhayek, M. Welk, A. Bruhn, S. Grewenig, K. B¨ ose, A. Kraegeloh, J. Weickert, Enhancing 3-D cell structures in confocal and

385

sted microscopy: a joint model for interpolation, deblurring and anisotropic

CE

smoothing, Meas. Sci. Technol. 24 (12) (2013) 125703.

AC

[13] M. M¨ akitalo, A. Foi, Optimal inversion of the anscombe transformation in

390

low-count poisson image denoising, Image Processing, IEEE Transactions on 20 (1) (2011) 99–109.

[14] M. M¨ akitalo, A. Foi, A closed-form approximation of the exact unbiased inverse of the anscombe variance-stabilizing transformation, Image Processing, IEEE Transactions on 20 (9) (2011) 2697–2698.

19

ACCEPTED MANUSCRIPT

[15] D. Chen, Regularized generalized inverse accelerating linearized alternating minimization algorithm for frame-based poissonian image deblurring, SIAM

395

J. Imaging Sci. 7 (2) (2014) 716–739.

CR IP T

[16] M. Carlavan, L. Blanc-F´eraud, Sparse poisson noisy image deblurring, IEEE Trans. Image Process. 21 (4) (2012) 1834–1846.

[17] Y. Zhang, B. Dong, Z. Lu, l0 minimization for wavelet frame based image restoration, Math. Comput. 82 (2013) 995–1015.

400

[18] Y. Jiao, B. Jin, X. Lu, A primal dual active set with continuation algorithm

39 (3) (2015) 400–426.

AN US

for the l0 -regularized optimization problem, Appl. Comput. Harmon. Anal.

[19] L. Shen, Y. Xu, X. Zeng, Wavelet inpainting with the l0 sparse regularization,

405

Applied

and

Computational

Harmonic

Analy-

sis(2015).http://dx.doi.org/10.1016/j.acha.2015.03.001.

M

[20] D.-Q. Dai, L. Shen, Y. Xu, N. Zhang, Noisy 1-bit compressive sensing:

ED

models and algorithms, Appl. Comput. Harmon. Anal. 40 (1) (2016) 1–32. [21] Z. Lu, Y. Zhang, Sparse approximation via penalty decomposition methods, SIAM J. Optim. 23 (4) (2013) 2448–2478.

410

PT

[22] B. Dong, Y. Zhang, An efficient algorithm for l0 minimization in wavelet frame based image restoration, J. Sci. Comput. 54 (2-3) (2013) 350–368.

CE

[23] M. Nikolova, M. K. Ng, C.-P. Tam, Fast nonconvex nonsmooth minimization methods for image restoration and reconstruction, IEEE Trans. Image Process. 19 (12) (2010) 3073–3088.

AC

415

[24] J. Xiao, M. K.-P. Ng, Y.-F. Yang, On the convergence of nonconvex minimization methods for image recovery, IEEE Trans. Image Process. 24 (5) (2015) 1587–1598.

20

ACCEPTED MANUSCRIPT

[25] M. Nikolova, Description of the minimizers of least squares regularized with l0 -norm. uniqueness of the global minimizer, SIAM J. Imaging Sci. 6 (2)

420

(2013) 904–937.

CR IP T

[26] R. Giryes, M. Elad, Sparsity-based poisson denoising with dictionary learning, Image Processing, IEEE Transactions on 23 (12) (2014) 5057–5069.

[27] G. Liuzzi, F. Rinaldi, Solving l0 -penalized problems with simple constraints via the frank-wolfe reduced dimension method, Optim. Lett. 9 (1) (2015)

425

57–74.

AN US

[28] T. Blumensath, M. E. Davies, Iterative thresholding for sparse approximations, J. Fourier Anal. App. 14 (5-6) (2008) 629–654.

[29] J. Nocedal, S. J. Wright, Numerical Optimization, 2nd Edition, Springer, 2006.

430

[30] R. Nishihara, L. Lessard, B. Recht, A. Packard, M. I. Jordan, A general

M

analysis of the convergence of admm, arXiv preprint.

ED

[31] C. Chen, Y. Shen, Y. You, On the convergence analysis of the alternating direction method of multipliers with three blocks, Abstract and Applied Analysis 2013 (r-2) (2013) 1–13.

435

PT

[32] D. Han, X. Yuan, Local linear convergence of the alternating direction method of multipliers for quadratic programs, SIAM Journal on numerical

CE

analysis 51 (6) (2013) 3446–3457. [33] M. Hong, Z.-Q. Luo, M. Razaviyayn, Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems, in:

AC

440

2015 ICASSP, IEEE, 2015, pp. 3836–3840.

[34] Y. Xu, W. Yin, A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion, SIAM Journal on imaging sciences 6 (3) (2013) 1758–1789.

21

ACCEPTED MANUSCRIPT

445

[35] C. Bao, B. Dong, L. Hou, Z. Shen, X. Zhang, X. Zhang, Extrapolated proximal iterative hard thresholding methods for wavelet frame based image restoration, submitted.

CR IP T

[36] H. Zhang, Y. Dong, Q. Fan, The convergence of alternating direction method of multipliers for an l0 regularized image restoration model, Technical Report, Wuhan Uiversity 1 (May) (2016) 1–12.

450

[37] Q. Fan, D. Jiang, Y. Jiao, A multi-parameter regularization model for image restoration, Signal Process. 114 (2015) 131–142.

AN US

[38] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed opti-

mization and statistical learning via the alternating direction method of multipliers, Found. Trends Mach. Learn. 3 (1) (2011) 1–122.

455

[39] M. K. Ng, R. H. Chan, W.-C. Tang, A fast algorithm for deblurring models with neumann boundary conditions, SIAM J. Sci. Comput. 21 (3) (1999)

M

851–866.

[40] B. He, X. Yuan, Convergence analysis of primal-dual algorithms for a saddle-point problem: from contraction perspective, SIAM J. Imaging Sci.

ED

460

5 (1) (2012) 119–149.

PT

[41] B. He, X. Yuan, On the o(1/n) convergence rate of the douglas-rachford alternating direction method, SIAM J. Numer. Anal. 50 (2) (2012) 700–709.

CE

[42] Y. Wang, W. Yin, J. Zeng, Global convergence of admm in nonconvex nonsmooth optimization, arXiv preprint arXiv:1511.06324.

465

AC

[43] Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, Image Processing, IEEE Transactions on 13 (4) (2004) 600–612.

22

CR IP T

ACCEPTED MANUSCRIPT

Table 1: The recovered results obtain from the FT-ADMM and CG-ADMM algorithms with different initialization conditions for u0 . The test image is Cameraman(size:256 × 256).

Initial u0

SNR(dB)

RelErr

SSIM

Iter

Time (s)

FT-ADMM

0

20.47

0.0947

0.7584

17

3.65

f

20.38

0.0958

0.7598

13

2.77

0

21.57

0.0836

0.8226

20

11.22

f

19.28

0.1086

0.5732

14

7.83

ED

M

CG-ADMM

AN US

Algorithm

PT

Table 2: The recovered results obtain from the FT-ADMM and CG-ADMM algorithms with different wavelet frame decomposition level L. The test image is Pepper(size:512 × 512).

CE

Algorithm

AC

FT-ADMM

CG-ADMM

L=1

L=2

L=3

L=4

SNR/SSIM/Time

SNR/SSIM/Time

SNR/SSIM/Time

SNR/SSIM/Time

22.68/0.7184/6.87

23.73/0.7846/12.41

23.88/0.7871/18.02

23.93/0.7860/23.69

23.50/0.7566/18.61

24.78/0.8237/25.93

24.85/0.8262/36.35

24.85/0.8253/48.86

23

ACCEPTED MANUSCRIPT

Table 3: Restoration results with our newly proposed FT-ADMM and CG-ADMM algorithms, GILAM[15], PIDAL[7] and ADMPNR [10] algorithms, respectively. All these quantities are the ten times average with each test image.

Cameraman

Pepper

RelErr

SSIM

l0 -FT

20.85

0.0906

0.8204

l0 -CG

21.17

0.0874

0.8539

l1

20.38

0.0957

0.8379

TV

18.99

0.1123

0.6882

4th-TV

19.81

0.1022

0.8170

l0 -FT

20.15

0.0983

0.7561

l0 -CG

21.59

0.0833

0.8228

l1

20.62

0.0931

TV

18.67

4th-TV

19.11

l0 -FT

23.34

l0 -CG

24.75

l1

23.28

1.91

20

6.28

8

3.52

17

1.09

19

1.27

13

1.46

20

6.24

0.7960

8

3.54

0.1165

0.6253

15

0.92

0.1108

0.7769

14

0.95

0.0681

0.7763

14

8.57

0.0579

0.8236

14

25.38

0.0685

0.8067

8

21.05

21.34

0.0857

0.6442

13

2.34

23.30

0.0684

0.8139

17

5.54

l0 -FT

18.20

0.1230

0.6259

12

7.31

l0 -CG

18.51

0.1187

0.6469

19

32.77

l1

17.46

0.1339

0.5196

7

19.54

TV

17.93

0.1265

0.6162

12

2.22

4th-TV

16.77

0.1451

0.4496

15

4.95

l0 -FT

11.67

0.2645

0.8021

19

2.08

l0 -CG

14.78

0.1823

0.8556

20

6.97

l1

13.33

0.1394

0.8328

18

7.40

TV

10.11

0.3129

0,7221

3

0.10

4th-TV

11.80

0.2569

0.8821

30

2.01

l0 -FT

24.47

0.0610

0.8512

12

1.37

l0 -CG

27.00

0.0447

0.8658

15

4.63

l1

25.18

0.0551

0.8328

8

3.36

TV

24.05

0.0627

0.7851

10

0.29

4th-TV

25.50

0.0531

0.8662

16

1.05

CE

PT

ED

4th-TV

AC

Blobs

Moon

Time (s)

17

TV

Baboon

Iter

CR IP T

SNR(dB)

AN US

Lena

Model

M

Image

24

(b)

(c)

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(d)

(e)

(f)

M

Figure 1: (a) Original image ”Cameraman”, size: 256 × 256. (b) Original image ”Lena”, size:256 × 256. (c) Original image ”Pepper”, size: 512 × 512. (d) Original image ”Baboon”, size: 512 × 512. (e) Original image ”Blobs”, size: 256 × 256. (f) Original image ”Moon”, size:

ED

256 × 256.

−1.5

log10(RelErr(uk,uk+1)

−1

FT−ADMM CG−ADMM −0.5

CE

log10(RelErr(uk,uk+1)

−0.5

0

FT−ADMM CG−ADMM

PT

0

−2

AC

−2.5

−3

0

5

−1

−1.5

−2

−2.5

10 Iteration

15

−3

20

(a)

0

5

10 Iteration

15

20

(b)

Figure 2: The relative difference curve of FT-ADMM and CG-ADMM for two test images: (a) Cameraman; (b) Lena.

25

(b)

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(e)

(f)

CE

PT

ED

(d)

(c)

(g)

(h)

(i)

Figure 3: Comparison of the denoised and deblurred images ”Lena”. (a) Degraded image with

AC

Poisson noise and Gaussian blur, (b) restored image with the l0 -FT, (c) restored image with

the l0 -CG, (d) restored image with the l1 model, (e) restored image with the TV model, (f) restored image with the high-order TV model.

26

(b)

(c)

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(e)

(f)

PT

ED

M

(d)

(h)

(i)

CE

(g)

Figure 4: Comparison of the denoised and deblurred images ”Baboon”. (a) Degraded image with Poisson noise and Gaussian blur, (b) restored image with the l0 -FT, (c) restored image

AC

with the l0 -CG, (d) restored image with the l1 model, (e) restored image with the TV model, (f) restored image with the high-order TV model. (g) zoomed part of the restored image with l0 -CG, (h)zoomed part of the restored image with l1 model, (i) zoomed part of the restored image with high-order TV model.

27

(b)

(c)

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(e)

(f)

PT

ED

(d)

(h)

(i)

CE

(g)

Figure 5: Comparison of the denoised and deblurred images ”Blobs”. (a) Degraded image with Poisson noise and Gaussian blur, (b) restored image with the l0 -FT, (c) restored image

AC

with the l0 -CG, (d) restored image with the l1 model, (e) restored image with the TV model, (f) restored image with the high-order TV model. (g) zoomed part of the restored image with l0 -CG, (h)zoomed part of the restored image with l1 model, (i) zoomed part of the restored image with high-order TV model.

28

(b)

(c)

M

AN US

(a)

CR IP T

ACCEPTED MANUSCRIPT

(e)

(f)

PT

ED

(d)

(h)

(i)

CE

(g)

Figure 6: Comparison of the denoised and deblurred images ”Moon”. (a) Degraded image with Poisson noise and Gaussian blur, (b) restored image with the l0 -FT, (c) restored image

AC

with the l0 -CG, (d) restored image with the l1 model, (e) restored image with the TV model, (f) restored image with the high-order TV model. (g) zoomed part of the restored image with l0 -CG, (h)zoomed part of the restored image with l1 model, (i) zoomed part of the restored image with high-order TV model.

29