1. INTRODUCTION

Most unstained biological specimens exhibit low optical scattering, i.e., they have poor contrast under bright-field microscopy. However, these transparent samples induce phase shifts in the transmitted light beam. Microscopy techniques, including phase contrast [1], differential interference contrast (DIC) [2], and Hoffman modulation contrast (HMC) [3], have been developed to translate these phase variations into 2D false-color intensity images of the optically transparent samples. However, these techniques cannot accurately image the 3D shape of live cells in complex, 3D spatial environments. This limitation results in a significant gap in our understanding of dynamic changes occurring in the 3D cell shape and behavior of thick unstained specimens. Digital Holographic Microscopy (DHM) overcomes the imaging depth limitations and distortions of the aforementioned microscopic methods [4–7]. DHM is based on optical interferometry to reconstruct both amplitude and phase maps of biological specimens, thereby providing functional and morphological sample information in a non-invasive way. DHM systems have become a powerful quantitative phase imaging (QPI) method due to their broad applicability in biomedical studies [8–17].

Because DHM systems can be utilized for quantitative analysis of biological systems and disease diagnostics, its accuracy for retrieving the phase information is with no doubt a factor to be followed closely. DHM heavily relies on computational processing to accurately reconstruct the phase of the biological sample. The mandatory numerical steps used in off-axis DHM technique are spatial frequency filtering [18,19] of the hologram spectrum and the compensation of the tilt between the reference and object waves. Depending on the optical configuration of the DHM system, numerical propagation and/or aberration correction may be required. The DHM computational processing should be performed automatically and be adaptable to different sample and imaging conditions to increase its applicability in biomedicine. An inaccurate compensation of the tilt between the reference and object waves could introduce errors in the quantitative phase measurements. One straightforward compensation method for precise extraction of the phase is manually selecting the tilt of the reference wave. This approach requires the manual generation of several reference waves until obtaining a phase reconstruction without sawtooth fringes. Selecting the best reconstructed phase image can be a challenging task for the end user in DHM since it depends significantly on the user’s expertise. Several research works have focused on developing computational approaches to automatically reconstruct phase images without any phase perturbations automatically. In 2004, Carl et al. proposed an automatic method for reconstructed phase images in an off-axis DHM system [20]. In their method, the parameters to compensate the tilt and spherical wave introduced by the non-telecentric configuration of the DHM system were estimated from a recorded hologram in which no sample was used (e.g., blank hologram). This method was implemented following an iterative process in which the optimized parameters were obtained by minimizing the standard deviation of the unwrapped reconstructed phase image in the blank hologram. The major inconveniences of this method are the computational time, the need for a blank hologram, the user’s input of the size of the object spectrum, and the need for an unwrapping method. The drawbacks of a blank hologram and the use of an unwrapping method were overcome in an upgraded version of this method [21], and the estimation of the size of object spectrum can be overcome by thresholding the spectrum of the hologram [19,22]. The main difference between these two approaches is the estimation of the tilt parameters. Whereas the authors in Ref. [19] used a centroid-based algorithm to determine the tilt parameters, Trujillo et al. generated multiple tilted plane waves with different parameters and searched for the reconstructed phase image without sawtooth fringes [22]. Whereas these algorithms aim to correct phase aberrations due to the interference tilt and the spherical wave introduced by the DHM system operating in non-telecentric configuration, automatic reconstruction phase methods based on polynomial curve fitting can correct high-order aberrations. In 2006, Colomb et al. proposed an automatic method in which the coefficients of the phase aberrations were estimated by fitting selected line profiles into polynomial curves [23]. Originally, these line profiles were extracted from regions in the unwrapped reconstructed phase image where the sample information was constant. To avoid the user’s input in selecting these flat regions, Nguyen et al. proposed the use of a convolutional neural network (CNN) to perform automatic background region detection [24]. The performance of these fitting-based automatic methods are significantly dependent on the size of the flat regions as well as the performance of the unwrapping algorithm. In 2014, Liu et al. avoided these issues by proposing a nonlinear optimization procedure to estimate the Zernike coefficients of the phase aberrations [25]. These coefficients were estimated by maximizing the spectral energy metric of a single recorded hologram.

In this paper, we present a different type of nonlinear optimization procedure to provide fast, accurate, and automatic compensation of the interference tilt in QPI-DHM. The proposed approach is based on the minimization of a cost function to provide precise phase measurements in off-axis DHM systems operating in telecentric regime. Telecentric-based DHM systems are linear shift-invariant imaging systems [26], leading to accurate QPI results [27] without the insertion of additional optical elements [23,28], the recording of an additional hologram [29], or the application of any computational compensation algorithm [23,28,30–32]. The proposed method has been implemented in a user-friendly tool, and it is available at https://oirl.github.io/tuDHM/. We intend to offer a more accessible reconstruction tool, which could increase the applicability of DHM systems. Preliminary results of this computational approach have been presented in the OSA Imaging and Applied Optics Congress [33]. The developed software tool will allow researchers in life and material sciences, even those without computational reconstruction knowledge, to analyze their results accurately, leading to new discoveries.

2. OFF-AXIS DIGITAL HOLOGRAPHIC MICROSCOPY IN A TELECENTRIC REGIME

Digital holographic microscopy is a hybrid imaging technique based on two stages: the optical recording of a hologram and its numerical reconstruction. DHM systems are optical interferometers. Figure 1 shows the optical configuration of a standard Mach–Zehnder interferometer used for imaging biological samples. The light beam emitted by a laser of wavelength $ \lambda $ impinges on a beam-splitter (BS). One of the split beams illuminates the sample, whose complex amplitude distribution, $o(x,y)$, is set at the front focal plane (FFP) of a microscope objective (MO) lens. The imaging system in Fig. 1 is composed of an infinity-corrected MO lens and a tube lens (TL). The sample image is then obtained at the back focal plane (BFP) of the TL. Commonly, this plane is named the image plane (IP) of the optical microscope. Without loss of generality, the aperture stop of the MO lens is set at the FFP of the TL. Therefore, the DHM system of Fig. 1 operates in the telecentric regime and follows the IP architecture (i.e., focused images of the sample are formed onto the sensor plane). The sensor records the interference pattern between a plane wave with amplitude distribution $r(x,y)$, named the reference wave, and the complex wavefield diffracted by the sample at the sensor plane, $o^\prime (x,y) = o(x,y) \otimes {\rm CPSF} (x,y)$, where ${\rm CPSF} (x,y)$ is the coherent (i.e., amplitude) point spread function of the optical microscope. This interference pattern, commonly called hologram, is given by

Equation (3) represents the filtered hologram spectrum, which is the spectrum of the sample displaced at the spatial frequencies, (${\sin}{\theta _x}/\lambda$, ${\sin}{\theta _y}/\lambda$). The amplitude distribution scattered by the sample can be obtained as the absolute value of the inverse Fourier transform of Eq. (3). However, the phase distribution requires the compensation of the lateral displacement introduced by the tilted reference beam. Thus, for QPI-DHM measurements, one should compensate the reference tilt. This compensation can be performed in both real and Fourier space. To compensate the reference angle using the Fourier space, the filtered hologram spectrum [Eq. (3)] should be centered on a new matrix. In the real space, the reference angle is compensated by multiplying the inverse Fourier transform of Eq. (3), ${h_F}({\boldsymbol x})$, and a digital replica of the reference wave, ${r_D}({\boldsymbol x})$, commonly called the digital reference wave. If the sensor plane coincides with the IP of the microscope, the reconstructed phase image (${\phi _{\hat o^\prime}}$) is given by the angle of

Otherwise, when the hologram is not recorded under IP conditions, one must numerically refocus the complex information of the object reconstructed by Eq. (4) via angular spectrum or Fresnel transform [4]. For simplicity, the proposed algorithm has been validated under IP conditions. In future work, we will expand it to defocused holograms and non-telecentric systems.

 

Fig. 1. Optical configuration of a Mach–Zehnder-based DHM system for imaging biological samples. BS, beam splitter; CL, converging lens; IP, image plane; M, mirror, MO, microscope objective; TL, tube lens.

Download Full Size | PDF

3. AUTOMATIC ALGORITHM TO RETRIEVE COMPLEX OBJECT INFORMATION

This section is devoted to describing our approach. If the hologram is recorded onto the surface of a discrete sensor with $M \times N$ square pixels of the ${\Delta _{xy}}$ side, the discrete digital reference wave ${r_D}(m,n)$ is expressed as

 

Fig. 2. Representation of the center position of the $+1$ term, $({{u_{{\max}}},{v_{{\max}}}})$. (a) The frequential position of the maximum peak coincides with an integer value, and (b) the frequential position of the maximum peak does not coincide with an integer position. The latter case is the most common in experimental DHM systems. The green circle denotes the position of the maximum peak to compute the digital reference wave and reconstruct accurate phase images.

Download Full Size | PDF

In Ref. [22], Trujillo et al. estimated the optimized values of $ u_{\max } $ and $ v_{\max } $ using an iterative approach based on double nested loops. In particular, Trujillo’s method was based on a mean-thresholding-and-intensity-summation metric. For each value of $ u_{\max } $ and $ v_{\max } $, including non-integer values, the authors compute the digital conjugated reference wave, reconstruct the phase image, generate the binary phase image using a fixed thresholding value (i.e., binary phase image), and sum all pixels in each binary image. For the optimal values of $ u_{\max } $ and $ v_{\max } $, no phase perturbations should appear in the reconstructed phase image (or minima phase perturbations in the presence of high-order aberrations such as astigmatism and coma). Therefore, its binary image should be white (i.e., all pixels in the binary phase image should be 1). Nonetheless, if the values of $ u_{\max } $ and $ v_{\max } $ are not the optimal ones, the reconstructed phase image is distorted by a residual plane tilt and, therefore, its binary image is a combination of black and white pixels. The higher the number of black pixels in the binary phase image, the higher the residual plane tilt. Figures 3(a)–3(c) show the reconstructed phase images of a section of the head of a Drosophila melanogaster fly using different values of the position of the maximum peak, $({{u_{\max}},{v_{\max}}})$, to generate the digital reference wave. Whereas panels (a) and (b) illustrate two non-optimal compensations, corresponding to erroneous positions of $({{u_{\max}},{v_{\max}}})$, Fig. 3(c) shows the optimal compensation of the reference wave, leading to the best-reconstructed phase image. The values of $ u_{\max } $, and $ v_{\max } $ are reported in the top left corner of the reconstructed phase images. As Fig. 3(f) shows, the best-reconstructed phase image should provide a binary image with a very minimum or no dark pixels at all. In other words, reconstructed phase images with perturbations (i.e., non-optimal phase reconstructions) have a higher number of dark pixels than the best-reconstructed phase image with no phase perturbations. These dark pixels are related to the off-center in the spectrum of the non-optimal phase reconstructions.

 

Fig. 3. Validation of the proposed method based on a minimization of a cost function. Panels (a)–(c) in the top row show reconstructed phase images of a section of the head of a Drosophila melanogaster fly for different values of ${u_{\max }}$ and ${v_{\max }}$. In the bottom row, panels (d)–(f), binary phase images are shown.

Download Full Size | PDF

 

Fig. 4. Block diagram of the proposed method.

Download Full Size | PDF

The difference between Ref. [22] and this work is the implementation of the search for the optimal values of $ u_{\max } $, and $ v_{\max } $. Trujillo et al. [22] implemented their approach using nested loops. The authors found the values of $ u_{\max } $, and $ v_{\max } $ by searching the maximum number of pixels equal to 1 in the binary image. Conversely, in this work, our implementation is based on unconstrained nonlinear optimization algorithms, which are intrinsically faster than nested loops. We find the values of ${u_{\max }}$, and ${v_{\max }}$ by searching the minimum number of pixels equal to 0 (e.g., black pixels) in the binary image. Therefore, the proposed cost function tracks the number of black pixels,

The block diagram of the proposed algorithm is shown in Fig. 4. The block diagram comprises four main stages: input parameters, spatial filtering, compensation of the tilt of the plane reference wave, and estimation of the best-reconstructed phase image. The input parameters of the proposed algorithm are an image-plane hologram recorded in a telecentric DHM system ($ h $), the source’s wavelength ($ \lambda $), and the pixel size ($\Delta xy$). The second stage is the spatial filtering of the $+1$ order. This filtering can be performed by estimating the integer position of the maximum peak of the $+1$ term in the hologram spectrum, $(u_{{\max}}^0,v_{{\max}}^0)$, and applying a circular mask. The dimensions of the circular mask filter depend on the magnification and numerical aperture of the DHM system [35]. Nonetheless, without knowing this information, the radius of the circular mask can be assumed to be equal to $1/3\sqrt {|{u_0} - u_{{\max}}^0{|^2} + |{v_0} - v_{{\max}}^0{|^2}}$, $({{u_0},{v_0}})$ being the integer position of the DC term. The next step is the compensation of the tilt of the plane reference wave. This stage is split into two parts: the setting of the cost function and its minimization of using the MATLAB built-in function fminunc. This function finds a local minimum of an unconstrained multivariable function using the quasi-Newton or trust-region algorithm. The scope of this research work is the investigation and validation of an automatic reconstruction approach for DHM that requires the minimum number of user inputs (e.g., hologram, wavelength, and pixel size). Because the trust-region algorithm requires a gradient function as the input of the minimization process, we have implemented the quasi-Newton algorithm. This work does not provide a rigorous study of the minimization solver and algorithm. Note that the selection of a different solver and algorithm, and the use of parallel computing may yield reduced response time; however, that investigation is beyond the scope of this paper.

Table 1. Comparison of Nested-Loops (NL), Centroid (C), and Proposed Methods

View Table

The input parameters of the fminunc function are the cost function [Eq. (8)], the integer position of the maximum peak in the filtered hologram spectrum, $(u_{{\max}}^0,v_{{\max}}^0)$, and a set of optimizations options. The optimization options are the maximum number of iterations allowed and a termination tolerance on the cost function. The tolerance has been set to ${{10}^{- 3}}$ since a lower tolerance value does not improve results. Regarding the number of iterations, the steps of the proposed method are repeated until convergence. The convergence is defined when the norm of the gradient of the cost function is smaller than the tolerance. For the tested holograms, convergence is achieved within the first five iterations. The number of iterations is highlighted dependent on the initial value in the minimization approach. Since the hologram has been recorded in a telecentric-based DHM system and the reference wave is a tilted plane wave, the center integer position of the $+1$ term, $(u_{{\max}}^0,v_{{\max}}^0)$, is the best candidate to start the search of the optimal non-integer spatial frequency $({{u_{\max}},{v_{\max}}})$ since the optimal values should be close to the integer position of the $+1$ term. Nonetheless, in the event of no achieving convergence, we have set 30 iterations as the maximum number of iterations allowed before the algorithm breaks the execution.

The estimation of the cost function [Eq. (8)] requires the binarization of the reconstructed phase image that depends significantly on the digital reference wave via the position $({u_{\max}},{v_{\max}})$. The binary image is obtained by applying the MATLAB built-in function imbinarize to the reconstructed phase image. The imbinarize function uses Otsu’s method, and its input parameter is the reconstructed phase image. In this study, we have used the default threshold value of Otsu’s method since it provides the best-reconstructed phase image with minimum phase perturbations. For each iteration of the minimization process, the algorithm finds a pair of possible values $(u_{{\max}},v_{{\max}})$, generates the corresponding digital reference wave, and computes the reconstructed phase image reconstruction. This resulting phase image is binarized, and the cost function J is computed. The proposed algorithm is an iterative method, starting with the integer position of the maximum peak in the $+1$ term, $(u_{{\max}}^0,v_{{\max}}^0)$ until finding the optimal non-integer position given by the $(u_{{\max}}^0,v_{{\max}}^0)$. Once the optimal values of $({{u_{\max}},{v_{\max}}})$ are estimated, the final phase image can be reconstructed without phase aberrations due to the interference tilt. It is important to remember that the amplitude image can be estimated after the spatial-filtering step.

For comparison purposes, we have shown the performance of the proposed method versus the one from nested loops [22] and centroid-based [19] algorithms. Table 1 shows the estimated values of $ u_{\max } $ and $ v_{\max } $, the accuracy of reconstructed phase values, and the processing time for three different optically thin samples. For each method, the reconstructed phase images of the phase USAF target of Benchmark Technologies, a transverse section of the head of a Drosophila melanogaster fly, and glioblastoma cells are shown in Fig. 5. For the phase USAF target, we have measured the average phase and its standard deviation of the elements in Group 6 to be equal to (${2.00} \pm {1.90}$) rad for the nested loops algorithm, (${1.32} \pm {0.1}$ 7) rad for the centroid-based algorithm, and (${1.83} \pm {0.16}$) rad for the proposed method. Regarding manufacturer specification, the nominal phase value of these elements is 1.84 rad ($= {2}\pi ({{n - 1}}){t}/\lambda$) for a wavelength of $\lambda = {532}\;{\rm nm}$, thickness of ${\rm t} = {300}\;{\rm nm}$, and refractive index of ${\rm n} = {1.52}$. The error difference between the experimental and nominal phase values is 0.54%, verifying the high accuracy of the proposed method to quantify phase measurements. The error difference for the nested loops and centroid-based algorithms is 8.7% and 28.3%, respectively, being higher due to the wrapped phase values. It is important to remember that, in biological imaging, phase measurements enable the estimation of biological parameters such as the integral intracellular refractive index [36]. Thereby, accurate phase measurements are imperative since variations in phase values are used as a diagnostic and measuring tool in biological research.

 

Fig. 5. Reconstructed phase images of the phase USAF target of Benchmark Technologies (first row), a section of the Drosophila melanogaster fly (second row), and glioblastoma cells (third row) for three different algorithms: nested-loops-based algorithm (first column), centroid-based algorithm (second column), and the proposed method (last column).

Download Full Size | PDF

Alternatively, we have estimated the performance of each compensation method by binarizing its reconstructed phase image (first row in Fig. 5) and counting the ratio between the total number of white pixels versus the total number of pixels in the reconstructed phase image. It is important to remember that the best reconstructed phase image’s binarized phase image has a very minimal or no dark pixels at all (i.e., the majority of the pixels in the binary image are white pixels). For example, in the USAF target, the image has ${1024} \times {1024} = {1.408} \times {{10}^6}$, and the binarized phase image has ${9.96} \times {{10}^5}$ white pixels, resulting in 95% of the total image (i.e., the performance of the method is 0.95). Table 1 compares the performance of each technique based on quantifying the quality of the binary image. Based on these values, the reconstructed phase images estimated by our proposed method provide the best reconstructed phase images without residual phases. This high performance is achieved because our proposed method achieves high precision in calculating the position of $ u_{\max } $ and $ v_{\max } $. We have analyzed the accuracy of the reconstructed phase images based on the precision of the values of $({{u_{\max}},{v_{\max}}})$, and measuring the mean square error (MSE) between them. Based on the MSE values, the highest the precision of the $({{u_{\max}},{v_{\max}}})$ values yields the minimum the error in the reconstructed phase image. Phase measurements are strongly dependent on ambient fluctuations and the experimental implementation. The maximum number of the decimal digits in the $({{u_{\max}},{v_{\max}}})$ values should be the one that generates a trustable reference wave and provides accurate phase measurements for any DHM implementation. Among the different DHM systems, common-path DHM systems are known to provide the least phase error, which is in the order of 0.0003 rad [37]. Based on this phase error and the estimated MSE values (not shown here), we need at least three decimal digits in the values of $ u_{\max } $, $ v_{\max } $ to provide accurate phase measurements. The processing time is reported based on a Windows-based i7-8700 K CPU (3.70 GHz) 16.0 Gbyte RAM desktop computer. The average processing time of our method is 1.55 s for holograms of ${1024} \times {1024}$ pixels, and 2.95 s for a hologram of ${2048} \times {2048}\;{\rm pixels}$. Note that the processing time has been reduced by a factor of ${2.3} \times$ compared with the centroid-based algorithm [19]. Regarding the nested-loops-based algorithm [22], which is the slowest algorithm in this comparison, our method is around ${40} \times {{- 90}} \times$ faster than the nested-loops-based method. The proposed method provides precise values of $ u_{\max } $ and $ v_{\max } $, and precise phase measurements without phase perturbations.

4. CONCLUSIONS

This work presents an automatic and fast method to reconstruct the quantitative phase distribution of unstained biological samples with minimal or no phase perturbation. The proposed method compensates for the tilt between the reference and object waves in off-axis DHM systems operating in telecentric configuration. The input parameters in our method are an in-focus hologram, the wavelength of the source, and the pixel size. The limitations of the proposed technique are: (1) one cannot reconstruct phase images of defocused holograms, (2) the approach only works for off-axis DHM systems operating in the telecentric regime, and (3) there is no compensation of high-order aberrations. Compensation of high-order aberrations such as astigmatism and coma is required to ensure accurate quantitative phase analysis. Future work will address these limitations to provide a more general reconstruction approach to ensure accurate quantitative phase analysis applicable to any off-axis DHM system. In addition, future work will investigate the applicability of the technique to optically thick samples. Among the hallmarks of the proposed approach is the high accuracy in estimating the parameters of the digital reference wave. The interference angle is calculated precisely, within three decimal places, without compromising the computation time. The proposed method performs ${40} \times$ faster than a previously reported automatic approach based on nested-loops and 2.3× faster than the centroid-based algorithm. The experimental results show that the proposed approach provides quantitative phase images without phase perturbations. Although the proposed method has been validated for transparent biological samples, it can also be applied for reflective-based DHM systems. A ready-to-use version of the MATLAB GUI plugin, the source code, example data sets, and a short user manual can be found online [38].