Automatic compensation of phase aberrations in digital holographic microscopy based on sparse optimization

In digital holographic microscopy, phase aberrations, which are usually caused by the imperfections of components and nontelecentric configuration of the optical system, severely affect the visuali...


I. INTRODUCTION
Digital holographic microscopy (DHM) is a powerful technique that can measure tiny structures such as MEMS circuits and biological specimens in real time due to its capability of recording the whole wavefront of a three-dimensional (3D) scene in a noninvasive and label-free way. By using a detector such as a CCD or a CMOS, the interference pattern between the object wave and the reference wave is recorded and stored in a digital manner. With the two-dimensional (2D) hologram, one can then numerically reconstruct the object's amplitude and quantitative phase distribution 1,2 using appropriate algorithms. By doing so, applications in tomography, 3 3D imaging, 4 and extended focused imaging 5 have been demonstrated in recent years.
Unfortunately, a serious distortion in the quantitative phasecontrast imaging of DHM is the undesired phase aberrations. On the one hand, in the off-axis configuration, the incidence between the object beam and the reference beam illuminating onto the detector produces a constant tilt aberration, which introduces a slant to the retrieved phase. On the other hand, the microscope objective (MO), if no additional component is added to form a telecentric configuration, inevitably causes a quadratic phase aberration to the extracted phase. Besides, due to the imperfect construction of the system and alignment of components, high-order aberrations, such as astigmatism and coma, may also exist in the final image. These aberrations appearing in the phase image not only distort the visualization of the object but also severely hinder the quantitative measurement of the object's true thickness/height. 6 In recent years, to eliminate the aberrations in the retrieved phase map, various methods have been proposed. By and large, they can be categorized into two groups: physical compensation in the optical recording stage and numerical compensation in the reconstruction stage. The former either requires additional/special optical components/configurations or demands multiple recordings. To ARTICLE scitation.org/journal/app name a few, by adding a tube lens, a telecentric configuration is formed, and therefore, the parabolic phase is physically compensated. [7][8][9] The tilt phase is further removed from the phase map by a spatial filtering. 9 Similarly, by inserting an electrically tunable/adjustable lens in the illumination path or in the reference arm, the parabolic phase factor or total phase distortions can be canceled out. [10][11][12] Besides, a common-path interferometer with one singlecube beam splitter is also demonstrated to physically compensate phase aberrations. 13 As for the latter scheme of multiple recordings, although the double exposure method is capable of eliminating all of the aberrations, a second object-free reference hologram has to be captured. 14,15 Under some situations, even two more laterally shifted holograms have to be collected. 16 As such, these methods are only appropriate for static scenarios. All in all, the physical compensation strategy complicates the holographic recording and makes it more bulky.
Another group of compensation methods, as mentioned above, is based on numerical postprocessing. Roughly speaking, fittingbased numerical methods involve in the use of different polynomials such as standard polynomials, 17,18 digital phase mask, 19,20 and Zernike polynomials. 17,21 By solving a least-squares problem 19,21 or a phase variation minimization problem, 20 the aberration surface is thus fitted and then subtracted from the phase map, resulting in a pure object phase finally. Comparatively, nonfitting methods employ various mathematical models and special operations on the hologram. For example, Zuo et al. used principal component analysis (PCA) to decompose the phase map into a set of values of uncorrelated variables, thus extracting aberration terms from the first principal component obtained. 22 A method that analyzes the frequency spectrum of a hologram for aberration removal is also demonstrated. 23,24 Besides, a hologram rotation-based method is proposed, in particular, for getting rid of the off-axis tilt as well as the parabolic phase aberration. 25,26 Despite the success and flexibility of current numerical approaches, manual operations of selecting a flat region as reference surfaces 17,19 are required, an assumption is made that the object is so thin that it is just a small perturbation in the aberration contribution to the overall reconstructed phase distribution, 18,21 or only low-order aberrations such as the phase curvature introduced by the MO can be removed. 18 It is worth noting that the polynomial-based method is a semiautomatic method since it requires background information to find the phase residual, which is detected by manually cropping a flat background area as Refs. 17 and 19. Besides, although PCA can automatically predict the phase residual by creating a self-conjugated phase to compensate aberrations and can be extended to deal with high-order aberrations, it assumes that only noncross terms exist, limiting the scope of this method. 22,27 Furthermore, manual filtering and centering 23 or hypothesis of a thin phase object 24 is inevitable in the spectrum-based analysis, and the self-hologram rotation is solely suitable for tilt aberration compensation. 25 Recently, a learning-based method is proposed to automatically detect the specimen-free background for the subsequent fitting using Zernike polynomials. 28 Aberration removal is further encapsulated in holographic reconstruction and thus can be automatically compensated during the feed-forward propagation along the network. 29 Despite its success and applicability to high-order aberrations, the learning-based method requires a lot of data and extensive computation for network training.
In this paper, we propose a purely numerical method that can remove phase aberrations of any orders in a totally automatic manner for DHM. By formulating the surface fitting as an inverse problem and by solving an ℓ 1 -norm-based optimization, the aberration surface can be accurately estimated and a high-quality aberrationfree quantitative phase-contrast image is then recovered. Synthetic and experimental results of multiple samples are demonstrated to verify the efficacy of the proposed method.

II. METHOD
On the reconstruction plane (x, y), based on the recovered complex wavefront Γ(x, y), the unwrapped and phase map is extracted by where UNWRAP[⋅] indicates the phase unwrapping operation and Re and Im denote the real and imaginary parts, respectively. Under situations that no physical compensation is implemented, taking aberrations into account, the recovered phase is the summation of the object phase and the aberration, which reads where ϕ obj (x, y) is the object phase distribution and ϕ abe (x, y) is the total aberration term overlapping with ϕ obj (x, y). The purpose of numerical aberration compensation is to find an approximation, ϕ abe (x, y), of the real aberration ϕ abe (x, y) as accurately as possible.
As such, the aberration-free object phase can be then fully recovered byφ In order to approximate the aberration, a theoretical fitting model is necessary. According to the principle of homogeneous polynomial fitting, the surface model we are using here is in secondorder and can be written aŝ where βi, i = 1, . . ., 6, denote the fitting coefficients to be estimated. Note that for mathematical simplification, the model we show here is second-order. This illustrates that the commonly occurred aberrations in DHM, such as tilt, astigmatism, and parabolic surface, can be successfully compensated using this fitting model. However, the extension of the fitting model to higher order is straightforward without any further efforts if high-order aberrations have to be removed. With the help of the fitting model, we can then proceed to estimate its coefficients, which determines the aberration surface. To do so, by assuming that the object phase is a small perturbation compared to the aberration (this is a reasonable hypothesis for the situations in DHM 18,21 ), we can rewrite ϕ(x, y) =φ abe (x, y) in the matrix notion as ARTICLE scitation.org/journal/app . . , bn) T denotes the phase ϕ(x, y), and n is the number of pixels of the phase image. Suppose that A denotes the matrix of the polynomials, Eq. (5) is written as This is a typical inverse problem, which is usually ill-posed. Conventionally, this problem is solved by the least squares method, in which the residual between the actual value and the value predicted by the model is minimized with the ℓ 2 -norm loss function. 18 Instead, to seek the estimation of ⃗ β, by incorporating regularization technique, which is a process of introducing additional information in order to prevent overfitting, the objective function we propose to minimize is where ∥⋅∥ 1 denotes the ℓ 1 -norm and λ is a regularization parameter determined by trial and error. Then, by solving Eq. (7), the fitting coefficients can be acquired bŷ Therefore, after having a proper approximation of ⃗ β, according to Eq. (4), an estimated phase aberrationφ abe (x, y) can be obtained. By using Eq. (3), an aberration-free object phase distribution can then be recovered for visualization and quantitative measurement. The process is totally automatic since no flat region has to be selected as the reference in our method. The only necessary information is the recovered phase image with aberration, which is always available.
To explain why changing from the ℓ 2 -norm to ℓ 1 -norm can improve the fitting accuracy, we turn to the well-known compressed sensing theory. As is well known, the least squares method for regression analysis is to minimize an ℓ 2 -norm loss function. As such, this optimization problem is equal to minimizing the amount of energy in the system and is mathematically simple to solve. However, such a formulation leads to poor fitting results for many practical applications for which the unknown coefficients have nonzero energy. Instead, to enforce the sparsity constraint when solving linear equations, one can minimize the number of nonzero components of the solution, which is the so-called ℓ 0 -norm in a broad sense. 30 This means that finding the sparsest solution x, which is feasible to sensing constraints given by the encoding matrix A and the observation vector b, can be summarized as an ℓ 0 minimization problem. Furthermore, in mathematics, it has been proved that for many problems such as the compressive imaging problem here, the ℓ 1 -norm can be used in lieu of the ℓ 0 -norm. Consequently, one can solve the ℓ 1 problem instead of the ℓ 0 problem since the former is less difficult and easy to solve mathematically. 31,32 Therefore, formulating an ℓ 1 -norm loss function has a better performance than the ℓ 2 -norm for our case.

A. Simulation
We first test the proposed method on synthesized data, which is implemented on a PC with Intel Xeon W-2145 at 3.70 GHz and 64 GB RAM. Figures 1(a) and 1(b) present a computer-generated phase map with a size of 256 × 256 on a flat background in 2D and 3D visualizations, respectively. In Fig. 1(c), a synthetic phase aberration is generated using the first 6 terms of the Zernike polynomials, 21 leading to tilt aberration in the x and y axes, astigmatism at 0 ○ and 45 ○ , and defocus. The coefficients used in the simulation are given in Table I. By adding up the object phase and aberration, an aberration-distorted phase is then created and shown in Fig. 1(d).
To remove the phase aberration, the proposed ℓ 1 -norm method and conventional least squares method are employed separately and solved using the SDPT3 4.0 algorithm, 33 and the compensated results are shown in Fig. 2. As can be seen, with the Astigmatism at y 0.5 5 x 2 − y 2 Astigmatism at x 0.0 ℓ 1 -norm-based optimization, the aberration is totally removed and the recovered object phase has exactly the same phase range to the original one, while the conventional method gives an inaccurate compensated background. This can be further observed by the quantitative comparison of phase profiles across the individual central lines in Fig. 2(e). As such, the superior performance of the proposed method can be verified.
Based on this model, we use the ℓ 1 -norm-based method and the least squares method for compensation. The results are given in Figs. 3(b) and 3(c), respectively. As can be clearly seen, the conventional method can basically work in this case. However, the bottom plane is not flat and contains fluctuation, which would introduce errors in the calculation of the object phase. While the proposed method can give a pretty flat bottom plane as well as a precise compensated result. The synthesized result here demonstrates an intuition of the convenience of extending the proposed method to higher-order aberration compensation.

B. Experiment
We further demonstrate the performance of this method using experimental data. The DHM system, as shown in Fig. 4(b), used to acquire data is an off-axis Mach-Zehnder interferometer with a 4X MO (GCO-2101, L = 45 mm, Daheng Optics) working in the reflection mode. The wavelength of the He-Ne laser source is 632.8 nm, and the pixel pitch of the detector (MER-130-30UM-L, 1280 × 1024 pixels, Daheng Imaging) is 5.2 μm. Two samples, a groove etched on an optical wafer shown in Fig. 4(c) and a microcircuit on a sapphire wafer shown in Fig. 4(d), are used as objects. 34 The collimated laser beam is split by the PBS and enters the interferometer along two paths. The reference beam, which is reflected by M1, directly goes into the detector, while another beam, which is reflected by M2, focuses onto the OBJ through the MO and illuminates the sample. Then, the object beam, which is the beam reflected by the sample, goes back to the detector. Two beams interfere and lead to an

ARTICLE
scitation.org/journal/app interference pattern, which is the hologram and is collected by using the detector. The PBS and two HWPs are designed to collaboratively change the polarization state, leading to the interference of two beams as well as the adjustment of the contrast of the interference pattern. Since it is challenging to intentionally create specific aberrations in the experiment, in this section, removal of two commonly occurred aberrations in DHM, tilt aberration and quadratic aberration, is demonstrated.

Removal of linear aberration
To ensure only tilt aberration is remained, the quadratic aberration introduced by the MO is optically corrected with the telecentric configuration. 8,35 In Fig. 5(a), the hologram of the groove is presented, and for a clear comparison, a true groove phase shown in Fig. 5(b) is acquired with a totally new DHM setup. All hardwares and configurations are different from our system used here. For a better visualization, the acquired phase map is cropped and magnified properly. Figure 5(c) shows the recovered phase image with aberration. As can be clearly seen, the groove is located on a tilt plane, even with some astigmatism at corners, while in Fig. 5(d), the tilt surface is successfully removed using the proposed method, and the absolute height of the groove can be thus easily measured according to the difference of the groove and the flat surface. Profiles of the individual diagonals in the aberrated and corrected phase images are also demonstrated in Fig. 5(e) for comparison. The true phase of the groove can therefore be measured as around 4 rad. Another example, as shown in Fig. 6, presents the tilt aberration removal for the sample circuit. Figure 6(a) shows the hologram of the circuit, and the true circuit phase is given in Fig. 6(b). Similarly, the recovered phase in Fig. 6(c) has a tilt distortion, while in Fig. 6(d), the slope is removed such that the true phase of the circuit can be measured. Figure 6(e) compares the phase profiles along central lines in the two phase images.

Removal of quadratic aberration
As for the compensation of quadratic aberration, we first capture a hologram of the groove without telecentric configuration, leading to the introduction of quadratic curvature. Figure 7 shows the recovered phase-contrast image in which the object phase is overlapped with a quadratic surface. The curvature, as can be seen clearly, introduced by the MO significantly hinders the measurement of the true thickness of the object. By using the proposed method, in Fig. 7(b), however, the quadratic surface is eliminated and the groove is thus located on a flat plane. From Fig. 7(c), we can see that removal of the notorious quadratic surface in phase imaging benefits the quantitative measurement of the object. Then, with the telecentric configuration, holograms of the groove and the circuit are collected. Note that although in theory, parabolic aberration should be removed in such configurations. However, due to the misalignment of components, the quadratic aberration is not totally removed. As such, the proposed method is utilized to compensate the second-order aberration. Figure 8 shows the amplitude reconstruction of the same groove. Figures 8(b) and 8(c) show the aberrated and aberration-free object phase maps, respectively. The existence and removal of quadratic surface are clearly demonstrated in Fig. 8(d). A similar example of imaging the microcircuit, as shown in Fig. 9, also supports our method.

IV. SUMMARY
To conclude, a numerical method is proposed in this paper for automatically compensating aberrations in phase-contrast imaging using DHM. By formulating the aberration fitting as an inverse problem and by adopting the ℓ 1 -norm as the objective function and regularization, an accurate aberration surface can be approximated. As such, by subtracting the estimated aberration surface from the original phase image, the true object phase image can then be recovered for visualization and quantitative measurement. Simulated and experimental results of multiple samples are demonstrated to verify the efficacy of the proposed method in quantitative phase-contrast imaging.
Despite the outperformance and automation of the proposed method here, it has a limitation that more efforts should be engaged in the future. Since solving an ℓ 1 -norm-based optimization problem (∼10 s) requires more time than that of the ℓ 2 -norm (∼1 s), this method has a lower computational efficiency (ten times slower) than the least-squares method. As such, for real-time applications such as particle tracking or live cell imaging, this method can only be employed in an offline fashion. Future directions of the reported approach include studies of the cause that makes the optimization slow. We will try faster optimization algorithms or software packages like the homotopy method or the MOSEK optimization suite to accelerate the computation of the proposed method.