Adaptive Cornea Modeling from Keratometric Data

Abstract

Purpose.: To introduce an iterative, multiscale procedure that allows for better reconstruction of the shape of the anterior surface of the cornea from altimetric data collected by a corneal topographer.

Methods.: The report describes, first, an adaptive, multiscale mathematical algorithm for the parsimonious fit of the corneal surface data that adapts the number of functions used in the reconstruction to the conditions of each cornea. The method also implements a dynamic selection of the parameters and the management of noise. Then, several numerical experiments are performed, comparing it with the results obtained by the standard Zernike-based procedure.

Results.: The numerical experiments showed that the algorithm exhibits steady exponential error decay, independent of the level of aberration of the cornea. The complexity of each anisotropic Gaussian-basis function in the functional representation is the same, but the parameters vary to fit the current scale. This scale is determined only by the residual errors and not by the number of the iteration. Finally, the position and clustering of the centers, as well as the size of the shape parameters, provides additional spatial information about the regions of higher irregularity.

Conclusions.: The methodology can be used for the real-time reconstruction of both altimetric data and corneal power maps from the data collected by keratoscopes, such as the Placido ring–based topographers, that will be decisive in early detection of corneal diseases such as keratoconus.

There is an increasing need for reliable and precise modeling of corneal surfaces, motivated both by technological and clinical applications. Given the significance of the shape of the front surface of the cornea to the refraction of the eye ¹ and the ability to correct refractive errors by laser ablation of the front surface of the cornea, a detailed wavefront error analysis of individual corneal topography data is crucial for clinicians as a basis for customizing treatment. It has been recognized that the corneal front surface generally provides the bulk of the ocular aberrations in the postsurgical or pathologic eye. ²  

Corneal modeling plays an essential role in diagnosing and managing keratoconus and other corneal diseases, with the goal of assessing the suitability of a subject for the treatment and to prevent improper refractive surgeries. ³ Also, the great role of the reliable visualization tools in clinical practice should not be underestimated. Modern techniques of design and fit of soft contact lenses can take into account features of the patient's eye, adapting the back surface of a lens to match the specific elevations of the cornea. These methods require again a detailed corneal topographic analysis of the anterior face of the cornea. 

With the introduction of high-speed videokeratoscopy ^4,5 in the study of the dynamics of corneal surface topography ⁶ and tear film stability, ⁷ data storage needs have become significant, motivating another important application of corneal surface modeling: data compression. ⁸  

Nearly all modern corneal topographers collect the data (elevation, curvature, mire displacement, or others) in a finite and structured set of data points. The data are contaminated by errors that stem from several sources, including the natural device noise, measurement and digitalization errors, algorithm errors (like those converting the mire displacement in elevation), rounding errors. Hence, we face the problem of the parsimonious fit of the actual surface data contaminated by noise, with a minimum number of coefficients or parameters, for clinical and technological applications. 

The solutions to this problem are usually classified as either zonal or modal methods. In the former group, the domain where the data are collected is subdivided into more elementary subdomains (e.g., triangles), and the surface is approximated in each subdomain with an additional requirement of global smoothness. ^{9 –11} In the modal methods of reconstruction, the approximation is found as a (typically linear) combination of functions from the given set or dictionary, defined by several parameters. Zonal methods are flexible and accurate, but they lack the simplicity of the modal approach, they are substantially more computer-intensive, and they encode the final shape in a larger amount of data. 

Crucial decisions to make in the modal reconstruction are the set of functions to use, the value of their parameters, and the number of functions needed to recover the relevant information without fitting (at least, in a large scale) the inevitable error (the so-called model selection problem). A standard functional basis, well-established in ophthalmology for expressing ocular wavefront error, is given by the Zernike polynomials. ¹² The coefficients of the expansions can be interpreted in terms of optical aberrations. As a fitting routine, Zernike polynomials are not limited to analysis of wavefront error surfaces, but can be applied to other ocular surfaces as well, including the anterior corneal surface. ^13,14 Corneal topography diagnostic tools can use the Zernike coefficients as inputs into corneal classification of neural networks, ^15,16 replacing or supplementing the currently used corneal indices embedded in many topography devices. 

However, potential limitations in this approach have been reported in the literature. ^12,17 There is a growing concern that the Zernike fitting method itself may be inaccurate in abnormal conditions. Furthermore, it is very difficult to assess a priori how many terms are necessary to achieve acceptable accuracy in the Zernike reconstruction of any given corneal shape. ¹⁸ It is known ¹⁷ that limiting Zernike analysis to only a few orders may cause incorrect assessment of the severity of the more advanced stages of keratoconus. ¹ This information is particularly needed in the discriminant analysis of the disease markers or when selecting the numerical inputs for neural network–based diagnostic software, such as corneal classification and utilities for grading condition severity. 

Several alternatives to the modal least-squares fit with Zernike polynomials have been suggested recently. Some of them intend to combine the modal and zonal approaches to preserve the best of both worlds ¹⁹ or to use nonlinear methods. ⁸ The possibility of achieving the accuracy and flexibility of the zonal methods within the framework of the linear modal approach by means of localized radial basis functions has been put forward, ²⁰ but without development of an actual implementation or procedure. In this follow-up paper we describe an adaptive, multiscale working algorithm for the parsimonious fit of the surface data, based on residual iteration with knot insertion that allows the adaptation of the number of functions used in the reconstruction to the conditions of each cornea. 

The residual iteration is well known in many branches of mathematics; it is related for instance to the iterative refinement methods of solution of systems of linear equations. In the context of purely radial basis functions, an adaptive greedy approximation algorithm using interpolation has been proposed by Schaback and Wendland. ²¹ The adaptive increase in the number of basis functions as a technique used to improve the quality of a given initial approximation is also standard. (See, for instance, Fasshauer ²² for a general discussion and references.) 

The method that we have developed allows building iteratively an approximation function as a linear combination of anisotropic Gaussian-basis functions, implementing also a dynamic selection of the parameters and the management of noise. It can be used to reconstruct altimetric data, corneal power maps, among other applications. Although it has been tuned for Placido ring–based keratoscopes, with the data nodes located in almost concentric rings, the technique is actually applicable to any scattered data approximation (Martínez-Finkelshtein A, et al. IOVS. 2010;51:ARVO E-Abstract 5690). 

Methods

The General Setting

The input data are a 3D cloud (x_k , y_k , z_k ), with k = 1, …, N, corresponding to either elevation or curvature z_k at the node P_k of the anterior corneal surface with Cartesian coordinates (x_k , y_k ). We will discuss the case in which z_k corresponds to elevation. A standard procedure is to “flatten” the data by fitting it with the best-fit sphere ²³ of the form S(x, y) = z ₀ +

$\sqrt{R^2-{\left(x-x_0\right)}^2-{\left(y-y_0\right)}^2}$

, where R and (x ₀, y ₀, z ₀) are its radius and the Cartesian coordinates of its center, respectively. 

As a result of the previous step, the residual errors ε_k⁽¹⁾ = z_k − S(x_k, y_k) contain the relevant information at different scales and noise. Our aim is to fit these residuals by a function E(x,y) in such a way that an analytic expression for the corneal height is given by  

 

where S accounts for the global shape of the cornea, while E captures the small irregularities on the surface. Function E is a linear combination of n functions from a given dictionary. In an ideal setting, n depends on the actual data and should be large enough to allow fitting all relevant information, but not too large to fit pure noise. Here, we use as basis functions the Gaussians of the form  

 

where the superscript T denotes the matrix transpose, Q = (Q_x, Q_y)^T is a certain point on the plane (“center”), and A is a positive-definite matrix in ℝ^2×2. For such a matrix, the A-norm of a point (column vector) P in R² is defined as  

 

for  

 

with α_x > 0 and α_xα_y > α_xy². 

In general, these are anisotropic radial basis functions that boil down to standard radial basis Gaussian functions (RBGF) when A is a positive multiple of the identity matrix I ₂. 

Hence, we seek the expression of the form  

 

A fitting routine should allow for an adequate selection of all parameters: centers Q ^(j), shape matrices A_j , scaling factors c_j , and number of terms n in the functional representation. 

We propose an iterative algorithm of reconstruction, such that in each step we fit partially the residual error by one anisotropic radial basis Gaussian function (A-RBGF and compute the new residuals, which will become the input for the next iteration (residual iteration with knot insertion). All the parameters will be chosen dynamically, depending on the residual data in each step. A key observation is the possibility of linearizing the problem by using logarithms. We take advantage of the parallelization that modern hardware allows, as well as discuss the appropriate stopping criteria. The details are presented in the following sections. 

Description of the Iterative Algorithm

Let E_{j −1} be already computed (we take E ₀ = 0). The input data for j's iteration (j = 1, 2, …) is the cloud (x_k , y_k , ε _k ^(j)) of nodes P_k = (x_k , y_k ) ^T and the corresponding residuals ε _k ^(j), k = 1, 2, …, N; recall that ε _k ⁽¹⁾ = z_k − S(x_k , y_k ) are the residual errors after the best fit sphere. We perform the following steps. 

Step 1: Selection of the Active Center Q^(j).

The problem of the selection of a center of a radial basis function has been discussed in Jamshidi and Kirby, ²⁴ where the criterion of maximum cross-correlation is used to choose the centers. Another criterion is based on the power function. ^21,22,25 Both methods, although computationally demanding, can be implemented to perform step 1. However, in our practice, we found the much simpler criterion of the maximal residual to be as satisfactory, at a minimum cost; it also correlates with the geometry of the A-RBGF. Hence, we chose Q ^(j) = (x _{k 0} , y _{k 0} )^T , where k ₀ = arg max _k |ε _k ^(j)|, and we denote m ^(j) = ε _{k 0} ^(j). 

Step 2: Dynamic Filtering.

Once the center Q ^(j) has been selected, we check the number, ℓ _k , of nodes P_k lying in the largest disc, centered at Q ^(j) and containing only nodes with the residues of the same sign as m ^(j). If ℓ _k < 20, we consider Q ^(j) to be an outlier, and it cannot be chosen as the center at this iteration. This can be done by setting ε _{k 0} ^(j) = 0 after which we return to step 1. Otherwise, we proceed to the next step. 

Step 3: Selection of the Shape Parameters.

We determine first the influence nodes _j (s), defined as the maximal set of, at most, s nodes P_k closest to Q ^(j), with residues of the same sign as m ^(j). It is convenient to parallelize the subsequent computations for several values of s: We performed experiments using the vector of values s = [s _min, 100, 150, 200, 300], where s _min = min(ℓ _k ,50), with ℓ _k defined in step 2. 

The interpolating conditions ε_k0^(j)h_j(x_k, y_k) = ε_k^(j), k ∈ _j(s), are equivalent to the overdetermined linear system  

 

in the three unknown entries of the shape matrix  

 

We can solve this system in the sense of weighted linear least squares (WLS), where the kth equation is multiplied by the weight η_k = (1 + ‖P_k − Q^(j)‖²)⁻¹. This solution is obtained by standard methods.²⁶ Observe also that  

 

 

As an additional step, we examine the matrix A_j built out of the solution of equation 3. If it is not positive definite, we force h_j to be a standard radial basis function, A_j = αI ₂. In this case, equation 3 is reduced to  

 

the solution of which is computed in the sense of the WLS. An explicit formula for α can be derived, but it is omitted here for the sake of brevity. 

Step 4: Selection of the Scaling Factor.

We can calculate the coefficient c_j from  

 

using the WLS with the same weights η _j as before. Again, c_j is computed by an explicit formula that can be easily derived, and in consequence, the formula is omitted here. It should be noted, however, that in many cases the simpler interpolation condition c_k = m ^(j) yields comparable results. 

Step 5: Computation of the New Residuals.

Having found the values of c_j and A_j, we update  

 

As mentioned, all the computations are performed in parallel for different values of s and hence, different nested sets of influence nodes _j(s). We now keep the value of s (and the corresponding values of c_j and A_j) that yields the smallest norm of the residue vector (ε_k^(j+1)). As a result, we find the new approximation E = E_j−1 + c_jh_j. As a final step, we check the stopping criterion, which will now be discussed. If this is not satisfied, we return to step 1. 

Stopping Criteria

In theory, the algorithm run indefinitely yields an interpolating function, and in consequence, a 0 residue vector. Measurement errors may influence the model selection: We want to capture all the relevant information without fitting the noise. Many solutions to this problem are described in the literature. The choice of the number of Zernike polynomials for the modal reconstruction of the altimetric data has been discussed, ^27,28 but the suggested algorithms are typically computationally intensive. ^28,29  

Less demanding methods use information theory criteria, such as the Akaike information criterion (AIC) or the efficient detection criterion (EDC).³⁰ However, we can gain information analyzing directly the behavior of the normalized mean square errors MSE_j defined by  

 

Typically, these errors start decaying at an exponential rate. After several iterations, we observe a stabilization in this rate of decay that becomes linear. This normally happens when the values of MSE_j are between 10⁻³ and 10⁻⁴ μm². Based on this experience, we have successfully used the following stopping criterion: We allow the algorithm to run for up to 50 iterations and calculate the weighted slopes  

 

The sequence δMSE, although oscillatory, is negative and tends to 0, and so we find the last iteration 1 ≤ J₁ ≤ 50 when δMSE_j ≥ −0.02. If MSE_j < 10⁻³ μm, we fix J₁ + 1 as the stopping iteration. Otherwise, we seek the last iteration 1 ≤ J₂ ≤ 50 when δMSE_j ≥ −0.01, and stop the algorithm at the (J₂ + 1)th iteration. 

Experimental Data

The altimetric and curvature data from in vivo corneas used for the experiments described below were collected by the Placido-based topography system (CSO, Firenze, Italy), which in ideal conditions digitizes up to 24 rings with 256 equally distributed points on each mire. All procedures were implemented in commercial software (MatLab 7; The MathWorks, Natick, MA) and run on standard platforms. 

Results

In this section, we present numerical results obtained by our method and by the standard Zernike reconstruction, applied both to simulated and real cornea surfaces. The choice of the Zernike-based method for comparison is motivated by its modal character and by being a standard de facto in industry and clinical practice. 

We analyze first the global error, measured by MSE_j , which tends to decrease monotonically with the exponential rate. A typical behavior for the reconstruction of the altimetric data for both synthetic and real corneas can be observed in Figure 1. In particular, in all cases we appreciate the clear transition from an overexponential to a linear rate of decay, which is used as the stopping criterion. For comparison, we have reconstructed the same data with the Zernike polynomials using linear least squares, which is the standard procedure in clinical practice, implemented in most topographers. The MSE_j for the Zernike reconstruction is plotted in the same Cartesian coordinate system, where j indicates the total number of Zernike polynomials used. Recall that j varies from 1 to (m + 1)(m + 2)/2, where m is the maximum radial order used. 

Figure 1.

 

Comparison of the MSE_j for altimetric data reconstruction by Zernike polynomials and by the A-RBGF algorithm. Top left: a synthetic cornea made of a linear combination of 10 Gaussian functions. Top right: a normal cornea. Bottom row: keratoconic corneas. Arrows: stopping time of the algorithm according to the criterion described in the section on Stopping Criterion. The horizontal axis (n) indicates the number of functions used.

Zernike polynomials easily capture the global shape of the surface, which is expressed in a typical fast decay of the corresponding error. However, smaller details on the surface (such as areas of localized steepening) are much less suited for this tool. It explains the saturation observed in the Zernike error behavior after a few orders, typically, around 25 polynomials. Such saturation is not the case in the reconstruction with A-RBGF, whose multiscale, adaptive character allows adjustment of the parameters adequately in each iteration. 

However illuminating these graphs may be, the global error is not the only way to compare both reconstruction approaches. Recall that the modal method with Zernike polynomials is tailored precisely to achieve a maximum reduction of the MSE_j for each j, whereas the iterative algorithm presented here has a totally different goal: locating the most salient feature of the residual surface and incorporating it into the analytic expression in equation 1. 

This observation can be illustrated by fitting a synthetic cornea with a simulated scar, used earlier in Martínez-Finkelshtein at al. ²⁰ Its contour plot is depicted in Figure 2, top left. The top right panel shows the contour plot of the surface reconstructed with the adaptive procedure described here using as few as 20 functions. The other two contour plots correspond to the same surface reconstructed with 36 (order ≤ 7) and 136 (order ≤ 15) Zernike polynomials. 

Figure 2.

 

Comparison of the reconstruction of a synthetic “scar” (top left) with 20 iterations of the adaptive algorithm (top right) and 36 (bottom left) and 136 (bottom right) Zernike polynomials.

The situation becomes even more apparent if we compare the residual errors along the eighth mire for both methods (Fig. 3, left). Although the Zernike polynomials work perfectly on smooth portions of the surface, they find difficulties in adapting to fast varying shapes, whereas the A-RBGF algorithm uses its multiscale and adaptive character to fit the surface almost perfectly after a few iterations. It takes a long time for the MSE errors of Zernike polynomials to progress, as illustrated in Figure 3, right. 

Figure 3.

 

Left: residual errors along the eighth mire for the synthetic “scar” reconstructed with 20 iterations of the adaptive algorithm (bold line) and up to 36 Zernike polynomials (dotted line). Right: MSE_j for the “scar” reconstruction by Zernike polynomials and by the A-RBGF algorithm.

Regarding the stopping criterion, the experiments performed with real and synthetic corneas show that the reasonable number of iterations lies between 20 and 40; there is no clear correlation between the number of iterations and the condition of the cornea, as Figure 1 shows. This is why we consider it appropriate to perform 50 iterations (taking advantage of the speed of the algorithm) to decide the correct value for n in equation 2. 

However, more correlation exists with the location and grouping of the centers Q ^(j). For the normal corneas, the centers typically cluster at the border of the area, where most of the oscillations occur, whereas for corneas affected by keratoconus, we observe how some centers match the deformation already at the first iterations. 

The adaptive algorithm described is suitable for a reliable reconstruction of any surface for the discrete set of data. In particular, it is possible also to reconstruct a corneal power map or the wavefront. Taking into account the typical shape of such a surface, it is convenient here to skip the fit by a sphere or Zernike polynomials of a low order, making S ≡ 0 in equation 1. 

An alternative way of assessing the quality of the approximation rendered by our method is by analyzing the resulting point spread function (PSF). Figure 4 shows the effect of modeling the corneal surface data of a simulated highly irregular eye on the estimated PSF. The original PSF corresponding to a wavefront described by an expansion in 136 Zernike polynomials, and a corneal diameter of 8 mm is shown (top left), together with a Zernike polynomial approximation of radial order ≤ 5 (top right) and ≤ 9 (bottom left) (21 and 55 functions, respectively), and the A-RBGF approximation with 21 functions (bottom right). Clearly, the latter leads to a PSF that more closely resembles the original PSF (capturing some finer features) than that derived from the Zernike polynomials. 

Figure 4.

 

Original PSF for a simulated cornea with a pupil size 8 mm in diameter and high wavefront error (top left), and the corresponding PSFs from a reconstruction with first 21 (order ≤ 5) Zernike polynomials (top right), first 55 (order ≤ 9) Zernike polynomials (bottom left), and with 21 A-RBGF functions (bottom right).

Discussion

In this work, we have developed an adaptive fitting method for corneal data, combining modal simplicity with the advantages of zonal reconstruction. It consists of a preliminary fit of the data with a sphere and an iterative procedure that adds terms to the analytic representation of the corneal data. Each term consists of a scaled anisotropic radial basis Gaussian function. The coefficients are computed dynamically and allow fitting the data in each iteration, independent of the scale. The method comprises also a filtering procedure that discards the outliers (data clearly corresponding to measurement noise) and a stopping criterion that chooses the final number of functions in the analytic representation in accordance with the evolution of the residual error. 

Although it is computationally more intense in comparison with the standard Zernike fit, the numerical implementation of this algorithm in a standard personal computer is very fast and its execution times are negligible, as illustrated by Tables 1 and 2. Experimental results allow us to draw the following conclusions:

•  
  The least-squares approximation by a linear combination of Zernike polynomials of a radial order up to six (which is the standard in modern aberrometers ³¹ ) fits adequately the altimetric data in the case of a normal cornea. It can be used also to capture the major features of the shape of the surface. However, for strongly aberrated corneas, the Zernike-based procedure saturates relatively early, and we need a high number of terms to achieve the desired accuracy at regions of localized steepening, at the price of overparameterizing the model and fitting the measurement noise.
•  
  In contrast, the iterative method presented herein exhibits a steady exponential error decay, independent of the complexity of the cornea. Its actual rate is basically influenced by the residue distribution. The fast decrease in the first iterations, when all salient features are reconstructed, is followed by a stable linear decay, when essentially errors are being fit. This can be used as a stopping criterion for the iterations. In this way, the minimal number of functions in the analytic representation of the cornea is used.
•  
  Unlike in the case of the Zernike polynomials, the complexity of each term in the functional representation with A-RBGF is the same, but their parameters vary to fit the current scale. This scale is determined only by the residual errors and not by the number of the iteration.
•  
  Because of the localized character of A-RBGF, the position and clustering of their centers, as well as the size of the shape parameters, provides an additional spatial information about the regions of higher irregularity.

Table 1.

 

Reconstruction Time for the Surfaces in Figure 1

Surface	Top Left	Top Right	Bottom Left	Bottom Right
Functions used, n	15	28	45	21
Computation time (A-RBGF), s	0.3169	0.3537	0.4819	0.2671
Computation time (Zernike), s	0.0632	0.0550	0.1183	0.0373

 In each case, the same number of functions is used for the A-RBGF algorithm and the standard Zernike fitting. This number is of the form (m + 1)(m + 2)/2, with m natural, and is selected as the closest to the recommended stopping value, indicated in each figure, with the purpose to include all Zernike polynomials of the maximal radial order.

Table 2.

 

Time of Processing for the Three Fitting Procedures Used in Figure 2

Method	A-RBGF Algorithm, with 20 Functions (Fig. 2, Top Right)	Zernike Fitting, with 36 Polynomials (Fig. 2, Bottom Left)	Zernike Fitting, with 136 Polynomials (Fig. 2, Bottom Right)
Computational time, s	0.5116	0.0990	0.7210

The iterative adaptive algorithm for the cornea modeling proposed here provides a method of obtaining a compact mathematical description of the shape or power map of the cornea. All information is ultimately encoded in the following set of values: the center and radius of the best-fit sphere, plus the center locations, shape parameters, and scaling factors, which, again, need a negligible amount of data storage space. This description can be used for global visualization of the cornea or of its portions, capturing smaller details than with standard procedures. It can serve also as the input data for resampling and computation of some other relevant values via ray tracing, numerical integration, and other calculations. 

Footnotes

 Supported in part by Junta de Andalucía Grant P06-FQM-01738 (AM-F, DRL); the Ministry of Science and Innovation of Spain Project MTM2008-06689-C02-01 (AM-F); Junta de Andalucía Grant P09-FQM-4643 (AM-F); and Carlos III Health Institute of the Ministry of Science and Innovation of Spain, Grants PI08/90519 and PI10/01843.

Footnotes

 Disclosure: A. Martínez-Finkelshtein, None; D.R. López, None; G.M. Castro, None; J.L. Alió, None