1.

Introduction

In the last few years, a large number of accurate devices for measuring corneal morphology have appeared. Among them, those based on Placido rings and Scheimpflug cameras are the most widely used.^{1, 2, 3, 4, 5, 6, 7} These devices can provide height data of the corneal surfaces that can be fitted to analytical surfaces through polynomial expressions (modal approach).^{8, 9, 10} These can be later used for aberration calculation and ray trace solving.^{11, 12, 13}

Zernike polynomials are often used as adjusting polynomials, but they have shown to be not precise when describing highly irregular corneas,^{14, 15, 16} since they do not fit properly to irregularities. In 2003, Smolek and Klyce¹⁴ concluded that Zernike polynomials did not fully characterize the surface features that affect vision. Later, Klyce, Karon, and Smolek¹⁵ found that using Zernike polynomials in normal eyes was acceptable, but for eyes after corneal surgery or eyes with corneal pathology such as keratoconus, the Zernike method fails to capture all clinically significant higher-order aberrations. Furthermore, it is very difficult to assess a priori how many terms are necessary to achieve an acceptable accuracy in the Zernike reconstruction for any given corneal shape.¹⁷

Common measuring devices are based in imaging techniques. An image reflected from a pattern—linear or annular—is projected onto a charge-coupled device (CCD) and is digitally processed.¹⁸ From this information, a height map is obtained. Surfaces with radial structure do not fit well with squared detectors or matrices,¹⁹ since density of samples in the azimuthal coordinate depends on the radial coordinate, due to the fact that the number of different angular samples in the center of the matrix is lower than in the periphery. Unfortunately, acquiring devices are based on squared CCD, and little can be done to change this fact. Reconstruction methods may consider such limitations and obtain the best possible results. Traditional reconstruction of the optical surface with Zernike polynomials²⁰ does not pay attention to the particular sampling distribution of the analyzed surfaces, and thus it does not provide optimum results in the central area around the apex, which is the most important zone under the optical point of view. One of the main limitations of the polynomial representation is that each term extends its influence over the entire pupil.

The goodness of fit of a surface model strongly depends on the number of samples and its distribution. Let us consider a spherical-like surface given by height samples in Cartesian coordinates. The number of points that contribute to the polynomial fit is proportional to the area of the considered surface. Thus the ratio between the number of samples included in an external annular region defined between radii $R_1$ and $R_2$ and the central part of aperture radius $R_1$ is quadratic, i.e., ${(R_2∕R_1)}^2-1$ . The effect is not very important for a spherical surface, but if the surface is affected by any irregularity in the periphery, it determines the coefficients in the central optical area, although the real effect of the irregularity on the retinal image is limited. Therefore, the fitting error is distributed over the entire surface, not just at those zones where the deformation is located.

B-spline polynomials are especially well suited for fitting complex-shaped surfaces.^{21, 22} They have the advantages of being locally defined and having great flexibility, which allows control of their smoothness and polynomial degree.^{23, 24} In Ref. 25, the authors propose an alternative technique for reconstructing the corneal shape from elevation data using sets of radial basis functions (RBFs) in a modal approach. Due to the fast decay of RBFs, the modal approach exhibits features of the zonal approach, eventually capturing small deformations of the surface, which are missed by the polynomial fitting. As Ref. 25 states, we think that there is no unique or best approach to corneal surfaces, and a combination of techniques can be a good strategy.

Light distributions propagated from the cornea can be numerically evaluated by means of a Fresnel propagation algorithm.²⁶ Numerical calculation of convergent Fresnel patterns through fast Fourier transform usually requires a large number of samples $(4\times {10}^6)$ to fulfil the Nyquist condition around the focus.²⁷ In Ref. 28, the authors present a simple method that permits subsampling the Fresnel pattern while maintaining the Nyquist condition, and thus preventing the appearance of aliasing effects in the calculation. Later, in Ref. 29, the authors proposed modifying the initial wavefront and relaxing the Nyquist condition, thus giving a more efficient numerical algorithm. Corneal height data sampling provided by corneal topographers is still lower than the one imposed by the Nyquist condition. Hence, we need a method that allows an accurate resample of corneal height data.

In this work, we propose a zonal Zernike fitting (combination of zonal and modal approaches) of corneal height data. It consists of obtaining Zernike coefficients over local areas. A local mask is applied on the surface, and polynomial fit is done for this part of the surface. The coefficients so obtained have no physical meaning due to the fact that Zernike polynomials are usually defined over a unit disk. However, we used the method as a reconstruction tool in a way similar to Lundström, Unsbo, and Gustafsson.³⁰ The mask is then displaced over the surface and the process is repeated on each local area. The surface reconstruction takes into account the fittings in all zones, diminishing the influence of smooth areas over irregular ones and vice versa. The method limits the influence of the peripheral irregularities over the central corneal area calculation, thus giving accurate reconstruction of the central optical zone. This fact will be of special interest in the evaluation of wavefront aberrations of irregular corneas (e.g., keratoconus, pellucid marginal corneal degeneration, corneal etc.). Thus, we join the advantage of a simpler analytical form of the modal representation together with the benefit of local definition of the zonal one. We have compared zonal Zernike fitting with the modal approach, evaluating both techniques over two different surfaces: a theoretical irregular surface and height data of a real keratoconic cornea.

2.

Methods

Zernike polynomials are often used as an expansion of corneal height data and to analyze of optical wavefronts, being low order terms directly related to classical Seidel aberrations.³¹ As a complete modal set, any surface can be approximated by a linear combination of circular polynomials as follows:³²

The modal fitting principle of the Zernike representation is shown in Eq. 1. The data $W(\rho ,\theta )$ are approximated by a polynomial function extended over the whole domain (Fig. 1 ). The best estimation of $c_j$ parameters is obtained by solving the linear least-squares problem described by the system of equations that can be deduced from Eq. 1. Our basic equation is $W=\mathbf{Z}c$ , which is a linear transformation on $c$ , where $c$ are the expansions coefficients, $W$ is a discrete set of elevation data, and $\mathbf{Z}$ is a matrix of discrete Zernike polynomials. The problem is choosing the coefficients $c$ to minimize, in some sense, the difference between the observed $W$ discrete elevations and the prescribed output of the system $\mathbf{Z}c$ . The function to be minimized is:^{33, 34}

Fig. 1

Area of the height data where modal fitting is obtained.

Zonal Zernike fitting is developed in the same way as a modal one, but the equation $W=\mathbf{Z}c$ is not solved over the whole domain. The surface is divided into zones that overlap each other (see Video 1 ). It takes advantage of increasing the influence of the central zone in the adjustment, thus providing a better description of the more interesting optical area. On the contrary, in the modal fitting, periphery data have a heavier influence than central ones due to the fact that they are more numerous.

Video 1

Areas of the height data where zonal Zernike fitting is obtained. Note the overlap of masks (QuickTime, 823.5 KB).

Let us call $\mathbf{W}(u,v)$ a surface described by an $\mathbf{N}\times \mathbf{N}$ matrix, with $u$ and $v$ being discrete indexes from 1 to $\mathbf{N}$ . A local area of $\mathbf{M}\times \mathbf{M}$ pixels, being $\mathbf{M}<\mathbf{N}$ , is selected in this matrix where we solve Eq. 1. Zernike decomposition for this region is obtained, thus having a local description of the surface. The calculation area is displaced along the entire surface, and a zonal Zernike polynomial fit is obtained each time. Reconstruction on one point is done by evaluating all Zernike decompositions, describing the local areas that overlap on this point and taking the mean value of the obtained height. Calculation windows for selecting local regions are implemented by constructing auxiliary matrices of size $\mathbf{N}\times \mathbf{N}$ , whose elements are equal to one in the region of interest and zero elsewhere. These matrices are described by:

As we said before, one point of the surface may belong to different local regions. Thus, it is necessary to take into account how many regions overlap on a single pixel. According to our previous description, the reconstructed surface can be written as

3.

Results

Obtained $l\mathrm{RMSD}$ for a spherical surface of curvature radius of $6.5\mathrm{mm}$ plus a scaled Franke’s function defined on an area with a diameter of $4\mathrm{mm}$ sampled every $0.1\mathrm{mm}$ are shown in Fig. 2 . We have applied modal and zonal Zernike fittings with two different masks of size $21\times 21\mathrm{px}$ and $11\times 11\mathrm{px}$ , this last one being optimum according to Eq. 9. The best fit is obtained for the maximum $l\mathrm{RMSD}$ . We find different maxima for each fitting: zonal Zernike fitting with a $11\times 11\mathrm{px}$ mask at 43 polynomials, 96 polynomials for a $21\times 21\mathrm{px}$ mask, and 114 polynomials for modal fitting. Notice that zonal Zernike fitting provides a better adjustment than the modal one for any number of considered polynomials.

Fig. 2

Reconstruction error between the Franke’s plus sphere surface and the reconstructed surfaces. $l\mathrm{RMSD}$ for modal and zonal Zernike fitting with a $21\times 21$ and $11\times 11\mathrm{px}$ mask.

Figures 3, 4, 5 present the differences between the original surface and the reconstructed ones using the polynomial decomposition that provides the best adjustment for each method. All differences are set to the same scale to emphasize the influence of zonal Zernike fitting. Notice that for both zonal Zernike fitting cases, the reconstruction error is lower than for the global case. Moreover, as the local mask becomes smaller, the central area is better reconstructed.

Fig. 3

Difference between the Franke’s plus sphere surface and the reconstructed surface obtained using the zonal Zernike fitting that provides the best adjustment for $11\times 11\mathrm{px}$ (43 polynomials).

Fig. 4

Difference between the Franke’s plus sphere surface and the reconstructed surface obtained using the zonal Zernike fitting that provides the best adjustment for $21\times 21\mathrm{px}$ (96 polynomials).

Fig. 5

Difference between the Franke’s plus sphere surface and the reconstructed surface obtained using the modal fitting that provides the best adjustment (120 polynomials).

Obtained $l\mathrm{RMSDs}$ for the height data of a real keratoconic cornea with a pupil diameter of $4\mathrm{mm}$ have been calculated and are shown in Fig. 6 . We have applied modal and zonal Zernike fittings using the same masks as before. Again, it is shown that zonal Zernike fitting provides a better adjustment than the global one, as we saw in the theoretical surface case.

Fig. 6

Reconstruction error between the height data of a real keratoconic cornea and the reconstructed surfaces. $l\mathrm{RMSD}$ for modal and zonal Zernike fitting with a $21\times 21$ and $11\times 11\mathrm{px}$ mask.

Contrary to what one may think, increasing the number of terms in the Zernike polynomial decompositions does not gives a more realistic approximation in numerical image analysis. In practice, the SVD method can always provide a unique solution, although the matrix ${\mathbf{Z}}^{+}$ is ill-conditioned and numerically rank deficient. This rank deficiency may be diagnosed through the condition number. The problem arises when the condition number is too large, because the SVD provides a bad least-squares solution.²² This fact can be seen in the falls of the $l\mathrm{RMSD}$ plotted in Fig. 6. There is a number of polynomials from which the solution obtained is degraded.

Figures 7, 8, 9 present the height differences between real keratoconic corneal data obtained from the Pentacam and the reconstructed ones using the polynomial decomposition that provides the best adjustment for each method. All differences are set to the same scale to emphasize the influence of zonal Zernike fitting. Results are quite similar to those obtained for the simulation. Again, results show that the keratoconic corneal height description through Zernike polynomials over the whole surface produces fitting errors in the whole area, whereas the zonal Zernike analysis proposed here provides a lower surface fitting error, mainly in the central region (see Fig. 7), which is optically the most important.

Fig. 7

Difference between the height data of a real keratoconic cornea and the reconstructed surface obtained using the zonal Zernike fitting that provides the best adjustment for $11\times 11\mathrm{px}$ (217 polynomials).

Fig. 8

Difference between the height data of a real keratoconic cornea and the reconstructed surface obtained using the zonal Zernike fitting that provides the best adjustment for $21\times 21\mathrm{px}$ (196 polynomials).

Fig. 9

Difference between the height data of a real keratoconic cornea and the reconstructed surface obtained using the modal fitting that provides the best adjustment (178 polynomials).

4.

Conclusion

In this work, we develop a zonal Zernike fitting of corneal height data. It consists of obtaining Zernike coefficients over local areas. A local displacing mask is applied on the surface, and polynomial fit is done for this part of the surface. The surface reconstruction takes into account the fittings in all the zones, thus diminishing the influence of smooth areas over irregular ones and vice versa. Zonal Zernike fitting has the same advantage of the simple analytical form of the modal representation, together with the benefit of local definition of the zonal one. If we compare for the same polynomial order (28 coefficients) our results with those presented in Ref. 24, they achieve a RMSD around $3\cdot {10}^{-3}\mathrm{mm}$ for an irregular surface (modeling a keratoconic cornea) using radial basis functions, while our analysis on a similar case (see Fig. 6) shows that zonal Zernike fitting reaches between 2 and 4 orders of magnitude lower, depending on the size of the zonal mask we used. Furthermore, we compare our technique with that proposed in Ref. 22 on the same complex surface (Franke’s function), and we find that the rms of the reconstructed surface is four times lower in our model. Regarding the computational cost of our method, it is at least 2 orders of magnitude higher than the traditional Zernike calculation over the entire surface.

Our goal is to propose a combination of techniques that could be a good strategy for fitting, reconstructing, and resampling corneal height data. Local fitting reduces the influence of regions far from the calculation point. The implementation of the method also shows that the central corneal surface part is better evaluated than the outer parts, since calculation is more intensive in this zone and is not affected by peripheral irregularities. Thus, the presented method allows for better reconstruction of optical surfaces, mainly of the optical central zone, than the traditional one based on polynomial fitting of the whole surface.