Poisson Image Denoising Using Best Linear Prediction: A Post-processing Framework
PPoisson Image Denoising Using Best LinearPrediction: A Post-processing Framework
Milad Niknejad, M´ario A.T. Figueiredo
Instituto de Telecomunicac¸ ˜oes, andInstituto Superior T´ecnico, Universidade de Lisboa, Portugal
Abstract —In this paper, we address the problem of denoisingimages degraded by Poisson noise. We propose a new patch-based approach based on best linear prediction to estimatethe underlying clean image. A simplified prediction formula isderived for Poisson observations, which requires the covariancematrix of the underlying clean patch. We use the assumption thatsimilar patches in a neighborhood share the same covariancematrix, and we use off-the-shelf Poisson denoising methods inorder to obtain an initial estimate of the covariance matrices.Our method can be seen as a post-processing step for Poissondenoising methods and the results show that it improves uponseveral Poisson denoising methods by relevant margins.
I. I
NTRODUCTION
Image denoising in the presence of Poissonian observationsis an important problem, appearing in several application areassuch as astronomy and medical imaging, where photon-limitedimages are very common.In Poisson denoising, each pixel in the noisy image isa realization of a Poisson random variable with expectedvalue equal to the true underlying pixel to be estimated.Some Poisson image denoising methods apply a nonlinear variance stabilization transform (VST), such as the Anscombetransform [1], to the noisy image, in order to approximatelytransform the noise into Gaussian-distributed. The resultingimage is then processed using a denoising method for Gaussiannoise; finally, an inverse of the transform is applied to theoutput of denoising method, in order to obtain the final imageestimate. Some methods have been proposed to improve theaccuracy of this inverse transform [2]. However, VST-basedmethods are less accurate in regimes of low signal-to-noiseratios (SNR), which led to the development of methods thatdirectly tackle the Poissonian image. Two recent examples ofsuch denoising methods are
Poisson non-local means (PNLM)[3] and non-local PCA (NL-PCA) [4] (more references to ear-lier work can be found in [3] and [4]). Approaches that directlydeal with the Poissonian nature of the observations have alsobeen proposed for image deconvolution under Poisson noise(see [5], [6], and references therein).The underlying assumption in many image restorationmethods is that the clean patches lie in a low dimensionalsubspace. Whereas in the Gaussian noise case, this assump-tion can be used efficiently [7], [8], in the Possionian case,applying this assumption as a prior leads to an intractableoptimization problem due to the non-quadratic log-likelihood.To circumvent this difficulty, some methods, such as NL-PCA and the sparsity-based method in [9], assume the patches are represented as the exponential of some data lying in alow-dimensional subspace. However, the exponential of low-dimensional data is not guaranteed to be low-dimensional.Recently, approaches based on sampling have been proposedto directly approximate the minimum mean squared error (MMSE) estimate [10]. Although these approaches are gener-ally computationally complex for generic images, it has beenshown that it can be effectively applied in a class-specificsetting [11].In this paper, we propose a new method based on the bestlinear prediction (BLP) principle, in order to estimate theclean patches as a linear transformation of the noisy ones.This estimator is independent of the noise distribution and itrequires merely the covariance matrix of the underlying cleanvector. We assume that similar patches in a neighborhoodof a reference patch share the same covariance matrix. Inorder to estimate this covariance matrix, we use an off-the-shelf Poisson denoising algorithm. Thus, our method can beseen as a post-processing step for Poisson denoising methods.The experimental results reported in this paper show that ourmethod improves over state-of-the-art denoising methods byrelevant margins.In the following sections, we first review the BLP principle;then, the proposed method is described. Finally, experimentalresults are reported and conclusions are drawn.II. B
EST L INEAR P REDICTION – BLPIn this paper, we use a uppercase normal font to indicaterandom vectors and uppercase bold font for deterministicmatrices. Consider a random vector X ∈ R n with unknownprobability density function is to be estimated from an ob-served random vector Y ∈ R m as an affine function of theform (cid:98) X = B Y + a , (1)where B ∈ R n × m and a ∈ R n are a fixed matrix andvector, respectively. Let µ y and µ x be the mean of the randomvariables Y and X , respectively. Based on the BLP theorem[12], the optimal choice for B and a in the MMSE sense iswell-known to be B = Σ xy Σ − yy , a = µ x − Σ xy Σ − yy µ y (2)where Σ xy and Σ yy are the cross-covariance matrix between X and Y , and the covariance matrix of Y , respectively. Aremarkable feature of BLP is that it is independent of the a r X i v : . [ c s . C V ] M a r istribution of the random variables involved, depending onlyon the aforementioned matrices and mean vectors. In thecase of Gaussian noise, it can easily be shown that BLPis equivalent to using a multivariate Gaussian prior withcovariance Σ xx , which is often used in image denosing [8],[13]; in that case, (2) is not only the best linear predictor, butthe best predictor overall in the MMSE sense [14].In this paper, we propose to use BLP (1)–(2) in the casewhere Y is a Poisson random vector with the underlying mean X ( E [ Y | X ] = X ), that is, for m = n and P ( Y = y | X = x ) = n (cid:89) j =1 e − x j x y j j y j ! , (3)where y ∈ N n and x ∈ R m + . The main challenge in thisapproach will be to obtain the aforementioned covariancematrices and mean vectors, as discussed in the next section.III. T HE PROPOSED METHOD
A. BLP from Poisson Observations
The Poisson conditional distribution of Y in (3) has animportant implication in the cross-covariance matrix Σ xy .Recall that the cross-covariance matrix is given by Σ xy = E X,Y [ XY T ] − µ x µ Ty , (4)where E X,Y indicates the expectation with respect to the jointdistribution of X and Y . Using the so-called law of iteratedexpectation , the above can be written as Σ xy = E X (cid:2) E Y | X [ XY T ] (cid:3) − µ x µ Ty = E X ( XX T ) − µ x µ Tx = Σ xx , (5)because µ y = µ x and E Y | X [ XY T ] = X E Y | X [ Y T ] = XX T .In conclusion (as is also true for any zero-mean additiveindependent noise), in the case of Poissonian observations(which is not additive noise), we have that Σ xy = Σ xx .A particular property of the Poisson distribution (namelythat its mean and variance are equal) underlies a simplerelationship between Σ yy and Σ xx . A similar relationshipexists in the case of additive independent zero-mean noise withknown covariance , but for the Poisson model (3), no furtherinformation is needed. Using iterated expectation again, Σ yy = E Y ( Y Y T ) − µ y µ Ty = E X (cid:2) E Y | X [ Y Y T ] (cid:3) − µ x µ Tx . (6)Conditioned on X , the components of Y are independent, thusthe off-diagonal elements of E Y | X [ Y Y T ] are simply ( i (cid:54) = j ) ⇒ (cid:0) E Y | X [ Y Y T ] (cid:1) i,j = E [ Y i ] E [ Y j ] = X i X j . Concerning the diagonal elements, (cid:0) E Y | X [ Y Y T ] (cid:1) i,i = E Y | X [ Y i ] = X i + X i , For example, for additive zero-mean white noise of known variance σ ,we have the well-known relationship Σ yy = Σ yy + σ I . because the mean and variance of a Poisson random variableare identical. Consequently, E Y | X [ Y Y T ] = XX T + diag ( X , ..., X n ) and E X (cid:2) E Y | X [ Y Y T ] (cid:3) = E X [ XX T ] + E X [ diag ( X , ..., X n )]= Σ xx + µ x µ Tx + diag ( µ x ) , (7)where diag ( µ x ) denotes a diagonal matrix with the compo-nents of µ x . Finally, plugging (7) into (6), yields Σ yy = Σ xx + diag ( µ x ) , (8)which is a natural result, since Poisson noise can be seenas additive independent noise with variance equal to theunderlying clean value.Finally, the BLP (1) of X from Poisson observations Y modeled by (3) is given by (cid:98) X = µ x + Σ xx ( diag ( µ x ) + Σ xx ) − ( Y − µ y ) . (9) B. Application to Patch-Based Poisson Denoising
As in classical patch-based methods, the image is dividedinto overlapping √ n ×√ n patches [3], [15], [16]. The patchesare denoised in collaboration with similar patches and thenreturned to the original position in the image, with the multipleestimates of each pixel (due to the overlapping nature of thepatches) averaged to reconstruct the final clean image estimate.Denoting the (vectorized) i -th noisy patch as y i , we estimatethe corresponding clean patch x i by BLP, i.e. , using (9), (cid:98) x i = µ x + Σ xx ( diag ( µ x ) + Σ xx ) − ( y i − µ y ) . (10)where µ x and Σ xx are estimated as describe in the nextparagraphs.Some methods for Gaussian denoising use the same un-derlying covariance matrix for similar patches in the image[8], [13]. These approaches can all be seen as inspired by thecollaborative filtering idea pioneered in the famous BM3Ddenoising method [16]. In this paper, we follow [8], byassuming the same covariance matrix Σ xx for patches thatare similar to a reference patch y r in a neighborhood thereof.Estimating µ x and Σ xx from a set of clean patches wouldsimply amount to computing the sample mean and sample co-variance. In the presence of additive zero-mean white Gaussiannoise of known variance σ , it would still be possible to obtainestimates µ x and Σ xx from noisy patches [17]. However, thisis not so simple in the case of Poissonian observations, becausethe variance of each observation depends on the underlyingclean value, which is of course unknown. Here, we proposeto use an off-the-shelf Poisson denoising methods to obtain aninitial estimate of the clean image, from which µ x and Σ xx can be estimated; these estimates are then plugged into (10)to obtain the final patch estimator. . Proposed Algorithm Let y denote the whole Poisson noisy image. The first stepis to use an off-the-shelf Poisson denoising method to obtain aso-called pilot estimate ˜ x . In the experiments reported below,we will use NL-PC [4], VST+BM3D, and the recent state-of-the-art method in [18].The remaining steps of the algorithm follow closely themacro structure of BM3D [16], [19]. Denoting a referencepatch of the pilot estimate as ˜ x r , this image is searched withina window of size N × N , centered at ˜ x r , for the k -nearestpatches in Euclidean distance; let X r denote the set of patchesthus obtained. From this set of patches, their sample mean andsample covariance are obtained and denote as µ x and Σ xx .Let Y r denote the set of patches in the noisy image atthe same locations as the patches in X r . Using µ x and Σ xx obtained as explained in the previous paragraph, all the patchesin Y r are denoised by BLP, according to (10). These denoisedpatches are then returned to their locations, and averagedwherever they overlap.The algorithm may be repeated L times, using the obtaineddenoised image as the next pilot estimate.It is clear from the description that, rather than a full self-contained denoising method, our algorithm can be seen asa post-processing step for other Poisson denoising methods.As shown below in the experiments, the proposed approachimproves over several state-of-the-art Poisson denoising algo-rithms by a non-negligible margin.IV. R ESULTS
In this section, we evaluate the performance of our algo-rithm using different methods to obtain the pilot estimate : NL-PCA [4], VST+BM3D [2], and a recent state-of-the-art methodin [18]. The parameters used in the proposed method are asfollows. Reference patches are selected with step-size alongboth the row and the columns of the image. The patches havesize × and the search window is × . The numberof selected similar patches to each reference patch is set to k = 30 . finally, we use L = 2 iterations of the algorithm.Table I, II, and III show the results of the proposed methodfor some benchmark images, in comparison to the methodsused to obtain the corresponding pilot estimates. It can be seenthat in all examples, the BLP improves the pilot methods forevery tested images. The improvements often exceed . dB.The average improvement of our method is also noticeable.In these tables, mainly the low SNR values with peak valuesfrom to are considered. Although the peak value of in images is sometimes not considered as a low SNR, it isincluded to show the capability of our method to improve thehigher SNR’s even when VST methods are quite accurate inthese ranges.Some examples of the proposed denoising methods usingdifferent initialization are shown and compared to the corre-sponding initialized images in Figures 1, 2, and 3. It can beseen that the proposed technique is able to reduce some ofthe artifacts produced by the underlying initialization method,yielding visually more pleasing results. (a)(b) (c) Fig. 1. Example of denoising results of a part of the image Cameraman (a)Noisy image (Peak= ), (b) VST+BM3D (PSNR=23.14) (c) Proposed methodinitialized by VST+BM3D (PSNR=23.51) (a)(b) (c) Fig. 2. Example of denoising of a part of the image Barbara, (a) Noisyimage (Peak= ), (b) VST+BM3D (PSNR=26.29), (c) VST+BM3D+BLP(PSNR=26.61)ABLE IR ESULTS OF THE PROPOSED METHOD IN COMPARISON WITH
NLPCApeak= peak= peak= NLPCA NLPCA+BLP NLPCA NLPCA+BLP NLPCA NLPCA+BLPCameraman 20.84
House 22.64
Barbara 21.47
Lena 23.64
Average 22.14
TABLE IIR
ESULTS OF OUR PROPOSED METHOD IN COMPARISON WITH
VST+BM3Dpeak= peak= peak= VST+BM3D VST+BM3D+BLP VST+BM3D VST+BM3D+BLP VST+BM3D VST+BM3D+BLPCameraman 22.04
House 24.04
Barbara 22.09
Lena 24.37
Average 23.13
TABLE IIIR
ESULTS OF OUR PROPOSED METHOD IN COMPARISON WITH THE METHOD IN [18]peak= peak= peak= peak= [18] [18]+BLP [18] [18]+BLP [18] [18]+BLP [18]+BLP [18]+BLPCameraman 21.17 House 23.56
Barbara 21.37
Lena 24.01
Average 22.53 (a)(b) (c)
Fig. 3. Example of denoising results of the House image, (a) Noisy image(peak= ), (b) NLPCA (PSNR= . ), (c) NLPCA+BLP (PSNR= . ) V. C
ONCLUSION
In this paper, we have proposed a new approach based on best linear prediction (BLP) for denoising images contami-nated with Poisson noise. BLP is independent of the noisedistribution and we showed that, in the Poisson noise case,the only required information is the mean and covariancematrix of the underlying vector being estimated. To imple-ment the BLP approach, our method relies on an off-the-shelf Poisson denoising method in order to obtain a pilotestimate, from which the required means and covariances areestimated. Consequently, the proposed method can be seenas a post-processing step for Poisson denoising methods. Theexperiments reported show that our method is able to obtaina noticable improvement over the initial denoised imagesobtained by several state-of-the-art methods.R
EFERENCES[1] F. Anscombe, “The transformation of Poisson, binomial and negative-binomial data,”
Biometrika , vol. 35, no. 3/4, pp. 246–254, 1948.[2] M. M¨akitalo and A. Foi, “Optimal inversion of the Anscombe transfor-mation in low-count Poisson image denoising,”
IEEE Transactions onImage Processing , vol. 20, no. 1, pp. 99–109, 2011.[3] C. Deledalle, F. Tupin, and L. Denis, “Poisson nl means: Unsupervisednon local means for poisson noise,” in
Image processing (ICIP), 201017th IEEE int. conf. on . IEEE, 2010, pp. 801–804.[4] J. Salmon, Z. Harmany, C. Deledalle, and R. Willett, “Poisson noisereduction with non-local PCA,”
Journal of Mathematical Imaging andVision , vol. 48, no. 2, pp. 279–294, 2014.[5] M.A.T. Figueiredo and J.M. Bioucas-Dias, “Restoration of Poissonianimages using alternating direction optimization,”
IEEE Transactions onImage Processing , vol. 19, pp. 3133 –3145, 2010.6] Z. Harmany, R. Marcia, and R. Willett, “This is SPIRAL-TAP: SparsePoisson intensity reconstruction algorithmstheory and practice,”
IEEETransactions on Image Processing , vol. 21, pp. 1084–1096, 2012.[7] S. Gu, L. Zhang, W. Zuo, and X. Feng, “Weighted nuclear normminimization with application to image denoising,” in
Proceedings of theIEEE Conference on Computer Vision and Pattern Recognition , 2014,pp. 2862–2869.[8] M. Niknejad, H. Rabbani, and M. Babaie-Zadeh, “Image restorationusing Gaussian mixture models with spatially constrained patch cluster-ing,”
IEEE Transactions on Image Processing , vol. 24, pp. 3624–3636,2015.[9] A. Rond, R. Giryes, and M. Elad, “Poisson inverse problems by theplug-and-play scheme,”
Journal of Visual Communication and ImageRepresentation , vol. 41, pp. 96–108, 2016.[10] S. Pyatykh and J. Hesser, “MMSE estimation for Poisson noise removalin images,” arXiv:1512.00717 , 2015.[11] M. Niknejad, J. Bioucas-Dias, and M. A. T. Figueiredo, “Class-specificPoisson denoising by patch-based importance sampling,” in
IEEE Intern.Conf. Image Processing (ICIP) , 2017.[12] S. Searle, G. Casella, and C. McCulloch,
Variance components , JohnWiley & Sons, 2009.[13] G. Yu, G. Sapiro, and S. Mallat, “Solving inverse problems withpiecewise linear estimators: From Gaussian mixture models to structuredsparsity,”
IEEE Transactions on Image Processing , vol. 21, pp. 2481–2499, 2012.[14] L. Scharf and C. Demeure,
Statistical signal processing: detection,estimation, and time series analysis , Addison-Wesley, 1991.[15] Antoni Buades, Bartomeu Coll, and J-M Morel, “A non-local algorithmfor image denoising,” in , 2005, vol. 2, pp. 60–65.[16] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising bysparse 3-D transform-domain collaborative filtering,”
IEEE Transactionson image processing , vol. 16, pp. 2080–2095, 2007.[17] A. Teodoro, M. Almeida, and M. Figueiredo, “Single-frame imagedenoising and inpainting using Gaussian mixtures,” in , 2015.[18] L. Azzari and A. Foi, “Variance stabilization for noisy+estimatecombination in iterative Poisson denoising,”
IEEE Signal ProcessingPetters , vol. 23, pp. 1086–1090, 2016.[19] M. Lebrun, “An analysis and implementation of the BM3D imagedenoising method,”