Roerdink, J.B.T.M.

English

ARTS repository - University of Groningen

   University of GroningenNonlinear Image Restoration in Confocal MicroscopyRoerdink, J.B.T.M.Published in:EPRINTS-BOOK-TITLEIMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite fromit. Please check the document version below.Document VersionPublisher's PDF, also known as Version of recordPublication date:1994Link to publication in University of Groningen/UMCG research databaseCitation for published version (APA):Roerdink, J. B. T. M. (1994). Nonlinear Image Restoration in Confocal Microscopy. In EPRINTS-BOOK-TITLE University of Groningen, Johann Bernoulli Institute for Mathematics and Computer Science.CopyrightOther than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of theauthor(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).The publication may also be distributed here under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license.More information can be found on the University of Groningen website: https://www.rug.nl/library/open-access/self-archiving-pure/taverne-amendment.Take-down policyIf you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediatelyand investigate your claim.Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons thenumber of authors shown on this cover page is limited to 10 maximum.Download date: 20-06-2022Nonlinear Image Restoration in Confocal Microscopy J.B.T.M. RoerdinkDept. of Computing Science, University of GroningenP.O. Box 800, 9700 AV Groningen, The NetherlandsTel. +31-50-633931; Fax +31-50-633800; Email roe@cs.rug.nl1 IntroductionA major problem in 3D imaging by a CSLM (confocal scanning laser microscope) in the (epi)uorescencemode is the darkening of the deeper layers due to scattering and absorption of excitation and uorescencelight. In [2] we developed a new restoration method, called the `FFT-method', to correct for these eectsby deriving a correction factor to the standard restoration in the form of a 3D convolution of the measuredsignal, which can be eÆciently computed by the use of the Fast Fourier Transform (FFT). In this way,the complexity of computation is reduced to O(Nz logNz), where Nz is the number of vertical layers.We also compared the computational eÆciency of our algorithm with the iterative layer method `withcondensation' developed in [4], which has complexity O(N2z ). For spatially varying image densities therestoration quality using our method was found to be a little poorer than in the layer method.This paper gives a succinct summary of an extension to this approach as published in [3]. There itis shown that the accuracy of the FFT-method can be improved by rst order moment and cumulantestimators leading to a nonlinear integral equation for the unknown uorescent density, which is solvedby an iterative method. The new estimators, the moment estimator in particular, are more accurate thanthe layer method. Since the computations involve only discrete 3D convolutions computable by the FFT,the advantage in computational eÆciency over the layer method is retained.2 The CSLM transformWhen a CSLM is operating in the uorescence mode a laser beam is focussed upon a point r = (x; y; z)in the object, see Fig. 1. Here the z-direction is chosen along the optical axis. The rays converging to theobject point are contained in a `light cone' with angle !, called the `semi-aperture angle'. The radiationabsorbed at the point in focus is uniformly reemitted as uorescent radiation and the part which travelsback the same route as the incoming radiation is detected. While travelling through the sample theradiation is attenuated with an attenuation factor which is spatially dependent. Then, assuming that (i)the attenuation of the excitation light equals that of the uorescence light, and (ii) that this attenuationcoeÆcient is proportional to the (unknown) uorescent density (r), one can derive the following nonlinearintegral transform (`CSLM-transform'):f(r) = (r) f (r)b(r); (1)Here f(r) is the measured uorescent intensity,f (r) := Cf Z !0 d Z 20 d sin  cos  exp  "Z z0 dz0cos  (r̂) (2)is the forward attenuation factor, andb(r) := Cb Z !0 d Z 20 d sin  exp "Z z0 dz0cos  (r̂) (3)In: Proc. Computing Science in the Netherlands, 21-22 november, Utrecht, 1994, pp. 253-259. Postscript versionobtainable at http://www.cs.rug.nl/~roe/.xyzddyzdxP(x,y,z)Rωφ      θFigure 1: CSLM geometry of light cone with apex at a point P (x; y; z) in the object. R: radius ofspherical bundle; !: semi-aperture angle; (; ): polar angles of light ray; dz: depth of the sample. Theoptical axis coincides with the z-axis.is the backward attenuation factor. In these equations " is a proportionality constant and Cf and Cb arenormalization constants. Here r̂ is the vectorr̂(r; ; ; z0) = (x+ (z   z0) tan  cos; y + (z   z0) tan  sin; z0): (4)The quantity " dz is a measure for the degree of attenuation, where dz is the depth of the sample.Inversion of the CSLM transformBy a perturbation expansion in the parameter ", which is tantamount to the assumption of weak atten-uation, we derived in [2] the following approximation ~(r) for the uorescent density,~(r) = f(r) f1 + " c(r)g ; (5)where c(r) is the convolution integralc(r) = Z 1 1 Z 1 1 Z 1 1 dr0 (r0) f(r  r0); (6)with (r) the space-invariant kernel given by(x; y; z) = (Cf z(x2+y2+z2)3=2 + Cb 1x2+y2+z2 ; 0  z  dz ; x2 + y2  (z tan!)20 elsewhere : (7)Numerical computationFor numerical computation, the integral (6) is discretized on a grid of Nx Ny Nz voxels, each voxelbeing a box of dimensions Æx; Æy; Æz in the x-, y- and z-directions. Then the approximation (5) is replacedby Rijk = Fijk(1 + " Cijk); (8)where Cijk = Mx2Xi0= Mx2 +1 My2Xj0= My2 +1 MzXk0=1Ki0j0k0 Fi i0 ;j j0 ;k k0 ; (9)with Rijk ; Cijk ;Kijk and Fijk the discrete counterparts of ~(x; y; z); c(x; y; z); (x; y; z) and f(x; y; z),respectively, with (i; j; k) in the index set II := f(i; j; k) : i = 1; ::; Nx; j = 1; ::; Ny; k = 1; ::; Nzg, andarray elements are dened to be zero when the indices are not in the index set II . The discrete convolution(9) can be computed eÆciently by FFT methods [2].3 Iterative algorithmsBy applying moment and cumulant expansions of characteristic functions [3] one obtains equations forthe `moment approximation' (m)(r) and the `cumulant approximation' (c))(r) in the form(c)(r) = f(r) exp h" (  (c))(r)i ; (10)and (m)(r) = f(r) h1  " (  (m))(r)i 1 ; (11)where   () denotes the convolution of the functions  and (),  = m; c, and the kernel  is identicalto that in (7). We assume that " is chosen small enough for the inverse in (11) to exist.After discretization, both (10) and (11) lead to a nite system of nonlinear equations of the formRijk = Fijk G(K R)ijk; (i; j; k) 2 II; (12)where G(x) = exp("x) for the cumulant approximation and G(x) = (1  "x) 1 for the moment approxi-mation, respectively, with K R the discrete convolution of the 3D arrays K and R.The equations (12) can be solved by Picard iteration [1], with Fijk as the initial estimate:R(0)ijk = Fijk ; (13)R(n)ijk = Fijk G(K R(n 1))ijk ; n = 1; 2; 3; ::: : (14)Each iteration step involves the computation of the discrete convolutionKR(n 1) of the estimate R(n 1)of the previous iteration (with the same convolution kernel K) which can be eÆciently computed by theFFT. The rst iterate of Eq. 14 with G(x) = 1 + "x coincides with the discrete analogon of Eq. 5 and isthe approximation used in [2].It can be shown that (for proofs, see [3]): The equations (12) have a unique solution; The iterates R(n)ijk of Eq. 14 converge in a nite number of Nz steps towards the unique solution ofEq. 12; The convergence is monotonous, that is, R(n)ijk  R(n 1)ijk ; The estimators (m) and (c) satisfy the inequality(c)  (m): (15)The restoration procedureThe complete procedure of reconstructing the uorescent density from measured data contains the fol-lowing steps:1. Read the measured data Fijk ; i = 1; ::; Nx; j = 1; ::; Ny; k = 1; ::; Nz.2. Determine the appropriate value of the attenuation constant ".3. Iteratively compute R(n)ijk = Fijk G(K R(n 1))ijk ; n = 1; 2; 3; :::; (16)where R(0)ijk = Fijk , and G(x) = (1   "x) 1 or G(x) = exp("x) for the moment and cumulantestimator, respectively.In each iteration the convolution of the previous estimate is computed by means of the FFT (using thesame kernel K as dened above).Input parameters of the algorithm are the dimensions Nx; Ny; Nz of the data array, the scanning stepsÆx; Æy; Æz, the semi-aperture angle ! and the attenuation constant ". For the determination of the correctvalue of " one may resort to a calibration experiment in which a homogeneous test sample is used [4]. Inthis paper we will consider test densities for which the value of " is known.4 Restoration of a test imageIn this section we consider a test image consisting of a number of circles and ellipses already used in[2]. Signal data Fijk were generated by numerically computing the integrals in (2){(3) for a number ofequidistant 3D positions. The parameters were chosen as follows: dx = dy = 1:0; dz = 0:1, Nx = Ny =128, Nz = 8, ! = 1:04719. We computed the relative root mean square errorE(z) :=  PNxx=1PNyy=1f(x; y; z)  ~(x; y; z)g2PNxx=1PNyy=1f(x; y; z)g2 ! 12 ;between original density  and restored density ~ at each plane z = constant. Computations wereperformed on a SPARC workstation (35 Mhz, 26 MIPS), taking about one minute per iteration step."z signal error iter = 1 iter = 2 iter = 30:0000 0:000 0:000 0:000 0:0000:0625 0:116 0:004 0:004 0:0040:1250 0:218 0:006 0:016 0:0160:1875 0:305 0:025 0:041 0:0450:2500 0:382 0:065 0:065 0:0900:3125 0:450 0:127 0:075 0:1530:3750 0:509 0:200 0:056 0:2360:4375 0:560 0:278 0:026 0:330Table 1: Signal error and restoration errors by the moment estimator (m) after one, two and threeiterations as a function of the eective depth "z."z signal error iter = 1 iter = 2 iter = 30:0000 0:000 0:000 0:000 0:0000:0625 0:116 0:010 0:010 0:0100:1250 0:218 0:035 0:022 0:0220:1875 0:305 0:080 0:047 0:0460:2500 0:382 0:141 0:089 0:0840:3125 0:450 0:213 0:146 0:1360:3750 0:509 0:286 0:213 0:1980:4375 0:560 0:357 0:283 0:265Table 2: Signal error and restoration errors by the cumulant estimator (c) after one, two and threeiterations as a function of the eective depth "z.Figure 2: Restoration of the circle image. First row: the attenuated test images; second row: restorationby the cumulant estimator (iter = 2); third row: restoration by the moment estimator (iter = 2). In eachrow, the rst, fourth and seventh layer is displayed from left to right. The original image in each layer isidentical to the rst image in row 1.Results are shown in Table 1 for the moment estimator (11) and in Table 2 for the cumulant estimator(10). In the case of the moment estimator the errors rst decrease and then start to grow again afterthe second iteration. This is due to the fact that the initial estimate f(r) is smaller then the exactdensity (r), so that at rst the iterates underestimate the true solution. Because of the monotonicityproperty mentioned in Section 3 the iterates always increase so that (if the solution (m) is larger thanthe true , which is apparently the case here) they will start to overestimate the true density. Thecumulant estimator (c) underestimates the true density, and the values were stable within an accuracyof three digits after the third iteration. For comparison we give in column 2 of the Tables the error beforerestoration, denoted by `signal error' and computed according to (4) with ~ replaced by f .Comparing with the numbers in Table 4 of [2], we conclude that both the moment estimator withiter = 1; 2 and cumulant estimator with iter  2 are more accurate than the layer method of [4] whichgives a restoration error of 0.301 at the deepest layer. From the Tables it is clear that the momentestimator, when run to convergence, overestimates the exact image densities. The rst iterate, however,underestimates the exact values. Therefore, in case of the moment estimator, we take the reconstructioncorresponding to the intermediate value iter = 2, which gives the best results. In Fig. 2 we show thecorresponding restored images. We rescale the values of the densities to make sure that they occupy thecomplete grey-scale, which consists of the set of integer values from 0 to 255. In order to avoid that a fewoutliers cause a large visual degradation of the resulting images, we constrained the approximate solutions~ to lie between the known lower and upper bounds, i.e. 0  ~  1. In each row, the rst, fourth andseventh layer is displayed from left to right, out of a total of 8 depth layers. Since the exact density (r)does not depend on z, the original image in each layer is identical to the rst image in row 1. The rstrow contains the attenuated test images f(r), the second row the restoration by the cumulant estimatorand the third row the restoration by the moment estimator, both after two iterations. The images whichare restored by the cumulant estimator are virtually identical to those of the layer method of Visser et al.[4], cf. Figure 4 of [2]. Clearly the largest improvement in restoration quality has been obtained by usingthe moment estimator. The reconstruction is not perfect, however: the central regions in the centers ofthe light circular regions are slightly overestimated. The calculations have been repeated for other testimages, leading to similar conclusions: the reconstruction errors are smallest when using the momentestimator with iter = 2, but the reconstructed images still show some dierences when compared tothe original images, see [3]. Nevertheless, a considerable improvement in restoration accuracy has beenobtained by using the estimators developed here, which in addition are eÆciently computable by usingFFT methods.5 Restoration of a real CSLM imageThe method of this paper was applied to a CSLM image (Nx = Ny = 256, Nz = 8) of a geological sampleconsisting of sandstone cavities lled with a uorescent oily substance, also considered in [4] and [2].The result is shown in Fig. 3. Again the best result is obtained by the moment estimator, which showssubstantially more detail than the original image. This shows the practical usefulness of the method. Wedetermined the value of " in this case so as to give a visually satisfactory restoration.Figure 3: Restoration of the sandstone CSLM image. First row: the attenuated CSLM images; secondrow: restoration by the cumulant estimator (iter = 2); third row: restoration by the moment estimator(iter = 2). In each row, the rst, fourth and seventh layer is displayed from left to right.6 SummaryIn this paper we describe the method of [3], which is a renement of the approach developed in [2],for attenuation correction in Fluorescence Confocal Microscopy using Fast Fourier Transform methods.This approach, valid for weak attenuation, consists in multiplying the measured uorescent intensityby a correction factor involving a convolution integral of the measured signal, which can be computedeÆciently by an FFT-based algorithm. By a statistical reformulation of the problem it is possible to deriverst order moment and cumulant estimators leading to a nonlinear integral equation for the unknownuorescent density, which can be solved by an iterative method. In each iteration the convolution of theprevious estimate is computed by means of the FFT. The rst iterate of the new estimators coincides withthe approximation used in [2] for very weak attenuation. It turns out that the moment estimator withtwo iterations gives the best results, which are more accurate than the layer method of [4]. Since onlytwo iterations are needed, the advantage in computational eÆciency over the layer method is retained.We conclude therefore that the combined results of [2] and [3] provide an eÆcient and accurate methodfor attenuation correction in confocal microscopy.References[1] Ortega J.M. and W.C. Rheinboldt (1970). Iterative Solution of Nonlinear Equations in SeveralVariables. Ac. Press, New York.[2] Roerdink J.B.T.M. and M. Bakker (1993). An FFT-based method for attenuation correction inuorescence confocal microscopy. J. Microscopy 169, 3{14.[3] Roerdink J.B.T.M. (1994). FFT-based methods for nonlinear image restoration in confocal mi-croscopy, J. Math. Imaging and Vision 4(2), 199{207.[4] Visser T.D., F.C.A. Groen and G.J. Brakenho (1991). Absorption and scattering correction inuorescence confocal microscopy. J. Microscopy 163, 189-200.

Nonlinear Image Restoration in Confocal Microscopy

https://pure.rug.nl/ws/files/3296008/1994ProcCompSciNLRoerdink.pdf

Nonlinear Image Restoration in Confocal Microscopy

Abstract

Similar works

Full text

Available Versions

ARTS repository - University of Groningen