July 2011
Volume 52, Issue 8
Free
Cornea  |   July 2011
Adaptive Cornea Modeling from Keratometric Data
Author Affiliations & Notes
  • Andrei Martínez-Finkelshtein
    From the Department of Statistics and Applied Mathematics, Almería University, Almería, Spain;
    the Institute Carlos I of Theoretical and Computational Physics, Granada University, Granada, Spain;
  • Darío Ramos López
    From the Department of Statistics and Applied Mathematics, Almería University, Almería, Spain;
  • Gracia M. Castro
    the Ophthalmological Institute VISSUM-Almería, Almería, Spain;
  • Jorge L. Alió
    the Department of Cornea and Refractive Surgery, Vissum Instituto Oftalmológico de Alicante, Alicante, Spain; and
    the Division of Ophthalmology, Miguel Hernández University, Medical School, Alicante, Spain.
  • Corresponding author: Andrei Martínez-Finkelshtein, Depto. de Estadística y Matemática Aplicada, Universidad de Almería, Carretera de Sacramento s/n, 04120 Almería, Spain; andrei@ual.es
Investigative Ophthalmology & Visual Science July 2011, Vol.52, 4963-4970. doi:10.1167/iovs.10-6774
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Andrei Martínez-Finkelshtein, Darío Ramos López, Gracia M. Castro, Jorge L. Alió; Adaptive Cornea Modeling from Keratometric Data. Invest. Ophthalmol. Vis. Sci. 2011;52(8):4963-4970. doi: 10.1167/iovs.10-6774.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
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 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  
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. 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. 
With the introduction of high-speed videokeratoscopy 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  
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. 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. 
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. 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. 
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 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. 
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. 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.) 
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 (xk , yk , zk ), with k = 1, …, N, corresponding to either elevation or curvature zk at the node Pk of the anterior corneal surface with Cartesian coordinates (xk , yk ). We will discuss the case in which zk corresponds to elevation. A standard procedure is to “flatten” the data by fitting it with the best-fit sphere 23 of the form S(x, y) = z 0 +
R 2 ( x x 0 ) 2 ( y y 0 ) 2
, where R and (x 0, y 0, z 0) are its radius and the Cartesian coordinates of its center, respectively. 
As a result of the previous step, the residual errors εk(1) = zkS(xk, yk) 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 = (Qx, Qy)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 R2 is defined as   for   with αx > 0 and αxαy > αxy2
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 2
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 Aj , scaling factors cj , 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 Ej −1 be already computed (we take E 0 = 0). The input data for j's iteration (j = 1, 2, …) is the cloud (xk , yk , ε k (j)) of nodes Pk = (xk , yk ) T and the corresponding residuals ε k (j), k = 1, 2, …, N; recall that ε k (1) = zk S(xk , yk ) 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, 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 )T , where k 0 = 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 Pk 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 Image not available j (s), defined as the maximal set of, at most, s nodes Pk 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)hj(xk, yk) = εk(j), kImage not availablej(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 + ‖PkQ(j)2)−1. This solution is obtained by standard methods.26 Observe also that    
As an additional step, we examine the matrix Aj built out of the solution of equation 3. If it is not positive definite, we force hj to be a standard radial basis function, Aj = α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. 
Step 4: Selection of the Scaling Factor.
We can calculate the coefficient cj from   using the WLS with the same weights η j as before. Again, cj 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 ck = m (j) yields comparable results. 
Step 5: Computation of the New Residuals.
Having found the values of cj and Aj, we update   As mentioned, all the computations are performed in parallel for different values of s and hence, different nested sets of influence nodes Image not availablej(s). We now keep the value of s (and the corresponding values of cj and Aj) that yields the smallest norm of the residue vector (εk(j+1)). As a result, we find the new approximation E = Ej−1 + cjhj. 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).30 However, we can gain information analyzing directly the behavior of the normalized mean square errors MSEj 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 MSEj are between 10−3 and 10−4 μm2. 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 ≤ J1 ≤ 50 when δMSEj ≥ −0.02. If MSEj < 10−3 μm, we fix J1 + 1 as the stopping iteration. Otherwise, we seek the last iteration 1 ≤ J2 ≤ 50 when δMSEj ≥ −0.01, and stop the algorithm at the (J2 + 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 MSEj , 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 MSEj 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 MSEj 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.
Figure 1.
 
Comparison of the MSEj 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 MSEj 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. 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. 
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.
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: MSEj for the “scar” reconstruction by Zernike polynomials and by the A-RBGF algorithm.
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: MSEj 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).
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 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.
Table 1.
 
Reconstruction Time for the Surfaces in Figure 1
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
Table 2.
 
Time of Processing for the Three Fitting Procedures Used in Figure 2
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
References
Schwiegerling JT Greivenkamp JE . Keratoconus detection based on videokeratoscopic height data. Optom Vision Sci. 1989;73;721–728. [CrossRef]
Artal P Berrio E Guirao A Piers P . Contribution of the cornea and internal surfaces to the change of ocular aberrations with age. J Opt Soc Am A. 2002;19:137–143. [CrossRef]
Maeda N Klyce SD Smolek MK . Comparison of methods for detecting keratoconus using videokeratography. Arch Ophthalmol. 1995;113(7):870–874. [CrossRef] [PubMed]
Nemeth J Erdelyi B Csakany B . High-speed videotopographic measurement of tear film build-up time. Invest Ophthalmol Vis Sci. 2002;43(6):1783–1790. [PubMed]
Iskander DR Collins MJ , Applications of high-speed videokeratoscopy. Clin Exp Optom. 2005;88(4):223–231. [CrossRef] [PubMed]
Zhu M Collins MJ Iskander DR . Dynamics of ocular surface topography. Eye. 2007;21:624–632. [PubMed]
Iskander DR Collins MJ Davis B . Evaluating tear film stability in the human eye with high-speed videokeratoscopy. IEEE Trans Biomed Eng. 2005;52(11):1939–1949. [CrossRef] [PubMed]
Schneider M Iskander DR Collins MJ . Modeling corneal surfaces with rational functions for high-speed videokeratoscopy data compression, IEEE Trans Biomed Eng. 2009;56(2):493–499. [CrossRef] [PubMed]
Halstead MA Barsky B Klein SA Mandell RB . A spline surface algorithm for reconstruction of corneal topography from a videokeratographic reflection pattern. Optom Vision Sci. 1995;72(11):821–827. [CrossRef]
Liu X Gao Y . B-spline based wavefront reconstruction for lateral shearing interferometric measurement of engineering surfaces. Adv Abras Technol. 2003;238–239:169–174.
Ares M Royo S . Comparison of cubic B-spline and Zernike-fitting techniques in complex wavefront reconstruction. Appl Opt. 2006;45:6945–6964. [CrossRef]
Klyce SD Karon MD Smolek MK . Advantages and disadvantages of the Zernike expansion for representing wave aberration of the normal and aberrated eye. J Refract Surg. 2004;20:S537–S541. [PubMed]
Smolek MK Klyce SD . Zernike polynomial fitting fails to represent all visually significant corneal aberrations. Invest Ophthalmol Vis Sci. 2003;44(11):4676–4681. [CrossRef] [PubMed]
Carvalho LA . Accuracy of Zernike polynomials in characterizing optical aberrations and the cornel surface of the eye. Invest Ophthalmol Vis Sci. 2005:46:1915–1926. [CrossRef] [PubMed]
Carvalho LA . Preliminary results of neural networks and Zernike polynomials for classification of videokeratography maps. Optom Vision Sci. 2005;82(2):151–158. [CrossRef]
Smolek MK Klyce SD . Current keratoconus detection methods compared with a neural network approach, Invest Ophthalmol Vis Sci. 1997;38(11):2290–2299. [PubMed]
Smolek MK Klyce SD . Goodness-of-prediction of Zernike polynomial fitting to corneal surfaces, J Cataract Refract Surg. 2005;31:2350–2355. [CrossRef] [PubMed]
Iskander R Collins MJ Davis B . Optimal modeling of corneal surfaces with Zernike polynomials. IEEE Trans Biomed Eng. 2001;48(1):87–95. [CrossRef] [PubMed]
Espinosa J Mas D Perez J Illueca C . Optical surface reconstruction technique through combination of zonal and modal fitting, J Biomed Opt 2010;15(2):026022.
Martínez-Finkelshtein A Delgado AM Castro GM Zarzo A Alió JL . Comparative analysis of some modal reconstruction methods of the shape of the cornea from the corneal elevation data, Invest Ophthalmol Vis Sci. 2009;50(12):5639–5645. [CrossRef] [PubMed]
Schaback R Wendland H . Adaptive greedy techniques for approximate solution of large RBF systems. Numer Algorithms. 2000;24(3):239–254. [CrossRef]
Fasshauer GE . Meshfree Approximation Methods with MATLAB. Singapore: World Scientific Publishers; 2007.
Ahn SJ Rauh W Warnecke H-J . Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola, Pattern Recogn. 2001;4(12)2283–2303. [CrossRef]
Jamshidi AA Kirby MJ . Towards a black box algorithm for nonlinear function approximation over high-dimensional domains. SIAM J Sci Comput 2007;29(3):941–963. [CrossRef]
De Marchi S Schaback SR Wendland H . Near-optimal data-independent point locations for radial basis function interpolation. Adv Comput Math. 2005;23(3):317–330. [CrossRef]
Trefethen LN Bau DIII . Numerical Linear Algebra. Philadelphia: Society for Industrial and Applied Mechanics (SIAM); 1997.
Iskander R Morelande MR Collins MJ Davis B . Modeling of corneal surfaces with radial polynomials, IEEE Trans Biomed Eng. 2002;49(4):320–328. [CrossRef] [PubMed]
Alonso-Caneiro D Iskander R Collins MJ . Estimating corneal surface topography in videokeratoscopy in the presence of strong signal interference. IEEE Trans Biomed Eng. 2008;55(10):2381–2387. [CrossRef] [PubMed]
Iskander DR Alkhaldi W Zoubir AM . On the computer intensive methods in model selection, in: Proceedings of the 33rd IEEE International Conference on Acoustic Speech Signal Process. (ICASSP). 2008:3461–3464.
Rao CR Wu Y . A strongly consistent procedure for model selection in a regression problem. Biometrika. 1989;76(2):369–374. [CrossRef]
Rozema JJ Van Dyck DEM Tassignon M-J . Clinical comparison of 6 aberrometers, Part 1: technical specifications J Cataract Refract Surg. 2005;31:1114–1127. [CrossRef] [PubMed]
Figure 1.
 
Comparison of the MSEj 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.
Figure 1.
 
Comparison of the MSEj 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.
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.
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.
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: MSEj for the “scar” reconstruction by Zernike polynomials and by the A-RBGF algorithm.
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: MSEj for the “scar” reconstruction by Zernike polynomials and by the A-RBGF algorithm.
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).
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).
Table 1.
 
Reconstruction Time for the Surfaces in Figure 1
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
Table 2.
 
Time of Processing for the Three Fitting Procedures Used in Figure 2
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
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×