purpose. Zernike polynomials have been successfully used for approximately 70 years in many different fields of optics. Nevertheless, there are some recent discussions regarding the precision and accuracy of these polynomials when applied to surfaces such as the human cornea. The main objective of this work was to investigate the absolute accuracy of Zernike polynomials of different orders when fitting several types of theoretical corneal and wave-front surface data.

methods. A set of synthetic surfaces resembling several common corneal anomalies was sampled by using cylindrical coordinates to simulate the height output files of commercial videokeratography systems. The same surfaces were used to compute the optical path difference (wave-front [WF] error), by using a simple ray-tracing procedure. Corneal surface and WF error was fit by using a least-squares algorithm and Zernike polynomials of different orders, varying from 1 to 36 OSA-VSIA convention terms.

results. The root mean square error (RMSE) ranged—from the most symmetric corneal surface (spherical shape) through the most complex shape (after radial keratotomy [RK]) for both the optical path difference and the surface elevation for 1 through 36 Zernike terms—from 421.4 to 0.8 μm and 421.4 to 8.2 μm, respectively. The mean RMSE for the maximum Zernike terms for both surfaces was 4.5 μm.

conclusions. These results suggest that, for surfaces such as that present after RK, in keratoconus, or after keratoplasty, even more than 36 terms may be necessary to obtain minimum accuracy requirements. The author suggests that the number of Zernike polynomials should not be a global fixed conventional or generally accepted value but rather a number based on specific surface properties and desired accuracy.

^{ 1 }Zernike was a Dutch physicist who worked in several fields of optics and won the Nobel Prize for inventing the phase-contrast microscope in 1935.

^{ 2 }This instrument is still used today to study biological specimens without the need for dyes. Dyes can help in the visualization process by emphasizing contrast but usually spoil the sample.

^{ 3 }Although they have only been applied more recently to the description of the optical aberrations of the human eye,

^{ 4 }

^{ 5 }they have become a standard in this field also.

^{ 6 }Nevertheless, there have been some recent discussions (Smolek MK, et al.

*IOVS*1997;38:ARVO Abstract 4298; Smolek MK, et al.

*IOVS*2002;43:ARVO E-Abstract 3943)

^{ 7 }

^{ 8 }

^{ 9 }

^{ 10 }among colleagues in the eye care community regarding the accuracy and even the usefulness of these polynomials, specifically for application in the visual sciences. There are arguments that ZPs are not sufficient to represent visually significant aberrations,

^{ 7 }

^{ 8 }and other investigators

^{ 9 }

^{ 10 }have stated that ZPs should be used carefully when fitting complex surfaces.

^{ 11 }on corneal elevation and/or wave-front (WF) data from in vivo eyes, but also to implement this method on theoretical surfaces, synthetically generated by computer algorithms, which eliminates the problem of videokeratography measurement errors

^{ 12 }

^{ 13 }and even avoids the limitations of the molding process of test surfaces on lathes, which are usually limited to spherical, ellipsoidal, or parabolic surfaces. In the analysis conducted in this study I applied ZPs in many different ways and to different theoretical surfaces, avoiding problems such as misalignment, data noise from image processing, and the other just-mentioned technical limitations. This may give an insight into the intrinsic relations between actual surface irregularities and accuracy when fitting them with ZPs.

^{ 6 }convention, which has recently also become an ANSI standard,

^{ 14 }determines, among other aspects, the nomenclature that should be used to refer to each Zernike coefficient, to avoid confusion when comparing results from different instrumentation, but does not suggest an upper limit to the number of terms. Nevertheless, many researchers and manufactures throughout the eye care community are currently using, at most, the first 36 terms, regardless of the application. The computations shown herein suggest that the ideal number of Zernike terms should not be a fixed standard for any general surface, but rather one that would allow for a minimum standard error for specific types of surfaces that are being represented by ZPs.

*wave aberration*(sometimes also referred to

*wave-aberration function*or

*wave-front error*—which will be abbreviated herein simply as

*W*).

^{ 15 }

^{ 16 }(although we do not know whether nature did this also for economic reasons!). Fortunately, today there are several very sophisticated optical design software programs that make this task much less empiric than it was years ago.

^{ 3 }

*R.*For each point on the exit pupil, the optical path difference (OPD), represented by the wave-aberration function (

*W*), is measured between the spherical reference surface and the aberrated WF along the radius of the reference surface. The wave aberration is usually described as a function of Cartesian or Polar coordinates (

*W*(

*x*,

*y*) or

*W*(ρ,θ)) obtained over the exit pupil, which is now used to describe the wave aberration in a continuous format.

*W*. The usual procedure is to represent it as a polynomial expansion. This is very useful because each term of the expansion may itself represent a specific type of aberration and also may determine how much of it is present on the entire WF. There are two sets of such polynomials that traditionally have been used for the description of the aberration function. In optical design the

*Seidel polynomials*are typically used. In optical testing and measurement,

*W*usually has to be deciphered (and not chosen or minimized, as happens in optical design), and the procedure that became common is to fit

*W*with a set of ZPs. I will briefly describe the SPs and some of their terms, to show why they are different from ZPs and why ZPs are a better solution for more complex applications.

*W*can be represented mathematically by a set of SPs by

*C*

_{ ijk }is the WF aberration coefficient,

*H*is the fractional image/object height for the chief ray, and the other terms are the radial order and polar frequency. The use of normalized pupil coordinates is a matter of convenience, and dimensionality is maintained by the

*C*

_{ ijk }coefficient.

*i + j*= 4 in the expansion given by equation 1 . The most familiar aberrations are spherical aberration, coma, astigmatism, field of curvature, and distortion. Figure 3shows two- and three-dimensional representations of some of the Seidel aberrations.

^{ 17 }When decentralization and tilt are present, there is no term in the SP expansion that can account for the wave aberration induced, because SPs are based on rotationally symmetric terms. When this problem arises, ZPs make a difference. As will be shown in the next section, they have all the invariance, normalization, and other interesting properties of the Seidel expansion, but they can also account for a more complete set of possible optical imperfections.

^{ 1 }

^{ 5 }

^{ 12 }

^{ 13 }

^{ 1 }This guarantees that one may fit any piece-wise continuous surface in space, given that a sufficiently large number of terms are used. This is certainly the case for most corneas

^{ 18 }that have not undergone physical trauma, early post-keratoplasty, in keratitis, or in any other severe disease that may cause abrupt changes in curvature.

*f*(

*x*) as an infinite summation of sines and cosines with different coefficients and arguments. Its general expression is given by

*f*(

*x*).) For signal type 2, the story is different. Because the signal is much more complex (has several frequencies associated with it) the truncation of the FS at different frequencies has significant influence on accuracy.

*f*(

*x*) with an FS we calculated the root mean square error (RMSE) associated with each function, given by

*IF*

_{ n }represents the FS up to the

*nth*order. Results for a successive number of terms for the square wave and sinusoidal functions are provided in Table 1 .

*f*(

*x*) and

*g*(

*x*), are said to be orthogonal over the interval

*a*≤

*x*≤

*b*with weighting function

*w*(

*x*) if

*n*and

*m*represent the indexes of each polynomial,

*C*

_{ n }is a constant and, when it assumes a value of 1, the polynomials are also normalized, and δ represents the Kronecker delta

^{ 19 }and assumes a value of 0 if

*m*=

*n*and 1 if

*m*≠

*n*. The general properties given by equations (7 8 9 -10)are also applicable to functions or polynomials defined in larger domains (such as the

*x–y*plane). The only difference is that a double integral should be implemented. The orthonormal functions are also said to be “complete” in the closed interval

*x*∈ (

*a*,

*b*) if, for every piece-wise continuous function

*f*(

*x*) in this interval, the squared error

*n*→∞. Symbolically, a set of functions is complete if

*x*in the considered interval. If this same concept is applied to a set of functions or polynomials for which the domain now is the

*x*,

*y*plane, the symbolic representation becomes

*x*-axis (

*a*,

*b*) and the

*y*-axis (

*c*,

*d*). Also, if these polynomials are in cylindrical coordinates—that is, φ(ρ,θ)—then equation 13can be written in terms of these new coordinates, and the limiting interval will also be determined by (ρ,θ). More specifically, it is useful to write these polynomials in cylindrical coordinates when the domain has polar symmetry, which often happens in the case of optical apertures and also the human pupil. This is why the VSIA standards for ZPs are given in cylindrical coordinates. ZPs also obey all the properties discussed thus far—a fact that is not proven herein, given that there are very good references that demonstrate these properties and also for the sake of brevity. The reader is directed to Born

^{ 1 }for a thorough discussion and demonstration of these properties and also an explanation of the recursive formulas for generating ZPs of any order. In this way, it can also be affirmed that ZPs form a “complete set of orthonormal polynomials” inside the unit circle domain (0 − 1, 0 − 2π).

*r*

_{ a }is the apical radius,

*p*= 1 −

*e*

^{2}is the “shape factor,” and

*e*is the eccentricity. When

*e*is set to 0, the shape factor is 1, and equation 14becomes a sphere with radius

*r*

_{a}. Another typical parameter is the “asphericity” (

*Q*) of the surface, which is equal to −

*e*

^{2}, so that the shape factor may also be written as

*p*= 1 +

*Q*. To model an astigmatic cornea, an ellipsotoric surface

^{ 20 }was used, where the apical radius is a function of the polar angle (θ)

*r*

_{v}and

*r*

_{h}are the vertical and horizontal apical radii, respectively.

*a*is keratoconus intensity. This surface has continuous curvature (proportional to the second derivative), an important feature to mimic in vivo corneal surfaces.

^{ 18 }Although decentralized keratoconus is also common in in vivo eyes, this surface mimics a centralized keratoconus—that is, one for which the apex coincides with the intersection of the optic axis with the anterior cornea. This simplifies the computational efforts of synthetic Placido image simulation. There is also the possibility of tilting this surface by using rotation and translation transformation matrices commonly used in computer graphics,

^{ 21 }

^{ 22 }to obtain decentralized keratoconus. This was not implemented in the present work for the sake of brevity, although it is certainly an important analysis to undertake in further work.

^{ 18 }a generalized version of a surface proposed by Rand et al.,

^{ 23 }was used for the post-RK cornea

*z*) is now a function of the axial distance (ρ) and the polar angle (θ).

*R*is a parametric factor that depends on the axial distance, so that

*a*here is the “wound intensity,” proportional to the depth of the scalped incisions. This surface is a sphere of radius

*r*

_{0}with an added sinusoidal corrugation. Table 2is a list of surfaces that were used in the simulations and their respective parameters, according to the description just provided.

*z*is the ZP fit of the corneal height,

*z*

_{c}is the theoretical surface height computed from equations 12 13 15 16 and 17 , and the ρ values are presented in their parameterized form—that is, as a function of polar angle (θ) and Placido disc (

*n*). An analogous procedure was implemented to compute the errors associated with the WF error fit.

*n*and

*n*′). An object point at position (P) localized at the object plane (O) has its image formed at point (P′) on the image plane (I). A marginal ray intersects the exit pupil at some point (ρ

_{c},θ

_{c}) on the corneal surface.

*W*along the marginal ray is calculated as the difference in optical path length from the chief ray. For a single ray

*s*) is chosen to be at infinity (>6 m), and the image distance

*s*′ is calculated by approximating the cornea to a lens, with radius (

*r*

_{m}) calculated from the mean value of all axial radii of curvature of the surface

*r*

_{a}is the apical radius (computed in the corneal elevation simulation phase). Alternatively, but with little loss in precision and generalization, the paraxial approximation may be applied and the

*r*

_{m}is then replaced by the apical radius of curvature. When the surface is astigmatic, the mean value of the vertical plus horizontal radii may be used [(

*r*

_{v}+

*r*

_{h})/2].

*s*′ is determined, the distances of the marginal and chief rays can be computed, and equation 20can be applied over the entire Placido image domain

*W*,

*z*

_{c},

*x*

_{c}, and

*y*

_{ c }are all parameterized functions of (ρ(θ,

*n*), θ). From this equation, the optical aberration (ρ,θ) for each corneal surface point

*z*

_{ c }(ρ,θ) obtained from the theoretical surface can be calculated.

^{ 12 }). Because of this implementation detail, all Placido image domains and exit pupil domains of all the synthetic images were limited to a total radial distance of 4 mm when computing the corneal elevation, optical aberrations, and errors. This radial distance limit guarantees that all Zernike fit errors for different surfaces are compared for domains having the same area. Aberrations for smaller pupils were not computed, because an 8-mm pupil maximizes aberration and at the same time serves as a reasonable approximation of maximum in vivo pupil sizes.

^{ 24 }

*n*,θ),

*C*

_{ j }are the Zernike coefficients, and

*Z*

_{ j }are the ZPs.

*N*data points. This procedure consists of minimizing the sum

*dS/dC*

_{ t }= 0 for

*t*= 1,…,

*k*, must be found, where

*k*is the total number of coefficients, so that

*AC = b*with 36 equations and 36 unknown values of

*C*. By solving this linear system through conventional procedures, such as the

*Gaussian elimination method*, the 36 Zernike coefficients are found for each surface.

^{ 18 }has obtained accuracy of up to 0.1 μm for elevation of synthetic surfaces with an algorithm that avoids the skew ray problem, and Guirao and Artal

^{ 13 }have obtained accuracy of 0.2 μm in practical measurements on test surfaces, using a commercial videokeratograph and synthetic ellipsoidal surfaces.

^{ 25 }On the contrary, for commercial diagnostic instrumentation used by eye care professionals, processing time should be a major factor. With most commercial videokeratography instruments available today, data processing takes only a few seconds,

^{ 26 }

^{ 27 }

^{ 28 }

^{ 29 }which means that algorithms that take on the order of minutes would be unacceptable for professional clinical use. I used commercial programming language (MatLab; The MathWorks, Natick, MA), input files which have the same number of data points as the Eyesys videokeratograph (5760), and three different IBM-compatible computers with processors and RAM, respectively, of 1.6 GHz with 1 GB, 1.7 GHz with 1 GB, and 1.7 GHz with 0.5 GB. Processing times for each number of Zernike coefficients are illustrated in Figure 11 .

*t*) as a function of number of Zernike coefficients (

*n*)

^{ 27 }

^{ 25 }

^{ 26 }

^{ 27 }

^{ 28 }

^{ 29 }).

**Figure 1.**

**Figure 1.**

**Figure 2.**

**Figure 2.**

**Figure 3.**

**Figure 3.**

**Figure 4.**

**Figure 4.**

n | Sine Function | Square-Wave Function |
---|---|---|

1 | 0.5671 | 0.8800 |

2 | 0.3108 | 0.5373 |

3 | 0.2962 | 0.5035 |

4 | 0.2965 | 0.4962 |

5 | 0.2854 | 0.4400 |

6 | 0.2854 | 0.4400 |

7 | 0.2817 | 0.3838 |

8 | 0.2897 | 0.3765 |

9 | 0.2849 | 0.3427 |

**Figure 5.**

**Figure 5.**

**Figure 6.**

**Figure 6.**

Surface | Description | Equation | Parameters |
---|---|---|---|

A | Sphere | 14 | p = 1, r _{a} = 7.80 |

B | Ellipsoid | 14 | p = 0.50, r _{a} = 7.80 |

C | Ellipsotoric | 14 15 | p = 0.75, r _{h} = 7.50, r _{v} = 8.00 |

D | Ellipsotoric | 14 15 | p = 0.50, r _{h} = 8.00, r _{v} = 7.50 |

E | Ellipsotoric | 14 15 | p = 0.30, r _{h} = 9.00, r _{v} = 7.00 |

F | Keratoconic | 16 | a = 0.0125, ρ_{1} = 1.00, ρ_{2} = 3.00, r _{0} = 7.00 |

G | Keratoconic | 16 | a = 0.0250, ρ_{1} = 1.50, ρ_{2} = 3.00, r _{0} = 8.00 |

H | Keratoconic | 16 | a = 0.0500, ρ_{1} = 2.00, ρ_{2} = 3.00, r _{0} = 9.00 |

I | Post-RK | 17 18 | a = 0.0250, ρ_{1} = 2.00, ρ_{2} = 4.00, r _{0} = 7.00 |

J | Post-RK | 17 18 | a = 0.0500, ρ_{1} = 2.00, ρ_{2} = 5.00, r _{0} = 8.00 |

**Figure 7.**

**Figure 7.**

**Figure 8.**

**Figure 8.**

**Figure 9.**

**Figure 9.**

**Figure 10.**

**Figure 10.**

**Figure 11.**

**Figure 11.**