1. Introduction

In this article, we characterize the lateral field distortions of a low numerical aperture (NA) and large field-of-view (FOV) fluorescence imaging system. We used a commercial fluorescent zooming microscope (also known as MACROscope) where the FOV changes with the optical zoom as it is well suited to illustrate the applicability of our characterization procedure. This setup from Leica™ (Fig. 1(a)) is a macro documentation system, combined with fluorescence techniques for visualization of sample preparations at a range of zooms [1].

 

Fig. 1 (a) Schematic of a simple wide-field fluorescence MACROscope (Reproduced from [1]); Best of two worlds: maximum intensity projection along the optical axis of a Convallaria majalis sample taken using a Leica™ ZAPO16, fit with a confocal scanning head, at (b) a minimum zoom setting with lateral pixel size of 1.09 μm and (c) a sub-region of the sample at the maximum zoom setting with lateral pixel size of 0.89 μm (Courtesy of INRA). The scale bars are 100 μm in length.

Download Full Size | PDF

As in the case of a microscope, the emitted fluorescence from the sample is collected by an objective lens, but a fine/coarse focusing can be obtained by using a macro lens. The apochromatic macro lens is combined with the objective lens to image large fields (about 20mm diagonal diameter) and to provide larger working distances (about 97mm). Existing wide-field microscopes offer high resolutions but with limited FOV, while stereomicroscopes offer larger FOV but compromise on the resolution. The MACROscope offers higher FOV and good lateral resolution (for NA between 0.12–0.50, a lateral resolution of 1.65–0.39 μm respectively). Variations of this setup is also available from other commercial vendors like the Axio Zoom V16 from Carl Zeiss or the AZC2 from Nikon. The principal difference between these commercial adaptations is the range of zooms that they work at.

A sample of the plant Convallaria Majalis is used to highlight the MACROscope’s imaging capabilities, under two different settings: minimum and maximum zoom. The three-dimensional (3-D) image volumes are shown in Fig. 1(b) and 1(c) as maximum intensity projections (MIP) along the optical axis.

1.1. Context

Most of the commercial MACROscopes are guaranteed to be telecentric. In telecentric images,

However, in the observation of the Convallaria Majalis specimen in Fig. 1(b), we noticed that with shifts in the focal plane, the specimen expands or contracts. The cell walls ‘appear’ tilted and thicker than expected. This effect is best illustrated in the observed image of a Haemocytometer (Fig. 2).

 

Fig. 2 Haemocytometer grid used for illustrating and measuring the distortion in the field. The square area indicated in red is of size 1 mm², in green is 0.0625 mm², in yellow is 0.04 mm² and finally the smallest in blue is 0.0025 mm². Reproduced from Wikimedia Commons

Download Full Size | PDF

In such a slide, the grid dimensions are calibrated, so that the image distortions can be quantified. For example, the square area indicated in red is of size 1 mm², in green is 0.0625 mm², in yellow is 0.04 mm² and finally the smallest in blue is 0.0025 mm². This slide was illuminated from below and the image was captured, from above, by a Photometrics CoolSNAP HQ² cooled CCD camera (6.45 μm × 6.45 μm pixel). The radial pixel size for the 12.7× zoom is 390nm, while the axial slice width was fixed at 50 μm to capture the entire volume. We have shown in Fig. 3 a single lateral focal plane and the projection along the y-direction for the transmitted image volume. The dimensions of the displayed volume is 343×343×6200 μm.

 

Fig. 3 The focal plane of the observed transmitted volume of the Haemocytometer (top) and the maximum intensity projection along the y-direction (bottom). The object is imaged using a 2x/air PlanApo objective fit to a Leica™ MacroFluo™ APOZ16. The zoom for this acquisition was set at 12.7x, the lateral sampling at 390nm, the slice thickness at 50 μm and the scale bar length is 50 μm. The total size of the displayed volume is 343×343×6200 μm.

Download Full Size | PDF

From the acquired data (see Media 1), we noticed the following:

1.2. Motivation and outline

The fundamental motivation underlying the study and characterization of distortions in any optical system is to identify the cause, understand the limits in system usage, and take necessary precautions to avoid or actions to correct it. The questions that we wish to answer as a result of this study are:

1.3. Notations

Scalar variables used in this article are denoted by lowercase letters (x), vectors are denoted by the boldface lowercase letters (x), and the matrices by the boldface uppercase letters (X). As the images are discrete, their spatial support is Ω_s = {(x,y,z) : 0 ≤ x ≤ N_x – 1, 0 ≤ y ≤ N_y – 1, 0 ≤ z ≤ N_z – 1}. By 𝒪(Ω_s) = {o = (o_xyz) : Ω_s ⊂ 𝕅³ ↦ 𝕉}, we refer to the possible observable objects, and we assign the function h : Ω_s ↦ 𝕉 to the microscope PSF. The observed intensities are denoted by i(x) : x ∈ Ω_s (bounded and positive), and a 3-D convolution operation between two functions is denoted by ‘*’. When the same symbol is used as a superscript over a given function, as in h^*(x), it represents the Hermitian adjoint operation on h(x). For a complex function h_A : Ω_s ↦ 𝔺, by |h_A(x)| and ∠h_A(x), we refer to its magnitude and phase respectively. While, by ℜ(h_A(x)) and ℑ(h_A(x)), we refer to its real and imaginary components. By Pr(·), we denote the probability density function.

The objective lenses of a microscope are defined by their magnification (M), numerical aperture (NA), and the medium in between the lens and the cover slip. For example, a lens of 5x magnification, 0.5 NA, and air as medium between the lens and cover slip is written as ‘5x/0.5 air’.

As mentioned earlier, we present 3-D images by their 2-D maximum intensity projection (MIP) along the optical axis in the 2-D XY plane or along the y-direction in the 2-D XZ plane.

2. Sensing the back aperture field

It is necessary to understand the conditions where our imaging system performs optimally, if we wish to extract the best from it. Nearly diffraction-limited (or aberration free) performance can be obtained, when the sources of these distortions are isolated. Once isolated, the objective would be to restore telecentricity in the images, by using post-acquisition computational methods.

2.1. PSF and role of the phase

The effective NA of the combined objective-zoom system is usually ≪ 0.7, and we work under near paraxial conditions. The effect of polarization is neglected, and the incoherent PSF can be modeled by using the scalar diffraction model. From the Kirchhoff-Fraunhofer approximation [4], we can write the near-focus amplitude PSF, h_A(x, y, z), in terms of the inverse Fourier transform of the two-dimensional (2-D) exit pupil function, P(k_x, k_y, z;NA), at each defocus z as

We observed that smooth variations in the amplitude of the pupil in Eq. (2) do not strongly affect the final PSF, while phase variations, such as defocus or aberrations, can produce an entirely different PSF. We use this as the basis for distortion characterization. We thus rephrase the question raised in Sect. 1.2: ‘Can the amplitude or the phase of the field at the back aperture of the optical system be an indication of the source of the distortions?’.

2.2. A Bayesian perspective

Although the phase information is not directly measurable in an incoherent imaging setup, it can be retrieved by choosing the PSF model that best fits the given bead image. This problem is also known as sensor-less wavefront sensing.

As the problem is both nonlinear and non-convex, for successful phase estimation, we require some prior knowledge about the field at the back aperture that we wish to estimate. It was shown in [6] that the out-focus-highlights (OOFHs) contains ‘partial’ information of the pupil or the back focal plane. There are several literature works like [7–9] that have studied wavefront reconstruction for adaptive optics (AO) control. The AO methods are based on the idea of phase aberration compensation by adding deformable mirror or phase modulating element in the optical system. A review of the recent trends in AO is given in [10]. The amplitude and phase of the pupil function can also be measured by using a fiber-optic interferometer, as was done in [11].

Wavefront sensing could also be accomplished computationally, for example, by using the GS algorithm [3]. Wavefront sensing [12] by phase retrieval is the process of estimating the amplitude and the phase of a pupil function from the observed 3-D intensities of an imaged point source. In the expression for h(x) in Eq. (3), as the only unknown is the phase φ_a from the aberrations, the problem of phase retrieval is a question of estimating the aberrated phase from the observed intensities. This problem of phase retrieval is normally under-determined. However, as the phase that is to be estimated does not change with defocus, it can be estimated if images of point source at multiple defocus positions are available. The only requirement is that these sections are sufficiently far from the focus. As the distance from the central focal plane grows to infinity, the intensity approaches that at the back pupil plane. However, in practice, the measurement of defocused beads becomes increasingly difficult for larger defocusing due to the decaying fluorescence intensities. We remark that the distortions that we observe are mainly amplitude aberrations. That is, they do not generate variations in the optical path difference (or the phase) of the light and so is not an aberration in the strict sense.

If we consider Poissonian photon counting statistics [13], the observed bead image can be written as:

For the GS algorithm, there is no intrinsic smoothness term on the solution. To compare our approach with the GS algorithm, we drop the prior energy term in Eq. (5) (by assuming a uniform distribution). The amplitude PSF can be estimated by the maximum likelihood (ML) algorithm:

For the fixed-point iterative algorithm, the step size, τ, was chosen to be 0.6 in all our experiments. The iterations are continued until either the mean-squared error (MSE) between the phase estimates for two successive iterations is below a pre-defined threshold ε or a pre-defined maximum number of iterations is reached by the algorithm.

Algorithm 1:. Proposed Algorithm.

View Table

2.3. Gerchberg-Saxton algorithm as a special case

The GS algorithm [3] is a technique to estimate the field at the back focal plane of the obejctive by following a forward and inverse Fourier transforms of the observation. The fixed-point iterative algorithm and the GS algorithm are initialized in the same manner for the pupil function. After initialization, a suitable curvature is added to the phase of the complex pupil function to obtain the defocus adjusted complex pupil function, P(k_x, k_y, z), at every defocus z [15]. This is inverse Fourier transformed (cf. Eq. (3)) to get the corresponding amplitude PSF (h_A(x)) intensities at the different defocus planes. The magnitude of h_A(x) is assigned to the corresponding measured intensities (after background subtraction) at the different defocus planes. A Fourier transform of this modified h_A(x) gives the new estimate of the defocus-adjusted complex exit pupil function, P̂(k_x, k_y, z), at the different defocus positions of z. The resulting defocus-adjusted complex pupil functions are readjusted back to zero and averaged to get a new estimate of the complex exit pupil function P̂(k_x, k_y, z = 0). This process is repeated until the MSE criterion or the maximum iteration is reached. Some constraints are introduced during the iterative algorithm that can aid in the convergence of the algorithm. The progress of the GS algorithm is the same as the fixed-point algorithm except for Step 6 in Algorithm 1. We see that when τ = 1 in Eq. (12), then the factor $i(x)/\left({\left|{\widehat{h}}_{\text{A}}^{(n)}(\mathbf{x})\right|}^2+b(\mathbf{x})\right)$, at each iteration, performs the assigning operation of Step 6. This ratio also has a physical significance. It has the role of replacing the incorrect amplitude of h_A(x) by the correct experimentally obtained magnitude i(x).

3. Experiments

In the previous section, we derived an iterative algorithm for phase retrieval from intensity data. In order to answer our original question in Sect. 2.1, we apply this proposed iterative algorithm on observed images of microbeads. For the imaging experiments, we chose InSpeck™Green (505/515) fluorescent microspheres of size 2.5 μm from Molecular Probes^®. This size was chosen because the resolution of the system, in the axial direction, for certain zooms are worse than this size. In addition, for lenses with low NA, the diversity defocus planes are not affected by the size. We diluted a 1 μl of this suspension in 20 μl of distilled water and dried a drop of this on to a coverslip. These dried beads were then imaged using a Leica™ MacroFluo™ fit with a 5× planapochromatic HR objective and 16 zoom positions. The maximum NA of the optical system is 0.5. Of the several data acquired and processed, we have chosen two settings for the objective and the zoom. When the objective is 2 × /air PlanApo, the zoom position is at 4.6 × /air (radial sampling 421nm and axial sampling 3000nm). While for the 5 × /air PlanApo, the zoom position is at 1.6 × /air (radial sampling 998.3nm and axial sampling 1000nm). As the beads are distributed randomly, we chose five locations in the lateral plane for cropping the bead images. These locations are shown in Fig. 4 with the positions of the PSFs marked out. Note that in Fig. 4, the PSFs are shown by their MIPs along the y-direction but their positions are shown in the lateral plane. For analysis, we chose only the beads from positions close to the periphery that exhibited apparent intensity chopping.

 

Fig. 4 Empirical PSFs are shown at the different positions (denoted by a cross) in the lateral field of the lens. Here, the PSFs are shown as the MIP along the y-direction.

Download Full Size | PDF

3.1. NA for convex relaxation

Given the fact that the MACROscope works under a variable zoom, the NA of the optical mount is variable. The setup that we performed our experiments with was a prototype that is not pre-calibrated by the vendor. So the effective NA for all the zooms of the mount are not directly available. However, if we consider the light as a cone, the apex of the cone is at the central focal intensity plane and the base of the cone is the observed diffraction ring at a defocus plane of distance H away from the center. D is the diameter of the largest concentric ring of the base (cf. Fig. 5). For example, at a zoom of 9.2× with the 2× objective, the maximum radius of the diffractive ring pattern at a distance of H = 61 μm away from the center was measured to be 32.46 μm. The measured radius D/2 is related to the angle α and the defocus distance H by tan(α) = D/2H = 32.46/61 = 0.53. The maximum subtended semi-cone angle will be α = arctan(0.53) = 0.49 radians. Since n_i = 1, the effective NA can be calculated to be 0.47 which is closer to the manufacture specified NA of 0.5. If we consider another set of images taken under a zoom of 1.6×, D was measured to be 32.70 μm for a H of 71 μm. In this case the NA was calculated as 0.22. We use such ‘loosely’ calculated NA values to limit the frequency bandwidth in the pupil plane of Eq. (1).

 

Fig. 5 The schematic to measure the maximum object spread for constraining the iterative algorithm.

Download Full Size | PDF

3.2. Results

As the problem is under determined, to introduce diversity, four defocus sections (M = 4 in Algorithm 1) were chosen, symmetrically, above and below the central focal plane of the 1.6× image. These sections lie at a distance of about 2–5 times the Rayleigh length from the focal plane. This bead was cropped from the periphery of the field and these individual sections approximate the OOFH that was mentioned earlier. One of the defocus sections is shown in Fig. 6(a). The unwrapped phase of the pupil, that was retrieved, φ̂_a, is shown in Fig. 6(b), after about 32 cycles of the fixed-point algorithm, with τ = 0.6. We allowed the algorithm to continue, although the solution for the electric field amplitude converged between 12–15 cycles. In order to reduce noisy estimates, at each iteration, it was suggested in [15] to filter the estimate obtained from the GS algorithm by a Gaussian filter. To avoid such an ad hoc method we propose, as future work, that at each cycle of the fixed-point iterative algorithm, the field amplitude be regularized by a total variation functional [13]. For reproducibility of the experiments, the complete source code in Matlab™, and the data are provided here: http://bioimageanalysis.org/praveen/code.zip.

 

Fig. 6 (a) The first section of the observed intensity, with z = −57 μm, and the scale bar is 10 μm, and (b) retrieved unwrapped pupil phase, φ̂_a. The bead image was cropped from the right peripherary intensity image of Fig. 4 at a zoom 1.6× (radial sampling of 998.3nm and axial sampling of 1000nm). τ = 0.6, the maximum number of iteration is 32 and the phase scale is between [−π, +π] radians.

Download Full Size | PDF

3.3. Discussion

From the estimated phase in Fig. 6(b), we see that the pupil function of the optical system is partially chopped. The chopping of the pupil is such that it resembles a ‘cat’s’ eye [6]. This could be the result of two limiting apertures (from the sizes of lenses in the objective and the zoom) creating the vignetting effect in peripheral regions of the field. In the schematic that is shown in Fig. 7, we illustrate such a distortion. The axial aperture is the complete circle while the oblique aperture is vignetted. Our reconstruction of the cat’s eye in Fig. 6(b) is not sharp. This could be explained as follows. Every lens creates a limiting pupil due to its physical dimension. The imaging of these pupils through the optical system on to the back aperture plane will make diffused circles when it is far from the conjugated back aperture planes.

 

Fig. 7 A schematic showing the effect of two limiting apertures (here zoom and objective lenses) at the back focal plane of the optical system. Here the on-axis and off-axis positions are shown.

Download Full Size | PDF

The retrieval of the phase allows us to understand the physics behind the field distortions. Although, some of these setups are claimed to be telecentric, we found that the system is not telecentric in the entire zoom range but only for high zooms. In addition, based on the amplitude of the retrieved electric field, we can also see the extent of overlap in the apertures. For example, from the defocus intensities in Fig. 8, the algorithm was able to retrieve the back aperture amplitude as shown in Fig. 9(b). Although the effective NA for this particular acquisition was calculated to be 0.17, the amplitudes were also retrieved with minor variations in the effective NA. We found that the estimation of the phase requires the exact NA value. In spite of erroneous input NA, the estimated amplitude on the other hand could still validate our hypothesis of vignetting (Fig. 9(a) and 9(c)). Even with an erroneous NA, the chopping could still be quantified to be between 84–89%.

 

Fig. 8 Diversity sections, i(x), taken at four symmetrical positions with defocus at (a) z = −36 μm, (b) z = −15 μm, (c) z = +15 μm and (d) z = +36 μm. The objective is a 2x/air PlanApo and the zoom is set at 4.6x. The slice width is fixed at 3 μm and the effective NA was calculated to be 0.17.

Download Full Size | PDF

 

Fig. 9 Retrieved back focal pupil amplitude, |P̂(k_x, k_y, z = 0)|, from the defocus sections in Fig. 8. The algorithm was run with variations in the effective NA (a) 0.05, (b) 0.17, (c) 0.20.

Download Full Size | PDF

The approach presented here can be thought of as adding a virtual wavefront or interferometric pupil plane sensor, and using it to retrieve the field information. By using fluorescence beads as ‘guide stars’, we can validate the causes of our distortions and quantify them. Although we have demonstrated our algorithm and experiments on a MACROscope, the field distortions studied here are present in any optical system working under low NA and low magnification conditions.

There are many ways to overcome the field distortions. Some of them are listed here:

In order to avoid the radial distortions that in turn can be produced as a result of such translations, the imaged object is moved in addition to the lens components. However, it is not clear what would be the defocus contribution of such an optical element translation to the final output image. It is likely that the imaging plane needs to be moved as well in addition to the optical elements.

Each of the above three methods have their advantages and difficulties. In the first case, the zoom system cannot be used at its full capacity while the other two are suggestions in redesigning the setup at the optical level. Our future work is therefore aimed at correcting these field distortions, computationally, after the images have been acquired.