Ophthalmic & Physiological Optics ISSN 0275-5408

Eigencorneas: application of principal component analysis to corneal topography Pablo Rodrıguez1, Rafael Navarro1 and Jos J. Rozema2,3 ICMA, Consejo Superior de Investigaciones Cientıficas-Universidad de Zaragoza, Facultad de Ciencias, Zaragoza, Spain, 2Department of Medicine and Health Sciences, University of Antwerp, Antwerp, Belgium, and 3Department of Ophthalmology, Antwerp University Hospital, Edegem, Belgium

1

Citation information: Rodrıguez P, Navarro R, Rozema JJ. Eigencorneas: application of principal component analysis to corneal topography. Ophthalmic Physiol Opt 2014; 34: 667–677. doi: 10.1111/opo.12155

Keywords: corneal model, corneal topography, principal component analysis, Zernike polynomials Correspondence: Rafael Navarro E-mail address: [email protected] Received: 7 March 2014; Accepted: 1 August 2014; Published Online: 14 September 2014

Abstract Purpose: To determine the minimum number of orthonormal basis functions needed to accurately represent the great majority of corneal topographies from a normal population. Methods: Principal Component Analysis was applied to the elevation topographies of the anterior and posterior corneal surfaces and central thickness of 368 eyes of 184 healthy subjects. PCA was applied directly to the input elevation data points and after fitting them to Zernike polynomials (up to 8th order, 8 mm diameter). The anterior and posterior surfaces, as well as right eye and left eye data, were analysed both separately and jointly. A threshold based on the amount of explained variance (99%) was applied to determine the minimum number of basis functions (eigencorneas) or degrees of freedom (DoF) in the population. Results: The eigenvectors directly obtained from elevation data resemble Zernike polynomials. The separate principal component analysis on the Zernike coefficients of anterior and posterior surfaces yielded 5 and 9 DoF, respectively. An additional reduction to 11 DoF (instead of 15 DoF) was achieved when performing a joint PCA that included both surfaces as well as central thickness. Finally, a further reduction was obtained by pooling right and left eye data together, to only 18 DoF. Conclusions: The combination of Zernike fit and Principal Component Analysis yields a strong reduction of dimensionality of elevation topography data, to only 19 independent parameters (18 DoF plus population average), which indicates a high degree of correlation existing between anterior and posterior surfaces, and between eyes. The resulting eigencorneas are especially well suited for practical applications, as they are uncorrelated and orthonormal linear combinations of Zernike polynomials.

Introduction The cornea is the first element of the optical system of the eye and as such it has a major effect on the optical and imaging properties of the eye. As it is part of the eye surface its topography can be measured with a reasonable reliability by a corneal topographer. Different types of basis functions were used in the literature to represent corneal elevation topography. Usually elevation is given by a matrix of sag z values sampled on a square grid of coordinates (x, y), which is not very practical mathematically. For reconstruction

from videokeratographic patterns (Placido disks) it is necessary to use local interpolating functions1 such as splines.2 For this reason modal representations such as Zernike,3–5 Fourier6 or Taylor7 analysis were introduced, which seem especially well suited for reconstructing elevation topography from Scheimpflug images.8 Other localised functions such as Bessel functions,9 adaptive arrays of Gaussians,10 or algorithms combining local and modal fitting11 were also proposed. But given the oversensitivity of most of these methods to irregular corneal shapes, the relatively stable and reproducible Zernike polynomials became the most

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

667

P Rodrıguez et al.

Application of PCA to corneal topography

popular methods for global representation despite their potential limitations.12 After topographical reconstruction from raw data, a comprehensive analysis is required by fitting the topography with optically significant shapes, such as spheres, conics, etc. From such analyses it is well known that the cornea is not spherical,13 lacks rotational symmetry and that surface models require biconics, three-axis ellipsoids14 or other alternative functions.15 Another experimental observation is that the conic constant becomes more negative as one increases the diameter of the cornea analysed. Several authors have proposed generalised conics16,17 to model this feature, but the most realistic models use a linear combination of a regular surface (conicoid, biconic, etc.) plus some polynomial functions (e.g. Zernike) to represent departures from that regular shape. Passing from the general to the canonical form of a three-axis ellipsoid it was possible to determine the position and orientation of the optical axis18 or higher order (4th and 6th) aspherical terms, as well as their changes with age.19 Generalizations to multizone models can be especially well suited to analyse postsurgical corneas.20 In summary, there is a wide variety of basis functions and models of increasing complexity. However, to some extent, the choice of the models seems somewhat arbitrary. The goal of this work was to define a new set of intrinsic basis functions to reconstruct corneal elevation topography. Secondly, we aimed to determine how many basis functions,21 or degrees of freedom (DoF), were required to accurately represent the great majority of corneas from a normal population. To this end principal component analysis (PCA) was applied to data from a previous study.19 PCA is a powerful method,22 especially suited to answer the two questions above. It produces a set of significant eigenvectors of the covariance matrix that are modes (or basis functions) intrinsic to the population analysed. PCA was used in many other fields including vision (eigenfaces proposed for face recognition task23) and optics (light polarisation24). It was also applied to the aberrations of the eye25,26 and to the analysis of corneal curvature maps as an illustrative example of robust statistical methods.27

Methods Materials and subjects This work uses previously published19 anterior and posterior corneal elevation (sagittal height) data, as well as the central corneal thickness, of 384 eyes of 192 healthy subjects recruited from the personnel of the Antwerp University Hospital and the University of Antwerp. All measurements were performed using a Pentacam Scheimpflug camera (www.pentacam.com). Subjects were between 22 and 70 years old, with a mean age of 42  (SD) 13 years, and without any previous ocular pathology or surgery. Other exclusion criteria were the wearing of hard contact lenses less than 1 month before testing and pregnancy. The study adhered to the tenets of the Declaration of Helsinki and received approval of the ethical committee of the Antwerp University Hospital (Ref. 7/6/24). Signed informed consent was obtained from participating subjects before testing. The Pentacam topographer provides a wide variety of data (elevation, power, pachymetry, etc.) In what follows we used the elevation z(x, y) of both anterior and posterior surfaces within a 8 mm corneal diameter, as well as the central thickness, as it is the minimum set of data providing a complete non-redundant description of the measured corneal geometry. Data analysis Left and right eyes were analysed independently, apart from an analysis on how the symmetry between both eyes of a subject may lead to the possibility of describing the topographies of both eyes simultaneously with only very few additional DoF to the single eye description. The analysis consists of several stages, which are summarized in Table 1. The first stage is to apply PCA to the raw elevation data points in order to decorrelate the data set by finding a complete orthogonal basis. PCA consists of computing the mean m and the covariance matrix C of the set of parameters over the population under study, and then diagonalising C. Then the first basis function is the mean m, whereas the following basis functions are the eigenvectors of matrix C, ordered according to decreasing magnitude of its

Table 1. Summary of data analyses

Input data & analysis

Number of parameters

Objective

Right eye surface elevation z(x,y) PCA Right eye Zernike coefficients PCA Right eye Zernike coefficients (re-ordering) Right eye Zernike coefficients & central thickness PCA Right and left eye Zernike coefficients & central thickness PCA

1257 45 45 91 182

Decorrelate the original data points, z(x,y) Get optimal basis on a variable space with n s Assess Zernike polynomials decorrelation capability Decorrelate anterior and posterior surfaces Decorrelate right and left eye data (mirror symmetry)

668

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

P Rodrıguez et al.

corresponding eigenvalue. In practice, it is common to consider only eigenvalues above a certain preset threshold to decrease the dimensionality. As a result PCA reduces the initial number of parameters n to n0 = 1 + e, with e the number of remaining eigenvectors so that ideally e  n. Here we apply a standard criterion in PCA in such a way that the first e eigenvectors account for 99% of the variance of the population. The resulting basis functions are orthogonal and quasi-complete and consist of linear combinations of the input basis functions. The main technical limitation of rigorous PCA treatment is that the number of parameters must be lower than the number of samples (i.e. n < s). As in our case s = 192 subjects, it will not be possible to compute more than 191 eigenvectors based on these data. Given that a typical elevation map consists of 75008000 points z(xi, yi), or 15 000–16 000 points if anterior and posterior surfaces are considered together, this condition is not met even though when a downsampling factor of 2 was applied to the original data in order to get covariance matrices of a reasonable (from a computational point of view) size. Nevertheless when n ≫ s, it is still possible to solve this as an underdetermined system and extract useful information. An alternative way to avoid these problems is to use a Zernike polynomial representation of the corneal elevation, as an intermediate stage, and then reduce the set of parameters (now Zernike coefficients) through PCA. This approach is especially warranted since we found that the eigenvectors of the raw elevation data resemble slightly distorted Zernike polynomials (see Results section). For this reason the elevation data of each corneal surface were fitted to a Zernike polynomials basis on a 8 mm diameter area. Zernike polynomials up to 8th order were used to guarantee a proper representation of the original data. Previous studies found an important correlation between the anterior and posterior surface of the cornea19. This means that a further decrease in dimensionality may be expected if PCA is applied to both surfaces simultaneously (joint PCA) rather than to each surface separately. To test this hypothesis the results of applying separate or joint PCA to both corneal surfaces were compared. In the separate PCA, each surface described by 45 Zernike coefficients was analysed separately. In the joint PCA, the corneal geometry is fully described by 45 9 2 (anterior and posterior surfaces) plus 1 (central thickness) parameters, so that the

Application of PCA to corneal topography

total number of parameters is n = 91. Thus, the joint PCA considers 91 parameters. For comparison purposes, we obtained another basis by re-ordering the Zernike polynomials according to the amount of population variance explained, using again the 99% variance threshold. A comparison between this basis and the PCA basis is useful to assess how far from optimality the Zernike polynomials are. Results PCA applied to raw elevation data The resulting basis functions, obtained directly from the initial elevation data z(x,y) are displayed on Figure 1 for the separate PCA applied to the anterior and posterior surfaces respectively. Note that after PCA the first basis function is always the average, followed by the first e orthonormal eigenvectors of the covariance matrix, ordered by the magnitude of their corresponding eigenvalues. The minimum and maximum values (blue and red in the colour maps), as well as the accumulated percentage of variance explained by the eigenvectors are given in Table 2. These colour maps suggest that the eigenvectors, especially the first ones, look like rotated and distorted Zernike polynomials such as tilt, defocus, astigmatism or coma. Nevertheless a number of artefacts are seen in the border areas of the higher eigenvectors were there were missing points in several subjects. The lack of data (s  n) also causes an artefact in the succession of the 400 most significant eigenvalues, normalised by the first (highest) eigenvalue (Figure 2). This artefact is clearly visible around the 200th eigenvalue (i.e. close to the number of subjects) as an abrupt drop to reach computing noise levels (109). This is a consequence of the insufficient number of subjects, making the system of equations to be solved underdetermined. Although it is fortunate that this cut occurs at very low levels (107), solving an underdetermined problem is risky in terms of artefacts and far from mathematically rigorous. Performing PCA directly on the elevation data z(x,y) poses important problems in our case, mainly due to the limitation in the total number of subjects in our data set. However, we have seen that the resulting eigenvectors resemble Zernike polynomials in both surfaces. Thus we decided fitting each anterior and posterior elevation topography to Zernike polynomials up to 8th order (45

Figure 1. Basis functions (eigencorneas) obtained for the anterior (top row) and posterior (bottom row) elevation topographies of all right eyes. The first basis function is the average, followed by the first 7 eigenvectors sorted according to decreasing eigenvalues.

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

669

P Rodrıguez et al.

Application of PCA to corneal topography

Table 2. Accumulated percentage of explained variance and limits of the colour maps in Figure 1 Average

#1

#2

#3

#4

#5

#6

#7

Anterior Variance (%) Maximum (mm) Minimum (mm)

– 1.116 0

87.2 0.053 0

92.8 0.095 0.068

96.1 0.050 0.099

98.1 0.070 0.095

99.1 0.071 0.081

99.3 0.084 0.147

99.5 0.109 0.079

Posterior Variance (%) Maximum (mm) Minimum (mm)

 1.400 0

76.3 0.056 0

86.7 0.062 0.066

92.3 0.059 0.083

94.9 0.075 0.084

96.7 0.054 0.120

97.9 0.062 0.121

98.3 0.086 0.159

for separate and joint PCA are very different. Differences tend to increase with the eigenvector number, and some of them change their rank in the order. The upper and lower limits of the colour maps, the accumulated percentage of variance explained by the eigenvectors and central thickness for joint PCA are summarised in Table 3. Zernike decomposition of the eigencorneas

Figure 2. Normalised eigenvalues corresponding to the eigenvectors of Figure 1: right eyes separate PCA, anterior (solid line) and posterior (dashed line) surfaces.

parameters). This is an intermediate step towards the final eigenvectors, allowing a drastic drop of dimensionality which guarantees the condition n  s. PCA applied to Zernike coefficients of surface elevation Figures 3 and 4 show the basis functions obtained for right eyes with separate (45 parameters per surface) and joint (91 parameters) PCA respectively. Only the 8 more significant basis functions are displayed: average and the 7 eigenvectors with the highest eigenvalues. As expected the results

The Zernike decomposition for the basis functions in Figure 3 is represented by bar plots in Figure 5. We have chosen the separate PCA as it allows for a simpler, more intuitive independent analysis of each surface. A quantitative description of the basis functions for both separate and joint PCA is given in Table 4. Both the average and the first eigenvector are simple for both surfaces. These are the two more relevant basis functions and both of them are mainly composed of a constant term z00 plus a revolution paraboloid z20 . If we take into account that the constant z00 is a relative position along the keratometric axis, we conclude that to a first approximation the mean and the first eigenvector e1 are revolution paraboloids. Nevertheless, other Zernike coefficients have small contributions: vertical and horizontal tilts (z11 and z11 ) confirm that the corneal surfaces are misaligned with the keratometric axis.18 In addition, the quadratic term z40 (first eigenvector) is a small correction to the parabolic shape which accounts for the well-known

Figure 3. Basis functions (eigencorneas) for anterior (top row) and posterior (bottom row) surfaces of the right eyes processed separately (each 45 parameters) using Zernike coefficients. The first basis function is the average, followed by the first 7 eigenvectors sorted according to decreasing eigenvalues.

670

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

P Rodrıguez et al.

Application of PCA to corneal topography

Figure 4. Basis functions for anterior (top row) and posterior (bottom row) surfaces obtained using joint processing of the Zernike coefficients of the anterior and posterior surfaces, as well as the central corneal thickness (91 parameters). The first basis function is the average, followed by the first 7 eigenvectors sorted according to decreasing eigenvalues.

Table 3. Accumulated percentage of explained variance, limits and central thickness of the eigenvectors in Figures 3 and 4 Average Anterior (separate PCA) Variance (%)  Maximum (mm) 1.116 Minimum (mm) 0 Posterior (separate PCA) Variance (%)  Maximum (mm) 1.401 Minimum (mm) 0 Anterior (joint PCA) Variance (%)  Maximum (mm) 1.116 Minimum (mm) 0 Posterior (joint PCA) Variance (%)  Maximum (mm) 1.401 Minimum (mm) 0 Central thickness (joint PCA) Value (mm) 0.549

#1

#2

#3

#4

#5

#6

#7

87.1 0.000 1.859

92.7 2.490 3.194

96.3 3.557 1.822

98.3 2.430 3.384

99.2 2.950 2.400

99.4 3.811 2.092

99.6 4.589 2.398

79.7 0 2.166

87.1 2.353 2.307

93.0 2.032 3.159

95.7 2.700 3.020

97.3 2.922 2.758

98.3 2.553 3.260

98.7 3.188 5.155

53.4 0.900 0

86.1 0.528 0

90.2 1.003 1.207

93.8 0.004 1.184

95.7 0.199 2.417

96.9 0.527 2.248

97.8 0.921 1.647

53.4 1.749 0

86.1 0.498 0.002

90.2 2.433 2.181

93.8 2.846 1.505

95.7 2.730 2.307

96.9 2.547 1.488

97.8 3.131 2.274

0.327

0.937

0.005

0.095

0.057

0.035

0.015

overall ellipsoidal shape of the cornea. The composition of the first eigenvector is quite similar to the mean: it is straightforward to realize that its role is the first and more important correction to the mean to account for the intersubject variability of the overall shape. To analyse the following more complex eigenvectors, we shall use the indexes n and m ðznm Þ as radial and angular frequencies respectively. Then using an analogy with music we will find the most significant term (i.e. the one with highest magnitude), which we shall call fundamental. The relative intensities (% contribution to the eigenvector RMS) of the fundamental components are listed in Table 4, where we can see that such percentage is higher than 50% in most cases. Therefore, the fundamental term in both the average and first eigenvector is n = 2, m = 0. In the same way, we will call harmonics those terms with the same m, and higher n. The 2nd eigenvector contains a variety of terms, most of them odd-symmetric (m = 1 or 1). The fundamental terms are tilts (z11 and z11 ) and their harmonics (coma z31

and z31 ). This is consistent with a relevant contribution of the corneal tilt and tip (and odd-symmetry) to explain the amount of variance (inter-subject variability) of the population. (Note that the relative contribution is given by the corresponding eigenvalue, so most of the variance is explained by the very first eigenvectors.) Again, tilts are the fundamental terms of 3rd eigenvector, but now the relative contribution of the horizontal tip is reversed in two ways: now z11 is higher in magnitude and negative. The fundamental contribution to eigenvectors 4th and 5th is toricity (astigmatism) which is a well-known geometrical feature of the cornea. The 4th fundamental is horizontal/vertical toricity, whereas the 5th fundamental is toricity along the diagonal. This order, according to their relative contribution to the variance is also consistent with the tendency observed experimentally of a vertical orientation (with the rule astigmatism). Again there are other minor contributions, especially for the posterior surface. The level of complexity increases considerably for the next eigenvectors, but a fundamental term can still be

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

671

P Rodrıguez et al.

Application of PCA to corneal topography

Figure 5. Bar plots with the distribution of normalised Zernike coefficients for the average and first 7 eigenvectors corresponding to Figure 3: separate PCA, right eyes, for the anterior (blue) and posterior surfaces (red).

determined. In the anterior cornea, the fundamental of the 6th eigenvector is z40 and smaller contributions by harmonics (z20 and z60 ) which contributes to modify the conic constant (more or less elongated ellipsoids). In the posterior 672

cornea, z40 is also a major contributor on the 6th eigenvector, but the fundamental mode is coma z31 . Additionally we can observe many other terms that increase the eigenvector complexity.

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

P Rodrıguez et al.

Application of PCA to corneal topography

Table 4. Fundamental components of anterior and posterior subsets of each eigenvector, and its percentage contribution to the eigenvector modulus of the eigenvectors depicted in Figures 3 and 4 Average Anterior (separate PCA) Fundamental z20 % RMS 50.8% Posterior (separate PCA) Fundamental z20 % RMS 51.1% Anterior (joint PCA) Fundamental z20 % RMS 50.8% Posterior (joint PCA) Fundamental z20 % RMS 51.1%

#1

#2

#3

#4

#5

#6

#7

z20 52.3%

z11 86.7%

z11 77.3%

z22 89.6%

z22 97.1%

z40 59.9%

z31 42.5%

z20 54.0%

z11 74.8%

z11 72.7%

z22 91.2%

z22 90.2%

z31 55.1%

z31 61.5%

z20 52.1%

z20 52.2%

z11 82.3%

z20 48.8%

z20 43.5%

z22 59.2%

z22 77.8%

z20 54.0%

z20 52.6%

z11 83.6%

z11 81.5%

z22 69.0%

z22 84.9%

z22 93.0%

Finally, eigenvector 7th is highly complex and, in the anterior cornea, the fundamental term (z31 ) is not much higher than the rest. The same fundamental term is observed in the posterior cornea, although in this case its preponderance is higher (see Table 4). We performed a similar analysis for the results of the joint PCA. The quantitative results are included in Table 4. In this case, the contribution of the fundamental component to the eigenvector RMS is calculated for the anterior and posterior subsets of each eigenvector independently. The strong correlations existing between the anterior and posterior surfaces, and the inclusion of corneal thickness to the analysis make the results more entangled and less intuitive to interpret. Even though the numbers are quite different from the separate PCA, the fundamental contributions are basically the same, at least for the initial eigenvectors. If we consider the complete eigenvectors (i.e. including central thickness, not displayed on Table 4 for the sake of clarity), the fundamental component for the mean and first eigenvectors is the parabolic coefficient z20 added with a constant z00 . However, the fundamental term in the 2nd eigenvector is central thickness (even though it is not displayed in Table 4) with a contribution of 90.1%, and lower contributions of z00 and z20 . The 3rd fundamentals are tilts (and harmonics) again. As the eigenvector number increases we observe a complicated pattern with a clear fun-

damental in one surface but a highly mixed distribution in the other surface. Re-ordered Zernike polynomials Table 5 shows the results of re-ordering the Zernike polynomials according to the amount of variance explained for right eye data (anterior and posterior corneal surfaces). The Table is divided in three rows, the first of which is the ordering number. The second row shows the Zernike polynomial and the third row contains the cumulative percent contribution of that polynomial to the global sample variance. The modes that account for most of the variance are: the constant z00 (average sag); then the revolution paraboloid z20 , followed by the tip/tilt terms; then even-symmetric (toricity z22 and z22 ) and odd-symmetric (z31 and z31 ) polynomials; then revolution symmetry aspherical term z40 ; trefoil, etc. Most of these Zernike polynomials were also the fundamental modes in the PCA eigenvectors. Dimensionality reduction The minimum number of eigenvectors needed to account for 99% of the database variance will be referred to as degrees of freedom (DoF). Table 6 summarizes the results obtained with this threshold criterion for the cases analysed

Table 5. Zernike polynomials ordered according to their relative contribution (from highest to lowest) to the total sample variance. Results correspond to right eyes, anterior and posterior corneal surfaces

Anterior Zernike % Variance Posterior Zernike % Variance

#1

#2

#3

#4

#5

#6

#7

#8

#9

#10

z00 63.2

z20 87.1

z11 91.9

z11 95.1

z22 97.4

z22 98.5

z31 99.0

z31 99.3

z40 99.4

z33 99.5

z00 55.8

z20 79.3

z11 85.9

z11 92.1

z22 95.0

z22 96.7

z31 97.4

z31 98.0

z40 98.5

z44 98.7

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

673

P Rodrıguez et al.

Application of PCA to corneal topography

Table 6. Summary of input parameters and number of eigenvectors needed to account for 99% of the sample variance (degrees of freedom)

Right eye elevation z(x,y) Right eye Zernike coefficients (PCA) Right eye Zernike coefficients (re-ordering) Right eye Zernike coefficients plus central thickness Right and left eye Zernike coefficients plus central thicknesses

in the previous sections, including the re-ordered Zernike polynomials. Under that criterion, 5 and 9 eigenvectors are needed to represent the anterior and posterior corneal surfaces respectively using the Zernike PCA basis, whereas 8 and 12 eigenvectors are needed with the re-ordered Zernikes. It is interesting to note that a similar result was obtained with the raw elevation data (5 and 10 eigenvectors, respectively) in spite of the mathematical and computing issues pointed out before. As for the joint PCA, the total number of joint eigenvectors needed to represent the geometry (both surfaces and thickness) was 11, which is substantially lower than the 15 eigenvectors required to represent the complete corneal geometry separately (i.e. 5 anterior, 9 posterior and 1 central thickness). This is consistent with the high correlation (redundancy) between both corneal surfaces,19 which is removed by the joint PCA. Discussion The results presented above show the feasibility of a strong reduction of the number of parameters needed to represent the corneal shape in a population of normal eyes. PCA provides the optimal linear solution in several aspects. The eigenvectors of the covariance matrix are not only orthogonal modes, but also fully decorrelated from each other. In the case that the input parameters have Gaussian distributions, then this guarantees a complete statistical characterization of that population, and that the coefficients on the new basis functions (eigencorneas) are statistically independent from each other. The whole population is described by a constant mean plus a linear combination of an orthonormal, decorrelated and reduced set of basis functions. When PCA is applied directly to the elevation data points, we obtain a drastic reduction of the dimensionality (number of independent parameters) of more than two orders of magnitude, and the resulting eigenvectors resemble distorted and rotated Zernike polynomials. The application of PCA directly to the elevation data points seems highly promising, but in our case its reliability is limited since the number of subjects in our study is lower than the number of parameters. Nevertheless, these results suggest 674

Input parameters

Anterior cornea

Posterior cornea

Complete cornea

1257 45 45 91 182

5 5 8 – –

10 9 12 – –

16 15 21 11 18

that Zernike polynomials represent a natural set of basic shapes and symmetries that are inherent to the corneal shape. Zernike fitting was highly powerful to represent elevation data and the dimensionality was dramatically reduced. The low-pass filtering effect of Zernike polynomials can potentially affect both noise and signal, but in normal corneas, regular and smooth, one can expect that the signal is almost unaffected so that most of the filtered component is noise. Furthermore, the results of Table 5 suggest that a simple re-ordering of the Zernike polynomials may provide nearly optimal results. The re-ordered Zernike polynomials provide a sub-optimal (nearly optimal) generic orthogonal basis quite interesting in practical terms. Note that this reordering was performed following the same criterion as in PCA: the only difference with PCA is that in this case the covariance matrix was only reordered but not diagonalised. Therefore, we can get a progressive reduction of dimensionality by applying successive stages: (1) Zernike fit acts as a powerful intermediate stage that permits us reducing the dimensionality (45 + 1 + 45 parameters), (2) a simple reordering of Zernike polynomials according to their associated population variance (as in PCA) and then applying a threshold permitted a further reduction (8 + 1 + 12 parameters), (3) finally, the combination of Zernike fit (1) and PCA yielded a stronger reduction of dimensionality (5 + 1 + 9 parameters.) The final number or parameters, in cases (2) and (3), depends on the eigenvalues threshold (99% of variance in our case). We only presented results for the right eye, but the same analysis was applied to left eye data. The results were similar with some numerical differences in the eigenvalues and eigenvectors. As a result, the DoF (number of basis functions) are 5 (right or left eye data) for the anterior surface and 9 (right or left eye data) for posterior surface. A central result of this study is that the joint PCA applied to both surfaces and corneal thickness further reduces dimensionality to a total number of 11 for right or left eye data. Note that this further reduction of dimensionality cannot be attained through a simple re-ordering of Zernike terms, since it is due to the high correlation between anterior and posterior surfaces. Thus this result is exclusive of the joint PCA method. This means that we only need 11

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

P Rodrıguez et al.

independent parameters to describe the overall geometry (both surfaces and thickness) in a population of corneas with a reasonable accuracy. Most of this reduction can be explained by the strong correlation existing between the anterior and posterior topographies.19 Note that the joint eigenvectors are more complex (91 components instead of 45) and hence they are able to carry more information. We also tested the possibility of a further reduction of dimensionality due to the mirror symmetry between left and right eyes by taking together left and right eyes to apply a conjoint PCA. This was implemented by duplicating the number of parameters (91 for right eye plus 91 for left eye) describing the complete set of data for each subject. The resulting DoF were 18 which means a more modest but still significant (compared to 22 if right and left eye data are analysed independently) reduction of dimensionality. Such further reduction (4 DoF less) can be realised as the removal of the redundancy existing due to the symmetry between right and left eyes. Note that for perfect mirror symmetry we expect to have no more than 12 DoF (only one more to account for horizontal flip). Then the modest reduction of DoF attained so far suggests that the redundancy (and hence symmetry) is moderate. These conjoint eigenvectors are vectors of 182 components (i.e. 91 per eye). This possibility of further reductions of dimensionality taking advantage of the strong correlations between anterior and posterior surfaces or between right and left eyes is probably the main advantage of PCA eigenvectors over the Zernike polynomials. PCA is an efficient method to reduce the dimensionality of large data sets such as the corneal topographies of the population of this study. The 99% variance of our database was explained by only 18 parameters. Nevertheless, the method has limitations that one has to take into account. A limitation inherent to the method is that the eigenvectors are specific for a particular population of eyes, and this imposes the same generalization problems as any other statistical approach. Another fundamental limitation for the direct application of PCA to the elevation data points is the lack of data, as we found artefacts both in eigenvectors and eigenvalues either when the total number of subjects is lower than the number of parameters (s < n), or when there were topographies with peripheral missing points. The analysis of the influence of the topographer is beyond the present work. Nevertheless, comparative studies in the literature28,29 suggest that Scheimpflug systems such as Pentacam show excellent repeatability and reasonable agreement with classic videokeratoscope instruments, except for higher angular frequencies m. Such specific difference is probably due to the coarser angular sampling (and finer radial sampling) of Scheimpflug system, as opposed to the finer angular (and coarser radial) sampling of Placido disks. Thus, it is reasonable to expect that these

Application of PCA to corneal topography

different sampling schemes (and different topography reconstruction algorithms) might produce differences in the eigenvector modes derived from data of different instruments. This could be a subject of future work. We believe that these eigenvectors, which are a set of orthonormal decorrelated vector modes, may be useful in a variety of applications. On the one hand, their low dimensionality facilitates a simple, fast and robust implementation, which consists of subtracting the mean and fitting the eigenvectors to the residual. The possibility of a complete characterization of the corneal shape with just a few parameters may be useful as well, such as in data compression.30 Furthermore, each eigenvector coefficient means a progressively finer correction to the initial rough approximation given by the mean. Tables 1 and 2 show that the first eigenvalue tends to be more than one order of magnitude higher than the second one. This suggests that, in many cases, only one coefficient (eigenvector #1) may be enough to provide a reasonable approach to the residual, which added to the mean allows a reasonable reconstruction of the corneal topography. The following eigenvectors provide finer approaches, but their relative contributions are progressively smaller. On the other hand, statistical independence (assuming Gaussian distribution) enables interesting applications, such as the numerical generation of synthetic, but realistic corneas from random numbers to perform Monte Carlo simulations, or building synthetic eye models.31 Finally, here we obtained eigencorneas for a normal population. It is reasonable to expect that deformations associated to pathologies such as keratoconus, or as a result of corneal surgery may yield significantly different eigencorneas. As an illustration of this potential application, we have performed a similar analysis on another population of

Figure 6. RMS fit residual of Zernike elevation data to the significant (99% variance) eigenvectors obtained with the normal population for normal (black circles) and keratoconic (red crosses) eyes.

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

675

P Rodrıguez et al.

Application of PCA to corneal topography

normal eyes. The basis functions obtained have then been used to fit the Zernike elevation data of these normal eyes as well as another population of keratoconic eyes, obtaining the results displayed in Figure 6. These results would be potentially useful in classification, early detection and hence it will be part of our future work. Acknowledgements This research was supported by the Spanish Ministry of Economıa y Competitividad and the European Union, grant FIS2011-22496, to R.N. and by the Flemish Agency for Innovation by Science and Technology, grant IWT/ 110684, to J.R. Disclosure The authors report no conflicts of interest and have no proprietary interest in any of the materials mentioned in this article. References 1. Klein SA. A corneal topography algorithm that produces continuous curvature. Optom Vis Sci 1992; 69: 829–834. 2. Halstead MA, Barsky BA, Klein SA & Mandell RB. A spline surface algorithm for reconstruction of corneal topography from a videokeratographic reflection pattern. Optom Vis Sci 1995; 72: 821–827. 3. Carroll JP. A method to describe corneal topography. Optom Vis Sci 1994; 71: 259–264. 4. Schwiegerling J, Greivenkamp JE & Miller JM. Representation of videokeratoscopic height data with Zernike polynomials. J Opt Soc Am A 1995; 12: 2105–2113. 5. Iskander DR, Collins MJ & Davis B. Optimal modeling of corneal surfaces with Zernike polynomials. IEEE Trans Biomed Eng 2001; 48: 87–95. 6. Hjortdal JØ, Erdmann L & Bek T. Fourier analysis of video-keratographic data. A tool for separation of spherical, regular astigmatic and irregular astigmatic corneal power components. Ophthalmic Physiol Opt 1995; 15: 171–185. 7. Howland HC, Glasser A & Applegate R. Polynomial approximations of corneal surfaces and corneal curvature topography. Tech Dig Ser Opt Soc Am 1992; 3: 34–37. 8. Sicam VA & Van der Heijde RG. Topographer reconstruction of the nonrotation-symmetric anterior corneal surface features. Optom Vis Sci 2006; 83: 910–918. 9. Trevino JP, G omez-Correa JE, Iskander DR & Chavez-Cerda S. Zernike vs. Bessel circular functions in visual optics. Ophthalmic Physiol Opt 2013; 33: 394–402. 10. Martınez-Finkelshtein A, Ramos-L opez D, Castro GM & Ali o JL. Adaptive cornea modeling from keratometric data. Invest Ophthalmol Vis Sci 2011; 52: 4963–4970.

676

11. Espinosa J, Mas D, Perez J & Illueca C. Optical surface reconstruction technique through combination of zonal and modal fitting. J Biomed Opt 2010; 15: 026022. 12. Smolek MK & Klyce SD. Zernike polynomial fitting fails to represent all visually significant corneal aberrations. Invest Ophthalmol Vis Sci 2003; 44: 4676–4681. 13. Kiely PM, Smith G & Carney LG. The mean shape of the human cornea. Opt Acta (Lond) 1982; 29: 1027–1040. 14. Burek H & Douthwaite WA. Mathematical models of the general corneal surface. Ophthalmic Physiol Opt 1993; 13: 68–72. 15. Kasprzak HT & Jankowska-Kuchta E. A new analytical approximation of corneal topography. J Mod Opt 1996; 43: 1135–1148. 16. Kasprzak HT & Iskander DR. Approximating ocular surfaces by generalised conic curves. Ophthalmic Physiol Opt 2006; 26: 602–609. 17. Rosales MA, Juarez-Aubry M, L opez-Olazagasti E, Ibarra J & Tepichın E. Anterior corneal profile with variable asphericity. Appl Opt 2009; 48: 6594–6599. 18. Navarro R, Gonzalez L & Hernandez JL. Optics of the average normal cornea from general and canonical representations of its surface topography. J Opt Soc Am A 2006; 23: 219–232. 19. Navarro R, Rozema JJ & Tassignon MJ. Optical changes of the human cornea as a function of age. Optom Vis Sci 2013; 90: 587–598. 20. Gonzalez L, Hernandez-Matamoros JL & Navarro R. Multizone model for postsurgical corneas. analysis of standard and custom LASIK outcomes. J Biomed Opt 2008; 13: 044035– 1,12. 21. Iskander DR, Morelande MR, Collins MJ & Buehren T. A refined bootstrap method for estimating the Zernike polynomial model order for corneal surfaces. IEEE Trans Biomed Eng 2004; 51: 2203–2206. 22. Jolliffe IT. Principal Component Analysis 2nd edition, Springer-Verlag: New York, 2002. 23. Turk M & Pentland A. Eigenfaces for recognition. J Cogn Neurosci 1991; 3: 71–86. 24. San JI & Gil JJ. Invariant indices of polarimetric purity: generalized indices of purity for n9n covariance matrices. Opt Commun 2011; 284: 38–47. 25. Porter J, Guirao A, Cox IG & Williams DR. Monochromatic aberrations of the human eye in a large population. J Opt Soc Am A 2001; 18: 1793–1803. 26. Thibos LN, Hong X, Bradley A & Cheng X. Statistical variation of aberration structure and image quality in a normal population of healthy eyes. J Opt Soc Am A 2002; 19: 2329– 2348. 27. Locantore N, Marron JS, Simpson DG, Tripoli N, Zhang JT & Cohen K. Robust principal component analysis for functional data. Test 1999; 8: 1–65. 28. Read SA, Collins MJ, Iskander DR & Davis BA. Corneal topography with Scheimpflug imaging and

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

P Rodrıguez et al.

videokeratography: comparative study of normal eyes. J Cataract Refract Surg 2009; 35: 1072–1081. 29. McAlinden C, Khadka J & Pesudovs K. A comprehensive evaluation of the precision (repeatability and reproducibility) of the Oculus Pentacam HR. Invest Ophthalmol Vis Sci 2011; 52: 7731–7737.

Application of PCA to corneal topography

30. Schneider M, Iskander DR & Collins MJ. Modeling corneal surfaces with rational functions for high-speed videokeratoscopy data compression. IEEE Trans Biomed Eng 2009; 56: 493–499. 31. Rozema JJ, Atchison DA & Tassignon MJ. Statistical eye model for normal eyes. Invest Ophthalmol Vis Sci 2011; 52: 4525–4533.

© 2014 The Authors Ophthalmic & Physiological Optics © 2014 The College of Optometrists Ophthalmic & Physiological Optics 34 (2014) 667–677

677

This document is a scanned copy of a printed document. No warranty is given about the accuracy of the copy. Users should refer to the original published version of the material.

Eigencorneas: application of principal component analysis to corneal topography.

To determine the minimum number of orthonormal basis functions needed to accurately represent the great majority of corneal topographies from a normal...
978KB Sizes 2 Downloads 5 Views