Aberration correction method based on double-helix point spread function

Abstract. Point spread function (PSF) engineering has met with lots of interest in various optical imaging techniques, including super-resolution microscopy, microparticle tracking, and extended depth-of-field microscopy. The intensity distributions of the modified PSFs often suffer from deteriorations caused by system aberrations, which greatly degrade the image contrast, resolution, or localization precision. We present an aberration correction method using a spiral-phase-based double-helix PSF as an aberration indicator, which is sensitive and quantitatively correlated to the spherical aberration, coma, and astigmatism. Superior to the routine iteration-based correction methods, the presented approach is iteration-free and the aberration coefficients can be directly calculated with the measured parameters, relieving the computing burden. The validity of the method is verified by both examining the intensity distribution of the conventional Gaussian PSF in three dimensions and observing muntjac skin fibroblast cells. This iteration-free correction method has a potential application in PSF engineering systems equipped with a spatial light modulator.


Introduction
Point spread function (PSF) engineering has been a routine technique in various applications, such as super-resolution microscopy, 1 microparticle tracking, 2 and extended depth-offield microscopy. 3,4 Apart from the conventional imaging systems, PSFs are regularly modified into a variety of shapes, including biplane PSF, 5 Bessel PSF, 6 cubic-phase PSF, 7,8 astigmatic PSF, 9 and double-helix PSF 10,11 (DH-PSF), by which improved imaging resolution, axial localization precision, or depth-of-field can be achieved. PSFs of imaging systems are often distorted by system aberrations, which stem from misaligned optical elements, surface-shape error of the optical components, and refractive-index mismatch. [12][13][14][15] As a consequence, imaging performance inevitably degrades, making aberration correction an essential procedure in these applications.
Generally, the premeasured aberration can be compensated by introducing a conjugated phase with adaptive optical elements, such as a spatial light modulator (SLM). Existing aberration measurement techniques can be categorized into direct wavefront sensing and indirect wavefront sensing. 16,17 Direct sensing is usually implemented with a Shack-Hartmann wavefront sensor, [18][19][20][21] which can provide high speed measurement; however, its requirement of both a wavefront sensor and a correction element increases the cost and complexity of the system. 17 In contrast, indirect sensing methods only require a wavefront correction element, which makes them more economical and convenient for the construction of microscopes. [22][23][24][25] Typically, images of fluorescent beads at different focusing depths are used to retrieve, the aberration, by transferring phases and amplitudes iteratively between the focal plane and the pupil plane. 22,26 However, the precision of the indirect sensing method is limited by the employed Gaussian PSF, due to its weak response to wavefront error. Superior to Gaussian PSF, optical vortices with symmetrical doughnut shapes are more sensitive to phase irregularity. Slight wavefront distortion will destroy its rotationally symmetric intensity distribution. For this reason, Jesacher et al. 27 utilized the doughnut PSF to determine the astigmatism caused by the flatness deviation of SLM; however, this approach was irresponsive for spherical aberration. In addition, the above indirect sensing methods are all modified from iteration-based algorithms, which are time-consuming and may not converge sometimes. [28][29][30] To circumvent the above problems, we propose an iterationfree aberration correction method with high sensitivity by employing a newly proposed double-helix PSF 31,32 as the aberration indicator. Double-helix PSF is originally generated with the superposition of Gaussian-Laguerre (GL) modes, 33 which is proved to be of high efficiency and has become the mainstream generation method in related applications. 10,34,35 The influence of aberrations on the double-helix PSF has been characterized by former researchers. 36,37 Nevertheless, none of them reveal the potential of the double-helix PSF in aberration correction, because the shape of aberration-free double-helix PSF based on GL modes somehow deviates from the ideal double-helix shape, which makes it difficult to depict the distortions. 31 Recently, an approach to generating the double-helix PSF based on spiral phase 38 was developed and proved to be more efficient and flexible. 32 In Ref. 31, we employed this doublehelix PSF to build a single-shot three-dimensional (3-D) fluorescent microscope and surprisingly found that the doublehelix PSF based on spiral phase is closer to an ideal doublehelix shape for its uniform rotation rate and interlobe distance, with which the distortions of the PSF can be easily quantified. Therefore, in this paper, we further adopt the spiral-phase-based double-helix PSF as the indicator of aberrations. The proposed double-helix PSF is sensitive to not only the coma and the astigmatism, but also the spherical aberration. 31 Moreover, the coefficients of the aberrations can be directly calculated by analyzing the shape of measured double-helix PSFs without the iterative procedure, which relieves the computing burden and takes no convergence risk.
As a demonstration, we corrected the primary spherical aberration, primary x-coma, primary y-coma, primary 0-deg astigmatism, and primary 45-deg astigmatism of a high numerical aperture fluorescent microscope with the developed approach. The intensity distribution of the corrected 3-D Gaussian PSF was measured to be in good agreement with well-known aberration-free distribution. In addition, we also observed muntjac skin fibroblast cells with the identical system, of which the brightness and the image contrast were significantly improved after the correction. We anticipate that the iteration-free aberration correction method will have a potential application in PSF engineering systems.

Methods
The computer-generated holograms (CGH) to generate the spiral-phase-based double-helix PSF can be simply described by the equation below E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 4 3 2 P N ðρ; θÞ ¼ where ρ represents the normalized radius of the selected aperture, P N ðρ; θÞ denotes the modulated phase at the polar coordinates ðρ; θÞ, and N is the total number of the Fresnel zones. The phase mask is composed of a sequence of radial sampling of the spiral phase into the Fresnel zones. The rotation rate and the interlobe distance of the PSF can be easily adjusted by changing the number of Fresnel zones. 31 Without loss of generality, we select N as 6 for all simulations. The complex amplitude of the beam at the exit pupil plane is assumed to be E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 2 6 0 U ex ðρ; θÞ ¼ E in ðρ; θ; z ¼ 0Þ expfi½P N ðρ; θÞ þ Wðρ; θÞg; (2) where E in ðρ; θ; z ¼ 0Þ denotes the incident wave field distribution at the exit pupil plane, Pðρ; θÞ represents the modulated phase to generate the double-helix PSF, and Wðρ; θÞ is the wavefront aberration of the system. For simplicity, the incident wave is treated as uniform plane wave, i.e., let E in ðρ; θ; z ¼ 0Þ ¼ 1.
The propagation from the exit pupil plane to the observation plane is described by Fresnel-Kirchhoff diffraction theory, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 1 2 9 U o ðρ 0 ; θ 0 ; zÞ ∝ where k represents the wave number and z is the distance between the two planes.
To keep consistent with the circular pupil of a routine optical system and the designed phase mask, we use Zernike circle polynomials to characterize the optical aberrations. For simplicity, we only consider primary spherical aberration, coma, and astigmatism here. Higher order aberrations can be analyzed in the same manner. Field curvature and distortion are out of our scope for they have no influence on the shape of the PSF. 39 Therefore, the wavefront aberration here can be described as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 6 4 2 where A s , A cx , A cy , A ax , and A ay are the corresponding coefficients of primary spherical aberration, primary x-coma, primary y-coma, primary 0-deg astigmatism, and primary 45-deg astigmatism in units of wavelength λ. Figure 1 shows the double-helix PSFs affected by different types of aberrations. We can draw a sketchy conclusion that spherical aberration causes focal plane shift and generates a nonuniform rotation rate, coma unbalances the energy of the main lobes, and astigmatism brings about varying interlobe distances at different depths. 31 In fact, not only the spherical aberration but also the astigmatism affects the rotation rate. In addition, astigmatism and coma in different directions will affect the PSF in different ways, as shown in Figs. 1(c) and 1(d), respectively. In the following section, quantitative results will be presented in detail. We confine the coefficient of every single aberration within 0.6 λ.

Quantitative Analysis of the Point Spread
Function Deformation Figure 2 shows the quantitative deformation of the double-helix PSFs with varying spherical aberrations. To measure the rotation angle and the interlobe distance, positions where the two main lobes reach their maximum intensity are first obtained, which can be denoted as ðx 1 ; y 1 Þ and ðx 2 ; y 2 Þ (supposing y 2 > y 1 ), respectively. Then, the rotation angle can be directly calculated as the angle deviation of the attachment of the two points from the y-axis, that is, arctan½ðx 2 − x 1 Þ∕ðy 2 − y 1 Þ. Similarly, the interlobe distance can also be measured by calculating the distance of the two maximum points. As shown in Fig. 2(a), the focal plane where the rotation angle reaches 0-deg shifts with an increase in the spherical aberration, and the uniformity of the rotation rate is severely broken as A s gets >0.4 λ. Meanwhile, the interlobe distance gradually deviates from the nearly uniform curve with an increase in the spherical aberration [ Fig. 2(b)]. Although the distorted interlobe distance curve cannot be approximated to a familiar function, its monotonicity could be helpful for the aberration correction, which will be discussed in Sec. 2.2.1. All deformations in the PSF become opposite as we choose an opposite A s . The quantitative deformation of the double-helix PSFs with varying 0-deg astigmatisms is shown in Figs. 3(a) and 3(b), respectively. Similarly, the focal plane shifts with the increasing 0-deg astigmatism and the rotation rate also turns into nonuniform. But the changing tendency is quite different. As the aberration coefficient increases, the rotation rate curve gradually  deviates from the original linear curve while keeping constant at two symmetric positions. The difference of the deviated curve from the linear curve can be approximated to a cosine function. The interlobe distance gradually turns nonuniform and the variation resembles a sinusoidal function, whose amplitude increases monotonously with the increasing coefficient. A similar conclusion can be drawn for 45-deg astigmatism. The constant point of the rotation rate curve turns to the focal plane, which implies the focal plane does not shift with increasing 45-deg astigmatism. The deviation of the rotation rate from the linear curve becomes a sinusoidal function, while the interlobe distance resembles a cosine function. Figure 4 shows how coma affects the shape of the doublehelix PSF. As shown in Figs. 4(a) and 4(c), neither the x-coma nor the y-coma impacts the uniformity of the rotation rate. In fact, not only the rotation rate, but also the interlobe distance is almost unaffected by coma. The interlobe distance curves are not provided here for brevity. Then, only the energy unbalance of the main lobes is involved for the increasing coma [ Fig. 1(d)]. To describe the extent of unbalance, we define a specific parameter that is calculated by the natural logarithm value of the lobe-contrast (LVC), as depicted by the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 2 0 4 where E 1 and E 2 , respectively, denote the energy of two main lobes, which are calculated by integrating the intensity of the corresponding main lobe in the region where the intensity is brighter than half the maximal intensity of the lobe. The LVC curves of x-coma and y-coma are, respectively, shown in Figs. 4(b) and 4(d), which can be approximated to sinusoidal curves and cosinoidal curves. The amplitudes of both curves increase with the absolute value of the coefficients A cx and A cy , respectively.

Single aberration correction
Based on the analysis above, straightforward correction can be readily carried out if the system only contains a single aberration. As shown in Fig. 2, the difference in the interlobe distance where the rotation angle reaches −90 deg and 90 deg monotonously increases with the increasing spherical aberration. Therefore, with premeasured mapping relation between the difference and the spherical aberration, the coefficient of the spherical aberration can be directly measured by subtracting the interlobe distances at the both ends.
Correction of astigmatism and coma is a bit more complicated, for both of them consist of two components. As mentioned above, the interlobe distance curves of 0-deg astigmatism and 45-deg astigmatism, respectively, resemble a sinusoidal function and a cosinoidal function with identical period. Only when the aberration of the system consists of either 0-deg astigmatism or 45-deg astigmatism can it be directly calculated by measuring the amplitude of the sinusoidal or cosinoidal function. Otherwise, the crosstalk between 0-deg astigmatism and 45-deg astigmatism needs to be considered. To investigate this issue, we simulate the combined effect of the two astigmatisms, as is shown in Fig. 5(a). No distinct crosstalk is found, i.e., the combined effect of the 0-deg astigmatism and the 45-deg astigmatism on the interlobe distance is approximate to the superposition of a sine function and a cosinoidal function with corresponding amplitudes. In other words, the mathematical relation of the combined effect can be described as the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 1 5 5 where p represents the period of the periodic interlobe distance. C x and C y , respectively, denote the weights of the 0-deg astigmatism and the 45-deg astigmatism, which are to be measured. By experimentally measuring the interlobe distance curve of the aberrated double-helix PSF, the weights C x and C y can be directly calculated. Further, the coefficient of the 0-deg astigmatism and the 45-deg astigmatism can be obtained according to the premeasured mapping relations. For x-coma and y-coma, the LVC curves are also approximated to a sinusoidal and a cosinoidal function with identical period, respectively. The combined effect of the coma in two directions is also investigated [see Fig. 5(b)], verifying the absence of the crosstalk between the coma in two directions. Therefore, the coefficients of x-coma and y-coma can be measured in a similar approach to astigmatism, and the mathematical relation in Eq. (6) is also applicable for coma.

Multiple aberration correction
As demonstrated above, three shape parameters of the doublehelix PSF, i.e., rotation rate, interlobe distance, and LVC, are significantly influenced by the aberrations. In turn, the variations of these parameters can provide clues for the determination of the aberration coefficients. Both of the former two are influenced by spherical aberration and astigmatism, while the LVC is only affected by coma. However, as shown in Figs. 2-4, variation in the rotation rate is not as remarkable and regular as the latter two. Thus, in the workflow below, we only use the interlobe distance to determine the spherical aberration and astigmatism. The flowchart of the aberration correction is shown in Fig. 6. The entire process includes three procedures and each procedure corrects one type of aberration. At the beginning of each procedure, the double-helix PSF is experimentally remeasured to minimize the disturbance among different types of aberrations.
Spherical aberration is corrected first. From the above interlobe distance curve in Figs. 2(b), 3(b), and 3(d), we found the discrepancy of the interlobe distance where the rotation angle equals −90 deg and 90 deg only varies with spherical aberration and increases monotonously with the coefficient A s . Thus, spherical aberration can be evaluated by remapping the measured difference to a premeasured lookup table. It is worth noting that only the axial points at both ends are used in this step, which means the measurement speed can be improved by capturing fewer images. Correction of astigmatism follows after the spherical aberration is corrected. The combined effect of the 0-deg astigmatism and the 45-deg astigmatism on the interlobe distance has been investigated in Sec. 2.2.1, which is essentially the overlay of a sine function and a cosine function with different amplitudes, as described in Eq. (6). Thus, the weights of each aberration C x and C y could be obtained with the remeasured interlobe distance curve. By remapping the weights to the simulative results, the coefficients of the 0-deg astigmatism and the 45-deg astigmatism are directly obtained. Similarly, coefficients of x-coma and y-coma can be acquired with the remeasured LVC curve.

Experimental Results and Discussion
As a demonstration, correction of the primary aberrations is carried out in an SLM-based adaptive fluorescence microscope, which is shown in Fig. 7. The 3-D PSF of the system is measured by axially scanning a commercial fluorescent nanoparticle with a diameter of 100 nm (F8803, Thermo Fisher Scientific Inc.) through the interested range, which is illuminated by the collimated beam from a solid-state laser (λ ¼ 491 nm, Calypso 491, Cobolt AB Inc., Sweden). The emission light is collected by the identical microscopy objective (100×, NA ¼ 1.25, Nikon Inc., Japan), of which the back focal plane is reimaged onto the SLM (1920 × 1080 pixels, Pluto II, HoloEye Photonics AG, Germany) with 1:1. The modulated light is then filtered by a long-pass emission filter (λ c ¼ 500 nm, Thorlabs Inc.) and collected by the zoom lens, and the specimen is eventually  imaged onto a scientific CMOS camera (Orca Flash4.0, 2048 × 2048 pixels, 16bits, Hamamatsu Inc., Japan). The double-helix PSF is generated by addressing the corresponding CGH (N ¼ 6) on the liquid-crystal panel. Meanwhile, the measured aberrations are compensated by the conjugate phaseerror onto the CGH. The axial scanning interval to capture the double-helix PSF is set as 100 nm. Even though the optical components are elaborately calibrated, primary aberrations including spherical aberration, astigmatism, and coma can be observed by the measured curves in Figs. 6(b) and 6(c), respectively. These aberrations will be corrected in sequence with the developed method. First, the spherical aberration coefficient is easily measured as −0.4 λ by calculating the difference of the interlobe distance at both ends of the curve in Fig. 6(b). The spherical aberration can be directly corrected by appending the measured conjugate phase-error on the CGH. The effect of the correction proved to be fairly good by Fig. 6(d), in which the remeasured difference value of the interlobe distance at both ends is close to zero. With the identical interlobe distance curve [ Fig. 6(d)], weights of the 0-deg astigmatism or 45-deg astigmatism can be, respectively, measured by fitting the measured curve to sinusoidal function. A ax and A ay are, respectively, calculated as 0 and 0.25 λ by remapping the calculated weights to the coefficients. As a consequence, the interlobe distance curve of the PSF after correction becomes uniform as is shown in Fig. 6(f), indicating the astigmatism is well corrected. Similarly, the coma coefficients A ax and A ay are, respectively, measured as −0.08 and 0 λ by fitting the measured LVC curve in Fig. 6(g). In the end, after all the primary aberrations are corrected, both the interlobe distance and the LVC turn out to be approximately uniform, which matches with that of the simulated aberration-free double-helix PSF. In total, the aberration of this system mainly consists of −0.4 λ spherical aberration, 0.25 λ 45-deg astigmatism, and −0.08 λ x-coma. The spherical aberration is mostly caused by the refractive-index mismatch of the immersion medium, whereas the astigmatism and the coma collectively stem from the misalignment of the system and the surface curvature of the SLM's reflective panel.
The above correction process takes about 4 s, which mostly consists of the acquisition time of the double-helix PSFs and the computing time of the employed curves. As the total number of acquired images is 64 (two for the first step; 31 for the second step, and the third step), the total acquisition time for the correction sums up to be 3.2 s with a exposure time of 50 ms. With an illumination power of 10 mW, the maximal grayscale level of the acquired double-helix PSF at the focal plane is about 3000 (full gray level is 65,535 with the sCMOS), which is adequate for our correction process. Further reduction in the acquisition time can be achieved by lowering the exposure time of the camera and the axial sampling rate of the PSF. In addition, it also takes about 0.8 s to calculate and fit the interlobe distance curve and the LVC curve in the second and the third step, respectively.
To confirm the validity of the correction, the 3-D Gaussian PSFs before and after correction are measured, as shown in  Fig. 8(a). Three axial positions are provided here. The slight asymmetry of the defocused spots implies the subtle astigmatism and coma of the system, while significant asymmetry below and above the focal plane reveals the severe spherical aberration. As expected, the symmetry of the PSF greatly improves after the correction, verifying the good performance of the correction method. To provide a more intuitional comparison, the rendered 3-D view of the PSF is shown in Fig. 8(b). In summary, the primary aberrations of the system are properly corrected with the developed method. Next, a section of muntjac skin fibroblast cells (F36925, Thermo Fisher Scientific Inc.) is tested with epi-illumination in the identical system. The prominent F-actin and mitochondria of the cell are, respectively, labeled with Alexa Fluor 488 phalloidin and an anti-OxPhos Complex V inhibitor protein mouse monoclonal antibody in conjunction with Alexa Fluor 555 goat antimouse IgG. Both of them are excited by the laser and eventually imaged by the camera. Improvement in the image brightness after correction can be instinctively observed by comparing the observation results in Figs. 9(a) and 9(b), respectively. To provide a clearer comparison, magnified views of the two images are also shown in Figs. 9(c) and 9(d). The mitochondria after correction are clearly discernible, while that in the original image remains dark and blurry. The intensity profiles along the solid lines in Figs. 9(c) and 9(d) are shown in Fig. 9(e) to make a quantified comparison of the image quality, in which the brightness of the image is calculated to be improved by a factor of 1.8. Owning to the degraded signal-to-noise rate (SNR) and the distorted PSF, only three of the five marked peaks are recognizable in the original image, while all of them are clear and sharp after correction. Collectively, brightness and contrast of the acquired image are significantly improved after the correction procedure.
It is not difficult to notice the irregular fluctuations in the curves of Fig. 6, which inevitably impacts the measurement precision. Such fluctuations are mainly attributed to the random noise in the measured PSFs that are captured in low light dose. For the fluorescent beads in our experiment, improving the illumination power will be helpful at the very start, but the subsequent images will suffer from poorer SNR because of the increasing photobleaching effect. Thus, the more effective approach to improve the signal level may be to adopt the brighter and antifade subresolution point source, such as quantum dots. In addition, as the difference value used in the first step is calculated by the data at only two points, the measured coefficient of spherical aberration will be more susceptible to the noise. Thus, it is possible to improve the performance of the correction by conducting a second run of the first step in some experiments.

Summary
We have quantitatively studied the response of the double-helix PSF to different aberrations. Both the spherical aberration and the astigmatism lead to the nonuniform rotation rate and interlobe distance, while the coma aberration only unbalances the energy of the main lobes. Crosstalk between different aberrations is also investigated, based on which we develop an aberration correction method. Superior to the routine iteration-based correction methods, this approach is iteration-free and the aberration coefficient can be directly calculated with the measured parameter, which is time-saving and takes no convergence risk. The validity of the promoted approach is verified by examining the 3-D Gaussian PSF after correction, which matches pretty well with the well-known aberration-free intensity distribution. The developed approach might provide an easy and economical tool to improve the performance in PSF engineering systems equipped with an SLM.

Disclosures
The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.