Back‐calculation of keratometer index based on OCT data and raytracing – a Monte Carlo simulation

This study aims to develop a raytracing‐based strategy for calculating corneal power from anterior segment optical coherence tomography data and extracting the individual keratometer index, which converts the corneal front surface radius to corneal power.


Background
Corneal power is one of the most important parameters for calculating the power of lens implants (Olsen 1986;Preussner, Wahl & Weitzel 2005), but it cannot be measured in situ.The power of a lens is always defined by the reciprocal of the focal length corrected for the refractive index of the optical media, and the focal length has to be referenced either to the image side principal plane (so-called equivalent power of the lens) or to the front or back vertex (so-called front vertex or back vertex power).Instead of measuring corneal power directly, clinicians measure the shape of the corneal surfaces and estimate corneal power based on a simplified optical model (Olsen & Jeppesen 2018).Most cases are restricted to a measurement of radius of curvature of the corneal front surface (R a ) using, for example a keratometer or a corneal topographer (Olsen & Jeppesen 2018;Langenbucher et al. 2020).Based on any schematic model eye, corneal front surface radius can be converted into corneal power using assumptions such as a fixed ratio of corneal front surface radius to back surface radius (R b ) and a fixed corneal thickness (CCT), and refractive indexes of cornea and aqueous humour.It is obvious that in situations where the proportions do not properly match the underlying model such a conversion results in incorrect values for corneal power, whichinserted in a lens power calculation schemeyields inappropriate results for the power of the lens and a refractive surprise after cataract surgery (Olsen 1986;Langenbucher, Haigis & Seitz 2004;Preussner, Wahl & Weitzel 2005;Haigis 2012).
With modern tomography techniques, the shape of the cornea can be assessed with a high precision (Ho et al. 2008;Olsen & Jeppesen 2018).In particular, Scheimpflug tomographers with or without integrated Placido topography and high-resolution optical coherence tomographers (OCT) permit measurements of corneal front and back surface curvature as well as asphericity and corneal thickness with a high degree of accuracy and reliability.The refractive indices of the cornea and aqueous humour cannot be directly measured in situ, but for both parameters the variation seems to be quite low, and most schematic model eyes describe congruent data with n C = 1.376 for the cornea and n AQ = 1.336 for the aqueous humour (Liou & Brennan 1997).
In clinical routine, conversion of corneal front surface radius R a to corneal power CP is performed with a keratometer index n K , which is an assumed parameter or calibration value (Olsen 1986;Fam & Kim 2007;Langenbucher et al. 2020).In many keratometers or topographers, the conversion is hard-coded meaning that the user cannot modify the keratometer index (Langenbucher, Haigis & Seitz 2004;Haigis 2012).The formula behind reads.
Several different values are used, including the so-called Javal index with n K = 1.3375 and the Zeiss index with n K = 1.332 which are the most popular calibrations (Olsen 1986;Fam & Lim 2007;Olsen & Jeppesen 2018).Interpreting these two keratometer indices based on the Gullstrand schematic model eye, the Javal index refers to the back vertex power and the Zeiss index to the front vertex power of the cornea.This results in a difference between R a converted to CP using the Javal, and the Zeiss index is around 0.75 dioptres on average (Langenbucher et al. 2020).This means that these two calibrations cannot be used interchangeably, and if a lens power calculation scheme is designed to use CP instead of R a the user should check carefully which calibration should be used.
Even if the radius of curvature of both corneal surfaces together with central corneal thickness is measured, and a thick lens formula used for calculation of corneal power, the calculation is mostly performed with restrictions of linear Gaussian optics (in the paraxial space), and the asphericity of the front (Q a ) and back surface (Q b ) and the aperture of the eye (PUP) are not taken into account.
The purpose of this study was to assess corneal power and back-calculated keratometer index based on a large dataset of high-resolution anterior segment OCT using raytracing techniques in terms of a Monte Carlo simulation, and to present the results for back-calculated keratometer index as a function of corneal front and back surface radius, asphericity of corneal front and back surface, corneal thickness and pupil size.

Measurement data
Out of a total of 11,276 measurements in patients without pathologies or a history of ocular surgery performed between January 2019 and July 2020, 10,218 OCT measurements with a complete set of data (10,218 eyes of 8,430 patients) were enrolled in this study.In all of these measurements, the internal quality check of the tomographer indicated a proper measurement and no movement artefacts.All measurements were carried out at Augenklinik Castrop-Rauxel.An ethics approval was not required for this study.This study is a retrospective evaluation of data which were collected during routine examinations.No extra examinations or measurements were performed.
The data were transferred to us in an anonymized fashion, which precludes back-tracing of the patient.
Data from the Casia 2 (Tomey, Nagoya, Japan) were exported in standard.csvdata format and imported to MATLAB (MathWorks, Natick, USA, Version 2019b) for subsequent data analysis.From the dataset, we used the central corneal front surface radius R a , the central corneal back surface radius R b , corneal eccentricity of the corneal front (ECC a ) and back surface (ECC b ) both derived in the central 6 mm zone, central corneal thickness CCT and projected (visible) pupil size.Corneal eccentricity with a positive value in the dataset indicating a prolate shape of the corneal surface was converted to corneal asphericity using.
Corneal surfaces with an oblate shape indicated in the dataset with a negative value of eccentricity were converted to We assumed a centred optical system without any tilt of the 2 corneal surfaces.Both surfaces were considered as quadric surfaces described by a central radius of curvature (R a or R b ) and an asphericity (Q a and Q b ).The axial symmetry of this model means that restriction to a 2-dimensional raytracing strategy was sufficient (Langenbucher et al. 2006;Langenbucher et al. 2011;Langenbucher et al. 2014;Langenbucher et al. 2016).The apex of the corneal front surface was assumed to be located at z = 0 and the corneal back surface apex at a distance of z = CCT.
A collimated bundle of 601 rays starting from a plano surface at z = 0 (apex plane of the cornea) was projected to each cornea where the radial ray spacing was adapted in a quadratic fashion from centre to periphery to realize an equally spaced sampling over the pupil size (area-correction).The diameter of the ray bundle was adjusted to the measured diameter of the pupil PUP.Then, the ray-surface intersection was calculated and the direction of the refracted rays was derived using Snell's law for the corneal front surface (Langenbucher et al. 2006;Langenbucher et al. 2011;Langenbucher et al. 2014;Langenbucher et al. 2016).After tracing all rays through the cornea, ray intersection with the corneal back surface was calculated and the direction of the refracted rays in the anterior chamber was derived using the Snell's law.The best focus plane z = z F was determined based on the criterion of least scatter (root mean squared error).A nonlinear search algorithm (Levenberg--Marquardt algorithm (Levenberg 1944;Marquardt 1963); maximum iterations: 100; termination tolerance on the function value: 1e-16; termination tolerance for the step size; 1e-14) was used for calculating the best focus (Langenbucher et al. 2006).
At best focus plane z = z F , the cumulative optical path length was calculated for each ray starting from object plane z = 0, and the standard deviation of optical path length for all rays in the ray bundle was used as a measure for optical aberration (Langenbucher et al. 2006).From the best focus plane z F , we assessed corneal power referenced to the corneal front vertex plane (as this surface is the reference for biometry prior to cataract surgery) and backcalculated the individual keratometer index from z F and R a .
Data from the N = 10,218 eyes were analysed descriptively using mean, standard deviation (SD), median, minimum and maximum.We then investigated the data pool using different multiple linear regression models (Model 1: Multivariate linear model using all 6 input parameters R a , R b , Q a , Q b , CCT and PUP; Model 2: Simplified linear model with three parameters R a , R b and CCT; Model 3: Simple linear model based on a single parameter (corneal front surface radius model R a ) to investigate the overall performance of the models as well as the regression coefficients, standard errors and the relevance (P-value) of the effect sizes.Finally, we plotted the back-calculated keratometer index as a function of R a , R b , Q a , Q b , CCT and PUP.

Results
The descriptive statistics of the 6 input data of R a , R b , Q a , Q b , CCT and PUP is shown in Table 1.
In Fig. 1, we have plotted the histogram of the distribution for all N = 10,218 input data.
Table 2 displays the descriptive statistics for the output variables: position of the best focus z F , corneal power CP, root mean squared error of optical path length differences and back-calculated keratometer index.
Figure 2 provides a histogram of the distribution of the best focus plane z F , the corneal power CP referenced to the front vertex plane of the cornea and the back-calculated keratometer index n K .
The results of the linear model for estimation of the keratometer index from all 6 input parameters R a , R b , Q a , Q b , CCT and PUP yield: The root mean squared error/R 2 / significance vs. a constant model is 5.54e-4/8.98e-1/<0.0001.The standard errors of the regression coefficients are 1.76e-4 for the intercept and 3.32e-5/ 2.37e-5/3.43e-5/2.73e-5/1.52e-7/7.23e-6for R a /Q a /R b /Q b /CCT/PUP.The significance level of all 6 effect sizes and the intercept is lower than 6.44e-15.
The simplified linear model for estimation of the keratometer index from the input data used for a paraxial backcalculation of the corneal power as a thick lens R a , R b , and CCT yields: The root mean squared error/R 2 / significance vs. a constant model is 1.69e-3/5.50e-2/<0.0001.This model is not appropriate for estimation of the keratometer index.The standard error of the regression coefficients are 5.30e-4 for the intercept and 9.91e-5/1.02e-4/4.61e-4 for R a /R b /CCT.The significance level of the 3 effect sizes and the intercept is lower than 1e-16.
If calculating the keratometer index with a linear model based only on the corneal front surface radius model R a , the model description reads as follows: The root mean squared error/R 2 / significance vs. a constant model is 1.66e-3/4.00e-3/0.229.This model is not appropriate for estimation of the keratometer index.The standard errors of the regression coefficients are 6.15e-4 for the intercept and 7.96e-5 for R a .The significance level of the R a as effect sizes and the intercept is lower than 1e-16.
Figure 3 shows the back-calculated keratometer index as a function of corneal front surface radius R a (Fig. 3 A), corneal front surface asphericity Q a (Fig. 3B), corneal back surface radius R b (Fig. 3C), corneal back surface asphericity Q b (Fig. 3D), central corneal thickness CCT (Fig. 3E) and pupil size PUP (Fig. 3F).If we restrict to a linear estimation model with one effect size, only the asphericity of the corneal front surface Q a (Fig. 3B) appears to be a good predictor for the keratometer index.In all graphs, we have included the linear fit in terms of minimizing the squared fit error.The respective slope of the fit is provided in the legend of each graph.Due to the special characteristics of the data distribution in Fig. 3B, we added a sigmoidal function which is defined by the lower/upper boundary of 1.3260/1.3345,an inflection point at Q a = −0.4895and a slope at the inflection point of 3.9537.The linear and the sigmoidal model both yielded a comparable fitting performance with a root mean squared fit error of 1.38e-3 and 1.36e-3, respectively.

Discussion
There is always discussion in the community of cataract surgeons about the proper keratometer index (Olsen 1986;Preussner, Wahl & Weitzel 2005;Fam & Lim 2007;Olsen & Jeppesen 2018).If we keep in mind that the basic measure is the shape of the corneal front surface, corneal power is always a result of a conversion from corneal front surface radius alone (e.g.keratometer), of corneal front and back surface radius together with central corneal thickness [e.g.thick lens formula with linear Gaussian optics Table 1.Descriptive statistics of the input data with mean, standard deviation (SD), median, minimum and maximum.R a refers to the corneal front surface radius, Q a to the corneal front surface asphericity, R b to the corneal back surface radius, Q b to the corneal back surface asphericity, CCT to the central corneal thickness and PUP to the pupil size considered at the corneal plane  (Liou & Brennan 1997), and if we do not use a full set of tomographic data we must make some assumptions which are typically derived from the schematic model eye.As an example, the Zeiss and the Javal index are both based on the classical Gullstrand model eye, and the derivation of both keratometer indices involves the use of the ratio between corneal front to back surface curvature and corneal thickness (Langenbucher, Haigis & Seitz 2004;Haigis 2012).The difference between the Zeiss and Javal index is primarily that the two indices refer to a different reference plane.
With modern optical measurement techniques, we get accurate and reproducible measures of the entire corneal shape including front surface and back surface architecture and central thickness.Therefore, it is obvious to use these measurement data to upgrade determination of corneal power and to question the traditional keratometer index.We extracted a large dataset of anterior segment OCT data from normal eyes measured in the last 18 months at one clinical centre in Germany.Based on this dataset and an estimate of the refractive index for the cornea and aqueous humour, we used raytracing techniques to calculate the best focus plane.As we used a concentric optical model based on the assumptions that both corneal surfaces are aligned and perpendicular to the fixation axis, we could restrict the analysis to a 2D raytracing strategy instead of 3D raytracing which is more complex and time-consuming.If we were to trace an equidistant bundle of rays (in one meridian) through the cornea, the data would be much denser in the centre and sparser in the periphery.Therefore, we decided to adapt the ray density in a quadratic fashion to resample a homogeneous distribution over the entire pupil.The further data processing was performed in a similar way to a classical Monte Carlo simulation.The back-calculated keratometer index for all eyes was analysed as a function of all potential effect sizes (R a , R b , Q a , Q b , CCT and PUP) separately, and after that multivariate linear models were defined which simplify the raytracing calculation and make a linear prediction for the appropriate keratometer index to be used for translating corneal front surface curvature to a corneal power.In total, 3 multivariate linear models were described, one of them considering all 6 potential effect sizes, one considering 3 effect sizes (R a , R b and CCT) and one as a simple linear model with one effect size (R a ) only).
For the model that we present in this paper, we used a collimated bundle of rays starting from the corneal front vertex plane and traced through both corneal surfaces to the best focus plane.Alternatively, a divergent bundle of rays starting, for example on-axis at a refractometry lane distance of 6 m could be used (Langenbucher et al. 2020).By accumulation of the optical path length from object to image (through air, cornea and aqueous humour), we were able to directly read out the optical path length differences for all the rays, which refer to the optical aberration (in µm).If we reference the best focal plane to the corneal front or back vertex plane, we can Table 2. Descriptive statistics of the output data with mean, standard deviation (SD), median, minimum and maximum.z F refers to the position of the best focus in terms of lowest ray scatter, CP to the corneal power with respect to the corneal front vertex plane, SD OPL to the standard deviation of optical path length differences from the object plane z = 0 to the best focus plane z = z F and n K to the back-calculated keratometer index easily derive the power of the cornea with respect to the front surface or back surface (Langenbucher et al. 2020).As we require corneal power with respect to the corneal front vertex plane for intraocular lens power calculation, we decided to use the term corneal power with the reciprocal of the distance from z = 0 to z = z F corrected for the refractive index of the aqueous humour.
In this paper based on a large dataset, we derived a focal length of around 31.2 mm behind the corneal front vertex, which refers to a corneal power of around 42.9 dioptres.On average, lateral ray scatter (root mean squared error) at the best focal plane was 3.32 µm, standard deviation of optical path length differences was 0.18 µm and the back-calculated keratometer index was 1.3317 on average.Considering that we used raytracing techniques and adapted the diameter of the ray bundle traced through the eye individually to the measured pupil size, the resulting keratometer index is very similar to that what was shown in the literature before (Fam & Lim 2007;Olsen & Jeppesen 2018;Langenbucher et al. 2020).In a recent study, we used paraxial calculation strategies based on a thick lens model and back-calculated the refractive index for clinical data after cataract surgery for a finite lane distance for refractometry and obtained a keratometer index of around 1.330 (Langenbucher et al. 2020).As the cornea shows some positive spherical aberration (in this study around 0.18 µm) which is not fully compensated by the typical prolate shape, it is obvious that the keratometer index here is somehow slightly higher compared to a paraxial assessment.
It was surprising to us that in the Monte Carlo simulation and the monovariate assessment of the keratometer index as shown in Figure 3, only the front surface asphericity, seems to affect the keratometer index.For instance, the mono-variate model (as shown in model 3 for the corneal front surface radius) does not show any significant difference from a constant model.Even using corneal front and back surface radius and central corneal thickness, the regression model does not show a significant difference compared to the constant model (model 2).But if we use a multivariate regression analysis using all input parameters (R a , Q a , R b , Q b , CCT and PUP; model 1) the model shows a very good performance and therefore it seems that a proper estimate of the keratometer index is possible with this model.
In order not to overload this simulation, we restricted the analysis to a rotationally symmetric optical model (without decentration and tilt).Using that model, we characterized our dataset with 6 input parameters.In general, such a Monte Carlo simulation would be much more difficult if astigmatic surfaces or even free-form surfaces and misalignments such as decentration and tilt were included: firstly this would require 3D raytracing, and secondly we would require a general description of both surfaces instead of quadric surfaces which are more easily handled.Finally, it would no longer be possible to describe the results with a linear multivariate model (e.g.model 1 in this paper).Therefore, interpretation might be very complex.
Clinicians have to keep in mind that we could not ignore the variation of the keratometer index in general.Based on our dataset with N = 10,218 data points, we observed for n K an individual variation of 1.3233 to 1.3390 and a standard deviation of 0.0017.This means, for example for an average corneal front surface curvature of 7.7 mm that considering the 2.5% and 97.5% quantile of the back-calculated keratometer index (1.3242 and 1.3371), in the worst-case scenario corneal power could be 324.2/7.7 = 42.17dioptres or 337.1/7.7 = 43.78dioptres (difference 1.67 dioptres).In reality, these extreme values are retrieved from extreme combinations of input parameters, and in normal cases the variation will be much smaller!A simple alternative to this raytracing-based calculation would be to use a straightforward paraxial model for thick lenses as shown in (Langenbucher et al. 2020), which considers at least the radius of curvature of both corneal surfaces and the central corneal thickness and could be applied with data from any modern optical biometer which measures the radii of both corneal surfaces (e.g.IOLMaster 700, Carl-Zeiss-Meditec, Jena, Germany).
In conclusion, we have shown with a raytracing-based calculation strategy applied to a large dataset of N = 10,218 measurements, how the back-calculated keratometer index varies with combinations of corneal front and back surface radius, corneal front and back surface asphericity, central corneal thickness and pupil size.Based on the best focus plane, we determined the corneal power and derived an individual keratometer index which translates the corneal front surface radius to the corneal power.We fitted a multiple linear regression model to this dataset, which shows a very good performance and which could be used to estimate the individual keratometer In addition to the data, the linear fit in terms of minimisation of the least squared fit error is displayed.In Figure 3B due to the distribution of the data a sigmoidal function was fitted.If we restrict to 1 effect size, only the asphericity of corneal front surface is usable for prediction of keratometer index.
index providing all input parameters are available.

Fig. 1 .
Fig. 1.Histogram of the distribution of input parameters (N = 10,218).R a refers to the corneal front surface radius, Q a to the corneal front surface asphericity, R b to the corneal back surface radius, Q b to the corneal back surface asphericity, CCT to the central corneal thickness and PUP to the pupil size.

Fig. 2 .
Fig.2.Histogram of the distribution of output parameters (N = 10,218).z = z F refers to the best focus plane (plane with the lowest ray scatter), CP to the corneal power with respect to the corneal front vertex plane z = 0 derived from z F and n K to the back-calculated keratometer index.

Fig. 3 .
Fig.3.Back-calculated keratometer index as a function of corneal front surface radius R a (A), corneal front surface asphericity Q a (B), corneal back surface radius R b (C), corneal back surface asphericity Q b (D), central corneal thickness CCT (E) and pupil size PUP (F).In addition to the data, the linear fit in terms of minimisation of the least squared fit error is displayed.In Figure3Bdue to the distribution of the data a sigmoidal function was fitted.If we restrict to 1 effect size, only the asphericity of corneal front surface is usable for prediction of keratometer index.