Standardization Techniques for Single-Shot Digital Holographic Microscopy

Abstract

Digital holographic microscopy (DHM) is a mature technology for quantitative phase imaging. Thousands of articles have been published on this topic over the last couple of decades. Our goal in this article is to emphasize that single-shot holographic microscopy systems offer several practical advantages and in principle capture the full diffraction-limited information of interest. Since phase cannot be measured directly, phase reconstruction is inherently a computational problem. In this context, we describe some traditional algorithmic ideas as well as newer sparse optimization-based methodologies for phase reconstruction from single-shot holograms. Robust operation of a DHM system additionally requires a number of auxiliary algorithms associated with fractional fringe detection, phase unwrapping, detection of focus plane, etc., that will be discussed in some detail. With the data-driven nature of applications of DHM being developed currently, the standardization or benchmarking of algorithmic ideas for DHM systems is important so that same sample imaged by different DHM systems provides the same numerical phase maps. Such uniformity is also key to establishing effective communication between DHM developers and potential users and thereby increasing the reach of the DHM technology.

1. Introduction

Digital holography (DH) is an interferometric imaging modality [1, 2, 3, 4] which enables recording and numerical reconstruction of phase map of a coherent light beam that has interacted with the object of interest. In DH, the interference pattern is recorded on an array sensor. As a result, unlike traditional film-based holography, the interference fringe data is now available as a matrix of numbers. This numerical nature of the data enables one to apply a number of algorithmic ideas for quantitative phase estimation and 3D image formation. While film-based holograms were typically reconstructed by re-illumination with the conjugate reference beam, DH does not necessarily have to mimic this physical hologram replay process. As we will see in this article, novel algorithmic ideas can in principle take DH beyond film based holographic imaging performance. Digital form of holography was first demonstrated in early work [5] where recording of holograms was performed using a vidicon tube array detector and a fast Fourier transform (FFT)-based numerical reconstruction was demonstrated. The DH modality gained much popularity starting early 1990s when charge-coupled device (CCD) and complementary metal oxide semiconductor (CMOS) array sensors became readily available with simultaneous increase in computational power at low cost. Since then, DH has has received a lot of attention in the Science and Engineering literature. A quick web search for the term “digital holography” yields over 150,000 documents, and interestingly almost half of them are related to “digital holographic microscopy” (DHM) suggesting the attention received by DHM as an imaging modality. Compared to intensity or amplitude-based microscopic imaging modalities commonly in use in Life Sciences, quantitative phase is known to be more sensitive modality and provides a much more complete information about a sample object under investigation. Quantitative phase imaging further allows the possibility of label-free live cell imaging. The minimal wet-lab processing required for label-free imaging also makes quantitative phase a potentially cost-effective imaging modality. A large number of publications on DHM continue to appear (mainly in the physics, optics, and engineering journals); however, it is still not used widely as a preferred imaging modality in the Life Sciences community. For example, if one walks into a typical bioscience research laboratory in a university setting or a pathology clinic in a hospital, the researchers and technicians there are most likely not using quantitative phase as an imaging modality in their work. While DHM can be considered as an emerging modality at present, its widespread use in future will require standardization and benchmarking of algorithmic methodologies as well as systems. In particular with a number of system configurations and algorithmic combinations being used, it is not clear if we will get identical numerical phase maps if the same sample is imaged by two different DHM systems. This aspect is very important in the light of current trend of data-driven imaging application development [6]. Otherwise, individual applications will likely become system or algorithm-specific, thus limiting their usage. This article summarizes a number of ideas from our work over the last several years that are relevant to standardization of algorithms for single-shot DHM systems that record image plane digital holograms.

The organization of this article is as follows. In Section 2, we discuss the importance of employing a single-shot DHM system, first from the perspective of simplicity of deploying such a system in practical settings. We additionally argue that the single-shot holographic image record actually contains all the diffraction-limited information about the object field that is of interest. In Section 3, we present algorithmic methodologies that may be used to demodulate a single-shot hologram record with emphasis on the sparse-optimization-based approaches and their noise and resolution advantage. Beyond the core requirement of good-quality phase reconstruction, additional auxiliary methodologies are required to make the DHM imaging performance robust. In Section 4, we discuss a few algorithmic methods for this purpose that make a DHM instrument more user-friendly. Section 5 describes some interesting issues related to computational 3D image formation from holographic data which is an important research topic in itself. We conclude the discussion in Section 6 by providing further insights gained from our experience of working with clinicians and emphasize the need for a DHM user consortium which will help popularization of DHM technologies within the Life Sciences community.

2. Single-shot DHM systems: motivation

A DHM system in bare-bones form is shown in Figure 1. The system as shown here mainly comprises of a coherent plane wave illuminating a sample of interest, which is then imaged using an afocal imaging system consisting of a microscope objective (MO) and a tube lens. The array sensor is placed in the image plane so that an image field $O\left(x,y\right)$ corresponding to the magnified version of the exit wave in the sample plane is incident on the sensor. Since DHM is an interferometric modality, a part of the illuminating beam is split beforehand and is recombined at the sensor. The reference beam field at the sensor may be denoted as $R\left(x,y\right)$. If sufficient temporal coherence is available and the path lengths in the interferometer are matched well, an interference pattern with good contrast may be recorded on the sensor, which may be denoted as:

Figure 1.

Single-shot digital holographic microscopy (DHM) in a bare-bones form. The basic imaging system is an afocal infinity-corrected microscope which forms image field $O\left(x,y\right)$ at the sensor plane containing full diffraction-limited information as permitted by the microscope objective (MO). The reference beam $R\left(x,y\right)$ is derived from the same illuminating source. The signal detected at the sensor consists of off-axis interference fringes.

$$H\left(x,y\right)={\left|R(x,y)+O(x,y)\right|}^2.$$E1

If future DHM systems are to be deployed in field, e.g. at a clinic, one may not always have access to research laboratory like environment with vibration isolation platforms. Luckily, the CCD or CMOS arrays available off-the-shelf at reasonable costs today are quite sensitive. As a result, a DHM system employing a few mW laser source requires just a few milliseconds of exposure time in order to record good-quality interference fringes with high contrast. Such short exposure times (much shorter compared to the timescale of typical ambient vibration) imply that vibration isolation is not required for recording a single frame of hologram data. The configuration in Figure 1 further suggests that the image field $O\left(x,y\right)$ contains full diffraction-limited information that one may wish to recover. So in principle, a single recorded data frame $H\left(x,y\right)$ contains interference fringes with good contrast and also embeds full diffraction-limited information in it. Multi-shot DHM systems, e.g. employing the phase shifting technique, on the other hand require stringent vibration isolation platforms. Additionally, there is a need of deploying piezo-transducers for introducing three or four phase shits in a repeatable manner for generation of a single phase image. The vibration isolation and phase shifting hardware can increase the system costs beyond what is affordable in a place like a primary healthcare center. The necessity of multiple interference records for a phase-shifting-based DHM system also suggests that it cannot perform phase imaging of dynamic samples as may be of interest to researchers engaged in live cell imaging. As described in the next section, traditionally the discussion regarding use of single-shot vs. multi-shot (or phase shifting) DHM systems has concentrated on achievable image resolution in the two methodologies. However, the possibility of full diffraction-limited imaging performance in single-shot DHM records via sparse optimization algorithms as we will present here addresses this debate to some extent. We can now therefore envision DHM systems that offer all practical advantages of a single-shot operation and at the same time do not suffer in terms of loss of resolution or phase accuracy.

3. Single-shot phase reconstruction algorithms

Having made a case for single-shot operation of DHM systems, we now turn to DHM configurations and corresponding algorithms for phase reconstruction. In recording of a digital hologram, two geometrical configurations are typically encountered. The off-axis configuration involves a plane reference beam $R\left(x,y\right)$ which nominally makes an angle with respect to the object beam $O\left(x,y\right)$ as shown in Figure 1. The second one is the in-line hologram recording configuration where the reference and object beam are nominally colinear. The in-line configuration is typically employed when we wish to image weak scattering objects like particulates or microorganisms floating in a fluid medium with a DHM system. Colinear nature of such systems makes them more compact which in itself offers several practical advantages. When morphology of more complex structures is to be imaged, the weak scattering approximation may not be fully valid, and in such cases, the off-axis configuration is preferred. The off-axis configuration has received a lot of attention in the literature. Let us for instance assume a reference beam for the off-axis case to be of the form: $R\left(x,y\right)=R_{0}exp\left(i2\pi f_{0x}x\right)$. The interference signal as in Eq. (1) can then be written more explicitly as:

$$H\left(x,y\right)={\left|O\right|}^2+{\left|R_0\right|}^2+2\mid O\Vert R_0\mid cos\left[2\pi f_{0x}x-{\varphi }_O\left(x,y\right)\right].$$E2

In the aforementioned relation, ${\varphi }_O\left(x,y\right)$ represents the object beam phase. The special form of the reference beam means that the Fourier transform structure of the hologram signal is quite simple. In particular, the first two terms of the hologram are located at the center of the Fourier space and the cosine term leads to two lobes in Fourier space centered on the carrier fringe frequency as shown in Figure 2. Note that if the spatial bandwidth of the function $O\left(x,y\right)$ is given by $2B_x$ in the $x-$ direction, the bandwidth of the ${\left|O\right|}^2$ is $4B_x$. The bandwidth of ${\left|R_0\right|}^2$ is practically well within the bandwidth of ${\left|O\right|}^2$ as it most contains some low-frequency variations. In the traditional film-based holography, the replay of the off-axis hologram would involve re-illuminating the developed film by conjugate reference beam. The digital processing equivalent of this procedure if the filter out the cross-term represented by the side lobe, place it in the center of the Fourier space, and evaluate inverse Fourier transform in order to get an estimation of the object beam $O\left(x,y\right)$ at the array sensor plane [7]. As seen from Figure 2, we however note that in order to recover the full information about the object wave, the three terms cannot overlap in the Fourier space. If the hologram can be considered to be sampled well on the digital array sensor, then the non-overlap requirement on the three energy lobes in Fourier space is highly restrictive as it implies that the allowable bandwidth of the object wave is much lower compared to the bandwidth of the array sensor itself. In other words, the filtering methodology inherently leads to lower-resolution object wave recovery, despite the fact that a higher resolution information is present in the hologram data. With the hologram being sampled on digital array sensors with typical pixel sizes of 2–5 μm already restricts the reference beam angle in order to get well-sampled fringes on the sensor and the non-overlap requirement on Fourier bands further worsens the situation. A relaxation in the overlap condition to enable higher resolution is achieved by two notable approaches [8, 9]. If the afocal microscope system as in Figure 1 was used with typical incoherent bright-field illumination (e.g. using LED or lamp), the user is expected to get full diffraction-limited image resolution. However, as soon as one switches to the digital holography mode with coherent illumination and introduction of off-axis reference beam, the recovered resolution is worse compared to the diffraction-limited resolution when linear filtering methodologies are used. This situation is highlighted in image plane holography as poor resolution implies loss of edge information which is critical to visual perception of the image field by users. We emphasize here that the loss of image resolution as discussed here is an artifact of the Fourier filtering-based algorithmic framework used for object field recovery which prevents us to recover the full pixel resolution. The single-shot object field recovery problem therefore needs to be examined further. We observe that the recorded hologram represented by Eq. (2) is a quadratic equation in terms of the amplitude $\mid O\mid$ of the object beam. A solution for the object beam amplitude can thus be expressed as:

Figure 2.

Fourier space representation of the off-axis hologram. Note that if the hologram is processed by filtering of the cross-term, then the allowable bandwidth for object wave is much lower compared to the bandwidth of the detector array.

$$\mid O\mid =-\mid R_0\mid cos\left(2\pi f_{0x}x-{\varphi }_O\right)\pm \sqrt{H-{\left|R_0\right|}^2{sin}^2\left(2\pi f_{0x}x-{\varphi }_O\right)}.$$E3

This relation is curious as it suggests that the magnitude $\mid O\mid$ explicitly depends on the phase ${\varphi }_O$ of the object beam itself. One may now choose any arbitrary phase map ${\varphi }_O$, compute $\mid O\mid$ as per Eq. (3), and the corresponding complex-valued object field will exactly satisfy the hologram data represented by $H\left(x,y\right)$. The solution of the single-shot holography problem is therefore ambiguous when we consider the problem in this numerical form. More specifically, since a single hologram frame $H\left(x,y\right)$ is to be used to estimate the complex-valued object wave function $O\left(x,y\right)$ (assuming the prior knowledge of the reference beam $R\left(x,y\right)$), the numerical problem of estimating $O\left(x,y\right)$ from a single hologram frame $H\left(x,y\right)$ can be considered to be an incomplete data problem. The particular form of the reference beam played no role in this incomplete data argument, and so this solution ambiguity exists for in-line holography case as well. In the in-line case, various terms of the hologram signal actually overlap completely in the Fourier space, and there is no possibility of obtaining a solution for $O\left(x,y\right)$ by means of a simple linear filtering operation. The off-axis case at least provides an approximate lower-resolution estimate of the unknown object wave $O\left(x,y\right)$ due to the particular Fourier structure of the hologram signal. The size of Fourier space window to be used to filter out the cross-term, however, still remains a subjective choice. In Figure 3(a), we show an off-axis image plane hologram for a red blood cell (RBC) sample. The recovered phase map of the object beam is shown in Figure 3(b)–(d) for three choices of filter windows of sizes of sizes $0.2{\rho }_0$, $0.5{\rho }_0$, and $0.7{\rho }_0$, respectively, where ${\rho }_0$ refers to the distance between the zero frequency and the carrier frequency peak positions.

Figure 3.

(a) Image plane digital hologram of RBCs. (b)–(d) Phase of object beam determined by Fourier filtering of the cross-term with filter window sizes $0.2{\rho }_0$, $0.5{\rho }_0$, and $0.7{\rho }_0$, where ${\rho }_0$ is the distance between the zero frequency and the carrier frequency peak positions. The variability of resolution and background artifacts makes this methodology highly subjective.

The three phase reconstructions differ in terms of their resolution and background artifacts which may have arisen due to the contribution of the dc or zero-frequency terms to the reconstruction. The choice of filter sizes is generally left to the user; however, this may lead to variability of numerical phase maps between different users which will generally not be acceptable to the user community. Overall, the aforementioned discussion suggests that the phase reconstruction problem for single-shot image plane holograms is not quite simple. A solution maybe designed for one setup for a particular sample by a given group of researchers. But translating this solution to a robust system to be used by third-party users is a nontrivial proposition. Solving such variability issues and obtaining a full pixel resolution seem to be possible with the phase shifting methodology which is however a multi-shot approach not suitable for practical deployment of DHM systems. Since the single-shot phase reconstruction is an incomplete data problem, the use of image sparsity ideas popularized by the compressive sensing community can be handy. This, however, requires one to change the reconstruction framework to an optimization procedure as we will explain in detail in the following subsection.

3.1 Single-shot phase reconstruction as an optimization problem

It is now a widely accepted fact in signal and image processing community that natural images have number of degrees of freedom that are much smaller than the number of visual pixels used to represent them. In fact, a number of image and video compression standards like JPEG or MPEG regularly exploit this redundancy in natural images that arise due to structures of the objects present in the typical images. The same ideas must apply to phase reconstruction as well. In particular, when observing objects like biological cells using a DHM, one should be able to exploit the representational redundancy of the desired complex-valued field $O\left(x,y\right)$ at the sensor plane corresponding to the objects being imaged. As explained before, the Fourier filtering solution for the off-axis holography essentially has origins in the way holograms used to be replayed in film-based holography. A different framework for complex object wave recovery is, however, possible in digital holography as the interference record is now available in the numerical form, and there is no need to mimic film-based holography in the numerical processing of this numerical data. In the following discussion, we describe an optimization framework [10] for recovering the object wave $O\left(x,y\right)$ from an image plane hologram $H\left(x,y\right)$ that may help in addressing the difficulties associated with the single-shot phase reconstruction problem. In particular, we seek to minimize the functional:

$$\begin{array}{c}C\left(O,O^{\ast }\right)=C_1+C_2\\ ={‖H-\left({\left|R\right|}^2+{\left|O\right|}^2+R^{\ast }O+{\mathit{RO}}^{\ast }\right)‖}^2+\alpha \psi \left(O,O^{\ast }\right).\end{array}$$E4

Here, ${‖\dots ‖}^2$ denotes the squared L2-norm of the quantity inside the norm. The first term $C_1$ of the cost function is, therefore, a least-square data fit to the interference model. The second term refers to a constraint or penalty term which encourages solution with some desired properties. The positive constant $\alpha$ decides the weight between the two terms of the cost function. Note that the overall cost function in Eq. (4) is real- and positive-valued, whereas the solution we are seeking is complex-valued. In such cases, the steepest descent direction needed for the purpose of iterative optimization algorithms may be evaluated in terms of the complex or Wirtinger derivatives. The Wirtinger derivative is evaluated with respect to the variable $O^{\ast }$ and makes the implementation almost as simple as a real-valued optimization problem. The Wirtinger derivative for the first term of the cost function is straightforward to evaluate and is given by:

$${\nabla }_{O^{\ast }}{C}_1=-2\left[H-{\left|R+O\right|}^2\right]\cdot \left(R+O\right),$$E5

A few suitable choices for the penalty term, the corresponding expressions for the Wirtinger derivatives, and the characteristic properties of the penalty functions are provided in Table 1. The first penalty is the squared gradient sum which encourages locally smooth solutions. The second penalty function is the total variation (TV) penalty which is known to be edge-preserving and encourages piece-wise constant solutions. The third is the modified Huber penalty which locally acts like the squared gradient sum or TV penalties depending on the parameter $\delta$. This parameter can be selected based on the statistics of gradients in a given image. For example, $\delta$ may be made proportional to the median of the gradient magnitudes. The Huber penalty will then act like the edge-preserving TV penalty at pixels where $\mid \nabla O\mid >>\delta$ and act like the smoothing quadratic penalty for pixels where $\mid \nabla O\mid <<\delta$. The expressions for the Wirtinger derivatives for various choices of the penalty function are provided in the second column of Table 1 for convenience of the readers. Implementation of the optimization solution can proceed with any of the gradient-based iterative schemes, the simplest being the gradient descent scheme. The $\left(n+1\right)-$th iteration of this scheme may be described as:

Penalty term	Wirtinger derivative	Remarks
${\sum }_k{\left|\nabla O_k\right|}^2$	$-{\left[{\nabla }^{2}O\right]}_k$	Smoothness penalty
${\sum }_k\mid \nabla O_k\mid$	$-{\left[\nabla \cdot \left(\frac{\nabla O}{\mid \nabla O\mid }\right)\right]}_k$	Total variation (edge-preserving)
${\sum }_k\left[\sqrt{1+\frac{{\left|\nabla O_k\right|}^2}{{\delta }^2}}-1\right]$	$-{\left[\nabla \cdot \left(\frac{\nabla O}{\sqrt{1+\frac{{\left|\nabla O\right|}^2}{{\delta }^2}}}\right)\right]}_k$	Huber penalty (locally adaptive)

Table 1.

Suitable of penalty functions for optimization-based image recovery in single-shot digital holography along with the Wirtinger derivatives.

$$O_{n+1}=O_n-\tau {\left[{\nabla }_{O^{\ast }}C\right]}_{O_n}.$$E6

The parameter $\tau$ represents a step size which may be derived using backtracking line search [11]. The gradient descent scheme can be replaced by alternative iterative methods such as conjugate gradient or Nesterov accelerated gradient. While several algorithmic choices are available, one of the important practical problems is that the solutions obtained by this procedure depend on the regularization parameter $\alpha$. Tuning of $\alpha$ for individual hologram data sets can become tedious. In Figure 4, we show object wave recovery for an image plane hologram corresponding to a step object. The TV penalty has been used in this case. It is observed that the quality of phase solution can change significantly depending on the value of $\alpha$. For a low numerical value of $\alpha$, the phase reconstruction shows fringe-like artifacts. On the other hand for a large numerical value of $\alpha$, the phase solution shows over-smoothing. A “good” value of $\alpha$ proves that, in principle, an excellent recovery of the edge (and thus full pixel resolution) is possible via the optimization approach. A methodology that does not require any free parameter like $\alpha$ is thus practically very useful. The optimization problem for this purpose may be restated as follows [12]. We wish to determine a solution $O\left(x,y\right)$ satisfying

Figure 4.

Phase recoveries from a single-shot hologram corresponding to a step phase object. The optimization method with TV penalty is used, and the solution is shown for three values of the free parameter $\alpha =\mathrm{0.1,1,10}$ in the cost function in Eq. (4).

$${‖H-{\left|R+O\right|}^2‖}_2^2<\epsilon ,$$E7

such that among all solutions satisfying the aforementioned criterion, the desired solution has the lowest numerical value of the penalty function. This objective may be achieved by progressing the iterative solution in the following manner:

• At the beginning of (n + 1)th iteration, the solution $O_n$ is updated so as to reduce the error objective $C_1$.

• The resulting intermediate solution may be denoted as $O_{n,\mathit{int}}$.

• The TV (or other penalty) for this iterative solution may be reduced with recursive steepest descent iterations of the form:

$$O_{n,k+1}=O_{n,k}-\tau {\left[{\nabla }_{O^{\ast }}\psi \left(O,O^{\ast }\right)\right]}_{O_{n,k}}.$$E8

At the end of $N_{\mathit{TV}}$ sub-iterations of this form starting with $O_{n,0}=O_{n,\mathit{int}}$, we set $O_{n+1}=O_{n,N_{\mathit{TV}}}$.

• Further, the reduction of $C_1$ and $\psi$ is performed in an adaptive such that the distances $d_1={‖O_n-O_{n,\mathit{int}}‖}_2$ and $d_2={‖O_{n+1}-O_{n,\mathit{int}}‖}_2$ are made approximately equal.

The aforementioned strategy is inspired by the Adaptive Steepest Descent Projection onto Convex sets (ASD-POCS) algorithm [13]. The main thought process behind the scheme aforementioned is that the solution is driven away from minima in $C_1$, where the numerical value of the penalty function $C_2$ may be large. Eventually, the solution is driven to an equilibrium point where both the reduction in $C_1$ and $C_2$ oppose each other. While the scheme aforementioned needs multiple sub-iterations in each iteration, another scheme that naturally achieves the goal of $d_1\approx d_2$ is the mean gradient descent (MGD) [14, 15]. In MGD, the solution is moved along the direction $\widehat{u}$ that bisects the steepest descent directions corresponding to the two objectives $C_1$ and $C_2$. The solution therefore proceeds as:

$$O_{n+1}=O_n-\tau \left[{\widehat{u}}_{O_n}\right],$$E9

where,

$$\widehat{u}=\frac{{\widehat{u}}_1+{\widehat{u}}_2}{2},$$E10

and ${\widehat{u}}_{1,2}$ are the unit vectors defined as:

$${\widehat{u}}_{1,2}=\frac{{\nabla }_{O^{\ast }}{C}_{1,2}}{{‖{\nabla }_{O^{\ast }}{C}_{1,2}‖}_2}.$$E11

Moving in the bisector direction typically increases the angle between the unit vectors ${\widehat{u}}_1$ and ${\widehat{u}}_2$ till the angle settles to a large obtuse value. At this point, simultaneous reduction of both the objectives $C_1$ and $C_2$ is not possible, and the solution thus reaches a balance point. Nominally, the two schemes aforementioned lead to a similar solution. Note that the goal of the optimization here is not to achieve a minimum of an overall cost function but rather to achieve an equilibrium between various objective functions (like $C_1$ and $C_2$) that represent the desirable properties of the solution $O\left(x,y\right)$. Figure 5 shows two illustrations of the adaptive optimization process applied to the single-shot image plane hologram data corresponding to a step phase object. The recovery of phase step is excellent independent of the off-axis or on-axis configuration. As long as the reference beam is known, the two cases are operationally the same for the optimization algorithm. The holograms in Figure 5 have been simulated with Poisson noise, and the RMS error achieved for the phase solution using the optimization procedure is observed to be better than the single-pixel-based shot noise limit. An experimental demonstration of this sub-shot noise phase recovery via the optimization procedure was also shown earlier with low-light-level interferograms. In particular, it was demonstrated that for interferograms with light level of the order of $10$ photons/pixel, the RMS phase accuracy for a simple lens phase object was 5× better that what is expected from the shot noise limit [16]. Achieving such phase accuracy with classical light interferometry is not surprising as we are implicitly constraining the solutions by making a suitable choice of the penalty function. The important message here is that optimization-based phase recovery is expected to provide phase measurements accuracy better than array sensor’s noise floor. The single-shot full pixel resolution capability along with improved phase accuracy makes this technique suitable for employing with a quantitative phase microscopy system. In the illustration in Figure 6, we demonstrate optimization-based phase recovery using an image plane hologram of cervical cell nucleus. The cervical cell sample was obtained from a clinical collaborator (All India Institute of Medical Sciences) and is a typical pap-smear slide sample. Phase solutions using Fourier filtering as well as optimization method are displayed. The Fourier filtering solution used a circular filter window of radius $0.6$ times the distance between the dc and cross-term peaks. An interesting aspect of the optimization procedure is that the iterations are fully in the image domain, and as a result, the iterative reconstruction may be performed over a selected region of interest (ROI) [17]. This makes it possible to implement reconstruction of individual cell regions in near real time and also allows for the possibility of parallelizing the phase reconstruction process. The iterative optimization here used the Huber penalty and clearly has a higher resolution compared to the Fourier filtering reconstruction. We emphasize here that the higher resolution information is already present in the image plane hologram and the optimization method that able to extract is much better compared to the Fourier filtering method. Improved resolution in single-shot holography operation offers multiple advantages. First of all, higher resolution enables one to observe fine textural features in a phase reconstruction, which may for example be important in diagnostics. For example, it is well known that cervical cell nuclei have higher “roughness” which show up in the phase profile [18, 19]. Secondly, the single-shot operation makes the DHM system simpler from hardware perspective, thus making it affordable for practical field deployment. The cost reduction here is not due to use of cheaper optics or other hardware components but due to the superior resolution and noise capabilities of the optimization approach to phase recovery which makes such a system possible.

Figure 5.

Adaptive optimization-based phase recoveries for image plane hologram of a step phase object. Both off-axis and on-axis configurations are shown.

Figure 6.

Adaptive optimization-based phase recoveries for image plane hologram of a cervical cell nucleus, (a) bright-field image of cell nucleus, (b) $256\times 256$ ROI of image plane hologram for the cell nucleus, phase reconstruction using (c) Fourier filtering method and (d) optimization method with Huber penalty. Optimization-based recovery clearly shows higher-resolution edge reconstruction. Fourier filtering-based reconstruction was performed with whole camera frame, and then the same ROI in phase reconstruction is cropped in (c).

4. Allied algorithmic requirements for robust phase reconstruction

Beyond the core phase reconstruction algorithm, the making of a robust DHM system requires a few other allied concepts so that the phase reconstructions are numerically repeatable for the same sample across different users and/or DHM systems. We discuss three such ideas in this section that are handy in building a robust DHM system.

4.1 Fractional fringe shift detection

The accurate knowledge of the reference beam $R\left(x,y\right)$ is essential for any phase reconstruction algorithm. While the amplitude $\mid R\left(x,y\right)\mid$ can be determined as by a calibration step by blocking the object beam, estimating the phase of $\mid R\left(x,y\right)\mid$ needs some additional discussion. Any error in phase estimate for $R\left(x,y\right)$ ends up as an additional phase factor in the object phase. If an afocal system is used for image formation in a DHM system, the aberrations are minimized, but the problem of accurate determination of carrier fringe frequency still remains an important problem. In off-axis holography configuration, the fringe frequency is typically estimated by locating the Fourier domain peaks corresponding to the cross-terms in the digital hologram. This step is commonly performed by using the fast Fourier transform (FFT) operation on the recorded off-axis hologram. For a hologram with $N\times N$ pixels, the spatial frequencies in the FFT domain are sampled at discrete intervals of $1/N$ pixel⁻¹. In a given setup, the true carrier fringe frequency peak may be located in between these discrete sample locations in the Fourier domain. If the fringe frequency determination has an error of $\left({\mathrm{\Delta }}_x,{\mathrm{\Delta }}_y\right)$ with $\mid {\mathrm{\Delta }}_x\mid ,\mid {\mathrm{\Delta }}_y\mid \le 0.5$ units in the $x$ and $y$ directions, this leads to a phase error at pixel location $\left(m,n\right)$ that may be given by:

$$\mathrm{\Delta }\varphi \left(m,n\right)=\frac{2\pi }{N}\left({\mathrm{\Delta }}_{x}m+{\mathrm{\Delta }}_{y}n\right).$$E12

This phase error causes a ramp phase background over the desired phase solution. For multiple DHM systems, it is not practically possible to ensure that the fringes will be oriented identically on an array sensor, and this will cause different ramp phase errors for the same sample. As demonstrated in [20], even a small fractional degree camera rotation can lead to different ramp phase backgrounds, making this fractional fringe problem quite sensitive for applications. Detection and elimination of this fractional fringe effects requires a two-step algorithmic approach as we describe here.

The main task at hand here is to get a higher-resolution representation of the cross-term peak in the Fourier domain compared to what is already available from the FFT of the hologram. The fractional fringe detection therefore needs the following steps:

• Compute 2D FFT of the hologram $H\left(x,y\right)$ and locate the integer pixel location $\left(u_0,v_0\right)$ corresponding to the cross-term peak,

• A local discrete Fourier transform computation is then performed as a three matrix product in the following manner in the neighborhood of $\left(u_0,v_0\right)$:

$$\tilde{H}\left(\mathbf{P},\mathbf{Q}\right)=exp\left(-i\frac{2\pi }{N}{\mathbf{PX}}^T\right)H\left(\mathbf{X},\mathbf{Y}\right)exp\left(-i\frac{2\pi }{N}{\mathbf{YQ}}^T\right).$$E13

Here $\mathbf{X},\mathbf{Y}$ denote coordinate vectors ${\left[-N/2,-N/2+1,\dots ,N/2-1\right]}^T$. If an upsampling factor of $\alpha$ is to be employed for a region of $\left(1.5\times 1.5\right)$ pixels surrounding $\left(u_0,v_0\right)$, then the vectors $\mathbf{P},\mathbf{Q}$ are defined as:

$$\left(\mathbf{P},\mathbf{Q}\right)=\left(u_0,v_0\right)+\left[-\frac{1.5}{2},-\frac{1.5}{2}+\frac{1}{\alpha },\dots ,\frac{1.5}{2}-\frac{1}{\alpha }\right].$$E14

The three matrices in Eq. (13) are of size $\left(1.5\alpha \times N\right),\left(N\times N\right),\left(N\times 1.5\alpha \right)$, respectively, and their multiplication has marginal additional computational cost when the upsampling factor $\alpha$ is much less than hologram size $N$. Following this procedure, an oversampled local Fourier transform $\mid \tilde{H}\left(\mathbf{P},\mathbf{Q}\right)\mid$ is computed and its sub-pixel peak shift is determined. The ramp phase error $\mathrm{\Delta }\varphi$ in Eq. (12) can then be corrected.

In Figure 7, we show phase reconstruction for a red blood cell sample with and without the fractional fringe correction, which clearly shows the effectiveness of this simple algorithmic method in removing any residual phase ramp background. A 20-fold upsampling of a $\left(1.5\times 1.5\right)$ pixel region centered on integer pixel peak location was performed in this case, making it possible to detect small fractional fringe shift. Such upsampling of Fourier transform of the hologram can also be performed by zero-padding of the hologram prior to computing its Fourier transform. However, a 20-fold upsampling will make the total image size after zero-padding impractical when high-resolution Fourier transform is required only in the neighborhood of the carrier-frequency peak. The methodology may be used routinely with off-axis holograms with marginal additional computational burden.

Figure 7.

Illustration of fractional fringe elimination technique, (a) image plane hologram of RBCs, (b), (c) phase reconstruction with and without sub-pixel localization of carrier-frequency peak in Fourier domain. The phase ramp along the diagonal direction in (b) has been removed in (c).

4.2 Focusing of unstained cells

The unique advantage of a DHM system is that it can image unstained cell samples, since the contrast mechanism for imaging is the natural refractive index of the cells relative to their surrounding medium. The second interesting aspect is that since the complex-valued object field in the image plane is being recovered from the hologram, it is possible to numerically propagate this field to achieve computational refocus. One of the difficulties with this computational refocusing is that the typical Life Sciences users are not accustomed to such a methodology and often want to physically focus the cells by moving the sample stage in z-direction. Further, since most users without prior Optics or Physics background cannot interpret the phase images, they usually prefer to first record a bright-field image in focused position to know what they are imaging. The quantitative phase image obtained at the same sample location makes it easier for them to interpret the morphology of the cell being imaged. For this purpose, a DHM system can be built with a dual illumination so that it can be easily switched between the bright-field and the quantitative phase modes of operation. Focusing of unstained samples is however not trivial in bright-field mode as there is minimal amplitude contrast. Computational refocus is also not possible in the bright-field mode as the recorded image does not contain any phase information. The focusing problem can be addressed in a clever manner by employing a focusing criterion directly in the hologram domain [21]. Unstained samples can be considered to be nearly pure phase objects in the focus plane, and this implies that the interference fringes in the focus-plane hologram will predominantly show phase modulation. In other words in the focus plane, fringes mainly show bending without any amplitude modulation at the location of the cell. It is well known that when the sample is defocused, the phase information gets transferred to amplitude which causes amplitude modulation of fringes in the form of dark or bright halo near the cell boundaries. This effect is clearly illustrated in Figure 8 where we show three through-focus holograms of the same red blood cell as the microscope stage is translated in z-direction in small steps. The bottom row of Figure 8 also shows the corresponding phase reconstructions for these sample positions. Interestingly, the total phase variation is highest for the hologram in Figure 8(b) where the amplitude modulation of fringes is minimal. One may define a hologram domain focusing criterion as follows: The focus plane may be defined as the one with minimal amplitude modulation on the interference fringes. This criterion is very important for standardization. Note that the phase profiles in the defocus planes in Figure 8 are quite distorted in comparison with that in the focus plane. This in turn suggests that numerical phase map obtained by different users may be different for the same sample unless a uniform focusing criterion is used. Such a situation is not desirable considering that many imaging applications including DHM imaging are moving toward data-driven applications where numerical values of phase maps are very important in statistical learning, for example, for the purpose of cell classification applications. The minimum amplitude modulation criterion described here is easily implementable visually or by means of an autofocus algorithm in the hologram domain. The methodology also provides a seamless experience to potential Life Sciences users of DHM systems who may prefer to focus the sample physically.

Figure 8.

Illustration of hologram domain focusing criterion, (a), (b), (c): through-focus holograms of the same cell; (d), (e), (f): their corresponding phase reconstructions.

4.3 Fast unwrapping of two-dimensional phase maps

The phase map ${\varphi }_O\left(x,y\right)$ recovered from the single-shot image plane hologram may be interpreted approximately (at least for thin samples like mono-layer of cells) as:

$${\varphi }_O\left(x,y\right)=\frac{2\pi }{\lambda }\int \mathit{dz}n\left(x,y,z\right).$$E15

Here, $\lambda$ is the illumination wavelength, $z$ is the nominal propagation direction, and $n\left(x,y,z\right)$ stands for the relative refractive index of the sample with respect to the surrounding medium. Typically, the reconstruction algorithm provides the complex-valued object wave $O\left(x,y\right)$, and the phase ${\varphi }_{O,w}\left(x,y\right)$ is defined as the arc-tangent of the ratio of imaginary and real parts of the complex-valued function:

$${\varphi }_{O,w}\left(x,y\right)=arctan\left(\frac{\mathit{Im}\left[O\left(x,y\right)\right]}{\mathit{Re}\left[O\left(x,y\right)\right]}\right).$$E16

Since the arc-tangent function is defined over the interval $\left(-\pi ,\pi \right]$, the phase map ${\varphi }_{O,w}\left(x,y\right)$ in Eq. (16) has jump discontinuities that are not expected in ${\varphi }_O$ as per Eq. (15). The $-\pi$ to $\varphi$ phase jumps may occur as per the morphology of the sample in a complicated manner in the phase image. The phase jumps are however not physical and need to be removed. The unwrapped phase may therefore be expressed as:

$${\varphi }_O\left(x,y\right)={\varphi }_{O,w}\left(x,y\right)+2m\left(x,y\right)\pi ,$$E17

where $m\left(x,y\right)$ are integers that make the phase map continuous without the $2\pi$ jumps. Phase unwrapping in two dimensions is not a trivial problem, and multiple solutions have been proposed to address it. The popular solutions have come from radar literature, and their path following nature makes phase unwrapping a challenging image processing problem. A fast and robust solution to the two-dimensional phase unwrapping problem is possible [22] via an approach based on the transport of intensity equation (TIE). The TIE is a curious relation that relates the longitudinal intensity derivative of a propagating field (in Fresnel zone) to its transverse phase gradient. The TIE may be stated in the context of our problem as:

$$-k\frac{\partial I}{\partial z}={\nabla }_{\mathit{xy}}\cdot \left(I{\nabla }_{\mathit{xy}}{\varphi }_O\right).$$E18

Here, the gradient operator ${\nabla }_{\mathit{xy}}$ on the right-hand side denotes the transverse gradient and $I\left(x,y\right)$ denotes field intensity. This is probably the only relation in Optics which expresses the phase ${\varphi }_O$ without the use of the arc-tangent function but rather as a differential equation. When the intensity $I\left(x,y\right)$ is nearly constant as for a pure phase object in focus plane, the TIE simplifies to the Poisson equation. We proposed a novel scheme based on the TIE for addressing the unwrapping problem which is fast and robust. In particular, one may use the wrapped phase map ${\varphi }_{O,w}\left(x,y\right)$ obtained from $O\left(x,y\right)$ to generate an auxiliary field $u\left(x,y\right)=exp\left[i{\varphi }_{O,w}\left(x,y\right)\right]$. To estimate the longitudinal intensity derivative corresponding to this auxiliary field, it is propagated by small distances $\pm \mathrm{\Delta }z$ and the intensity derivative is obtained as:

$$\frac{\partial I}{\partial z}\approx \frac{{\left|u\left(x,y,\mathrm{\Delta }z\right)\right|}^2-{\left|u(x,y,-\mathrm{\Delta }z)\right|}^2}{2\mathrm{\Delta }z}+O\left({\left(\mathrm{\Delta }z\right)}^2\right).$$E19

Note that the distance $\mathrm{\Delta }z$ aforementioned is not determined by any experimental considerations like detector noise but is purely selected as a computational parameter which may be made small to improve the accuracy of the longitudinal intensity derivative. The TIE may be solved to get the unwrapped phase ${\varphi }_O\left(x,y\right)$. A popular approach for solving the TIE involves the use of inverse Laplacian operator:

$${\varphi }_O\left(x,y\right)=-\frac{k}{I}{\nabla }_{\mathit{xy}}^{-2}\frac{\partial I}{\partial z}.$$E20

The inverse Laplacian may be implemented readily in Fourier domain as follows:

$${\nabla }_{\mathit{xy}}^{-2}g\left(x,y\right)={\mathrm{ℱ}}^{-1}\left[\frac{D^{2}G\left(f_x,f_y\right)}{D^2+{\epsilon }^2}\right].$$E21

Here, $D^2=-4{\pi }^2\left(f_x^2+f_y^2\right)$ and ${\epsilon }^2$ is a small positive constant which avoids division by zero. This methodology offers a fast FFT-based solution to the phase unwrapping problem whose processing time is independent of the phase structure to be unwrapped. The phase unwrapping operation can thus be implemented with marginal additional computational efforts. In Figure 9, we illustrate unwrapping of a phase map obtained from image plane hologram of a pollen grain using the TIE-based solution. This technique is useful and works well except in cases where the transverse gradient of the phase map has a circulating (or curl) component.

Figure 9.

Illustration of TIE-based phase unwrapping; (a) wrapped phase map recovered from image plane hologram of a pollen grain and (b) unwrapped version of (a).

5. Does digital holography truly provide 3D information?

The readers may have noticed that in this article so far, we have avoided showing the phase maps associated with various objects in a 3D surface-rendered form. This is because the recovered phase ${\varphi }_O\left(x,y\right)$ using image plane digital holography is purely a 2D function and does not in itself contain any tomographic information. The 3D surface rendering is usually justified, since the phase map is related to the optical thickness of the sample under consideration. In this context, we provide a short discussion on the topic of whether phase imaging via DHM even provides a true 3D information. Holography if often naturally associated with “3D imaging.” In traditional film-based display holography, the holographic replay is performed by re-illumination of the recorded hologram by conjugate reference beam. When a human observer inspects the replayed wavefront visually, there is a perception of 3D objects which mainly arises because of the capability of eyes to focus onto high-contrast objects while ignoring any diffuse and blurred background. In digital holography, this ability of human eye to focus onto sharp objects has no important role to play, since the reconstruction is performed numerically. The problem of 3D reconstruction over $\left(N_x\times N_y\times N_z\right)$ voxels using object field $O\left(x,y\right)$ recovered from a digital hologram of size $\left(N_x\times N_y\right)$ also poses a problem of mismatch of degrees-of-freedom [23]. The 3D reconstruction problem has been attempted in the prior literature in two different ways. In one configuration, multiple (typically a few hundred) digital holograms of a transparent object are recorded from different views and the individual phase reconstructions are them combined using the Fourier diffraction theorem [24]. This method requires dedicated hardware where in order to change illumination angle or sample orientation or both. Such tomographic approaches also have relevance to the 3D reconstruction problems in cryo-electron microscopy of viruses or macromolecules [25]. In a second approach, a much smaller number of holograms are used and the reconstruction is performed via a sparse optimization method [26, 27]. The problem of dimensionality mismatch is attempted in such cases by exploiting the object sparsity. Tomographic imaging remains to be a challenging problem for digital holography, and much work needs to be done to model forward (object to hologram) and backward (hologram to object volume) propagation of light fields. The problem is particularly difficult when the sample is thick. Due to the complexity of the problem, some recent works have also explored machine learning-based approaches to tomographic phase (or refractive index) reconstruction [28]. This is likely to be a research problem of interest in future.

6. Popularization of DHM technology among life sciences researchers

Widespread usage of DHM technology among Life Sciences researchers and practitioners is somewhat lacking at present, and a lot needs to be done by DHM researchers to correct this status. The multiple advantages of the quantitative phase imaging modality suggest that it can be routinely used in practice just like bright-field or fluorescence microscopy. Here, I lay out a few steps that are needed going forward.

• Education about quantitative phase imaging: A quick look at biosciences and medical school curriculum suggests that quantitative phase as a modality is never introduced to them. This means that most Life Sciences researchers are not even aware of quantitative phase imaging. They might hear about it at advanced level as some exotic microscopic methodology but are then not inclined to study or use it for their own work. Even with very large number of publications by DHM researchers in Optics and Engineering journals, the situation may not change unless concrete efforts are made in this direction.

• Standardization of protocols: Life Sciences community is very particular about sample preparation and imaging protocols. The lack of standardized protocols for phase microscopy prevents multiple laboratories to communicate with each other effectively. The present article, in a limited sense, describes standardization issues associated with algorithmic methodologies. However, wide ranging efforts must be made in terms of benchmarking DHM hardware configurations, associated reconstruction algorithms, and their capabilities such as spatial resolution and phase accuracy. This will allow potential users to understand the modality much better before they decide to use it seriously. Such standardization is even more important if phase images are to be used with machine learning for designing diagnostic applications.

The popularization of DHM technology thus seems to be a task beyond what individual research groups can achieve in isolation. A DHM users’ consortium will be needed in order to move toward the aforementioned goals.

7. Conclusions

It this article, we summarized algorithmic methodologies that enable robust phase reconstruction performance for a single-shot digital holographic microscope system. First, we presented a sparse optimization-based phase reconstruction approach that provides single-shot full-resolution imaging performance. Further allied methodologies for fractional fringe correction, hologram domain focusing of unstained samples, and transport-of-intensity-equation-based fast 2D-phase unwrapping were discussed. Finally, we highlighted the importance and need of a DHM users’ consortium for benchmarking of DHM hardware and reconstruction algorithms, if this technology is to be popularized as a routine microscopy modality among Life Sciences researchers.

Conflict of interest

The authors declare no conflict of interest.