**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.

^{ 1 }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.

^{ 2 }

^{ 3 }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.

^{ 4,5 }in the study of the dynamics of corneal surface topography

^{ 6 }and tear film stability,

^{ 7 }data storage needs have become significant, motivating another important application of corneal surface modeling: data compression.

^{ 8 }

^{ 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.

^{ 12 }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.

^{ 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.

^{ 18 }It is known

^{ 17 }that limiting Zernike analysis to only a few orders may cause incorrect assessment of the severity of the more advanced stages of keratoconus.

^{ 1 }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.

^{ 19 }or to use nonlinear methods.

^{ 8 }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,

^{ 20 }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.

^{ 21 }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

^{ 22 }for a general discussion and references.)

*IOVS*. 2010;51:ARVO E-Abstract 5690).

*x*,

_{k}*y*,

_{k}*z*), with

_{k}*k*= 1, …,

*N*, corresponding to either elevation or curvature

*z*at the node

_{k}*P*of the anterior corneal surface with Cartesian coordinates (

_{k}*x*,

_{k}*y*). We will discuss the case in which

_{k}*z*corresponds to elevation. A standard procedure is to “flatten” the data by fitting it with the best-fit sphere

_{k}^{ 23 }of the form

*S*(

*x*,

*y*) =

*z*

_{0}+

*R*and (

*x*

_{0},

*y*

_{0},

*z*

_{0}) are its radius and the Cartesian coordinates of its center, respectively.

_{k}

^{(1)}=

*z*−

_{k}*S*(

*x*,

_{k}*y*) contain the relevant information at different scales and noise. Our aim is to fit these residuals by a function

_{k}*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}*is a certain point on the plane (“center”), and*

^{T}*A*is a positive-definite matrix in ℝ

^{2×2}. For such a matrix, the

*A*-norm of a point (column vector)

*P*in

*R*

^{2}is defined as for with α

*> 0 and α*

_{x}*α*

_{x}*> α*

_{y}_{xy}

^{2}.

*A*is a positive multiple of the identity matrix

*I*

_{2}.

*E*

_{j}_{−1}be already computed (we take

*E*

_{0}= 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}*and the corresponding residuals ε*

^{T}_{ k }

^{(j)},

*k*= 1, 2, …,

*N*; recall that ε

_{ k }

^{(1)}=

*z*−

_{k}*S*(

*x*,

_{k}*y*) are the residual errors after the best fit sphere. We perform the following steps.

_{k}^{(j)}.

^{ 24 }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 }

*)*, where

^{T}*k*

_{0}= arg max

*|ε*

_{k}_{ k }

^{(j)}|, and we denote

*m*

^{(j)}= ε

_{ k 0 }

^{(j)}.

**Q**

^{(j)}has been selected, we check the number, ℓ

*, of nodes*

_{k}*P*lying in the largest disc, centered at

_{k}**Q**

^{(j)}and containing only nodes with the residues of the same sign as

*m*

^{(j)}. If ℓ

*< 20, we consider*

_{k}**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.

*(*

_{j}*s*), defined as the maximal set of, at most,

*s*nodes

*P*closest to

_{k}**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(ℓ

*,50), with ℓ*

_{k}*defined in step 2.*

_{k}_{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

*k*th equation is multiplied by the weight η

*= (1 + ‖*

_{k}*P*−

_{k}*Q*

^{(j)}‖

^{2})

^{−1}. This solution is obtained by standard methods.

^{26}Observe also that

*A*built out of the solution of equation 3. If it is not positive definite, we force

_{j}*h*to be a standard radial basis function,

_{j}*A*

_{j}= αI_{2}. 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.

*c*from using the WLS with the same weights η

_{j}*as before. Again,*

_{j}*c*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

_{j}*c*=

_{k}*m*

^{(j)}yields comparable results.

*c*and

_{j}*A*, we update As mentioned, all the computations are performed in parallel for different values of

_{j}*s*and hence, different nested sets of influence nodes

*(*

_{j}*s*). We now keep the value of

*s*(and the corresponding values of

*c*and

_{j}*A*) that yields the smallest norm of the residue vector (ε

_{j}_{k}

^{(j+1)}). As a result, we find the new approximation

*E*=

*E*

_{j}_{−1}+

*c*. As a final step, we check the stopping criterion, which will now be discussed. If this is not satisfied, we return to step 1.

_{j}h_{j}^{ 27,28 }but the suggested algorithms are typically computationally intensive.

^{ 28,29 }

^{30}However, we can gain information analyzing directly the behavior of the normalized mean square errors

*MSE*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

_{j}*MSE*are between 10

_{j}^{−3}and 10

^{−4}μm

^{2}. 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*

_{1}≤ 50 when δ

*MSE*≥ −0.02. If

_{j}*MSE*< 10

_{j}^{−3}μm, we fix

*J*

_{1}+ 1 as the stopping iteration. Otherwise, we seek the last iteration 1 ≤

*J*

_{2}≤ 50 when δ

*MSE*≥ −0.01, and stop the algorithm at the (

_{j}*J*

_{2}+ 1)th iteration.

*MSE*, 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

_{j}*MSE*for the Zernike reconstruction is plotted in the same Cartesian coordinate system, where

_{j}*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.

*MSE*for each

_{j}*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.

^{ 20 }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.

*n*in equation 2.

**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.

*S*≡ 0 in equation 1.

- 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
^{ 31 }) 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.

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 |

*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.