Prediction model for best focus, power, and spherical aberration of the cornea: Raytracing on a large dataset of OCT data

Purpose To analyse corneal power based on a large optical coherence tomography dataset using raytracing, and to evaluate corneal power with respect to the corneal front apex plane for different definitions of best focus. Methods A large OCT dataset (10,218 eyes of 8,430 patients) from the Casia 2 (Tomey, Japan) was post-processed in MATLAB (MathWorks, USA). Using radius of curvature, corneal front and back surface asphericity, central corneal thickness, and pupil size (aperture) a bundle of rays was traced through the cornea. Various best focus definitions were tested: a) minimum wavefront error, b) root mean squared ray scatter, c) mean absolute ray scatter, and d) total spot diameter. All 4 target optimisation criteria were tested with each best focus plane. With the best-fit keratometer index the difference of corneal power and keratometric power was evaluated using a multivariate linear model. Results The mean corneal powers for a/b/c/d were 43.02±1.61/42.92±1.58/42.91±1.58/42.94±1.59 dpt respectively. The root mean squared deviations of corneal power from keratometric power (nK = 1.3317/1.3309/1.3308/1.3311 for a/b/c/d) were 0.308/0.185/0.171/0.209 dpt. With the multivariate linear model the respective RMS error was reduced to 0.110/0.052/0.043/0.065 dpt (R² = 0.872/0.921/0.935/0.904). Conclusions Raytracing improves on linear Gaussian optics by considering the asphericity of both refracting surfaces and using Snell’s law of refraction in preference to paraxial simplifications. However, there is no unique definition of best focus, and therefore the calculated corneal power varies depending on the definition of best focus. The multivariate linear model enabled more precise estimation of corneal power compared to the simple keratometer equation.

However, there is no unique definition of best focus, and therefore the calculated corneal power varies depending on the definition of best focus. The multivariate linear model enabled more precise estimation of corneal power compared to the simple keratometer equation.

Background
In most clinical situations dioptric corneal power is not really decisive and does not have any diagnostic or therapeutic consequence [1]. However, corneal power does have a high impact in refractive surgery and cataract surgery, and clinicians rely on accurate values for procedures such as the determination of an appropriate replacement lens power during cataract surgery, implantation of an additional lens into a phakic or pseudophakic eye, or determination of an ablation profile for corneal refractive surgery. In normal eyes an overestimation / underestimation of corneal power by 1 dioptre significantly affects the required power of a replacement lens during cataract surgery, resulting in an underestimation / overestimation of around 1.5 dioptres [2].
Currently, corneal power is most commonly derived from curvature measurement of the corneal front surface using a manual or automatic keratometer or topographer. The radius of curvature is measured in the mid-periphery at a distance of around 1.5 to 2.0 mm from the corneal apex, and radius of curvature is converted to corneal power using a keratometer index [3]. In the last decade more and more tomographers have been launched to the market with the ability to measure height data of the corneal front and back surface as well as corneal thickness (CCT). But in most cases derivation of corneal power has not changed as a result of these new measurement options [1].
Both corneal surfaces can be approximated to good accuracy by quadric surfaces. For stigmatic surfaces, a 2-axis ellipsoid characterized by a central radius of curvature and asphericity is sufficient [4]. In the case of astigmatism, a fit with a 3-axis ellipsoid is used, characterised by: 2 radii of curvature (in both cardinal meridians), an asphericity, and the orientation of the ellipsoid [5]. In the most general case, a biconic surface fit could be used to extract the radii of curvature and asphericity in both cardinal meridians together with the orientation of the steep or flat meridian [6]. The fewer degrees of freedom used to describe the parametric surface, the more robust the parameters extracted from the fit.
Converting corneal front surface radius of curvature (R f ) to corneal power using a keratometer index involves a simplified optical model of the cornea, in which the refractive index of cornea and aqueous humour, central corneal thickness, and the ratio of corneal front to back surface curvature (back surface radius R b ) from a schematic model eye are used. Using these assumptions, mostly taken from the Gullstrand model eye, in the normal eye corneal power referenced to the front (keratometer index n K = 1.332) or to the back surface (n K = 1.3375) may be derived [7].
In reality, the ratio of corneal front to back surface curvature and corneal thickness varies individually from eye to eye, and more modern schematic model eyes incorporate a more curved back surface curved than that of the Gullstrand eye from 1909. Additionally, this conversion from corneal front surface curvature to power is based on paraxial optics (linear Gaussian optics) without properly considering the asphericity of either surface [2].
Raytracing is a standard strategy for analysing optical systems [8]. A bundle of rays is traced through the refractive surfaces using Snell's law (instead of the paraxial simplification) and restricted by an aperture stop. This strategy considers the curvature and asphericity of both corneal surfaces (front and back surface asphericity Q f and Q b ) in addition to the individual aperture diameter (PUP, pupil size) [9]. Tracing rays from a point light source at a finite (e.g. 4-6 m in front of the eye) or infinite distance, we evaluate the best focus position. But there are complementary definitions of best focus: e.g. the plane with the lowest wavefront (WF) error, the plane with the lowest root mean squared (RMS) scatter, the lowest mean absolute (MA) scatter of rays, or the smallest total spot diameter (TSD).
The objective of this paper is to • use 2D raytracing based on a rotationally symmetric optical model of the cornea defined with central curvature (R f and R b ) and asphericity (Q f and Q b ) of the corneal front and back surface and to evaluate the best focus, • derive best focus in terms of minimising WF error, RMS scatter, MA scatter, and TSD, • describe corneal spherical aberration (SA), mean WF error, RMS scatter, MA scatter, and TSD for all 4 focal planes, • calculate corneal power from the best focus position with respect to the corneal front surface as required for lens power calculation for cataract surgery.
• provide multivariate linear models for corneal power as a function of R f , Q f , R b , Q b , CCT and PUP in a large dataset of > 10,000 measurements from an OCT device.

Measurement data
This retrospective study included 10,218 eyes of 8430 patients examined between January 2019 and July 2020 at Augenklinik Castrop-Rauxel. All data are from a cataract population, obtained during pre-cataract examination. The 10,218 datasets from the Casia 2 (Tomey, Nagoya, Japan) were batch exported in standard .csv format and imported to MATLAB (MathWorks, Natick, USA, Version 2019b) for subsequent data analysis. Incomplete datasets were filtered out prior to exporting the data. From this data set, we used the central corneal front surface radius R f , the central corneal back surface radius R b , corneal eccentricity of the corneal front and back surface both derived in the central 6 mm zone, central corneal thickness CCT, and projected (visible) pupil size. Corneal eccentricity was converted to corneal asphericity, and quadric surfaces were defined for the corneal front and back surface based on central radius of curvature and asphericity. The sag z of these quadric surfaces is defined by: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Where z 0 refers to the axial position of the surface apex, R f,b to the central radius of curvature of the front and back surface, and Q f,b to the asphericity of the front and back surfaces respectively.

Raytracing scheme
For the cornea we assumed a centred optical system without tilt. Both surfaces were considered as quadric surfaces as described above, and the apex of the corneal front surface was considered to be located at the origin, with the apex of the back surface at z = CCT. The axial symmetry of this model means that a 2-dimensional raytracing strategy was sufficient [10,11]. The respective refractive indices of cornea and aqueous humour were taken from the Liou & Brennan schematic model eye [12].
A collimated bundle of 1001 rays starting from a plano surface at z = 0 was projected to the corneal front surface. Rays were equally spaced in the radial direction over the projection of the individual pupil size (PUP) to the corneal front surface from-PUP/2 to PUP/2. To ensure an equidistant spacing over the pupil (area correction), rays were weighted in a quadratic manner with respect to the distance from the optical axis. For each ray we calculated the intersection with the corneal front surface, the direction of the refracted ray (within the cornea) using Snell's law, and the intersection with the corneal back surface [8]. Again, applying Snell's law we derived the direction of the refracted ray in the anterior chamber filled with aqueous humour.
Behind the corneal back surface we calculated the wavefront curvature and fitted a circle to obtain a preset value for the focus (centre of the circle). Starting from this centre we applied a nonlinear annealing algorithm to extract the best focus planes with: • the minimal RMS WF error (BFP-WF), • the minimum RMS ray scatter (BFP-RMS), • the minimum MA ray scatter (BFP-MA), • the minimum total spot diameter TSD (BFP-TSD).
For all 4 best focus planes we analysed the RMS wavefront error (WFE), the RMS ray scatter (RMSS), the MA ray scatter (MAS), plus the TSD, together with the Zernike coefficient for spherical aberration Z40 referenced to a central 6 mm zone at the corneal front surface plane [13][14][15].
From the location of the best focus planes we calculated corneal power with respect to the corneal front apex plane. Defining a best fit keratometer index for each best focus plane, we investigated the dataset of 10,218 OCT measurements to define multivariate linear models describing corneal power as functions of keratometric power (based on the best fit keratometer index) and linear correction terms for corneal front and back surface radius R f and R b , corneal asphericity Q f and Q b , central corneal thickness CCT and pupil size PUP.

Results
The descriptive statistics of the 6 input data of R f , R b , Q f , Q b , CCT, and PUP are shown in Table 1 in terms of mean, standard deviation, median, and 90% and 99% confidence level. Table 2 displays the descriptive statistics for the position of the best focus planes BFP-WF, BFP-RMS, BFP-MA, and BFP-TSD alongside with the refractive power referenced to the corneal front apex plane. Fig 1 shows the RMS wavefront error WFE (row 1), the RMS ray scatter RMSS (row 2), the mean absolute ray scatter MAS (row 3), as well as the total spot diameter TSD (row 4) for each of the 4 best focus planes BFP-WF, BFP-RMS, BFP-MA, and BFP-TSD.
On the diagonal of the plot we find those conditions where the target criterion fits to the respective best focal plane (e.g. BFP-WF and WFE, row 1 / column 1 or BFP-RMS and RMSS, row2 / column 2). Fig 2 displays the ratio of RMS wavefront error to RMS wavefront error at BFP-WF, ratio of RMS ray scatter to RMS ray scatter at BFP-RMS, ratio of MA ray scatter to MA ray scatter at BFP-MA, and ratio of total spot diameter TSD to TSD at BFP-TSD. Table 3 provides the descriptive statistics of the corneal spherical aberration in terms of Zernike coefficient Z40 extracted from the wavefront error at corneal front apex plane for a central optical zone of 6 mm, plus the keratometer index derived from corneal power and corneal front surface radius R f . for BFP-WF, BFP-RMS, BFP-MA, and BFP-TSD.
The customised keratometer index for BFP-WF, BFP-RMS, BFP-MA, and BFP-TSD calculated with a least squares fit yields: n K = 1.3317 forBFP-WF, n K = 1.3309 for BFP-RMS, n K = 1.3308 for BFP-MA, n K = 1.3311 for BFP-TSD.  Table 4 shows the respective descriptive data.
Modelling corneal power with respect to best focus planes 1 to 4 results in the following equations: Model 1: For BFP-WF corresponding to optimisation for root mean squared wavefront error WFE: The Description is given in terms of mean, standard deviation (SD), median, minimum and maximum, 0.5%, 5%, 95% and 99.5% quantiles (90% and 99% confidence intervals). Best focal planes BFP-WF, BFP-RMS, BFP-MA, and BFP-TSD refer to the focal plane with the lowest RMS wavefront error, the lowest RMS ray scatter, the lowest mean absolute scatter, and the lowest total spot diameter respectively.

Discussion
In many ophthalmological situations the dioptric corneal power is of minor relevance, without direct consequence to diagnosis or therapy [1]. However, in some situations such as planning a cataract surgery or refractive intervention at the cornea or lens, the corneal power has an  The description is given in terms of mean, standard deviation (SD), median, minimum and maximum, 0.5%, 5%, 95% and 99.5% quantiles (90% and 99% confidence intervals).
https://doi.org/10.1371/journal.pone.0247048.t003 extremely high impact. In modern biometry for cataract surgery, misinterpretation of corneal power is one of the most important error sources [1,2,8]. In spite of this, in most clinical cases  Corneal power was derived from the best focus plane with respect to the front apex plane, and keratometer power was calculated using (n K -1)/R f with a keratometer index specified in Table 3. The 4 best focus planes indicating lowest RMS WF error (BFP-WF), lowest RMS spot size (BFP-RMS), lowest mean absolute scatter (BFP-MA), and lowest overall spot diameter (BFP-TSD). The description is given in terms of mean, standard deviation (SD), median, minimum and maximum, 0.5%, 5%, 95% and 99.5% quantiles (90% and 99% confidence intervals). https://doi.org/10.1371/journal.pone.0247048.t004 ophthalmologists tend to restrict themselves to measurement of corneal front surface radii and interpret corneal power with a simplified thin lens model using a keratometer index for conversion. However, even if the shape of the corneal front and back surface is measured with a tomographer and corneal thickness and the diameter of the aperture stop is known, raytracing does not solve all problems: even using Snell's law, it is necessary to specify the best focus plane of the cornea. There are several options in optics to define the best focus, and in professional raytracing tools we have for instance, the options of minimising the root mean squared wavefront error, the root mean squared ray scatter, the mean absolute ray scatter, or the overall spot size at the focus plane. Besides these options, we could also optimise for the Strehl ratio or some characteristics of the modulation transfer function.
In this paper we wish to show that there is no general rule for defining best focus, and depending on the selection of a target criterion such as optimising for the RMS wavefront error WFE or RMS ray scatter RMSS we obtain different values for the best focus plane. And with different locations of the best focus plane we obtain different values for the corneal power [7,13], for the keratometer index [3], or for spherical aberration [13][14][15]. We restricted the study in this paper to 4 different options for definition of the best focus plane: minimisation of the RMS wavefront error (BFP-WF), the RMS ray scatter (BFP-RMS, spot size), the mean absolute ray scatter (BFP-MA, mean spot size), and the total spot size (BFP-TSD, overall diameter of the spot). Fig 1 shows the behaviour of these measures for all 4 best focus planes. Clearly, if we optimise e.g. BFP-WF, the RMS wavefront error WFE at BFP-WF will be lower compared to BFP-RMS, BFP-MA, or BFP-TSD. Therefore, the lowest measures are always found on the diagonal in Fig 1. By referencing to those optimised conditions (e.g. calculating the ratio of RMS wavefront error at best focus planes 2 to 4 to RMS wavefront error at best focus plane 1) we get a measure of robustness as shown in Fig 2. Based on our dataset these ratios range between 1 and 5, where 1 refers to a situation with no loss in our target criterion if we optimise for a condition A and evaluate condition B. For instance, if we consider the ratio of RMS wavefront error at BFP-RMS, BFP-MA, or BFP-TSD to RMS wavefront error at BFP-WF we see that the RMS wavefront error increases by around 100-200% (the ratio ranges between 2 and 3). It is interesting to note that if we optimise for RMS ray scatter RMSS and analyse BFP-MA or optimise for MA ray scatter MAS and analyse for BFP-RMS the ratios are more or less 1, which means that based on our dataset, both optimisation criteria could be used interchangeably.
Another interesting finding is that even if we customise the keratometer index to represent corneal power with respect to BFP-WF, BFP-RMS, BFP-MA, and BFP-TSD, we find some deviation of corneal power from keratometric power in the range -0.56 to 0.44 dpt for BFP-WF, -0.31 to 0.27 dpt for BFP-RMS, -0.28 to 0.24 dpt for BFP-MA, and -0.36 to 0.30 dpt for BFP-TSD for the 95% confidence interval ( Table 4). This means that considering the 95% confidence interval, even in the simplified case with centred quadric surfaces, the corneal power derived with raytracing techniques deviates from the keratometric power by typically ±¼ to ±½ dioptre. This deviation has to be taken into account for lens power calculation prior to cataract surgery.
We attempted to characterise the deviation of corneal power from keratometric power (as the dependent target variable) with a multivariate linear model. As effect sizes we selected the central curvature of the corneal front and back surface, the asphericity of the corneal front and back surface, central corneal thickness, and diameter of the pupil. The root mean squared fit errors of models 1 to 4 are 0.110 / 0.052 / 0.043 / 0.065 mm respectively. Comparing these root mean squared fit error of models 1 to 4 to the root mean squared error of the deviation of corneal power minus keratometric power (0.308 / 0.185 / 0.171 / 0.209 mm), we see that with models 1 to 4 the deviation is significantly reduced. The performance of model 1 seems to be slightly lower (R 2 = 0.872) compared to models 2 to 4 (R 2 = 0.921 / 0.935 / 0.904). But overall, the deviation is reduced with the multivariate linear models to around 1/3.
In conclusion, corneal power can be derived from corneal front and back surface shape including central radius of curvature and asphericity, central corneal thickness, and pupil diameter. If we restrict our analysis to a centred optical system without tilt, the aspherical corneal front and back surface can easily be represented with quadric surfaces and raytracing simplified from 3D to 2D. We should be aware that, when using raytracing techniques to evaluate the focusing properties of the cornea to overcome the limitations of paraxial optics, there is no generally accepted definition regarding the best focus. Depending on the best focus definition we obtain different best focus planes and therefore different corneal powers. Keratometric power does not fully reflect the refractive properties of the cornea, and the deviation between corneal power and keratometric power ranges between ±¼ and ±½ dpt depending on the selection of the best focus criterion.