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

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

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

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

^{ 6 }

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

^{ 8 }It is known

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

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

*C*(

*P*) is the corneal elevation at the point

*P*of the plane,

*f*is the a priori chosen basis function, and

_{j}**a**= (

*a*

_{1},…,

*a*is the expansion coefficient (the superindex

_{v})^{T}*T*means matrix transposition). In this setting, fitting the equation to a discrete set of elevation data

**Z**= (

*Z*

_{1},…,

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

^{T}M**a***M*, which is unique if the collocation matrix

^{T}**Z***M*is of maximum rank

*v*.

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

_{j}*Z*(corresponding to the radial non-negative integer

_{n}^{m}*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.

^{ 10 }

*n*,

*m*) into a 1-D array

*j*. The most widely acceptable one is Choosing in equation 1 as

*f*the Zernike polynomials of radial order ≤

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

^{ 11 }corresponding to values

*v*= 15 and

*v*= 28 for

*n*= 4 and

*n*= 6, respectively. It was shown in Iskander et al.

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

^{ 12 }

*B*, whose fitting properties have been analyzed in Iskander et al.

_{n}^{m}^{ 13 }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)

^{2}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*are not algebraic polynomials in the Cartesian variables

_{n}^{m}*x*and

*y*, but they expand in series of monomials

*x*, ρ = (

^{i}y^{j}ρ^{k}*x*

^{2}=

*y*

^{2})

^{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*by

_{j}*j*=

*n*(

*n*= 1) =

*m.*

^{ 14 }and the Sobolev orthogonal polynomials on the disc.

^{ 15,16 }

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

*f*in equation 1 sets of RBFs, defined in their simplest form by translates of a given function Φ: where points

_{j}*Q*, called centers of the RBFs, are conveniently chosen, and ‖·‖ denotes the Euclidean distance on the plane. The general theory of

_{ j }*interpolation*by RBFs is developing rapidly, and several criteria for Φ can be found in the literature.

^{ 22 }In particular, standard options are the so-called Gaussians and inverse multiquadrics, corresponding, respectively, to Φ(

*t*) = exp(−

*at*) and Φ(

*t*) = (

*t*=

*c*

^{2})

^{−β}, with positive parameters

*a*,

*c*, and β. However, we are unaware of any deep theoretical analysis of the LSF with RBF. This, according to Buhmann,

^{ 23 }is a highly relevant and interesting field of research.

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

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

_{j}^{ 24 }

*Q*are fixed a priori, the values of with

_{j}*f*given by equation 3, can be computed and stored so that the Zernike coefficients are easily recovered by the scalar product of two vectors.

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

*P*, s = 1,…,6144, with polar coordinates (ρ

_{s}*, θ*

_{i}*) where ρ*

_{j}*=*

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

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

^{T}M**a***M*directly. The use of the Moore-Penrose pseudoinverse of

^{T}**Z***M*computed by its singular value decomposition (SVD), complemented with regularization, is preferable instead (see e.g., Ref. 25).

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

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

^{ 21 }

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

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

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

^{ 20 }; see also Iskander et al.

^{ 8 }); beyond this point an increase in the number of terms in the expansion (equation 1) becomes counterproductive.

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 |

*n*= 1 or

*n*= 2), subtracting the fit from the original elevations, and approximating the new data by LSF with the RBF.

*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*in equation 1 grows with the index

_{j}*j*.

*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*. Higher values of

_{j}*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*in equation 1 for different values of the index

_{j}*j*are computationally identical, simply “aimed” at different points of the disc.

*v*of terms used for approximation with radial polynomials or with RBF should be compared with care.

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

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

*J Refract Surg*. 2004;20:S537–S541. [PubMed]

*Invest Ophthalmol Vis Sci*. 2003;44(11):4676–4681. [CrossRef] [PubMed]

*Invest Ophthalmol Vis Sci*. 2005;46:1915–1926. [CrossRef] [PubMed]

*Optom Vis Sci*. 2005;82:151–158. [CrossRef] [PubMed]

*Optom Vis Sci*. 1989;73:721–728. [CrossRef]

*J Opt Soc Am A*. 2002;19:137–143. [CrossRef]

*J Cataract Refract Surg*. 2005;31:2350–2355. [CrossRef] [PubMed]

*IEEE Trans Biomed Eng*. 2001;48(1):87–95. [CrossRef] [PubMed]

*Appl Opt*. 2006;45:6945–6964. [CrossRef]

*Vision Science and its Applications*,

*OSA Technical Digest.*Washington, DC: Optical Society of America; 2000;paper SuC1.

*J Cataract Refract Surg*. 2005;31:1114–1127. [CrossRef] [PubMed]

*Proc Cambridge Philos Soc*. 1954;50:40–53. [CrossRef]

*IEEE Trans Biomed Eng*. 2002;49(4):320–328. [CrossRef] [PubMed]

*SIAM J Appl Math*. 1966;14(3):476–489. [CrossRef]

*J Approx Theory*. 2006;138:232–241. [CrossRef]

*J Approx Theory*. 2008;152:52–65. [CrossRef]

*Appl Opt*. 1991;30:501–503. [CrossRef]

*J Refract Surg*. 2006;22(11):943–948. [PubMed]

*J Cataract Refract Surg*. 2007;33:999–1004. [CrossRef] [PubMed]

*J Refractive Surg*. 2008;24:582–590.

*J Refract Surg*. 2009;25:9–11. [PubMed]

*Cambridge Monographs on Applied and Computational Mathematics*. Cambridge, UK: Cambridge University Press; 2005.

*Radial Basis Functions: Theory and Implementations*. Cambridge, UK: Cambridge University Press; 2003.

*Opt Commun*. 1999;163:259–269. [CrossRef]

*Numerical Linear Algebra*. Philadelphia: Society for Industrial and Applied Mathematics (SIAM); 1997.

*IEEE Trans Biomed Eng*. 2008;55(10):2381–2387. [CrossRef] [PubMed]

*J Refract Surg*. 2006;22:539–545. [PubMed]