A L0 sparse analysis prior for blind poissonian image deconvolution Xiaojin Gong,1,∗ Baisheng Lai,1 and Zhiyu Xiang1 Dept. of Information Science and Electronic Engineering, Zhejiang University, Hangzhou, 310027, China ∗ [email protected]
Abstract: This paper proposes a new approach for blindly deconvolving images that are contaminated by Poisson noise. The proposed approach incorporates a new prior, that is the L0 sparse analysis prior, together with the total variation constraint into the maximum a posteriori (MAP) framework for deconvolution. A greedy analysis pursuit numerical scheme is exploited to solve the L0 regularized MAP problem. Experimental results show that our approach not only produces smooth results substantially suppressing artifacts and noise, but also preserves intensity changes sharply. Both quantitative and qualitative comparisons to the specialized state-of-the-art algorithms demonstrate its superiority. © 2014 Optical Society of America OCIS codes: (100.1455) Blind deconvolution; (100.1830) Deconvolution; (100.3020) Digital image processing; (100.2000) Inverse problem.
References and links 1. D. Kundur and D. Hatzinakos, “Blind image deconvolution revisited,” IEEE Signal Proc. Mag. 11, 61–63 (1996). 2. P. Campisi and K. Egiazarian, Blind Image Deconvolution: Theory and Applications (CRC Press, 2007). 3. M. Demenikov and A. R. Harvey, “Parametric blind-deconvolution algorithm to remove image artifacts in hybrid imaging systems,” Opt. Express 18(17), 18035–18040 (2010). 4. S. V. Vorontsov, V. N. Strakhov, S. M. Jefferies, and K. J. Borelli, “Deconvolution of astronomical images using SOR with adaptive relaxation,” Opt. Express 19(14), 13509–13524 (2011) 5. F. Dupe, J. M. Fadili, and J. Starck, “A proximal iteration for deconvolving poisson noisy images using sparse representations,” IEEE Trans. Image Process. 18(2), 310–321 (2009). 6. W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. 62, 55–59 (1972). 7. L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. 79, 745–754 (1974). 8. D. A. Fish, A. M. Brinicombe, E. R. Pike and J. G. Walker, “Blind deconvolution by means of the richardson-lucy algorithm,” J. Opt. Soc. Am. 12, 58–65 (1995). 9. Z. Xu and E. Lam, “Maximum a posteriori blind image deconvolution with huber-Markov random-field regularization,” Opt. Lett. 34, 1453–1455 (2009). 10. L. Yan, H. Fang, and S. Zhong, “Blind image deconvolution with spatially adaptive total variation regularization,” Opt. Lett. 37, 2778–2780 (2012). 11. S. Nam, M. E. Davies, M. Elad, and R. Gribonval, “The cosparse analysis model and algorithms,” Appl. Comput. Harmon. A. 34, 30–56 (2013). 12. J.-L. Starck, E. J. Candes, and D. L. Donoho, “The curvelet transform for image denoising,” IEEE Trans. Image Process. 11, 670–684 (2002). 13. J. Cai, H. Ji, C. Liu, and Z. Shen, “Framelet-based blind motion deblurring from a single image,” IEEE Trans. Image Process. 21, 562–572 (2012). 14. L. Xu, S. Zheng and J. Jia, “Unnatural l0 sparse representation for natural image deblurring,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 1107–1114 (2013). 15. http://www.imageprocessingplace.com/root files V3/image databases.htm 16. H. Fang, L. Yan, H. Liu, and Y. Chang, “Blind poissonian images deconvolution with framelet regularization,” Opt. Lett. 38, 389–391 (2013). 17. T. Chan and C. Wong, “Total variation blind deconvolution,” IEEE Trans. Image Process. 7, 370–375 (1998).
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3860
1.
Introduction
Blind Poissonian image deconvolution is a challenging task, to which much research has been devoted as in microscopy, medical and astronomical imaging [1–5]. Among the extensively developed techniques, the most widely used one is the Richardson-Lucy (RL) algorithm [6–8]. It relies on a Bayesian model and formulates deconvolution as a maximum likelihood estimation (MLE) problem, so that it is able to take into account Poisson statistics of the noise. However, due to the ill-posed nature of deconvolution, RL cannot prevent from noise amplification, leading to ’ringing’ artifacts. Therefore, extra regularizations such as the Huber-Markov random field [9], total variation (TV) and spatially adaptive TV [10] have been incorporated. These regularized MLE methods, equivalent to the maximum a posteriori (MAP) estimation with the incorporation of priors, successfully achieve better performance. But unfortunately, these methods are commonly prone to over-smooth intensity edges, especially when images are suffer from significant noise and heavy blur. In this paper, we propose to incorporate a new regularizer, that is the L0 sparse analysis prior [11], together with the total variation regularization into the MAP framework for deconvolution. This sparse analysis prior implies that the analyzed result Ωx of a signal x ∈ Rn is sparse. Here, Ω ∈ Rm×n is a redundant analysis operator with m ≥ n. Common examples of Ω include the finite difference operator, the wavelet transform, curvelet [12], framelet [13] and so on. The sparsity of the analysis result in such a transform domain is commonly represented by a L1 regularization and has been successfully applied for various signal processing tasks in recent years. More recently, the L0 scheme has started to be investigated [11, 14]. It, as a high-sparsity-pursuiting regularizer, assures to preserve salient intensity changes better. Therefore, we explore its performance in blind image deconvolution. Moreover, a simple but effective greedy-like approach [11] is further exploited to solve the L0 regularized MAP problem. Benefitted from the joint use of the L0 sparse analysis prior and the TV regularization, our approach can not only produce smooth results that substantially suppress artifacts and noise, but also preserve intensity changes sharply. 2.
Proposed algorithm
Let us consider an imaging system that is linear and spatial invariant. Then, an observed image y ∈ Rn can be viewed as a true image x ∈ Rn convolved with a point spread function (PSF) h and corrupted by a Poisson noise process P. In matrix-vector notation, this image formation model can be represented by the following two equivalent forms: y = P(Hx) = P(Xh),
(1)
where H ∈ Rn×n denotes the matrix notation of the convolution of the PSF; X is the convolution matrix of x; both x and y stand for the vectorized images. Assuming that the intensities of all pixels in the observed image are random variables independent and identically distributed, and each yi obeys a Poisson distribution, we get the likelihood in the form of n (Hx)yi exp(−Hx) i i . (2) p(y|x, h) = ∏ y ! i i=1 Further, we model the joint regularization as a prior defined below: p(x) ∝ exp(−α||x||TV ) exp(−β ||Ωx||0 ),
(3)
where α and β are two parameters regularizing the weights of the TV and the L0 sparse priors respectively.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3861
Using the Bayesian model, we can reconstruct the image x via maximizing the a posteriori probability p(x, h|y). This MAP estimation is equivalent to the following minimization problem: arg min J(x, H, y) + α||x||TV + β ||Ωx||0 , {x,h}
(4)
subject to x ≥ 0 where J(x, H, y) = 1T Hx − yT log(Hx), and 1 is a n-size vector whose all entries are 1. The constraint x ≥ 0 is placed to ensure the non-negativity of x. The problem defined in Eq.(4) is NP-complete due to the L0 sparse analysis term. In this work, we exploit a simple but effective approach, the Greedy Analysis Pursuit (GAP) algoˆ 0 , which is defined to be rithm [11], to solve it iteratively. GAP starts from a cosupport set Λ the indexes of the zero entries of Ωx and is initialized to be the whole set. An initial solution is estimated via {ˆx0 , hˆ 0 } = arg min J(x, H, y) + α||x||TV + β ||ΩΛˆ 0 x||2 {x,h}
(5)
subject to x ≥ 0. Then, elements outside the cosupport set are detected according to the values of Ωˆx0 and removed. With the updated cosupport set, x and h are re-estimated for removing more elements. By this means, GAP iteratively reduces the cosupport set till a stop criterion is reached. The details of GAP-based blind image deconvolution is illustrated in Algorithm 1. Algorithm 1: GAP for Blind Image Deconvolution Input: Ω, y, N, t ∈ (0, 1]. Output: xˆ = xˆ k ; hˆ = hˆ k Initialization: k = 0; ˆ 0 = {1, 2, ..., m}; Initialize cosupport: Λ Initialize solution xˆ 0 and hˆ 0 by solving arg min J(x, H, y) + α||x||TV + β ||ΩΛˆ 0 x||2 , {x,h}
subject to x ≥ 0. while k < N do k = k + 1; Compute z = Ωˆxk−1 ; Find largest entries: Γk = {i : |zi | ≥ t max j |z j |} ˆk =Λ ˆ k−1 \Γk ; Update cosupport: Λ ˆ Update solution xˆ k and hk by solving arg min J(x, H, y) + α||x||TV + β ||ΩΛˆ k x||2 , {x,h}
subject to x ≥ 0. end The problem defined in Eq.(5) and those within the GAP iteration are solved by the
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3862
expectation-maximization approach. That is, h is obtained iteratively via y hk+1 = hk X∗ Xhk hk+1 hk+1 = , ||hk+1 ||2
(6)
and the image x is estimated by xk+1 =
xk 1 − αdiv
∇x ||∇x||2
y ∗ H , Hxk + β ΩTΛ ΩΛ x
(7)
xk+1 = max {xk+1 , 0} where X∗ , H∗ denote the adjoint operators of X and H. 3.
Experiments
In order to evaluate the performance of our approach, we first conduct experiments on four test images: a 255 × 252 Cameraman image [15], a 250 × 250 satellite image [16], a 300 × 300 brain image [16], and the 230 × 230 spiral nebula NGC 7027 image [15]. Each degraded image is simulated by convolving the original one with a 15 × 15 Gaussian blur kernel of a 2.5 standard deviation, and further contaminated with Poisson noise. We compare our approach to a classical method and a state-of-the-art approach, which are, respectively, the Richardson-Lucy total variation (RLTV) [17] and the blind poissonian images deconvolution with framelet regularization (BPIDFR) methods [16]. Experimental results are quantitatively assessed in terms of the normalized mean square error (NMSE), which is defined by NMSE =
||ˆx − x||22 , ||x||22
(8)
where x and xˆ are, respectively, the original and the restored image. In all experiments, in order to achieve small NMSEs, we manually set the parameters as follows: α = 0.0015, β = 0.00001, t = 0.6 and 30 GAP iterations. Moreover, for the analysis operator Ω, a candidate could be the wavelet, curvelet, framelet, or other pre-defined transforms. Here, we simply employ the level 4 Daubechies wavelets (db2-wavelets), which is of less computational complexity than the framelet transform used in BPIDFR. The NMSE of each experimental result is listed in Table 1. From it we can observe that our approach performs much better than RLTV in all cases. When comparing to BPIDFR, our approach achieves quantitatively the same results on the Cameraman and the brain images, but outperforms BPIDFR in a considerable extent for the other two cases. Image Restoration by Image Cameraman Satellite Brain Nebula
Degraded Image 0.0242 0.0741 0.0567 0.0951
RLTV 0.0107 0.0500 0.0293 0.0615
BPIDFR 0.0064 0.0393 0.0166 0.0510
TVL0 0.0063 0.0372 0.0164 0.0430
Table 1. NMSE of the degraded images and the restored images obtained by different algorithms. Here, TVL0 refers to our proposed approach.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3863
In order to achieve visual comparisons on the performance, Fig. 1 presents all the degraded images and the restored results. From the Cameraman scenario, we can obviously observe that RLTV yields artifacts near object boundaries, while both BPIDFR and our approach suppress such artifacts quite well. Moreover, our approach suppresses noise better than BPIDFR. Such observations are also validated in the remaining cases. More generally speaking, the proposed approach produces smoother results while preserving intensity changes, including both object boundaries the and isolated points such as stars in the nebula image, sharper than RLTV and BPIDFR.
(a) Degraded Image
(b) RLTV
(c) BPIDFR
(d) TVL0
Fig. 1. Restoration of the degraded images. (a) is the degraded image. (b), (c), (d) are the results restored by the RLTV, BPIDFR, and TVL0 approaches respectively. From top to bottom, the scenarios are Cameraman, satellite, brain, and nebula.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3864
We further validate the capability of our approach on a real Saturn image [16] and a real neuron image [5]. In this experiment, the same parameters are used as previous, but the PSF is initialized as a 25 × 25 Gaussian kernel of a 1.5 standard deviation. Fig. 2 illustrates the restored results. From it we observe that RLTV yields piecewise smooth artifacts in the Saturn case and fails to restore part of the thin structures for neuron. In contrast, both BPIDFR and our approach suppress the artifacts and restore thin structures better. The restored images obtained by our approach preserves intensity changes slightly sharper.
(a) Degraded
(b) RLTV
(c) BPIDFR
(d) TVL0
Fig. 2. Restoration of the degraded images. (a) is the degraded image. (b), (c), (d) are the images restored by the RLTV, BPIDFR, and TVL0 approaches respectively. The top is the Saturn and the middle is neuron. The bottom shows a zoomed part of the neuron image.
4.
Conclusion
In conclusion, this paper has presented a new blind Poissonian image deconvolution method, which incorporates the L0 sparse analysis prior, together with the total variation constraint within the MAP framework. A greedy analysis pursuit algorithm has been exploited to solve the formulated L0 regularized MAP problem. The proposed approach substantially suppresses artifacts and noise, producing smooth results while preserving intensity changes sharp. Comparative results on both simulated and real images show that our approach outperforms the classical RLTV and the state-of-the-art framelet regularized (BPIDFR) methods. Acknowledgments This work was supported in part by the National Natural Science Foundation of China via grants 61001171, 61071219, and Chinese Universities Scientific Fund. The authors are grateful to Prof. Luxin Yan and Mr. Houzhang Fang (the authors of BPIDFR [16]) for providing the satellite and brain test images. We are also grateful to the reviewers for the valuable suggestions.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3865
Abstract: This paper proposes a new approach for blindly deconvolving images that are contaminated by Poisson noise. The proposed approach incorporates a new prior, that is the L0 sparse analysis prior, together with the total variation constraint into the maximum a posteriori (MAP) framework for deconvolution. A greedy analysis pursuit numerical scheme is exploited to solve the L0 regularized MAP problem. Experimental results show that our approach not only produces smooth results substantially suppressing artifacts and noise, but also preserves intensity changes sharply. Both quantitative and qualitative comparisons to the specialized state-of-the-art algorithms demonstrate its superiority. © 2014 Optical Society of America OCIS codes: (100.1455) Blind deconvolution; (100.1830) Deconvolution; (100.3020) Digital image processing; (100.2000) Inverse problem.
References and links 1. D. Kundur and D. Hatzinakos, “Blind image deconvolution revisited,” IEEE Signal Proc. Mag. 11, 61–63 (1996). 2. P. Campisi and K. Egiazarian, Blind Image Deconvolution: Theory and Applications (CRC Press, 2007). 3. M. Demenikov and A. R. Harvey, “Parametric blind-deconvolution algorithm to remove image artifacts in hybrid imaging systems,” Opt. Express 18(17), 18035–18040 (2010). 4. S. V. Vorontsov, V. N. Strakhov, S. M. Jefferies, and K. J. Borelli, “Deconvolution of astronomical images using SOR with adaptive relaxation,” Opt. Express 19(14), 13509–13524 (2011) 5. F. Dupe, J. M. Fadili, and J. Starck, “A proximal iteration for deconvolving poisson noisy images using sparse representations,” IEEE Trans. Image Process. 18(2), 310–321 (2009). 6. W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Am. 62, 55–59 (1972). 7. L. B. Lucy, “An iterative technique for the rectification of observed distributions,” Astron. J. 79, 745–754 (1974). 8. D. A. Fish, A. M. Brinicombe, E. R. Pike and J. G. Walker, “Blind deconvolution by means of the richardson-lucy algorithm,” J. Opt. Soc. Am. 12, 58–65 (1995). 9. Z. Xu and E. Lam, “Maximum a posteriori blind image deconvolution with huber-Markov random-field regularization,” Opt. Lett. 34, 1453–1455 (2009). 10. L. Yan, H. Fang, and S. Zhong, “Blind image deconvolution with spatially adaptive total variation regularization,” Opt. Lett. 37, 2778–2780 (2012). 11. S. Nam, M. E. Davies, M. Elad, and R. Gribonval, “The cosparse analysis model and algorithms,” Appl. Comput. Harmon. A. 34, 30–56 (2013). 12. J.-L. Starck, E. J. Candes, and D. L. Donoho, “The curvelet transform for image denoising,” IEEE Trans. Image Process. 11, 670–684 (2002). 13. J. Cai, H. Ji, C. Liu, and Z. Shen, “Framelet-based blind motion deblurring from a single image,” IEEE Trans. Image Process. 21, 562–572 (2012). 14. L. Xu, S. Zheng and J. Jia, “Unnatural l0 sparse representation for natural image deblurring,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, 1107–1114 (2013). 15. http://www.imageprocessingplace.com/root files V3/image databases.htm 16. H. Fang, L. Yan, H. Liu, and Y. Chang, “Blind poissonian images deconvolution with framelet regularization,” Opt. Lett. 38, 389–391 (2013). 17. T. Chan and C. Wong, “Total variation blind deconvolution,” IEEE Trans. Image Process. 7, 370–375 (1998).
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3860
1.
Introduction
Blind Poissonian image deconvolution is a challenging task, to which much research has been devoted as in microscopy, medical and astronomical imaging [1–5]. Among the extensively developed techniques, the most widely used one is the Richardson-Lucy (RL) algorithm [6–8]. It relies on a Bayesian model and formulates deconvolution as a maximum likelihood estimation (MLE) problem, so that it is able to take into account Poisson statistics of the noise. However, due to the ill-posed nature of deconvolution, RL cannot prevent from noise amplification, leading to ’ringing’ artifacts. Therefore, extra regularizations such as the Huber-Markov random field [9], total variation (TV) and spatially adaptive TV [10] have been incorporated. These regularized MLE methods, equivalent to the maximum a posteriori (MAP) estimation with the incorporation of priors, successfully achieve better performance. But unfortunately, these methods are commonly prone to over-smooth intensity edges, especially when images are suffer from significant noise and heavy blur. In this paper, we propose to incorporate a new regularizer, that is the L0 sparse analysis prior [11], together with the total variation regularization into the MAP framework for deconvolution. This sparse analysis prior implies that the analyzed result Ωx of a signal x ∈ Rn is sparse. Here, Ω ∈ Rm×n is a redundant analysis operator with m ≥ n. Common examples of Ω include the finite difference operator, the wavelet transform, curvelet [12], framelet [13] and so on. The sparsity of the analysis result in such a transform domain is commonly represented by a L1 regularization and has been successfully applied for various signal processing tasks in recent years. More recently, the L0 scheme has started to be investigated [11, 14]. It, as a high-sparsity-pursuiting regularizer, assures to preserve salient intensity changes better. Therefore, we explore its performance in blind image deconvolution. Moreover, a simple but effective greedy-like approach [11] is further exploited to solve the L0 regularized MAP problem. Benefitted from the joint use of the L0 sparse analysis prior and the TV regularization, our approach can not only produce smooth results that substantially suppress artifacts and noise, but also preserve intensity changes sharply. 2.
Proposed algorithm
Let us consider an imaging system that is linear and spatial invariant. Then, an observed image y ∈ Rn can be viewed as a true image x ∈ Rn convolved with a point spread function (PSF) h and corrupted by a Poisson noise process P. In matrix-vector notation, this image formation model can be represented by the following two equivalent forms: y = P(Hx) = P(Xh),
(1)
where H ∈ Rn×n denotes the matrix notation of the convolution of the PSF; X is the convolution matrix of x; both x and y stand for the vectorized images. Assuming that the intensities of all pixels in the observed image are random variables independent and identically distributed, and each yi obeys a Poisson distribution, we get the likelihood in the form of n (Hx)yi exp(−Hx) i i . (2) p(y|x, h) = ∏ y ! i i=1 Further, we model the joint regularization as a prior defined below: p(x) ∝ exp(−α||x||TV ) exp(−β ||Ωx||0 ),
(3)
where α and β are two parameters regularizing the weights of the TV and the L0 sparse priors respectively.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3861
Using the Bayesian model, we can reconstruct the image x via maximizing the a posteriori probability p(x, h|y). This MAP estimation is equivalent to the following minimization problem: arg min J(x, H, y) + α||x||TV + β ||Ωx||0 , {x,h}
(4)
subject to x ≥ 0 where J(x, H, y) = 1T Hx − yT log(Hx), and 1 is a n-size vector whose all entries are 1. The constraint x ≥ 0 is placed to ensure the non-negativity of x. The problem defined in Eq.(4) is NP-complete due to the L0 sparse analysis term. In this work, we exploit a simple but effective approach, the Greedy Analysis Pursuit (GAP) algoˆ 0 , which is defined to be rithm [11], to solve it iteratively. GAP starts from a cosupport set Λ the indexes of the zero entries of Ωx and is initialized to be the whole set. An initial solution is estimated via {ˆx0 , hˆ 0 } = arg min J(x, H, y) + α||x||TV + β ||ΩΛˆ 0 x||2 {x,h}
(5)
subject to x ≥ 0. Then, elements outside the cosupport set are detected according to the values of Ωˆx0 and removed. With the updated cosupport set, x and h are re-estimated for removing more elements. By this means, GAP iteratively reduces the cosupport set till a stop criterion is reached. The details of GAP-based blind image deconvolution is illustrated in Algorithm 1. Algorithm 1: GAP for Blind Image Deconvolution Input: Ω, y, N, t ∈ (0, 1]. Output: xˆ = xˆ k ; hˆ = hˆ k Initialization: k = 0; ˆ 0 = {1, 2, ..., m}; Initialize cosupport: Λ Initialize solution xˆ 0 and hˆ 0 by solving arg min J(x, H, y) + α||x||TV + β ||ΩΛˆ 0 x||2 , {x,h}
subject to x ≥ 0. while k < N do k = k + 1; Compute z = Ωˆxk−1 ; Find largest entries: Γk = {i : |zi | ≥ t max j |z j |} ˆk =Λ ˆ k−1 \Γk ; Update cosupport: Λ ˆ Update solution xˆ k and hk by solving arg min J(x, H, y) + α||x||TV + β ||ΩΛˆ k x||2 , {x,h}
subject to x ≥ 0. end The problem defined in Eq.(5) and those within the GAP iteration are solved by the
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3862
expectation-maximization approach. That is, h is obtained iteratively via y hk+1 = hk X∗ Xhk hk+1 hk+1 = , ||hk+1 ||2
(6)
and the image x is estimated by xk+1 =
xk 1 − αdiv
∇x ||∇x||2
y ∗ H , Hxk + β ΩTΛ ΩΛ x
(7)
xk+1 = max {xk+1 , 0} where X∗ , H∗ denote the adjoint operators of X and H. 3.
Experiments
In order to evaluate the performance of our approach, we first conduct experiments on four test images: a 255 × 252 Cameraman image [15], a 250 × 250 satellite image [16], a 300 × 300 brain image [16], and the 230 × 230 spiral nebula NGC 7027 image [15]. Each degraded image is simulated by convolving the original one with a 15 × 15 Gaussian blur kernel of a 2.5 standard deviation, and further contaminated with Poisson noise. We compare our approach to a classical method and a state-of-the-art approach, which are, respectively, the Richardson-Lucy total variation (RLTV) [17] and the blind poissonian images deconvolution with framelet regularization (BPIDFR) methods [16]. Experimental results are quantitatively assessed in terms of the normalized mean square error (NMSE), which is defined by NMSE =
||ˆx − x||22 , ||x||22
(8)
where x and xˆ are, respectively, the original and the restored image. In all experiments, in order to achieve small NMSEs, we manually set the parameters as follows: α = 0.0015, β = 0.00001, t = 0.6 and 30 GAP iterations. Moreover, for the analysis operator Ω, a candidate could be the wavelet, curvelet, framelet, or other pre-defined transforms. Here, we simply employ the level 4 Daubechies wavelets (db2-wavelets), which is of less computational complexity than the framelet transform used in BPIDFR. The NMSE of each experimental result is listed in Table 1. From it we can observe that our approach performs much better than RLTV in all cases. When comparing to BPIDFR, our approach achieves quantitatively the same results on the Cameraman and the brain images, but outperforms BPIDFR in a considerable extent for the other two cases. Image Restoration by Image Cameraman Satellite Brain Nebula
Degraded Image 0.0242 0.0741 0.0567 0.0951
RLTV 0.0107 0.0500 0.0293 0.0615
BPIDFR 0.0064 0.0393 0.0166 0.0510
TVL0 0.0063 0.0372 0.0164 0.0430
Table 1. NMSE of the degraded images and the restored images obtained by different algorithms. Here, TVL0 refers to our proposed approach.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3863
In order to achieve visual comparisons on the performance, Fig. 1 presents all the degraded images and the restored results. From the Cameraman scenario, we can obviously observe that RLTV yields artifacts near object boundaries, while both BPIDFR and our approach suppress such artifacts quite well. Moreover, our approach suppresses noise better than BPIDFR. Such observations are also validated in the remaining cases. More generally speaking, the proposed approach produces smoother results while preserving intensity changes, including both object boundaries the and isolated points such as stars in the nebula image, sharper than RLTV and BPIDFR.
(a) Degraded Image
(b) RLTV
(c) BPIDFR
(d) TVL0
Fig. 1. Restoration of the degraded images. (a) is the degraded image. (b), (c), (d) are the results restored by the RLTV, BPIDFR, and TVL0 approaches respectively. From top to bottom, the scenarios are Cameraman, satellite, brain, and nebula.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3864
We further validate the capability of our approach on a real Saturn image [16] and a real neuron image [5]. In this experiment, the same parameters are used as previous, but the PSF is initialized as a 25 × 25 Gaussian kernel of a 1.5 standard deviation. Fig. 2 illustrates the restored results. From it we observe that RLTV yields piecewise smooth artifacts in the Saturn case and fails to restore part of the thin structures for neuron. In contrast, both BPIDFR and our approach suppress the artifacts and restore thin structures better. The restored images obtained by our approach preserves intensity changes slightly sharper.
(a) Degraded
(b) RLTV
(c) BPIDFR
(d) TVL0
Fig. 2. Restoration of the degraded images. (a) is the degraded image. (b), (c), (d) are the images restored by the RLTV, BPIDFR, and TVL0 approaches respectively. The top is the Saturn and the middle is neuron. The bottom shows a zoomed part of the neuron image.
4.
Conclusion
In conclusion, this paper has presented a new blind Poissonian image deconvolution method, which incorporates the L0 sparse analysis prior, together with the total variation constraint within the MAP framework. A greedy analysis pursuit algorithm has been exploited to solve the formulated L0 regularized MAP problem. The proposed approach substantially suppresses artifacts and noise, producing smooth results while preserving intensity changes sharp. Comparative results on both simulated and real images show that our approach outperforms the classical RLTV and the state-of-the-art framelet regularized (BPIDFR) methods. Acknowledgments This work was supported in part by the National Natural Science Foundation of China via grants 61001171, 61071219, and Chinese Universities Scientific Fund. The authors are grateful to Prof. Luxin Yan and Mr. Houzhang Fang (the authors of BPIDFR [16]) for providing the satellite and brain test images. We are also grateful to the reviewers for the valuable suggestions.
#204012 - $15.00 USD Received 7 Jan 2014; revised 31 Jan 2014; accepted 2 Feb 2014; published 11 Feb 2014 (C) 2014 OSA 24 February 2014 | Vol. 22, No. 4 | DOI:10.1364/OE.22.003860 | OPTICS EXPRESS 3865
