Richardson–Lucy Algorithm With Total Variation Regularization for3D Confocal Microscope Deconvolution NICOLAS DEY,1LAURE BLANC-FERAUD,1*CHRISTOPHE ZIMMER,2PASCAL ROUX,3ZVI KAM,4
JEAN-CHRISTOPHE OLIVO-MARIN,2AND JOSIANE ZERUBIA1
1Ariana Group,INRIA/I3S,2004route des Lucioles–BP93,06902Sophia Antipolis,France
2Quantitative Image Analysis Group,Institut Pasteur,25–28rue du Dr.Roux,75015Paris,France
3Dynamic Imaging Platform,Institut Pasteur,25–28rue du Dr.Roux,75015Paris,France
4Department of Molecular Cell Biology,Weizmann Institute of Science,Rehovot,Israel76100
KEY WORDS image deconvolution;total variation regularization;Poisson noi;fluorescence confocal microscopy
ABSTRACT Confocal lar scanning microscopy is a powerful and popular technique for3D imaging of biological specimens.Although confocal microscopy images are much sharper than standard epifluorescence ones,they are still degraded by residual out-of-focus light and by Poisson noi due to
photon-limited detection.Several deconvolution methods have been propod to reduce the degradations,including the Richardson–Lucy iterative algorithm,which computes maximum likelihood estimation adapted to Poisson statistics.As this algorithm tends to amplify noi,regu-larization constraints bad on some prior knowledge on the data have to be applied to stabilize the solution.Here,we propo to combine the Richardson–Lucy algorithm with a regularization con-straint bad on Total Variation,which suppress unstable oscillations while prerving object edges.We show on simulated and real images that this constraint improves the deconvolution results as compared with the unregularized Richardson–Lucy algorithm,both visually and quanti-tatively.Microsc.Res.Tech.69:260–266,2006.V C2006Wiley-Liss,Inc.
INTRODUCTION
The confocal lar scanning microscope(CLSM)is an opticalfluorescence microscope that us a lar to scan and image point-by-point a sample in3D,and where a pinhole is ud to reject most out-of-focus light.The ability of the CLSM to image optical ctions explains its widespread u in biological rearch. Nevertheless,the quality of confocal microscopy images still suffers from two basic physical limitations. First,out-of-focus blur due to the diffraction-limited nature of the optics remains substantial,even though it is reduced compared to epifluorescence microscopy. Second,the
confocal pinhole drastically reduces the amount of light detected by the photomultiplier,lead-ing to Poisson noi.The images produced by the CLSM can therefore benefit from postprocessing by deconvolution methods designed to reduce blur and noi.
Several deconvolution methods have been propod for3D microscopy.In wide-field microscopy,Agard and Sedat(1983)and many others achieved out-of-focus deblurring by constrained iterative deconvolution.For confocal images,Tikhonov–Miller inverfilter(van Kempen et al.,1997),Carrington(van Kempen and van Vliet,2000a)and Richardson–Lucy(RL)algo-rithms(Lucy,1974;Richardson,1972)have also been applied.The latter has been ud extensively in astro-physical or microscopic imaging(van Kempen et al., 1997),and is of particular interest for confocal micros-copy,becau it is adapted to Poisson noi.An impor-tant drawback of RL deconvolution,however,is that it does not converge to the solution becau the noi is amplified after a small number of iterations.This n-sitivity to noi can be avoided with the help of regula-rization constraints,leading to much improved results. van Kempen et al.have applied a RL algorithm using energy-bad regularization to biological images(van Kempen and van Vliet,2000a).The term bad on the Tikhonov–Miller expression suppress noi amplifi-cation,but it has the drawback to smooth edges.
In this study,we propo to u the total variation (TV)as a regularization term.The main advantages of the TV are that it prerves the edges in the image and smoothes homogeneous areas.It has been propod for 2D optical gray level image restoration in Rudin et al. (1992).This kind of nonlinear image restoration has been extensively ud in image processing inver problems,for the last10years,and has generated mathematical studies on the Bounded Variation(BV) function space(Ambrosio et al.,2000).In this study,we propo to apply this nonlinear regularization to the deconvolution of3D confocal microscopy images.Con-sidering the statistical Poisson nature of the noi,this leads to a RL algorithm with nonlinear regularization. The paper is organized as follows.Thefirst ction prents the image formation model,by describing the point spread function and the noi distribution.In the cond ction,wefirst recall the RL algorithm in
*Correspondence to:Laure Blanc-Feraud,Ariana Group,INRIA/I3S,2004 route des Lucioles–BP93,06902Sophia Antipolis,France.
E-mail:Laure.Blanc_Feraud@sophia.inria.fr
Contract grant sponsors:INRIA ARC De´MiTri,P2R Franco-Israeli Program. Received28July2005;accepted in revid form13October2005
DOI10.1002/jemt.20294
Published online in Wiley InterScience(www.).
V C2006WILEY-LISS,INC.
MICROSCOPY RESEARCH AND TECHNIQUE69:260–266(2006)
the context of microscope imaging,and then describe the nonlinear TV regularized RL algorithm.Next,we prent the results of the propod algorithm on syn-thetic and real data,and compare them with the stand-ard algorithm without regularization,both in terms of visual image quality and a numerical criterion.The last ction concludes the paper.
IMAGE FORMATION MODEL
The image formation can be modeled with the Point Spread Function(PSF)and the statistical distribution of the detection noi.Such a general model can be written as:
i¼}ðoÃhÞð1Þ
where i is the obrved(measured)image,o is the object we want to retrieve(corresponding to the distri-bution offluorescent markers inside the specimen),h is the PSF,and}models the noi distribution.
Noi梅条肉
In a confocal microscope,a pinhole is ud to reject most out-of focus light.The amount of light reaching the detector is therefore low,and the noi statistics is well described by a Poisson process.The distribution of i at point s,i(s),knowing(oÃh)(s),is a Poisson distri-bution of the parameter(oÃh)(s):
PðiðsÞ=ðoÃhÞðsÞÞ¼½ðoÃhÞðsÞ iðsÞeÀðoÃhÞðsÞ
iðsÞ!
ð2Þ
Assuming that the noi is spatially uncorrelated,the statistics of the obrved image i,knowing the object o, (and the PSF h which is considered here as a parame-ter)is the likelihood distribution given by:
Pði=oÞ¼
Y
s e S ½ðoÃhÞðsÞ iðsÞeÀðoÃhÞðsÞ
iðsÞ!
!
ð3Þ
where s is the coordinate of a3D voxel and S is the total t of voxel coordinates in the imaged volume.
Confocal Point Spread Function
Under the assumptions of a translation invariant PSF,incoherent imaging,and monochromatic light (excitation and emission),we u a well-known image formation model of the CLSM,which has been pro-pod by van Kempen and van Vliet(1999),and is derived from Sheppard and Cogswell(1990).For a defocusing Z,the3D PSF(Sheppard and Cogswell, 1990)is given by:
hðX;Y;ZÞ¼j A RðX;YÞÃ^P k emðX;Y;ZÞj2
3j^P k exðX;Y;ZÞj2
ð4Þ
where k ex and k em are the excitation and emission wavelengths,respectively.The pinhole(Born and Wolf, 1999;Goodman,1996)is reprented by A R(X,Y),and ^P
幼儿便秘k is the2D Fourier transform of the pupil function P k given by:
P kðU;V;ZÞ¼P q
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
U2þV2
领导核心
p
e2i p k WðU;V;ZÞð5Þ
Here W(U,V,Z)¼1
2
Z(1Àcos2a)is the aberration pha derived from Goodman(1996).It gives a measure of the wavefront path difference between a focud and a defo-cud beam(Born and Wolf,1999;Stokth,1969).q¼NA
2k
is the lateral cutoff frequency,and NA¼n o sin a the numerical aperture,which is related to the cone of light collected by the objective and the immersion medium re-fractive index n o.The model takes into account thefinite size of the pinhole,which is ud for the image acquisi-tion.This PSF model will be ud both for the creation of the simulated images and for the deconvolution.
In van Kempen and van Vliet(1999)and van Kempen and van Vliet,(2000a),the authors u an image forma-tion model for confocal microscopy,which is slightly dif-ferent from Eq.(1):namely i¼q(oÃhþb),where b models the background signal.In a previous work(Dey et al.,2004a,b),we showed that the u of a background term in the model is not necessary if a well suited regu-larization term is included.In fact,it happens that using a background model is equivalent to imposing a support constraint to regularize the solution.We also showed that when using a regularization term,
no im-provement is reached by adding the background term, and that it is therefore not needed.
CONFOCAL MICROSCOPY IMAGE
DECONVOLUTION
The Multiplicative RL Algorithm
The RL algorithmfinds o from the obrvation i, knowing the PSF h,by maximizing the likelihood dis-tribution(Eq.(3))with respect to o.The RL algorithm can be derived by using the Expectation–Maximization (EM)algorithm(Dempster et al.,1977),or by minimiz-ing the functionalÀlog p(i|o),using a multiplicative type algorithm.Let us briefly describe the cond approach.MinimizingÀlog p(i|o)is equivalent to mini-mizing J1(o)given by:
J1ðoÞ¼
X
s
ðÀiðsÞlog½ðoÃhÞðsÞ þðoÃhÞðsÞÞð6Þ
Since the functional J1(o)is convex in o,arching for a minimum is equivalent to arching for a zero of the gradient of J1(o)which results in solving
iðsÞ
ðo kÃhÞðsÞ
ÃhðÀsÞ¼1ð7Þ
Let o k(s)be the estimate of o at iteration k.Then one iteration of a multiplicative gradient type algorithm (for more details e Lucy(1974)and Richardson1972) or Appendix C of Dey et al.(2004a))is given by:
o kþ1ðsÞ¼
iðsÞ
ðo kÃhÞðsÞ
ÃhðÀsÞ
o kðsÞð8Þ
261
TOTAL VARIATION BASED DECONVOLUTION
which defines one iteration of the RL algorithm.An important property of this scheme is that it ensures nonnegativity if thefirst estimate(here taken as the measured3D image,scaled by a value equal to the mean of the image stack)is nonnegative.However, becau the inversion problem is ill-pod and the Max-imum Likelihood(ML)estimator is nonregularized,the solution given by Eq.(8)when k?þ?consists of noi only,becau of noi amplification.To obtain good results,the algorithm is stopped before the noi is amplified.In order to improve the trade-off between a good deconvolution with sharp edges and noi amplifi-cation,we propo to regularize the solution by mini-mizing its TV.
Total Variation Regularization
In a Bayesian approach,a regularization is obtained by Maximizing the a posteriori distribution(MAP esti-mator)P(o/i)¼P(i/o)3P(o)/P(i)with a suitable a priori distribution P(o).Tikhonov–Miller(TM)(van Ke
mpen and van Vliet,1999;van Kempen and van Vliet,2000b)regularization is the most popular algo-rithm in3D image deconvolution.TM regularization is designed to avoid noi amplification by smoothing, causing smear of the object edges.de Monvel et al., (2003)ud a RL algorithm with maximum entropy regularization.This regularization gives better re-sults than that of TM regularization,in the n that the edges are less smoothed.However,the trade-off between good regularization with a good reduction of the noi,and edge prervation still exists.This compromi can be improved by using wavelet coeffi-cient thresholding,as propod by de Monvel et al. (2003).
Here,we introduce an other a priori constraint on the object by adding to the energy J1,a regularization term J reg,defined by the TV of the solution o as in Rudin et al.(1992):
J regðoÞ¼k TV
X
s
jr oðsÞjð9ÞThe total functional to be minimized is then:
J1ðoÞþJ regðoÞ¼
X
s
ðÀiðsÞlog½ðoÃhÞðsÞ
þðoÃhÞðsÞÞþk TV
X
s jr oðsÞj
ð10Þ
Using the L1norm over!o rather than the L2norm as in TM regularization(van Kempen and van Vliet, 2000b)allows to recover a smooth solution with sharp edges.It can be shown(Rudin et al.,1992)that the smoothing process introduced by J reg acts only in the direction tangential to the image level lines and not in the orthogonal direction,so that edges are prerved. The derivative of J reg with respect to o is a nonlinear
term@
@o J reg¼Àk TV div
谷物大脑
r o
jr o j
where div stands for the
divergence.We minimize Eq.(10)using the multiplica-tive gradient-bad algorithm(or equivalently by using EM algorithm for the penalized criterion of Eq.(10)),and we adopt an explicit scheme,as in Green (1990),defined by:
o kþ1ðsÞ¼
iðsÞ
ðo kÃhÞðsÞ
ÃhðÀsÞ
3o k
ðsÞ
1Àk TV div r okðsÞ
jr okðsÞj
ð11Þ
Numerically we noticed that the regularization param-eter k TV should be neither too small nor too large:if k TV<10À6,RL is dominated by the data model;if k TV $1,RL is dominated by the regularization term.For larger k TV,the denominator of Eq.(8)can become zero or negative.This must be avoided becau small de-nominators create points of very high intensity,which are amplified at each iteration.A negative value vio-lates the nonnegativity constraint.We u a param-eter k TV¼0.002for computations,and iterations are stopped at convergence,defined when the relative dif-ference between two images is less than a threshold (10À5in the experiments prented below).
RESULTS
Evaluation Criteria
藤县中学
In the following,we report results on three synthetic (simulated)objects and two types of real specimens.To quantify the quality of the deconvolution we ud the I-divergence criterion,which is the best measure in the prence of a nonnegativity constraint(Csisza´e´r,1991). The I-divergence between two3D images A and B is defined by:
I A;B¼
X
s e S
A s ln
A s
B s
我的家园ÀðA sÀB sÞ
ð12Þ
分钟换算小时
The I-divergence is ,I A,B=I B,A.In the experiments with synthetic data,image A is the known undegraded initial image(the one we want to retrieve)and image B is the restored one.The I-diver-gence then quantifies the difference between the ideal image and the reconstructed one.Ideally,a perfect deconvolution should end with I-divergence equal to zero.This numerical criterion cannot be computed on results from real data becau the undegraded object is not known.
In order to display the visual quality,we prent cross ctions of the objects before and after deconvolu-tion with the RL algorithm and with the TV regular-ized RL algorithm.
Results on Synthetic Data
We show results on three different synthetic objects, shown in Figures1–3(column a).The objects are artificially blurred and corrupted by Poisson noi using the image formation model described in the c-tion of image formation model.The resulting degraded images are shown in Figures1–3(column b).Deconvo-lution results of the degraded images using standard RL are shown in Figures1–3(column c).For standard RL,we ud the number of iterations that achieves minimum I-divergence.The deconvolution results
262N.DEY ET AL.
using RL with TV regularization (RL–TV)are shown in Figures 1–3(column d).
The first simulated object is a full cylinder (Fig.1a).It is a single object,prenting no corner (except at the top and bottom)and which is homogeneous (there are only two intensity values:the background and the cylinder).Figure 1c shows the result of the deconvolution using the standard RL algorithm:in lateral and axial views there is still some blur,and intensity oscillation artifacts are clearly apparent.The RL–TV algorithm greatly improves the results (Fig.1d):the cylinder is sharp (the edges are sharp in lateral and axial views)and there are no intensity oscillations (e the corresponding in-tensity profile in Figure 6a).The axial view shows that the heights of the deconvolved cylinder (Fig.1d)and of the initial cylinder (Fig.1a)are similar.The final I-divergence is 0.766for standard RL and 0.220for RL–TV ,resulting in a more than 3-fold improvement.
Figure 2a prents another 3D synthetic image com-pod of five geometrical shapes.The are more com-plex objects with corners and homogeneous regions of different intensity levels.The deconvolution result,with standard RL,is shown in Figure 2c.The different shapes are still blurred and appear thinner than before degradations.Many intensity oscillations appear inside the shapes.All of the artifacts are avoided when using the RL–TV algorithm (Fig.2d):all shapes are sharp,the thickness is prerved,and there are no more intensity oscillations (e the intensity profile in Fig.6b).
The cross-ction displays a slight effect of rounding corners.In the axial view,the borders are not perfectly sharp.Nevertheless,the final I-divergence is 0.691for the pro-pod method,which reprents nearly a 2-fold im-provement compared with standard RL (1.365).
In Figure 3a,we prent a 3D synthetic object com-bining texture and fine structures.The object is inho-mogeneous and the fine structure exhibits sharp tran-sitions.The degradation (Fig.3b)is sufficiently strong to almost hide the details.The deconvolution ob-tained with standard RL still shows strong axial blur.With RL–TV ,this blur is almost completely removed,and the fine structures are discernable again.However,a drawback of TV regularization is apparent on the XY ction:a stair-casing effect.It is due to the
nonderiv-
Fig.1.Deconvolution of a binary cylinder.(a )the original 3D syn-thetic image;(b )the same image,degraded by blur and Poisson noi according to the image formation model;(c )deconvolution by stand-ard RL (I-divergence is 0.766);(d )deconvolution by RL–TV with k TV ¼0.002(I-divergence is 0.220).The top row shows the lateral (XY )
ction through the center of the image stack.The bottom row shows the axial (YZ )ction corresponding to the dotted line shown in (a).The image size is 1283128364with voxel size 30330350nm in X ,Y ,and Z
.
Fig.2.Deconvolution of a composite image.(a )original image;(b )degraded by blur and Poisson noi;(c )deconvolved with standard RL;(d )deconvolved with RL–TV.The objects exhibit different inten-sities:255for the cylinder,221for the annulus,170for the cross,102for the triangle,238and 102for the two bars,10for the background.Final I-divergence is 1.365for standard RL (c )and 0.691for RL–TV with k TV ¼0.002(d ).The size of the stack is 1283128364with voxel size 30330350nm in X ,Y ,and Z .
263
TOTAL VARIATION BASED DECONVOLUTION
ability of the absolute value in zero,and appears mainly in parts of the image which contain oscillations of small amplitudes,as textures.The intensity profiles in Figure 6c do not show significant differences be-tween standard RL and RL–TV results.Final I-diver-gence is 0.462for standard RL and 0.404for RL–TV ,which is an improvement by 12.6%.
Results on Real Data
To evaluate the performances of our method on real data,we cho two types of objects,beads and shells,who size and shape are well characterized and which can therefore be ud as bona fide references.The
beads are FocalCheck TM F-24634(Molecular Probes)fluorescent microspheres,with a diameter of 6l m.The shells are 15l m in diameter,with a fluorescent layer of thickness between 0.5and 0.7l m.The fluorescence excitation wavelength is k ex ¼488nm;the emission wavelength is k em ¼520nm;the pinhole size was t to 1Airy unit;in object space,the XY dimension of each voxel was 89nm and the Z dimension 230nm.The microscope was a Zeiss LSM 510,equipped with an immersion oil Apochromat 633with numerical aper-ture NA ¼1.4objective,and an internal magnification of 3.33.Image acquisition was performed using the Zeiss Vision
software.
泪珠儿Fig.3.Deconvolution of a textured image with fine structure.(a )original image;(b )degraded by blur and Poisson noi;(c )deconvolved with standard RL;(d )deconvolved with RL–TV.Final I-divergence is 0.462for standard RL (c )and 0.404for RL–TV with k TV ¼0.002(d ).Note that axial blur is much reduced with RL–TV compared to RL.The image size is 12831283
32.
Fig.4.Deconvolution of a spherical shell image.(a )shows the raw image,(b )the result of standard RL deconvolution,(c )the result of RL–TV deconvolution,with k TV ¼0.002.RL–TV deconvolution allows a better estimate of the shell thickness than RL (e text).The image size is 25632563128with voxel size 893893230nm in X ,Y ,and
Z .First row reprents the central XY image,and cond line the cen-tral axial ction of this stack corresponding to the dotted line.There is no resampling to obtain the same scale in X ,Y ,and Z ,only a resiz-ing of the image.
264
N.DEY ET AL.
