A non-stationary model for simulating the dynamics of ocular aberrations

The time-evolution of ocular aberrations has been the subject of many studies, but so far there has been little discussion involving the modelling of the underlying temporal statistics. This paper presents a non-stationary modelling approach based on a coloured-noise generator, which can be applied to ocular aberration dynamics. The model parameters are computed from measured ocular aberration data. A custom-built aberrometer based on a Shack-Hartmann sensor was used for measurement. We present simulations based on our modelling approach, and validate them through comparison to real data. This work could be useful in areas such as the testing of ophthalmic devices and the development of improved algorithms for laser refractive surgery. © 2010 Optical Society of America OCIS codes: (330.4060) Vision modeling; (330.0330) Vision, color, and visual optics. References and links 1. E. N. Bruce, Biomedical Signal Processing and Signal Modeling (Wiley Series in Telecommunications and Signal Processing, 2001). 2. C. Roberts, “Future challenges to aberration-free ablative procedures,” J. Refract. Surg. 16, 623-629 (2000). 3. A. Mira-Agudelo, L. Lundström, and P. Artal, “Temporal dynamics of ocular aberrations: Monocular vs binocular vision,” Ophthal. Physiol. Opt. 29, 256-263 (2009). 4. R. Navarro, E. Moreno-Barriuso, S. Bará, and T. Mancebo, “Phase plates for wave-aberration compensation in the human eye,” Opt. Lett. 25, 236-238 (2000). 5. S. O. Galetskiy, T. Yu. Cherezova, and A. V. Kudryashov, “Adaptive optics in ophthalmology: human eye wavefront generator,” Proc. SPIE 6849, 68490918 (2008). 6. D. R. Iskander, M. Collins, M. Morelande, and M. Zhu, “Analyzing the dynamic wavefront aberrations in the human eye,” IEEE Trans. Biomed. Eng 51, 1969-1980 (2004). 7. H. Hofer, P. Artal, and D. R. Williams, “Dynamics of the eyes wave aberration,” J. Opt. Soc. Am. A. 18, 497-506 (2001). 8. K. Hampson, E. Mallen, and J. C. Dainty, “Coherence function analysis of the higher-order aberrations of the human eye,” Opt. Lett. 31, 184-186 (2006). 9. K. Hampson, I. Munro, C. Paterson, and J. C. Dainty, “Weak correlation between the aberration dynamics of the human eye and the cardiopulmonary system,” J. Opt. Soc. Am. A. 22, 1241-1250 (2005). 10. M. Zhu, M. J. Collins, and D. R. Iskander, “Microfluctuations of wavefront aberrations of the eye,” Opthal. Physiol. Opt. 24, 562-571 (2004). 11. D. R. Iskander, M. R. Morelande, and M. J. Collins, “Estimating the dynamics of aberration components in the human eye,” In IEEE Signal Process. Workshop Statist. pages 241-244 (2001). 12. K. Billah and M. Shinozuka, “Numerical method for colored noise generation and its application to a bistable system,” Phys. Rev. A 42, 7492-7495 (1990). 13. N. Davies, L. Diaz-Santana, and D. Lara-Sucedo, “Repeatability of ocular wavefront measurement,” Optom. Vis. Sci. 80, 142-150 (2003). 14. L. Diaz-Santana, V. Guériaux, G. Arden, and S. Gruppetta, “Newmethodology to measure the dynamics of ocular wave front aberrations during small amplitude changes of accommodation,” Opt. Express 15, 5649-5663 (2007). 15. L. Diaz-Santana, C. Torti, I. Munro, P. Gasson, and C. Dainty, “Benefit of higher closed-loop bandwidths in ocular adaptive optics,” Opt. Express 11, 2597-2605 (2003). #131250 $15.00 USD Received 6 Jul 2010; revised 9 Sep 2010; accepted 22 Sep 2010; published 23 Sep 2010 (C) 2010 OSA 27 September 2010 / Vol. 18, No. 20 / OPTICS EXPRESS 21386 16. S. Gruppetta, F. Lacombe, and P. Puget, “Study of the dynamic aberrations of the human tear film,” Opt. Express 13, 7631-7636 (2005). 17. K. Y. Li and G. Yoon, “Changes in aberrations and retinal image quality due to tear film dynamics,” Opt. Express 14, 12552-12559 (2006). 18. L. Rankine, N. Stevenson, M. Mesbah, and B. Boashash, “A nonstationary model of newborn EEG,” IEEE Trans. Biomed. Eng 54, 19-28 (2007). 19. B. Boashash and M. Mesbah, “A time-frequency approach for newborn seizure detection,” IEEE Eng. Med. Biol. Mag. 20, 54-64 (2001). 20. M. Muma, D.R. Iskander, and M. J. Collins, “The role of cardiopulmonary signals in the dynamics of the eye’s wavefront aberrations,” IEEE Trans. Biomed. Eng. 57, 373-383 (2010). 21. N. Kasdin, “Discrete simulation of colored noise and stochastic processes and 1/f power law noise,” Proc. IEEE 83, 802-827 (1995). 22. B. M. Hill “A simple general approach to inference about the tail of a distribution,” Ann. Statist. 3, 1163–1174 (1975). 23. A. Clauset, C. R. Shalizi, and M. E. J. Newman, “Power-law distributions in empirical data,” arXiv:0706.1062v1, URL http://arxiv.org/abs/0706.1062v1 (2007). 24. N. R. Lomb, “Least-squares frequency analysis of unequally spaced data,” Astrophys. & Space Sci. 39, 447-462 (1975). 25. J. D. Scargle, “Studies in astronomical time series analysis ii. statistical aspects of spectral analysis of unevenly spaced data,” Astrophys. J. 263, 835-853 (1982). 26. P. Celka and P. Colditz, “Nonlinear nonstationary Wiener model of infant seizures,” IEEE Trans. Biomed. Eng. 49, 556-564 (2002). 27. L.R. Stark and D.A. Atchison, “Pupil size, mean accommodation response and the fluctuations of accommodation,” Opthal. Physiol. Opt. 17, 316-323 (1997). 28. 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, 26682681 (2010). 29. M. B. Priestley, Non-Linear and Non-Stationary Time Series Analysis (Academic Press, London, United Kingdom, 1988). 30. L. N. Thibos, A. Bradley, and X. Hong, “A statistical model of the aberration structure of normal, well-corrected eyes,” Ophthal. Physiol. Opt. 22, 427-433 (2002). 31. 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, 2329-2348 (2002). 32. L. Cohen, Time-Frequency Analysis, (Prentice-Hall, New York, 1995).


Introduction
Signal models of physiological processes are commonplace [1]; however there has been very little research regarding models of dynamic ocular aberrations other than defocus.Roberts [2] asserted that measuring and understanding the dynamic performance of the eye is critical to the future development of advanced customised refractive surgeries, which is in itself important for generating an overall improvement of refractive surgery results.Therefore, any improvements in the measurement and modelling of dynamic aberrations could be beneficial in areas such as the development of customised ablation algorithms.Such models could also be of interest in the study of the relationship between dynamics and visual performance [3], the development of customised contact lenses [4], and in the design of ocular adaptive optics systems with real-time correction [5].
Many dynamic properties of ocular aberrations have been established by previous authors [6,7].In particular, the time evolutions of some Zernike aberration modes have been shown to be correlated with each other [8], as well as with the cardiopulmonary system [9,10].Previous work has shown that the power spectral density (PSD) of aberration signals tends to follow a power-law relationship for frequencies below approximately 10 Hz, when the subject is not actively accommodating [3,7].The dynamics of Zernike modes have also been shown to exhibit features such as non-stationarity, harmonic distortion, and harmonic amplitude modulation [6,11].
Iskander et al. [11] proposed a methodology for modelling the dynamics of higher-order ocular aberrations using a parametric, amplitude-modulated, frequency-modulated (AM-FM) approach.This study focused in particular on the Zernike spherical aberration and coma coefficients, however measurements of just 5 s in duration (without taking into account the effects of subjects' blinking) were considered.Other work has demonstrated a potential application of a model of dynamic aberrations, by employing a deformable mirror to generate aberrations in real time [5].
In this work, we present a methodology for constructing power-law models of the temporal dynamics of ocular aberrations.The modelling approach is based on a numerical method for coloured-noise generation, developed by Billah and Shinozuka [12].We describe the connection between power-law relationships in the frequency domain and temporal stationarity, and how this may be exploited in the modelling of dynamic aberrations.We attempt to validate the performance of our model by comparing real and simulated data.

Data collection
All data was collected with a custom-built aberrometer based on a Shack-Hartmann wavefront sensor.The subject was a 26 year-old male.For each measurement, data was collected from the subject's dominant eye, whereas the non-dominant eye was covered with an eye patch.Eye dominance was assessed through an eye examination by a qualified optometrist prior to our experiments.Data was collected over 40 s, with an active pupil diameter of 5.4 mm.The sampling frequency was 100 Hz, giving a total of 4,000 samples per trial.The defocus term had the greatest magnitude, but there were significant contributions from astigmatism, coma, spherical aberration, and other higher-order terms.The system noise was white, with an RMS value of approximately 3 nm, as evaluated with an artificial eye (consisting of an 18 mm singlet coupled with a screen in the focal plane).
Analysis of dynamic ocular wavefront sensor measurements often takes the form of analysis of the Zernike coefficients as time series [3,[6][7][8][9][13][14][15][16][17].The temporal reconstruction of ocular wavefronts using Zernike polynomials can be represented as given in Eq. (1) [6]: where W (ρ, θ ;t) is the dynamic wave aberration function, m and n are the radial and azimuthal orders respectively, Z m n (ρ, θ ) is the Zernike polynomial expansion term, c m n (t) gives the corresponding Zernike coefficients, and ε(ρ, θ ;t) represents the modelling error.For our studies of aberration signals, we used an expansion of Zernike terms up to and including the 8 th radial order.
Measurements of dynamic ocular aberrations of more than a few seconds in duration are inevitably degraded by the subject's blinking.This leads to discontinuities and/or missing values in the measured signal [6,28].Though the blink interval is typically on the order of 250 ms [6], the transient components of the measured aberration signal associated with the blink can last for up to a second or more [28].For all measurements, we systematically removed any part of the signal that was related to the blink interval: tear film break-up and build-up, eye movement, and increased noise in the reconstruction of the Zernike modes of the wavefront due to a vignetted pupil.To incorporate these phenomena would significantly increase the complexity of any modelling approach.Therefore, as an approximation, we did not include them as part of the modelling objectives in this study.

Modelling approach
The modelling method is based on the assumption of a linear or piecewise linear power-law structure of the power spectral density (PSD) of the aberration signals.Synthesis of signals is based on a modification of an algorithm for the generation of power-law processes (or "coloured noise") [12].This algorithm has also been utilised as the basis for modelling approaches used in the field of EEG research [18], which deals with the study of signals that have many similarities to ocular signals in terms of their non-stationarity, power-law form of the spectral density, and dynamic range [19].
We initially assume that the power spectral density of the process to be modelled can be approximated by a power-law of the form given in Eq. ( 2): where f is the cyclic frequency, γ is the power-law exponent, and c is a constant.In general, for a finite-length, continuous-time signal x(t), the power spectral density can be written as shown in Eq. ( 3) [1]: where the angular brackets denote ensemble average.X( f ) is the Fourier transform of x(t), and can be written in the complex exponential form of Eq. ( 4): where φ ( f ) denotes the phase spectrum.To synthesise a time-domain signal from X( f ), we can specify the magnitude and phase spectra, and then take the inverse Fourier transform.We assume that the magnitude spectrum follows a power-law of slope γ/2.
In practice, the method described above can be implemented on a computer using an FFT algorithm.A detailed discussion of the discrete implementation of such an algorithm is given by Billah and Shinozuka [12].Rankine et al. [18] suggested that the "roughness" of real-world power spectra (as opposed to the "smooth" power-law form of the model) can be approximated by performing the synthesis for several different realisations of φ ( f ), and then summing together the resultant signals in the time-domain.

Parameter estimation
For estimating the slope value to fit to the PSD, the well-known Hill estimator was used [22].This method is simple to implement and gives an asymptotically normal and consistent estimate [23].For the phase spectrum, φ ( f ), one may use the phase spectrum of a measured signal as a "surrogate" phase.Alternatively, Billah and Shinozuka [12] as well as Kasdin [21] asserted that the phase spectrum of a power-law process can be assumed to be composed of random phase angles uniformly distributed over the interval [0, 2π].Artificial dynamic aberration signals can thus be generated by randomly choosing phase angles from a uniform distribution, to go with a PSD slope value arbitrarily chosen or empirically selected based on estimates from real data.

Validation of model
To justify the use of the model, we must assess its validity.Unfortunately, for non-stationary signal models, no well-established validation methodology exists [26].Rankine et al. [18] suggested using the Pearson product-moment correlation coefficient to validate simulated nonstationary data, in part due to its widespread familiarity.The application of this measure to processes that exhibit non-stationarity and non-linearity is however questionable.Muma et al. [20] Fig. 1.Dynamics of selected Zernike aberrations (left) and the associated PSD estimates (right).described a method of using time-frequency coherence estimates to obtain a normalised estimate of the coupling between two non-stationary signals.We are grateful to the authors for providing us with scripts for the implementation of this method.In the absence of a more rigorous validation approach, we use the time-frequency coherence between real and simulated signals as a means to assess model performance.We do not claim this to be a conclusive or quantitative method of validation (our use of the computed phase spectrum of measured data as a model parameter means that a significant level of coupling is to be expected), but merely an indicator that the simulated data does indeed resemble measured data to a reasonable degree.

Results
Figure 1 (left) shows the measured time-evolution of several Zernike coefficients; Fig. 1 (right) shows the corresponding PSD estimates (via the Lomb-Scargle periodogram [24,25]).As mentioned earlier, the temporal behaviour of ocular aberrations is in general non-stationary.This is illustrated by Fig. 2, which gives a time-frequency representation of a measured Zernike astigmatism (Z −2 2 ) signal.The time-frequency plot was obtained using the ZAM distribution, a member of Cohen's class of distribution functions [32].The plot shows low frequency content modulated in the < 0.5 Hz range over the full duration of the 40 s trial.There is significant power at approximately 1 Hz, which also varies in frequency and is likely to be related to cardiopulmonary activity [10,20].
To assess the validity of our model, we generated signals from the model using parameters obtained directly from real data.From a measured signal, we determined the measured phase spectrum and an estimate γ of the power-law slope, and used these as the model parameters φ ( f ) and γ respectively.We then synthesised a time-domain representation using these model parameters.For all measurements, our estimates of the power-law slope on a doubly logarithmic plot were greater than unity, and were similar to the PSD slope estimates of previous authors [3,7].Figure 3 compares the dynamics of measured Zernike spherical aberration to a simulated signal generated using our modelling method, while Fig. 4 shows the corresponding PSD estimates.Figure 5 is a plot of the measured phase spectrum φ ( f ) that is used as the surrogate phase.
Figure 6 shows the time-frequency coherence between the real and modelled signals of     that of the original measured signal from which the model parameters were estimated.

Performance of models
Rankine et al. [18] suggested that greater flexibility in the modelling of non-stationary data can be achieved by allowing the estimated power-law exponent γ to be time-varying, i.e. γ = γ(t).We have observed that γ can vary significantly over a single epoch of data.However, by allowing γ to vary over time in our model, we observed only marginal improvement in the model performance at best.This could be partly due to the difficulties in validating model performance.However, we feel that the potential benefits in allowing γ to be time-varying must be weighed against some drawbacks.Firstly, the time-scales upon which γ varies significantly are difficult to determine.Secondly, by attempting to estimate γ over shorter time-scales, one must take into account the additional bias in the estimate of γ.

Two-slope model
The shape of the PSD varies with the level of accommodative effort [27,28].When the eye actively accommodates, the shape of the PSD in this range often exhibits two distinct slopes with a "breaking frequency" between them [28].This effect is not observed when the eye is in its unaccommodated state, such as in Fig. 1.For the cases where the estimated PSD is better fitted by two slopes rather than one, the model can be adapted accordingly.A bilinear fit can be performed by adopting a regression approach and imposing a constraint on the regression function.In this case, γ = γ( f ), with γ( f ) given by Eq. ( 5): where γ 1 and γ 2 are constants representing estimates of the PSD slope for lower and higher frequencies respectively.The two slope regions are separated by the breaking frequency f br .

Relationship Between Power-Law Slope and Stationarity
Kasdin [21] stipulated that for power-law processes, the stationarity of the process is related to the slope γ of the PSD.In particular, the slope of the PSD at low frequencies is indicative of whether the process is stationary or not.This can be summarised as follows: • A value of γ > 1 at low frequencies implies inherent non-stationarity of the underlying process.
• Conversely, a value of 0 ≤ γ ≤ 1 at these lower frequencies implies that the process is stationary.
We illustrate this notion in Fig. 7, with some examples of signals synthesised using the twoslope model.The higher frequency slope was fixed to γ 2 = 2.0.We then generated signals using five different values of γ 1 .In each case, an identical random phase spectrum (composed of uniform random variates distributed over the interval [0, 2π]) was used.To employ a random phase spectrum (rather than a measured one), is potentially more useful as a general approach to modelling the temporal behaviour of the eye's aberrations.Both methods can be used to generate either stationary or non-stationary signals with frequency content shaped by γ( f ), as required.As the value of γ 1 exceeds unity, the effect on the stationarity of the signals becomes clearly visible (for example, the mean of the process wanders significantly).We can consider the value of γ 1 to relate to the "statistical stability" of the signal, rather than to stationarity in a more rigorous sense.It is useful in this context to think of stationarity in a relative sense, as truly stationary signals are mathematical constructs that do no not typically occur in real-world data [29].

Application of models
The model for the dynamics of ocular aberrations was originally conceived as a way to extend static models of aberrations over a large population, such as the model developed by Thibos et al. [30].If the dynamic model parameters could be similarly estimated over a large population, the authors' idea of creating a database of "virtual eyes" could be extended to include the dynamics of aberrations.A notable feature of the aberration signals is their dependence on the mean accommodative effort [28].The modelling technique described in this work could similarly be used to model dynamic accommodation, as measures of accommodation can be computed from Zernike terms [31].A model of dynamic accommodation could be adapted to include the state-dependence by modelling the relationship between the power-law slope values and the accommodative effort.
It should be noted that in our modelling efforts we have not attempted to compensate for the movements of the eye during measurement.Given that the ultimate aim of this work is to simulate the temporal behaviour of the eye's intrinsic aberrations (rather than those measured by the wavefront sensor), it follows that these movements should ideally be compensated for, e.g., by using some active tracking of the pupil.
Another property that we have not directly addressed in our modelling approach is the nonlinearity of aberration and accommodation signals.We previously discussed the presence of harmonic distortion in our measured signals, which suggests the presence of non-linearity.This is not surprising for a physiological process.Non-linearity has been directly addressed in EEG signal models [26], and we feel that proper study of non-linearity in ocular signals is also warranted (e.g., using bispectrum or bicoherence methods) for the improvement of the modelling effort.

Fig. 1 .Fig. 2 .
Figure 1 (left) shows the measured time-evolution of several Zernike coefficients; Fig. 1 (right)shows the corresponding PSD estimates (via the Lomb-Scargle periodogram[24,25]).As mentioned earlier, the temporal behaviour of ocular aberrations is in general non-stationary.This is illustrated by Fig.2, which gives a time-frequency representation of a measured Zernike astigmatism (Z −2 2 ) signal.The time-frequency plot was obtained using the ZAM distribution, a member of Cohen's class of distribution functions[32].The plot shows low frequency content modulated in the < 0.5 Hz range over the full duration of the 40 s trial.There is significant power at approximately 1 Hz, which also varies in frequency and is likely to be related to cardiopulmonary activity[10,20].To assess the validity of our model, we generated signals from the model using parameters obtained directly from real data.From a measured signal, we determined the measured phase spectrum and an estimate γ of the power-law slope, and used these as the model parameters φ ( f ) and γ respectively.We then synthesised a time-domain representation using these model parameters.For all measurements, our estimates of the power-law slope on a doubly logarithmic plot were greater than unity, and were similar to the PSD slope estimates of previous authors[3,7].Figure3compares the dynamics of measured Zernike spherical aberration to a simulated signal generated using our modelling method, while Fig.4shows the corresponding PSD estimates.Figure5is a plot of the measured phase spectrum φ ( f ) that is used as the surrogate phase.Figure6shows the time-frequency coherence between the real and modelled signals of Fig.1.The value of the time-frequency coherence function is considerably close to unity for most values of t and ω, which indicates that the simulated signal has spectral content that is close to

Fig. 3 .
Fig. 3. Comparison of measured and simulated Zernike spherical aberration signals in the time-domain.

Fig. 4 .
Fig. 4. Comparison of measured and simulated Zernike spherical aberration signals via their estimated power spectral densities.

Fig. 5 .
Fig.5.Discrete Fourier phase spectrum computed from measured data and used as the "surrogate" phase to produce the simulated aberration signal of Fig.3.This Fourier phase spectrum representation was computed using an interpolation algorithm to replace missing data points.

Fig. 7 .
Fig. 7. Illustration of the relationship between statistical stability about a mean and the PSD slope at low frequencies.All signals were generated using the two-slope model, with an identical phase spectrum generated from a uniform distribution.