Deconvolution of confocal microscopy images using proximal iteration and sparse representations

DECONV OLUTION OF CONFOCAL MICR OSCOPY IMA GES US ING PRO XIMAL ITERA TION AND SP AR SE REPRESENT A TIONS F .-X. Dup ´ e a , M.J. F adili a and J .-L. Star ck b a GREYC UMR CNRS 6072 b D APNIA/SEDI-SAP CEA-Saclay 14050 Caen France 91191 Gif-sur-Yv ette France ABSTRA CT W e propo se a deconv olution algorith m for images blurred and degraded by a Poisson noise. The algorith m uses a f ast proxim al backward-fo rward splitting iteration. This it er- ation minimizes an energy wh ich combines a non -linear data fidelity term, adap ted to Poisson noise, and a non - smooth sparsity-pro moting regularization ( e.g ℓ 1 -norm ) over the image representatio n co efficients in some dictio - nary of transfo rms ( e.g. wav elets, cu rvelets). Our results on simulated m icroscopy images o f neuro ns an d cells are confro nted to some state-of- the-art algorith ms. They show that our approach is very com petitiv e, an d as expected, the imp ortance of the non-lin earity du e to Poisson noise is mo re salient at low an d med ium in tensities. Finally an experiment on real fluorescen t conf ocal m icroscopy data is reported . Index T erms — Deco n volution, Poisson noise, Confo- cal m icroscopy , Itera ti ve thresh olding, Sparse representa- tions. 1. INTRODUCTION Fluorescent m icroscopy suffers fro m two main sources o f degradation: th e optica l system and th e acquisition noise. The optical system has a finite resolution introducin g a blur in the observation. This degradation is modeled as conv olution with the point spread functio n (PSF). The sec- ond source o f image degradation is due to Poisson coun t process ( shot noise). In p resence o f Poisson no ise, sev- eral deconv olu tion algorithms hav e b een prop osed such as the well-k nown Richardson-Lucy (RL) algorithm or T ikho nov-Miller in verse filter, to name a few . RL is ex- tensiv ely used for its good adaptation to Poisson noise, but it tends to amplify noise after a few iterations. Reg- ularization can be in troduc ed in order to avoid this issue. In b iological imag ing deconv olutio n, many kind s of reg- ularization have been su ggested: to tal variation with RL [1] which gives staircase artifacts, T ikho nov with RL (see [2] for a revie w), etc. W av elets h ave also been used as a regularization scheme when deconvolving biome dical im- ages; [3] presents a version of RL com bined with wa velets denoising , and [4] use the th resholded Landweb er itera- tion of [5]. The latter approach implicitly assumes that the contamina ting noise is Gau ssian. In the context o f de conv olution with Gaussian white noise, sp arsity-pro moting regularization over o rthog onal wa velet co efficients ha s been recently propo sed [5, 6]. Generalization to frames was prop osed in [7, 8]. In [9], the auth ors presented an image deconv olutio n algorithm by iterative th resholdin g in an overcomplete d ictionary of transforms. Howe ver , all sparsity-b ased approaches published so far ha ve mainly focused on Gaussian noise. In this paper, we propose an image dec onv olution al- gorithm for d ata blurred and co ntaminated by Poisson noise. Th e Poisson n oise is han dled prop erly by using th e Anscombe variance stabilizing tr ansform (VST), lead ing to a non-lin ear degrad ation equation with additive Gaus- sian noise, see (2). The deconv olution pro blem is then formu lated as the minim ization of a con vex functional with a no n-linear data-fidelity ter m reflecting the noise proper ties, and a n on-smoo th sparsity-prom oting penalty over representation coefficients of the ima ge to r estore, e.g. wav elet or cur velet coefficients. Inspired b y the work in [6], a fast pro ximal iterative a lgorithm is p ropo sed to solve the minimization prob lem. Exper imental results ar e carried ou t to compare our app roach on a set of simu- lated an d real confocal m icroscopy images, and show the striking benefits gained from taking into account the Pois- son nature o f the noise and th e morph ological structure s in volved in the imag e. Notation Let H a real Hilbert spac e, here a finite dimensional vector subspace of R n . W e denote by k . k 2 the no rm associated with the inne r product in H , and I is the ident ity operator on H . x and α are respecti vel y reordered vec tors of image samples and transform coeffic ients. A func tion f is co- erci ve, if lim k x k 2 → + ∞ f ( x ) = + ∞ . Γ 0 ( H ) is the class of all proper lo wer semi-continuous con vex functions from H to ] − ∞ , + ∞ ] . 2. PROBLEM ST A TEMENT Consider the image for mation model whe re an input image x is blur red b y a poin t spread function (PSF) h an d con- taminated by Poisson n oise. The observed image is the n a discrete collection of counts y = ( y i ) 1 6 i 6 n where n is the number of pixels. Each coun t y i is a realization o f an in- depend ent Poisson rando m variable with a mean ( h ⊛ x ) i , where ⊛ is the circular conv olution o perator . Formally , this writes y i ∼ P (( h ⊛ x ) i ) . A naive solution to this deconv olutio n problem would be to ap ply trad itional approach es design ed for Gaussian noise. But this would be awkward as (i) the n oise tend s to G aussian on ly for large me an ( h ⊛ x ) i (central limit theorem) , and (ii) the noise variance depen ds on the mean anyway . A more adapted way would be to ado pt a b ayesian framework with an appropr iate a nti-log-likeliho od sco re reflecting th e Poisson statistics of the n oise. Unfor tunately , doing so, we would end-up with a fu nctional which do es not satisfy s om e key prop erties (th e Lipschitzian p rop- erty stated a fter (3)), hen ce preventing us fr om using th e backward-fo rward splitting pr oximal algo rithm to solve the op timization p roblem . T o circumvent this difficulty , we propo se to handle the noise statistical p roperties by using the Anscombe VST defined as z i = 2 q y i + 3 8 , 1 6 i 6 n. (1) Some previous authors [10] have alread y suggested to use the Anscombe VST , an d then deconv olve with wa velet- domain regula rization as if the stabilized ob servation z i were linearly degraded by h and contaminated by additi ve Gaussian no ise. But this is no t valid as stan dard asym p- totic results of the Anscombe VST state that z i = 2 q ( h ⊛ x ) i + 3 8 + ε, ε ∼ N (0 , 1) (2) where ε is an additive white Gaussian noise of u nit vari- ance 1 . In words, z is non-lin early related to x . In Sec- tion 4, we p rovide a n elegant optimization pro blem and a fixed p oint algorithm tak ing into accoun t such a non- linearity . 3. SP ARSE IMA GE REPRESENT A TION Let x ∈ H be an √ n × √ n image . x can be wr itten as the sup erposition of elem entary atoms ϕ γ parametrize d by γ ∈ I such that x = P γ ∈ I α γ ϕ γ = Φ α, |I | = L, L > n . W e deno te by Φ the dictionar y i.e. the n × L matrix whose co lumns are the gen erating wa veforms ( ϕ γ ) γ ∈ I all n ormalized to a unit ℓ 2 -norm . The fo rward transform is then defin ed by a non-n ecessarily square ma- trix T = Φ T ∈ R L × n . In the rest of th e paper, Φ will be an orthob asis or a tig ht frame with constant A . 4. SP ARSE ITE RA TIVE DECONV OLUTION 4.1. Optimizat ion pr oblem The class of min imization prob lems we are inter ested in can be stated in the general form [6]: arg min x ∈H f 1 ( x ) + f 2 ( x ) . (3) where f 1 ∈ Γ 0 ( H ) , f 2 ∈ Γ 0 ( H ) and f 1 is differentiable with κ -Lipschitz gr adient. W e denote by M the set o f solutions. From (2), we immediately deduce the data fidelity term F ◦ H ◦ Φ ( α ) , with (4) F : η 7→ n X i =1 f ( η i ) , f ( η i ) = 1 2  z i − 2 q η i + 3 8  2 , 1 Rigorously speaking, the equation is to be understood in an asymp- totic sense. where H d enotes the conv olution operato r . From a statisti- cal p erspective, (4) corresp onds to the anti-log-likelihood score. Adopting a bayesian framework an d using a standard maximum a posteriori ( MAP) ru le, our go al is to minimize the following functional with respect to the representation coefficients α : ( P λ,ψ ) : arg min α J ( α ) (5) J : α 7→ F ◦ H ◦ Φ ( α ) | {z } f 1 ( α ) + ı C ◦ Φ ( α ) + λ L X i =1 ψ ( α i ) | {z } f 2 ( α ) , where we implicitly assumed that ( α i ) 1 6 i 6 L are indepen- dent an d ide ntically distributed. The penalty fu nction ψ is chosen to enforc e sp arsity , λ > 0 is a regularizatio n parameter an d ı C is the indicator fu nction of a conve x set C . In o ur case, C is the p ositiv e orthan t. W e r emind that the positivity constraint is becau se we are fitting Poisson intensities, which are positiv e by nature. 4.2. P roximal itera tion W e now presen t o ur m ain proximal iter ativ e algorithm to solve the minimization problem ( P λ,ψ ) : Theorem 1. ( P λ,ψ ) h as at least one solutio n ( M 6 = ∅ ). The solution is uniq ue if ψ is strictly conve x or if Φ is a orthobasis a nd Ker( H ) = ∅ . F or t > 0 , let ( µ t ) t be such that 0 < inf t µ t 6 sup t µ t <  3 2  3 / 2 /  2 A k H k 2 2 k z k ∞  . F ix α 0 ∈ C ◦ Φ , for every t > 0 , set α t +1 = prox µ t f 2 ( α t − µ t ∇ f 1 ( α t )) , (6) wher e ∇ f 1 is the gr ad ient of f 1 and pro x µ t f 2 is computed using the following iteration: let P t ν t (1 − ν t ) = + ∞ , take γ 0 ∈ H , and defin e the sequence of iter ates: γ t +1 = γ t + ν t  rprox µ t λ Ψ+ 1 2 k . − α k 2 ◦ rprox ı C ′ − I  ( γ t ) , (7) wher e prox µ t λ Ψ+ 1 2 k . − α k 2 ( γ t ) =  prox µ t λ 2 ψ (( α i + γ t i ) / 2)  1 6 i 6 L , P C ′ = prox ı C ′ = A − 1 Φ T ◦ P C ◦ Φ +  I − A − 1 Φ T ◦ Φ  , rprox ϕ = 2 pr ox ϕ − I and P C is the pr ojector on to the positive orthant ( P C η ) i = max( η i , 0) . Th en, γ t ⇀ γ and prox µ t f 2 ( α ) = P C ′ ( γ ) . (8) Then ( α t ) t > 0 conver ges (weakly) to a solution of ( P λ,ψ ) . A proof can be found in [11]. prox δψ is giv en by , Theorem 2. Suppose that (i) ψ is conve x even-symmetric , non-n e gative and non-d ecr easing on [0 , + ∞ ) , a nd ψ (0) = 0 . (ii) ψ is twice differ entiab le on R \ { 0 } . (iii) ψ is c on- tinuous on R , it is n ot necessarily smooth at zer o a nd ad- mits a p ositive right d erivative at zer o ψ ′ + (0) > 0 . The n, the pr oximity operator prox δψ ( β ) = ˆ α ( β ) h as exactly one continuo us solutio n decoupled in each coor dinate β i : ˆ α i ( β i ) = ( 0 if | β i | ≤ δ ψ ′ + (0) β i − δ ψ ′ ( ˆ α i ) if | β i | > δ ψ ′ + (0) (9) See [ 9]. Among the most po pular pen alty fu nctions ψ satisfy ing the above req uiremen ts, we have ψ ( α i ) = | α i | , in which case the associated prox imity operato r is soft-thresho lding. Th erefor e, (6) is e ssentially an iterative thresholdin g algorithm with a positivity constraint. 5. RESUL TS The performanc e of our approach has been assessed on se veral datasets of b iological images: a neu ron phanto m and a cell. Our algor ithm was compar ed to RL with total variation regularization (RL-TV [1]), RL with mu lti- resolution supp ort wa velet regularizatio n (RL- MRS [12]), the naive proxim al metho d that would assum e the noise to be Gaussian (NaiveGauss [4]), and the ap proach of [10] (AnsGauss). For all results pr esented, each algorithm was run with 200 it era tions, enough to reach con vergence. Simulations were carried out with an approximated but re- alistic PSF [13] whose parameters are obtained from a real confoc al micro scope. As usual, the choice o f λ is cru cial to balance b etween regularization and deco n volution. For all the situation s below , λ was adjusted in order to reach the best compr omise. In Fig. 1(a), a ph antom of a neuron with a m ushroo m- shaped spine is depicted . The maximum intensity is 30. Its blurred and blur red+no isy versions are in (b) and (c). W ith this n euron , a nd for NaiveGauss, AnsGau ss and our ap- proach , the diction ary Φ contained the wa velet orthog onal basis. The de conv olution r esults are shown in Fig. 1(d)- (h). As expe cted the worst results a re fo r the AnsGau ss version, as it does n ot take ca re o f the non-lin earity of the Anscombe VST . RL-TV shows rather good results b ut the backgr ound is full of artifacts. Our appr oach provid es a vi- sually ple asant deco n volution result on this example. It ef- ficiently restores the spine, altho ugh t he backg round is not fully cleaned . RL-MRS also exhibits goo d d econv olution results. Qualitative visual resu lts are c onfirmed by quan- titati ve m easures of the qu ality of decon volution, where we used b oth the ℓ 1 -error (adapted to Poisson noise), a nd the traditional MSE criteria. Th e ℓ 1 -errors fo r this imag e are shown by T ab. 1 (similar results were obtain ed for the MSE). In genera l, our appr oach pe rforms we ll. It is com- petitiv e compared to RL-MRS which is designed to d i- rectly handle Poisson n oise. At low intensity r egimes, our approa ch an d RL-MRS are the two algo rithms that giv e th e best results. At h igh intensity , RL-TV per forms very well, although RL-MRS, Nai veGauss an d our appr oach are very close to it. NaiveGauss perform s po orly at lo w i nte nsity as it does not correspo nd to a d egradation mo del with Pois- son noise. AnsGauss giv es the worst results probably be- cause it do es not hand le the non-linear ity o f the degrada- tion m odel (2) after the VST . T o assess the computatio nal burden o f the compared alg orithms, T ab . 2 summarizes the execution times with an Intel PC Core 2 Duo 2GHz, 2Gb RAM. Ex cept RL-MRS which is written in C++, all other algorithm s were imp lemented in MA TLAB. The same exp eriment as above was carried ou t with a microscopy image of the endothelial cell of the blood mi- crovessel walls; see Fig. 2. For th e NaiveGauss, An sGauss and our approach, th e dictionar y Φ contained th e wa velet orthog onal basis and the curvelet tight fr ame. The An s- Gauss a nd the Naiv eGau ss r esults ar e spoiled by artifacts and suf fer from a loss of photometry . RL-TV result shows a goo d restoration of small isolated details but with a dom- inating staircase-like artifacts. RL-MRS and our approach giv e very similar results althoug h an extra-effort could be made to better restor e tiny details. The quantitative mea- sures dep icted in Fig. 3 confirm this qualitati ve discussion. Finally , we app lied ou r algo rithm o n a real con focal microscopy image of neu rons. Fig. 4(a) depicts the ob - served image 2 using the GFP fluorescent pro tein. Fig. 4(b ) shows the restored image using our alg orithm with the or- thogon al wa velets. The ima ges ar e shown in log-scale fo r visual pu rposes. W e can no tice that the b ackgro und has been c leaned and some structures have reappea red. The spines are well resto red and part of the dendr itic tree is re- constructed , howe ver so me info rmation can be lost (see tiny holes). This can be imp roved using mor e re lev ant transform s. (a) (b) (c) (d) (e) (f) (g) (h) Fig. 1 . De con volutio n of a simulate d neuron (Intensity 6 30). (a) Original , (b) Blurred, (c) Blurred&noisy , (d) RL-TV , (e) Naiv eGauss, (f) AnsGauss, (g) RL-MRS, (h) Our Algorithm. 6. CONCLUSION In th is paper, we presented a sparsity- based fast iterative thresholdin g deconv olu tion algor ithm that take accou nts of the presence of Poisson noise. Com petitiv e results on con- focal micro scopy images with state- of-the-a rt algorithm s 2 Courtesy of the GIP Cyc ´ eron, Caen France. (a) (b) (c) (d) (e) (f) (g) (h) Fig. 2 . Decon vol ution of the ce ll image (Intensity 6 30). (a) Origi- nal, (b) Blurre d, (c) Blurred&noi sy , (d) RL-TV , (e) Nai veGa uss, (f) Ans- Gauss, (g) RL-MRS, (h) Our Algorithm. 50 100 150 200 250 10 −1 10 0 10 1 Intensity l 1 −error Our method NaiveGauss AnsGauss RL−MRS RL−TV Fig. 3 . Mean ℓ 1 -error of all algorithms as a function of the intensity le vel for the decon volut ion of the cell are shown. The combination of se veral transfor ms leads to some a dvantages, as we can easily adapt the d ictionary to the kind of imag e to restore. The param eter λ can be tr icky to find, and we are developing a m ethod helping to solve this issue. 7. REFERENCES [1] N. Dey et al., “ A decon volu tion method for conf ocal microscopy with total v ariation regulariz ation, ” in IEE E ISBI , 2004. [2] P . Sarder and A. Nehorai, “Dec onv olution Method for 3-D Fluo- rescenc e Micr oscopy Images, ” IEEE Sig. Pr o. , vol. 23, pp. 32–45, 2006. [3] J. Boutet de Mon vel et al, “Image Restora tion for Confocal Mi- croscop y: Improving the Limits of Decon v olution, with Applica - tion of the V isualiza tion of the Mammal ian Hearing Or gan, ” Bio- physical J ournal , vol. 80, pp. 2455–2470, 2001. [4] C. V onesch and M. Unser, “ A fast iterati ve thresholding algorithm for wa velet-r egulariz ed decon vol ution, ” IEEE ISBI , 2007. [5] I. Daubechies, M. Defri se, and C. De Mol, “ An iterati ve thresh- olding algorit hm for linear in ve rse problems with a sparsity con- straints, ” CP A M , vol . 112, pp. 1413–1541, 2004. Intensit y regime Method 6 5 6 30 6 100 6 255 Our method 0.21 0.76 2.39 4.79 Nai veGauss 0. 41 0.89 2.31 5.15 AnsGauss 2.07 7.56 21.25 50.41 RL-MRS 0.21 0.96 2.2 4.51 RL-TV 0.73 1.47 2.72 4.28 T able 1 . Mean ℓ 1 -error of all algorit hms as a function of the intensit y le vel for the decon volut ion of the neuron phantom. Method T ime (in s) Our method 2.7 Nai veGauss 1.7 AnsGauss 1.7 RL-MRS 15 RL-TV 2.5 T able 2 . Executi on time for the simula ted neuron using the o rthogonal wa velet transform [6] P . L . Combett es and V . R. W ajs, “Signal reco very by proximal forward -backward splitting, ” SI AM MMS , vol. 4, no. 4, pp. 1168– 1200, 2005. [7] G. T eschke , “Mu lti-frame representatio ns in linear in verse prob- lems with mixed m ulti-co nstraints, ” ACHA , vol. 22, no. 1, pp. 43– 60, 2007. [8] C. Chaux, P . L. Combettes, J. -C. Pesquet , and V . R. W ajs, “ A vari- ationa l formulati on for frame-based in verse problems, ” In v . P r ob . , vol. 23, pp. 1495 –1518, 2007. [9] M. J . Fadili and J .-L. Starck, “Sparse representati on-based image decon v olution by iterati ve thresholding, ” in ADA IV , Franc e, 2006, Else vier . [10] C. Chaux, L . Blanc-F ´ eraud, and J. Zerubia, “W av elet-based restora - tion met hods: applicati on to 3d confoca l m icroscop y images, ” in SPIE W avelets X II , 2007. [11] F-X Dup ´ e, M.J. Fadili, and J.-L. Starck, “ A proximal iteratio n for decon v olving poisson noisy ima ges using sparse representati ons, ” IEEE T rans. Im. Pr oc. , 2008, submitted. [12] J.-L. Starck and F . Murta gh, Astron omical Image and Data Analy- sis , Springer , 2006. [13] B. Zhang, J. Zerubia, and J.-C. Oliv o-Marin, “Gaussian approxi- mations of fluorescenc e microscope PSF models, ” OSA , 2006. (a) (b) Fig. 4 . Decon volu tion of a real neu ron. (a) Origi nal noisy , (b) Resto red with our algorit hm