Comparative Analysis of Some Modal Reconstruction Methods of the Shape of the Cornea from Corneal Elevation Data

Abstract

Purpose.: A comparative study of the ability of some modal schemes to reproduce corneal shapes of varying complexity was performed, by using both standard radial polynomials and radial basis functions (RBFs). The hypothesis was that the correct approach in the case of highly irregular corneas should combine several bases.

Methods.: Standard approaches of reconstruction by Zernike and other types of radial polynomials were compared with the discrete least-squares fit (LSF) by the RBF in three theoretical surfaces, synthetically generated by computer algorithms in the absence of measurement noise. For the reconstruction by polynomials, the maximal radial order 6 was chosen, which corresponds to the first 28 Zernike polynomials or the first 49 Bhatia-Wolf polynomials. The fit with the RBF was performed by using a regular grid of centers.

Results.: The quality of fit was assessed by computing for each surface the mean square errors (MSEs) of the reconstruction by LSF, measured at the same nodes where the heights were collected. Another criterion of the fit quality used was the accuracy in recovery of the Zernike coefficients, especially in the case of incomplete data.

Conclusions.: The Zernike (and especially, the Bhatia-Wolf) polynomials constitute a reliable reconstruction method of a nonseverely aberrated surface with a small surface regularity index (SRI). However, they fail to capture small deformations of the anterior surface of a synthetic cornea. The most promising approach is a combined one that balances the robustness of the Zernike fit with the localization of the RBF.

Zernike analysis is used commonly in ophthalmology to express ocular wavefront error in the form of a polynomial function. ¹ The coefficients of these expansions have interpretation in terms of the basic aberrations such as defocus, astigmatism, coma, trefoil, and spherical aberrations, along with higher order 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. ^2,3 It has been suggested that Zernike analysis may be applicable in the development of corneal topography diagnostic tools (e.g., Zernike coefficients as inputs into corneal classification of neural networks [Smolek MK, et al. IOVS 1997;38:ARVO Abstract 4298] ⁴ ), replacing or supplementing the currently used corneal indices included with many topography devices. 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 corneal topography data is clinically useful and important. It has been recognized that the corneal front surface generally provides the bulk of the ocular aberrations in the postsurgical or pathologic eye. ⁶  

However, several potential limitations in this approach have been reported in the literature. ^1,7 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 decease markers, or when selecting the numerical inputs for neural network–based diagnostic software such as corneal classification and grading utilities for condition severity. 

In this sense, several alternatives to Zernike polynomials have been recently suggested. 

This is a report of a comparative study of the ability of some modal approaches to reproduce corneal shapes of varying complexity. Rather than dwelling further on the shortcomings of the Zernike fit, we compare several techniques in some “model” situations, ignoring on purpose all sources of noise that exist in any real system. In this study we avoided experiments using third-party software on corneal elevation from in vivo eyes, but implement the fitting methods on theoretical surfaces, synthetically generated by computer algorithms. This method gives an insight into the intrinsic accuracy of each approach. 

It should be emphasized that our primary goal was to assess the behavior of some methods in different situations. As a result of our study, it may be concluded that there is no unique and best approach to corneal surface reconstruction that can be considered preferable in every scenario, and that a combination of techniques may therefore be the optimal strategy. 

In corneal topography, the elevation of the corneal surface is collected on a discretely sampled grid, which is typically a polar grid for Placido-ring–based systems. These raw data are used to reconstruct the corneal shape, by applying either zonal (see e.g., Ref. 9) or modal algorithms. The modal approach is taken most often because it is easy to use and it offers better noise-suppressing properties. Within the modal approach, the anterior surface of the cornea can be modeled by a linear combination of some basis functions,  

 

where C(P) is the corneal elevation at the point P of the plane, f_j is the a priori chosen basis function, and a = (a₁,…,a_v)^T is the expansion coefficient (the superindex T means matrix transposition). In this setting, fitting the equation to a discrete set of elevation data Z = (Z₁,…,Z_N), N ≥ v, at the nodes P_i, i = 1,…,N, can be restated in terms of solution of the overdetermined linear system  

 

The LSF corresponds formally to the solution of the normal equations M^TMa = M^TZ, which is unique if the collocation matrix M is of maximum rank v. 

Different basis functions f_j in equation 1 can be used, such as radial polynomials (e.g., Zernike, Bhatia-Wolf), Fourier series, and radial basis functions (RBFs). 

Zernike polynomials Z_n ^m (corresponding to the radial non-negative integer n and azimuthal integer m indices, with |m|≤n and n − m even) exhibit special properties that make them an interesting expansion set for the description of general surfaces in the fields of optical engineering and physiological optics. They form a complete set of orthonormal polynomials on the unit disc with respect to the Lebesgue (plane) measure. Since they are well-known, we omit their explicit description, referring the reader to the standard. ¹⁰  

There are several methods of mapping the double indices (n, m) into a 1-D array j. The most widely acceptable one is  

 

Choosing in equation 1 as f_j the Zernike polynomials of radial order ≤ n yields the value v = (n = 1)(n = 2)/2. Most clinical aberrometers use Zernike expansion up to the 6th (typically, the 4th) radial order to reconstruct wavefront data or a corneal surface, ¹¹ corresponding to values v = 15 and v = 28 for n = 4 and n = 6, respectively. It was shown in Iskander et al. ⁸ that, at least for the normal and astigmatic corneas, the optimal value is v = 11 so that even n = 4 leads in most cases to overparameterization of the model. 

Closely related with the Zernike polynomials are the Bhatia-Wolf polynomials, ¹² B_n ^m , whose fitting properties have been analyzed in Iskander et al. ¹³ They also satisfy the orthonormality condition with respect to the unit Lebesgue measure on the disc. One difference between Bhatia-Wolf and Zernike polynomials is that the only constraint |m|≤n on the radial and azimuthal indices results in the generation of (n = 1)² linearly independent polynomials for a given radial degree n instead of the (n = 1)(n = 2)/2 for the Zernike polynomials. It should be noted that B_n ^m are not algebraic polynomials in the Cartesian variables x and y, but they expand in series of monomials xⁱy^jρ^k , ρ = (x ² = y ²)^1/2 for i ≥ 0, j ≥ 0, and k ≥ 0. The double indices (n, m) are easily converted into the polynomial order j of B_j by j = n(n = 1) = m.  

There are other possible choices of radial polynomials, such as the generalized Zernike polynomials ¹⁴ and the Sobolev orthogonal polynomials on the disc. ^15,16  

A special note is owed to another well-known fitting method based on the (bidimensional) Fourier transform, ¹⁷ which reconstructs wavefront data by decomposing the image into spatial frequency components (see e.g., Refs. 18–21). Standard Fourier methods build the surface as a combination of the trigonometric basis whose coefficients can be computed via the FFT algorithm. In some situations, the input information is the set of slopes and not the elevations, in which case an additional step (reduction to the laplacian) is needed. The typical Gibbs phenomenon (high oscillation at the boundary) is handled via a Gershberg-type iterative method (see the literature mentioned thus far). 

In this study, we investigated an alternative meshless technique for reconstructing the corneal shape from the elevation data using as f_j in equation 1 sets of RBFs, defined in their simplest form by translates of a given function Φ:  

 

where points Q _j , called centers of the RBFs, are conveniently chosen, and ‖·‖ denotes the Euclidean distance on the plane. The general theory of interpolation by RBFs is developing rapidly, and several criteria for Φ can be found in the literature. ²² In particular, standard options are the so-called Gaussians and inverse multiquadrics, corresponding, respectively, to Φ(t) = exp(−at) and Φ(t) = (t = c ²)^−β, with positive parameters a, c, and β. However, we are unaware of any deep theoretical analysis of the LSF with RBF. This, according to Buhmann, ²³ is a highly relevant and interesting field of research. 

There are several advantages in the use of equation 1 with the RBF. Because of the fast decay of the Gaussians or multiquadrics, functions f_j in equation 3 are practically locally supported. Hence, equation 1 exhibits features of the zonal approach, eventually capturing small deformations of the surface, which are missed by the polynomial fitting. The rate of decay or the size of the effective support of f_j can be controlled with the parameters of the RBFs, endowing the model with a flexibility that lacks in other modal schemes described earlier. The correct selection of these parameters depends on several factors, such as the frequency of the sampling data, the separation between centers of the RBFs and the grade of variation of the surface. As far as we are aware of, the only work in which such a use of the Gaussians has been discussed, but in the context of the wavefront fitting, is that of Montoya-Hernández et al. ²⁴  

We want to point out that the choice of the RBF in equation 1 does not imply renouncement of the Zernike coefficients as the output information of the reconstructed surface. On the contrary, since the centers Q_j are fixed a priori, the values of  

 

with f_j given by equation 3, can be computed and stored so that the Zernike coefficients are easily recovered by the scalar product of two vectors. 

Material and Methods

It is well known that in a situation close to ideal, when the corneal surface presents only small and smooth deviations from a sphere, almost any reasonable fitting scheme renders good results. In particular, in such a situation, the use of the Zernike polynomials is perfectly justified. Hence, to assess the fitting properties of the different approaches, we have chosen the following three model surfaces with high surface regularity indices (SRIs):

1.  
   Surface A: a “flat sphere”, roughly simulating a surgically altered cornea and a surface with a gradient discontinuity (see Fig. 4).
2.  
   Surface B: a sphere with a radial deformation, or “scar” (Fig. 1).
3.  
   Surface C: a cornea with topographic asymmetry and decentered corneal apex (keratoconus), but with an incomplete set of data (see Fig. 6).

In cases A and B we obtain the elevation data at a discrete set of points P_s , s = 1,…,6144, with polar coordinates (ρ _i , θ _j ) where ρ _i = i/24, i = 1,…,24, whereas the meridians θ _j , j = 1,…,256 are equidistributed in (0, 2π). For the surface C, we collect the elevations at a subset of the nodes just described (Fig. 2), simulating the standard situation in clinical practice, when some of the measurements are obstructed by the eyelashes or other obstacles. A common procedure in such cases is to discard the elevations corresponding to incomplete rings, which may imply an unnecessary loss of information. One of the advantages of the RBF is that they are not bound intrinsically to circular domains. This fact gives an additional interest to the analysis of the situation modeled by Surface C.

Figure 1.

 

A 3-D representation of surface B.

Figure 2.

 

Elevation data for surface C.

Another important observation is related to the units of measurement. Since the elevations are obtained from synthetic surfaces where the scaling is irrelevant, we chose to fit the data on a unit disc. Hence, the plots appearing in Figures 2 3 4 5 6–7 and the data in Table 1 are given in universal units, whose choice does not affect the results. 

We gathered the discrete elevation data into a vector Z without adding any noise, and solved the overdetermined system (equation 2) in the sense of the LSF. In practice, the collocation matrix M can be very ill-conditioned and numerically rank deficient, so we have to avoid solving the normal equations M^TMa = M^T Z directly. The use of the Moore-Penrose pseudoinverse of M computed by its singular value decomposition (SVD), complemented with regularization, is preferable instead (see e.g., Ref. 25). 

The method can be easily adapted to include the weighted least square fit (WLSF) by left-multiplying equation 2 by a diagonal positive matrix representing the weights. In real-life computation, these weights can reflect the reliability of the data (e.g., portions of the cornea obstructed by eyelashes, poor quality of the tear film; see Ref. 26 for the algorithms that allow separation and identification of the regions of a strong interference). 

All the computations were performed in commercial software (MatLab ver. 7.6, MathWorks Inc., Natick, MA). The vectorization capabilities of MatLab have been extensively exploited, and highly efficient algorithms have been achieved that reduce the computation time drastically (to ∼3 seconds or less), even for the most time-consuming Zernike polynomials fit (compare e.g., with Refs. 1, 3). 

We performed a comparison between families of radial polynomials (Zernike and Bhatia-Wolf) and radial basis functions (Gaussians and inverse multiquadrics). For the reconstruction by polynomials, we normally choose the maximal radial order 6 which corresponds to the first 28 Zernike polynomials, or the first 49 Bhatia-Wolf polynomials, which is the standard in modern aberrometers. ¹¹ The fit with the RBF was performed with a regular grid of centers, like those represented in Figure 3. Observe that to avoid high oscillations on the edge we must use centers situated outside of the cornea, although we omit those located too far from the nodes. 

Figure 3.

 

Dots denote centers of the RBF.

Experiments have been performed also with other functions, such as Sobolev orthogonal polynomials on the disc or multiquadric RBF; however, the results obtained do not differ significantly from those corresponding to other members of the same class, and we decided to omit that discussion for the sake of brevity. 

We have left out of the comparison the Fourier-based techniques for several reasons. First, these methods can be implemented in different ways. If we choose a number of terms v in equation 1 that is smaller than the size of the dataset, then the behavior of the truncated Fourier expansion is very similar to that of the Zernike polynomials (take note that these bases differ only in the radial coordinate). Alternatively, we can take the maximum possible size of v, which endows the Fourier methods with the maximum resolution capacity, but deprives them at the same time of the smoothing ability of the other modal approaches (see Refs. 18–20). Last but not least, the implementation of the Fourier methods is still far from clear and reliable, as a recent discussion shows. ²¹  

Subject A (Fig. 4, top left) is given analytically by  

 

simulating a sphere with a cap removed by a flat cut. The main goal is to check the goodness of detection of the fast variations of the gradient. 

Figure 4.

 

A 3-D representation of surface A. Top left: the original surface. Top right: reconstruction with Zernike polynomials of radial order 6 (v = 28). Bottom left: reconstruction with Bhatia-Wolf polynomials of radial order 6 (v = 49). Bottom right: reconstruction with the inverse multiquadric RBFs (β = −1.5, c = 0.6) with 177 centers.

For subject B we used a sphere with a radial slit (Fig. 1); its level curves are represented in Figure 5, upper left. Its analytic expression is cumbersome, and thus we avoid presenting it here. 

Figure 5.

 

Contour plot for the surface B (top left) and its reconstruction with Zernike polynomials up to order 6 (top right) and 18, v = 190 (bottom left). Bottom right: reconstruction with the inverse multiquadrics with 177 centers, using the parameters (c = 1, β = 5).

The data of subject C were collected from measurements by a corneal topographer (CM02; CSO, Florence, Italy) of the corneal elevations of an actual patient with keratoconus. To approximate the situation to real-life scenarios, we retained the nodes where the elevations were obtained and considered a simulated keratoconus corneal surface modeled with a series of the first 136 Zernike polynomials (radial order 15). This yields an analytic expression for the surface in Figure 6, first row, for which we already know the exact values of the corresponding Zernike coefficients. We represent it only over the domain where the reliable information is available. 

Figure 6.

 

Contour Surface C. Up: original surface. Center: distributions of errors of the model with Zernike (left) and Bhatia-Wolf (right). Down: distributions of errors of the model with Gaussians (left) and multiquadrics (right).

We assess the quality of fit by computing in each case the mean square errors (MSEs). For that purpose, after obtaining vector a in equation 2 we reconstruct the surface by formula in equation 1 and evaluate it at the same nodes where the heights were collected. This gives us the vector of fitted elevations Z̃. Then  

 

Another criterion of the fitting quality is the accuracy in recovery of the Zernike coefficients, especially in the case of incomplete data, which becomes crucial for the discriminant analysis of the decease markers, or for the neural network-based diagnostic software ⁴ such as corneal classification and condition severity-grading utilities. In this sense, for subject C we performed the discrete LSF with Zernike polynomials directly from the raw input data and alternatively fitting the surface previously reconstructed by the Gaussian RBF (see Fig. 7). 

Figure 7.

 

Relative error of the reconstruction of Zernike coefficients for subject C directly by the LSF or fitting previously the surface with Gaussian RBF. The horizontal axis represents the 1-D index j of the Zernike polynomial Z_j , whereas the vertical axis represents the (dimensionless) relative errors.

Results

In this section, we present a comparison of the numerical results obtained with the different methods applied to the three simulated surfaces. 

There are two aspects related to the numerical side of the problem. One is the computational cost, and the other is the sensitivity of the scheme to data perturbations. In the former, the RBF clearly outperform the radial polynomials. The computational time is several times higher for the Zernike polynomials—even for a highly optimized vector algorithm implemented in MatLab—independently if we generate them by their recurrence relation or by the explicit formula. Nevertheless, in our implementation, the reconstruction by the slowest Zernike fit with as many as 231 polynomials (radial order ≤20), took approximately 2.5 seconds (compare with Refs. 1, 3), so that execution time becomes less and less of a concern with the progress of the software optimization and computing power. 

On the other hand, the condition number of the collocation matrices M in equation 2 for the RBF grows very rapidly with the dimension of the problem (number of centers). Despite this undesirable feature, they are still better fit to capture the small local variations in the shape of the cornea than are the standard Zernike polynomials (see e.g., the analysis for subject B in the next section). This problem can be easily addressed though using a Tikhonov-type regularization combined with an SVD computation of the pseudoinverse. On the contrary, the numerical condition of the collocation matrix for Zernike polynomials initially grows slowly with the size (value of v in equation 1), until undersampling sets in, causing an exponential growth of the condition number (phenomenon nicely described in Yoon et al. ²⁰ ; see also Iskander et al. ⁸ ); beyond this point an increase in the number of terms in the expansion (equation 1) becomes counterproductive. 

Let us discuss in more detail the main results in the three cases. 

Surface A: In the portions of the plane where the surface is smooth, all approximation schemes work very well, so that the error is localized in a neighborhood of the ridge formed by the cut. This explains the minimal deviation in the MSE for all methods (see Table 1). However, if our primary goal is the accurate reconstruction of the shape of the cornea, the visual analysis (Fig. 4) shows that the fitting with RBFs outperforms the fitting with radial polynomials. Clearly, the flexibility of the RBF approach, built into the scaling parameters, allows us to capture more easily the rapid variations (or discontinuities) of the gradient of the surface. 

Table 1.

 

Comparison of the MSE Obtained with the Different Methods

Method	Surface A	Surface B	Surface C
Zernike	8.6337e-6	7.0416e-7	2.1226e-4
Bhatia-Wolf	1.9262e-6	6.3971e-7	1.5979e-4
Gaussians	4.6234e-7	2.9088e-7	8.6161e-6
Inverse multiquadrics	4.5003e-7	2.9042e-7	8.9577e-6

Surface B: the purpose was to detect the relatively small details on the surface by fitting the elevation data. As it follows from Figure 5, Zernike polynomials (and in general, radial polynomials) are less well suited for reflecting the small deformation, even if we allow unusually high orders in equation 1. 

Surface C: in this case the deviations of the surface from an “ideal one” are more global, and this fact is immediately reflected in the MSE in Table 1, where we gain at least two orders of precision with the RBF. However, the experiment with the reconstruction of the Zernike coefficients of the surface, either fitting it first by Gaussians or directly by Zernike polynomials, is not conclusive (see Fig. 7). Still, some coefficients are clearly better fit with the former approach, which can be significant for an early detection of keratoconus (see Ref. 27). 

These results can be improved further by using a combined approach: fitting the raw data with Zernike polynomials with a very low order (n = 1 or n = 2), subtracting the fit from the original elevations, and approximating the new data by LSF with the RBF. 

Discussion

The first important observation concerns the number of terms v in the modal reconstruction (equation 1). When we use radial polynomials (Zernike, Bhatia-Wolf), the number of terms corresponds to the maximum order of aberrations or frequencies that can be captured or represented by the right-hand side in equation 1. The computational complexity of the basis functions f_j in equation 1 grows with the index j. 

Experiments show a saturation phenomenon: although for a low v, an addition of a new term renders a substantial improvement in the goodness of the fit, higher orders have less and less impact. Moreover, a small change in a localized subset of data may imply a substantial modification of all entries of the coefficient vector a. On the contrary, the number of terms v used in the fitting with the RBF is given by the numbers of centers Q_j . Higher values of v imply in this case more flexibility in the approximating scheme. The localization property of the RBF used implies also that a small local variation in the data has only a “local” impact on the coefficients of a. Basis functions f_j in equation 1 for different values of the index j are computationally identical, simply “aimed” at different points of the disc. 

Hence, the amount v of terms used for approximation with radial polynomials or with RBF should be compared with care. 

As it was observed previously in Iskander et al., ¹³ Bhatia-Wolf polynomials achieve higher precision in surface approximation than their classic Zernike counterparts. Nevertheless, a clear conclusion of this research is that the Zernike polynomials still work perfectly well as a reconstruction method of a nonseverely aberrated surface with a small SRI. They also are an appropriate tool for recovering the lower Zernike coefficients. 

However, these coefficients fail to capture small deformations of the anterior surface of the cornea. In particular, if such deformations turn out to be markers of an eye disease, it is reasonable to complement the Zernike coefficients with additional input parameters for the neural network–based diagnostic software (see the pioneering work, ⁴ in which corneal cases with no surface singularities were considered). When severe curvature changes are present, the accuracy of the fit (taking into account the small features of the surface) can become a priority, since it allows extracting reliably other shape indices of the approximated surface. In such a situation, the flexibility of the RBF functions, combining some properties of a zonal reconstruction (localization) with the simplicity of a modal scheme, can become relevant. 

Thus, a combined approach seems promising: using Zernike or Bhatia-Wolf polynomials of a low degree to obtain the fundamental part of the shape of the cornea, with a subsequent refinement by RBF. 

However, additional research is needed to address some computational concerns such as an automatic selection of the scaling parameters of the RBF, or better control of the condition numbers of the corresponding collocation matrices. These aspects will be subject of a further investigation. 

Footnotes

 Supported in part by Junta de Andalucía and Grants FQM-229 and P06-FQM-01738 (AM-F, AMD); Ministry of Science and Innovation of Spain Project MTM2008-06689-C02-01 (AM-F, AZ); Spanish Ministry of Education and Science Project MTM2006-13000-C03-02 (AMD); Junta de Andalucía Grant FQM-481 (AM-F); and Carlos III Health Institute of the Ministry of Science and Innovation of Spain Grant PI08/90519 (all authors).

Footnotes

 Disclosure: A . Martinez-Finkelshtein, None; A. M. Delgado, None; G. M. Castro, None; A. Zarz o, None; J. L. Alió, None

Footnotes

 The publication costs of this article were defrayed in part by page charge payment. This article must therefore be marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

The authors thank the anonymous referees, whose helpful comments clearly improved the text.