1. 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 adjacent centroid positions, and can thus be considered as local extrapolation [11–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–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.

2. 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 _ref, of size 2N×1 (with at most 2n valid components). X _ref 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 _ref 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 ₀ (of size 18×1) is defined to initialize our algorithm. A first reduced set of Zernike coefficients Z ₀ (with typically five terms: tilts, defocus, and astigmatism) is then computed from δX ₀ using the least square method, which is classically described as a modal reconstruction of the wavefront [17–19]. The five components of Z ₀ 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 ₀ will only differ slightly from Z because of the measurement noise on the centroid positions δX ₀. From the approximated set of Zernike coefficients Z ₀, 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.

3. 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.

Table 1. Parameter of the aberrometer

View Table

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 ₀[i] the component of x that is closer to X _ref [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 ₀, 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 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 ₀ to form a vector X ₁ (of size 50×1). A set of 14 Zernike coefficients Z ₁ is estimated from the available centroid displacement δX ₁, 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.

 

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

Download Full Size | PDF

 

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.

Download Full Size | PDF

4. 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 z2,-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.

 

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

Download Full Size | PDF

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).