Optical imaging featuring both long working distance and high spatial resolution by correcting the aberration of a large aperture lens

High-resolution optical imaging within thick objects has been a challenging task due to the short working distance of conventional high numerical aperture (NA) objective lenses. Lenses with a large physical diameter and thus a large aperture, such as microscope condenser lenses, can feature both a large NA and a long working distance. However, such lenses suffer from strong aberrations. To overcome this problem, we present a method to correct the aberrations of a transmission-mode imaging system that is composed of two condensers. The proposed method separately identifies and corrects aberrations of illumination and collection lenses of up to 1.2 NA by iteratively optimizing the total intensity of the synthetic aperture images in the forward and phase-conjugation processes. At a source wavelength of 785 nm, we demonstrated a spatial resolution of 372 nm at extremely long working distances of up to 1.6 mm, an order of magnitude improvement in comparison to conventional objective lenses. Our method of converting microscope condensers to high-quality objectives may facilitate increases in the imaging depths of super-resolution and expansion microscopes.

Scientific REPORtS | (2018) 8:9165 | DOI: 10.1038/s41598-018-27289-1 as mentioned above, the critical drawback of using lenses such as condensers for imaging is their strong aberration. Aberrations tend to vary steeply with respect to the spatial frequency, particularly at high frequencies, which significantly degrades spatial resolution. Thus, correcting the aberration of large aperture lenses is essential for the development of microscopes featuring both long working distance and high spatial resolution.
Various techniques have been proposed in the field of adaptive optical microscopy for the correction of both system and sample-induced aberrations. Adaptive optical methods make use of spatial light modulators or deformable mirrors to control the wavefront of illumination and/or detection beams. Depending on the type of aberration measurement, these methods can be classified into direct [25][26][27][28] and indirect [29][30][31][32][33] wavefront sensing schemes. These methods have proven to be effective in various applications. For example, the adaptive optics was applied to the super-resolution techniques such as STED and STORM to acquire extremely high-resolution images from samples inducing aberrations [34][35][36][37] . Most such schemes are effective at identifying sample-induced aberrations and designed to work at high speed using a relatively small number of aberration correction elements in wavefront shaping devices. Therefore, they are poorly suited to coping with steep variations in aberrations spanning the entire pupil because this will require a large number of aberration correction elements. We recently developed an aberration correction approach employing indirect wavefront sensing methods called closed-loop accumulation of single scattering (CLASS) microscopy working in the reflection geometry 38 . This method measures and corrects sample-induced aberration by maximizing the total intensity of synthetic aperture images. One important feature of CLASS microscopy is that the angular resolution of the aberration measurements can be extremely high because the method relies on finely angularly scanned images rather than feedback control by the wavefront shaping devices for the correction of aberrations. Thus, it is well suited to dealing with the steep angular variations in aberration associated with the lenses with a large diameter.
In this study, we constructed a transmission microscope composed of two commercial microscope condensers and applied the CLASS algorithm to identify and numerically correct the aberrations of the individual condensers. The innate aberrations of these lenses cause waves containing high spatial frequency components of the object spectrum to be inaccurately accumulated. Consequently, resolving power is attenuated, and PSF is significantly broadened in comparison with ideal diffraction-limit PSF. The CLASS algorithm was applied to correct such strong aberrations by maximizing the total intensity of the synthetic aperture image. By performing the maximization operation on both forward and phase-conjugated waves, the aberration maps of the two condensers could be distinguished. At a source wavelength of 785 nm, we demonstrated wide-field transmission imaging with a spatial resolution of 372 nm to the extremely long working distance of 1.6 mm. The performance of the proposed method was evaluated when imaging both a resolution target and yeast cells. Our method will help to improve the working depth of super-resolution microscopes, which require high NA objectives for maximal performance.

Experimental setup and data acquisition
We constructed a transmission-mode coherent imaging system based on a Mach-Zehnder interferometer (Fig. 1). To physically guarantee both long working distance and high spatial resolution, we employed two microscope condensers (MBL78700, Nikon; NA: 1.4, working distance: 1.6 mm) as illumination and imaging objective lenses. In this configuration, since the working distance of each condenser was 1.6 mm, the space between the two condensers where the sample could be inserted was 3.2 mm thick. A laser diode (LP785-SF100, Thorlabs; wavelength: 785 nm) was used as a light source. The output beam from the laser was split into the sample (red) and reference (blue) beams at the beam splitter (BS1), and the direction of illumination of the sample beam was controlled by a pair of galvanometer mirrors (GV011/M, Thorlabs). The sample beam was then introduced to the sample plane (SP) through an input condenser (OL i ), and the wave transmitted through the sample was captured by the output condenser (OL o ). The sample image was relayed to a scientific complementary metal oxide semiconductor (sCMOS) camera (pco.edge 4.2, PCO). The angle of the reference beam was tilted by selecting the first-order diffraction from a diffraction grating (G) with a spatial filter (SF), and then combined with the sample beam at the camera to acquire off-axis interference images at the camera. A 2D Hilbert transform was applied to the recorded interferograms to obtain the complex field maps, i.e., the phase and amplitude maps, of the transmitted waves 39,40 . The overall magnification from the sample plane to the camera was 80× and the view field was 50 × 50 µm 2 .
With a target in place, we took a set of complex field maps of transmitted waves, → E x y k ( , ; ) o i , by varying the illumination angle, or the transverse wavevector → k i of the incident wave. The → k i was chosen in such a way to cover all the orthogonal free modes in the given field of view and the numerical aperture. For the view field of 50 × 50 µm 2 and the numerical aperture of 1.2, the required number of images amounted to19,861 images. Typical image acquisition rate was 20 fps in the present study such that the total data acquisition takes approximately 15 minutes. The image acquisition time can potentially be reduced by using a high-speed camera and by minimizing the number of images required for the aberration correction. We took Fourier transform of each recorded image with respect to the spatial coordinate (x, y) and obtained the spatial frequency spectrum is the transverse wavevector of the transmitted wave conjugate to (x, y). We used the  to identify the system aberration and to reconstruct a target image (see details in the following sections).

Effect of aberration on synthetic aperture imaging
The condenser lens exhibits strong aberration, especially at large incidence angles. To describe the detrimental effects of this aberration on the formation of the image, consider a unit-amplitude plane wave at the sample plane, whose transverse wavevector is k k k ( , ) This plane wave propagates from the light source to the target object through the illumination condenser OL i . Aberration along this illumination pathway can be described by an input pupil function describes the phase retardation of the plane wave as a function of the transverse momentum → k . For an ideal lens, the pupil function P i will be unity for α | → | ≤ k k 0 , where k 0 is the magnitude of the wavevector in free space and α the NA of the lens, or zero otherwise 41 . In the presence of aberration, however, the plane wave may experience additional phase retardation depending on its propagation direction. Due to this aberration, the plane wave impinging on the pupil plane of OL i cannot form a clean focus at the sample plane, as indicated by the red curves in Fig. 2(a). Similarly, the optical transfer function from the target object through OL o to the detector can be described by an output pupil function can be regarded as the aberration induced by OL o . Due to this output aberration, the light scattered from a point target at the sample plane forms a distorted plane wave at the pupil plane of OL o (blue curves in Fig. 2(a)). Our analysis assumes no loss at the pupil planes throughout, i.e.
The effect of aberration from OL i to SP is relatively straightforward as the incident wavevector → k i does not change direction. On the other hand, the effect of the output aberration is more complicated since the incoming wave to SP experiences alteration of the wavevector due to the target object. Let us consider the target object whose object spectrum is written as  ∆ → k ( ). Due to the conservation of momentum, the transverse wavevector of the wave transmitted through the target object is modified to for each spatial frequency component ∆ → k of the target object. For each plane wave illumination, which is a Dirac delta function in Fourier space, the angular spectrum of the scattered field will simply be the translation of the object spectrum by → k i , which can be written as . Therefore, the angular spectrum of the transmitted wave is written as is the uncontrolled phase drift occurring during image acquisition. The phase drift originates from fluctuations in the optical path length difference between the sample waves and reference waves over the multiple measurements taken when scanning → k i . The synthetic aperture image can be computed by demodulating the target spectra from the illumination wavevector. In other words, the angular spectrum of the synthetic aperture image can be obtained by adding the electric fields of the transmitted waves that have the same momentum difference ∆ → k acquired for various Note that if the optics are free from aberration and uncontrolled drift, the angular spectrum of the synthetic aperture image will simply be proportional to . The term with summation is the cross-correlation of output pupil function P o a and the input pupil function multiplied by uncontrolled drift P ig function is complex-valued, compared to the ideal case the cross-correlation will have a smaller amplitude and stronger dependence on ∆ → k : The system aberration has two major detrimental effects on the synthetic aperture image: the absolute value of the cross-correlation becomes smaller, which causes a reduction in signal intensity, and the phase value of the cross-correlation becomes dependent on the momentum difference ∆ → k . This distorts the PSF of the optical system, thereby distorting the target structure in the reconstructed image.

Algorithm for aberration correction
Due to the inequality in Eq. (3), when optical system aberrations exist the total intensity of the synthetic aperture image becomes weaker than the aberration-free case. This concept is fundamental when correcting aberration in post-processing. We search for both the angle-dependent input correction function θ → k ( ) i i and the output correc- that maximize the total intensity of synthetic aperture image: Since it is difficult to find the correction functions θ θ and i o immediately, a tentative correction function θ that maximizes the total intensity of the synthetic aperture image while leaving θ = 0 o is applied. This step is equivalent to placing a wavefront shaping device at the input pupil function and controlling its wavefront so as to maximize the total intensity of the reconstructed image (the graph with black bars in Fig. 2(b)). The maximization operation represented by Eq. (4) is satisfied when the cross-term due to the absolute square is real. In other words, it is satisfied when the following scalar product of the two images taken at two arbitrary incident wavevectors  In the transmission side (blue curves), the ideal point source appears to be a distorted plane wave after transmitting through the output condenser lens. (b) By placing a virtual SLM at the Fourier plane to compensate for the additional phase shift at different angles (graph with black bars), the back aperture of the input condenser is illuminated by a distorted plane wave. Due to this initial phase shift, the light focused at the sample plane forms a clean PSF. (c) Here, the direction of propagation (optical axis) is reversed numerically by phase-conjugation. In a reversal of the previous process, a virtual SLM (graph with black bars on the right) is placed at the output pupil to focus light at the sample plane to compensate for output aberration.
Here the term However, the output pupil also adds aberration. Therefore, i m i n i is nonzero and represents the error of the aberration correction for the input, so additional independent steps are required to correct the output aberrations. In any case, this first round of correction serves to sharpen the PSF.
After finding the first-step input correction function θ i , the output aberration can be estimated in a similar way, but by reversing the process. As shown in Fig. 2(c), a virtual SLM (the graph with black bars on the right in Fig. 2(c)) can be placed at the output pupil plane and control its wavefront so as to maximize the intensity of the synthetic aperture image for the case when sending the waves from the opposite side. Note that for this operation the direction of the optical axis is inverted to use the same algorithm as for the input correction operation, but that additional data measurement is not required since applying phase conjugation to the original data set can supply the equivalent information. In the phase conjugation process, the complex fields in the momentum space can be written as After conducting same search process as before, the first-step correction to the output θ This correction process preferentially corrects the output aberration such that the PSF of the illumination from output side is made sharper. With the output correction in place, the input correction was repeated and the entire process performed iteratively. Then the accumulation of the correction functions converged to its corresponding aberration maps, i.e. θ φ ∑ → + g

Results
Test target imaging. The correction algorithm was demonstrated by imaging a custom-made Siemens star test target. The target was fabricated on a gold film coated on a slide glass. The target was covered by a second slide glass so that the total thickness of the sample became 2 mm. The configuration of measurement is depicted in Fig. 3(a). The inset in Fig. 3(b) shows an atomic force microscope (AFM) image of the finest element marked by a red square on the main image. The central 50 × 50 μm 2 region is imaged while scanning → k i for 19,861 different illumination angles up to 1.2 NA. The number of illumination angles depends on the size of the field of view, which is the inverse of the angular resolution, and the NA of the optical system. Figure 3(c) is a single-shot holographic image acquired with normal illumination. One can observe image distortion in the fine structures caused by the optical system's aberration. For example, Fig. 3 and applied our aberration correction algorithm to this sets of the measured images. The input correction functions (Fig. 3(e-h)) and the output correction functions (Fig. 3(j-l)) were acquired after 1, 3, 5, and 15 rounds of iteration, respectively. Note that the input correction functions appear to be meaningless because of the coupling between input aberration and uncontrolled drift. On the other hand, the output correction functions, which are decoupled from the uncontrolled phase drift, converge to spherically symmetrical phase maps. This is a reasonable outcome considering the spherically symmetric properties of the condenser lens. The synthetic aperture images corrected by each correction function are displayed in Fig. 3(m-p). As the iterations continue, fine target structures progressively resolve. To verify the convergence of the correction algorithm, we monitored the total intensity of the synthetic aperture image at each step ( Fig. 4(a)). The intensity converges to 3.5 × 10 16 (arbitrary units) after nine or ten iterations. The effect of aberration correction on the spatial resolution of the synthetic aperture images was investigated by observing the periodicity of a line profile. The image at the left-hand side of Fig. 4(b) shows the central 5 × 5 μm 2 area of a synthetic aperture image without aberration correction. In this image, only uncontrolled drift and the input aberration were compensated. The line profile along the white dashed line is presented on the right-hand side of Fig. 4(b). Its periodicity was 819 nm, which is equivalent to the diffraction limit of the aberration-free optical system with 0.4 NA objective lenses. After compensating for aberration, the finest structure of the target becomes resolvable as shown in the left-hand side of Fig. 4(c). Two line profiles along the white dashed line are obtained from the synthetic aperture images with (blue curve) and without (red curve) aberration correction. Note that the intensity of the line profile of the uncorrected image was multiplied by 10 3 for better visualization. Although there was slight dependence on the direction of the pattern, the structure of 372 nm periodicity was resolved with aberration correction. Thus, spatial resolution was enhanced from 819 nm to 372 nm. Considering that the diffraction-limit spatial resolution ∆ = = λ +

nm (NA N A )
ill c ol in our system, this result confirms that our aberration correction algorithm almost perfectly compensates for system aberration and allows microscope condensers to be used as objective lenses. The detrimental effect of the system aberration and the result of the aberration correction can be observed from the single shot images. Figure 4(d) is an uncorrected image captured for the normal illumination. While the third and fourth groups of patterns were resolvable, the second and the first groups remained unresolved. After applying the output correction function obtained from the CLASS algorithm, the second finest patterns became resolvable (Fig. 4(e)). Here, the input correction function needs not be applied because the sample was illuminated with a plane wave of a single illumination angle.
The difference between the images with and without the application of the correction function was more pronounced for the case of oblique illumination. In this case, the target was illuminated with a plane wave of a lateral momentum ≈ . k 1 0 0 . Without output correction, the single-shot image ( Fig. 4(f)) was more severely distorted than that of the normal illumination. Because of the oblique illumination, relatively large portion of the scattered light propagated through the outer region of the output pupil where the aberration was steep. However, the target structure was retrieved with the output correction function in place (Fig. 4(g)). These results support the reliability of the CLASS algorithm.
Yeast cell imaging. We applied the developed method to the imaging of biological cells. We prepared budding yeast Saccharomyces cerevisiae (S. cerevisiae) which was purchased from the Korean collection for type cultures (KCTC, Korea). The yeast cells were mixed with a 3% gelatin solution (G1393, Sigma) and then poured into a Petri dish (10090, SPL; bottom thickness ≈ 1 mm) to make a gelatin block of about 150 μm thickness at 4 °C. The block was covered by a slide glass, making a total sample thickness of 2.1 mm. This was too thick for an ordinary oil-immersion objective lens to reach the sample. The geometry of the sample is illustrated in Fig. 5(a).
First, the synthetic aperture image was computed without aberration correction (Fig. 5(b)). From the same dataset, the optical system's aberration including sample-induced aberration was measured by optimizing the total intensity of the synthetic aperture image. Figure 5(c) shows the aberration-corrected phase map of the yeast cells after ten rounds of iteration. The corresponding input and output correction functions are plotted in Fig. 5(d,e), respectively. Unlike the image reconstructed without aberration correction, in the aberration-corrected image granular structures within the cells were clearly visible (see white arrows in Fig. 5(b,c)). Figure 5(f) presents the phase retardation along the white dashed lines in each of the phase images, and shows that the spatial phase variation was sharpened by aberration correction.
The input correction function in Fig. 5(d), like the functions measured from the Siemens star target, has a complicated structure because of uncontrolled phase drift. On the other hand, the output correction function for the yeast cell sample shown in Fig. 5(e) has a different structure to that of the Siemens star dataset. This highlights the effect of the sample-induced aberration caused by the petri dish, which has a different refractive index from the glass. Separation of the input aberration and the uncontrolled drift. The input and output aberration correction maps identified by our algorithm are presented in Fig. 3(e-l) and Fig. 5(d-e). At first glance, the output correction maps were accurately measured, as judged by their circularly symmetrical structure. The condenser lens for detection presents strong aberration, especially at NA higher than 0.8, suggesting the likelihood of the pronounced degradation of spatial resolution. However, the input correction map presented rather irregular patterns and looked quite different from the output correction map despite the same model of condenser lens being used. In fact, this is due to uncontrolled phase drift during measurement. As discussed in Eq.
. The mixing of aberration and uncontrolled drift is not an issue for the present study since the correction of their summation is sufficient for the reconstruction of synthetic aperture images. But the separation of other scanning-type microscopes such as confocal or two-photon microscopes in such a way that our method provides aberration map to be corrected for the focus scanning microscopes.
This uncontrolled drift cannot be separated from input aberration without additional measurements. To address this issue, we propose a way to separate the input aberration from uncontrolled drift by measuring one more set of images, which contain the combined aberration of the two condensers. Firstly, the two correction functions are computed from the target images at the sample plane. The output correction function θ φ compensates for output aberration alone as shown in Fig. 6(a). In contrast, the input correction function corresponds to the sum of input aberration and uncontrolled drift. One more set of images was then taken by placing a test target at the conjugate image plane before the input condenser lens (Fig. 6(b)). In this configuration, we can consider the system as having no input aberration and a stronger output aberration than before as the two condensers' pupil functions are multiplied. Therefore, the complex field ′ SA for the synthetic aperture can be written as ) describes the angular spectrum of the test target. In this case, the correction functions of our aberration correction algorithm will yield θ′ , respectively. Therefore, the aberration of the input condenser φ i can be computed from Figure 6(c-e) show the experimental results of θ φ  To verify whether the input aberration map identified in this way is correct, we separated the three correction i on the images, and computed PSF at the sample plane ( Fig. 6(f)). In this case, PSF distortion is solely due to the aberration φ i caused by the input condenser. The correction for φ i was then additionally applied to form the PSF shown in Fig. 6(g). The resulting width of the PSF in Fig. 6(g) is much narrower than that in Fig. 6(f). From the line profiles along the white dashed lines (Fig. 6(h)), we could acquire the FWHMs of the uncorrected and corrected PSFs, which were 770 nm and 444 nm, respectively. The shape of the focus became an Airy function pattern after correcting for input aberration. We confirm this with a curve fitting of an Airy function model (yellow dashed line in Fig. 6(h)). The squared norm of the residual, which quantifies the variation of the data from the model, is 2.5 × 10 −3 . Note that the intensity of the uncorrected PSF was multiplied by ten for better visualization and that the peak intensity was enhanced 34.7 times by compensating for input aberration. Ordinary measurement configuration. The target is placed at the sample plane. The measured output aberration is equal to the aberration of the imaging condenser (marked by the red box). (b) A test target is placed at the image plane before the illumination condenser lens. The measured output aberration in this configuration corresponds to the combined aberration of the two condenser lenses. (c) Measured output aberration by placing the target at the sample plane, i.e., the configuration in (a). (d) Measured output aberration from a test target placed at the image plane before the illumination condenser, i.e. the configuration in (b). (e) Subtracted aberration map, (c) from (d). This corresponds to the input aberration originating from the illuminating condenser lens. (f) Intensity of PSF simulated from the transmission matrix with input aberration in place. To single out the effect of input aberration, only uncontrolled drift and output aberration are corrected numerically. (g) Intensity of PSF simulated from the aberration-free transmission matrix. Scale bar, 1 µm for both images. The two images are normalized to make the maximum intensity of aberration-free PSF unity. (h) Line profiles along the white dashed lines in the PSF images. The intensity of the PSF with input aberration is amplified 10-fold for better visualization. The PSF with input aberration correction is well described by an Airy function model (yellow dashed line). Note that the peak intensity is enhanced 34.7 times (inverse of Strehl ratio) by correcting for input aberration.

Discussion and Conclusion
In this study, we developed a method to measure and correct the aberration of an optical system in both its illumination and detection paths. By using this method, a pair of commercial microscope condensers turned into diffraction-limited objectives with an order-of-magnitude longer working distance than the conventional oil-immersion objective lenses. This allowed us to image specimens prepared between thick slide glasses with a spatial resolution of approximately 372 nm, a value close to the diffraction limit spatial resolution set by 1.2 NA at the wavelength of 785 nm. Furthermore, we separately identified the aberrations of the illumination and collection lenses. In particular, we could separate out the uncontrolled phase drift innate to coherent imaging from the input aberration to exclusively identify the aberration of the input condenser lens.
Relatively slow data acquisition speed is one of the current technical limitations. To capture a synthetic aperture image, one needs to capture multiple wide-field holographic images. The required number of images, although it depends on the field of view and the numerical aperture, is 19,861 in our measurements if one wants to sample the illumination angles fully. In our experiment, the camera frame rate was 20 fps such that it takes about 1,000 s for a set of measurements. Although this is not a critical issue for the present study where the characterization of the imaging system's aberration is a major concern, there is room to shorten the data acquisition time. Since there are plenty of photons to make use of in the transmission measurement, we can use a high-frame-rate camera to speed up the measurement. The number of images required for the aberration correction can also be reduced depending on the angular variation of aberrations, which will further enhance overall speed of the proposed method.
The proposed method ameliorates the image quality in the presence of steeply varying aberration at the pupil planes, which typically occurs with lenses of a large diameter, by the fine stepping of angular scanning and post-processing of the acquired set of coherent images. In comparison with the correction methods using feedback control and wavefront shaping devices 29,31 , the proposed method can be faster, especially when a large number of the correction modes matters. Even at the current implementation of slow image acquisition, we could correct 20 different angular modes per second. This is equivalent to correcting the same number of Zernike modes per second, which is high speed in the conventional adaptive optics. In the present study, we restricted our interest to deal with optical system's aberrations in transmission geometry, but the algorithm would be equally valid for correcting sample-induced aberrations as we have demonstrated in the reflection geometry. Since elastic scattering, not fluorescence, is used for aberration correction, the light intensity level is low enough to induce photo-bleaching. Therefore, our technique can be effectively combined with fluorescence-based imaging, by providing aberration correction maps and helping to increase the working distances at which high spatial resolution can be maintained.

Methods
Data processing. Custom scripts by Matlab R2016b (MathWorks Inc.) were used to CLASS optimization to the acquired set of images was conducted by using on a desktop computer (Intel Core i7-6700K CPU, 4.0 GHz, 64 GB RAM). A set of raw interference images were converted into complex fields with amplitude and phase by applying 2D Hilbert transformation. Hilbert transformation of 19,861 images takes 498 s. From these image sets, , which took 595 s. After that, the CLASS optimization process takes 440 s per iteration. The required number of iterations varies depending on the target objects. The iterative process was terminated by monitoring convergence of total synthetic aperture intensity. Including 15 iterations of CLASS optimization, it takes 8,289 s to reconstruct an aberration-free SA image from raw images. 2D Hilbert transform. 2D Hilbert transform is a process that reconstructs the phase and amplitude maps of a sample beam from an off-axis hologram. The transform consists of Fourier transform, low-pass filtering, and shift of the spectrum. To acquire an off-axis hologram, a reference beam of lateral spatial frequency → k R is introduced at the detector. The camera measures the intensity map, = + I E E S R 2 , where E S and E R are the sample and reference waves, respectively. The off-axis hologram can be decomposed into three different spectra centered at → 0 , → k R , and − → k R in the Fourier space. The first one is an autocorrelation of the sample complex field (Fourier transform of sample intensity). The second and the third terms are the spatial frequency spectrum of E S and its complex conjugate, respectively. The circular spectrum centered at → k R is selected by low-pass filtering with the bandwidth (the radius of the circular low-pass filter) set by the field of view and numerical aperture, and shifted by − → k R to eliminate the carrier frequency. By taking the inverse Fourier transform of the resulting spectrum, we could obtain the complex field map of E S , which completes the 2D Hilbert transform.