1. Introduction

Quantitative Phase Imaging (QPI) techniques enable non-invasive analysis of in vitro biological samples. In contrast to qualitative phase imaging (e.g. Zernike Phase Contrast [1] and Differential Interference Contrast [2]), QPI methods recover the phase delay caused by the sample, decoupled from absorption effects. Having full complex-field information makes it possible to digitally remove known aberrations [3] or to digitally refocus in post-processing [4]. System aberrations are typically measured once as a pre-calibration step, by conducting QPI with a guide star or ground-truth object. Then the aberrations can be computationally removed from subsequent measurements by a deconvolution process [3]. This strategy, however, requires precisely known calibration objects and will be sensitive to calibration drift and sample-induced aberrations.

Here, we propose a method of algorithmic self-calibration for Differential Phase Contrast (DPC) microscopy [5–10], where we avoid the need for pre-calibration by jointly recovering both the sample’s complex-field and the spatially-varying aberrations of the system, directly from raw images. The method is an extension of illumination-based DPC microscopy [6, 7], where images are captured with different source patterns, then a reconstruction algorithm recovers the complex-field. Experiments are implemented in a commercial brightfield microscope with a low-cost programmable LED array light source, such that patterns can be switched quickly with no moving parts [7, 11–13]. Among the wide variety of QPI methods, those which use partially coherent illumination, like DPC, are advantageous since they provide 2× better resolution, more light throughput and reduced speckle [7, 14–16], as compared to coherent methods.

The concept of algorithmic self-calibration – solving jointly for the reconstructed image and the calibration parameters – has proven particularly useful in coherent computational imaging. Examples include probe retrieval in Ptychography [17–20], source recovery for through-focus phase imaging [16, 21], pupil and source recovery for Fourier Ptychography Microscopy (FPM) [22–24], and calibration-invariant inverse scattering models [25]. The standard approach to self-calibration uses alternating projections (AP), which optimizes multiple variables serially, keeping other parameters fixed during each sub-iteration. The non-convexity of AP provides no guarantee of global convergence, but in practice it works with sufficiently diverse data. Self-calibration of aberrations has been demonstrated previously in an LED array microscope for the cases of through-focus phase [26] and FPM [13, 22, 27–29], but required a large number of images to be captured. For example, a typical FPM setup [30] uses approximately 5× as many images as our system to achieve the equivalent resolution.

Here we propose an AP framework that uses 4 captured images for simultaneous phase retrieval and digital correction of spatially-varying aberrations (Fig. 1). Three of the measurements are partially-coherent conventional DPC images (with rotated half-circle sources); these provide good phase contrast, but poor aberration contrast. The fourth image uses single-LED (spatially coherent) illumination; this provides aberration contrast, but alone cannot resolve the ambiguity between complex-field and pupil aberration [31]. Using both partially-coherent and coherent images together improves sensitivity to aberrations without sacrificing the benefits of partial coherence. We model the aberrations parametrically (with a Zernike basis [26, 32]) to dramatically reduce the number of unknowns. In addition, by segmenting the field-of-view (FoV), we are able to recover and digitally correct for spatially-varying aberrations across the FoV.

 

Fig. 1 Our LED array microscope captures 4 images with different illumination source patterns (three half-circles and one single LED). The intensity images are used to simultaneously reconstruct both amplitude and phase of the sample, and to estimate the pupil aberrations at each spatial location, which are then digitally corrected for. We show reconstructions for 4 regions with different spatially-varying aberrations.

Download Full Size | PDF

2. Joint estimation of complex-field and aberrations

Our LED array microscope and data capture scheme is shown in Fig. 1. Quantitative DPC (without aberration correction) requires a minimum of 3 intensity measurements to reconstruct phase [10], though only two quantities (amplitude and phase) are reconstructed at each pixel. Hence, there is significant redundancy in the data which may, in principle, be used to solve for aberrations. Unfortunately, intensity images formed by partially-coherent illumination do not exhibit significant aberration contrast. Hence, we modify the capture scheme to add one additional measurement with spatially-coherent illumination (a single on-axis LED). The example in Fig. 1 has a different aberration in each quadrant of the FoV, and only the single-LED image displays visible differences in contrast for each region. An off-axis LED would also provide the necessary coherent contrast, but the on-axis LED provides higher intensity and better signal-to-noise ratio (SNR). Having achieved both phase and aberration contrast with 1 on-axis and 3 half-circle sources, we use this dataset to jointly recover both spatially-varying aberrations and the complex-field of the sample (with resolution set by the incoherent diffraction limit).

The process of estimating the sample’s complex-field and the system’s aberrations simultaneously is a joint estimation large-scale nonlinear non-convex problem. To simplify, DPC typically makes a weak scattering approximation which linearizes the forward model. This approximation is generally valid for optically thin or index-matched samples, like biological cells. Intensity measurements can then be related to absorption (µ) and phase (ϕ) of the sample [7, 10] by simple convolutions. As derived in Appendix A, the DPC forward model in matrix form is:

Here, diag(x) denotes a diagonal matrix with diagonal values x, F and F⁻¹ represent the DFT and inverse DFT matrices, I_n is the vectorized normalized measured intensity (with the DC term subtracted), and H_µ and H_ϕ are vectorized transfer functions for absorption and phase. If the source satisfies the Köhler illumination configuration, the transfer functions can be numerically evaluated based on Eq. (13) and the cross-correlation property of Fourier transforms:

Our alternating projections algorithm initializes c with zero. At the start of one iteration, the Zernike coefficients c are fixed and a DPC deconvolution sub-procedure (Section 2.1) is performed in order to update the estimates of amplitude, exp(µ), and phase, ϕ. This new complex-field estimate is then held fixed while an aberration estimation procedure (Section 2.2) is performed. After updating the aberration estimate based on this procedure, a new iteration begins. Eventually, the objective function converges to a stationary point, giving the final estimates of amplitude, phase, and aberration coefficients. In general, this optimization strategy works as long as there exist enough diversity and redundancy in the measurements.

2.1. DPC phase retrieval sub-procedure

The phase retrieval sub-procedure amounts to solving the conventional DPC inverse problem [7, 10], except that it incorporates the current estimate of the Zernike coefficients c_k at iteration k:

The regularization term, R, should be chosen based on a priori information about the sample. For instance, Tikhonov regularization can mitigate noise, and the solution of Eq. (6) can then be found using a non-iterative deconvolution [10]. If the gradients of the object are relatively sparse, the Total Variation (TV) regularizer can be used to reduce noise without degrading edges. Since the regularization term for TV is not differentiable, iterative algorithms for solving Eq. (6) are needed. In this paper, we use the Alternating Direction Method of Multipliers (ADMM) to implement TV regularization [33], typically requiring about 20 iterations.

2.2. Aberration recovery sub-procedure

The aberration recovery sub-procedure uses a nonlinear optimization algorithm to update the aberration estimate based on the newly updated complex-field, which is held fixed. The sub-procedure is initialized with the Zernike coefficients estimate from the previous iteration, c_k. Mathematically, the sub-procedure problem is written as:

Equation (7) may be solved by a gradient descent (first-order optimization) approach, or more sophisticated second-order optimization routines (e.g. Newton’s method [21]). All of these require computation of the gradient of the objective function with respect to c. If we define the cost function as $f={\displaystyle {\sum }_{s=1}^{N_s}{\Vert {\epsilon }_s\Vert }_2^2}$, in which ε_s = FI_n,s − diag(H_µ,s)Fµ_k+1 − i · diag(H_ϕ,s)Fϕ_k+1 is the residual vector, the gradient becomes ${\nabla }_{\mathbf{c}}f={\displaystyle {\sum }_{s=1}^{N_s}{\left[\partial {\epsilon }_s/\partial \mathbf{c}\right]}^{\mathrm{H}}{\epsilon }_s}$, where H denotes Hermitian transpose. Using Eqs. (2)–(4), the gradient can be calculated analytically as:

In this gradient, T denotes transpose and the Zernike basis, Z^T = [Z₀, Z_1,…, Z_M ]^T, contains a finite number of modes where Z₀ − Z_M are the vectorized Zernike modes. For efficient computation, we adopt the L-BFGS algorithm [34] and use the gradient in Eq. (8) to solve this nonlinear optimization problem, which generally takes ∼ 10 iterations to converge.

3. Simulation results

To verify the performance of our joint estimation framework, we show simulation results in Fig. 2. The system parameters were chosen to match our experimental setup (0.4NA, wavelength 514nm, 177 source LEDs), with the LED array placed sufficiently far away from the sample such that the illumination from each LED is effectively spatially coherent (plane wave) [11, 13, 28].

 

Fig. 2 Performance of joint phase and aberrations estimation on a simulated dataset. (a) Simulated FPM and DPC measurements. Red dashed circles indicate the NA of the objective lens. (b) Joint estimation of optical field and pupil aberrations, comparing ground truth, FPM and DPC measurements. (c) Errors for complex-field and aberrations at each iteration.

Download Full Size | PDF

We compare our results with joint phase and aberration recovery FPM in Fig. 2. FPM captures a separate image for each of the 177 LEDs, whereas DPC requires only 4 images to reconstruct the same quantities. FPM intensity images are simulated by using different tilted plane wave illuminations corresponding to each single LED. All the intensity images contain the same pupil aberration, which is a weighted sum of the first 21 Zernike modes. Our simulated DPC measurements are the sum of the intensity images from half-circle blocks of LEDs on the top, bottom, or right regions of the LED array [7]. For all measurements, we added synthesized noise using a Poisson distribution with a mean of ∼3000 photons per pixel. Equation (5) was then solved with ℓ₂ regularization using the 4 DPC images, while we implemented the same algorithm in [35] to recover complex-field and aberrations using FPM with the full 177 image dataset. In this setup, FPM could use as few as 32 images (illuminations from the outer-most annular LEDs only) to achieve the same spatial frequency coverage as our DPC method. Therefore, we also include the results of FPM with 32 measurements for comparison.

Reconstructions from both FPM and our DPC algorithm match the ground truth (Fig. 2(b)); however, our method only requires 4 measurements, reducing acquisition time and memory requirements. Figure 2(c) plots the normalized root-mean-square error (║x − x_true ║₂/║x_true║₂) at each iteration for both complex-field and pupil aberrations. FPM incurs lower complex-field error than DPC, likely due to both the weak scattering approximation and the larger dataset. As for the pupil aberration error, FPM performance varies significantly with dataset size. While FPM with the full dataset outperforms the proposed DPC framework, our method provides a better result than FPM with 32 images. This is because FPM with fewer measurements has lower effective SNR, adding noise to the recovered pupil aberration (inset of Fig. 2(b)). Our method requires significant computation since it solves two optimization problems at each iteration, but the computation time is comparable to a sequential FPM reconstruction. Both methods were implemented in MATLAB on a desktop computer (Intel Core i7 CPU, Nvidia Tesla C2075 GPU). With a 650×584 pixel object, each iteration took 2.2s for FPM and 2.5s for our DPC algorithm.

4. Experimental results

Experimentally, we use an LED array microscope with the illumination module replaced with a custom-built LED array (λ = 0.514µm) [7, 13]. A phase target (Benchmark Technologies), which contains periodic patterns of continuous spatial frequencies, is imaged by a 20× 0.4NA objective lens (Nikon, CFI Plan Achro) in a Nikon TE300 microscope and images are recorded by a PCO.edge 5.5 sCMOS camera on the front port of the microscope (which adds 2× magnification). To test the weak phase gradient assumption, phase images of 6 resolution targets of different heights are recovered using DPC. After validating the reconstructed phase values against theoretical ones, we find that DPC provides accurate results when the phase of the sample is below 0.64 radians, then underestimates the phase values due to breakdown of the approximation (see Visualization 1). We collect 177 measurements by scanning individual LEDs within a maximum 0.4NA illumination angle. To provide a fair comparison, we use the same measurements to synthesize the 4 images for our method. As shown in Fig. 3(a), DPC measurements with half-circle source patterns have high resolution and qualitatively reveal the phase gradients of the sample, while measurements with single-LED illumination have lower resolution. Two single-LED image zoom-ins that contain the same structure at different orientations are shown in Fig. 3(a). One has high contrast, while the other does not, because of directional aberrations. By processing the images using the FPM algorithm and the proposed method with TV regularization, we recover the phase of the sample as shown in Fig. 3(b). The FPM and DPC reconstructions are similar, and can resolve features with period as small as λ/(2 × NA) = 0.643µm. In addition, both results provide reliable quantitative phase of the sample. The refractive index of the binary phase target is 1.52, with height of 100nm, resulting in ∼0.64 radians peak-to-valley. Looking at the 1D cut-lines (taken along dashed lines) in Fig. 3(b), after subtracting the mean of each, the reconstructions show good agreement with the ideal height.

 

Fig. 3 (a) Experimental FPM and DPC measurements for different LED source patterns. Zoomed regions at different orientations for coherent illumination are marked in cyan and pink boxes, respectively. (b) Quantitative phase of a star target using FPM and DPC, along with 1D cutlines for FPM (red) and DPC (blue) along the dashed lines. (c) Reconstructed wavefront error function and the weights of each Zernike mode up to the 4^th radial degree.

Download Full Size | PDF

Pupil aberrations recovered by both FPM and DPC algorithms are in Fig. 3(c). While we have no ground truth, aberrations estimated from FPM and DPC match well within the 4^th radial degree of Zernike modes. The dominant aberration in the objective lens is the 8^th Zernike mode (horizontal coma); this agrees with the evidence of directional aberration mentioned above. One reason causing a difference between the reconstructed pupils is the high-frequency fluctuation shown in the pupil aberration from FPM, which does not exist in the low-order Zernike modes. Although high-order aberrations might be estimated with more measurements to avoid overfitting, it’s usually enough to improve image quality by correcting the low-order aberrations.

To quantify the performance of our method, we introduced known defocus aberrations using an axial motion stage (Thorlabs, MZS500-E). We translated the phase test target across a range of known axial steps over a total range of 40µm. Each 2µm of translation, the 4 images for our method were acquired using a 10× 0.25NA objective lens (Nikon, CFI Plan Achro) at an acquisition rate of 20Hz. Quantitative phase reconstructions assuming zero aberration, aberration-corrected phase images and the recovered aberrations are shown in Fig. 4(a). With increased defocus, the uncorrected phase results degrade. Though our method also suffers from resolution degradation at large defocus distances (±20µm), it provides better resolution than the uncorrected phase retrieval. For example, in Fig. 4(b) the numbers in group 9 are resolved with pupil correction, while they are not readable in the uncorrected result.

 

Fig. 4 (a) Quantitative phase reconstructions of a USAF 1951 resolution target at various defocus distances with and without aberration correction, and the corresponding recovered aberrations (see Visualization 2). (b) Zoomed-in reconstructions at defocus of 10µm. (c) Known and experimentally-estimated defocus values from the 4^th Zernike mode over time.

Download Full Size | PDF

We can also recover the focus distances from the defocus term of the Zernike modes and compare the estimated defocus values to known values (Fig. 4(c)). Since the Zernike basis is normalized in a range from −1 to 1, the defocus value, d, at each time point, t, can be evaluated as d(t) = c₄(t)×λ/π/(1−(1−NA²)^0.5). As expected, the predicted pupil aberrations have quadratic form, which indicates that defocus dominates. Experimentally, our method overestimates defocus values by a factor of 1.16. This discrepancy might originate from mis-calibration of the experiment or parameters used in computation (e.g. wavelength, NA of the objective lens, precision of the motion stage). This linear relationship between the predicted positions and the true focus holds when the magnitude of defocus is less than 16 µm, ∼2× Depth of Field (DoF). However, the accuracy of aberration correction drops as the defocus value exceeds 2× DoF, when the maximum phase difference in the pupil becomes larger than 2π and the algorithm becomes more likely to converge to a local minimum. From these experiments, we conclude that our method accurately estimates time-varying aberrations within a total range of 4× DoF.

Our method is easily extended to account for spatially-varying aberrations, simply by solving the joint estimation problem separately over different patches of the FoV. In order to visualize spatially-varying aberrations, we prepared oil (Cargille, n_D = 1.58) immersed 10µm polystyrene beads (Sigma-Aldrich) as our sample, which is assumed to be nearly spatially invariant. The sample was imaged by a 4× 0.2NA objective lens (Nikon, CFI Plan Apo Lambda) that has a FoV of 1.7 × 2.1mm. Four images were captured at 12.5Hz using the same source patterns as in Fig. 3(a). Field-dependent aberrations are primarily caused by two types of system imperfections. First, the incident angle of individual LEDs changes from one FoV to another, since the relative position between the LEDs and the sample at each FoV varies. Second, there exist spatially-varying aberrations native to the optical system. To accommodate this spatial variance, we break the full FoV into small patches, and apply the joint estimation algorithm on each patch independently. In this case, the incident angle of individual LEDs and the aberrations are assumed invariant within the patch, which is usually valid when the region is much smaller than the total FOV. In Fig. 5(a), the full FoV is divided into 100 patches. Within each patch, the absorption, phase and pupil aberration were recovered independently before being stitched together to form the full FoV images. In practice, optimization for each patch converges in 20 – 30 iterations.

 

Fig. 5 (a) Reconstructed absorption, phase and spatially-varying aberrations (recovered pupil wavefronts for different regions of the field-of-view). (b) Comparison of results with and without pupil estimation for central and edge regions of the field-of-view.

Download Full Size | PDF

Looking at the aberrations recovered, the center of the FoV is essentially aberration-free, which is consistent with the fact that optical systems are usually optimized there. Therefore, the absorption and phase images with pupil recovery are similar to that without aberration correction. Aberrations are much stronger along the edges of the FoV, where the field curvature blurs the image. Consequently, obvious differences occur between the normal DPC and aberration-corrected DPC reconstructions at those patches. The recovered amplitude (absorption) at the bottom right corner is dramatically changed after applying aberration correction. Qualitatively, the uncorrected absorption at the edge of the field appears to contain some phase information, leading to a shadow-like appearance; these artifacts are removed with pupil recovery. In addition, while aberrated phase evaluated with Tikhonov deconvolution suffered from high-frequency noise, phase with pupil correction shows reduced noise due to the use of TV regularization.

5. Conclusion

In computational imaging systems, aberrations can significantly degrade complex-field reconstructions if not properly compensated for. In order to correct aberrations without pre-calibration, we propose joint estimation of the sample and the pupil aberration. Our DPC-based phase retrieval technique simultaneously recovers the system pupil aberration and quantitative phase and absorption of the sample from only 4 intensity images. By combining 3 half-circle source pattern measurements and 1 made under spatially coherent (single-LED) illumination, diffraction-limited complex-field and aberration Zernike coefficients were recovered through an alternating projections non-convex optimization. This method not only reduces the data acquisition time, but also increases the SNR due to higher light throughput compared to single-LED acquisition, which should aid in imaging dynamic biological samples. The recovered system aberrations are general to the system, and may be used for image correction in other subsequent imaging modalities on the same microscope, such as deconvolution of fluorescence images [29].