Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Fourier space aberration correction for high resolution refractive index imaging using incoherent light

Open Access Open Access

Abstract

An aberration correction method is introduced for 3D phase deconvolution microscopy. Our technique capitalizes on multiple illumination patterns to iteratively extract Fourier space aberrations, utilizing the overlapping information inherent in these patterns. By refining the point spread function based on the retrieved aberration data, we significantly improve the precision of refractive index deconvolution. We validate the effectiveness of our method on both synthetic and biological three-dimensional samples, achieving notable enhancements in resolution and measurement accuracy. The method's reliability in aberration retrieval is further confirmed through controlled experiments with intentionally induced spherical aberrations, underscoring its potential for wide-ranging applications in microscopy and biomedicine.

© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Achieving theoretical imaging resolution is a fundamental goal in optical systems, which depend on the precise alignment of finely polished lenses for high-resolution imaging [1]. Microscopes are susceptible to image quality degradation due to even minor misalignments in their components. Apart from mechanical alignment issues, microscopy frequently contends with sample-induced aberrations. These aberrations, originating either from the sample itself or its substrate, can substantially impair image quality. While Fourier space aberrations represent only one category of aberration and do not encompass other aberrations like image deformation or field curvature, they are of particular interest in microscopy due to their profound impact on resolution. Furthermore, many complex aberrations can often be represented as spatially variable Fourier space aberrations [2]. Consequently, effectively characterizing and rectifying these Fourier space aberrations is imperative for enhancing imaging quality.

Recent advancements in quantitative phase imaging (QPI) techniques have been transformative in various applications [37], particularly due to their label-free and quantitative analysis capabilities [810]. Among these, phase deconvolution microscopy stands out for its simplified optical setup using an incoherent light source [1115]. This setup is notably resilient to noise and vibration, especially when compared to other QPI methods that rely on interferometry with coherent light sources [1416]. This simplicity and robustness have facilitated its application in label-free 3D quantitative imaging of various thick biological samples, including live Caenorhabditis elegans [17] and renal cancer tissues [18].

The significance of aberration correction in microscopy has led to the development of numerous correction methods for QPI techniques. These encompass algorithms for imaging thin samples, including 2D differential phase contrast microscopy [19] and ptychography [20,21]. Other approaches to aberration correction have been implemented, such as background subtraction [22], lateral shifting [23], the 3D angle scanning method [24], and deep-learning-based strategies [25,26]. Additionally, some techniques involve pre-imaging retrieval of system aberrations under specific experimental conditions [27,28]. While aberration correction and adaptive optics methods are well-established in fluorescence microscopy [2,29], their direct application to QPI is not straightforward.

In this paper, we propose a versatile method for aberration correction in 3D phase deconvolution microscopy. Our approach, inspired by techniques used in 2D differential phase contrast microscopy [19] or ptychography [20,21], retrieves aberration from redundant information obtained under varied illumination patterns. This method addresses critical limitations in current microscopy techniques, particularly in analyzing thick tissue specimens where aberrations can significantly distort imaging results. By leveraging redundant information from varied illumination patterns, our approach ensures more accurate aberration correction, crucial for detailed and reliable analysis in biological and medical research. It is uniquely suited for complex samples where traditional methods fall short, enabling precise and deep insights into specimen structure and dynamics.

2. Principle and methods

2.1 Aberration in 3D phase deconvolution microscopy

Our approach to aberration correction within 3D phase deconvolution microscopy initiates with the sample being illuminated via an LED source. This illumination is precisely modulated by a spatial light modulator (SLM) (Fig. 1). The acquisition of a 3D intensity stack is accomplished by the precise axial scanning of the sample stage. In the domain of phase deconvolution microscopy, the computation of the refractive index (RI) is executed through deconvolution with point spread functions (PSFs) that are inferred from these specific illumination patterns.

 figure: Fig. 1.

Fig. 1. (a) Schematic representation of the optical arrangement, detailing the path of light through the system components. (b) Aberration extraction process, where a gradient descent algorithm iteratively refines the aberration model to align with the correlation between simulated and empirical intensity profiles. (c) Enhancement sequence demonstrating aberration mitigation and the consequent recovery of the high-fidelity tomographic image.

Download Full Size | PDF

To augment this methodology, we engineered an algorithm designed to extract Fourier space aberrations from the observed intensity under various illumination patterns (Fig. 1(b)). Within this framework, Fourier space aberrations are conceptualized as the product of the transmitted light field E and an aberration function A(kx, ky) at the microscope's pupil plane. The Fourier transform, symbolized by ∼, facilitates the mathematical representation of aberrations. The transmitted field can then be expressed as $\widetilde E \to \widetilde E{e^{iA({k_x},{k_y})}}$.

The normalized transmitted intensity S is approximated as the convolution of the scattering potential’s real and imaginary parts ${F_{Re }},{F_{{\mathop{\rm Im}\nolimits} }}$ with two point spread functions ${H_A},{H_P}$, where ${{\cal F}^{ - 1}}$ denotes the inverse Fourier transform [12,30,31].

$$\begin{array}{c} S = {F_{Re }} \ast {H_P} + {F_{{\mathop{\rm Im}\nolimits} }} \ast {H_A}\\ {H_A} ={-} 16{\pi ^2}Re [{{{\cal F}^{ - 1}}[{({\widetilde {{G_{\textrm{obj}}}}} )\ast ({{k_z}{{|\rho |}^2}\widetilde {{G^\ast }_{\textrm{obj}}}} )} ]} ]\\ {H_P} ={-} 16{\pi ^2}{\mathop{\rm Im}\nolimits} [{{{\cal F}^{ - 1}}[{({\widetilde {{G_{\textrm{obj}}}}} )\ast ({{k_z}{{|\rho |}^2}\widetilde {{G^\ast }_{\textrm{obj}}}} )} ]} ]\end{array}.$$

The aberration in this model is included in the Green function, ${G_{obj}}$

$$\widetilde {{G_{\textrm{obj}}}} = \widetilde {{P_{obj}}}\widetilde G,\textrm{ }\widetilde {{P_{obj}}} = \frac{{\sqrt {{k_z}} }}{k}\delta ({k_x^2 + k_y^2 < \textrm{N}{\textrm{A}^2}/{\lambda^2}} )\cdot {e^{iA({{k_x},{k_y}} )}}.$$

In phase deconvolution microscopy, the intensity of transmitted light is captured under various patterned illuminations, denoted as |ρi|2, where i represents the sequential number of the illumination. These illuminations provide overlapping data regarding the specimen's characteristics. Typically, a single illumination pattern suffices for reconstructing transparent objects, while a minimum of two patterns are required for simultaneous retrieval of the refractive index and absorption properties. Employing additional illumination patterns can be beneficial, enhancing the overall image quality [15]. A detailed comparative analysis of the effects of using different numbers of illumination patterns can be found in Supplement 1.

Our optimization algorithm identifies the optimal Fourier space aberration that aligns the redundant information obtained from the sample under diverse illumination conditions. The process involves simulating the anticipated intensity based on the reconstructed RI and calculating the mean squared error relative to the experimentally captured intensity. Subsequently, we employ a gradient descent technique to pinpoint the aberration that minimizes this discrepancy, as described in Fig. 1(b).

2.2 Retrieval of aberration information

Our method is tailored for transparent objects, as they yield more redundant information for a set number of illumination patterns, facilitating the identification of the optimal aberration. This scenario is mathematically represented as FIm = 0. Although most biological samples meet this criterion, it is important to acknowledge that the phase transfer function exhibits smaller values at lower frequencies compared to the absorption transfer function. Consequently, our approach capitalizes on high-frequency information to determine aberrations, thereby ensuring applicability to a broad spectrum of biological samples.

Mathematically, the optimization problem can be expressed as follows:

$$\begin{array}{c} {A_{optimal}} = \mathop {\arg \min }\limits_A \{{C(A )} \},\textrm{ }C(A )= \sum\limits_i {||{{\varepsilon_i}} ||_2^2} ,\textrm{ }{\varepsilon _i}(A )= {S_i} - {F_{Re }} \ast {H_{P,i}}(A )\\ {H_{P,i}}(A )={-} 16{\pi ^2}{\mathop{\rm Im}\nolimits} [{{{\cal F}^{ - 1}}[{({\widetilde {{G_{\textrm{obj}}}}(A)} )\ast ({{k_z}{{|{{\rho_i}} |}^2}\widetilde {{G_{\textrm{obj}}}^\ast }(A )} )} ]} ]\end{array}$$
$$\begin{array}{c} \widetilde {{G_{\textrm{obj}}}}(A) = \widetilde {{P_{obj}}}(A)\widetilde G,\textrm{ }\widetilde {{P_{obj}}}(A) = \frac{{\sqrt {{k_z}} }}{k}\delta (k_x^2 + k_y^2 < \textrm{N}{\textrm{A}^2}/{\lambda ^2}) \cdot {e^{iA({k_x},{k_y})}}\\ {F_{Re }} = {{\cal F}^{ - 1}}[\sum\limits_i {\widetilde {{S_i}}} {\widetilde {{H_{P,i}}}^\ast }(A)/\sum\limits_i {{{|{\widetilde {{H_{P,i}}}(A)} |}^2}} ] \end{array}$$

Here, i represents the illumination number, Aoptimal is the optimal aberration that best describes the observed transmitted intensity, and C(A) is the cost function to be optimized. We solve this optimization problem using a gradient descent method, specifically the FISTA method [32,33]. To facilitate gradient computation, we discretize our problem by representing the field and scattering potential on a discrete grid, converting these 3D quantities into vectorized arrays. Consequently, the Fourier transform is treated as a discrete Fourier transform, expressed as matrix multiplication. For a given discretized vector B, such as the field or the scattering potential its Fourier transform is written as $\widetilde B = UB$ with $U$ the matrix composed of the Fourier basis vectors. Equivalently the inverse Fourier transform is given by $B = {U^{ - 1}}\widetilde B$. The optimization problem can now be written as:

$$\scalebox{0.9}{$\begin{array}{l} {A_{optimal}} = \mathop {\arg \min }\limits_A \{ C(A)\} ,\textrm{ }C(A) = \sum\limits_i {||{\widetilde {{\varepsilon_i}}} ||_2^2} ,\textrm{ }\widetilde {{\varepsilon _i}}(A) = \widetilde {{S_i}} - \textrm{diag}(\widetilde {{H_{P,i}}})\widetilde {{F_{Re }}}\\ \widetilde {{H_{P,i}}}(A) = 16i{\pi ^2}U({\textrm{diag}({{U^{ - 1}}\widetilde {{G_{\textrm{obj}}}}} )U\textrm{diag}({k_z})\textrm{diag}({{|{{\rho_i}} |}^2}){{\widetilde {{G_{\textrm{obj}}}}}^\ast } - \textrm{diag}({U{{\widetilde {{G_{\textrm{obj}}}}}^\ast }} ){U^{ - 1}}\textrm{diag}({k_z})\textrm{diag}({{|{{\rho_i}} |}^2})\widetilde {{G_{\textrm{obj}}}}} )\\ \widetilde {{G_{ {\pm} \textrm{obj}}}}(A) = \textrm{diag}(\widetilde {{P_{obj}}})\widetilde G,\textrm{ }\widetilde {{P_{obj}}}(A,{k_x},{k_y}) = \frac{{\sqrt {{k_z}} }}{k}\delta (k_x^2 + k_y^2 < \textrm{N}{\textrm{A}^2}/{\lambda ^2}) \cdot {e^{iA({k_x},{k_y})}}\\ \widetilde {{F_{Re }}}(A) = {{\sum\limits_i {\textrm{diag}({{{\widetilde {{H_{P,i}}}}^\ast }} )\widetilde {{S_i}}} } / {\sum\limits_i {\textrm{diag}({{{|{\widetilde {{H_{P,i}}}} |}^2}} )} }} \end{array}$}$$

Note that here we use the diag(B) function to form a diagonal matrix whose diagonal elements are the elements of B. In contrast to previous equations where the multiplications were element wise multiplications, here all multiplication are matrix multiplications. In order to perform the optimization, we need to compute the gradient of the cost function with respect to the aberration vector: ${\nabla _A}(C(A))$. Our discretized formulation now enables us to use matrix calculus identity to compute the gradient. Using the chain rule formula for the gradient of the norm [34] we can express the gradient of the cost function in following form. See Supplement 1 for detailed derivation.

$$\begin{array}{c} {\nabla _A}C(A) = 2\textrm{Re} \left[ { - \sum\limits_j {{{\left( {\frac{{\partial \widetilde {{H_{P,j}}}}}{{\partial A}}} \right)}^H}\textrm{diag}({{{\widetilde {{F_{\textrm{Re} }}}}^\ast }} ){\varepsilon_j}} - {{\sum\limits_i {\left( {\frac{{\partial {{\widetilde {{H_{P,i}}}}^\ast }}}{{\partial A}}} \right)} }^H}\textrm{ diag}({{{\widetilde {{S_i}}}^\ast }} )\frac{{\sum\limits_j {\textrm{diag}({H_{P,j}}^\ast ){\varepsilon_j}} }}{{\sum\limits_k {\textrm{diag}({{{|{\widetilde {{H_{P,k}}}} |}^2}} )} }}} \right.\\ \left. { + \sum\limits_i {\left( {{{\left( {\frac{{\partial \widetilde {{H_{P,i}}}}}{{\partial A}}} \right)}^H}\textrm{diag}({\widetilde {{H_{P,i}}}} )+ {{\left( {\frac{{\partial {{\widetilde {{H_{P,i}}}}^\ast }}}{{\partial A}}} \right)}^H}\textrm{diag}({{{\widetilde {{H_{P,i}}}}^\ast }} )} \right)} \textrm{ diag}({{\widetilde {{F_{\textrm{Re} }}}}^\ast })\frac{{\sum\limits_j {\textrm{diag}({{\widetilde {{H_{P,j}}}}^\ast }){\varepsilon_j}} }}{{\sum\limits_k {\textrm{diag}({{{|{\widetilde {{H_{P,k}}}} |}^2}} )} }}} \right] \end{array}$$
$$\scalebox{0.7}{$\begin{aligned}&\left( {\displaystyle{{\partial \widetilde{{H_{P,j}}}} \over {\partial A}}} \right)^H = 32\pi ^2P_z{\mathop{\rm Re}\nolimits} \left[ {{\rm diag}\left( {\widetilde{{G_{{\rm obj}}}}} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right){\rm diag}\left( {k_z} \right)U^{-1}{\rm diag}\left( {U{\widetilde{{G_{{\rm obj}}}}}^*} \right)-{\rm diag}\left( {{\widetilde{{G_{{\rm obj}}}}}^*} \right)U{\rm diag}\left( {U^{-1}{\rm diag}\left( {k_z} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right)\widetilde{{G_{{\rm obj}}}}} \right)} \right]U^{-1}\\& \left( {\displaystyle{{\partial {\widetilde{{H_{P,j}}}}^*} \over {\partial A}}} \right)^H = 32\pi ^2P_z{\mathop{\rm Re}\nolimits} \left[ {-{\rm diag}\left( {{\widetilde{{G_{{\rm obj}}}}}^*} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right){\rm diag}\left( {k_z} \right)U{\rm diag}\left( {U^{-1}\widetilde{{G_{{\rm obj}}}}} \right) + {\rm diag}\left( {\widetilde{{G_{{\rm obj}}}}} \right)U^{-1}{\rm diag}\left( {U{\rm diag}\left( {k_z} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right){\widetilde{{G_{{\rm obj}}}}}^*} \right)} \right]U\end{aligned}$}$$

To recover the phase of the aberration, we implement FISTA, without the use of a regularization term. When starting the first iteration the aberration was initialized to zero. The omission of regularization is due to its tendency to confine the optimization within local minima, as a result of the phase's modulo nature. Post-optimization, a low-pass filter is applied to the aberration to mitigate high-frequency noise, effectively reducing the resolution by a factor of four. Given the circular configuration of the aberration image, we perform iterative low-pass filtering and enforce circular symmetry over 2000 iterations. This process is designed to precisely align the low-frequency components with the observed aberration pattern on the circle.

After aberration retrieval, we proceed with the reconstruction of the RI, considering both its real and imaginary components, employing the algorithm detailed in Ref. [15]. To mitigate artifacts stemming from the ‘missing cone problem'—an inherent issue in tomographic reconstruction—we incorporate total variation regularization into our process

To boost computational efficiency, we carried out the optimization process on a graphics processing unit (GPU) utilizing MATLAB's parallel computing toolbox. For example, 100 FISTA iterations on a 701 × 701 × 79 pixel volume were completed in just 202 seconds using an NVIDIA RTX 3090 GPU. This approach synergizes computational algorithms with advanced imaging methods, optimizing aberration correction in 3D phase deconvolution microscopy, especially for complex or dense samples.

2.3 Experimental setup

To demonstrate our method, we utilized a custom-built phase deconvolution tomography setup. The illumination pattern was controlled at the pupil plane of the microscope using a ferroelectric spatial light modulator (Forth dimension display M150). Subsequently, the illumination was then projected on the sample using a high numerical aperture (NA) condenser lens (Olympus UPLSAPO60XW, NA 1.2). Scattered light was collected by an objective lens (Olympus UPLSAPO60XW, NA 1.2) and projected to form an image on the camera using an achromatic tube lens. The imaging system was equipped with a monochrome camera (Oryx ORX-10G-71S7M-C) with pixel size 4.50 µm, and illumination was provided by an LED (Thorlabs Inc., M595L4) emitting at a nominal wavelength of 595 nm. A schematic of the setup is presented in Fig. 1(a).

In the experiments related to Fig. 4, the objective and condenser were changed to a slightly lower NA lens (Olympus LUMFLN60XW, NA 1.1), which exhibited reduced coma aberration. Once the aberrated point spread function was found, the final refractive index reconstruction was carried out as described in [15], employing total variation regularization as described in [35].

2.4 Sample preparation

Formalin fixed-paraffin embedded pancreatic tissue slides were obtained after approval from the Institutional Review Board (IRB) with waiver of informed consent from patients [2021-1698]. The tissues were processed according to standard histopathological procedures except for staining process. Sections were cut to a thickness of 10 µm and 20 µm, respectively. Then, the samples were mounted on cover glass slides for imaging.

3. Results

3.1 Validation with a microsphere

Our methodology was initially verified using a synthetic reference sample, specifically an 8-µm-diameter Poly(Methyl Methacrylate) (PMMA) microbead. Figure 2(a) presents the RI tomograms of an 8-µm-diameter PMMA microbead embedded in UV-cure resin, before and after applying the aberration correction. The RI tomogram, after the aberration adjustment, shows a better restoration of the microbead’s geometric integrity and enhanced RI contrast. The tomogram post-aberration correction also more accurately reflects the bead’s expected spherical shape. Additionally, it demonstrates a substantial reduction in the ‘missing cone’ artifacts, recognizable as X-shaped distortions at the periphery and base of the bead in the vertical cross-section. This improvement suggests that applying total variation regularization to a more precise signal enhances the mitigation of these artifacts.

 figure: Fig. 2.

Fig. 2. Comparative analysis of a PMMA microbead, utilizing 8 illumination patterns. (a) RI tomograms of an 8-µm-diameter PMMA microbead embedded in UV-cure resin, before and after the aberration correction. The central inset displays the retrieved aberration. (b) RI profiles in axial (represented by slanted lines) and radial (represented by horizontal lines) directions, from the center-slices of the RI tomograms. (c) Experimental point spread functions (PSFs) obtained by deconvolving the RI tomograms with the known spherical shape.

Download Full Size | PDF

Aberration correction markedly improved both contrast and resolution, as evidenced by Figs. 2(b) and 2(c). In Fig. 2(b), the refractive index (RI) profiles in both axial and radial directions exhibit more pronounced distinctions between the PMMA microbead and its surrounding resin mounting medium following aberration correction. In contrast, the original measurements underestimated the RI values, leading to a blurred boundary, especially in the axial direction. Further evidence of resolution enhancement is detailed in Fig. 2(c), through the analysis of experimental point spread functions (PSFs). These PSFs, derived by deconvolving the RI tomograms using the known spherical shape of the bead (as introduced in Ref. [36]), illustrate that aberration correction notably diminished the axial elongation of the PSF. This finding unequivocally supports the improvement in axial resolution.

3.2 Application to histology samples

To evaluate our method's effectiveness in more complex samples with pronounced aberrations, we applied our aberration correction algorithm to pancreatic tissue sections, sliced to a thickness of 10 µm. Figure 3 highlights the substantial enhancement in image quality post-correction. In the initial tomogram, distortion artifacts due to aberrations significantly obscured the tissue's upper and lower boundaries, making them nearly indistinguishable. However, post-correction, these boundaries are sharply defined and accurately reflect the 10 µm thickness of the tissue.

 figure: Fig. 3.

Fig. 3. Enhanced visualization of pancreatic tissue via aberration correction, employing eight illumination patterns. (a) The left and right panels display the original and corrected tomograms of a selected region within pancreatic tissue, respectively. The aberration corrected RI tomograms emphasize the enhanced resolution and clarity achieved through aberration correction. The RI scale is provided, with detailed insets offering a magnified view (5 µm scale) of the tissue's microstructure. (b) The lower panels replicate this comparative format for a different pancreatic tissue section, with corresponding aberration patterns presented in the central images, underscoring the improvements in structural delineation post-correction.

Download Full Size | PDF

A notable enhancement is also observed in image contrast after correction. The RI values of the tissue relative to the background in the corrected image are nearly 50 percent higher than that in the original. This increase suggests that aberration in phase deconvolution microscopy can lead to destructive interference of the scattered signal, resulting in diminished intensity contrast and an underestimation of RI. Such underestimation has been reported in earlier studies [15], and it is likely that this phenomenon contributes to these findings. The correction process not only clarifies features in the vertical plane but also significantly improves the overall resolution of the image. Additionally, as evidenced in Fig. 3, the uniformity of aberrations across different tissue regions underscores the consistency and reliability of our aberration correction method.

3.3 Systematic analyses of aberration correction

To further evaluate our algorithm's performance under conditions of controlled aberration, we conducted experiments to determine its precision in detecting artificially induced spherical aberrations and its effectiveness in improving RI reconstruction. We deliberately introduced specific degrees of spherical aberration using the objective lens's correction collar, thereby mimicking potential real-world challenges. The correction collar settings ranged from 0.3 mm to 1.1 mm in target focus depth, and measurement where carried at a fixed axial position of the pancreatic tissue sample.

Figure 4 presents the algorithm's response to varying levels of spherical aberration, with each row representing a different correction collar setting. Although substantial spherical aberration initially compromised image quality, our algorithm effectively restored it. The aberrations identified by our algorithm corresponded closely with the theoretical radial phase profiles. For improved visualization, the 2D phase profile is presented in an unwrapped format. As anticipated, the 1st and 2nd order Zernike terms displayed a linear relationship with the levels of spherical aberration applied.

 figure: Fig. 4.

Fig. 4. Systematic analysis of aberration correction, utilizing eight illumination patterns. The outcomes of employing a correction collar to adjust for various spherical aberrations during RI reconstruction are shown. The leftmost panels display RI tomographic slices at progressive depths, marked with a 10 µm scalebar for reference. The central column exhibits the phase maps corresponding to different levels of spherical aberration, while the rightmost graphs plot the Zernike coefficients for the first and second-order spherical aberrations, offering a quantitative insight into the aberration retrieval process.

Download Full Size | PDF

Interestingly, the minimal aberration was noted at a setting of d = 0.7 mm. This finding is attributed to the disparity between the refractive index of the pancreatic tissue mounting medium and cover glass (1.5) and the 1.36 RI for which the lens's correction collar is calibrated, resulting in additional spherical aberration and is consistent with theoretical predictions (Supplement 1). For an in-depth explanation of how theoretical spherical aberrations are calculated based on adjustments to the correction collar, refer to Supplement 1.

We emphasize that our method is adept at correcting a wide range of aberrations, extending beyond merely spherical aberration. To substantiate this claim, we conducted further experiments using a lower quality objective lens (Olympus UPLSAPO60XW, NA 1.2), as depicted in Figure S5. This particular lens exhibits coma aberration, and notably, the optimal correction collar setting does not align with the cover glass thickness. Despite these challenges, our aberration correction algorithm successfully identified and corrected the aberration. It's pertinent to note that while modern objective lenses typically exhibit minimal aberrations, identifying the optimal correction collar position for each sample can be laborious. This issue is more pronounced with older, less expensive, or lenses designed for a larger field of view, where residual aberrations might persist.

4. Discussion and summary

In this study, we successfully developed and validated a method for aberration correction in 3D phase deconvolution microscopy. Demonstrated to be effective on both synthetic and biological thick specimens, our method markedly improved axial resolution and addressed the issue of RI underestimation. Precision in aberration retrieval was further confirmed through experiments with controlled spherical aberrations, where the retrieved aberration phase profiles closely aligned with theoretical expectations.

The present method offers several advantages in optical imaging systems. Firstly, it characterizes both system-induced and sample-induced aberrations, even when a sample is present. This is particularly beneficial for large histopathological samples, where aberrations may vary based on the sample's position, challenging conventional correction methods. Additionally, aberration removal can be performed post-sample data measurement, enhancing application potential in varied samples.

However, the method also faces challenges. It necessitates extra measurements for aberration retrieval, using redundant information. This can be mitigated by employing the method intermittently. Moreover, its iterative nature demands substantial computational resources. We have significantly improved computational efficiency by utilizing GPU capabilities for faster data processing. To further accelerate this process, we propose two strategies: First, aberrations, which are often consistent within a single sample as seen in pancreatic tissue imaging, can be retrieved once and applied across all fields-of-view for that sample. Second, for increased computational speed, aberrations can be identified in a smaller sub-volume and extrapolated to correct the entire field of view.

Our method also shows considerable promise for wide field-of-view microscopes, commonly employed in histology and medical imaging. Its applicability extends to volumetric histopathology, cytology, and the emerging field of 3D biology, including the study of organoids [3739]. This versatility underscores the method's promise in advancing microscopy techniques across a broad spectrum of biological research and diagnostic applications. These microscopes, typically featuring low-magnification and high numerical aperture lenses, are susceptible to spatially non-uniform aberrations, especially near the field's periphery. Adaptable to this challenge, our method could correct aberrations in individual sub-volume positions, effectively compensating for spatially varying aberrations [2,19].

Funding

Institute for Information and Communications Technology Promotion (2021-0-00745); National Research Foundation of Korea (2015R1A3A2066550, 2022M3H4A1A02074314); KAIST Institute of Technology Value Creation, Industry Liaison Center (G-CORE Project) (N10240002).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. J. Sasián, Introduction to Aberrations in Optical Imaging Systems (Cambridge University Press, 2013).

2. M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light: Sci. Appl. 3(4), e165 (2014). [CrossRef]  

3. K. Kim, W. S. Park, S. Na, et al., “Correlative three-dimensional fluorescence and refractive index tomography: bridging the gap between molecular specificity and quantitative bioimaging,” Biomed. Opt. Express 8(12), 5688–5697 (2017). [CrossRef]  

4. G. Kim, Y. Jo, H. Cho, et al., “Learning-based screening of hematologic disorders using quantitative phase imaging of individual red blood cells,” Biosens. Bioelectron. 123, 69–76 (2019). [CrossRef]  

5. Y. Jo, H. Cho, W. S. Park, et al., “Label-free multiplexed microtomography of endogenous subcellular dynamics using generalizable deep learning,” Nat. Cell Biol. 23(12), 1329–1337 (2021). [CrossRef]  

6. M. Esposito, C. Fang, K. C. Cook, et al., “TGF-β-induced DACT1 biomolecular condensates repress Wnt signalling to promote bone metastasis,” Nat. Cell Biol. 23(3), 257–267 (2021). [CrossRef]  

7. S.-A. Lee, L.-C. Chang, W. Jung, et al., “OASL phase condensation induces amyloid-like fibrillation of RIPK3 to promote virus-induced necroptosis,” Nat. Cell Biol. 25(1), 92–107 (2023). [CrossRef]  

8. Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” Nat. Photonics 12(10), 578–589 (2018). [CrossRef]  

9. T. L. Nguyen, S. Pradeep, R. L. Judson-Torres, et al., “Quantitative phase imaging: recent advances and expanding potential in biomedicine,” ACS nano 16(8), 11516–11544 (2022). [CrossRef]  

10. D. Jin, R. Zhou, Z. Yaqoob, et al., “Tomographic phase microscopy: principles and applications in bioimaging,” J. Opt. Soc. Am. B 34(5), B64–B77 (2017). [CrossRef]  

11. L. Tian and L. Waller, “Quantitative differential phase contrast imaging in an LED array microscope,” Opt. Express 23(9), 11394–11403 (2015). [CrossRef]  

12. M. H. Jenkins and T. K. Gaylord, “Three-dimensional quantitative phase imaging via tomographic deconvolution phase microscopy,” Appl. Opt. 54(31), 9213–9227 (2015). [CrossRef]  

13. M. Chen, L. Tian, and L. Waller, “3D differential phase contrast microscopy,” Biomed. Opt. Express 7(10), 3940–3950 (2016). [CrossRef]  

14. J. M. Soto, J. A. Rodrigo, and T. Alieva, “Label-free quantitative 3D tomographic imaging for partially coherent light microscopy,” Opt. Express 25(14), 15699–15712 (2017). [CrossRef]  

15. H. Hugonnet, M. Lee, and Y. Park, “Optimizing illumination in three-dimensional deconvolution microscopy for accurate refractive index tomography,” Opt. Express 29(5), 6293–6301 (2021). [CrossRef]  

16. H. Fujii and T. Asakura, “A contrast variation of image speckle intensity under illumination of partially coherent light,” Opt. Commun. 12(1), 32–38 (1974). [CrossRef]  

17. J. Li, A. Matlock, Y. Li, et al., “High-speed in vitro intensity diffraction tomography,” Adv. Photonics 1(06), 1–066004 (2019). [CrossRef]  

18. J. Park, S.-J. Shin, J. Shin, et al., “Quantification of structural heterogeneity in H&E stained clear cell renal cell carcinoma using refractive index tomography,” Biomed. Opt. Express 14(3), 1071–1081 (2023). [CrossRef]  

19. M. Chen, Z. F. Phillips, and L. Waller, “Quantitative differential phase contrast (DPC) microscopy with computational aberration correction,” Opt. Express 26(25), 32888–32899 (2018). [CrossRef]  

20. X. Ou, G. Zheng, and C. Yang, “Embedded pupil function recovery for Fourier ptychographic microscopy,” Opt. Express 22(5), 4960–4972 (2014). [CrossRef]  

21. G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photonics 7(9), 739–745 (2013). [CrossRef]  

22. T. Colomb, E. Cuche, F. Charrière, et al., “Automatic procedure for aberration compensation in digital holographic microscopy and applications to specimen shape compensation,” Appl. Opt. 45(5), 851–863 (2006). [CrossRef]  

23. I. Choi, K. Lee, and Y. Park, “Compensation of aberration in quantitative phase imaging using lateral shifting and spiral phase integration,” Opt. Express 25(24), 30771–30779 (2017). [CrossRef]  

24. Y. Baek, H. Hugonnet, and Y. Park, “Pupil-aberration calibration with controlled illumination for quantitative phase imaging,” Opt. Express 29(14), 22127–22135 (2021). [CrossRef]  

25. B. Seong, I. Kim, T. Moon, et al., “Untrained deep learning-based differential phase-contrast microscopy,” Opt. Lett. 48(13), 3607–3610 (2023). [CrossRef]  

26. T. Chang, D. Ryu, Y. Jo, et al., “Calibration-free quantitative phase imaging using data-driven aberration modeling,” Opt. Express 28(23), 34835–34847 (2020). [CrossRef]  

27. G. Gunjala, S. Sherwin, A. Shanker, et al., “Aberration recovery by imaging a weak diffuser,” Opt. Express 26(16), 21054–21068 (2018). [CrossRef]  

28. S. Lee, K. Lee, S. Shin, et al., “Generalized image deconvolution by exploiting the transmission matrix of an optical imaging system,” Sci. Rep. 7(1), 8961 (2017). [CrossRef]  

29. N. Ji, “Adaptive optical fluorescence microscopy,” Nat. Methods 14(4), 374–380 (2017). [CrossRef]  

30. N. Streibl, “Three-dimensional imaging by a microscope,” J. Opt. Soc. Am. A 2(2), 121–127 (1985). [CrossRef]  

31. H. Hugonnet, M. J. Lee, and Y. K. Park, “Quantitative phase and refractive index imaging of 3D objects via optical transfer function reshaping,” Opt. Express 30(8), 13802–13809 (2022). [CrossRef]  

32. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]  

33. Y. E. E. Nesterov, “A method of solving a convex programming problem with convergence rate O\bigl(k^2\bigr),” in Doklady Akademii Nauk (Russian Academy of Sciences, 1983), pp. 543–547.

34. K. B. Petersen and M. S. Pedersen, “The matrix cookbook,” Technical University of Denmark 7, 510 (2008).

35. A. Beck and M. Teboulle, “Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems,” IEEE Trans. on Image Process. 18(11), 2419–2434 (2009). [CrossRef]  

36. C. Park, S. Shin, and Y. Park, “Generalized quantification of three-dimensional resolution in optical diffraction tomography using the projection of maximal spatial bandwidths,” J. Opt. Soc. Am. A 35(11), 1891–1898 (2018). [CrossRef]  

37. H. Hugonnet, Y. W. Kim, M. Lee, et al., “Multiscale label-free volumetric holographic histopathology of thick-tissue slides with subcellular resolution,” Adv. Photonics 3(02), 026004 (2021). [CrossRef]  

38. Y. K. Lee, D. Ryu, S. Kim, et al., “Machine-learning-based diagnosis of thyroid fine-needle aspiration biopsy synergistically by Papanicolaou staining and refractive index distribution,” Sci. Rep. 13(1), 9847 (2023). [CrossRef]  

39. M. J. Lee, J. Lee, J. Ha, et al., “Long-term three-dimensional high-resolution imaging of live unlabeled small intestinal organoids using low-coherence holotomography,” bioRxiv, 2023.2009.2016.558039 (2023). [CrossRef]  

Supplementary Material (1)

NameDescription
Supplement 1       Supplementary document

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)

Fig. 1.
Fig. 1. (a) Schematic representation of the optical arrangement, detailing the path of light through the system components. (b) Aberration extraction process, where a gradient descent algorithm iteratively refines the aberration model to align with the correlation between simulated and empirical intensity profiles. (c) Enhancement sequence demonstrating aberration mitigation and the consequent recovery of the high-fidelity tomographic image.
Fig. 2.
Fig. 2. Comparative analysis of a PMMA microbead, utilizing 8 illumination patterns. (a) RI tomograms of an 8-µm-diameter PMMA microbead embedded in UV-cure resin, before and after the aberration correction. The central inset displays the retrieved aberration. (b) RI profiles in axial (represented by slanted lines) and radial (represented by horizontal lines) directions, from the center-slices of the RI tomograms. (c) Experimental point spread functions (PSFs) obtained by deconvolving the RI tomograms with the known spherical shape.
Fig. 3.
Fig. 3. Enhanced visualization of pancreatic tissue via aberration correction, employing eight illumination patterns. (a) The left and right panels display the original and corrected tomograms of a selected region within pancreatic tissue, respectively. The aberration corrected RI tomograms emphasize the enhanced resolution and clarity achieved through aberration correction. The RI scale is provided, with detailed insets offering a magnified view (5 µm scale) of the tissue's microstructure. (b) The lower panels replicate this comparative format for a different pancreatic tissue section, with corresponding aberration patterns presented in the central images, underscoring the improvements in structural delineation post-correction.
Fig. 4.
Fig. 4. Systematic analysis of aberration correction, utilizing eight illumination patterns. The outcomes of employing a correction collar to adjust for various spherical aberrations during RI reconstruction are shown. The leftmost panels display RI tomographic slices at progressive depths, marked with a 10 µm scalebar for reference. The central column exhibits the phase maps corresponding to different levels of spherical aberration, while the rightmost graphs plot the Zernike coefficients for the first and second-order spherical aberrations, offering a quantitative insight into the aberration retrieval process.

Equations (7)

Equations on this page are rendered with MathJax. Learn more.

$$\begin{array}{c} S = {F_{Re }} \ast {H_P} + {F_{{\mathop{\rm Im}\nolimits} }} \ast {H_A}\\ {H_A} ={-} 16{\pi ^2}Re [{{{\cal F}^{ - 1}}[{({\widetilde {{G_{\textrm{obj}}}}} )\ast ({{k_z}{{|\rho |}^2}\widetilde {{G^\ast }_{\textrm{obj}}}} )} ]} ]\\ {H_P} ={-} 16{\pi ^2}{\mathop{\rm Im}\nolimits} [{{{\cal F}^{ - 1}}[{({\widetilde {{G_{\textrm{obj}}}}} )\ast ({{k_z}{{|\rho |}^2}\widetilde {{G^\ast }_{\textrm{obj}}}} )} ]} ]\end{array}.$$
$$\widetilde {{G_{\textrm{obj}}}} = \widetilde {{P_{obj}}}\widetilde G,\textrm{ }\widetilde {{P_{obj}}} = \frac{{\sqrt {{k_z}} }}{k}\delta ({k_x^2 + k_y^2 < \textrm{N}{\textrm{A}^2}/{\lambda^2}} )\cdot {e^{iA({{k_x},{k_y}} )}}.$$
$$\begin{array}{c} {A_{optimal}} = \mathop {\arg \min }\limits_A \{{C(A )} \},\textrm{ }C(A )= \sum\limits_i {||{{\varepsilon_i}} ||_2^2} ,\textrm{ }{\varepsilon _i}(A )= {S_i} - {F_{Re }} \ast {H_{P,i}}(A )\\ {H_{P,i}}(A )={-} 16{\pi ^2}{\mathop{\rm Im}\nolimits} [{{{\cal F}^{ - 1}}[{({\widetilde {{G_{\textrm{obj}}}}(A)} )\ast ({{k_z}{{|{{\rho_i}} |}^2}\widetilde {{G_{\textrm{obj}}}^\ast }(A )} )} ]} ]\end{array}$$
$$\begin{array}{c} \widetilde {{G_{\textrm{obj}}}}(A) = \widetilde {{P_{obj}}}(A)\widetilde G,\textrm{ }\widetilde {{P_{obj}}}(A) = \frac{{\sqrt {{k_z}} }}{k}\delta (k_x^2 + k_y^2 < \textrm{N}{\textrm{A}^2}/{\lambda ^2}) \cdot {e^{iA({k_x},{k_y})}}\\ {F_{Re }} = {{\cal F}^{ - 1}}[\sum\limits_i {\widetilde {{S_i}}} {\widetilde {{H_{P,i}}}^\ast }(A)/\sum\limits_i {{{|{\widetilde {{H_{P,i}}}(A)} |}^2}} ] \end{array}$$
$$\scalebox{0.9}{$\begin{array}{l} {A_{optimal}} = \mathop {\arg \min }\limits_A \{ C(A)\} ,\textrm{ }C(A) = \sum\limits_i {||{\widetilde {{\varepsilon_i}}} ||_2^2} ,\textrm{ }\widetilde {{\varepsilon _i}}(A) = \widetilde {{S_i}} - \textrm{diag}(\widetilde {{H_{P,i}}})\widetilde {{F_{Re }}}\\ \widetilde {{H_{P,i}}}(A) = 16i{\pi ^2}U({\textrm{diag}({{U^{ - 1}}\widetilde {{G_{\textrm{obj}}}}} )U\textrm{diag}({k_z})\textrm{diag}({{|{{\rho_i}} |}^2}){{\widetilde {{G_{\textrm{obj}}}}}^\ast } - \textrm{diag}({U{{\widetilde {{G_{\textrm{obj}}}}}^\ast }} ){U^{ - 1}}\textrm{diag}({k_z})\textrm{diag}({{|{{\rho_i}} |}^2})\widetilde {{G_{\textrm{obj}}}}} )\\ \widetilde {{G_{ {\pm} \textrm{obj}}}}(A) = \textrm{diag}(\widetilde {{P_{obj}}})\widetilde G,\textrm{ }\widetilde {{P_{obj}}}(A,{k_x},{k_y}) = \frac{{\sqrt {{k_z}} }}{k}\delta (k_x^2 + k_y^2 < \textrm{N}{\textrm{A}^2}/{\lambda ^2}) \cdot {e^{iA({k_x},{k_y})}}\\ \widetilde {{F_{Re }}}(A) = {{\sum\limits_i {\textrm{diag}({{{\widetilde {{H_{P,i}}}}^\ast }} )\widetilde {{S_i}}} } / {\sum\limits_i {\textrm{diag}({{{|{\widetilde {{H_{P,i}}}} |}^2}} )} }} \end{array}$}$$
$$\begin{array}{c} {\nabla _A}C(A) = 2\textrm{Re} \left[ { - \sum\limits_j {{{\left( {\frac{{\partial \widetilde {{H_{P,j}}}}}{{\partial A}}} \right)}^H}\textrm{diag}({{{\widetilde {{F_{\textrm{Re} }}}}^\ast }} ){\varepsilon_j}} - {{\sum\limits_i {\left( {\frac{{\partial {{\widetilde {{H_{P,i}}}}^\ast }}}{{\partial A}}} \right)} }^H}\textrm{ diag}({{{\widetilde {{S_i}}}^\ast }} )\frac{{\sum\limits_j {\textrm{diag}({H_{P,j}}^\ast ){\varepsilon_j}} }}{{\sum\limits_k {\textrm{diag}({{{|{\widetilde {{H_{P,k}}}} |}^2}} )} }}} \right.\\ \left. { + \sum\limits_i {\left( {{{\left( {\frac{{\partial \widetilde {{H_{P,i}}}}}{{\partial A}}} \right)}^H}\textrm{diag}({\widetilde {{H_{P,i}}}} )+ {{\left( {\frac{{\partial {{\widetilde {{H_{P,i}}}}^\ast }}}{{\partial A}}} \right)}^H}\textrm{diag}({{{\widetilde {{H_{P,i}}}}^\ast }} )} \right)} \textrm{ diag}({{\widetilde {{F_{\textrm{Re} }}}}^\ast })\frac{{\sum\limits_j {\textrm{diag}({{\widetilde {{H_{P,j}}}}^\ast }){\varepsilon_j}} }}{{\sum\limits_k {\textrm{diag}({{{|{\widetilde {{H_{P,k}}}} |}^2}} )} }}} \right] \end{array}$$
$$\scalebox{0.7}{$\begin{aligned}&\left( {\displaystyle{{\partial \widetilde{{H_{P,j}}}} \over {\partial A}}} \right)^H = 32\pi ^2P_z{\mathop{\rm Re}\nolimits} \left[ {{\rm diag}\left( {\widetilde{{G_{{\rm obj}}}}} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right){\rm diag}\left( {k_z} \right)U^{-1}{\rm diag}\left( {U{\widetilde{{G_{{\rm obj}}}}}^*} \right)-{\rm diag}\left( {{\widetilde{{G_{{\rm obj}}}}}^*} \right)U{\rm diag}\left( {U^{-1}{\rm diag}\left( {k_z} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right)\widetilde{{G_{{\rm obj}}}}} \right)} \right]U^{-1}\\& \left( {\displaystyle{{\partial {\widetilde{{H_{P,j}}}}^*} \over {\partial A}}} \right)^H = 32\pi ^2P_z{\mathop{\rm Re}\nolimits} \left[ {-{\rm diag}\left( {{\widetilde{{G_{{\rm obj}}}}}^*} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right){\rm diag}\left( {k_z} \right)U{\rm diag}\left( {U^{-1}\widetilde{{G_{{\rm obj}}}}} \right) + {\rm diag}\left( {\widetilde{{G_{{\rm obj}}}}} \right)U^{-1}{\rm diag}\left( {U{\rm diag}\left( {k_z} \right){\rm diag}\left( {{\left| {\widetilde{{\rho _i}}} \right|}^2} \right){\widetilde{{G_{{\rm obj}}}}}^*} \right)} \right]U\end{aligned}$}$$
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.