Inverse optical design of the human eye using likelihood methods and wavefront sensing

We are developing a method for estimating patient-specific ocular parameters, including surface curvatures, conic constants, tilts, decentrations, thicknesses, refractive indices, and index gradients. The data consist of the raw detector outputs from one or more Shack-Hartmann wavefront sensors, and the parameters in the eye model are estimated by maximizing the likelihood. A Gaussian noise model is used to emulate electronic noise, so maximum likelihood reduces to nonlinear least-squares fitting between the data and the output of our optical design program. The Fisher information matrix for the Gaussian model was explored to compute bounds on the variance of the estimates for different system configurations. In this preliminary study, an accurate estimate of a chosen subset of ocular parameters was obtained using a custom search algorithm and a nearby starting point to avoid local minima in the complex likelihood surface.


Introduction
The availability of a complete patient-specific eye model would provide many key advantages in vision science and ophthalmology. For example, the model could be used as a theoretical basis for vision correction of higher-order aberrations in laser refractive surgery or with corrective lenses, since classical devices only measure and correct for basic refractive errors (i.e., defocus and astigmatism). Accurate determination of optical parameters in normal and abnormal eyes could also be valuable in developing data bases for clinical diagnosis of pathologies, while measurements of ocular surface misalignments would be useful after implantation of intraocular lenses in cataract surgery [1][2][3]. Moreover, knowledge of the refractive index distribution in the lens may be beneficial in optical coherence tomography, which is based on interferometric reflectometry and index changes [4,5]. It could additionally lead to a substantial improvement in retinal imaging. For instance, current adaptive-optics ophthalmoscopes incorporating a Shack-Hartmann wavefront sensor (WFS) and wavefront corrector conjugated to a single surface of the eye offer high resolution [6,7], but over a very limited field-of-view (FOV) [8] due to a form of anisoplanatism involving aberrations of the eye. Aberrations collected over different field positions on the retina result from the passage through different parts of the ocular media so that the adaptive-optics correction is valid only over a certain field area, referred to as the isoplanatic patch. One solution is to conjugate multiple wavefront sensors and correctors to various refractive surfaces in the eye, thereby increasing the isoplanatic patch size and enabling wide-field measurements, but choice of the optimal planes at which to conjugate the correctors would be facilitated by knowing the real eye structure of the individual.
In addition to improvements in vision correction and retinal imaging, the availability of patientspecific parameters could facilitate a broad range of ongoing vision science studies. Of significant interest is the in vivo gradient-index (GRIN) distribution and lenticular geometry of the human crystalline lens as a function of both age and accommodation [3], but this information has been difficult to obtain, and reliable measurements are scarce [4,5,[8][9][10]. While previous studies suggested that aspheric surfaces in the anterior segment and an effective refractive index for the lens are sufficient to model spherical aberration, lack of knowledge regarding the GRIN distribution precludes both the prediction of off-axis aberrations and study of dispersion in the lens, so that experimental data are limited [9]. A complete mapping of the human eye could also be used to evaluate intersubject variability and statistical variations, as well as vision performance and image quality in the central and peripheral visual fields [11,12], which could be enhanced by accurate measurement of the retinal curvature [13,14]. Another fundamental study in physiological optics is how individual ocular components factor into the overall performance of the human eye [15] and how such performance would change if one or more surfaces are altered, a critical element in surgical procedures. While schematic eyes have been extremely useful for that purpose, they often lack asymmetry such as decentration of the lens or pupil, which manifest in the fovea as aberrations of nonaxisymmetrical systems (e.g., coma, astigmatism, and transverse chromatic aberration) and may have a significant impact on ocular performance [16,17]. We believe that a patient-specific mapping of the entire eye, including nonaxially symmetric components, would enable further investigations that have been previously unapproachable.
We are developing a new method for studying the human eye, either for clinical ophthalmology or for basic research. In our approach, we estimate patient-specific parameters of the eye, including surface curvatures, conic constants, tilts, decentrations, thicknesses, refractive indices, and graded-index distributions. The data consist of the raw detector outputs from a Shack-Hartmann wavefront sensor (WFS), a device that measures the phase distortions of a wavefront and provides useful information about aberrations in an optical system. Once we have the data, the parameters in the eye model are estimated via maximum-likelihood (ML) estimation using an optical design program. This process is the opposite of traditional optical design, where the lens parameters at each iteration are fixed and data on one or more output planes are computed. Here the WFS data are given, and the lens parameters are estimated, so we refer to the method as inverse optical design. It should be noted that this is a nonlinear estimation problem, since the data depend nonlinearly on the parameters. Similar methods have been applied in astronomy and optical shop testing [18,19]. In astronomy, image data are used to estimate optical prescription parameters (i.e., first-order geometrical parameters), a process referred to as prescription retrieval. In optical shop testing, coefficients in an arbitrary wavefront expansion (e.g., Zernike coefficients) are estimated to characterize aberrations in a test element. However, this is a linear estimation problem in contrast to our method.
Our method of inverse optical design using wavefront sensing could provide an in vivo, noninvasive, and complete mapping of the human eye, including dozens of parameters that are essential to an accurate representation of the eye and its aberrations. Existing in vivo methods supply a small subset of ocular parameters. For example, a common technique in phakometry uses Purkinje images of the back reflections from the anterior and posterior surfaces of both the cornea and crystalline lens, providing basic curvatures, tilts, and decentrations [2,3]. However, one difficulty in this approach is that insufficient knowledge of the refractive index distribution of the lens leads to significant measurement errors in the lens posterior radius [20]. Scheimpflug slit imaging is increasingly being used to obtain sharp cross-sectional images of the anterior eye segment, imparting surface shapes, misalignments, and intraocular distances, although accurate determination of these parameters relies on the correction of optical distortions in the imaging system and within the eye itself [2]. Distortion due to the geometry of the Scheimpflug camera can be corrected analytically with relative ease, but correction of distortion due to refraction at intermediate ocular surfaces is much less approachable. Measurements of a particular surface are subjected to refraction at all successive surfaces [3,10,21] and traversal through media of individually varying thickness and curvature. Hence, arbitrary quantification errors in one surface are propagated throughout the system [22]. Conversely, magnetic resonance imaging has recently been used for in vivo visualization of structures in the anterior segment, which eliminates the distortion dilemma [21], but suffers from low resolution, signal-to-noise ratio (SNR) constraints, and eye motion artifacts due to longer acquisition times [23]. On the other hand, corneal topography is a rapidly developing technique that provides very detailed and reliable measurements regarding corneal curvature [24][25][26][27], including astigmatism and surface irregularities, although it does not provide information about the remaining ocular surfaces. However, such accurate corneal information could be used to supplement or validate the parameter estimates acquired with our system, or even used as input to inverse optical design to narrow the high-dimensional parameter space.
In this paper, we will discuss some of the obstacles that must be overcome to make the system practical and reliable for clinical use. For instance, we discovered strong coupling between the ocular parameters and a complex likelihood surface. Therefore, the crux of future research will involve the development of an efficient global search algorithm as well as rapid processing techniques.

Data-acquisition
The data consist of the raw detector outputs from a Shack-Hartmann WFS for measurement of aberrations in human eyes, as illustrated in Fig. 1. In the clinical system developed by Straub et al. [28], a narrow collimated beam from a laser diode centered at 780 nm with a 30-nm bandwidth produces a point source on the retina, which acts as a laser beacon and fills the dilated pupil upon reflection. In this preliminary study, we used a single wavelength of 780 nm in our computerized model for simplicity and to reduce computation time. We also illuminated the eye with a Gaussian beam (2 mm at FWHM) at multiple angles of 0, 6, and 12 degrees in the vertical direction to assess both on-axis and off-axis aberrations such as coma and astigmatism. The center of the beam was coincident with the intersection of the optic axis and anterior corneal surface. We treated the aberrated retinal spot as a perfect diffuse scatterer, but did not account for scattering within the ocular media or internal reflections. Furthermore, we excluded relay optics as used in the clinical configuration to conjugate the exit pupil of the eye to the WFS, and simply placed the lenslet array 10 mm away from the corneal apex. Global wavefront tip and tilt were discarded by rotating the WFS for off-axis angles. While we did not incorporate a beamsplitter into the computational model, we will later consider a pellicle beamsplitter to reduce unwanted back-reflections in the clinical system. The lenslet array sampled the aberrated wavefront and produced a grid of focused spots on a detector placed in the focal plane. Generally, the amount of displacement of each spot from its ideal, on-axis location is proportional to the local average slope of the wavefront at the respective lenslet, while finer details of the aberrated wavefront manifest in the spot profiles. Therefore, both the spot positions and profiles provide useful information about aberrations in the eye. Our method does not involve centroid estimation as in conventional wavefront sensing, and the data consist of all pixel intensity values in the focal plane of the WFS, which we refer to as the raw detector outputs. For this simulation, we considered a commercially available micro-WFS, consisting of an array of 600 × 600 μm 2 lenslets with a focal length of 24 mm. The detector pixels had a pitch of 10 μm, and the peak SNR was 10 3 .
We generated an optical design program in MATLAB, which performs non-paraxial raytracing through quadric surfaces, and simulated a trial set of data for ocular parameters comparable to those in the Navarro wide-angle schematic eye model [13], as provided in Table  1. We used the chromatic dispersion model developed by Atchison and Smith [29] to calculate refractive indices at 780 nm. Therefore, we are currently utilizing an effective refractive index for the crystalline lens, but will later incorporate a GRIN distribution according to models suggested in the latest vision science studies [30][31][32]. In addition, we decentered the lens by 0.20 mm in the horizontal direction and -0.10 mm in the vertical direction, which agrees with the experimental range of crystalline lens misalignments [2]. We generated WFS data by tracing a 256 × 256 grid of rays, although about 60% of the rays survived the double-pass after the random retinal reflection and vignetting at the pupil.
Next, we assumed that the system is not photon-starved and used Gaussian statistics to represent electronic noise. For detector elements that are i.i.d. (independent and identically distributed) and noise that is independent of the illumination level, the probability density function (PDF) from which the data are drawn is given by (1) where g is the M × 1 vector set of random data, θ is the set of estimable parameters, and σ 2 is the variance in each detector element [33]. Furthermore, overbars are used to denote means, such that ḡ m (θ) represents the average value of the m th detector element evaluated at θ. Since the noise is zero-mean Gaussian, ḡ m (θ) is simply the output of the optical design program.
We calculated the variance to obtain a peak SNR of 10 3 and applied noise in the data using a Gaussian random number generator.

Maximum-likelihood estimation
After generating the data, we disregarded knowledge of the values of selected parameters and estimated them by maximizing the likelihood. Essentially, ML estimation returns the θ argument which maximizes the probability of occurrence of the observed data, defined as (2) where θ̂ represents an estimate of the vector set of parameters. Since the logarithm increases monotonically with its argument, Eq. (2) can be rewritten as (3) For a purely Gaussian noise model with i.i.d. detector elements, it can be shown that ML estimation according to Eq. (3) reduces to nonlinear least-squares fitting between the data and the output of the optical design program: (4) Immediate advantages are that ML estimation is efficient (i.e., unbiased and yields the best possible variance) if an efficient estimator exists, and that it is asymptotically unbiased as more photons are acquired. However, one challenge is that an accurate probability model must be used that includes all sources of randomness (e.g., photon noise and electronic noise) [33] as well as fluctuations of ocular aberrations associated with live imaging of the eye, including the tear film effect.

Fisher information
The performance of the ML estimator can be analyzed with the Fisher information matrix (FIM), denoted F, because it describes the ability to estimate a vector set of parameters and is related to the sensitivity of the data to changes in those parameters. For a vector parameter of P real components, the FIM is a P × P symmetric matrix with components given by (5) where the angle brackets denote the average. Mathematically, it is the covariance matrix of the gradient of the logarithm of the likelihood. Therefore, off-diagonal entries in the FIM indicate coupling between different pairs of parameters; stronger coupling leads to greater noise in the parameter estimates. For Gaussian i.i.d. data, the FIM components reduce to (6) which depend on the average data vector evaluated at the set of parameters and are inversely proportional to the variance in the detector elements [33].
Inversion of the FIM provides the Cramér-Rao lower bound (CRB) on the variance (i.e. the theoretical minimum possible variance) of the parameter estimates. It is well known that the variance of any unbiased estimate obeys the Cramér-Rao inequality: (7) For an efficient estimator, the inverse of the FIM is the covariance matrix of the parameter estimates; thus, its off-diagonal elements represent coupling between these estimates. We computed the FIM, its inverse, and the CRBs for the chosen subset of parameters using our optical design program.
In general, the FIM can be computed for any system configuration and therefore used to design and optimize the system that acquires the WFS data to be used as input to the inverse optical design, prior to clinical application. We refer to that system as the inverse design system.
Adjustable system parameters to increase Fisher information include the beam size, beam angle, lenslet array geometry, detector element spacing, and variance in the detector elements. Additionally, data from multiple beam angles may serve to decouple pairs of ocular parameters in the FIM, which could improve parameter estimability. We will investigate this more thoroughly in future studies.

Forward model and trial data
The geometrical eye model corresponding to the ocular parameters provided in Table 1 and crystalline lens decentrations stated earlier is illustrated in Fig. 2, which was generated using our optical design program. Sample rays for the multiple beam angles (i.e. 0, 6, and 12 degrees) used in this study are also plotted. The trial set of WFS data used as input to the inverse optical design problem are provided in Fig. 3. The respective focal spots on the retina are shown in Fig. 4 with the coordinate system centered on the optical axis, such that the position on the retina can be read from the axes.

Fisher information and Cramér-Rao bounds
In this proof-of-principle study, we estimated a reduced set of ocular parameters including the posterior radius, thickness, and refractive index of the cornea, the thickness and index of the anterior chamber, and the anterior and posterior radius, thickness, and equivalent index of the crystalline lens. The FIM of the system including data from all beam angles is displayed on a logarithmic scale shown in Fig. 5(a), with values ranging from 10.145 to 17.375. This matrix is order-specific, and the indices of the estimated parameters are listed in Table 1. Note that the jk th entry has units of [units of j th parameter] -1 [units of k th parameter] -1 . The large values in the FIM indicate that the data are very sensitive to changes in the parameters, but that the magnitudes of off-diagonal entries reveal a high degree of undesired coupling between pairs of parameters. This is intuitive since first-order geometrical parameters combine to form various optical quantities, such as refractive index and thickness in optical path length, and curvatures and indices in optical power. An adverse effect of parameter coupling is that an error in one parameter estimate will cause errors in the estimates of all other parameters, while it introduces many local minima in the likelihood surface.
We inverted the FIM as shown in Fig. 5(b), also on a logarithmic scale, and read off the diagonal entries to determine the CRB for each parameter estimate. The corresponding standard deviations (i.e., square-roots) are given in Table 2. These diminutive values permit the estimation of parameters to high precisions even under pessimistic noise levels since we implemented a peak SNR of only 10 3 , and the theoretical lower bound on the variances can immediately be improved by decreasing the variance in the detector pixels. However, the performance of the ML estimator relies greatly on an accurate forward model of the system, particularly when we work with real data. In later studies, we will add complexity to our basic model by incorporating other noise sources, coherence properties of the source as well as the retinal reflection, a GRIN distribution in the crystalline lens according to the latest models in vision science, the Stiles-Crawford effect, scattering within the ocular media, Fresnel reflections between the ocular surfaces, topography of the corneal surface, and the optical effect of the tear film. We also anticipate practical issues such as stray light, misalignments, or other effects that may contribute to estimation errors. Once we include sufficient detail, we can investigate the effect of a model-mismatch by quantifying how much estimation error results from deficiencies in the forward model.

Parameter estimation
Algorithms that implement a likelihood approach to estimating parameters perform a search through parameter space to find a point that maximizes the probability of generating the observed data. Virtually all search algorithms locate an extremum via iterative methods that execute in a variable number of steps, depending on the starting location, the complexity of the probability surface, and values selected for convergence factors. One significant limitation to overcome in our method of parameter estimation is that the estimation step is very time consuming, especially since we are dealing with a complicated surface with many local features and strong coupling between parameters. Making our system practical will require dedicated computer hardware and the development of a proficient global search algorithm tailored to the probability surface. The algorithm may involve prior information obtained from statistical studies or by other reliable modalities such as corneal topography, which can be used to select a starting point or put limits on the parameter values to narrow the search space. An approximate method based on optimization with reverse ray-tracing [34] might give a better starting point for our likelihood approach. However, we conducted the pilot study described in this research by using a custom simplex algorithm for the optimization and a nearby starting point to avoid stagnation in local basins of the surface.
The estimated parameters are listed in Table 2, including the true values used to generate data, the starting point in the ML search, the final estimates, and the fractional error in each estimate. The parameters were estimated to two or three decimal places of accuracy for radii and thicknesses and four to five decimal places for refractive indices. All errors were within a fraction of a percent, indicating very good agreement between the true and estimated values. The reconstructed eye model shows precise overlap, as illustrated in Fig. 6. Although the corneal anterior radius was not estimated assuming that its value is known from direct measurements, we plotted it in the reconstruction to indicate the estimated corneal thickness.
When normalized by the detector variance and total number of elements, the sum-of-squares described in Eqn. 4 is 401.779 at the starting point, 1.65236 at the termination point in the optimization, and 0.998794 at the true global minimum. (As the number of pixels increases, the true minimum should approach unity using the normalization stated above.) We will investigate other ML search algorithms of greater performance in finding the true minimum and achieving the Cramér-Rao bounds.

Conclusion
Patient-specific parameters of the eye were determined using simulated Shack-Hartmann WFS data, a Gaussian noise model, and ML estimation. A system configuration was employed that gave a feasible lower CRB, while an accurate estimate of the parameters was obtained using a nearby starting point. Future analyses will consider other system configurations, noise models, and more complexities in the forward model, such as source coherence, the GRIN distribution in the crystalline lens, and irregularities of the corneal surface. Global optimization algorithms will be investigated using prior information about the eye and characteristics of the likelihood surface to overcome the need for a nearby starting point. Rapid processing techniques will be explored to improve computation time, allowing an increase in model robustness and an exhaustive global search algorithm. This will be accomplished with dedicated computer hardware and parallelization of the optical design program. The ultimate goal is to make the system practical in the clinical setting and to obtain reliable estimates of the full set of ocular parameters for a given patient. Shack-Hartmann WFS system for the eye. In our computational model, we exclude relay optics and place the center of the lenslet array 10 mm from the corneal apex. The angle between the laser beam and the optical axis of the eye is denoted α. Geometrical eye model used to generate WFS data, corresponding to ocular parameters in Table 1. Sample rays for 0, 6, and 12 degrees are shown. The y-axis corresponds to the vertical direction. Corresponding WFS data used as input to inverse optical design, for beam angle: (a) 0 degrees, (b) 6 degrees, and (c) 12 degrees. Focal spots on the retina due to the source beam for angle: (a) 0 degrees, (b) 6 degrees, and (c) 12 degrees.  Reconstructed eye model of the estimated parameters, superimposed with the true values in the data.