1.

Introduction

Uncorrected refractive error is the leading cause of vision impairment worldwide.^1–3 There is a need for inexpensive solutions, particularly in the developing world, which lacks sufficient resources and specialists. Furthermore, the prevalence of myopia has been rapidly increasing in recent years in Asia,⁴ putting a strain on eyecare systems. Wavefront sensors provide one potential avenue for automation to help address this problem. Unlike the phoropter,⁵ which is a complex apparatus that requires operation by a trained expert to determine the best subjective correction, a wavefront sensor-based device can directly measure the optical performance of the patient’s eyes. Conventional autorefractors⁶ used older technologies and were generally unsuccessful at automating this evaluation.

The advantage of the wavefront sensor is that it also provides measurements of higher-order aberrations (where low order refers to defocus and astigmatic error) that have been shown to be critical for determining an accurate correction.^7,8 Given these measurements, a more accurate correction may then be performed with conventional spectacle lenses by optimizing a sharpness metric that estimates the net effect on the retinal point spread function.⁹ Unfortunately, wavefront sensors still tend to be very sophisticated systems, primarily designed for expensive applications such as refractive surgery¹⁰ when used in ophthalmic practice; in such an application, the cost of errors is very high and clearly an expert operator must be presumed to be available. For the application of determining the needed refractive correction, however, we can potentially trade off much of the complexity and cost to address a much wider population. When the goal is correction of refractive errors, an estimate of the full wavefront error is not needed. It is primarily a subset of third- and fourth-order aberrations, particularly the more central terms in the Zernike expansion such as coma and spherical aberration, which dominate the effects of aberrations on visual acuity.^7,8 Furthermore, these are the largest high-order aberrations seen in the population.¹¹ A better reduction in complexity would, therefore, be achieved with a method that focuses on robustly estimating just these important terms.

The most common ophthalmic wavefront sensor utilizes a Shack–Hartmann lenslet array.^10,12 This is a relatively expensive device due to fabrication costs,¹³ particularly for high-density arrays, though prices might be reduced significantly given the development of high-volume fabrication techniques targeting the smartphone market.¹⁴ Indeed, higher density arrays are desirable to improve accuracy;¹⁵ lower density arrays suffer more from the curvature of the wavefront across each lenslet, as well as edge effects and related distortions caused by localized artifacts. A potentially cheaper high-density option¹⁵ is grating-based sensors, such as those that utilize the Talbot effect.^16,17 These can easily use grating patterns with a very small pitch (equivalent to a high-density array), which can be produced by simply etching the grating onto a glass slide attached to the detector.

In addition to the sensor device, wavefront sensing requires a specialized algorithm to estimate the wavefront from the detected intensity pattern. Figure 1 shows an overview of common methods utilized for Shack–Hartmann and Talbot devices. With higher density lenslet arrays (and with Talbot gratings), it becomes increasingly difficult to determine the displacements for individual elements to perform the classical methods of Fig. 1(a), due to diffraction and reduced signal-to-noise ratio,²⁰ which suggests Fourier methods may be more attractive. As depicted in Fig. 1(b), Fourier methods typically must unwrap the phase after demodulation of the grating frequency.²¹ Unfortunately phase unwrapping remains an open problem,²² as it is ill-posed and NP-hard;²³ even in a noise-free scenario, one cannot guarantee an optimal result without testing a combinatoric number of possible combinations of $2\pi$ steps. In practical terms, phase unwrapping algorithms are particularly sensitive to localized image artifacts, as are the spot detection methods of Fig. 1(a), and common issues such as eyelashes or corneal scars can cause catastrophic failures in reconstruction.²⁴ Figure 2 shows examples of eyelashes occluding the pupil, a common artifact that can result in significant difficulties. Research on phase unwrapping continues, utilizing sophisticated nonlinear optimization techniques such as total-variation denoising,²⁵ unscented Kalman filters,²⁶ and sparse modeling.²⁷

Fig. 1

(a) In classical Shack–Hartmann methods,¹⁸ focal spot displacements are estimated to get samples of the wavefront $x$- and $y$-gradients; the wavefront is computed from these gradients via an inversion algorithm, and Zernike coefficients can be computed from the reconstructed wavefront. (b) Fourier-based methods¹⁹ replace the spot detection step with a demodulation that directly extracts components containing the $x$- and $y$-gradients in their phase; the phase values are modulus $2\pi$, however, so phase unwrapping is required. (c) The proposed method operates directly on the demodulated signal components, by estimating the peak locations of products of signals; no phase unwrapping or inversion is needed.

Fig. 2

Examples of eyelash shadows partly occluding pupil signal; the bright regions are reflections of the laser source from the cornea and subsequent structures.

Researchers have also proposed techniques to improve or extend other components of the Fourier methods of Fig. 1(b). Demodulation can be performed by a fast filtering in the image domain.²⁸ Gradient inversion is itself a nontrivial inverse problem, and a variety of techniques have been developed²⁹ including fast noniterative approaches.³⁰ If the wavefront phase is small, it may also be extracted directly after demodulation in the spatial frequency domain, by taking the anti-hermitian component after centering,³¹ which can then be directly used in Fourier-based inversion methods. However, this approach requires symmetry of the signal aperture, and the small-angle approximation leads to increasing errors for larger aberrations.

As noted earlier, there are only a limited number of major aberration components of interest for refractive correction, so a more specialized approach may sidestep the major difficulties of these conventional reconstruction algorithms. This paper formulates a direct approach to estimate major aberrations, generally described in Fig. 1(c), which does not need phase unwrapping or iterative inversion. The approach uses a numerical technique mathematically similar to a spatial heterodyne, to create a diversity of images with different aberrations. At each polynomial order, the algorithm estimates each aberration coefficient by finding the peak frequency component. This technique is employed in a successive fashion to compute and correct for successively lower aberrations until only the low-order terms remain to be calculated much more accurately. Furthermore, the fact that this approach is not reliant on high-quality estimates of spot locations or pixel phase allows us to take advantage of higher pitch gratings, which reduces the relative effect of edges and localized artifacts, and also allows for a very compact and easily produced sensor. We will demonstrate the performance of the method with simulations depicting the robustness against severe artifacts, as well as over a range of realistic aberration magnitudes.

2.

Method

An ophthalmic wavefront sensor typically utilizes a monochromatic laser source scattered from a point on the subject’s retina. Viewing this as a point source, the scattered light passes through the eye’s optics in reverse, and deviations of the resulting signal from a planar wavefront provide an estimate of the eye’s optical aberrations. Most commonly, the returning signal is reimaged from the eye onto a wavefront sensor, such as in Fig. 3, which depicts a sensor utilizing the Talbot effect (a) and a Shack–Hartmann array (b). The resulting image at the camera may be described as the intensity of a interference pattern of the following equation:

Fig. 3

Wavefront sensor: (a) employing grating with pitch $P$ at appropriate distance from a camera to utilize Talbot effect and (b) employing lenslet array with pitch $P$ at appropriate distance from camera to produce focused spot pattern.

Fig. 4

Simulated grating (top row) and Shack–Hartmann (bottom row) result; detector image (left) and spatial-frequency image (right). Detector images are $1024\times 1024\text{pixels}$ covering 10 mm, with a pitch of $42\mu \mathrm{m}$ and a magnification of 2.0; here images are zoomed to central $200\times 200\text{pixel}$ region for visibility of pattern.

The two-dimensional Fourier transform of the detected image is

The frequency content of the image is depicted in Fig. 5. The important information can be retained with just the first-linear harmonic terms where $(n,m)=(\mathrm{1,0})$ and (0,1) denoted by the boxes in Fig. 5. For a plane wave, the peaks of the first terms will be located at $(0,\pm k_0)$ and $(\pm k_0,0)$, where $k_0=(1/P)$. Assuming the information is sufficiently concentrated near these points, we can extract the subimages depicted in Fig. 6, where the origins correspond to the $(0,k_0)$ and $(k_0,0)$ points, respectively, from the original spatial frequency image. By taking the inverse Fourier transform of these subimages, we get complex images modulated by the phase gradients

Fig. 5

Example of $\tilde{s}(k_x,k_y)$, two-dimensional spatial Fourier transform of sensor image, depicting locations of subimage selection about harmonic terms of interest.

Fig. 6

Subimages ${\tilde{s}}_{\alpha }$ and ${\tilde{s}}_{\beta }$ extracted from Fig. 5, depicting peak locations resulting from aberrations.

2.2.

High-Order Aberrations

To address high-order aberrations, we will take advantage of the fact that ophthalmic wavefronts are dominated by a small number of Zernike terms, while retaining the robustness and simplicity of Fourier-transform-based methods such as in Sec. 2.1. We address high-order terms by numerically computing shifted products of the demodulated images and their conjugates, and thereby convert higher-order terms to lower orders by a process similar to differentiation. This will allow us to subsequently utilize a low-order estimation method like that of Sec. 2.1. While the approach appears quite involved mathematically, the derivation only requires manipulations of small matrices, and the implementation consists of phase adjustments and fast Fourier transforms (FFTs) of the relatively small subimages around the peaks.

The first step is to process the subimages as follows, for both the $s_{\alpha }$ and $s_{\beta }$ subimages,

Eq. (10)

$$\begin{gathered} s_{\alpha }^{(x)}(x,y)=s_{\alpha }(x-\mathrm{\Delta },y)s_{\alpha }{(x+\mathrm{\Delta },y)}^{*}, \\ {s}_{\alpha }^{(xx)}(x,y)=s_{\alpha }^{(x)}(x-\mathrm{\Delta },y)s_{\alpha }^{(x)}{s}_{\alpha }{(x+\mathrm{\Delta },y)}^{*}, \\ {s}_{\alpha }^{(y)}(x,y)=s_{\alpha }(x,y-\mathrm{\Delta })s_{\alpha }{(x,y+\mathrm{\Delta })}^{*}, \\ {s}_{\alpha }^{(yy)}(x,y)=s_{\alpha }^{(y)}(x,y-\mathrm{\Delta })s_{\alpha }^{(y)}{(x,y+\mathrm{\Delta })}^{*}, \\ {s}_{\alpha }^{(xy)}(x,y)=s_{\alpha }^{(y)}(x-\mathrm{\Delta },y)s_{\alpha }^{(y)}{(x+\mathrm{\Delta },y)}^{*}. \end{gathered}$$

Recall that the demodulated subimages are complex, unlike the original detected image. We will label the phase of these subimages with $\alpha$ for those derived from the (1,0) term and $\beta$ for those derived from the (0,1) terms. For example, the phase of $s_{\alpha }^{(x)}$ is ${\alpha }^{(x)}$, and the phase of $s_{\beta }^{(x)}$ is ${\beta }^{(x)}$. The parameter $\mathrm{\Delta }$ is similar conceptually to the distance used in a finite-difference differentiation estimate (instead of first and second derivatives, we refer to them as first and second products). The precise choice is not critical, except that a large $\mathrm{\Delta }$ discards more information (as depicted in Fig. 7), while a smaller $\mathrm{\Delta }$ will result in a smaller frequency shift for our subsequent calculation, making the result more susceptible to noise and quantization errors. Section 3 gives a simulation demonstrating that the output is relatively insensitive to $\mathrm{\Delta }$ except for extreme values.

Fig. 7

An excessively large $\mathrm{\Delta }$ parameter reduces the amount of pupil information utilized (a) by $2\mathrm{\Delta }$ for the first product $s^{(y)}$ [and similar for $s^{(x)}$], (b) by $4\mathrm{\Delta }$ for the second products $s^{(yy)}$ [and similar for $s^{(xx)}$], and (c) by $2\mathrm{\Delta }$ in both directions for the remaining second product $s^{(xy)}$. So, a smaller $\mathrm{\Delta }$ is desirable to retain the most pupil area.

The result of such combinations is most easily understood in terms of transformations of monomials. Recall that a vector $\mathbf{z}$ of Zernike coefficients describing all aberrations up to some order can be described equivalently by a vector $\mathbf{p}$ of coefficients for monomial terms.³² The vector $\mathbf{p}$ is related to $\mathbf{z}$ by a simple linear transformation (i.e., $\mathbf{z}={\mathbf{M}}_{m2z}\mathbf{p}$). A fourth-order wavefront can be described as the monomial expansion

Eq. (11)

$$w(x,y)=\sum _{i=0}^4\sum _{j=0}^{4-i}{p}_{i,j}{x}^i{y}^j,$$

for which the coefficients may be easily visualized in the matrix form (using base-zero indexing)

Eq. (12)

$$\mathbf{P}=\left(\begin{array}{ccccc}{p}_{\mathrm{0,0}}& p_{\mathrm{0,1}}& p_{\mathrm{0,2}}& p_{\mathrm{0,3}}& p_{\mathrm{0,4}}\\ p_{\mathrm{1,0}}& p_{\mathrm{1,1}}& p_{\mathrm{1,2}}& p_{\mathrm{1,3}}& 0\\ p_{\mathrm{2,0}}& p_{\mathrm{2,1}}& p_{\mathrm{2,2}}& 0& 0\\ p_{\mathrm{3,0}}& p_{\mathrm{3,1}}& 0& 0& 0\\ p_{\mathrm{4,0}}& 0& 0& 0& 0\end{array}\right),$$

with $\mathbf{p}=\mathrm{vec}(\mathbf{P})$ using the vectorization operator, an operator that reforms matrix elements into a vector (purely for convenience of notation). In matrix form, separable linear operations may be performed using a matrix $\mathbf{M}$, such as ${\mathbf{M}}^T\mathbf{P}$ to transform the $x$-coordinate and $\mathbf{PM}$ to transform the $y$-coordinate. See Ref. 32 for a more thorough introduction.

The coefficients of the gradient wavefronts $\alpha (x,y)$ and $\beta (x,y)$ can be computed using a differentiation matrix ${\mathbf{M}}_D$, to get

Eq. (13)

$${\mathbf{P}}_{\alpha }={\mathbf{M}}_D^T\mathbf{P}=\left(\begin{array}{ccccc}{p}_{\mathrm{1,0}}& p_{\mathrm{1,1}}& p_{\mathrm{1,2}}& p_{\mathrm{1,3}}& 0\\ 2p_{\mathrm{2,0}}& 2p_{\mathrm{2,1}}& 2p_{\mathrm{2,2}}& 0& 0\\ 3p_{\mathrm{3,0}}& 3p_{\mathrm{3,1}}& 0& 0& 0\\ 4p_{\mathrm{4,0}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),$$

Eq. (14)

$${\mathbf{P}}_{\beta }={\mathbf{PM}}_D=\left(\begin{array}{ccccc}{p}_{\mathrm{0,1}}& 2p_{\mathrm{0,2}}& 3p_{\mathrm{0,3}}& 4p_{\mathrm{0,4}}& 0\\ p_{\mathrm{1,1}}& 2p_{\mathrm{1,2}}& 3p_{\mathrm{1,3}}& 0& 0\\ p_{\mathrm{2,1}}& 2p_{\mathrm{2,2}}& 0& 0& 0\\ p_{\mathrm{3,1}}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right).$$

In this section, we will assume the device constants $P=d=1$ to simplify the notation; the scale factor $P/d$ may be applied when done, the same as it appears in Sec. 2.1 for low-order coefficient estimates. Note that the locations (row and column, counting from zero) in the matrix ${\mathbf{P}}_{\alpha }$ and ${\mathbf{P}}_{\beta }$ give the power of the new term via an expansion similar to Eq. (11) for $\alpha (x,y)$ and $\beta (x,y)$, respectively. Hence we can easily see that the wavefronts are now third order, as would be expected from differentiation, since the highest row or column index (when counting from zero) for a nonzero element is three.

To determine the effect of the shifts by $\mathrm{\Delta }$, we use translation matrices of the form

Eq. (15)

$${\mathbf{M}}_{\mathrm{\Delta }}=\left(\begin{array}{cccc}1& 0& 0& 0\\ \mathrm{\Delta }& 1& 0& 0\\ {\mathrm{\Delta }}^2& 2\mathrm{\Delta }& 1& 0\\ {\mathrm{\Delta }}^3& 3{\mathrm{\Delta }}^2& 3\mathrm{\Delta }& 1\end{array}\right).$$

So, e.g., to determine the monomial coefficients for ${\alpha }^{(x)}(x,y)=\alpha (x-\mathrm{\Delta },y)-\alpha (x+\mathrm{\Delta },y)$ given ${\mathbf{P}}_{\alpha }$ we would simply compute ${\mathbf{M}}_{-\mathrm{\Delta }}^T{\mathbf{P}}_{\alpha }-{\mathbf{M}}_{\mathrm{\Delta }}^T{\mathbf{P}}_{\alpha }$. The resulting coefficients are

Eq. (16)

$${\mathbf{P}}_{{\alpha }^{(x)}}=\left(\begin{array}{cccc}4\mathrm{\Delta }{p}_{\mathrm{2,0}}+8{\mathrm{\Delta }}^3{p}_{\mathrm{4,0}}& 4\mathrm{\Delta }{p}_{\mathrm{2,1}}& 4\mathrm{\Delta }{p}_{\mathrm{2,2}}& 0\\ 12\mathrm{\Delta }{p}_{\mathrm{3,0}}& 12\mathrm{\Delta }{p}_{\mathrm{3,1}}& 0& 0\\ 24\mathrm{\Delta }{p}_{\mathrm{4,0}}& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$

A similar computation for ${\alpha }^{(xx)}(x,y)={\alpha }^{(x)}(x+\mathrm{\Delta },y)-{\alpha }^{(x)}(x-\mathrm{\Delta },y)$ gives

Eq. (17)

$${\mathbf{P}}_{{\alpha }^{(xx)}}=\left(\begin{array}{cccc}24{\mathrm{\Delta }}^2{p}_{\mathrm{3,0}}& 24{\mathrm{\Delta }}^2{p}_{\mathrm{3,1}}& 0& 0\\ 96{\mathrm{\Delta }}^2{p}_{\mathrm{4,0}}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$$

The wavefront described by these coefficients is linear. Writing out the polynomial gives

Eq. (18)

$$\begin{gathered} {\alpha }^{(xx)}(x,y)=\sum _{i=0}^4\sum _{j=0}^{4-i}{[{\mathbf{P}}_{{\alpha }^{(xx)}}]}_{i,j}{x}^i{y}^j \\ =24{\mathrm{\Delta }}^2{p}_{\mathrm{3,0}}+96{\mathrm{\Delta }}^2{p}_{\mathrm{4,0}}x+24{\mathrm{\Delta }}^2{p}_{\mathrm{3,1}}y. \end{gathered}$$

Hence, the processed image will be of the form (where, again, we assume $P=d=1$)

Eq. (19)

$$\begin{gathered} s_{\alpha }^{(xx)}(x,y)=C_{\mathrm{1,0}}\exp [-2\pi i{\alpha }^{(xx)}(x,y)] \\ =C_{\mathrm{1,0}}\exp [-2\pi i(96{\mathrm{\Delta }}^2{p}_{\mathrm{4,0}}x+24{\mathrm{\Delta }}^2{p}_{\mathrm{3,1}}y)], \end{gathered}$$

and its two-dimensional Fourier transform will have a peak at a location

Eq. (20)

$$(k_x,k_y)=(96{\mathrm{\Delta }}^2{p}_{\mathrm{4,0}},24{\mathrm{\Delta }}^2{p}_{\mathrm{3,1}}).$$

So by finding the location of this peak, we can estimate the $p_{\mathrm{4,0}}$ and $p_{\mathrm{3,1}}$ coefficients of $w(x,y)$. With similar logic applied to the other processed images, we are able to estimate all fourth-order terms, to get

Eq. (21)

$$\begin{gathered} p_{\mathrm{4,0}}=\frac{1}{96{\mathrm{\Delta }}^2}X[s_{\alpha }^{(xx)}], \\ {p}_{\mathrm{3,1}}=\frac{1}{72{\mathrm{\Delta }}^2}\{Y[s_{\alpha }^{(xx)}]+X[s_{\alpha }^{(xy)}]+X[s_{\beta }^{(xx)}]\}, \\ {p}_{\mathrm{2,2}}=\frac{1}{64{\mathrm{\Delta }}^2}\{Y[s_{\alpha }^{(xy)}]+X[s_{\alpha }^{(yy)}]+Y[s_{\beta }^{(xx)}]+X[s_{\beta }^{(xy)}]\}, \\ {p}_{\mathrm{1,3}}=\frac{1}{72{\mathrm{\Delta }}^2}\{Y[s_{\alpha }^{(yy)}]+Y[s_{\beta }^{(xy)}]+X[s_{\beta }^{(yy)}]\}, \\ {p}_{\mathrm{0,4}}=\frac{1}{96{\mathrm{\Delta }}^2}Y[s_{\beta }^{(yy)}], \end{gathered}$$

where the terms that appear in multiple images are averaged, and the coordinates for the peak of the Fourier transform of an image $s$ are denoted as $[X(s),Y(s)]$.

The next step is to remove these components from the lower order images, which we achieve by applying a phase adjustment based on the quantity of the term that exists in each (see Appendix). For example, from Eq. (16), we see that there are quadratic components resulting from the fourth-order terms of $w$, namely, $24\mathrm{\Delta }{p}_{\mathrm{4,0}}$ in the (2,0) entry, $12\mathrm{\Delta }{p}_{\mathrm{3,1}}$ in the (1,1) entry, and $4\mathrm{\Delta }{p}_{\mathrm{2,2}}$ in the (0,2) entry. As a result we would apply

Eq. (22)

$${\delta }_{\alpha }^{(x)}(x,y)=\exp [2\pi i\mathrm{\Delta }(24p_{\mathrm{4,0}}{x}^2+12p_{\mathrm{3,1}}xy+4p_{\mathrm{2,2}}{y}^2)]$$

to $s_{\alpha }^{(x)}(x,y)$. Similarly for the other three quadratic images, we would apply

Eq. (23)

$$\begin{gathered} {\delta }_{\alpha }^{(y)}(x,y)=\exp [2\pi i\mathrm{\Delta }(6p_{\mathrm{3,1}}{x}^2+8p_{\mathrm{2,2}}xy+6p_{\mathrm{1,3}}{y}^2)], \\ {\delta }_{\beta }^{(x)}(x,y)=\exp [2\pi i\mathrm{\Delta }(6p_{\mathrm{3,1}}{x}^2+8p_{\mathrm{2,2}}xy+6p_{\mathrm{1,3}}{y}^2)], \\ {\delta }_{\beta }^{(y)}(x,y)=\exp [2\pi i\mathrm{\Delta }(4p_{\mathrm{2,2}}{x}^2+12p_{\mathrm{1,3}}xy+24p_{\mathrm{0,4}}{y}^2)]. \end{gathered}$$

Now we can calculate the cubic terms from $w(x,y)$ using the peaks of the Fourier transform of these corrected images, as

Eq. (24)

$$\begin{gathered} p_{\mathrm{3,0}}=\frac{1}{12\mathrm{\Delta }}X[{\delta }_{\alpha }^{(x)}{s}_{\alpha }^{(x)}], \\ {p}_{\mathrm{2,1}}=\frac{1}{12\mathrm{\Delta }}\{Y[{\delta }_{\alpha }^{(x)}{s}_{\alpha }^{(x)}]+X[{\delta }_{\alpha }^{(y)}{s}_{\alpha }^{(y)}]+X[{\delta }_{\beta }^{(x)}{s}_{\beta }^{(x)}]\}, \\ {p}_{\mathrm{1,2}}=\frac{1}{12\mathrm{\Delta }}\{Y[{\delta }_{\beta }^{(x)}{s}_{\beta }^{(x)}]+X[{\delta }_{\beta }^{(y)}{s}_{\beta }^{(y)}]+Y[{\delta }_{\alpha }^{(y)}{s}_{\alpha }^{(y)}]\}, \\ {p}_{\mathrm{0,3}}=\frac{1}{12\mathrm{\Delta }}Y[{\delta }_{\beta }^{(y)}{s}_{\beta }^{(y)}]. \end{gathered}$$

Finally, we correct the third- and fourth-order terms in the original gradient subimages, to accurately compute the second-order aberrations, by forming

Eq. (25)

$$\begin{gathered} {\delta }_{\alpha }(x,y)=\exp [2\pi i\mathrm{\Delta }(4p_{\mathrm{4,0}}{x}^3+3p_{\mathrm{3,1}}{x}^{2}y+2p_{\mathrm{2,2}}x{y}^2 \\ +p_{\mathrm{1,3}}{y}^3+3p_{\mathrm{3,0}}{x}^2+2p_{\mathrm{2,1}}xy+p_{\mathrm{1,2}}{y}^2)], \\ {\delta }_{\beta }(x,y)=\exp [2\pi i\mathrm{\Delta }(p_{\mathrm{3,1}}{x}^3+2p_{\mathrm{2,2}}{x}^{2}y+3p_{\mathrm{1,3}}x{y}^2 \\ +4p_{\mathrm{0,4}}{y}^3+p_{\mathrm{2,1}}{x}^2+2p_{\mathrm{1,2}}xy+3p_{\mathrm{0,3}}{y}^2)]. \end{gathered}$$

With these, we can calculate

Eq. (26)

$$\begin{gathered} p_{\mathrm{2,0}}=\frac{1}{2}X({\delta }_{\alpha }{s}_{\alpha }), \\ {p}_{\mathrm{1,1}}=\frac{1}{2}[Y({\delta }_{\alpha }{s}_{\alpha })+X({\delta }_{\beta }{s}_{\beta })], \\ {p}_{\mathrm{0,2}}=\frac{1}{2}Y({\delta }_{\beta }{s}_{\beta }). \end{gathered}$$

Finally, we form the vector of monomial coefficients that can be transformed to Zernike coefficients for application

Eq. (27)

$$\mathbf{p}={(p_{\mathrm{2,0}},p_{\mathrm{1,1}},p_{\mathrm{0,2}},p_{\mathrm{3,0}},p_{\mathrm{2,1}},p_{\mathrm{1,2}},p_{\mathrm{0,3}},p_{\mathrm{4,0}},p_{\mathrm{3,1}},p_{\mathrm{2,2}},p_{\mathrm{1,3}},p_{\mathrm{0,4}})}^T.$$

To summarize, we have the following algorithm:

• 1. Demodulate detected image, producing complex subimages: $s_{\alpha }$ and $s_{\beta }$.

• 2. Compute first product images: $s_{\alpha }^{(x)}$, $s_{\alpha }^{(y)}$, $s_{\beta }^{(x)}$, $s_{\beta }^{(y)}$.

• 3. From first product images, compute second product images: $s_{\alpha }^{(xx)}$, $s_{\alpha }^{(xy)}$, $s_{\alpha }^{(yy)}$, $s_{\beta }^{(xx)}$, $s_{\beta }^{(xy)}$, $s_{\beta }^{(yy)}$.

• 4. From FFT peaks of second product images, compute fourth-order coefficients: $p_{\mathrm{4,0}}$, $p_{\mathrm{3,1}}$, $p_{\mathrm{2,2}}$, $p_{\mathrm{1,3}}$, $p_{\mathrm{0,4}}$.

• 5. Using fourth-order coefficients, generate phase corrections: ${\delta }_{\alpha }^{(x)}$, ${\delta }_{\alpha }^{(y)}$, ${\delta }_{\beta }^{(x)}$, ${\delta }_{\beta }^{(y)}$.

• 6. From FFT peaks of first product images with phase corrections applied, compute third-order coefficients: $p_{\mathrm{3,0}}$, $p_{\mathrm{2,1}}$, $p_{\mathrm{1,2}}$, $p_{\mathrm{0,3}}$.

• 7. Using third-order coefficients, generate phase corrections: ${\delta }_{\alpha }$, ${\delta }_{\beta }$.

• 8. From FFT peaks of original subimages with phase corrections applied, compute second-order coefficients: $p_{\mathrm{2,0}}$, $p_{\mathrm{1,1}}$, $p_{\mathrm{0,2}}$.

• 9. Transform monomial coefficients to Zernike coefficients: $\mathbf{z}={\mathbf{M}}_{m2z}\mathbf{p}$.

3.

Results

Next, the method was tested with a variety of simulations of a Talbot sensor. The simulations used a pupil size of 5 mm with a detected image size of 10 mm for $1024\times 1024\text{pixels}$. The grating pitch was $42\mu \mathrm{m}$, the wavelength was 850 nm, and the magnification was $m_s=2.0$ (this simply required that a scaling be applied to the final terms based on their order). Figure 8 shows the estimates versus true coefficients for the example from Fig. 4. Here, the result of the simple second-order estimation of Eq. (9) is also demonstrated, which does not use the subsequent high-order estimation and their correction. Figure 9 shows the Fourier transforms of the subimages calculated in the method, for one of the two linear harmonic terms used [specifically for just the $(n,m)=(\mathrm{1,0})$ term]. First note that the peak in (a), the original subimage, exhibits a large amount of high order, as do the subimages in (b) and (c), computed from a single-shifted product. However, the second-shifted product images in (d), (e), and (f) result in Airy disks, as there is no remaining high order at this point. Similarly, the correction of the previous subimages, depicted in (g), (h), and (i), result in Airy disks in the Fourier transforms, as the high order has been accurately estimated and removed. Figures 10 and 11 show results for the same case with severe artifacts intentionally added to the image; again, the high-order coefficients are accurately estimated.

Fig. 8

True Zernike coefficients used in simulation, and result of algorithm to estimate fourth-order Zernikes (Fast4), as well as a version that stops after estimating second-order coefficients (Fast2). Net rms error (the root of the sum of the errors for each coefficient squared) is $0.043\mu \mathrm{m}$-rms for Fast4, and $2.65\mu \mathrm{m}$-rms for Fast2.

Fig. 9

Example of Fourier transforms of subimages used in processing, showing aberrations of peak: (a) ${\tilde{s}}_{\alpha }$, (b) ${\tilde{s}}_{\alpha }^{(x)}$, (c) ${\tilde{s}}_{\alpha }^{(y)}$, (d) ${\tilde{s}}_{\alpha }^{(xx)}$, (e) ${\tilde{s}}_{\alpha }^{(xy)}$, (f) ${\tilde{s}}_{\alpha }^{(yy)}$, (g) ${\delta }_{\alpha }^{(x)}{\tilde{s}}_{\alpha }^{(x)}$, (h) ${\delta }_{\alpha }^{(y)}{\tilde{s}}_{\alpha }^{(y)}$, and (i) ${\delta }_{\alpha }{\tilde{s}}_{\alpha }$.

Fig. 10

Image artifact examples, with random spots or lines to demonstrate algorithm robustness.

Fig. 11

Resulting estimates for the three artifact examples. Net rms errors (the root of the sum of the errors for each coefficient squared) for Fast4 are 0.056, 0.18, and $0.19\mu \mathrm{m}$-rms. As compared to an error of $0.043\mu \mathrm{m}$-rms for the example with no artifacts. The net errors for Fast2 are 2.58, 2.56, and $2.74\mu \mathrm{m}$-rms for the three cases, comparable to the $2.65\mu \mathrm{m}$-rms error in the artifact-free case.

The previous examples simulate a purely fourth-order wavefront, meaning no fifth- or higher-order aberrations are present. This presumption is the basis for the simplifications we were able to make. To test the effect of realistic levels of higher-order aberrations, a large set of images were simulated based on a secondary analysis of real sixth-order Zernike coefficients for measurements from 1500 subjects collected at multiple clinics in North America and Asia utilizing multiple different aberrometers. Figure 12 shows the standard deviations for each Zernike coefficient at 5 mm diameter versus the residual (true minus estimate), which generally agrees with the statistics found in other studies.¹¹ We see that the magnitude of the terms drops off quickly after fourth order (corresponding to coefficients above 14). Typically, values in this region are $<0.02\mu \mathrm{m}$-rms, which is comparable to the repeatability of the measurement due to biological variation. In Fig. 12, we see that the error after correction is reduced to this range as well.

Fig. 12

Standard deviations of coefficients for each Zernike term in simulated data (top) and in residual error after estimate (bottom). Standard deviations of true Zernike coefficients 3, 4, and 5, are 0.27, 1.40, and 0.4274, respectively.

Finally, this large simulation was used to test the parameter choice $\mathrm{\Delta }$. This parameter directly results in the shift of the peaks within the processed subimages, hence larger $\mathrm{\Delta }$ makes this estimate more accurate. Though note also that interpolation of the Fourier transform, peak locations will generally be needed to accurately measure small coefficients. The choice of larger $\mathrm{\Delta }$ also comes at a price, however, of utilizing less data from the subimages as described in Fig. 7. In our initial testing, a moderate choice of seven pixels was sufficient, and the result was not very sensitive to this choice until extremes were reached. This agrees with the simulation results shown in Fig. 13, where error was found to be small and relatively constant between the extremes of small shifts and a maximum of about 13.

Fig. 13

Average wavefront error versus shift parameter $\mathrm{\Delta }$ for large population of simulated wavefronts. For this pupil size and sample spacing, the result is relatively insensitive to choices of $\mathrm{\Delta }$ in the range from three to twelve.