Accurate unsupervised estimation of aberrations in digital holographic microscopy for improved quantitative reconstruction

Dylan Brault; Thomas Olivier; Ferréol Soulez; Sachin Joshi; Nicolas Faure; Corinne Fournier

doi:10.1364/OE.471638

1. Introduction

Optical microscopy can be used to extract several characteristics from a biological sample, such as morphological parameters, birefringence or a phase shifts introduced by an unstained sample. For quantitative measurement of these properties, an accurate optical model is required [1]. Accounting for the characteristics of the optical system is an essential component of reconstruction algorithms in optical microscopy. For example, in fluorescence microscopy, accurate modeling of the Point Spread Function (PSF) is a way to improve the deconvolution step [2–5]. It can be performed using either a dedicated calibration step (by directly measuring the PSF on “point-like” objects [2]) or by estimating the PSF directly on an image that presents aberrations [4,6,7]). In the literature, estimating aberrations or PSF have been widely addressed using various microscopy methods (fluorescence, single-molecule localization, wide-field microscopy, holography, etc.), with different measurement or reconstruction approaches and models of the PSF. These models can be very simple (e.g., Gaussian model), more realistic, like the Gibson-Lanni model [5,7,8], or more versatile and general, like the Zernike polynomials modelization of the pupil function [6,9,10]. In the two latter cases, the coherent PSF is modeled as a phase error function in the exit pupil plane of the objective.

In the particular case of digital holographic microscopy, the issues of aberrations estimation and correction have been widely studied for off-axis configuration (e.g., [11–14]). However, it concerns essentially the wavefront mismatch between the object and the reference beams, which creates distortions of the interference fringes, thus inducing errors in the reconstruction.

In-line digital holographic microscopy requires a simpler setup involving a single beam. It simply consists in recording the intensity pattern diffracted by a sample. It is less bulky and less sensitive to vibrations than off-axis holographic setups [15,16]. Image processing makes it possible to reconstruct the optical properties of the sample including its absorption and its phase shift. These can be discriminant in a classification or characterization task in the context of microbiology or medical diagnoses [17].

The aberrations of an in-line holographic optical system can have different causes, such as non standard uses of the objective, tilts or collimation errors in the illumination. These aberrations are dependent on the setup, its alignment and vary in the field of view. They lead to reconstruction errors, not only in the quantitative estimation of the modulus and the phase but also in the geometrical properties of the reconstructed objects. Thus, the repeatability as well as the reproductibility of the reconstructions is affected. However, the aberrations of the optical system are usually not considered in the reconstruction step. Accounting for the aberrations in the image formation model makes it possible to reduce the bias introduced in the reconstructions. These aberrations are an important issue to overcome in applications such as medical diagnoses that require reconstructions to be as accurate as possible to make the decision as robust as possible. To our knowledge, it is only recently that the influence of optical aberrations has been studied in the context of in-line digital holographic microscopy [1,18,19]. These studies underlined the need for a fine estimation of aberrations in order to improve the quantitativity and the repeatability of the phase reconstructions as well as the axial positioning, by reducing the aberration-driven biases.

In the present paper, we first address the problem of estimating aberrations in the context of in-line digital holographic microscopy. To that end, we use calibration beads to estimate an aberrated forward model. Using an Inverse Problems Approach (IPA), we simultaneously fit Zernike coefficients and calibration beads parameters, which are parameters of the forward model, on data. Unlike many PSF estimation studies, our approach does not require axial stacks of images, i.e., only one hologram is needed. Moreover, we made no assumption of an aberration-free PSF in the center of the field, like in Zheng’s et al. study [10]. Finally, this model of aberration is more general than the Gibson-Lanni model [5,7,18]. As a forward model, we use a Lorenz-Mie model of the calibration beads that has been extended to account for the aberrations of the optical system using Zernike polynomials [9]. To jointly estimate the calibration beads and aberration parameters, we choose a parametric IPA as it is known to be accurate in estimating the parameters of simple shape objects [20–22] and of the experimental parameters required for calibration. It has already been successfully applied in the context of autofocusing [23], for the estimation of the spectral crosstalk on a Bayer sensor [24] and to estimate the parameters of an astigmatic reference wave [25].

Once Zernike coefficients estimated locally for each bead, they can be used to perform aberration free reconstruction of the sample. These reconstructions can be performed using regularized IPA algorithm [26,27] or Fienup algorithm [28,29]. To test the proposed methodology, we use the experimental procedure of Martin et al. in [18], i.e., the use of a water immersion microscope objective with a correction collar that causes aberrations when not set correctly.

In the following section, we describe the method to estimate aberration parameters (Zernike coefficients) and use them to refine the PSF model of our holographic setup in order to reconstruct aberration-free images. In the third section, we detail the setup used to validate the proposed methodology. In the fourth section, to demonstrate the robustness of the approach to reconstruct various kinds of aberrations, we first present the estimation of both aberrations and beads parameters on simulated holograms and on experimental holograms. Finally, to illustrate the relevance of our approach on experimental data, the experimental data are reconstructed with phase retrieval algorithms (Fienup and regularized IPA algorithm) that take into account the estimated aberrations.

2. Estimation of the aberration parameters and reconstruction

Inverse problems are a general class of problems where unknowns are linked to measurements through a known image formation model (simulating the measurements is referred to as the “forward problem”). In this framework, reconstructions are based on minimizing the discrepancy between the hologram (the data) ${{{\boldsymbol {d}}}}$ and an image formation model (forward model) ${\boldsymbol {m}}$. In a general case, such phase retrieval problem is ill-posed as it has many degeneracies (more unknowns than data, twin image ,…). To solve it, it is necessary to inject some a priori on the solution into the minimization problem by adding regularization terms and/or constraints. Another way to solve this degeneracies is to use a model of the measurement that depends on only a few parameters. The problem can be then solved using the parametric IPA framework [30]. This framework is well suited to calibrate the aberrations using holograms of spherical objects as the image formation model depends only on the parameters of the objects (position, diameter and refractive index) and on the aberrations that can be modeled with a complex pupil function described by few parameters. Once these aberrations are estimated, they can be used in a regularized reconstruction method to reconstruct any sample without any aberration artefacts. Fig. 1 shows a flowchart representing the two main steps, the calibration and the reconstruction, that are detailed here after.

Fig. 1. Flowchart representing the two main steps of the proposed method: calibration and reconstruction.

Download Full Size | PDF

2.1 Calibration: aberration parameters estimation

The diffraction pattern $\underline {\boldsymbol {a}}^{\textrm {Mie}}$ of a spherical bead is accurately modeled by the Lorenz-Mie model [31] which depends on the set of bead parameters ${\boldsymbol {\vartheta }}$ $=\{x,y,z,r,n\}$, where ${x,y,z}$ corresponds to the 3D position, ${r}$ is the radius and $n$ is the refractive index (complex-valued variables are underlined). The Lorenz-Mie model has been successfully used to reconstruct spherical objects from holograms by fitting methods [20,22] or, in a more general framework, by parametric IPA [21,32,33]. In the presence of aberrations, the new image formation model of the diffraction pattern of the beads $\boldsymbol {m}^{\textrm {P}}$ ($\textrm {P}$ stands for Parametric) also depends on the aberration parameters $\boldsymbol {\alpha }$ of the optical system that can be included in the model by mean of a complex pupil plane as follows:

(1)$$\boldsymbol{m}^{\textrm{P}}({\boldsymbol{\vartheta}},\boldsymbol{\alpha})= \left| \mathscr{F}^{{-}1}\left[\widetilde{\underline{p}}(\boldsymbol{\alpha})\odot\widetilde{\underline{\boldsymbol{a}}}^{\textrm{Mie}}({\boldsymbol{\vartheta}}) \right] \right|^{2}$$

where $\mathscr{F}^{-1}$ is the inverse Fourier Transform, $\widetilde {\underline {p}}(\boldsymbol {\alpha })$ is the pupil function in Fourier domain that depends on $(\kappa _x,\kappa _y)$, the spatial frequency coordinates. For the sake of compactness, Fourier space coordinates and spatial coordinates are omitted in the equations when they are not required. $\boldsymbol {\alpha }=\left \{ \alpha _n^{m} \right \}_{(m,n) \in \mathbb {Z}^{2}}$ is a vector of aberration parameters, that will be referred as Zernike coefficients in this work. $\widetilde {\underline {\boldsymbol {a}}}^{\textrm {Mie}}$ is the Fourier Transform of $\underline {\boldsymbol {a}}^{\textrm {Mie}}$, and $\odot$ is the Hadamard product.

As described in [9,34], Zernike polynomials $\left \{\boldsymbol {Z}_n^{m} \right \}_{m,n}$ provide a suitable basis to describe the pupil function $\widetilde {\underline {p}}$ (see Appendix A for details):

(2)$$\widetilde{\underline{p}}(\kappa_x,\kappa_y,\boldsymbol{\alpha})=e^{i\left[\sum\limits_{n,m} \alpha_n^{m} \boldsymbol{Z}_n^{m}(\kappa_x,\kappa_y)\right]}$$

To characterize the aberration effects of the optical system, the Zernike coefficients $\boldsymbol {\alpha }$ have to be estimated. Assuming a white and Gaussian noise, the maximum likelihood estimation of model parameters $\{{\boldsymbol {\vartheta }}, \boldsymbol {\alpha }\}$ of the bead and the aberrations corresponds to a least squares fitting problem [20,32]:

(3)$$\left\{{\boldsymbol{\vartheta}}^{{\dagger}},\boldsymbol{\alpha}^{{\dagger}}\right\}=\mathop{\mathrm{argmin}}_{{\boldsymbol{\vartheta}} \in \mathbb{P}, \boldsymbol{\alpha} \in \mathbb{D}} \lVert {{{\boldsymbol{d}}}} -\boldsymbol{m}^{\textrm{P}}({\boldsymbol{\vartheta}},\boldsymbol{\alpha})\rVert ^{2} _2$$

where $\left \{\mathbb {P},\mathbb {D}\right \}$ are optimization constraints and $\lVert \cdot \rVert _2$ is the $L_2$-norm. Note, taking the weighted version of the $L_2$-norm (i.e., the squared Mahalanobis distance) makes possible to consider non-stationary and correlated noise into account [32,35,36].

To numerically solve this optimization problem (Eq. (3)), only the first 15 Zernike coefficients are estimated in the following. As a phase piston has no effect on the image formation model (intensity image formation model), $\alpha _0^{0}$ is set to $0$. As varying Zernike coefficients $\alpha _{1}^{-1}$ and $\alpha _{1}^{1}$ simply amounts to shift parameters $x$ and $y$, these Zernike coefficients are also set to $0$. In these conditions, seventeen parameters are studied:

x, y, z, r, n , \alpha^{{-}2}_2, \alpha^{0}_2, \alpha^{2}_2, \alpha^{{-}3}_3, \alpha^{{-}1}_3, \alpha^{1}_3, \alpha^{3}_3, \alpha^{{-}4}_4, \alpha^{{-}2}_4, \alpha^{0}_4, \alpha^{2}_4, \alpha^{4}_4.

A study of the correlations between the estimated parameters is presented in Appendix B. It shows some high correlations in the correlation matrix. All the parameters $\{{\boldsymbol {\vartheta }},\boldsymbol {\alpha }\}$ should therefore be estimated simultaneously. An iterative detection/local optimization scheme [21] is used to guarantee the rapid and accurate reconstruction of a set of objects. Since the beads are monodispersed, a narrow parameter research domain $\mathbb {P}$ can be chosen depending on the size and refractive index of the beads used experimentally.

Since the aberration can differ depending on the location of the beads in the field of view, the aberration parameters have to be estimated for several different bead locations.

2.2 Reconstruction: including aberration model

Once the aberrations are modeled, they are taken into account to better reconstruct the modulus and the phase of the objects of interest. These samples are modeled by a 2D transmittance plane ${\underline { t}}(x,y)$. In that case, the image model will be referred as a non-parametric. For an infinite aperture and aberration free imaging system, this model is the squared modulus of the convolution between the Rayleigh-Sommerfeld propagation kernel $\underline {\boldsymbol {h}}^{\textrm {RS}}_z$ and the transmittance plane ${\underline {{\boldsymbol {t}}}}$ where:

(4)$${\underline{h}}_z^{\textrm{RS}}(x,y)=\frac{z}{i\lambda}\frac{\exp\left(i \frac{2\pi}{\lambda}\sqrt{x^{2}+y^{2}+z^{2}}\right)}{x^{2}+y^{2}+z^{2}}$$

and $\lambda$ is the wavelength of the illumination [37]. In order to account for aberrations in the image formation model, an aberrated PSF model should be used. Assuming a shift invariance of the pupil function with $z$, the Optical Transfer Function (OTF), which is equal to the Fourier Transform of the complex-valued PSF $\underline {\boldsymbol {h}}_z$, can be expressed as follows:

(5)$$\widetilde{\underline{\boldsymbol{h}}_z}(\boldsymbol{\alpha})=\widetilde{\underline{p}}(\boldsymbol{\alpha})\odot\widetilde{\underline{\boldsymbol{h}}}^{\textrm{RS}} _z$$

where $\widetilde {\underline {\boldsymbol {h}}}^{\textrm {RS}} _z$ is also called the angular spectrum.

In addition, the estimated depths of the beads provide reliable estimations of the propagation distance $z^{\dagger }$ [23]. Thus, the aberration corrected non-parametric model $\boldsymbol {m^{\textrm {NP}}}$ can be expressed as:

(6)$$\boldsymbol{m^{\textrm{NP}}}({{{\underline{\boldsymbol{t}}}}}, \boldsymbol{\alpha})=\left| \underline{\boldsymbol{h}}_z(\boldsymbol{\alpha}) \ast {\underline{{\boldsymbol{t}}}} \right|^{2}$$

Unlike the parametric case (section 2.1), minimizing the discrepancy between data and model is not sufficient to solve this ill-posed problem. A priori information about the sample must be added in the form of constraints on the optimization space $\mathbb {S}$ and in the form of a regularization term $\mathscr{R}^{\textrm {NP}}$ [27,38,39]:

(7)$${\underline{{\boldsymbol{t}}}}^{{\dagger}}=\mathop{\mathrm{argmin}}_{{\underline{{\boldsymbol{t}}}} \in \mathbb{S}} \ \lVert {{{\boldsymbol{d}}}} -\boldsymbol{m^{\textrm{NP}}}({\underline{{\boldsymbol{t}}}},\boldsymbol{\alpha})\rVert ^{2} _2 +\mu \mathscr{R}^{\textrm{NP}}({\underline{{\boldsymbol{t}}}})$$

where $\mu$ is an hyperparameter.

This reconstruction is called a regularized reconstruction. The knowledge of the propagation distance $z$ is crucial because the image formation model depends on it. This distance is chosen according to the parameters of the beads previoulsy estimated with the parametric IPA reconstructions [23].

In the following, the regularization term is a hyperbolic total variation term [40]. The hyperparameter is chosen empirically. The optimization domain is restricted to the unitary disk corresponding to a non-emissive object hypothesis. A FISTA algorithm is used to perform this minimization [41].

3. Experimental study

3.1 Principle

High quality microscope objectives are supposed to be diffraction limited as long as they are used in the standard conditions for which they have been optimized (coverslip thickness, refractive indices of the immersion medium, the sample medium and the coverslip and position of the sample relative to the coverslip) [8,42]. Yet, in some applications, these golden rules may be broken (wrong coverslip thickness, for instance). In inset A of Fig. 2, the refraction of the beam in the coverslip is shown before entering the objective. This illustrates the origin of the possible wavefront errors that may occur between the paraxial rays and the high angle rays when the standard conditions of use are not met. This wavefront error has been described by several authors [8,42] in on-axis situations, but it may vary with the position in the field of view. Finally, even when the rules are strictly applied, residual aberrations may still exist, especially out of the optical axis, and may differ from one objective to another. To experimentally study the influence of such aberrations, we used a water immersion objective with a coverslip correction collar. Thus, for a given coverslip thickness, a wrong correction collar setting will give rise to aberrations. This idea was recently proposed by Martin et al. [18].

Fig. 2. Experimental setup. F: monomode fiber coupled laser source, CO: collection optics, P: 200$\mu$m-pinhole, M: mirror, L: lens, Sa: sample that can be precisely moved in XYZ-directions, z: defocus distance of the sample from the focus plane, FP: objective focal plane, MO: microscope objective, BFP: objective back focal plane, TL: tube lens, Se: sensor. Inset A: Zoom on the sample and the objective showing the refraction of the rays occurring through the coverslip. SM: sample medium, C: coverslip, IM: immersion medium. Inset B: Picture of the setup showing the imaging system and the precision piezo-stage (ZS) and the XY-translation stage (XYS).

Download Full Size | PDF

3.2 Setup

Our home-made experimental setup [19] is presented in Fig. 2. The setting of a spatially coherent illumination may be difficult in a microscopy setup as it is very sensitive to any stray reflections or dust particles and leads to complex, sometime unstable speckle patterns. In this setup, the coherent illumination is set by illuminating a 200$\mu m$-pinhole (P) and a lens (L) set in a $2f$ configuration. Thus, an airy pattern illuminates the sample, with a large enough central peak to illuminate the whole field of view, but without inserting too much stray light in the imaging system. This leads to moderate vignetting which is corrected by dividing the holograms by a background intensity image.

In the present study, the sample was composed of 1$\mu$-diameter polystyrene beads diluted in glycerol. The use of transparent beads is interesting here, because we aim at reconstructing phase objects. The diameter of the beads is chosen to mimic biological objects such as bacteria. Usually, sub-resolution objects are used for PSF calibrations. However, in our context, with sub-resolution beads, the contrast of holograms would be too low and Rayleigh-Sommerfeld based models would fail to reconstruct correctly the beads [43].

As polystyrene beads float in glycerol and thanks to its high viscosity, the beads were located just below the coverslip and did not move during the exposure of one hologram (typically, few milliseconds). According to the Gibson-Lanni model of the aberrations [8] induced by wrong coverslip thicknesses and/or refractive indices, the fact that the sample medium was glycerol instead of water should not induce additional aberrations as the beads were just below the coverslip (i.e., $t_s = 0$ with the Gibson-Lanni notations). A coherent illumination with a laser at 637.6nm was used. The illumination power was sufficient to keep exposure times as short as 5ms with our Thorlabs-S805MU1 camera. The sensor pixel size was 5.5$\mu m$. With 22.6mm diagonal, the sensor covers an important part of the field of view of the image (the objective field number is 26.5mm). The microscope objective was a water immersion microscope objective (Olympus PlanSApo, 60$\times$, 1.2NA) with a coverslip correction collar. The tube lens was a 200mm-focal length apochromatic TTL200MP from Thorlabs that was used in a telecentric configuration. The measured magnification was 66.5, and not 60, as the tube lens has a greater focal length than the Olympus standard (180mm).

3.3 Experimental protocol

Five cases of aberration were tested in this experiment with the correction collar at different settings (0.13, 0.15, 0.17, 0.19 and 0.21mm). The coverslip thickness was measured to be 0.170mm with a digital indicator (with a resolution of $\pm 1 \mu m$). Thus, the 0.17mm setting of the correction collar is assumed to be the aberration free situation. A single bead was tracked through the whole field of view in regular steps in the X and Y directions. A total of 35 images ($7 \times 5$) were acquired in order to regularly cover the whole field of view ($273 \times 204 \mu m$). For each XY-position in the field, an axial stack was recorded with defocus positions ranging from $-10 \mu m$ to $+20 \mu m$ from the focus position with a step size of $0.5\mu m$. This stack is used for the illustration of Fig. 3, but only one axial position will be reconstructed in the next section. It should be noted that the sample is the only moving part, which is important for recording a background image by calculating the median value of the 35 XY-shifted images recorded at focus.

Fig. 3. Example of a mosaic of holograms (top) of 1$\mu$m-diameter polystyrene beads in glycerol for an approximate defocus of 12$\mu$m under 5 different settings of the correction collar (from left to right: 0.13 (green), 0.15 (yellow), 0.17 (red), 0.19 (blue) and 0.21mm (magenta)). XZ-views of the hologram stacks for the different correction collar settings (bottom).

Download Full Size | PDF

A view of a typical hologram is shown in the top part of Fig. 3. XZ-views of the stack along the vertical axis of the bead are represented at the bottom of Fig. 3. As shown in the figure, a change in the focus position is observed as a function of the setting of the correction collar, as well as modifications in the XY-profiles. The radial symmetry of the PSF is not always valid, as can be seen, for example, for the 0.13mm setting of the correction collar (green). This asymmetry is due to aberration effects that may break the radial symmetry of the holograms (e.g., coma, astigmatism, etc.). All aberrations may originate from the objective, but also from the tube lens or from misalignment of the illumination. Moreover, aberrations can also originate from inhomogeneities of the slide and the coverslip.

4. Results

In this section, we first apply the proposed method to simulated holograms to demonstrate the robustness of our approach for several kind of aberrations, especially in cases of difficult optimizations, i.e., with highly correlated Zernike coefficients. We then apply it to experimental holograms of beads. We compare our results with state-of-the-art parametric reconstruction algorithms in both simulated and experimental cases and finally evaluate and discuss the effects of aberration on regularized reconstructions.

4.1 Reconstructions on simulated data

A mosaic of $7 \times 5$ in-line holograms was simulated with aberrations varying in the field of view (see Fig. 4). Each hologram is a 512$\times$512 pixels sub-image simulated with the experimental parameters described in Table 1 and with the aberrated Lorenz-Mie model (see Eq. (1)). The defocus is set to 12$\mu$m. This distance was chosen to improve the accuracy of the estimation of the Zernike coefficients, as indicated by CRLB analysis of this parameter (see Appendix B Fig. 11).

Fig. 4. Top: 35 holograms simulated with variable Zernike coefficients depending on the position in the field of view. Bottom: magnifications of 3 holograms from different areas (first line), estimated model accounting for aberrations (C) (second line), residuals, i.e., difference between the first line and the second one (third line).

Download Full Size | PDF

Table 1. Experimental parameters (*from manufacturer, ThermoFisher Scientific, Inc.)

View Table | View all tables in this article

To simulate a varying PSF in the field of view, the aberrated pupil function was considered to depend on the position of the bead in the field of view. This pupil function corresponds to a linear combination of oblique astigmatism ($\boldsymbol {Z}_2^{-2}$), vertical coma ($\boldsymbol {Z}_3^{-1}$), horizontal coma ($\boldsymbol {Z}_3^{1}$), spherical aberration ($\boldsymbol {Z}_4^{0}$) and oblique secondary astigmatism ($\boldsymbol {Z}_4^{-2}$) (see Appendix A). This linear combination is weighted by the corresponding Zernike coefficients $\boldsymbol {\alpha }$ (see Section 2.1). We arbitrarily chose to set a linear behavior along $y$ for $\alpha _2^{-2}$ and $\alpha _3^{1}$, a linear behavior along $x$ for $\alpha _3^{-1}$ and $\alpha _4^{-2}$, and we set $\alpha _4^{0}$ constant in the field of view. This set of coefficients was chosen to demonstrate the performance of the proposed method in difficult cases, i.e., we chose Zernike coefficients that were highly correlated in the corrected model (see Appendix B, Table 6).

Finally, a white Gaussian noise $\bf {\epsilon }{}$ was added to the simulated holograms, which led to a Signal-to-Noise Ratio (SNR) of 4 in the holograms ($\textrm {SNR}=\frac {\Delta {\boldsymbol {m}}}{2\sigma _{\bf {\epsilon }{}}}$, where $\Delta {\boldsymbol {m}}$ is the peak-to-peak amplitude of the model).

For each simulated hologram, the reconstruction was performed using parametric IPA with or without aberration corrections in the model. The abbreviations C (standing for corrected), and UC (standing for uncorrected) will be used in the following. The optimization algorithm we used was the LINCOA algorithm [44]. To perform the reconstructions with the corrected model (C), the first step implies an exhaustive search in a 17 parameters space, which can be really demanding in terms of computational time. To reduce this exhaustive search, it can be fairly convenient to have at least a coarse knowledge of the Zernike coefficients. As our aberrations were quite low, we performed this step by considering no aberration, i.e., all Zernike coefficients were set to zero.

Then, the optimization step was performed with the fully corrected model (eq. (3)), with the constraints on parameters described in Table 2. The optimization domains $\left \{\mathbb {P},\mathbb {D}\right \}$ were chosen quite large in order to check the robustness of the proposed method.

Table 2. Optimization constraints for each estimated parameters ($z$ and $r$ are in micrometers)

View Table | View all tables in this article

Table 3 shows the bead parameters reconstructed without (UC) or with (C) taking the aberration into account in the model. It shows the biases introduced by geometrical aberrations. When using an unaberrated model (UC), the reconstructions converge either on a local optimization minimum or to the constraint domain bounds. Conversely, when using an aberrated model, the reconstructions always converge to the global minimum with low bias and a standard deviation close to the theoretical lower bound given by Cramér-Rao analysis.

Table 3. Statistical results on the estimated bead parameters with aberration corrected (C) and uncorrected (UC) models^a

View Table | View all tables in this article

Residuals between the data and the model are very low, indicating that the model fits the data accurately (see the bottom line in Fig. 4). On the upper part of Fig. 5 are presented the phase of the pupil functions that were simulated in each part of the field of view. This gives another view, in Fourier space, of the type of phase errors that aberrations may imply. On the lower part of Fig. 5 the residuals of the estimated pupil functions are presented (from the simulated ground truth). From these residuals, we see that our estimations of the Zernike coefficients are accurately describing the phase function introduced by aberrations in Fourier space.

Fig. 5. Simulated phase correction in the pupil planes $\widetilde {\underline {p}}(\boldsymbol {\alpha }^{\dagger })$ (a) and residuals of the estimated pupil functions from the ground truth (b). The white dashed circles correspond to the disk in which 95% of the energy of the power spectrum of the object (c) is contained. The black dashed circles correspond to the aperture (calculated from the numerical aperture of the objective).

Download Full Size | PDF

In the most difficult cases (upper part and lower part of the field on Fig. 5), the residuals are not negligible for the highest spatial frequencies, close to the cutoff frequency imposed by the numerical aperture of the objective (represented by a black dashed circle). Indeed, as we did not use a sub-resolution object, the power spectrum of the object is not filling the entire pupil. In the inset of Fig. 5, the typical power spectrum of the object is presented and a white dashed circle shows the part of the spectrum including 95% of its energy. In this white dashed circle, the residuals remain low. Actually, this is an unsurprising limitation of this approach: as the object spectrum does not cover the whole aperture of the objective, the pupil phase function can not be estimated precisely for the highest frequencies. However, the pupil function is correctly estimated for the spatial frequencies corresponding to the spectrum of the object, which ensures that a similar object will be correctly reconstructed. If the aberrations are important, this effect must be considered for the choice of the calibration objects: the size of the beads chosen for aberration estimation must be at least equal or smaller than the smallest detail of interest.

4.2 Reconstructions on experimental data

The experimental parameters are given in Table 1 and were the same as those used in the simulations. Once again, since the accuracy of the estimated parameters is better in a specific range of defocus $z$ (see Appendix B), the holograms to be reconstructed were located approximately 12$\mu$m from the focus position, as in the simulations. They were reconstructed using parametric IPA, with the same workflow that was described in the reconstructions of the previous subsection. Again, to compare the effect of aberrations on the estimation of the beads parameters, both corrected (C) and uncorrected (UC) models are used for the reconstructions. As illustrated in Section 3, the position of the focus varied with the setting of the correction collar. Parametric IPA provides an estimation of the defocus distance $z$ between the sample and the focal plane of the objective.

Table 4 presents a list of the mean values and the standard deviations of all 35 positions in the field for parameters $z$, $r$ and $n$ and for both (UC) and (C) reconstructions. According to the comparison of standard deviations for each collar setting, the dispersion over the field was only moderately modified by the model (UC) or (C). However, the mean values changed, especially that of the estimated defocus $\hat {z}$. A maximum difference of 1.68$\mu$m in the estimated defocus was found between the two models (UC) and (C).

Table 4. For the 5 correction collar settings, averages $<\hat {\vartheta _i}>$ and standard deviations $\sigma _{\hat {\vartheta _i}}$ of the estimated parameter $\hat {z}$, $\hat {r}$ and $\hat {n}$ using (UC) uncorrected model and (C) aberration corrected model.^a

View Table | View all tables in this article

Moreover, the estimated defocus highly depends on the correction collar setting, which varied from 8.8 to 13.9$\mu$m (UC) and from 10.5 to 12.5$\mu$m (C). Thus, this dispersion was reduced by taking the aberrations into account, indicating a correction of the bias in the evaluation of the defocus. Since regularized reconstruction algorithms rely on a precise knowledge of the image formation model (including the defocus distance), any misestimation of the axial position of the sample would bias the reconstructions. Finally, it must be noted that the remaining dependence of the estimated defocus with the correction collar setting may have a physical origin. Indeed, wrong settings of the correction collar may really change the focus position as it changes the properties of the objective.

For the estimated radii $\hat {r}$ and refractive indices $\hat {n}$, the dispersion over the field was reduced when the aberrations were taken into account. The averages were also less dispersed, but to a lesser extent. Indeed, some biases that depend on the correction collar setting appeared to remain.

Figure 6 presents the estimated bead parameters as scatter plots. This makes it possible to visualize the correlations between the estimated parameters $z$, $r$ and $n$.

Fig. 6. Scatter plots showing the biases and correlations between the estimated defocus $\hat {z}$, radius $\hat {r}$ and refractive index $\hat {n}$ for a single bead, for the 35 positions in the field, for the 5 settings of the coverslip correction collar and with corrected models (C) and uncorrected models (UC). With correction of the aberrations, the bias and the dispersion of the estimations due to aberrations are reduced.

Download Full Size | PDF

Moreover, for each collar setting (one color for one collar setting), taking the aberrations into account improved the repeatability of the parameter estimation independently of the introduced aberrations. Indeed, the aberration corrections not only reduce the biases between the different collar settings (differences from one color point cloud to another) but also reduce correlations coefficients between parameter estimations (correlations within one color point cloud). This is presented quantitatively on Table 5, for both models (C) and (UC) and for the less aberrated case (0.17mm). According to Table 5, the decorrelation is particularly important between $r$ and $n$.

Table 5. Correlation coefficients between the estimated parameters without aberration correction (left) and with aberration correction (right) for a correction collar setting of 0.17mm (less aberrated case)

View Table | View all tables in this article

According to the manufacturer’s specifications, the radius should be (0.5$\pm$0.03)$\mu$m and the refractive index should be around 1.587. The estimated parameters obtained with or without an aberration model were within the manufacturer’s confidence interval (0.47-0.53$\mu$m). It is important to note that the fit with the Mie model is constrained by the spherical hypothesis and thus may be quite robust to errors in the model, contrary to the case of regularized reconstruction that have more degrees of freedom, and will be more sensitive to aberrations, especially to non-radially symmetric ones, as it will be seen later on regularized reconstructions.

With the 35 recorded holograms corresponding to 35 bead positions in the field of view, we were able to check that the Zernike coefficients vary in the field of view, following continuous evolutions similar to those described in another work [10]. The Fig. 7, illustrates the evolution of the Zernike coefficients associated with oblique astigmatism, defocus, vertical coma, horizontal coma and spherical aberrations. These appeared to be the main components of the aberrated pupil function $\widetilde {\underline {p}}$. The evolution of these coefficients is continuous and, not surprisingly, increases with increasing errors in the correction collar setting. Vertical coma increases from the left to the right whereas horizontal coma increases from the top to the bottom of the field of view. Spherical aberration and defocus do not depend on the location in the field of view but change with the correction collar setting, with almost no spherical aberration and defocus for the less aberrated case (0.17mm). This is quite logical as a coverslip thickness error is known to induce spherical aberrations [18]. On the contrary, oblique astigmatism varies in the field of view without depending too much on the correction collar setting.

Fig. 7. Estimated Zernike coefficients $\alpha _n^{m}$ as a function of the position in the field of view and for 3 settings of the coverslip correction collar (0.13 mm, 0.17 mm, 0.21 mm). The evolution of the Zernike coefficients is continuous in the field of view. The coma and astigmatism coefficients depend on the position in the field of view and on the correction collar setting whereas defocus and spherical aberration only depend on the correction collar setting.

Download Full Size | PDF

Figure 8 illustrates the evolution of the phase correction for the 35 positions in the field of view and for a correction collar of 0.17mm. For this supposedly aberration-free case, the setup still suffer from aberrations that change in the field of view. These phase functions show significant aberration effects but, as expected, lower than for the other correction collar settings (not represented). This indicates the necessity of taking aberrations into account for hologram reconstruction even when the optical system is supposed to be compensated for aberrations. Indeed, these aberrations may come from residual aberrations of the objective, but also from other sources, like thickness inhomogeneities of the slide and the coverslip, as well as alignment issues.

Fig. 8. Evolution of the phase (in radians) of the pupil function correction in the field of view for a setting of the coverslip correction collar of 0.17mm and for the 7$\times$5 positions in the field where the aberrations were estimated. The black and white dashed circles are defined on Fig. 5

Download Full Size | PDF

From the numerical point of view, the detection of all 35 beads in the mosaic takes around 30 seconds on a 3296$\times$2472 pixels image. The local optimization step for each bead takes around 10 seconds when not considering aberrations while it takes 45 seconds when considering them. These estimations have been realized using an Intel Core i9-11950H CPU 2.60GHz with 16GBytes of RAM.

4.3 Reconstructions on experimental data using regularized algorithms

The evaluation of aberration’s effects on reconstructions is performed using beads holograms. This allows us to compare quantitatively the reconstructed transmittance with a ground truth (assumed to be the transmittance of the bead whose parameters are estimated by parametric IPA). However, since the non-parametric model is very general (not limited to spherical objects), similar results will be obtained with an aspherical sample. The reconstruction is performed with (C) and without (UC) the previously estimated aberration pupil function $\widetilde {\underline {p}}$ and the $\hat {z}$ parameters.

A Fienup phase retrieval algorithm [28,29], as well as a regularized IPA (as presented in 2.2) are used to reconstruct the data. These reconstructions are performed using the uncorrected propagator $\underline {\boldsymbol {h}}^{\textrm {RS}}_z$ (UC) or the corrected propagator $\underline {\boldsymbol {h}}_z$ (C) in the model (eq. (6)). Figure 9 illustrates the reconstructions results for both algorithms. The estimated aberrated Mie model that fits the data has been back propagated at the center plane of the bead (BPMie-C) and is considered as the ground truth here because it is the most accurate model. Similarly, a back propagation of the Mie model estimated without aberration has also been computed (BPMie-UC). Because of the coma aberrations, the bead position $(x,y)$ is not the same for (BPMie-C) and (BPMie-UC) parametric inversions, as mentioned in Appendix B. For comparison purpose the beads have then been centered in Fig. 9.

Fig. 9. Non-parametric reconstructions using regularized IPA and Fienup algorithm with (C) or without (UC) aberrations correction. The reconstructions are presented in real part an imaginary part. A reconstruction is compared with the back-propagation of the estimated Mie model without aberration estimation (BPMie-UC) and the back-propagation of the Mie model with aberration estimation (BPMie-C). Profiles of the real part and imaginary part at the center of the bead are presented.

Download Full Size | PDF

When aberrations are not considered in the reconstruction model, the morphological properties and quantitativity of the reconstructions are compromised. Indeed, either with Fienup or regularized IPA, the bead does not show a circular shape. As the aberrations vary in the field of view, the same bead does not have the same shape for each lateral position. The back-propagation of the Mie model without aberration illustrates the model error when the aberrations are not considered, but the radial symmetry is maintained as the Mie model is based on a spherical model. The regularized reconstructions without aberrations do not match with this model indicating bias in the estimation of the bead parameters. However, with aberrations correction the reconstructions fit the corresponding back-propagated Mie model and have the expected geometrical and quantitative properties. It demonstrates that whatever the reconstruction algorithm, aberrations should be taken into account to restore accurately the morphological and quantitative properties of the sample.

Taking into account the aberrations in regularized reconstructions has no effect on the computational time as the aberrated forward model has the same complexity as angular spectrum propagation. In the example of Fig. 9, reconstructing a whole field of view (2472 $\times$ 3296) and considering the spatial evolution of the PSF takes less 10 minutes. These estimations have been realized using an Intel Core i9-11950H CPU 2.60GHz with 16GBytes of RAM. This computational time can be reduced using GPU.

5. Conclusion

In this article, we present a method to estimate the aberrations and thus reduce reconstruction errors, by using a more accurate image formation model, in in-line holographic microscopy. This method is based on the use of calibration beads. We show that the rigorous and highly constrained Mie model can be used to estimate bead parameters and Zernike coefficients at the same time with a good precision and repeatability. Moreover, this approach requires only one hologram and does not require any assumption on the PSF evolution in the field of view. This calibration step could be done sequentially, like standard calibrations or in-situ by inserting calibrated beads in the biological sample itself. However, this may depend on the application or on the main origin of the aberrations (from the optical setup or from the sample itself). Actually, adding calibration beads in the sample has already proven to be useful for autofocusing [23]. In this context, with the present method of correction of aberrations, this autofocusing would be even more accurate.

Once the Zernike coefficients have been estimated, it is then possible to use them in a regularized approach framework to reconstruct any biological objects (spherical or not), as long as the sparsity constraint required in in-line digital holography is fullfilled. This methodology of aberration estimation was applied for the improvement of regularized reconstruction of holograms with the in-line holographic microscopy configuration. However, it is also applicable to off-axis holography or other coherent imaging techniques or simply used as a calibration method for microscopy systems.

The method proposed here offers interesting perspectives for reconstructing more accurately and with more quantitativity the absorption and the phase of the objects of interest, even with poorly corrected or misaligned optical systems, non-standard optical configurations (various sample media, variable axial position of the objects below the coverslip) and more generally, for any non-standard microscopy configurations that may introduce aberrations.

In this study, we estimated aberrations parameters on a discrete grid. The next step could be to interpolate the spatially varying PSF. This continuous PSF can then be used in the image reconstruction step, but with a higher computational cost. Nevertheless, fast algorithms can be used [45,46].

Appendix A: Zernike polynomials

Zernike polynomials depend on two parameters: the azimuthal angle $\phi =\arctan \left ( \frac {\kappa _y}{\kappa _x}\right )$ and the normalized radial distance $\rho =\frac {\lambda }{\textrm {NA}}\sqrt {\kappa _x^{2}+\kappa _y^{2}}$ and are defined as follows :

Z^{m}_n(\rho,\phi)= \left\{ \begin{array}{ll} R_n^{\lvert m \rvert}(\rho)\sin(m\phi)\ \textrm{if}\ m>0\\ R_n^{\lvert m \rvert}(\rho)\cos(m\phi)\ \textrm{otherwise} \end{array} \right.

where $n\in \mathbb {N}$, $m \in \mathbb {Z}$ and $R^{m}_n(\rho )$ is defined as :

R^{m}_n(\rho)=\sum_{k=0}^{\frac{n-m}{2}}\frac{({-}1)^{k}(n-k)!}{k!\left[\frac{n+m}{2}-k\right]!\left[\frac{n-m}{2}-k\right]!}\rho^{n-2k}

with $n\geq \lvert m \rvert$ and $n-\lvert m \rvert$ even.

Because of the numerical aperture, the pupil function is zero out of the disk defined by $\rho \leq 1$. An illustration of the polynomials is given on Fig. 10 [47].

Fig. 10. Illustration of the 15 first Zernike polynomials (adapted from [47]).

Download Full Size | PDF

Appendix B: Theoretical study of the aberration parameters accuracy

In this appendix, we aim at estimating the achievable precision on each estimated parameter and to study the correlation between these parameters. For these purposes, Cramér-Rao Lower Bounds (CRLB) and the correlation matrix are computed [48] using our aberrated model $\boldsymbol {m}^{\textrm {P}}$ presented in section 2.1 Eq. (1). According to Cramér-Rao inequality, the variance of any unbiased estimator $\hat {{\boldsymbol {\zeta }}}=\{\zeta _i\}_i=\{{\boldsymbol {\vartheta }}, \boldsymbol {\alpha } \}$ of the unknown vector parameter ${\boldsymbol {\zeta }}^{\dagger }$ is bounded from below by the i-th diagonal coefficient of the inverse of the Fisher information matrix:

(8)$$\textrm{Var}\left( \hat{\zeta_i} \right) \geq \left[\boldsymbol{I}^{{-}1}({\boldsymbol{\zeta}}^{{\dagger}})\right]_{i,i} = {\sigma^{\textrm{CRLB}}_{\zeta_i}}^{2}$$

where $\boldsymbol {I}({\boldsymbol {\zeta }}^{\dagger })$ is the Fisher information matrix. It is linked to the curvature of the cost function in the parameters space:

(9)$$\left[\boldsymbol{I}({\boldsymbol{\zeta}}) \right]_{i,j} = E \left[ \left| \frac{\partial^{2} \mathscr{D}^{\textrm{P}}({{{\boldsymbol{d}}}},\boldsymbol{m}^{\textrm{P}}({\cdot}))}{\partial \zeta_i \partial \zeta_j} \right|_{{\boldsymbol{\zeta}}} \right]$$

In the case of white Gaussian noise of standard deviation $\sigma _\epsilon$, neglecting quantization effect and considering a centered model [49] :

(10)$$\left[\boldsymbol{I}({\boldsymbol{\zeta}}) \right]_{i,j} =\frac{1}{\sigma_\epsilon^{2}} \sum_k\left( \frac{\partial \boldsymbol{m}^{\textrm{P}}(x_k,y_k,{\boldsymbol{\zeta}})}{\partial \zeta_i}\frac{\partial \boldsymbol{m}^{\textrm{P}}(x_k,y_k,{\boldsymbol{\zeta}})}{\partial \zeta_j}\right)$$

These bounds are computed for a bead at the center of the field of view and for several defocus distances with parameters of Table 1 (${\boldsymbol {\vartheta }}\left (x=0\; \mu m, y=0\; \mu m , z, r=0.5\; \mu m, n=1.58\right )$).

As the aberrations happen to be quite low in our case, the accuracy study has been performed with Zernike coefficient set to zero. Thus, the accuracy on the Zernike coefficents has been studied around a zero value.

Figure 11 illustrates the evolution of the CRLB with the propagation distance $z$ (i.e., the lower bound variance of each parameter versus $z$ value).

Fig. 11. Evolution of Cramér-Rao Lower Bounds on each parameter as a function of $z$. The experimental parameters of the model are given in Table 1. These CRLB have been computed for a hologram without aberrations.

Download Full Size | PDF

These CRLB have been computed considering $\sigma _\epsilon$ constant and using numerical derivatives. For most parameters the best accuracy is obtained for defocus distances between 10 and 15 $\mu m$. In this study, the defocus distance $z=12 \mu m$ was considered.

The correlation matrix $\boldsymbol {\Sigma }$ is obtained by:

(11)$$\boldsymbol{\Sigma}_{i,j}\left({\boldsymbol{\zeta}}\right)=\frac{\left[\boldsymbol{I}^{{-}1}({\boldsymbol{\zeta}})\right]_{i,j}}{\sigma^{\textrm{CRLB}}_{\zeta_i}.\sigma^{\textrm{CRLB}}_{\zeta_j}}$$

Table 6 shows this correlation matrix for the selected seventeen parameters. Coefficients below 0.05 are set to zero for a better visualization.

Table 6. Correlation matrix of the 5 beads parameters and 12 Zernike coefficients.^a

View Table | View all tables in this article

The correlation matrix indicates strong correlations between several parameters. Unsurprisingly, $r$ and $n$ are highly correlated as the phase shift induced by a an object depends on the product of these two parameters and the phase shift has a strong effect on the propagation. It is interesting to notice that coma coefficients represented by $\alpha ^{-1}_3$ and $\alpha ^{1}_3$ are highly correlated with $x$ and $y$. Therefore, ignoring the coma aberration could lead to lateral shifts in the reconstructions. Correlations between $\alpha ^{-4}_2$ and $\alpha ^{-2}_2$, $\alpha ^{0}_4$ and $\alpha ^{0}_2$ or $\alpha ^{2}_4$ and $\alpha ^{2}_2$, may lead to misestimations of these coefficients. This is studied in section 4.1 on simulation experiments. However, it would be probably worse not to take them into account because that would systematically introduce errors in the model. Most of the other coefficients of the correlation matrix are low or null and the corresponding parameters can then be considered as decorrelated. Because of the high correlation values in the correlation matrix, all parameters must be estimated at the same time to prevent estimation errors.

Funding

Agence Nationale de la Recherche (ANR-11-IDEX-0007, ANR-11-LABX-0063, ANR-18-CE45-0010).

Acknowledgment

The algorithmic tools presented in this work have been implemented within the framework of the Matlab library GlobalBioIm [50,51]. The Zernike polynomials models have been computed using Fricker’s implementation [52]. LINCOA optimization strategies have been computed using the PDFO library [53]. This work has been funded by the Auvergne-Rhône-Alpes region, France, under project DIAGHOLO. It was also performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program “Investissements d’Avenir” (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR) and ANR HORUS (ANR-18-CE45-0010).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request. Information related to GlobalBioIm can be found at [54].

References

1. R. Alexander, B. Leahy, and V. N. Manoharan, “Precise measurements in digital holographic microscopy by modeling the optical train,” J. Appl. Phys. 128(6), 060902 (2020). [CrossRef]

2. J. G. McNally, T. Karpova, J. Cooper, and J. A. Conchello, “Three-Dimensional Imaging by Deconvolution Microscopy,” Methods 19(3), 373–385 (1999). [CrossRef]

3. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Process. Mag. 23(3), 32–45 (2006). [CrossRef]

4. F. Soulez, L. Denis, Y. Tourneur, and E. Thiebaut, “Blind deconvolution of 3D data in wide field fluorescence microscopy,” in 9th International Symposium on Biomedical Imaging (ISBI) (IEEE, 2012), pp. 1735–1738.

5. J. Li, F. Xue, and T. Blu, “Fast and accurate three-dimensional point spread function computation for fluorescence microscopy,” J. Opt. Soc. Am. A 34(6), 1029–1034 (2017). Publisher: Optica Publishing Group. [CrossRef]

6. A. Aristov, B. Lelandais, E. Rensen, and C. Zimmer, “ZOLA-3D allows flexible 3D localization microscopy over an adjustable axial range,” Nat. Commun. 9(1), 2409 (2018). [CrossRef]

7. J. Li, F. Xue, F. Qu, Y.-P. Ho, and T. Blu, “On-the-fly estimation of a microscopy point spread function,” Opt. Express 26(20), 26120–26133 (2018). Publisher: Optica Publishing Group. [CrossRef]

8. S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc. Am. A 8(10), 1601–1613 (1991). [CrossRef]

9. V. Lakshminarayanan and A. Fleck, “Zernike polynomials: a guide,” J. Mod. Opt. 58(7), 545–561 (2011). [CrossRef]

10. G. Zheng, X. Ou, R. Horstmeyer, and C. Yang, “Characterization of spatially varying aberrations for wide field-of-view microscopy,” Opt. Express 21(13), 15131 (2013). [CrossRef]

11. P. Ferraro, S. D. Nicola, A. Finizio, G. Coppola, S. Grilli, C. Magro, and G. Pierattini, “Compensation of the inherent wave front curvature in digital holographic coherent microscopy for quantitative phase-contrast imaging,” Appl. Opt. 42(11), 1938–1946 (2003). Publisher: Optica Publishing Group. [CrossRef]

12. T. Colomb, E. Cuche, F. Charrière, J. Kühn, N. Aspert, F. Montfort, P. Marquet, and C. Depeursinge, “Automatic procedure for aberration compensation in digital holographic microscopy and applications to specimen shape compensation,” Appl. Opt. 45(5), 851–863 (2006). [CrossRef]

13. J. Min, B. Yao, S. Ketelhut, C. Engwer, B. Greve, and B. Kemper, “Simple and fast spectral domain algorithm for quantitative phase imaging of living cells with digital holographic microscopy,” Opt. Lett. 42(2), 227–230 (2017). [CrossRef]

14. L. Xu, X. Peng, J. Miao, and A. K. Asundi, “Studies of digital microscopic holography with applications to microstructure testing,” Appl. Opt. 40(28), 5046–5051 (2001). Publisher: Optica Publishing Group. [CrossRef]

15. J. Garcia-Sucerquia, W. Xu, S. K. Jericho, P. Klages, M. H. Jericho, and H. J. Kreuzer, “Digital in-line holographic microscopy,” Appl. Opt. 45(5), 836–850 (2006). [CrossRef]

16. T. Kreis, Handbook of Holographic Interferometry: Optical and Digital Methods (John Wiley & Sons, 2006).

17. G. Popescu, Quantitative Phase Imaging of Cells and Tissues (McGraw-Hill Education, 2011).

18. C. Martin, B. Leahy, and V. N. Manoharan, “Improving holographic particle characterization by modeling spherical aberration,” Opt. Express 29(12), 18212 (2021). [CrossRef]

19. T. Olivier, D. Brault, S. Joshi, T. Brard, A. Brodoline, L. Méès, and C. Fournier, “Effects of some model approximations in the reconstructions of digital in-line holograms: simulations, experiments on calibrated objects and model refinement assessment,” Proc. SPIE 12136, 1213603 (2022). [CrossRef]

20. S.-H. Lee, Y. Roichman, G.-R. Yi, S.-H. Kim, S.-M. Yang, A. Van Blaaderen, P. Van Oostrum, and D. G. Grier, “Characterizing and tracking single colloidal particles with video holographic microscopy,” Opt. Express 15(26), 18275–18282 (2007). [CrossRef]

21. F. Soulez, L. Denis, C. Fournier, É. Thiébaut, and C. Goepfert, “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” J. Opt. Soc. Am. A 24(4), 1164–1171 (2007). [CrossRef]

22. F. C. Cheong, B. J. Krishnatreya, and D. G. Grier, “Strategies for three-dimensional particle tracking with holographic video microscopy,” Opt. Express 18(13), 13563–13573 (2010). [CrossRef]

23. D. Brault, C. Fournier, T. Olivier, N. Faure, S. Dixneuf, L. Thibon, L. Mees, and L. Denis, “Automatic numerical focus plane estimation in digital holographic microscopy using calibration beads,” Appl. Opt. 61(5), B345 (2022). [CrossRef]

24. O. Flasseur, C. Fournier, N. Verrier, L. Denis, F. Jolivet, A. Cazier, and T. Lépine, “Self-calibration for lensless color microscopy,” Appl. Opt. 56(13), F189–F199 (2017). [CrossRef]

25. N. Verrier, C. Fournier, L. Méès, and T. Fournel, “In-line particle holography with an astigmatic beam: setup self-calibration using an “inverse problems” approach,” Appl. Opt. 53(27), G147–G156 (2014). [CrossRef]

26. L. Denis, D. Lorenz, E. Thiébaut, C. Fournier, and D. Trede, “Inline hologram reconstruction with sparsity constraints,” Opt. Lett. 34(22), 3475–3477 (2009). [CrossRef]

27. F. Jolivet, F. Momey, L. Denis, L. Méès, N. Faure, N. Grosjean, F. Pinston, J.-L. Marié, and C. Fournier, “Regularized reconstruction of absorbing and phase objects from a single in-line hologram, application to fluid mechanics and micro-biology,” Opt. Express 26(7), 8923–8940 (2018). [CrossRef]

28. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef]

29. T. Latychevskaia and H.-W. Fink, “Solution to the twin image problem in holography,” Phys. Rev. Lett. 98(23), 233901 (2007). [CrossRef]

30. A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, vol. 89 (SIAM, 2005).

31. F. Slimani, G. Gréhan, G. Gouesbet, and D. Allano, “Near-field lorenz-mie theory and its application to microholography,” Appl. Opt. 23(22), 4140–4148 (1984). [CrossRef]

32. F. Soulez, L. Denis, É. Thiébaut, C. Fournier, and C. Goepfert, “Inverse problem approach in particle digital holography: out-of-field particle detection made possible,” J. Opt. Soc. Am. A 24(12), 3708–3716 (2007). [CrossRef]

33. L. Méès, N. Grosjean, D. Chareyron, J.-L. Marié, M. Seifi, and C. Fournier, “Evaporating droplet hologram simulation for digital in-line holography setup with divergent beam,” J. Opt. Soc. Am. A 30(10), 2021–2028 (2013). [CrossRef]

34. R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66(3), 207–211 (1976). [CrossRef]

35. O. Flasseur, L. Denis, É. Thiébaut, T. Olivier, and C. Fournier, “Expaco: detection of an extended pattern under nonstationary correlated noise by patch covariance modeling,” in 27th European Signal Processing Conference (EUSIPCO) (IEEE, 2019), pp. 1–5.

36. D. Brault, L. Denis, S. Dixneuf, T. Olivier, N. Faure, and C. Fournier, “Fast and robust pattern detection: Application to spherical bead localization in holographic microscopy,” accepted in IEEE European Signal Processing Conference (EUSIPCO) (2022).

37. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Co. Publishers, 2004).

38. S. Sotthivirat and J. A. Fessler, “Penalized-likelihood image reconstruction for digital holography,” J. Opt. Soc. Am. A 21(5), 737–750 (2004). [CrossRef]

39. F. Soulez, M. Rostykus, C. Moser, and M. Unser, “A COnstrained Method for lensless Coherent Imaging of thin samples,” Appl. Opt. 61(9), F34–F46 (2022). [CrossRef]

40. P. Charbonnier, L. Blanc-Féraud, G. Aubert, and M. Barlaud, “Deterministic edge-preserving regularization in computed imaging,” IEEE Trans. on Image Process. 6(2), 298–311 (1997). [CrossRef]

41. A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. on Image Process. 18(11), 2419–2434 (2009). [CrossRef]

42. O. Haeberle, “Focusing of light through a stratified medium: a practical approach for computing fluorescence microscope point spread functions. Part II: confocal and multiphoton microscopy,” Opt. Commun. 235(1-3), 1–10 (2004). [CrossRef]

43. J. Dohet-Eraly, L. Méès, T. Olivier, F. Dubois, and C. Fournier, “Analysis of three-dimensional objects in quantitative phase contrast microscopy: a validity study of the planar approximation for spherical particles,” Proc. SPIE 11351, 113510Y (2020). [CrossRef]

44. M. J. Powell, “On fast trust region methods for quadratic models with linear constraints,” Math. Prog. Comp. 7(3), 237–267 (2015). [CrossRef]

45. L. Denis, E. Thiébaut, and F. Soulez, “Fast model of space-variant blurring and its application to deconvolution in astronomy,” in 18th International Conference on Image Processing (IEEE, 2011), pp. 2817–2820.

46. L. Denis, E. Thiébaut, F. Soulez, J.-M. Becker, and R. Mourya, “Fast approximations of shift-variant blur,” Int. J. Comput. Vis. 115(3), 253–278 (2015). [CrossRef]

47. Y. Hsieh, Y. Yu, Y. Lai, M. Hsieh, and Y.-F. Chen, “Integral-based parallel algorithm for the fast generation of the zernike polynomials,” Opt. Express 28(2), 936–947 (2020). [CrossRef]

48. S. Kay, “Fundamentals of statistical signal processing: estimation theory,” Technometrics 37, 465–466 (1993).

49. C. Fournier, L. Denis, and T. Fournel, “On the single point resolution of on-axis digital holography,” J. Opt. Soc. Am. A 27(8), 1856–1862 (2010). [CrossRef]

50. M. Unser, E. Soubies, F. Soulez, M. McCann, and L. Donati, “Globalbioim: A unifying computational framework for solving inverse problems,” in Computational Optical Sensing and Imaging (Optical Society of America, 2017), paper CTu1B–1.

51. E. Soubies, F. Soulez, M. T. McCann, T.-a. Pham, L. Donati, T. Debarre, D. Sage, and M. Unser, “Pocket guide to solve inverse problems with globalbioim,” Inverse Probl. 35(10), 104006 (2019). [CrossRef]

52. P. Fricker, “Zernike polynomials,” MATLAB Central File Exchange (2022), https://www.mathworks.com/matlabcentral/fileexchange/7687-zernike-polynomials.

53. T. M. Ragonneau and Z. Zhang, “PDFO: cross-platform interfaces for Powell’s derivative-free optimization solvers (version 1.1),” Zenodo (2021), https://doi.org/10.5281/zenodo.3887569.

54. Biomedical Imaging Group (EPFL), “GlobalBioIm Library,” Github (2017), https://biomedical-imaging-group.github.io/GlobalBioIm/index.html.

Wavelength	637.6nm
Magnification	66.5
Pixel pitch	83nm
Total field of view	273 $\times$ 204 $μ$ m
Beads diameter*	(1.0 $\pm$ 0.06) $μ$ m
Beads refractive index*	1.587 (polystyrene)
Refractive index of immersion medium	1.47 (glycerol)
Coverslip thickness	(0.170 $\pm$ 0.001)mm
Typical defocus	12 $μ$ m

	$z$	$r$	$n$	$α_{2}^{- 2}$	$α_{2}^{0}$	$α_{2}^{2}$	$α_{3}^{- 3}$	$α_{3}^{- 1}$
Lower bound	10	0.2	1.52	−10	−10	−10	−10	−10
Upper bound	14	0.7	1.63	10	10	10	10	10
	$α_{3}^{1}$	$α_{3}^{3}$	$α_{4}^{- 4}$	$α_{4}^{- 2}$	$α_{4}^{0}$	$α_{4}^{2}$	$α_{4}^{4}$
Lower bound	−10	−10	−10	−10	−10	−10	−10
Upper bound	10	10	10	10	10	10	10

$ϑ_{i}$	$ϑ_{i}^{G T}$	$< \hat{ϑ_{i}} >^{U C}$	$< \hat{ϑ_{i}} >^{C}$	$σ_{ϑ_{i}}^{CRLB}$	$σ_{\hat{ϑ_{i}}}^{U C}$	$σ_{\hat{ϑ_{i}}}^{C}$
$z (μ m)$	12	11.048	12.001	0.002	0.579	0.004
$r (μ m)$	0.5	0.267	0.500	0.001	0.030	0.001
$n$	1.58	1.619	1.5798	0.0006	0.0311	0.0007

	UC	C	UC	C
Collar	$< \hat{z} >$	$< \hat{z} >$	$σ_{z}$	$σ_{z}$
0.13	13.913	12.579	0.802	0.817
0.15	13.049	12.542	0.483	0.440
0.17	11.706	12.046	0.548	0.546
0.19	10.496	11.556	0.510	0.484
0.21	8.838	10.525	0.442	0.495
Collar	$< \hat{r} >$	$< \hat{r} >$	$σ_{r}$	$σ_{r}$
0.13	0.526	0.501	0.008	0.005
0.15	0.519	0.505	0.006	0.005
0.17	0.502	0.513	0.007	0.004
0.19	0.495	0.519	0.007	0.003
0.21	0.497	0.522	0.007	0.012
Collar	$< \hat{n} >$	$< \hat{n} >$	$σ_{n}$	$σ_{n}$
0.13	1.5733	1.5901	0.0041	0.0029
0.15	1.5773	1.5882	0.0027	0.0024
0.17	1.5856	1.5837	0.0044	0.0022
0.19	1.5902	1.5809	0.0036	0.0023
0.21	1.5878	1.5798	0.0032	0.0046

Uncorrected (UC)
$ϑ_{i}$	$z$	$r$	$n$
$z$	1	0.02	0.51
$r$	0.02	1	−0.56
$n$	0.51	−0.56	1

Accurate unsupervised estimation of aberrations in digital holographic microscopy for improved quantitative reconstruction

Abstract

1. Introduction

2. Estimation of the aberration parameters and reconstruction

2.1 Calibration: aberration parameters estimation

2.2 Reconstruction: including aberration model

3. Experimental study

3.1 Principle

3.2 Setup

3.3 Experimental protocol

4. Results

4.1 Reconstructions on simulated data

4.2 Reconstructions on experimental data

4.3 Reconstructions on experimental data using regularized algorithms

5. Conclusion

Appendix A: Zernike polynomials

Appendix B: Theoretical study of the aberration parameters accuracy

Funding

Acknowledgment

Disclosures

Data availability

References

Data availability

Cited By

Figures (11)

Tables (6)

Equations (14)

Optics Express

	$x$	$y$	$z$	$r$	$n$	$α_{2}^{- 2}$	$α_{2}^{0}$	$α_{2}^{2}$	$α_{3}^{- 3}$	$α_{3}^{- 1}$	$α_{3}^{1}$	$α_{3}^{3}$	$α_{4}^{- 4}$	$α_{4}^{- 2}$	$α_{4}^{0}$	$α_{4}^{2}$	$α_{4}^{4}$
$x$	1	0	0	0	0	0	0	0	0	−0.89	0	0	0	0	0	0	0
$y$	0	1	0	0	0	0	0	0	0	0	−0.89	0	0	0	0	0	0
$z$	0	0	1	0.24	−0.13	0	−0.49	0	0	0	0	0	0	0	−0.39	0	0
$r$	0	0	0.24	1	−0.85	0	−0.13	0	0	0	0	0	0	0	0	0	0
$n$	0	0	−0.13	−0.85	1	0	0.22	0	0	0	0	0	0	0	0.05	0	0
$α_{2}^{- 2}$	0	0	0	0	0	1	0	0	0	0	0	0	0	0.72	0	0	0
$α_{2}^{0}$	0	0	−0.49	−0.13	0.22	0	1	0	0	0	0	0	0	0	0.92	0	0
$α_{2}^{2}$	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0.72	0
$α_{3}^{- 3}$	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0	0
$α_{3}^{- 1}$	−0.89	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0	0
$α_{3}^{1}$	0	−0.89	0	0	0	0	0	0	0	0	1	0	0	0	0	0	0
$α_{3}^{3}$	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0	0
$α_{4}^{- 4}$	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	0
$α_{4}^{- 2}$	0	0	0	0	0	0.72	0	0	0	0	0	0	0	1	0	0	0
$α_{4}^{0}$	0	0	−0.39	0	−0.05	0	0.92	0	0	0	0	0	0	0	1	0	0
$α_{4}^{2}$	0	0	0	0	0	0	0	0.72	0	0	0	0	0	0	0	1	0
$α_{4}^{4}$	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1