## Abstract

This work presents a compact statistical model of the retinal image quality in a large population of human eyes following two objectives. The first was to develop a general modal representation of the optical transfer function (OTF) in terms of orthogonal functions and construct a basis composed of cross-correlations between pairs of complex Zernike polynomials. That basis was not orthogonal and highly redundant, requiring the application of singular value decomposition (SVD) to obtain an orthogonal basis with a significantly lower dimensionality. The first mode is the OTF of the perfect system, and hence the modal representation, is highly compact for well-corrected optical systems, and vice-versa. The second objective is to apply this modal representation to the OTFs of a large population of human eyes for a pupil diameter of 5 mm. This permits an initial strong data compression. Next, principal component analysis (PCA) is applied to obtain further data compression, leading to a compact statistical model of the initial population. In this model each OTF is approximated by the sum of the population mean plus a linear combination of orthogonal eigenfunctions (eigen-OTF) accounting for a selected percentage (90%) of the population variance. This type of models can be useful for Monte Carlo simulations among other applications.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Among the different optical and image quality functions, wavefront aberration W is the most extensive way to objectively characterize the optical performance of the human eye. The robustness and reliability of ocular aberrometers, as well as the popularization of commercial systems [1], caused the publication of a high number of both basic [2,3] and clinical studies [4–8]. A number of researchers studied the relationship between W and the most relevant magnitudes used in the clinical practice, such as visual acuity VA [9,10] or refractive error Rx [11]. Despite the complex, non-linear relationship between wavefront aberration and retinal image, these studies showed a reasonable overall correlation between wavefront-based metrics, typically calculated with the coefficients of a Zernike polynomial expansion [12], and clinical Rx or visual acuity. A more direct estimation of visual performance, however, would require analyzing the retinal image quality by means of the frequency response given by the Optical Transfer Function OTF [13] or the impulse response given by the Point Spread Function PSF [14,15]. Different authors proposed sets of analytical functions for OTF expansions before [16–18], but each suffered from important drawbacks, such as the lack of orthogonality, etc. Another work proposed an orthogonal expansion in which the first function was the OTF of the perfect system [19], causing OTFs close to the diffraction limit to be highly compressed, whereas OTFs of low-quality optical systems would require a large number of modes (basis functions).

To this end, the current work aims to develop a statistical model of the OTF image quality metric based on a large population of normal human eyes. Individual OTFs depend on spatial frequency, which can be described as a linear expansion of two-dimensional complex sinusoids (plane waves) that form a complete orthogonal basis. This would form a highly compact model to describe a large population of OTFs, able to maximally compress the initial individual data.

To avoid these issues, the departure point of this work is different by using a numerical approach to develop an orthogonal basis able to fit individual OTFs. Once the individual OTFs can be represented as a set of expansion coefficients, the second part consists of statistically modeling the whole population by applying standard Principal Component Analysis (PCA). This technique has been used before in the field of visual optics, in particular to model the statistics of wavefront aberrations [2,20,21] or corneal topography [22] of a population. In the case of the impulse-response Point Spread Function (PSF), PCA has been applied before to assess the performance of astronomical instrumentation [23]. As for the ocular OTF, previous statistical models were restricted to the computation of the mean OTF [24,25], even though the large inter-subject variability makes these mean-based models too simplistic. PCA has also shown to be a powerful tool enabling the implementation of realistic synthetic eye models [26].

## 2. Methods

The main stages of our method are outlined in Fig. 1. It starts with the theoretical and numerical steps required to compose the general model in the form of an orthogonal basis functions to represent individual OTFs in a compact way. Next, the model is applied to a population database by projecting each OTF on the basis functions to obtain a new compact representation, and subsequently perform the statistical model through standard PCA.

#### 2.1 Orthogonal functions

### 2.1.1 Initial basis functions

The original wavefront database consists of a vector of 44 Zernike coefficients for each eye (i.e. up to 8th order, excluding piston). Assuming pupil apodization by the Stiles-Crawford effect (SCE) remains small for the 5 mm pupil considered here, and assuming constant transmittance inside the pupil, the 2D complex pupil function P can be constructed as:

^{th}order; $k = {\raise0.7ex\hbox{${2\pi }$} \!\mathord{\left/ {\vphantom {{2\pi } \lambda }} \right.}\!\lower0.7ex\hbox{$\lambda $}}$ is the wave number. Using the fact that complex Zernike polynomials from a complete orthogonal basis on the unit circle able to represent any complex pupil function

*P*, we can rewrite it as: where ${C_k}$ are complex Zernike polynomials CZP [27]: with $N_n^m$ and $R_n^{|m|}$ the normalization factors and radial Zernike polynomials, respectively, as defined in the ISO24157:2008 standard [28]. Once the complex pupil function is expressed in the orthogonal basis (Eq. 3), the Optical Transfer Function OTF is obtained by cross-correlating the CZP coefficients [27]:

Although this basis is complete, it is not orthogonal and highly redundant. It is also
very large, requiring K^{2} functions to describe an OTF,
or the square of the number of basis functions needed for the
pupil function *P*. This number of
basis functions remains high, even after removal of all duplicate
cross-correlations. Such a number was to be expected, given the
strong non-linearity of the complex exponential in Eq. (1),
requiring about double the number of coefficients for the complex
pupil function (91 CZPs, 12th order) compared to the wavefront
*W* (44 ZPs, 8th order given by the
aberrometer) in order to have a sufficiently accurate expansion in
Eq. (3). Starting with a 12^{th} order expansion of
*P*, and considering all
cross-correlations between terms while removing simple
redundancies, this leads the huge number of 1184 basis functions.
We could expect to get additional reduction by removing possible
non-Hermitian functions, since the OTF is Hermitian by definition.
But since the SVD already removes all redundancies and provides an
orthogonal basis. Hermiticity testing of the functions was not
necessary.

### 2.1.2. Singular value decomposition

Singular Value Decomposition (SVD) is a standard powerful method of linear algebra that is especially well-suited to solve problems with high redundancy and a lack of orthogonality. SVD splits a *m* x *n* matrix **A** into a product of three matrices **A** = **UΣV ^{T}**, where

**U**(

*m*x

*m*) and

**V**(

*n*x

*n*) are unitary orthogonal matrices and

**Σ**is a

*m*x

*n*diagonal matrix containing the positive singular values in rapidly decreasing order. The non-zero singular values form the fundamental column subspace with a dimension

*r*that would be ideally much lower than the initial dimension

*N*.

Conversely, the null subspace corresponds to the subset of singular values equal to zero. This property is important in this case of a huge, highly redundant matrix. In this numerical implementation the OTF frequency domain is sampled in a square grid of 201 × 201 = 40401 points, which can be reduced to 31417 points (by a factor equal to π/4) considering only points below the cut-off frequency. Given the 1184 initial functions, this means that the size of the initial matrix A is 31417 (samples) by 1184 (functions). Another reduction can be made by imposing a threshold to the singular values slightly greater than zero by computing the sum of the singular values and selecting the set of values adding up a given percentage of the total sum. For 80% of the total this reduce the number from 1184 to 114 basis functions, or about 10% of the initial size and only slightly higher than the initial number of 91 CZPs. Meanwhile, a more conservative threshold of 95% corresponds with 196 basis function, or 16.5% of the initial size.

It is worth remarking that the resulting orthogonal basis functions are general, meaning they can be used for modelling the OTF of any optical system with a circular pupil, such as human eyes. This basis consists of numerical samples, rather than analytical functions reported before [19], but both approaches are equivalent. The coefficients can be computed through standard methods, by either determining the projections on each basis function (inner product), or through a least-squares fit.

#### 2.2. Database

This work uses the monochromatic wavefront aberration data [21] of 333 healthy right eyes and 330 healthy left eyes of West-European Caucasian subjects measured at the Antwerp University Hospital. The study adhered to the tenets of the Declaration of Helsinki and received approval of the ethical committee of the Antwerp University Hospital. Signed informed consent was obtained from participating subjects before testing. The cohort consisted of 42.9% men and 57.1% women with an average age of 40.8 ± 11.0 years (range: 20.2 - 59.6 years). Exclusion criteria were previous ocular pathology or surgery, an intraocular pressure higher than 22 mmHg, wearing rigid contact lenses less than 1 month before testing and pregnancy. In addition, only eyes with refractive errors within the range of [−5 D, + 5 D] for sphere and [−1.5 D, + 1.5 D] for astigmatism were considered. Finally, 10 potential outliers with wavefront aberrations larger than the mean, plus 3 times the standard deviation of the higher order aberrations HOA (excluding defocus and astigmatism) RMS wavefront error were removed. All measurements were performed using an iTrace ocular aberrometer (www.traceytechnologies.com) under natural pupil dilation with a measurement diameter of 5 mm. Left and right eyes were analyzed independently. In the following only the right eye data was presented, given that the left eye results are basically equivalent due to mirror symmetry.

#### 2.3. Principal component analysis

Principal Component Analysis (PCA) [29] was performed to find the simplest global model of the population with the lowest dimensionality (or degrees of freedom, DoF). This model permits expressing the image quality of any eye in the population as a linear combination of the mean *M* plus a limited number of *J* orthogonal eigenvectors ${E_j}$. This means that for eye number *i* the OTF in polar coordinates can be reasonably approximated by:

*M*and the OTF eigenfunctions (or eigenOTFs) ${E_j}$ are constant for a given population, each individual OTF can be fully described by a low-dimensional coefficient vector ${{\bf c}^i} = \{{c_j^i} \}$. If the variables in the database are normally distributed, PCA provides a complete compact statistical model of the specific population analyzed. To our knowledge, the combination of physically meaningful OTF basis functions and PCA is a novel approach to model the retinal image quality.

The initial set of 2D OTFs is given in a 201 × 201 Cartesian grid. The first step is to
compress the modal representation to a limited set of coefficients
(114 or 196 coefficients for 80% or 95% SVD threshold respectively).
Next, this modal representation is used for the statistical analysis
of the OTF population. PCA is subsequently applied for both further
data compression and obtaining a global statistical model of our
population Gaussian distributions fully characterized by the mean
**M** (average vector of coefficients) and the covariance
matrix **C**. The PCA consists of the diagonalization of
**C** in the form
**U ^{−1}CU** =

**Λ**, with

**Λ**a diagonal matrix formed by eigenvalues ${\lambda _j}$ in decreasing order. The unitary matrix

**U**is formed by the eigenvectors of the covariance and represents a change of basis through rotation to align the distribution along the axes of the Gaussian distribution. The eigenvalues are the variances along these axes, and their sum is the total variance of the distribution. Applying the rotation

**U**to the basis functions lead to the new orthogonal eigenfunctions (‘eigenOTFs’) of the population. Again a threshold to the cumulative variance (sum of eigenvalues) is applied to further reduce dimensionality. Finally, Eq. (5) is applied to represent all OTFs in the database. Note that this leads to another modal representation in which the first mode is the mean. To obtain a realistic model, the distribution should not only be reasonably Gaussian, but the number of samples (303 individuals) should be greater than the size of the covariance matrix. The initial SVD threshold must therefore be high enough to reduce the number of basis functions below the limit imposed by the dataset size.

## 3. Results

#### 3.1 Orthogonal OTF basis

The first part focuses on general results as the orthogonal modes (basis functions) can be used to obtain a compact modal representation of any OTF function.

Since SVD orthogonalization is numerical, the final orthogonal modes are arrays of discrete numerical matrices (201 × 201 elements). Some of the resulting basis functions are depicted in Fig. 3, with the upper half showing the modulus and the lower half the phase of these complex OTF functions. The first function again represents the OTF of a perfect optical system, but the rest changed substantially compared to Fig. 2 and can no longer be easily related to stand-alone Zernike modes.

A fundamental aspect to obtaining a compact representation is to apply a threshold to the singular values, as illustrated in Fig. 4. The curve shows the cumulative sum of singular values as a function of the number of modes (basis functions) included. This demonstrates the power of the SVD to reduce dimensionality, as well as the possibility to retain the most significant singular values and their associated basis functions by applying a suitable threshold. In this way, a high threshold (99%) will provide a high fidelity representation but at the cost of needing 271 modes, and conversely saving modes by a lower threshold will yield a lower fidelity (80% requiring only 114 modes).

#### 3.2. EigenOTFs for normal eyes

Applying PCA to the population data for a 5** **mm pupil diameter leads to the
average OTF (first mode) and first 10 eigenfunctions for our
population (Fig. 5),
consisting of both the modulus (Modulation Transfer Function, MTF) and
phase (Phase Transfer Function, PTF). Here too, a threshold may be
applied to reduce the dimensionality even further, this time in the
form of the eigenOTF functions (PCA modes; Fig. 6). Using similar criteria, the
number of modes can be reduced to about 20 (for 80% cumulative
variance), or more for a high-fidelity model. The important advantage
here is that the model can be adapted to the fidelity required. The
most simplistic model would be to only use the mean (0% variance, 100%
deterministic), or increase the number of eigenOTFs if the application
requires it.

Figure 7 shows the numbers of total modes for different combinations of SVD and PCA thresholds, which ranges from 20 modes for the most compact model (80% and 80% thresholds) to more than 120 for the highest fidelity model (99%, 99% thresholds). Therefore, a strong dimensionality reduction can be achieved with an adequate combination of both SVD and PCA thresholds.

#### 3.3. OTF reconstruction

One of the most interesting applications of such statistical models is generating a virtually unlimited number of synthetic individuals through realizations Monte Carlo simulations [26]. To this end the fidelity of the OTF representation was verified in terms for a low number of eigenOTFs by analyzing both average RMS fit error for the whole population, as well as a couple of randomly picked examples of the reconstructed OTFs. To obtain the modal representation a standard lineal least squares fit was applied, and a weighting function proportional to 1/ρ (radial frequency) was introduced to consider the fact that the number of the samples to fit for each radial frequency is proportional to 2πρ. Without weighting, the higher frequencies dominate the fit and the lower frequencies tend to have higher fit errors. This weighting strongly improves the reconstruction fidelity of the low and mid frequencies, which are highly relevant in human vision. Furthermore, by definition the OTF is normalized at zero frequency, and hence it is essential to obtain extremely high reconstruction fidelity for frequencies close to zero. Figure 8 compares the original versus reconstructed OTF modulus (MTF) for two subjects (upper and lower panels respectively).

As the number of modes necessary to approximate OTFs increases, the OTF departs from the diffraction limit [19]. The wavier the OTF, the harder the reconstruction. For this reason, this is a highly demanding problem and the reconstructions are not perfect, especially for high frequencies. Regardless, the model is able to reproduce the overall trends for this population of normal eyes. Note that the OTFs are far from diffraction limited since human eyes are strongly aberrated.

The reconstruction errors were analyzed for the whole population and different combinations of SVD and PCA thresholds. The metric used was the volume (computed in the same way as the Strehl ratio (i.e. normalized by the volume under the diffraction-limited OTF) of the difference between the original and the reconstructed OTFs. Figure 9 compares the reconstruction error histograms for the whole population for the two extreme (lowest and highest) applied thresholds. The statistical median for different combinations of SVD and PCA thresholds are listed in Table 1. As expected, as we increase both thresholds, both the mean (histogram shift) and the dispersion (histogram width) are reduced. When we pass from (80%, 80%) to (99%, 99%) thresholds, then the mean, median and standard deviation are roughly halved.

It is important to point out that increasing the accuracy of the OTF reconstruction is no longer useful once the reconstruction error is well below the noise level of the ocular wavefront or direct MTF measurements. Several sources of wavefront variability have been studied in the past, including residual accommodation fluctuations [30], pupil alignment [31], tear-film dynamics [32] and inter-device variability [1,33]. We have estimated our noise level by performing 6 repeated measurements on a single normal subject with the same aberrometer. The OTF of the average Zernike coefficients is then considered our “true OTF” and the OTF from the individual measurements the “reconstructed OTFs”. The average of the calculated reconstruction errors was 0.113 (range: 0.067 to 0.151), meaning that our lowest thresholds are already in the limit of the expected natural fluctuations.

From the above analysis we conclude that for human eyes (and possibly other strongly aberrated optical systems) we would need to choose among a minimum number of modes (see Fig. 7) compatible with the highest fidelity possible (Table 1). For example, we can attain a reconstruction error below 0.1 with 38 modes (that is 90% SVD and 90% PCA thresholds).

## 4. Discussion

This study presents a compact modal representation of the optical transfer function OTF in terms of a complete orthogonal basis, which differs from previous approaches [16–19]. It starts from an initial basis formed by cross-correlations between complex Zernike polynomials, considering the physics of optical image formation, to form an orthogonal basis on a circular pupil. Like in other approaches before [18], the resulting functions are non-orthogonal and highly redundant. To address this issue a Singular Value Decomposition SVD was applied so that the final basis is no longer analytical. This difference with other approaches is mostly of theoretical importance since in practical use one would have to discretize the basis functions. Other purely mathematical approaches [19] select the first function and subsequently build a basis with functions orthogonal to it. Interestingly, in most of these approaches one of the basis functions is the diffraction-limited OTF, allowing a strong compression rate for the OTF of well-corrected optical systems.

Besides the OTF defined on a circle of unit radius (e.g. the pupil) it would also be possible to develop a modal representation of the PSF. This would require the use of Bessel or related functions with infinite support, which seems a considerable drawback. Instead, one could apply image formation theory, such as the extended Nijboer-Zernike theory [34,35], to obtain the Amplitude Spread Function (ASF), defined as the Fourier transform of *P*. In this case the basis functions ${C_k}$ in Eq. (2) must be replaced by their respective Fourier transforms (i.e. products of Bessel functions by angular frequencies), while keeping the same coefficients ${b_k}$. Since $PSF = {|{ASF} |^2}$ the appropriate basis functions would have the form [36]:

*N*= 14) and $m ={-} ({n_1} + {n_2}), - ({n_1} + {n_2}) + 2,\ldots + ({n_1} + {n_2})$. The latter equation suggests that, here too, there would be a high dimensionality and redundancy due to a lack of orthogonality. Although it is possible to build an orthogonal basis through SVD, it would prove difficult to develop given the lack of finite support. Since OTF and PSF are related by a linear Fourier transform, it seems therefore reasonable to use the OTF in practice.

The generalization to other geometries (e.g. off-axis image formation) can be accomplished by either repeating the process of substituting Zernike polynomials by a related orthogonal basis [37] or applying a geometrical transformation to the final basis functions, such as a numerical 2D affine transformation to pass from a circular to an elliptical support.

The second part of this work considered building a statistical model of the OTFs for a large population of human eyes. The OTF-specific modal representation permits us to obtain an important data compression by reducing the dimensionality of the individual OTFs through a threshold applied to the singular values. This threshold represents a trade-off between the lowest dimensionality and the highest fidelity of the representation. The subsequent PCA provides a population-specific model consisting of the mean plus a linear combination of orthogonal modes that account for the population variance. A second threshold is imposed on these eigenvalues according to the cumulative population variance, meaning that the variance of the statistical model is lower than the actual variance in the population. Each eigenvalue is the contribution of its corresponding eigenOTF to the total variance. If the statistical distribution is Gaussian, the model is exact, but only for the specific population analyzed. Thus one has to be careful when attempting to extrapolate to other groups or populations. One of the most interesting applications is to generate a virtually unlimited number of realizations (individual OTFs) to obtain a synthetic population with nearly the same values as the original [26], for a wide variety of visual processing modelling and simulations.

## Funding

Agentschap voor Innovatie door Wetenschap en Technologie (TBM 110684); Ministerio de Economía y Competitividad (PGC2018_095795_B_I00).

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **J. J. Rozema, D. E. M. Van Dyck, and M.-J. Tassignon, “Clinical comparison of 6 aberrometers. Part 2: Statistical comparison in a test group,” J. Cataract Refractive Surg. **32**(1), 33–44 (2006). [CrossRef]

**2. **L. N. Thibos, X. Hong, A. Bradley, and X. Cheng, “Statistical variation of aberration structure and image quality in a normal population of healthy eyes,” J. Opt. Soc. Am. A **19**(12), 2329–2348 (2002). [CrossRef]

**3. **T. O. Salmon and C. van de Pol, “Normal-eye Zernike coefficients and root-mean-square wavefront errors,” J. Cataract Refractive Surg. **32**(12), 2064–2074 (2006). [CrossRef]

**4. **N. Maeda, T. Fujikado, T. Kuroda, T. Mihashi, Y. Hirohara, K. Nishida, H. Watanabe, and Y. Tano, “Wavefront aberrations measured with Hartmann-Shack sensor in patients with keratoconus,” Ophthalmology **109**(11), 1996–2003 (2002). [CrossRef]

**5. **S. Shah, S. Naroo, S. Hosking, D. Gherghel, S. Mantry, S. Bannerjee, K. Pedwell, and H. S. Bains, “Nidek OPD-scan analysis of normal, keratoconic, and penetrating keratoplasty eyes,” J. Refract. Surg. **19**(2), S255–S259 (2003).

**6. **Y. Oie, N. Maeda, R. Kosaki, A. Suzaki, Y. Hirohara, T. Mihashi, Y. Hori, T. Inoue, K. Nishida, T. Fujikado, and Y. Tano, “Characteristics of ocular higher-order aberrations in patients with pellucid marginal corneal degeneration,” J. Cataract Refractive Surg. **34**(11), 1928–1934 (2008). [CrossRef]

**7. **E. Moreno-Barriuso, J. M. Lloves, S. Marcos, R. Navarro, L. Llorente, and S. Barbero, “Ocular aberrations before and after myopic corneal refractive surgery: LASIK-induced changes measured with laser ray tracing,” Invest. Ophthalmol. Visual Sci. **42**(6), 1396–1403 (2001).

**8. **T. Kuroda, T. Fujikado, N. Maeda, T. Oshika, Y. Hirohara, and T. Mihashi, “Wavefront analysis in eyes with nuclear or cortical cataract,” Am. J. Ophthalmol. **134**(1), 1–9 (2002). [CrossRef]

**9. **X. Cheng, L. N. Thibos, and A. Bradley, “Estimating visual quality from wavefront aberration measurements,” J. Refract. Surg. **19**(5), S579–S584 (2002).

**10. **J. D. Marsack, L. N. Thibos, and R. A. Applegate, “Metrics of optical quality derived from wave aberrations predict visual performance,” J. Vision **4**(4), 8–328 (2004). [CrossRef]

**11. **X. Cheng, A. Bradley, X. Hong, and L. N. Thibos, “Relationship between refractive error and monochromatic aberrations of the eye,” Optom. Vis. Sci. **80**(1), 43–49 (2003). [CrossRef]

**12. **L. N. Thibos, R. A. Applegate, J. T. Schwiegerling, and R. Webb, “Standards for reporting the optical aberrations of eyes,” J. Refract. Surg. **18**(5), S652–S660 (2002).

**13. **P. Artal, S. Marcos, D. R. Williams, and R. Navarro, “Odd aberrations and double-pass measurements of retinal image quality,” J. Opt. Soc. Am. A **12**(2), 195–201 (1995). [CrossRef]

**14. **R. Navarro and M. A. Losada, “Phase transfer and point-spread function of the human eye determined by a new asymmetric double-pass method,” J. Opt. Soc. Am. A **12**(11), 2385–2392 (1995). [CrossRef]

**15. **R. Navarro and M. A. Losada, “Shape of stars and optical quality of the human eye,” J. Opt. Soc. Am. A **14**(2), 353–359 (1997). [CrossRef]

**16. **E. C. Kintner and R. M. Sillitto, “A new analytic’ method for computing the optical transfer function,” Opt. Acta **23**(8), 607–619 (1976). [CrossRef]

**17. **J. Schwiegerling, “Relating wavefront error, apodization, and the optical transfer function: on-axis case,” J. Opt. Soc. Am. A **31**(11), 2476–2483 (2014). [CrossRef]

**18. **J. A. Díaz, “Relating wavefront error, apodization, and the optical transfer function: Comment,” Opt. Soc. Am A **33**(8), 1622–1625 (2016). [CrossRef]

**19. **C. Ferreira, J. L. López, R. Navarro, and E. Pérez-Sinusía, “Orthogonal basis for the optical transfer function,” Appl. Opt. **55**(34), 9688–9694 (2016). [CrossRef]

**20. **J. Porter, A. Guirao, I. G. Cox, and D. R. Williams, “Monochromatic aberrations of the human eye in a large population,” J. Opt. Soc. Am. A **18**(8), 1793–1803 (2001). [CrossRef]

**21. **J. J. Rozema, P. Rodriguez, R. Navarro, and C. Koppen, “Bigaussian wavefront model for normal and keratoconic eyes based on PCA,” Optom. Vis. Sci. **94**(6), 680–687 (2017). [CrossRef]

**22. **P. Rodríguez, R. Navarro, and J. J. Rozema, “Eigencorneas: application of principal component analysis to corneal topography,” Ophthalm. Physiol. Opt. **34**(6), 667–677 (2014). [CrossRef]

**23. **M. J. Jee, J. P. Blakeslee, M. Sirianni, A. R. Martel, R. L. White, and H. C. Ford, “Principal component analysis of the time-and position-dependent point-spread function of the advanced camera for surveys,” Publ. Astron. Soc. Pac. **119**(862), 1403–1419 (2007). [CrossRef]

**24. **P. Artal and R. Navarro, “Monochromatic modulation transfer function of the human eye for different pupil diameters: an analytical expression,” J. Opt. Soc. Am. A **11**(1), 246–249 (1994). [CrossRef]

**25. **A. B. Watson, “A formula for the mean human optical modulation transfer function as a function of pupil size,” J. Vision **13**(6), 18 (2013). [CrossRef]

**26. **J. J. Rozema, P. Rodriguez, R. Navarro, and M.-J. Tassignon, “SyntEyes: a higher-order statistical eye model for healthy eyes,” Invest. Ophthalmol. Visual Sci. **57**(2), 683–691 (2016). [CrossRef]

**27. **M. Born and E. Wolf, * Principles of Optics*(Pergamon Press, 1983).

**28. **ISO (International Organisation for Standardization), “ISO 24157:2008 - Ophthalmic optics and instruments: Reporting aberrations of the human eye,” (2008).

**29. **I. T. Jolliffe, * Principal Component Analysis*, 2nd edition (Springer-Verlag, 2002).

**30. **C. Leahy, C. Leroux, C. Dainty, and L. Diaz-Santana, “Temporal dynamics and statistical characteristics of the microfluctuations of accommodation: dependence on the mean accommodative effort,” Opt. Express **18**(3), 2668–2681 (2010). [CrossRef]

**31. **J. Arines, E. Pailos, P. Prado, and S. Bará, “The contribution of the fixational eye movements to the variability of the measured ocular aberration,” Ophthalm. Physiol. Opt. **29**(3), 281–287 (2009). [CrossRef]

**32. **S. Koh, N. Maeda, Y. Hirohara, T. Mihashi, S. Ninomiya, K. Bessho, and Y. Tano, “Serial measurements of higher-order aberrations after blinking in normal subjects,” Invest. Ophthalmol. Visual Sci. **47**(8), 3318–3324 (2006). [CrossRef]

**33. **P. Rodriguez, R. Navarro, L. González, and J. L. Hernández, “Accuracy and Reproducibility of Zywave, Tracey, and Experimental Aberrometers,” J. Refract. Surg. **20**(6), 810–817 (2004).

**34. **A. J. E. M. Janssen, “Extended Nijboer–Zernike approach for the computation of optical point-spread functions,” J. Opt. Soc. Am. A **19**(5), 849–857 (2002). [CrossRef]

**35. **S. van Haver and A. J. E. M. Janssen, “Advanced analytic treatment and efficient computation of the diffraction integrals in the extended Nijboer-Zernike theory,” J. Europ. Opt. Soc. Rap. Public. **8**, 13044 (2013). [CrossRef]

**36. **A. A. Ramos and A. L. Ariste, “Image reconstruction with analytical point spread functions,” Astron. Astrophys. **518**, A6 (2010). [CrossRef]

**37. **R. Navarro, J. L. López, J. A. Díaz, and E. Pérez Sinusía, “Generalization of Zernike polynomials for regular portions of circles and ellipses,” Opt. Express **22**(18), 21263–21279 (2014). [CrossRef]