A Non-Local Means Filter for Removing the Poisson Noise

Noname man uscript No. (will be inserted b y the editor) A Non-Lo cal Means Filter for Remo ving the P oisson Noise Qiyu JIN · Ion Grama · Quansheng Liu the date of receipt a nd acceptanc e should be inserted later Abstract A new image denoising a lgorithm to deal with the Poisson noise mo del is given, which is based on the idea of No n- Loca l Mea n. By using the ”Orac le ” co ncept, w e esta blish a theorem to show that the Non- L o ca l Means Filter can effectively deal with P oisso n noise with some mo difica tion. Under the theoretical result, we construct our new algor ithm called Non-Lo cal Me a ns Poisson Filter and demo nstrate in theor y that the filter conv erges at the usual optimal rate. The filter is as simple as the c la ssic Non-Local Means and the simulation r esults sho w that our filter is very comp etitive. Keyw ords Non-Lo cal Mea ns · Mean Sq uare Er ror · Poisson noise · ”O r acle” estimator Qiyu JIN Institute of Image Pr ocessing and Pat tern Recognition, Shanghai Jiao T ong Universit y , No. 800 Dongch uan Road, Mi nhang District, Shanghai 200240, China UMR 6205, Lab oratoire de Mathmatiques de Br etagne Atlan tique, Universit´ e de Bretagne- Sud, Campus de T ohaninic, BP 573, 56017 V annes, F rance Unive rsit´ e Europ´ eenne de Bretagne, F rance Jiangsu Engineering Cen ter of Netw ork Monitoring, Nanjing Unive rsity of Information Sci- ence & T ech nology , Nanjing 210044, China E-mail: qiyu. j in2008@gmail.com Ion Grama UMR 6205, Lab oratoire de Mathmatiques de Br etagne Atlan tique, Universit´ e de Bretagne- Sud, Campus de T ohaninic, BP 573, 56017 V annes, F rance Unive rsit´ e Europ´ eenne de Bretagne, F rance E-mail: i on.grama@univ-ubs.fr Quansheng Li u UMR 6205, Lab oratoire de Mathmatiques de Br etagne Atlan tique, Universit´ e de Bretagne- Sud, Campus de T ohaninic, BP 573, 56017 V annes, F rance Unive rsit´ e Europ´ eenne de Bretagne, F rance Sc hool of Mathematics and Computing Sciences, Changsha Uni versity of Science and T ech - nology , Changsha 410076, China E-mail: quansheng.liu@univ-ubs.f r 2 Qiyu JIN et al. 1 In tro duction Noise is inevitable in an y image device. A digital imaging sys tem consists of an optical sys tem follow ed by a photo detector and asso ciated electrical filters. The photo detector conv e rts the incident optical in tensity to a detector current, i.e. photons to electrons. During the pro cess, the true sig nals are contaminated b y many differ e n t sour c es of noise. The Poisson noise appea rs in low-light conditions w he n the n umber o f collec ted photons is s mall, such as night v is ion, medica l imaging, underwater ima g ing, micr oscopic ima ging, optical microsco py imaging and astrono m y imaging . Such a noise is signal- depe ndent, and requires to adapt the usual denoising approa c hes. The key challenge in Poisso n int ensity estimation pro blems is that the v ar i- ances of the observed counts are different. As a res ult, ma n y metho ds are in- tro duced to transfor m the Poisso n distributed nois e to the data approximately Gaussian and homo scedastic. These metho ds ar e called V ariance Stabilizing T ransformations (VST), suc h as Anscombe ro ot tra nsformation (1 948 [4 ], and 1993 [6]), multiscal VSTs (200 8 [3 9]), conditional v a riance stabilizatio n (CVS) (2006 [18]), or Haa r-Fisz transformation (2 004 [15] and 2007 [14]). T he n we can dea l with these data as Gaussian nois e. Second, the noise is removed using a co n ven tio nal denoising algorithm for additiv e white Gaussian no ise, see for example B ua des, Coll and Morel (2 005 [8]), Ker vrann (2006 [22]), Aharon and Elad a nd Bruckstein (20 06 [2]), Hammond and Simoncelli (2008 [16]), Polzehl and Sp okoin y (2006 [34]), Hirak aw a a nd Parks (20 0 6 [17]), Mairal, Sapir o and Elad (2008 [2 7]), Portilla, Stre la , W ainwrigh t and Simoncelli (200 3 [35]), Roth and Blac k (2009 [36]), Katk ovnik, F oi, E giazaria n, and Astola (20 10 [2 1]), Dabov, F oi, K atk ovnik and E giazaria n (2006 [9]), Abraham, Abraham, Deso l- neux and Li-Thiao-T e (2007 [1]), and Jin, Grama and Liu (2011 [20]). Aft er denoising, s ome inv erse tra nsformations, like Exa ct Un biased In verse (EUI) (2009 [28] a nd 201 1 [29]), are a pplied to the denoised signal, obtaining the estimate of the s ignal of interest. Man y authors restore the Poisson nois e b y this type of methods with a three - step pro cedure (see [7, 25, 26, 3 9] ). Maxmum Likeliho od (ML) estimation (1996 [38], 2009[32]) and Similarity Measure (SM) (2006 [3 ]) are als o found to b e effectiv e since they can a ccount for the sp ecial proper ties of the Poisso n distribution. O thers methods lik e as Complexity-P enalized Likeliho od Estimation (CPLE) (2000 [33], 20 05 [23]) and T otal V ariation (TV) seminorm (20 09 [5]), have been in tro duced to deal with the Poisson noise. Le et al. ([24]) have adapted the success ful ROF mo del for total v ar iation regularizatio n to deal with Poisson noise. The gradient de- scent iteration for this mo del r eplaces the regular ization pa rameter with a function. The Non-Lo cal Mea ns Filter has b een propo sed by B uades et al (2 0 05 [8]) to denoise images damage d by additiv e white Gaussian noise. It is based o n the similarity phenomenon exis ting very often in natural images, and assumes that there is enough r edundan t information (pixels having identical noise-fre e v alue) in the image to r educe the noise significantly . This filter is known to efficiently reduce the no ise a nd to preserve structure s . Some a uthors (see 2008 A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 3 [7], 2010 [10]) co m bine the No n-Lo c al Mea ns metho d with other methods to restore the Poisson noise. Deledalle et al. (2010 [1 0]) prop osed an ex tens io n of the Non-Lo cal Means for images dama ged by Poisson noise. It is based on pro babilistic s imila rities to compar e nois y pa tches and patches of a pre- estimated image. In this pap er, a new image denoising alg orithm to deal with the Poisson noise model is given, which is based o n the idea of Non-Lo cal Mean. Our main idea is as follows: we first obtain an ”Oracle ” estimator by minimized a very tigh t upper b ound of the Mean Square Err or with c hanging the size of search window. The ”Oracle” estimator depends on the unknown ta rget func- tion (or ig inal image), whose co ncept is developed in Donoho and Johnstone [11]. So the ”Oracle” estimator is not computable, but it can help us to find an a v ailable a lgorithm in mathematic theory . W e second es tablish a theorem by the concept of the ”Oracle” to show that the Non-Lo cal Means Filter can effectively deal with Poisson noise with some mo dification. Finally , replacing the unknown ta rget function by some es timators, we construct our new a l- gorithm called No n-Loc al Means Poisson Filter and demons tr ate in statistic theory that the filter converges at the usua l optimal r a te. The filter is as simple as the cla ssic No n-Loc al Mea ns a nd the simulation r esults show that our filter is v ery comp etitive. The remainder of this pap er is orga nized a s follow: w e first introduce an ”Orac le ” estimator for P oisso n noise based o n the idea of Non-Lo cal Means, and pre s en t a theo r em to show the r ate of conv ergence of the ”Or acle” es tima- tor in Section 2. W e second constr uc t an ada ptiv e estimator according to the ”Orac le ” estimator a nd obtain so me convergence theorems o f the estimato r in Section 3. Finally , we demonstr ate in Sectio n 4 the ability of a pproach at restoring image contaminated b y P oisson noise with a brief analysis. 2 The ”Oracle” estimator 2.1 Some notations W e suppo se that the orig inal imag e of the ob ject be ing photogr aphed is a int egra ble tw o-dimensional function f ( x ), x ∈ (0 , 1] × (0 , 1]. L e t the mean v alue of f in a set B x be Λ ( B x ) = N 2 Z B x f ( t ) dt. Typically w e obser v e a discr ete dataset of counts Y = {N ( B x ) } , wher e N ( B x ) is a Poisson random v a riable of int ensity Λ ( B x ). W e consider that if B x ∩ B y = ∅ , then N ( B x ) is independent of N ( B y ). Supp o se that x = ( x (1) , x (2) ) ∈ I = { 1 N , 1 N , · · · , 1 } 2 , and B x = ( x (1) − 1 / N , x (1) ] × ( x (2) − 1 / N , x (2) ]. Then { B x } x ∈ I is a partition of the square (0 , 1] × (0 , 1]. Using this partition w e get a disc rete function f ( x ) = Λ ( B x ), x ∈ I . The denoising algor ithm aims 4 Qiyu JIN et al. at estimating the underlying intensit y profile discrete function f ( x ) = Λ ( B x ). The image function f is considered to b e co nstan t on ea ch B x , x ∈ I . Ther efore f ( x ) = N ( B x ), x ∈ I . F urthermore , w e can estimate the integrable function p by the discre te function f . Let Y ( x ) = N ( B x ) , x ∈ I . (1) This model has b een used effectiv ely in many co n texts. The Poisson noise mo del can be re wr itten in the regression form Y ( x ) = f ( x ) + ǫ ( x ) , x ∈ I , (2) where ǫ ( x ) = N ( B x ) − f ( x ). It is easy to see that E ( ǫ ( x )) = 0 and V ar ( ǫ ( x )) = f ( x ). Let us set some notations to b e us ed throughout the pap er. The Euclidea n norm of a vector x = ( x 1 , ..., x d ) ∈ R d is denoted by k x k 2 =  P d i =1 x 2 i  1 2 . The supremum norm of x is denoted by k x k ∞ = sup 1 ≤ i ≤ d | x i | . The cardinalit y of a set A is denoted card A . F or a po sitiv e integer N the uniform N × N gr id of pixels on the unit square is defined b y I =  1 N , 2 N , · · · , N − 1 N , 1  2 . (3) Each ele men t x of the grid I will b e called pixel. The n umber of pixels is n = N 2 . F or any pixel x 0 ∈ I a nd a g iv en h > 0 , the square window of pixels U x 0 ,h = { x ∈ I : k x − x 0 k ∞ ≤ h } (4) will b e called se ar ch window at x 0 . W e naturally take h as a m ultiple o f 1 N ( h = k N for some k ∈ { 1 , 2 , · · · , N } ). The size of the squa r e sear c h window U x 0 ,h is the pos itive integer num ber M = (2 N h + 1) 2 = ca rd U x 0 ,h . (5) F o r any pixel x ∈ U x 0 ,h and a given η > 0 a second square window of pixels U x,η will b e called p atch at x . Like h , the parameter η is also ta k en as a m ultiple of 1 N . The size of the patch U x,η is the pos itive integer m = (2 N η + 1) 2 = ca rd U x 0 ,η . (6) 2.2 The Non-Lo cal Means algorithm The Non-Lo cal Mea ns algo r ithm (2005 [8]) can b e describ ed as fo llows. F or any x ∈ I , e f = X x ∈ I w ( x ) Y ( x ) , (7) A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 5 where the weigh ts w ( x ) are given by w ( x ) = e − e ρ 2 x 0 ( x ) /H 2  X x ′ ∈ I e − e ρ 2 x 0 ( x ′ ) /H 2 , (8) with e ρ 2 x 0 = X y ∈ U x 0 ,η κ ( y ) | Y ( y ) − Y ( T x y ) | 2 P y ′ ∈ U x 0 ,η κ ( y ′ ) . Here H is a bandwidth parameter, U x 0 ,η is given by (4), κ ( y ) > 0 are so me fixed kernel, and T x is the translation mapping: T x : x 0 + y → x + y (9) In prac tice the bandwidth parameter H is often taken as a linear funct ion of σ (see [8]). 2.3 Oracle estimator In o rder to ada pt the Non-Lo cal Mea ns algor ithm to the Poisson noise, we int ro duce an ”Ora cle” estimator (for details on this concept see Donoho and Johnstone (1994 [11])). Denote f ∗ h = X x ∈ U x 0 ,h w ∗ h Y ( x ) , (10) where w ∗ h ( x ) = e − ρ 2 f,x 0 ( x ) H 2 ( x 0 )  X x ′ ∈ U x 0 ,h e − ρ 2 f,x 0 ( x ′ ) H 2 ( x 0 ) (11) with ρ f ,x 0 ( x ) ≡ | f ( x ) − f ( x 0 ) | , (12) H ( x ) is a control function sub ject to γ = inf { H ( x ) : x ∈ I } > 0 . (13) It is ob vious that X x ∈ U x 0 ,h w ∗ h ( x ) = 1 and w ∗ h ( x ) ≥ 0 . (14) Note that the function ρ f ,x 0 ( x ) ≥ 0 c haracterizes the similar it y of the image brightness at the pixel x with r espect to the pixel x 0 , therefore w e shall call 6 Qiyu JIN et al. ρ f ,x 0 similarity function. The usual bias-v ar iance decomp osition (cf. e.g. [12, 30, 37]) of the Mean Squared Error (MSE) E ( f ( x 0 ) − f ∗ h ( x 0 )) 2 (15) =   X x ∈ U x 0 ,h w ∗ h ( x ) ( f ( x ) − f ( x 0 ))   2 + X x ∈ U x 0 ,h w ∗ h ( x ) 2 f ( x ) ≤   X x ∈ U x 0 ,h w ∗ h ( x ) | f ( x ) − f ( x 0 ) |   2 + X x ∈ U x 0 ,h w ∗ h ( x ) 2 f ( x ) . The inequality (1 5) com bining with (12) implies the following upp er b ound E ( f ( x 0 ) − f ∗ h ( x 0 )) 2 ≤ g ( w ∗ h ( x )) , (16) where g ( w ) =   X x ∈ U x 0 ,h w ( x ) ρ f ,x 0 ( x )   2 + X x ∈ U x 0 ,h w ( x ) 2 f ( x ) . (17) W e sha ll define a family of es timates by minimizing the function g ( w h ) by changing the width of the sea r c h window. With a Poisson noise in low-light co n- ditions, the upp er bo und of signal function is small, so we let Γ = sup { f ( x ) : x ∈ I } . Acco rding to the similarity phenomenon existing very o ften in natural images, we supp ose that the function f satisfies the lo cal H¨ older condition | f ( x ) − f ( y ) | ≤ L k x − y k β ∞ , ∀ x, y ∈ U x 0 ,h + η , (18) where β > 0 and L > 0 are constants, h > 0, η > 0 and x 0 ∈ I . The following theorem gives the r ate o f conv ergence of the ”Or acle” estima to r and the pr oper width h of the search window. Theorem 1 Assume that h =  Γ 4 β L 2  1 2 β + 2 n − 1 2 β + 2 and γ > √ 2 Lh β . Supp ose that the fun ction f satisfies t he lo c al H¨ older c ondition (18) and f ∗ h ( x 0 ) b e given by (10) . Then E ( f ∗ h ( x 0 ) − f ( x 0 )) 2 ≤ c 0 n − 2 β 2 β + 2 , (19) wher e c 0 = 2 2 β + 6 2 β + 2 Γ 2 β 2 β + 2 L 4 2 β + 2 β 2 β 2 β + 2 . (20) F o r the pr oo f of this theorem see Section 6.1. This theorem shows that at least from the pr actical p oint of view, it is justified to optimize the upp er b ound g ( w ) instea d of o ptimizing the ris k E ( f ∗ h ( x 0 ) − f ( x 0 )) 2 itself. The theorem also justifies that we can choose a small search window in pla c e of the whole o bserved ima ge to estimate a point, without loss of visua l qualit y . That is wh y w e only co nsider small s earch win- dows for the simulations of our algorithm. A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 7 3 Non-Lo cal Means P o isson Filter 3.1 Construction of Non-Lo cal Means Poisson Filter With the theory of ”Ora c le ” es timator, we c o nstruct the No n- Loca l Means Poisson Filter. Let h > 0 and η > 0 be fixed num b ers. Since | f ( x ) − f ( x 0 ) | 2 = E | Y ( x ) − Y ( x 0 ) | 2 − ( f ( x 0 ) + f ( x )), a n ob vious estimator o f E | Y ( x ) − Y ( x 0 ) | 2 is given by 1 M X y ∈ U x 0 ,η | Y ( y ) − Y ( T x y ) | 2 , where T x is given by (9), and ( f ( x 0 ) + f ( x )) is estimated by 2 f ( x 0 ), where f ( x 0 ) = 1 M X x ∈ U x 0 ,h f ( x ) . Define an estimated similarity fun ction b ρ x 0 by b ρ 2 x 0 ( x ) =   1 M X y ∈ U x 0 ,η | Y ( y ) − Y ( T x y ) | 2 − 2 f ( x 0 )   + . (21) The following theorem implies that it is r easonable to let b ρ x 0 ( x ) b e the estimator of ρ f ,x 0 ( x ). Theorem 2 Assume t hat h =  Γ 4 β L 2  1 2 β + 2 n − 1 2 β + 2 and η = c 1 n − α  (1 − β ) + 2 β +2 < α < 1 2  . Supp ose that the function f s atisfi es t he lo c al H ¨ older c ondition (18) and b ρ 2 x 0 ( x ) is given by ( 21). Then ther e i s a c onstant c 2 such tha t P  max x ∈ U x 0 ,h   b ρ 2 x 0 ( x ) − ρ 2 f ,x 0 ( x )   ≥ c 2 n α − 1 2 √ ln n  ≤ O  n − 1  . (22) F o r the pr oo f of this theorem see Section 6.2. As a result, it is natural to define an adaptive estimator b f h by b f h ( x 0 ) = X x ∈ U x 0 ,h b w h ( x ) Y ( x ) , (23) where b w h = e − b ρ 2 x 0 ( x ) H 2  X x ′ ∈ U x 0 ,h e − b ρ 2 x 0 ( x ′ ) H 2 . (24) and U x 0 ,h given by (4). 8 Qiyu JIN et al. 3.2 Conv ergence theorem of Non-Loca l Means Poisson Filter Now, we turn to the study of the conv ergence o f the Non-Lo cal Means Poisson Filter. Due to the difficulty in dealing with the dep endence of the weigh ts we shall consider a slightly mo dified version of the prop osed algor ithm: we divide the set of pixels into tw o indep enden t parts, so that the weigh ts are constructed f rom the one par t, and the estimation of the target function is a weigh ted mean along the other part. More precisely , we split the set of pixe ls int o tw o par ts I = I ′ x 0 ∪ I ′′ x 0 for any x 0 ∈ I wher e I ′ x 0 =  x 0 +  i N , j N  ∈ I : i + j is even  , and I ′′ x 0 = I  I ′ x 0 . Define an estimated similarity function b ρ x 0 by b ρ ′ 2 x 0 ( x ) =   1 card U ′′ x 0 ,η X y ∈ U ′′ x 0 ,η | Y ( y ) − Y ( T x y ) | 2 − 2 f ′ ( x 0 )   + , x ∈ U ′ x 0 ,h (25) where f ′ ( x ) = 1 card U ′′ x 0 ,h X y ∈ U ′′ x 0 ,h Y ( y ) , (26) and U ′′ x 0 ,η = U x 0 ,η ∩ I ′′ x 0 with U x 0 ,h given b y (4). The adaptive estimator b f ′ h is denoted b y b f ′ h ( x 0 ) = X x ∈ U ′ x 0 ,h b w ′ h ( x ) Y ( x ) , (27) where U ′ x 0 ,h = U x 0 ,h ∩ I ′ x 0 and b w ′ h = e − b ρ ′ 2 x 0 ( x ) H 2 ( X 0 )  X x ′ ∈ U ′ x 0 ,h e − b ρ ′ 2 x 0 ( x ′ ) H 2 ( x 0 ) . (28) In the next theorem we prove that the Mea n Squar ed E r ror of the esti- mator b f ′ h ( x 0 ) conv erges at the rate n − 2 β 2 β + 2 which is the usual optimal r ate of conv ergence for a given H¨ older smo othness β > 0 (s e e e .g. F an a nd Gijb els (1996 [13])). Theorem 3 L et η = c 3 n − α , h =  Γ 4 β L 2  1 2 β + 2 n − 1 2 β + 2 , H ( x 0 ) > 4 c 2 n α − 1 2 √ ln n and γ > max { √ 2 Lh β , 4 c 2 n α − 1 2 √ ln n } . Supp ose tha t the fu n ction f satisfies the H¨ older c ondition (18) and b f ′ h ( x 0 ) is given by (25). Then E  b f ′ h ( x 0 ) − f ( x 0 )  2 ≤ c 4 n − 2 β 2 β + 2 , (29) A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 9 wher e c 4 = 8 2 2 β + 6 2 β + 2 Γ 2 β 2 β + 2 L 4 2 β + 2 β 2 β 2 β + 2 ! 2 . F o r the pr oo f of this theorem see Section 6.2. 4 Sim ulation results 4.1 Computational algor ithm Throughout the simulations, w e use the following algorithm to compute the Non-Lo cal Mea ns Poisson estimator b f h ( x 0 ) . The input v alues o f the alg orithm are Y ( x ) , x ∈ I (the image), and tw o num be r s m and M ar e given by (6) and (5) res pectively . In order to improv e the r esults, we int ro duce a smo othed version of the es timated similarit y distance b ρ 2 κ,x 0 ( x ) =    X y ∈ U x 0 ,η κ ( y ) | Y ( y ) − Y ( T x y ) | 2 − 2 f ( x 0 ) P y ′ ∈ U x 0 ,η κ ( y ′ )    + . (30) As smo othing k ernels κ we use the Gaussia n k ernel κ g ( y , h g ) = exp  − N 2 k y − x 0 k 2 2 2 h g  , (31) where h g is a constant, a nd the follo wing k ernel: for y ∈ U x 0 ,η , κ 0 ( y ) = N η X k =max(1 ,j ) 1 (2 k + 1) 2 , (32) if k y − x 0 k ∞ = j N for so me j ∈ { 0 , 1 , · · · , N η } . W e shall als o use the r ectangular kernel κ r ( y ) =  1 card U x 0 ,η , y ∈ U x 0 ,η , 0 , otherwise . (33) F o r the simulation we use the k ernel κ 0 ( y ) defined by (32 ). W e ha ve seen exp erimen tally that when we ta k e the filtering function H 2 ( x 0 ) as µ · q f ( x 0 ), where µ is a co ns tan t dep ending o n the character of the image, to obtain a denoising of hig h visual quality . W e mention that throug hout the pape r we symmetrize images near the bo undary . Algorithm Non-Local Means Poisson Filter (NLMPF) Let { M , m, h g } be the pa rameters. Repe a t for e ac h x 0 ∈ I 10 Qiyu JIN et al. NMISE=0.1009 NMISE=0.0593 NMISE=0.0368 NMISE=0.1150 NMISE=0.0870 M = 19 × 19 M = 13 × 13 M = 9 × 9 M = 15 × 15 M = 7 × 7 m = 13 × 13 m = 3 × 3 m = 21 × 21 m = 21 × 21 m = 13 × 13 (a) Sp ots (b) Galaxy (c) Ridges (d) Barbara (e) Cell s Fig. 1 These i mages r estored by the first step of our algorithm. - compute Step 1 b w ( x ) = exp ( − b ρ 2 κ,x 0 ( x ) /H 2 ( x 0 )) b f 1 ( x 0 ) = P x ∈ U x 0 ,h b w ( x ) Y ( x )  P x ∈ U x 0 ,h w ( x ) Step 2 If 1 (2 d +1) 1 P k x − x 0 k≤ d/ N b f 1 ( x ) < δ compute b f ( x 0 ) = P k x − x 0 k≤ d/ N κ g ( x, h g ) b f 1 ( x ) / P k x − x 0 k≤ d/ N κ g ( x, h g ) else b f ( x 0 ) = b f 1 ( x 0 ) . Note: we take δ = 15. 4.2 Numerical pe r formance of the Non-Lo cal Mea ns P oisson Filter By sim ulations we found that the images with brig h tness b et ween 0 and 255 (like Bar bara) are w ell denoised b y the first step, but for the low co un t lev els images (with br igh tness less than µ ), the re s tored images by NLMPF are not smooth enough (see Figure 1 ). This explains why f or the low co un t le vel images, we smo oth the restor ed imag e s by step 2. Our experiments are done in th e same wa y a s in [39] and [28] to produce comparable results; we also use the same set o f test ima ges (all of 25 6 × 256 in size): Sp o ts [0 . 08 , 4 . 99], Galaxy [0 , 5], Ridges [0 . 05 , 0 . 85], Ba rbara [0 . 93 , 1 5 . 73], and Cells [0 . 53 , 16 . 93]. The authors of [39] and [2 8] kindly pr o vided us with their prog rams and the tes t ima ges. A matlab implementation of the a lgo- rithms deriv ed in this paper is a v aila ble online 1 . This unoptimized implemen- tation pro cesses the set of 256 × 256 test images 145 seconds with a sea rch window of s ize 15 × 1 5 and patc hes of size 21 × 21, 52 seconds with a searc h window o f size 9 × 9 and patc hes of size 21 × 21 . The computational time is of ab out 10s p er iteration on a 256 × 25 6 ima ge a nd Matlab o n an Intel Pent ium Dual CPU T3200 32-bit @ 2.00GHz CPU 3.00GHz. 1 h ttp://www.pami.sjtu.edu.cn/people/jinqy/ A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 11 T able 1 A comparison of the denoising perf ormance (NMISE) of sev eral denoising algo- rithms. Algorithm Our OW PNF Poisson EUI+ MS-VST MS-VST algorithm NLM BM3D +7 / 9 +B3 Spots[0 . 08 , 4 . 99] 0 . 0260 0 . 0259 0.0790 0 . 0358 0 . 0602 0 . 0810 Galaxy[0 , 5] 0 . 0284 0 . 0285 0.0346 0 . 0297 0 . 0357 0 . 0338 Ridges[0 . 05 , 0 . 85] 0 . 0140 0 . 0162 0.0154 0 . 0121 0 . 0193 0 . 0416 Barbara[0 . 93 , 15 . 73] 0 . 1150 0 . 1061 0.1207 0 . 0863 0 . 2391 0 . 3777 Cells[0 . 53 , 16 . 93] 0 . 0785 0 . 0794 0.0816 0 . 0643 0 . 0909 0 . 1487 T able 2 A comparison of the denoising p erformance (PSNR, DB ) of sev eral denoising algorithms. Algorithm Our OW PNF Poisson EUI+ MS-VST MS-VST algorithm NLM BM3D +7 / 9 +B3 Spots[0 . 08 , 4 . 99] 31 . 45 31 . 31 3 1.12 31 . 95 31 . 64 30 . 35 Galaxy[0 , 5] 28 . 09 27 . 80 27.77 28 . 04 27 . 57 27 . 94 Ridges[0 . 05 , 0 . 85] 24 . 69 23 . 90 24.94 25 . 89 2 4 . 49 24 . 03 Barbara[0 . 93 , 15 . 73] 24 . 71 24 . 60 24.72 25 . 92 21 . 81 20 . 22 Cells[0 . 53 , 16 . 93] 29 . 08 29 . 91 29.40 30 . 18 28 . 87 26 . 69 T able 1 shows the NMISE v alues of images r e constructed b y NLMPF, OWPNF [19], Poisson NLM [10], EUI+BM3D [2 9], MS-VST+ 7 / 9 [39] and MS- VST+B3 [39]. Our algor ithm reach the best in the ca se of Ga laxy[0 , 5], while OWPNF reac h the b est in the case of Spots[0 . 08 , 4 . 99]; for Ridges[0 . 05 , 0 . 85], Barbar a[0 . 93 , 15 . 7 3], and Cells[0 . 53 , 16 . 93 ], the method EUI+BM3D g iv es the bes t results, but our metho d is a lso very c o mpetitive. T able 2 shows the PSNR v alues of imag es reconstructed. O ur algor ithm a lso reach the b est in the c ase of Galaxy[0 , 5 ]. The metho d E UI+BM3D have the hig hes t P SNR v alue. How e ver, the most imp ortant ev aluation criteria is the visual qualit y of r e stored imag e . Figures 2- 6 illustr ate the visual quality of these deno ised images. It is obvious that The visua l quality o f the outputs of our metho d hav e high visual qua lit y and many details Remained. F or example, in the case of restor ed images o f Spo ts (cf. Figures 2), our algorithm and OWPNF remain mos t sp ots. W e can see clearly 7 spots at the third column (from left) in Figures 2 (c), while EUI+BM3D just rema ins 4 s pots, Poisson NLM Makes sev eral spots stic k ing together, the images restored by MS-VST + 7 /9 and MS-VST + B3 are not smo oth enough. In the case of Galaxy (cf. Fig ures 3), vis ually , our algorithm bes t preser v es the fine textures. In the other case, o ur metho d also lead to go o d result visually . 12 Qiyu JIN et al. 5 Conclusion In this paper , we ha ve presen t a new image denois ing algorithm to deal with the Poisson noise mo del, whic h is based on the idea of Non-Lo cal Mean. The ”Orac le ” estimator is obtained b y minimized a very tight upper bo und of the Mean Squar e Err o r with changing the size of sear c h window. It help to establis h a theorem to sho w that the Non-Lo c a l Means Filter can effectiv ely deal with Poisson noise with some mo dification. As a result, we successfully cons truct the new algor ithm ca lle d Non-Lo cal Mea ns Poisson Filter and demonstr ate in statistic theory that the filter conv erges at the usual o ptimal ra te. The filter is as s imple as the classic Non-Lo cal Mea ns and the simulation r esults show that our filter is v ery competitive. The idea of how to construct an algorithm for Poisson noise model is creative. With our idea, many algorithms to Remov e Gaussian Noise could deal with the Poisson nois e with some modification. 6 App endix: Pro ofs of the main results 6.1 Pro of of Theorem 1 Denoting for brevity I 1 =   X x ∈ U x 0 ,h w ∗ h ( x ) ρ f ,x 0 ( x )   2 =        X k x − x 0 k ∞ ≤ h e − ρ 2 f,x 0 ( x ) H 2 ( x 0 ) ρ f ,x 0 ( x ) X k x − x 0 k ∞ ≤ h e − ρ 2 f,x 0 ( x ) H 2 ( x 0 )        2 , (34) and I 2 = f ( x 0 ) X x ∈ U x 0 ,h ( w ∗ h ( x )) 2 = f ( x 0 ) X k x − x 0 k ∞ ≤ h e − 2 ρ 2 f,x 0 ( x ) H 2 ( x 0 )   X k x − x 0 k ∞ ≤ h e − ρ 2 f,x 0 ( x ) H 2 ( x 0 )   2 , (35) then we hav e g ( w ∗ h ) = I 1 + I 2 . (36) The conditions (13) and γ > √ 2 Lh β imply that for x ∈ U x 0 ,h , w e ha ve L 2 k x − x 0 k 2 β ∞ H 2 ( x ) ≤ L 2 h 2 β γ 2 ≤ 1 2 . (37) A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 13 Noting that e − t 2 H 2 ( x 0 ) , t ∈ [0 , γ / √ 2) is decrea sing, and using one term T aylor expansion, the inequality (37) implies that X k x − x 0 k ∞ ≤ h e − ρ 2 f,x 0 ( x ) H 2 ( x 0 ) ≥ X k x − x 0 k ∞ ≤ h e − L 2 k x − x 0 k 2 β ∞ H 2 ( x 0 ) ≥ X k x − x 0 k ∞ ≤ h  1 − L 2 k x − x 0 k 2 β ∞ H 2 ( x 0 )  ≥ 2 h 2 n. (38) Considering that te − t 2 H 2 ( x 0 ) , t ∈ [0 , γ / √ 2) is increasing function, X k x − x 0 k ∞ ≤ h e − ρ 2 f,x 0 ( x ) H 2 ( x 0 ) ρ f ,x 0 ( x ) ≤ X k x − x 0 k ∞ ≤ h L k x − x 0 k β ∞ e − L 2 k x − x 0 k 2 β ∞ H 2 ( x 0 ) ≤ X k x − x 0 k ∞ ≤ h L k x − x 0 k β ∞ ≤ 4 Lh β +2 n. (39) The ab ov e three inequalities (34), (38) and (39) imply that I 1 ≤ 4 L 2 h 2 β . (40) T aking in to account the inequality X k x − x 0 k ∞ ≤ h e − 2 ρ 2 f,x 0 ( x ) H 2 ( x 0 ) ≤ X k x − x 0 k ∞ ≤ h 1 = 4 h 2 n, (35) and (38), it is easily seen that I 2 ≤ Γ h 2 n . (41) Combining (36), (40 ), a nd (41), we give g ( w ∗ h ) ≤ 4 L 2 h 2 β + Γ h 2 n . (42) Let h minimize the latter term of the ab ov e inequality (4 2). Then 8 β L 2 h 2 β − 1 − 2 Γ h 3 n = 0 from whic h we infer that h =  Γ 4 β L 2  1 2 β + 2 n − 1 2 β + 2 . (43) Substituting (43) to (42) leads to g ( w ∗ h ) ≤ 2 2 β + 6 2 β + 2 Γ 2 β 2 β + 2 L 4 2 β + 2 β 2 β 2 β + 2 n − 2 β 2 β + 2 . Therefore (16) implies (19). 14 Qiyu JIN et al. 6.2 Pro of of Theorem 2 W e shall use following lemma to finish the Pro of o f Theo rem 2. The lemma can be deduced form the results in Borovk ov [6], see also Merlev ede, P eligra d and Rio [31] Lemma 1 If, for some δ > 0 , γ ∈ (0 , 1 ) and K > 1 we have sup E ex p ( δ | X i | γ ) ≤ K, i = 1 , ..., n, then ther e ar e two p ositive c onstants c 1 and c 2 dep ending only on δ, γ and K such tha t, for any t > 0 , P n X i =1 X i ≥ t ! ≤ exp  − c 1 t 2 /n  + n exp ( − c 2 t γ ) . Pro of of Theorem 2 . Recall that M and m is giv en b y (5) and (6) resp ectively . Consider that | b ρ 2 x 0 ( x ) − ρ 2 f ,x 0 ( x ) | ≤     1 m S ( x ) + 1 m R ( x )     , (44) where S ( x ) = X y ∈ U x 0 ,η Z ( y ) , Z ( y ) = ( Y ( y ) − Y ( T x y )) 2 − ( f ( y ) − f ( T x y )) 2 − f ( y ) − f ( T x y ) , (4 5) and R ( x ) = X y ∈ U x 0 ,η  ( f ( y ) − f ( T x y )) 2 + f ( y ) + f ( T x y )  − ( f ( x 0 ) − f ( x )) 2 − 2 f ( x 0 ) . Since Y ( x ) has the Poisson distribution, with mea n f ( x ) and v ar iance f ( x ), E  e Y ( x )  = + ∞ X k =0 e k f k ( x ) e − f ( x ) k ! = ef ( x ) e ( e − 1) f ( x ) ≤ eΓ e ( e − 1) Γ . (46) F r om the inequalit y (46), we ea sily deduce sup E  e | Z ( y ) | 1 / 2  ≤ sup E  e Y ( y )+ Y ( T x y )+2 Γ +2 √ Γ  ≤ ( eΓ ) 2 e 2 eΓ +2 √ Γ (47) By Lemma 1, we s e e that there are tw o positive constants c 5 and c 6 such that for any z > 0 , P  1 m | S ( x ) | ≥ z √ m  ≤ exp( − c 5 z 2 ) + m e x p( − c 6 ( √ mz ) 1 2 ) . (48) A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 15 Considering m = (2 N η + 1) 2 and η = c 1 n − α ( (1 − β ) + ) 2 β +2 < α < 1 2 , w e ha ve m = c ′ 1 n 1 − 2 α (1 + o (1)). Ther efore, substituting z = q 1 c 5 ln n 2 int o the inequality (48), w e see that for n large enough, P   1 m | S ( x ) | ≥ q 1 c 5 ln n 2 √ m   ≤ 2 ex p  − ln n 2  = 2 n 2 . F r om this inequalit y we easily deduce that P   max x ∈ U x 0 ,h 1 m | S ( x ) | ≥ q 1 c 5 ln n 2 √ m   ≤ X x ∈ U x 0 ,h P   1 m | S ( x ) | ≥ q 1 c 5 ln n 2 √ m   ≤ 2 n . W e arrive at P ( B ) ≤ c 7 n − 1 , (49) where B = { max x ∈ U x 0 ,h 1 m | S ( x ) | < c 8 n α − 1 2 √ ln n } and c 8 is a constant de- pending only o n β and L . It is easy to see that R ( x ) = O  n α − 1 2  . (50) In the set B , the inequality (50) implies that max x ∈ U x 0 ,h | b ρ 2 x 0 ( x ) − ρ 2 f ,x 0 ( x ) | ≤ c 8 n − β 2 β + 2 √ ln n + O  n α − 1 2  = O ( n − β 2 β + 2 √ ln n ) . (51) Combining (49) and (5 1), we obtain (22). 6.3 Pro of of Theorem 3 T aking in to account (25), (27), and the indep endence of ǫ ( x ), we have E {| b f ′ h ( x 0 ) − f ( x 0 ) | 2   Y ( x ) , x ∈ I ′′ x 0 } ≤ g ′ ( b w h ) , (52) where g ′ ( w ) =   X x ∈ U ′ x 0 ,h w ( x ) ρ f ,x 0 ( x )   2 + f ′ ( x 0 ) X x ∈ I ′ x 0 w 2 ( x ) . By the pro of of Theorem 1, w e obtain g ′ ( w ∗ h ) ≤ 3 2 2 2 β + 6 2 β + 2 Γ 2 β 2 β + 2 L 4 2 β + 2 β 2 β 2 β + 2 n − 2 β 2 β + 2 ! . (53) By Theorem 2 and its pro of, for b ρ ′ x 0 is defined by (25), there is a cons tan t c 2 such that P ( max x ∈ U ′ x 0 ,h    b ρ ′ 2 x 0 ( x ) − ρ 2 f ,x 0 ( x )    ≥ c 2 n α − 1 2 √ ln n ) = O  n − 1  . (54) 16 Qiyu JIN et al. Let B = n max x ∈ U ′ x 0 ,h    b ρ ′ 2 x 0 ( x ) − ρ 2 f ,x 0 ( x )    ≤ c 2 n α − 1 2 ln n o . On the s et B , w e hav e ρ 2 f ,x 0 ( x ) − c 2 n α − 1 2 √ ln n < b ρ ′ 2 x 0 ( x ) < ρ 2 f ,x 0 ( x ) + c 2 n α − 1 2 √ ln n , from which we infer that b w ( x ) = e − b ρ ′ 2 x 0 ( x ) H 2 ( x 0 ) P x ′ ∈ U ′ x 0 ,h e − b ρ ′ 2 x 0 ( x ′ ) H 2 ( x 0 ) ≤ e − ρ 2 f,x 0 ( x ) − c 2 n α − 1 2 √ ln n H 2 ( x 0 ) P x ′ ∈ U ′ x 0 ,h e − ρ 2 f,x 0 ( x ′ )+ c 2 n α − 1 2 √ ln n H 2 ( x 0 ) ≤ e − ρ 2 f,x 0 ( x ) H 2 ( x 0 )  1 + 2 c 2 n α − 1 2 √ ln n H 2 ( x 0 )  P x ′ ∈ U ′ x 0 ,h e − ρ 2 f,x 0 ( x ′ ) H 2 ( x 0 )  1 − c 2 n α − 1 2 √ ln n H 2 ( x 0 )  =   1 + 2 c 2 n α − 1 2 √ ln n H 2 ( x 0 ) 1 − c 2 n α − 1 2 √ ln n H 2 ( x 0 )   e − ρ 2 f,x 0 ( x ) H 2 ( x 0 ) P x ′ ∈ U ′ x 0 ,h e − ρ 2 f,x 0 ( x ′ ) H 2 ( x 0 ) =   1 + 2 c 2 n α − 1 2 √ ln n H 2 ( x 0 ) 1 − c 2 n α − 1 2 √ ln n H 2 ( x 0 )   w ∗ h ( x ) . This implies that g ′ ( b w h ) ≤   1 + 2 c 2 n α − 1 2 √ ln n H 2 ( x 0 ) 1 − c 2 n α − 1 2 √ ln n H 2 ( x 0 )   2 g ′ ( w ∗ h ) . The condition 4 c 2 n α − 1 2 √ ln n < γ ≤ H 2 ( x 0 ) implies that   1 + 2 c 2 n α − 1 2 √ ln n H 2 ( x 0 ) 1 − c 2 n α − 1 2 √ ln n H 2 ( x 0 )   2 ≤ 2 . Consequently , (52) b e c omes E  | b f ′ h ( x 0 ) − f ( x 0 ) | 2   Y ( x ) , x ∈ I ′′ x 0 , B  ≤ 2 g ′ ( w ∗ h ) . (5 5) Since the function f satisfies the H¨ o lder condition, E  | b f ′ h ( x 0 ) − f ( x 0 ) | 2   Y ( x ) , x ∈ I ′′ x 0  < g ′ ( b w h ) ≤ c 9 , (56) A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 17 for a constan t c 9 > 0 depending only on β and L . Combining (22), (55), and (56), w e hav e E  | b f ′ h ( x 0 ) − f ( x 0 ) | 2   Y ( x ) , x ∈ I ′′ x 0  = E  | b f ′ h ( x 0 ) − f ( x 0 ) | 2   Y ( x ) , x ∈ I ′′ x 0 , B  P ( B ) + E  | b f ′ h ( x 0 ) − f ( x 0 ) | 2   Y ( x ) , x ∈ I ′′ x 0 , B  P ( B ) ≤ 2 g ′ ( w ∗ h ) + O  n − 1  Now, the a s sertion of the theorem is obtained easily if we take into account (53). References 1. I. Abraham, R. Abraham, A. Desolneux , and S. Li-Thiao -T e. Significant edges in the case o f non-stationar y gaus s ian noise. Pattern r e c o gnition , 40(11):32 77–329 1, 2007. 2. M. Aharon, M. Elad, and A. Bruckstein. r mk -svd: An a lgorithm for de- signing overcomplete dictionaries for sparse r epresentation. IEEE T r ans. Signal Pr o c ess. , 54(11):4 3 11–432 2, 2006. 3. F. Alter, Y. Matsushita, and X. T ang. An in tens it y similarity measure in low-ligh t conditions . Computer Vision–E CCV 2006 , pages 267–2 80, 2 006. 4. F.J. Anscombe. The tr ansformation of p oisson, binomial and neg ative- binomial data. Biometrika , 35(3/4 ):246–25 4, 194 8. 5. A. Beck and M. T eb oulle. F ast gradient-based algor ithms for constrained total v aria tion imag e denoising and deblurring problems. IEEE T r ans. Image Pr o c ess. , 18(1 1):2419–2 434, 200 9. 6. A.A. Borovko v. E stimates for the distributio n of sums a nd maxima of sums of ra ndo m v ariables without the cramer c o ndition. Sib erian Mathematic al Journal , 41(5):811– 848, 200 0. 7. J. Boulanger , J.B. Sibarita, C. Kervrann, and P . Bouthemy . Non- parametric r egression for patch-based fluores cence micros c op y image se- quence denoising. In in Pr o c. of IEEE Int. Symp. on Biome dic al Imaging: F r om Nano to M acr o, ISBI2008 , pages 748– 7 51. IEE E, 2 008. 8. A. Buades, B. Coll, and J.M. Morel. A review o f imag e denoising al- gorithms, with a new one. SIAM Journal on Mult isc ale Mo deling and Simulation , 4(2):490 – 530, 200 5. 9. T. Buades, Y. Lou, JM Morel, and Z. T ang. A note on m ulti-image de- noising. In Int. workshop on L o c al and Non-Lo c al Appr oximation in Image Pr o c essing , pag es 1 – 15, August 2009. 10. C.A. Deledalle, F. T upin, and L. Denis. Poisson nl means: Unsup ervised non loc a l means for poiss on noise. In IEEE Int. Conf. on Image Pr o c ess. (ICIP), 2010 17th , pages 801–8 0 4. IEEE, 2010. 11. D.L. Dono ho and J.M. Jo hnstone. Ideal spa tial ada ptation b y wav ele t shrink a g e. Biometrika , 81(3 ):4 25, 19 94. 18 Qiyu JIN et al. 12. J. F an. Lo cal linear r egression smo other s and their minimax efficiencies. The Annals of Statistics , pages 196–21 6, 199 3 . 13. J.Q . F an and I. Gijbels. Lo cal po ly nomial mo delling and its applications. In Cha pman & Hall, L ondon , 19 96. 14. P . F ryzlewicz, V. Delouille, a nd G.P . Naso n. Goes - 8 x-ray s e ns or v ar iance stabilization using the multiscale data - driven haa r–fisz transfor m. J. R oy. Statist. So c. ser. C , 56(1):9 9 –116, 20 07. 15. P . F ryzlew ic z and G.P . Nason. A haar-fisz alg orithm for pois son intensit y estimation. J. C omp. Gr aph. Stat. , 13(3):621–63 8, 2 004. 16. D.K. Hammond and E.P . Simoncelli. Imag e mo deling and denoising with orientation-adapted gaussian scale mixtures. IEEE T r ans. Image Pr o c ess. , 17(11):20 89–210 1, 2008. 17. K . Hirak aw a and T.W. Parks. Image denoising using total least s q uares. IEEE T r ans. Image Pr o c ess. , 15 (9 ):2730–27 42, 20 06. 18. M. J a nsen. Mult iscale p oisson data s moothing . J. R oy. S t atist. So c. B , 68(1):27– 48, 2006 . 19. Q iy u Jin, Ion Grama, a nd Q ua nsheng Liu. A new po isson noise filter based on w eights o ptimization. arXiv pr eprint arXiv:12 01.596 8 , 2012 . 20. Q .Y. Jin, I. Gra ma, a nd Q.S. Liu. Removing ga ussian no ise b y optimiza- tion of w eights in non-lo cal means. http://arxiv.or g/abs/110 9.5640 . 21. V. Ka tk ovnik, A. F oi, K . Egia z arian, and J. Astola. F rom lo ca l k e rnel to nonlo cal multiple-model image denoising. Int. J. Comput. V is. , 8 6(1):1– 32, 2010. 22. C. Kervrann a nd J. Boulang er. Optimal spatia l adaptation for patch- based ima ge deno ising. IEEE T r ans. Image Pr o c ess. , 15(10 ):2866–2 8 78, 2006. 23. E .D. K olaczyk and R.D. Now ak. Multiscale generalised linear mo dels for nonparametric function estimation. Biometrika , 92(1 ):119, 2005. 24. T. Le, R. Chartrand, and T. J. Asaki. A v ar iational approach to recon- structing images corrupted b y p oisson noise. Journal of Mathematic al Imaging and Vision , 27(3):257 –263, 2007 . 25. S. Lefkimmiatis, P . Maragos, a nd G. Papandreou. Bayesian inference on m ultiscale mo dels for po isson intensit y e s timation: Applications to photon- limited ima g e denoising. IEEE T r ans. Image Pr o c ess. , 1 8(8):1724– 1741, 2009. 26. F. Luisier, C. V onesch, T. Blu, and M. Unser. F a st interscale wa velet denoising of p oisson-cor rupted ima ges. Signal Pr o c ess. , 90(2):415 –427, 2010. 27. J. Mair a l, G. Sa pir o, and M. Elad. Le arning m ultiscale sparse represen- tations for image and video r e storation. SIAM Multisc ale Mo deling and Simulation , 7(1):214 – 241, 200 8. 28. M. Makitalo and A. F oi. On the inversion o f the ansco m be trans fo rmation in low-count p oisson ima g e deno ising. In Pr o c. Int. Workshop on L o c al and Non-L o c al App r ox. in Image Pr o c ess., LNLA 20 09, T uusula, Finla nd , pages 26–32. IEEE, 2009. A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 19 29. M. Makitalo and A. F o i. Optimal in version of the a nscombe transforma - tion in low-coun t p oisson image denoising. IEEE T r ans. Image Pr o c ess. , 20(1):99– 109, 201 1. 30. J. Mandel. Use of the singular v alue decomp osition in regress io n analysis . The Americ an Statistician , 36(1):15–24 , 1982. 31. F. Merlev` ede, M. P eligrad, and E. Rio. A bernstein type inequality a nd mo derate deviations for w eakly dep endent seque nc e s . Pr ob ab. The ory R e- late d Field s , 2 010. 32. I. Mo on and B. Javidi. Three dimensional imaging and reco gnition using truncated photon counting mo del and par a metric maximum lik eliho o d estimator. Optics ex pr ess , 17(18 ):15709– 1 5715, 2009 . 33. R.D. Nowak and E.D. Kolaczy k. A statistical m ultiscale framework for po isson inv erse pro blems. IE EE T r ans. Info. The ory , 46(5):1811 –1825, 2000. 34. J. Polzehl and V. Sp okoin y . Propag ation-separa tion approach for lo cal likelihoo d estimation. Pr ob ab. The ory R el. , 135(3):335– 362, 2 006. 35. J. Portilla, V. Strela , M.J. W ain wright, and E.P . Simoncelli. Image de- noising using scale mixtur e s of ga ussians in the wa v elet domain. IEEE T r ans. Image Pr o c ess. , 12(11):1338– 1351, 20 03. 36. S. Roth and M.J. Black. Fields of exp erts. Int. J. Comput. Vision , 82(2):205 –229, 2009 . 37. G. R. T errell and D. W. Scott. V ariable kernel densit y estimation. Th e Annals of Statistics , pag es 1 236–12 65, 1992 . 38. G.M.P . v an K empen, H.T.M. v an der V o ort, J.G.J . Bauman, and K.C. Strasters. Compar ing maximum lik eliho od estimation and constr a ined tikhonov-miller restora tion. IEEE Engine ering in Me dicine and Biolo gy Magazine , 15(1):76– 8 3, 199 6. 39. B. Zhang , J.M. F a dili, and J.L. Starc k. W av e lets, ridgelets, a nd c urvelets for poisso n noise remov al. IEEE T r ans. Image Pr o c ess. , 17(7):10 9 3–1108 , 2008. 20 Qiyu JIN et al. (a) O riginal image (b) Noisy image (c) N LM PF N M I S E = 0 . 0260 P S N R = 31 . 45 (d)O WPNF (e) Poisson NLM (f ) EUI+BM3D N M I S E = 0 . 0259 N M I S E = 0 . 0790 N M I S E = 0 . 0358 P S N R = 31 . 31 P S N R = 31 . 12 P S N R = 31 . 95 (g) M S-VST + 7 / 9 (h) MS-VST + B3 N M I S E = 0 . 0602 N M I S E = 0 . 0810 P S N R = 31 . 64 P S N R = 31 . 35 Fig. 2 D enoising an image of simulated s p ots of different radii (image size: 256 × 256). (a) simulated sources (amplitudes ∈ [0 . 08 , 4 . 99]; bac kground = 0 . 03); (b) observ ed coun ts; (c) NLMPF ( M = 19 × 19, m = 1 3 × 13, d = 3, σ H = 2 . 5, µ = 1 N M I S E = 0 . 0260) ; (d) Optimal W eigh ts Filter ( M = 19 × 19, m = 13 × 13, d = 2 and H = 1 , N M I S E = 0 . 0259); (e) Po isson NLM ( N M I S E = 0 . 0790); (f ) Exact un biased inv erse + BM3D ( N M I S E = 0 . 0358) ; (g) MS-V ST + 7 / 9 biorthogonal wa v elet ( J = 5, F P R = 0 . 01, N max = 5 iterations, N M I S E = 0 . 0602); (h) MS-VST + B3 isotropic wa velet ( J = 5, F P R = 0 . 01, N max = 5 iterations, N M I S E = 0 . 0810). A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 21 (a) O riginal image (b) Noisy image (c) N LM PF N M I S E = 0 . 0284 P S N R = 28 . 09 (d)O WPNF (e) Poisson NLM (f ) EUI+BM3D N M I S E = 0 . 0285 N M I S E = 0 . 0346 N M I S E = 0 . 0297 P S N R = 27 . 80 P S N R = 27 . 77 P S N R = 28 . 04 (g) M S-VST + 7 / 9 (h) MS-VST + B3 N M I S E = 0 . 0357 N M I S E = 0 . 0338 P S N R = 27 . 57 P S N R = 27 . 94 Fig. 3 D enoising a galaxy image (image size: 256 × 256). (a) galaxy image (inten sity ∈ [0 , 5]); (b) observ ed count s; (c) NLMPF ( M = 13 × 13, m = 3 × 3, d = 2, σ H = 1, µ = 0 . 6 N M I S E = 0 . 0284) ; (d) Optimal W eigh ts Filter ( M = 15 × 15, m = 5 × 5, d = 2 and H = 1, N M I S E = 0 . 0285); (e) Poisson NLM ( N M I S E = 0 . 0346); (f ) Exact unbiased inv erse + BM3D ( N M I S E = 0 . 0297 ) ; (g) MS-VST + 7 / 9 biorthogonal w av elet ( J = 5, F P R = 0 . 0001, N max = 5 iterations, N M I S E = 0 . 0357); (h) MS-VST + B3 isotropic wa velet ( J = 3, F P R = 0 . 0001, N max = 10 iterations, N M I S E = 0 . 0338). 22 Qiyu JIN et al. (a) O riginal image (b) Noisy image (c) N LM PF N M I S E = 0 . 0140 P S N R = 24 . 69 (d)O WPNF (e) Poisson NLM (f ) EUI+BM3D N M I S E = 0 . 0 . 0162 N M I S E = 0 . 0154 N M I S E = 0 . 0 . 0121 P S N R = 23 . 90 P S N R = 24 . 94 P S N R = 25 . 89 (g) M S-VST + 7 / 9 (h) MS-VST + B3 N M I S E = 0 . 0 . 0193 N M I S E = 0 . 0416 P S N R = 24 . 49 P S N R = 24 . 03 Fig. 4 P oisson denoising of sm ooth ri dges (image size: 256 × 256). (a) int ensity im age (the peak int ensities of the 9 v ertical ridges v ary progressively from 0 . 1 to 0 . 5; the i nclined ri dge has a maximum inte nsity of 0 . 3; backg round = 0 . 05); (b) Poisson noisy image; (c) NLM PF ( M = 9 × 9, m = 21 × 21, d = 4, σ H = 0 . 5, µ = 0 . 4, N M I S E = 0 . 0140) ; (d) Optimal W eigh ts Filter ( M = 9 × 9, m = 19 × 19, d = 3 and H = 2, N M I S E = 0 . 0162); (e) Poisson NLM ( N M I S E = 0 . 0154); (f ) Exact unbiased inv erse + BM 3D ( N M I S E = 0 . 0121); (g) MS-VST + 7 / 9 biorthogonal wa velet ( J = 5, F P R = 0 . 001, N max = 5 iterations, N M I S E = 0 . 0193); (h) MS-VST + B3 i sotropic wa velet ( J = 3, F P R = 0 . 00001, N max = 10 i terations, N M I S E = 0 . 0416). A Non-Lo cal Means Fi lter for Remov ing the Poisson Noise 23 (a) O riginal image (b) Noisy image (c) N LM PF N M I S E = 0 . 1150 P S N R = 24 . 71 (d)O WPNF (e) Poisson NLM (f ) EUI+BM3D N M I S E = 0 . 1061 N M I S E = 0 . 1207 N M I S E = 0 . 0683 P S N R = 24 . 60 P S N R = 24 . 72 P S N R = 25 . 92 (g) M S-VST + 7 / 9 (h) MS-VST + B3 N M I S E = 0 . 2391 N M I S E = 0 . 3777 P S N R = 21 . 81 P S N R = 20 . 22 Fig. 5 Poisson denoising of the Barbara image (image s i ze: 256 × 256). (a) intensit y im age (in tensit y ∈ [0 . 93 , 15 . 73]); (b) Poisson noisy image; (c) NLM PF ( M = 15 × 15, m = 21 × 21, d = 0, µ = 1, N M I S E = 0 . 1150); (d) Optimal W eights Filter ( M = 15 × 15, m = 21 × 21 and d = 0, N M I S E = 0 . 1061); (e) P oisson NLM ( N M I S E = 0 . 1207); (f ) Exac t un biased inv erse + BM3D ( N M I S E = 0 . 0863) (h) MS-VST + 7 / 9 biorthogonal w a ve let ( J = 4, F P R = 0 . 001, N max = 5 iterations, N M I S E = 0 . 2391); (h) MS-VST + B3 isotropic wa velet ( J = 5, F P R = 0 . 001, N max = 5 i terations, N M I S E = 0 . 3777). 24 Qiyu JIN et al. (a) O riginal image (b) Noisy image (c) N LM PF N M I S E = 0 . 0785 P S N R = 29 . 08 (d)O WPNF (e) Poisson NLM (f ) EUI+BM3D N M I S E = 0 . 0794 N M I S E = 0 . 0816 N M I S E = 0 . 0643 P S N R = 29 . 91 P S N R = 29 . 40 P S N R = 30 . 18 (g) M S-VST + 7 / 9 (h) MS-VST + B3 N M I S E = 0 . 0909 N M I S E = 0 . 1487 P S N R = 28 . 87 P S N R = 26 . 69 Fig. 6 Poisson denoising of fluorescent tubules (image si ze: 256 × 256). (a) int ensity image (in tensit y ∈ [ 0 . 53 , 16 . 93]); (b) P oisson noisy image; (c) NLMPF ( M = 7 × 7, m = 13 × 13, d = 2, σ H = 2, µ = 1, N M I S E = 0 . 0785); ( d) Optimal W eights Fil ter ( M = 11 × 11, m = 17 × 17, d = 1 and H = 0 . 6, N M I S E = 0 . 0794); (e) Poisson NLM ( N M I S E = 0 . 0816); (f ) Exact unbiase d inv erse + BM3D ( N M I S E = 0 . 0643) (g) MS-VST + 7 / 9 biorthogonal wa velet ( J = 5, F P R = 0 . 0001, N max = 5 iterations, N M I S E = 0 . 0909); (h) MS-VST + B3 i sotropic wa velet ( J = 5, F P R = 0 . 001, N max = 10 iterations, N M I S E = 0 . 1487).