Numerical refocusing in digital holographic microscopy with extended-sources illumination Matˇej T´ycˇ ,∗ Luk´asˇ Kvasnica, Michala Slab´a, and Radim Chmel´ık Institute of Physical Engineering, Faculty of Mechanical Engineering, Brno University of Technology, Technick´a 2, 616 69 Brno, Czech Republic CEITEC — Central European Institute of Technology, Brno University of Technology, Technick´a 10, 616 00 Brno, Czech Republic ∗ [email protected]

Abstract: Numerical refocusing can be seen as a method of compensating the defocus aberration based on deconvolution by inverse filtering [1] in digital holographic microscopy (DHM). It is well-understood in cases when a coherent (ie point and monochromatic) light source such as a collimated laser beam is used [2]. This paper extends the theory to the case of illumination by a quasi-monochromatic extended (spatially incoherent) source. Refocusing methods for spatially incoherent illumination are derived and benefits of this type of illumination are demonstrated. We have proved both theoretically and experimentally that coherent-based refocusing gives incorrect results for extended-source illumination, while results obtained using the newly derived method are correct. © 2013 Optical Society of America OCIS codes: (090.1995) Digital holography; (110.1830) Microscopy; (100.1830) Deconvolution; (180.9600) Three-dimensional microscopy; (110.4980) Partial coherence in imaging.

References and links 1. M. T´ycˇ and R. Chmel´ık, “Numerical refocusing of planar samples unlimited,” Proc. SPIE 7746 774620 (2010). 2. F. Dubois, L. Joannes, and J. Legros, “Improved three-dimensional imaging with a digital holography microscope with a source of partial spatial coherence,” Appl. Opt. 38, 7085–7094 (1999). 3. L. Lovicar, J. Komrska, and R. Chmel´ık, “Quantitative-phase-contrast imaging of a two-level surface described as a 2D linear filtering process,” Opt. Express 18, 20585–20594 (2010). 4. R. Barer, “Interference microscopy and mass determination,” Nature (London) 169, 366–367 (1952). 5. Y. Cotte, F. Toy, C. Arfire, S. Kou, D. Boss, I. Bergo¨end, and C. Depeursinge, “Realistic 3D coherent transfer function inverse filtering of complex fields,” Biomed. Opt. Express 2, 2216–2230 (2011). 6. F. Dubois, C. Yourassowsky, N. Callens, C. Minetti, and P. Queeckers, “Applications of digital holographic microscopes with partially spatial coherence sources,” JPCS, 139 (IOP Publishing, 2008), p. 012027. 7. T. Kozacki and R. J´oz´ wicki, “Digital reconstruction of a hologram recorded using partially coherent illumination,” Opt Commun 252, 188–201 (2005). 8. T. Slab´y, P. Kolman, Z. Dost´al, M. Antoˇs, M. Loˇst’´ak, and R. Chmel´ık, “Off-axis setup taking full advantage of incoherent illumination in coherence-controlled holographic microscope,” Opt. Express 21, 14747–14762 (2013). 9. F. Dubois, M. Novella Requena, C. Minetti, O. Monnom, and E. Istasse, “Partial spatial coherence effects in digital holographic microscopy with a laser source,” Appl. Opt. 43, 1131–1139 (2004). 10. F. Dubois, O. Monnom, C. Yourassowsky, and J–C. Legros, “Border processing in digital holography by extension of the digital hologram and reduction of the higher spatial frequencies,” Appl. Opt. 41, 2621–2626 (2002). 11. F. Dubois, C. Schockaert, N. Callens, and C. Yourassowsky, “Focus plane detection criteria in digital holography microscopy by amplitude analysis,” Opt. Express 14, 5895–5908 (2006). 12. O. Carriere, J. Hermand, and F. Dubois, “Underwater microorganisms observation with off-axis digital holography microscopy using partially coherent illumination,” in “OCEANS 2011,” (IEEE, 2011), pp. 1–7.

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28258

13. P. Kolman and R. Chmel´ık, “Coherence-controlled holographic microscope,” Opt. Express 18, 21990–22003 (2010). 14. R. Chmel´ık, “Three-dimensional scalar imaging in high-aperture low-coherence interference and holographic microscopes,” J. Mod. Opt. 53, 2673–2689 (2006). 15. J. W. Goodman, Introduction to Fourier optics (McGraw-Hill, 1996), international editions 1996 16. Y. Cotte, M. Toy, N. Pavillon, and C. Depeursinge, “Microscopy image resolution improvement by deconvolution of complex fields,” Opt. Express 18, 19462–19478 (2010). 17. R. Chmel´ık, Z. Harna, “Parallel-mode confocal microscope,” Opt. Eng. 38, 1635–1639 (1999). 18. L. Lovicar, L. Kvasnica, and R. Chmel´ık, “Surface observation and measurement by means of digital holographic microscope with arbitrary degree of coherence,” Proc. SPIE , 7141 p. 71411S (2008). 19. Y. Geerts, M. Steyaert, and W. Sansen, Design of multi-bit delta-sigma A/D converters, vol. 686 (Springer, 2002). 20. F. Harris, “On the use of windows for harmonic analysis with the discrete Fourier transform,” Proc. IEEE 66, 51–83 (1978). 21. E. Olson, “On computing the average orientation of vectors and lines,” in “Robotics and Automation (ICRA)”, Proc. IEEE, pp. 3869–3874, (2011). 22. C. Agostinelli and U. Lund, R package circular: Circular Statistics (version 0.4-3), CA: Department of Environmental Sciences, Informatics and Statistics, Ca’ Foscari University, Venice, Italy. UL: Department of Statistics, California Polytechnic State University, San Luis Obispo, California, USA (2011).

1.

Introduction

Digital holographic microscopy (DHM) is a powerful, non-invasive light-microscopy technique that enables retrieval of not only the amplitude of light scattered by the observed sample, but also of its phase. The detected phase shift has physical interpretation: in reflected-light DHM (RDHM), it is proportional to the displacement of the observed surface [3], and in transmittedlight DHM (TDHM) it is, to a large extent, proportional to dry-mass density distribution of observed cells [4]. Apart from having this significance, it provides us with more complete knowledge of the light wave scattered by a sample, thus allowing application of advanced methods of image processing. These methods can, for instance, compensate optical aberrations. The aberrations compensation through inverse filtering has been successfully applied to highly aberrated images of samples acquired using objectives with high numerical aperture (NA) [5]. Another of these methods is numerical refocusing (see section 2 for the definition). Numerical refocusing of reconstructed complex amplitudes acquired using a quasi-point, monochromatic (coherent) source has been demonstrated [2]. Additionally, limits for refocusing were derived for cases when images acquired with an extended-source are refocused using formulas that were derived for point-source illumination. According to [2] (after transformation of variables), the refocusing is accurate if the refocusing distance ∆z satisfies p λ 1 − NAi 2 |∆z|  , (1) NA NAi where NA is the numerical aperture of the objective and NAi is the effective numerical aperture of illumination — NAi = n sin ϕmax , where n is the index of refraction of immersion and ϕmax is the greatest angle of incidence of the illumination rays in the object space. The benefit of the extended-source illumination is twofold: it removes high-coherence effects which give rise to speckles that degrade the image [6], [7] and it also increases the transversal resolution [8]. The image acquired using extended-source illumination has its frequency spectrum (see (3) for the definition and (5) for the meaning) filtered, which also affects numerical refocusing. The function that expresses the attenuation of the image frequency spectrum has been derived in [9] for the case when the source is assumed to be partially spatially coherent and NAi  NA: O(Qt ) ∝ S (c |∆z| Qt ) , (2)

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28259

where c is a constant depending on the experimental setup, Qt is a particular spatial frequency vector and S is the Fourier transform of the source intensity distribution. In the referenced case, the source was a laser beam focused onto a rotating ground glass, so the intensity distribution was assumed to have the Gaussian shape. The refocusing distance then ranges from tens to hundreds of micrometers for objectives with NA = 0.25. Another limitation of numerical refocusing follows from the discrete nature of the data [2]. This limitation actually relates to whether the field of view is sufficiently large in the sense that surroundings of all points we want to refocus correctly is recorded, too. Typically, this limitation is not significant and it can be mitigated by using a border processing technique (see eg [10]). Advanced applications of numerical refocusing are known and include focal plane detection [11] and 3D motion tracking of moving objects [12]. Numerical refocusing can be a step in the complex amplitude reconstruction process [7]. We have derived correct formulas (propagators) for numerical refocusing for the case of spatially incoherent large extended sources (ie NAi > NA/2) that are based on 3D coherent transfer function of the system. These propagators significantly differ from those for coherent sources. Our arrangement uses extended-source primarily to increase both the transversal and longitudinal resolution, the speckle reduction is regarded as a side-effect. Inseparable from that, the presence of out-of-focus sample features is suppressed even further, so refocusing distances are in our case lower. We have experimentally proved that when using a large extended-source, propagators based on coherent source produce wrong results and our derived propagators produce accurate results. We present a comparison of the two propagators and evaluate their accuracy using statistical methods. 2.

Imaging in DHM

Let w(x, y, z) be a value of the reconstructed complex amplitude [13] at point [x, y, z] that we can express alternatively as wt (qt ; z), qt = (x, y). The image point is determined by the x, y coordinates of the optically conjugated point of the object plane [14]. The coordinate z refers to the axial position of the sample relative to the object plane. Numerical refocusing of reconstructed complex amplitudes is a process of recovering wt (qt ; z0 ) from wt (qt ; z0 + ∆z), where we call ∆z the refocusing distance. We define the frequency spectrum Wt (Qt ; z) of the acquired 2D reconstructed complex amplitude wt (qt ; z) in the form Wt (Qt ; z) =

ZZ R2

wt (qt ; z) exp (−2πi qt · Qt ) d2 qt ,

(3)

ie Wt is the Fourier transform of wt along x and y for a given z. Q = (X, Y, Z) denotes a spatial frequency vector and Qt = (X, Y ) then analogously denotes a transverse spatial frequency vector. Linear systems theory can be used to describe the image formation in DHM [14] if we restrict ourselves to the 1st Born approximation of the scattering theory, if we assume that the DHM consists of aplanatic imaging system and if we assume that the extended-source is spatially incoherent. Then, according to [14], the reconstructed complex amplitude w(q) = wt (qt ; z) can be expressed [1] as w(q) =

ZZZ R3

H(Q)T (Q) exp(2πi q · Q) d3 Q,

(4)

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28260

ie the inverse Fourier transform of H ×T , where H is the system’s 3D coherent transfer function (CTF) and T is the 3D Fourier transform of the sample’s 3D scattering potential t(q) [14]. Means of calculating the H have been described in [14] for sources of arbitrary coherence state. Cross-sections of supports of CTFs for monochromatic point-source illumination and monochromatic large extended-source illumination are shown in Figs. 1(a) and 1(b) respectively. We assume rotational symmetry of H with respect to the optical axis meaning that Qt 1 = Qt 2 ⇒ H(Qt 1 , Z) = H(Qt 2 , Z). We restrict ourselves to planar samples — typically reflective surfaces observed in RDHM, whose scattering potential t can be expressed [1] as t(q) = tt (qt ) δ (z), where δ refers to the Dirac distribution. Then we have T (Qt , Z) = Tt (Qt ) for all Z, where Tt is the 2D Fourier transform of tt . If we combine (4) with (3) while using the fact that T is constant with respect to its Z-axis, we obtain simpler expression: Wt (Qt ; z) = Tt (Qt ) Ht (Qt ; z),

(5)

where Ht (Qt ; z) can be regarded as a 2D CTF for a planar sample defocused by z: Ht (Qt ; z) =

Z

H(Q) exp (2πi z Z) dZ.

(6)

R

Z 0

b) Qt

HT

Qt

Qt =

NA λ

HT

Z 0

− λ1

− λ1

HR HR

Qt = 2 NA λ

a)

− λ2 + BL

− λ2

− λ2

Fig. 1. Cross-sections of support of CTFs for an objective with NA = 0.5 in the reflectedlight (HR ) and transmitted-light mode (HT ). Calculated for (a) point-source illumination and (b) extended-source illumination with NAi = NA, BL denotes the longitudinal frequency bandwidth (see section 4.2). CTFs are rotationally symmetric with respect to the Z axis.

3.

Refocusing theory

Numerical refocusing for planar samples from z to z+∆z can be expressed in terms of frequency spectra using (5) for DHM setups that satisfy (4): Wt (Qt ; z + ∆z) = Wt (Qt ; z) P(Qt ; ∆z; z),

(7)

where we call the generally complex-valued function P propagator and

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28261

P(Qt ; ∆z; z) =

Ht (Qt ; z + ∆z) , Ht (Qt ; z)

for all Qt so that Ht (Qt ; z) 6= 0.

(8)

We will discuss cases when Ht (Qt ; z) equals zero (or it is very low) in the section 3.2. We use the notation P(Qt ; ∆z) further in the text for cases when z = 0 (this corresponds to the case when recovery of focused sample is in question) and the function H the propagator is based upon clearly follows from the context. We can compute the 3D CTF according to [14], then the corresponding 2D CTFs according to (6) and the propagator according to (8). As a result, we obtain an invertible model of imaging for aberration-free imaging systems with extended non-monochromatic sources. By expressing the complex function P in a polar form as P = |P| exp(i arg P), we denote PA = |P| the amplitude and Pϕ = arg P the phase of the propagator P. Then filtering the image frequency spectrum Wt using only Pϕ is called phase deconvolution and the use of both the amplitude and the phase is referred to as to complex deconvolution [5]. Further in the text, we are going to focus only on imaging in RDHM, mainly because numerical refocusing of planar samples is more applicable in this type of microscopy. 3.1.

Coherent source

When illumination by an axial monochromatic point-source is used and observation of reflective surfaces is assumed, the corresponding RDHM 3D CTF can be expressed in a closed-form: 1

H(Q) = Po (Qt ) δ [K + Z + (K 2 − Q2t ) 2 ] for Qt ≤ NA/λ ,

(9)

where K = n/λ is the reduced wave vector in the object space and λ is the vacuum wavelength. The exact form of the term Po weakly depends on the sample type and can be looked up in [14]. This CTF support is shown in Fig. 1(a). The frequency transmission |Ht | corresponding to this CTF does not change with ∆z. Also note that it is virtually constant (except the cut-off) for low numerical apertures and reflective conductive surfaces, as it is shown in Fig. 2 (dashed line). Then, after calculating the respective Ht (Qt ; ∆z) and using it in (8), we obtain the appropriate propagator: 1 P(Qt ; ∆z) = exp[2πi∆z(K 2 − Q2t ) 2 ] for Qt ≤ NA/λ . (10) 1

We extend the definition of the propagator to exp[2πi∆z(K 2 − Q2t ) 2 ] for all spatial frequencies, because it is applied to Wt , which is equal to zero for Qt > NA/λ nevertheless. If we consider the inversion of the imaging process based on (5) and denote the result (ie the recovered frequency spectrum of the sample) by Tf∗ (Qt ; ∆z), we get Tf∗ (Qt ; ∆z) =

Wt (Qt ; ∆z) P(Qt ; ∆z)

(11)

for all Qt and ∆z. This is another expression of the fact that phase numerical refocusing produces an accurate image of the sample in case of point-source illumination. Finally, if we take the 2D inverse Fourier transform of Tf∗ (Qt ; ∆z), we obtain the corresponding sample’s scattering potential tf∗ (qt ; ∆z). 3.2.

Extended source

In this section, we focus on the extended monochromatic source, when generally neither the 3D CTF H(Q) nor its corresponding Ht (Qt ; z) can be expressed in a closed-form. Considering the #194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28262

phase deconvolution, Pϕ (Qt ; ∆z) is very sensitive to the size of the source. Behavior of some phase propagators that were obtained by numerical evaluation of their corresponding CTFs is shown in Fig. 3. Notice that in the TDHM case (Fig. 3(a)), the propagator’s phase variation decreases with the growing source size. In other words, phase refocusing is not relevant in TDHM with large extended-source (ie when NAi → NA). 2.0

NAi → 0

NAi = 0.24

any z

z = 0 µm z = 8 µm z = 10 µm

|Ht (0, Qt ; z)|

1.5 1.0 0.5 0.0

−2

−1

0 Qt / NA λ

1

2

Fig. 2. Dependence of |Ht (0, Qt ; z)| (ie transmission of spatial frequencies) on defocus. Computed for RDHM, NA = NAi = 0.25 (full lines) and NAi = 0.05 (dashed lines). Narrow-band filter with central λ = 550 nm, full width in the half maximum (FWHM) = 10 nm.

We extend the definition of P to all spatial frequencies analogously to the point-source case. This time, |Ht | is not constant on its support as it can be observed in Fig. 2. In order to prevent spurious amplification of spatial frequencies in the frequency spectrum we deconvolve, we introduce measures similar to ones introduced in [16] or [5]. We introduce the operator G{·}: G{ f ; c} = exp(i arg f )

|f|+c , 1+c

c ∈ [0, ∞).

(12)

Selecting c = 0 yields G{Ht ; 0} = Ht , conversely, we define G{Ht ; ∞} = exp(i arg Ht ). As it is shown in Fig. 2, transmission for defocus z = 10 µm equals zero both for high and low spatial frequencies. The operator G{·} prevents amplification of noise equally well when filtering the focused signal, where only the highest spatial frequencies are problematic, and the mentioned signal defocused by z = 10 µm, when possible division by zero becomes an issue for low spatial frequencies, too. The constant c depends on the noise, however the acquisition noise is not the main limiting factor. Artifacts present along the optical path cause disturbances that appear in reconstructed complex amplitudes. They are not band-limited and they are the main source of the image degradation when applying the deconvolution and they determine the minimum value of c. Given the defocus z, refocusing distance ∆z the system’s CTF H and the constant c, we define the propagator as Ht (Qt ; ∆z) . (13) P(Qt ; ∆z) = G{Ht (Qt ; 0); c}

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28263

0.0

3.0

arg Ht (0, Qt ; 8 µm)

a −0.5

2.0

−1.0

1.5

−1.5

NAi

−2.0

0.24 0.12 →0

−2.5 −3.0 0.0

b

2.5

0.5

1.0 1.5 NA Qt / λ

2.0

1.0

NAi 0.24 0.12 →0

0.5 0.0 −0.5 0.0

0.5

1.0 1.5 NA Qt / λ

2.0

Fig. 3. Phase of the 2D CTF for defocus of 8 µm) for TDHM (a) and RDHM (b). Monochromatic extended-source illumination for selected values of NAi . NA = 0.25, λ = 550 nm. Shifted so that arg Ht = 0 for Qt = 0.

We now introduce the function T ∗ (Qt ; ∆z) presenting the inversion of the imaging process of complex amplitudes defocused by ∆z for extended-source: T ∗ (Qt ; ∆z; c1 ) =

Wt (Qt ; ∆z) , G{Ht ; c1 max |Ht |}

(14)

where Ht ≡ Ht (Qt ; ∆z). The recovered sample’s scattering potential t ∗ (qt ; ∆z; c1 ) is again the inverse 2D Fourier transform of T ∗ (Qt ; ∆z; c1 ). Since further in the text, we use either the complex deconvolution (14) with c1 = 0.4 or the phase deconvolution, in order to keep the notation simple, we define tp∗ (qt ; z) = t ∗ (qt ; z; ∞) and tc∗ (qt ; z) = t ∗ (qt ; z; 0.4). 4.

Experiment

From now on, we use the term amplitude image to refer to the absolute value of reconstructed complex amplitude and we use intensity image to refer to the classical microscope output signal. 4.1.

Sample and methods

Experiments were conducted on a reflective amplitude planar sample — a golden film with small pits on SiO2 (Fig. 4). The optical setup referred to as RDHM is the one described in [17]. No immersion was used for any observation, so the index of refraction n = 1. Two fundamentally different devices were used for the following three series of observations marked as W, P and E: W Single reference observation was done on Nikon Eclipse L-150 classical reflected-light widefield microscope in white light with a 100 × / NA = NAi = 0.95 objective. P z-stack (2 µm step) using the RDHM with a 10 × / NA = 0.25 objective, central λ = 550 nm, filter FWHM 10 nm using a quasi-point-source. Value of NAi = 0.05 was determined by analysis of transmitted frequencies in the frequency spectra. #194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28264

0.75 0.50 0.25 20µm

|wt (qt ; 0)|/a.u.

1.00

0.00

Fig. 4. Visualized |wt (qt ; 0)| of the field of view recorded by RDHM in the extended-source imaging setup (setup E, see the list in section 4.1). Edges of the field of view were apodized (see section 4.4). The stroke in the upper right corner is a diffraction grating artifact, the white rectangle points out the section of the sample shown further in the paper in more detail.

E z-stack 1 µm step) using the RDHM with a 10 × / NA = 0.25 objective, central λ = 550 nm, filter FWHM 10 nm using an extended-source. The source image has filled condenser pupil so the theoretical NAi = 0.25. Analysis of transmitted frequencies in the frequency spectrum of the defocused complex amplitude has revealed the effective NAi we use for the processing is lower, NAi = 0.15. Since the non-monochromatic RDHM with an extended-source has the optical sectioning property [18], the refocusing length has an upper bound — the reconstructed complex amplitude’s frequency spectrum gets progressively attenuated (see the relative transmission decrease in Fig. 2). We have found out that we can acquire data with acceptable signal-to-noise ratio for |∆z| < 8 µm, so we analyze numerical refocusing for two extreme values of defocus ±8 µm for both quasi-point-source and large extended-source illumination. 4.2.

Sampling

Since the complex amplitude is a band-limited signal with the highest transmitted transverse spatial frequency equal to (NA + NAi )/λ ≤ 2 NA/λ (see Fig. 1), it can be sampled. The minimal sampling distance is the inverse of the signal’s bandwidth B [19]. Given a coordinate axis and the length of field of view along that axis D, N/B = D, where N is the number of samples. We can express N = B×D. The transversal frequency bandwidth of the reconstructed complex amplitude is twice the highest transmitted transversal spatial frequency, the longitudinal bandwidth BL is shown in Fig. 1(b). So for the number of transversal and longitudinal samples, Nt and NL respectively, we have p NA 1 − 1 − NA2 Dt , NL = 2 DL , (15) Nt = 4 λ λ where Dt refers to the length or width and DL to the height of the area in which we want to be able to refocus. The conclusion is that in our experimental setup, the sufficiently sampled #194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28265

a

4

b

Y / NA λ

1

2

0

0

−1 −2

Pϕ /rad

2

−2 −2

−1

0 X/ NA λ

1

2

−2

−1

0 X/ NA λ

1

2

Fig. 5. The unwrapped phase difference Pϕ (Qt ; 8 µm) between Wt (Qt ; 0) and Wt (Qt ; 8 µm). Both phase shifts are unwrapped in such way that the zero frequency has phase equal to 0 in both cases. (a) extended-source, setup E and (b) quasi-point-source P.

CTF 3D array contains 268 × 200 × 30 elements. Therefore it is feasible to perform calculations involving it on ordinary office PCs. Concerning the sample itself, spans beyond the lower edge of the field of view (see Fig. 4), so we have used the Hanning function [20] to apodize it. This measure suppresses spurious spatial frequencies in the sample’s frequency spectrum and its success is apparent in Fig. 5, since the visualized phase difference is pure noise for any frequency Qt > 2NA/λ . 4.3.

Observations

Observations of the section of the in-focus sample (see Fig. 4) are shown in Fig. 6. The black circle • at the bottom left corner of images has radius of 1.22 λ /(NA + NAi ), which corresponds to the Rayleigh resolution criterion in the sense that images of two points where one is located at the center of the circle and the other one outside of it are considered resolved [15]. Since the reconstructed complex amplitude wt alone is not an accurate representation of the sample because of the different transmission of spatial frequencies in the case of extendedsource (shown in Fig. 2), corresponding tc∗ was calculated using (14). We can see the difference between the unprocessed |wt (qt ; 0)| and |tc∗ (qt ; 0)| in Figs. 6(c) and 6(b) respectively and Figs. 7(c) and 7(b) respectively. Visual comparison of Fig. 7 (a)–(c) entitles us to say that the imaging inversion (14) produces a higher-contrast and more accurate representation of the sample without adding any kind of artifacts. Figure 6(d) shows the quasi-point-source amplitude image, that has clearly lower resolution than the extended-source one. The difference between Figs. 6(d) and 6(c) reveals the increase of resolution due to the extended-source illumination. Examining Wt of the two observations reveals that the frequency spectrum of the extended-source-illuminated amplitude image contains transmitted frequencies at least 1.6× higher than the other one (the theoretical increase is 2×, compare the difference in Fig. 5). There are two main reasons that can explain why the observed resolution increase 1.6× is lower than the one predicted by the theory 2×: • Higher spatial frequencies transmitted from the sample may be present in the frequency spectra, but they are overpowered by noise due to the fact that their transmission is pro#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28266

1.00 a

1.00 b

0.81

0.67

0.63

0.33

0.44

0.00

1.00 c

1.00 d

0.74

0.71

0.48

0.41

10µm 0.21

0.12

Fig. 6. Comparison of in-focus sample images (see the full sample in Fig. 4): (a) reference intensity image (setup W), (b) RDHM (setup E) filtered extended-source amplitude image |tc∗ (qt ; 0)|, (c) RDHM (setup E) extended-source amplitude image |wt (qt ; 0)| and (d) RDHM (setup P) point-source amplitude image |wt (qt ; 0)|. The black circle has radius of 1.22 (NA + NAi )/λ , see section 4.3. The colormap is linear and adjusted separately for each image to provide the highest contrast possible.

gressively reduced — see Fig. 2. This applies particularly to the defocused complex amplitude, where the signal-to-noise ratio is essentially lower. • Certain degree of illumination misalignment is inevitable, rendering the effective source size is lower than the nominal one as it is described in [13]. 4.4.

Propagator measurement

Having a z-stack allows us to measure the propagator introduced in (8) directly using the frequency spectra and relate it to the model. Specifically, assuming |Wt (Qt ; 0)| 6= 0, Pϕ (Qt ; ∆z) = arg

Wt (Qt ; ∆z) . Wt (Qt ; 0)

(16)

The experimental data are discrete and affected by noise, so we introduce a mapping that maps function of two discrete variables to a function of one (radial) continuous variable. This facilitates analysis of data that are assumed to be rotationally symmetric, such as Pϕ (Qt ; ∆z). The mapping is based on taking averages of the examined function over circles. Having a discrete function of two-dimensional variable such as the phase difference Pϕ (Qt ) (we omit the parameter ∆z now for simplicity), we focus on its domain, the space of spatial frequencies. We choose a set of spatial frequencies that lie on a ring of a given central radius Qt and width equal to the distance between adjacent spatial frequencies. We denote this set as S(Qt ). Finally, we define the set of all values of the function Pϕ (Qt ) on the ring as F(Pϕ , Qt ) = Pϕ (Qt ) for all Qt ∈ S(Qt ) Since we assume the function Pϕ (Qt ) to be rotationally symmetric, we could expect elements of F(Pϕ , Qt ) to be all the same. However, due to the noise or aberrations, it is not the case, so we define the function Pϕ 0 (Qt ) of the radial variable Qt as the average of values of F(Pϕ , Qt ).

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28267

1.0 a

b

c

0.8 0.6 0.4 0.2

5µm

0.0 1.0 |w|/ max |w|

0.8 0.6 0.4 0.2 0.0 reference (a) 30

31

32

filtered (b) 33

34

raw (c) 35

36

37

x/µm

Fig. 7. Comparison of close-ups corresponding to Fig. 6 (center of the upper edge): (a) reference intensity image (setup W), (b) RDHM (setup E) filtered extended-source amplitude image |tc∗ (qt ; 0)|, (c) RDHM (setup E) extended-source (raw) amplitude image |wt (qt ; 0)|. The colormap is common for all three images. The cross-section at the bottom is taken along the path marked by the white bar in individual images above.

As it was explained below (12), artifacts present along the optical path disturb the frequency spectrum most. In our case, those disturbances were caused by dust particles present on a glass just in front of the camera’s CCD chip. Such disturbances are not band-limited in the way the signal is and they are the main reason we have to set c1 in (14) as high as 0.4. Setting the right value of c1 is based on preventing accentuation of these disturbances, so they don’t distort their surroundings. Unfortunately, this type of noise can’t be easily separated neither from the image nor from the spectrum. The experimentally measured PA (Qt ; ∆z) = |Wt (Qt ; ∆z)/Wt (Qt ; 0)| differs significantly from the theoretical PA = |Ht (Qt ; ∆z)/Ht (Qt ; 0)|. As a result, we don’t use complex deconvolution (that involves it) for other purposes than to filter in-focus complex amplitudes. Analysis of the spectrum was not conducted on the part of the sample shown for example in Fig. 6, but on the whole sample — see Fig. 4. Figure 5 shows the measured Pϕ (Qt ; 8 µm) — ie the change of phase in the frequency spectra corresponding to the defocus of 8 µm for the two illumination source sizes. The measured Pϕ (Qt ; 8 µm) for extended-source is less symmetric, but the area where it is well-defined (ie not noisy) spans beyond NA/λ . Conversely, in the case of coherent illumination, Pϕ (Qt ; ∆z) is well-defined only for Qt < NA/λ . The experimental Pϕ 0 (Qt ; 8 µm) corresponding to Fig. 5 is visualized in Fig. 8(a) (solid lines). The experimental data are noisy, therefore 95% confidence intervals for the mean value of Pϕ 0 (Qt ) were added. Since the examined variable is phase, a statistical model that can work with directional data had to be used. Therefore we use statistical tools based on von Mises #194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28268

4.5

4.0 a

b

3.0 2.5

3.0 2.5 2.0 NAi = 0.15 NAi = 0.05

2.0

NAi = 0.15 NAi = 0.05

Pϕ 0 (Qt )/rad

3.5

3.5

Pϕ 0 (Qt )/rad

4.0

1.5

1.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 NA Qt / λ Qt / NA λ Fig. 8. Phase difference Pϕ 0 (Qt ) corresponding to defocus (a) 8 µm (see also Fig. 5) and (b) −8 µm. Solid lines correspond to experimental data, dashed lines to the theory and the bands show 95% confidence interval for the experimental data. The experimental data has been smoothed out to reduce the quantization error by applying convolution with a Gaussian kernel of σ = 0.05 NA/λ .

statistical distribution, which itself can be described as being very close to wrapped Gaussian distribution [21]. The 95% confidence intervals were obtained by fitting the von Mises distribution using [22] on the sample of values F(Pϕ , Qt ) for each given Qt . The number of values in the data sample |F(Pϕ , Q)| > 50 for Qt > 0.16 NA/λ . In other words, the data sample for any spatial frequency greater than the 0.08 multiple of the theoretical highest transmitted spatial frequency has more than 50 elements and therefore is large enough to draw conclusions from it. The theoretical Pϕ 0 (Qt ) is visualized in Fig. 8 (dashed lines). It has been computed using numerically evaluated CTF and it is described in more detail in section 3.2. Figure 8 shows that the theory is most inaccurate for low spatial frequencies. Otherwise the theoretically computed Pϕ 0 mostly falls in the confidence interval of the mean of experimental data. This indicates numerical refocusing using the theoretical values should yield good results. 4.5.

Refocusing

It is clear that there are features on the sample that are imaged in utterly different way between the classical and holographic microscopes (such as the features by the center of the upper edge in Fig. 6 or the feature above the cross-section line in Fig. 7 that yield much higher contrast in DHM images). This is caused by the fact that opposed to classical microscopes, the DHM imaging process itself is coherent even if an incoherent illumination source is used [14]. Consequently, we are going to compare cross-sections of data acquired only by RDHM. We have applied the filtering (14) to obtain |tc∗ (qt ; 0)| from in-focus complex amplitudes first. Then we have taken the defocused complex amplitude wt (qt ; 8 µm) from the z-stack, and computed the properly refocused |tp∗ (qt ; 8 µm)| and the coherent approximation |tf∗ (qt ; 8 µm)|. Then we compare the corresponding RDHM reference |tc∗ (qt ; 0)| with |tp∗ (qt ; 8 µm)|, |tf∗ (qt ; 8 µm)| and the originally acquired defocused amplitude image |wt (qt ; 8 µm)|. If we apply the refocusing formula (1) to our quasi-point-source setup P, we obtain the

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28269

b

c

d

|w|/ max |w|

a

reference (a)

refocused (c)

recorded (b)

coherent (d)

10µm 34.0 34.5 35.0 35.5 36.0 36.5 37.0 37.5 x/µm Fig. 9. Refocusing in the setup P (quasi-monochromatic, quasi-point-source) amplitude images: (a) in-focus deconvolved |tc∗ (qt ; 0)|, (b) defocused by 8 µm |wt (qt ; 8 µm)|, (c) phase refocused |tp∗ (qt ; 8 µm)| using the correct propagator and (d) phase refocused |tf∗ (qt ; 8 µm)| using the approximative (coherent) propagator. See Fig. 4 for the colorbar. The cross-section on the right is took along the path marked by the white bar in individual images on the left.

expression 8  40 that indicates that coherent propagation should yield acceptable results. Actual results presented in Fig. 9 confirm this. While the top view in Fig. 9 illustrates that it is possible recover the focused amplitude image of the sample even if the defocused amplitude image is severely distorted, the cross-section confirms that the recovered function is accurate in the sense that it differs very little from the reference (see Fig. 9(c) and 9(d)). Results for the extended-source are shown in Fig. 10. If we apply the refocusing formula (1) the setup E as we did in the previous case, we obtain the expression 8  14, which is not likely to hold. Figure 10(d) illustrates that phase refocusing using the coherent propagator expectedly produces wrong results, while the correctly derived phase propagator is able to invert the imaging process accurately. Although the difference between the focused and defocused images is not as apparent as in the previous case, now there are much finer details present and also the contrast is higher (compare Fig. 9 with Fig. 10). Consequently, the cross-section plot in Fig. 10 (right) becomes much more important because it shows that that the pit the cross-section traverses is not localized very well in the defocused amplitude image, while it is accurately recovered by the refocusing process. 5.

Conclusion

We have described the problem of numerical refocusing for large (NAi > NA/2) extended, spatially incoherent light sources, introducing the image inversion concept that is specific to the case of large extended sources. We have described the refocusing theory in terms of spatial frequencies spectra and developed methods how to verify the theory on experimental data. We have conducted a refocusing experiment with two quasi-monochromatic RDHM arrangements — one with quasi-point-source and one with large extended-source. The processed data have shown that the theory corresponds to the experiment in terms of phase of spatial frequency spectra (see Fig. 8) as well as in terms of complex amplitudes. Specifically, we have confirmed that the quasi-point-source setup can be accurately approximated by a point-source monochromatic setup (see Fig. 9) and the coherent propagation can #194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28270

b

c

d

|w|/ max |w|

a

reference (a)

refocused (c)

recorded (b)

coherent (d)

10µm 34.0 34.5 35.0 35.5 36.0 36.5 37.0 37.5 x/µm Fig. 10. Refocusing in the setup E (quasi-monochromatic, extended-source) amplitude images: (a) in-focus deconvolved |tc∗ (qt ; 0)|, (b) defocused by 8 µm |wt (qt ; 8 µm)|, (c) phase refocused |tp∗ (qt ; 8 µm)| using the correct propagator and (d) phase refocused |tf∗ (qt ; 8 µm)| using inappropriate (coherent) propagator. See Fig. 4 for the colorbar. The cross-section on the right is took along the path marked by the white bar in individual images on the left.

be used for refocusing. We have demonstrated that this is not the case when the large extendedsource is used (see Fig. 10(d)). We have confirmed that when the condition (1) doesn’t hold, usage of the coherent propagation for refocusing is inadequate. We have shown that the phase numerical refocusing can yield accurate results (see Fig. 10(c)), although it is only an approximation of the complex description of the phenomenon, formally similar to the coherent propagation. Moreover, this simplified approach described in section 4.5 is numerically stable since it doesn’t amplify the amplitude of spatial frequency spectra, which is often affected by disturbances that can disrupt the imaging process inversion. The usage of extended source of illumination makes the setup more susceptible to various aberrations and the assumption of spatial incoherence of the extended-source may be incorrect. Despite these issues, we have demonstrated that it is possible to accurately recover the function t ∗ (qt ; 0) representing the sample using solely the phase deconvolution. The refocused images retain high spatial frequencies that are not present in images acquired using the quasi-pointsource. Acknowledgments The authors thank to colleagues from Experimental biophotonics group (CEITEC) for helpful discussions, especially to Martin Loˇsˇta´ k. This work was supported by Technology Agency of the Czech Republic (TE01020229), Ministry of Industry and Trade of the Czech Republic (FRTI4/660), CEITEC – Central European Institute of Technology (CZ.1.05/1.1.00/02.0068) from European Regional Development Fund and by the grant “Specific research” FSI-S-11-22 of Brno University of Technology.

#194878 - $15.00 USD Received 13 Aug 2013; revised 18 Oct 2013; accepted 23 Oct 2013; published 11 Nov 2013 (C) 2013 OSA 18 November 2013 | Vol. 21, No. 23 | DOI:10.1364/OE.21.028258 | OPTICS EXPRESS 28271

Numerical refocusing in digital holographic microscopy with extended-sources illumination.

Numerical refocusing can be seen as a method of compensating the defocus aberration based on deconvolution by inverse filtering [1] in digital hologra...
4MB Sizes 2 Downloads 0 Views