A simple and robust method to extend the dynamic range of an aberrometer

We present an algorithm to extend significantly the dynamic range of a Shack-Hartmann wavefront sensor. With this method, the recorded Shack-Hartmann spots are not constrained to stay in the field of view of their lenslet. The proposed algorithm is computationally effective, robust to a high level of noise on the measured centroid positions and also to missing centroid values. The principle is closely related to the description of wavefronts using Zernike polynomials, which makes optimization for a given sensor and application achievable thanks to numerical simulation. These features make it useful for the measurements of highly aberrated eyes. © 2009 Optical Society of America OCIS codes: (120.5050) Phase measurement; (120.3940) Metrology. References and links 1. J.Primot, “Theoretical description of Shack-Hartmann wave-front sensor,” Opt. Commun. 222, 81-92 (2003). 2. A.Burvall, E. Daly, S. R. Chamot, and C. Dainty, “Linearity of the pyramid wavefront sensor,” Opt. Express 14, 11925–11934 (2006). 3. D.Malacara, Optical Shop Testing (Wiley-Interscience,2007). 4. R. Navarro and E. Moreno-Barriuso, “Laser ray-tracing method for optical testing,” Opt. Lett. 24, 951-953 (1999). 5. N. Lindlein, J. Pfund, and J.Schwider, “Algorithm for expanding the dynamic range of a Shack-Hartmann sensor by using a spatial light modulator array,” Opt. Eng. 40, 837–840 (2001). 6. G.Yoon, S. Pantanelli, and L. J. Nagy, “Large-dynamic-range Shack-Hartmann wavefront sensor for highly aberrated eyes,” J. Biomed. Opt. 11, 030502–030504 (2006). 7. L. E. Schmutz and B. M. Levine, “Hartmann sensors detect optical fabrication errors,” Laser Focus World 32, (1996). 8. X. Levecq and S. Bucourt, ”Method and device for analysing a highly dynamic wavefront,” US patent 6,750,957B1 (2004) 9. N. Lindlein, J. Pfund, and J. Schwider, “Expansion of the dynamic range of a Shack-Hartmann sensor by using astigmatic microlenses,” Opt. Eng. 39, 2220-2225 (2000). 10. N. Lindlein and J. Pfund, “Experimental results for expanding the dynamic range of a Shack-Hartmann sensor using astigmatic microlenses,” Opt. Eng. 41, 529-533 (2002). 11. J. Pfund, N. Lindlein, and J.Schwider, “Dynamic range expansion of a Shack-Hartmann sensor by use of a modified unwrapping algorithm,” Opt. Lett. 23, 995–997 (1998). 12. S. Groening, B. Sick, K. Donner, J. Pfund, N. Lindlein, and J. Schwider, “Wave-Front Reconstruction with a Shack-Hartmann Sensor with an Iterative Spline Fitting Method,” Appl. Opt. 39, 561–567 (2000). 13. D. G. Smith and J. E. Greivenkamp Schwider, “Generalized method for sorting Shack-Hartmann spot patterns using local similarity,” Appl. Opt. 47, 4548–4554 (2008). 14. L. Lundstrm and P. Unsbo, “Unwrapping Hartmann-Shack images from highly aberrated eyes using an iterative B-spline based extrapolation method,” Optom. Vis. Sci. 81, 383-388 (2004). 15. N. Lindlein and J. Pfund, “Phase retrieval by demodulation of a Hartmann-Shack sensor,” Opt. Commun. 215, 285-288 (2003). 16. M. C. Roggemann and T. J. Schulz, “Algorithm to increase the largest aberration that can be reconstructed from Hartmann sensor measurements,” Appl. Opt. 37, 4321–4329(1998). 17. R. K. Tyson, Principles of Adaptive Optics (Academic Press, 1998). #115839 $15.00 USD Received 17 Aug 2009; revised 2 Oct 2009; accepted 3 Oct 2009; published 7 Oct 2009 (C) 2009 OSA 12 October 2009 / Vol. 17, No. 21 / OPTICS EXPRESS 19055 18. L. Diaz-Santana, G. Walker, and S.Bará, “Sampling geometries for ocular aberrometry: A model for evaluation of performance,” J. Opt. Soc. Am. A 13, 8801-8818 (2005). 19. S. Bará, “Characteristic functions of Hartmann-Shack wavefront sensors and laser-ray-tracing aberrometers,” J. Opt. Soc. Am. A 24, 3700–3707 (2007). 20. L. N. Thibos, R. A. Applegate, J. T. Schwiegerling, and R. Webb, “Standards for reporting optical aberrations of eyes,” J. Ref. Surg. 18, 652-660 (2000). 21. N.Maeda, T.Fujikado, T.Kuroda, T.Mihashi, Y.Hirohara, and K.Nishida, “Wavefront aberrations measured with Hartmann-Shack sensor in patients with keratoconus,” Ophthalmology 109, 1996-2003 (2002). 22. S.Pantanelli, S.MacRae, T.M.Jeong, G.Yoon, “Characterizing the wave aberration in eyes with keratoconus or penetrating keratoplasty using a high-dynamic range wavefront sensor,” Ophthalmology 114, 2013-2021 (2007).


Introduction
The Shack-Hartmann wavefront sensor has some fundamental limitations in terms of spatial sampling [1], and is thus rarely the most suitable technique to measure wavefronts with high spatial frequencies.The dynamic range of the Shack-Hartmann is however not limited by any physical principle, which is not the case for some other wavefront sensors.The pyramid wavefront sensor, for example, has a limited range of linearity, which is defined by the amplitude of modulation of the focussed beam [2].Any measurement system based on interferometry has a dynamic range limited by the coherence properties of the light source [3].The dynamic range of a Shack-Hartmann wavefront sensor is only limited by the ability of the system to process the recorded data.For highly distorted wavefronts, two major difficulties occur for the processing of the Shack-Hartmann data.
The first problem occurs when two adjacent spots are so close that they cannot be processed independently.The "crossover" of adjacent spots happens when the curvature of the wavefront is locally too large.The main method to overcome this limitation of the Shack-Hartmann wavefront sensor is to record the spots sequentially, which is what a laser ray tracing wavefront sensor does [4].This can be done using a spatial light modulator that switches on and off the lenslets, as it was implemented by Lindlein et al. [5].Yoon et al. [6] successfully implemented a similar method using a translatable obstruction mask for the measurement of ocular aberrations in a population of keratoconus eyes.
A second difficulty occurs when the Shack-Hartmann spots leave the field of view of their lenslets.This effect is considered to be the first limitation in the operation of a Shack-Hartmann wavefront sensor, because the crossover of adjacent spots usually happens for a higher amount of aberration.This limitation of the Shack-Hartmann can be overcome with methods that are robust to large displacements of the spots.A simple hardware solution is to measure the centroid positions for different longitudinal positions of the detector [7].Some algorithmic solutions have also been suggested in the literature.Extending the range of a Shack-Hartmann with an algorithmic solution permits the measurement of highly aberrated wavefronts instantaneously, which is one of the key features of a Shack-Hartmann wavefront sensor.For most algorithmic approaches, the problem consists in finding the lenslet that corresponds to each measured centroid position.Various methods have been suggested in the literature.
One approach is to manufacture a lenslet array that has a calibrated periodic change in the structure of the lenslet array (for example the pitch of the lenslet), and to identify this structure in every measured map of centroid positions [8].A similar method was also introduced by Lindlein et al. [9,10], who suggested using a lattice of astigmatic lenslets with different (and calibrated) principal axes, and to identify the corresponding direction in the pattern of each Shack-Hartmann spot.
A powerful idea to retrieve the corresponding lenslet of each Shack-Hartmann spot is to apply an iterative extrapolation of the measured centroid positions, using some previously indexed centroid positions.Previously proposed methods to implement this idea rely on the use of adja-cent centroid positions, and can thus be considered as local extrapolation [11][12][13].One of these methods, which relies on a spline extrapolation [12], has been successfully implemented for the measurement of ocular aberrations at large field angles [14].
An interesting alternative to measure a highly aberrated wavefront is to directly extract the gradient of the wavefront from the raw CCD frame, using a Fourier demodulation technique [1,15].This approach is simple and has a low computation time.Roggemann et al. [16] also introduced a method to extend the dynamic range of a Shack-Hartmann wavefront sensor using the raw CCD frame and an additional image of the point spread function.
We propose in this paper a method based on the iterative extrapolation of the measured centroid positions.The main difference with the above methods [11][12][13][14] is that the extrapolation required to gradually analyze the whole pupil is performed on the wavefront itself, using Zernike polynomials.Doing this permits an understanding of our method in terms of Zernike coefficients of the incoming wavefront.The smoothing effect of Zernike coefficients makes this method robust to common features of the data recorded by the aberrometer, such as missing centroid values and noisy measurements of the centroid positions.We introduce the principle of our method, and discuss its implementation for the measurements of ocular aberrations.

Principle of the algorithm
We assume that the raw frame recorded by the Shack-Hartmann wavefront sensor has been accurately processed by an algorithm that yields a (2n × 1) set x of estimated centroid positions that cannot be easily associated with their reference counterpart (n ≤ N, N being the number of lenslets over the measured pupil).The problem we propose to solve is to retrieve the lenslet associated to each centroid position.The output of the algorithm is thus a set of centroid displacements δ X = X − X re f , of size 2N × 1 (with at most 2n valid components).X re f is the set of centroid positions measured with a reference wavefront, and X the set of sorted centroid positions.Before invoking the algorithm presented in this paper, each component of X re f is associated with a position in the pupil plane.
Assuming that a certain number of spots remain in the field of view of their lenslet (typically a 3 × 3 grid), a vector of centroid displacements δ X 0 (of size 18 × 1) is defined to initialize our algorithm.A first reduced set of Zernike coefficients Z 0 (with typically five terms: tilts, defocus, and astigmatism) is then computed from δ X 0 using the least square method, which is classically described as a modal reconstruction of the wavefront [17][18][19].The five components of Z 0 can only be considered as a rough approximation to their true values Z, because they were estimated using 9 centroid displacements only.In the case where these modes (tilts, defocus, and astigmatism) are actually the only significant modes of the true wavefront, Z 0 will only differ slightly from Z because of the measurement noise on the centroid positions δ X 0 .From the approximated set of Zernike coefficients Z 0 , it is straight-forward to compute an approximation X of all the centroid positions X of the sensor, using the geometry matrix A [17]: We recall that each component of A is the shift of a given Shack-Hartmann spot (row index) induced by a given Zernike aberration (column index).
For each lenslet i, we then search among the components of x for the computed centroid position that is the nearest to X[i].This allows us to form the set of centroid positions X (of size 2N × 1, but with invalid data that mark the missing centroid positions) that corresponds to the Shack-Hartmann spots that have been associated to their reference counterpart.A simple way to identify a Shack-Hartmann spot is to use a threshold value d (in pixels) for the maximum distance allowed between the approximated centroid position X[i] and the nearest component of measured centroid positions x.

Implementation
We briefly describe some details of the algorithm we implemented with an aberrometer whose specification is given in Table 1.The dynamic range of this aberrometer in terms of curvature of the wavefront (crossover of adjacent spots) is relatively large.For a 10 dioptres myopic eye, adjacent spots are still separated by 16 pixels (instead of the 18.5 pixels pitch), but some of the spots are shifted by more (24 pixels) than the dimension of a lenslet.Without any dedicated algorithm, the maximum power that can be measured with this aberrometer over a 6 mm pupil is approximatively 2 dioptres.The aberrometer uses a very narrow probing beam, which is not corrected for the refractive error of the subject and parallel to the axis of the instrument.Due to reciprocity of propagation of light, the local tilt of the measured wavefront is zero at the position where the probing beam enters the eye, when the probing beam is parallel to the optical axis of the wavefront sensor.We thus initialize the algorithm with a 3 × 3 grid of lenslet centred on the probing beam.For each lenslet i of this 3 × 3 grid, our algorithm attributes to X 0 [i] the component of x that is closer to X re f [i] than a distance equal to half the pitch of a lenslet.
Once the initialization is done, we compute the approximated set of Zernike coefficients Z 0 , and then a first set of approximated centroid positions X.To improve robustness, the vector δ X is not reconstructed at once.In many cases, the vector X is a good approximation of X, and all centroid positions can be retrieved at once.It is however safer to do the extrapolation in several steps, especially if there is a significant amount of coma in the measured wavefront.Figure 1 shows how we typically extrapolate over the whole pupil: each number corresponds 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14   to the number of steps used before attributing a centroid position to a given lenslet.From the first calculation of X, only 16 centroid positions (marked by a "1" in Figure 1) are identified and concatenated to X 0 to form a vector X 1 (of size 50 × 1).A set of 14 Zernike coefficients Z 1 is estimated from the available centroid displacement δ X 1 , and the process is repeated 19 times until the whole pupil has been analyzed.The number of Zernike terms we use for this algorithm is 14 for any further iteration, because we assume that a priori only these modes will be significantly large.Using a normal computer with a 2 GHz processor and a Matlab implementation of the algorithm, the full process takes approximatively 0.4 seconds and is thus only useable offline.For a real time computation of ocular aberration, it is not necessary to invoke the algorithm at each frame.Each centroid position can be tracked from the previous frame.

Results
Figure 2(a) shows the spot displacements reconstructed from the double-pass measurement of an ophthalmic cylindrical lens.The lens was placed in the measurement plane of our aberrometer, just in front of an artificial eye (a 40 mm doublet that has an opaque screen in its focal plane).The position of the probing beam in the measurement plane is marked by a black circle.
The measured astigmatism is (z 2,−2 = −5.1 µm, z 2,2 = 0.5 µm), using the convention adopted by the OSA [20].Figure 2(b) shows the spot displacements reconstructed from a simulation of the centroid positions obtained with a similar astigmatic wavefront z 2,−2 = −5.5 µm.The position of the centroid positions of the spots are computed using the geometry matrix of the Shack-Hartmann, and permuted before invoking our algorithm, in order to simulate the detection of the centroid positions with software windows that are not bound to the lenslet array.The global tilt of the measured wavefront is adjusted so that the spot displacement of the lenslet centred on the probing beam is zero, in order to satisfy the principle of reciprocity of the propagation of light.For this simulation and the one shown in Figure 3, we consider our algorithm as successful if all centroid positions are correctly paired to their reference counterpart.1, with and without our algorithm, for combinations of astigmatism z 2,−2 (vertical axis) and various higher order aberrations (horizontal axis).
The dynamic range of our algorithm, which is suitable for double-pass measurements, is primarily limited by the successful initialisation of the 3 × 3 grid.If the wavefront is locally too aberrated, the assumption that the spots of this grid remain in the field of view of their lenslet might not hold.(But the spot in the middle of this grid will remain unshifted, if the probing beam is correctly centred on the grid.)We show in Figure 3 the extension of the dynamic range for combinations of astigmatism z 2,−2 (vertical axis of each graph) and higher order aberrations (horizontal axis): trefoil z 3,−3 (a), coma z 3,−1 (b), secondary astigmatism z 4,−2 (c), and spherical aberration z 4,0 (d).The black region corresponds to the combinations of aberrations for which each spot remains in the field of view of its lenslet.The gray region corresponds to the region for which all the spots are correctly assigned to their reference counterpart with our algorithm.It is important to mention that these simulations are obtained without performing the full detection of the CCD frame with floating software windows, and that the magnitude of aberration in Figure 3 is much larger than that obtained in real eyes.The centroid positions are computed, as in Figure 2(b), with the geometry matrix of the sensor.According to the simulations shown in Figure 3, the dynamic range of our aberrometer is primarily limited by the crossover of adjacent spots, which is not taken into account in our simulations.For example, with a spherical aberration z 4,0 = 4 µm (and no astigmatism), some adjacent spots are separated by 6 pixels only.Assuming that their centroid positions are processed accurately, our algorithm retrieves their reference positions.
These results show that our algorithm is suitable for the measurements of patients with keratoconus, who have on average around 2.5 µm of root mean square higher order aberrations over a 6 mm pupil [21,22].A Matlab code and demonstration of the algorithm is available on our website (http://optics.nuigalway.ie/people/charlie).

Fig. 1 .
Fig. 1.Steps for the extrapolation of the centroid positions across the measured pupil.

Fig. 2 .
Fig.2.Successful pairing of the Shack-Hartmann spots (o) with their reference counterparts (x), for an astigmatic wavefront and a probing beam located 2.25 mm superior to the pupil centre.Left: double-pass measurement with an ophthalmic lens.Right: simulation.

Fig. 3 .
Fig. 3. Dynamic range of the aberrometer of Table1, with and without our algorithm, for combinations of astigmatism z 2,−2 (vertical axis) and various higher order aberrations (horizontal axis).