1.

Introduction

Hartmann-Shack aberrometry is one of the most common methods used to characterize the wavefront aberrations of an optical system,¹ most especially the aberrations of the human eye.² In its most common implementation, an array of lenslets is used to form an array of images of a point source that has passed through some optical media. Each subimage, or “spot,” is deviated away from the axis of its lenslet by an amount proportional to the average wavefront slope across the lenslet. Measuring these deviations across the pupil enables the wavefront shape to be reconstructed and/or an appropriate adaptive optics correction to be determined.¹

Locating the position of each spot is relatively trivial using standard image processing techniques. The more challenging task is to ensure that each spot is assigned to the correct lenslet, a process referred to as “sorting” the Hartmann-Shack pattern. The easiest way to accomplish this is to assume that each spot image is formed on the detector plane within the bounds of its lenslet. The problem is then reduced to a series of small search windows of fixed location and size. We refer to this as the “conventional” algorithm. With this algorithm, when the local wavefront slope is too steep, a spot will appear outside the bounds of its lenslet and therefore be ascribed to the wrong lenslet. This causes errors in the wavefront fit and in any adaptive optics correction.

The range of this conventional algorithm is sufficient for the majority of normal human eyes when measuring aberrations close to the optical axis. However, the algorithm can fail when measurements are sufficiently off-axis, when refractive error is moderate to high, or when pathology is present that distorts the refractive components of the eye.

Several hardware changes to Hartmann-Shack systems have been proposed to circumvent the limitations of the conventional approach. Larger magnification of the pupil can be used to reduce wavefront slope and hence spot deviation, but this requires a wider detector.² Lenslet focal length can be shortened to again provide less deviation at the detector plane, but this decreases accuracy due to pixelation effects.² Lenslet size can be made larger, but at the expense of reduced lateral spatial resolution of the wavefront.² Yoon suggested a rapidly translating plate that could block certain lenslets in alternating frames, although this sacrifices temporal resolution.³ Lindlein and Pfund suggested an array of astigmatic lenslets with different orientations to uniquely identify each spot in a local area.⁴ Ares described the use of an array of astigmatic lenslets together with a software algorithm that was able to trace the resultant line foci outside the bounds of the conventional search area.⁵ For aberrometry on the human eye, appropriate trial lenses are often put into place. This can inject unwanted reflections and is time consuming in that a refraction must be performed first, and it does not help with imaging the peripheral fundus where refraction changes markedly,⁶ or for higher order effects such as in patients with corneal disease.

Software-based approaches can also increase the range, while avoiding alterations to existing systems and minimizing the cost of system construction. For example, Pfund described a sorting algorithm that was inspired by a common approach used in unwrapping phase images in interferometry.⁷ Groening ⁸ and Lundström and Unsbro⁹ described more complex algorithms where a spline surface was fit to central spots, and then iteratively extrapolated to predict the locations of subsequent spots. Leroux and Dainty used a similar approach but fit Zernike polynomials to the spot displacement data instead of a spline, producing a smoothing effect that suppresses noise.¹⁰ Lee ⁴ and Smith and Greivenkamp¹² each described algorithms in which spots are located in a predetermined order in such a way as to reduce ambiguity as to lenslet-spot correspondence.

These software methods generally reported dramatic improvements in the range of aberrometry, but it is difficult to compare algorithms directly due to differences in system parameters, aberrations considered, and criteria for failure. In particular, no range data at all appeared with the Smith and Greivenkamp¹² algorithm. The major aim of this review is to facilitate the direct comparison of several of these algorithms by providing further range data obtained under consistent conditions. We also explore the robustness of the considered algorithms in the context of ocular aberrometry by considering the effects of poorly sensed spots, decentration of the pupil (in the absence of a pupil monitor), and decentration of the input beam.

2.

Methods—The Algorithms

Three of the preceding algorithms were selected to reflect the most modern tactics that provide large increases in range over the conventional approach, and to cover a range of complexity in terms of processing required.^{9, 10, 12} For each algorithm, it is assumed that the spot positions are already known via some image processing method, and that the only remaining problem is the assignment of the spots to their corresponding lenslets. Figure 1 depicts the order in which spots are located for each algorithm. For a more thorough treatment of a given algorithm, the reader should consult the original paper.

Fig. 1

Order in which lenslets of a square lenslet array are unwrapped for each of the algorithms explored here. The zero entries indicate the initial location of central spots before the main part of the algorithm takes effect. Subsequent numbers indicate iterations of each algorithm. The quoted terms are used to reference each algorithm in the text.

It is noteworthy that each of the algorithms explored here approach the spot allocation process in a generally similar way. An initial assumption is made in each case that a small number of spots nearest to the pupil center correspond to their directly overlying lenslets. This is a valid assumption since aberrations are smallest at the pupil center. From this anchor point, each algorithm expands iteratively outward from the pupil center. On a given iteration only spots that are adjacent to previously allocated spots are sorted, which enables the algorithms to proceed with minimum uncertainty since the wavefront changes minimally from one lenslet to the next.

3.

Methods—Analysis

To initially compare each of the preceding algorithms, custom MATLAB code was used to simulate double-pass Hartmann-Shack images obtained from a human eye. Monochromatic aberrations were simulated at 830 nm using a phase plate located in the pupil plane, with the retina located 20.2 mm behind this. The retinal space refractive index was 1.33. Fourier optics was used to propagate a 1.0-mm-diam entry beam to form a realistically blurred point source on the retina. Light from a point on this retinal object was then propagated back out of the “eye” with an exit pupil diameter of 6.0 mm. The pupil plane wavefront was sectioned into 0.4-mm-diam squares to simulate a 0.4-mm-pitch square lenslet array. For each section, the point spread function was computed for a lenslet with a 24.0-mm focal length, and convolved with the retinal spot to simulate the double-pass image of the spot on the detector.

Spot centroids were located by standard image processing functions in MATLAB in two steps. First, a binary image was generated by thresholding at intensity levels of 20% and the MATLAB function “bwlabel” was used to label the centers of the resulting binary spots. Second, these approximate labels were supplied as arguments to the MATLAB “regionprops” function, along with the original gray-scale image, to determine the center of mass for each spot.

To quantify the range for a given aberration, the aberration was varied in magnitude with a trial-and-error approach until the cut-off point was reached (to within 0.05 μm) at which spots began to be incorrectly assigned to a lenslet. Only negative Zernike coefficients were considered—e.g., the results for defocus represent a myopic eye. Negative units were chosen because the range with the improved algorithms is much larger for hyperopic defocus; since spots do not overlap here it is actually the size of the detector array that limits the range.

As well as determining range for a well-centered eye, we modeled the effects on the achievable range of some common sources of errors in ocular aberrometry. First, without a separate pupil monitor, there may be as much as a 0.2-mm error in lateral pupil alignment for our 0.4-mm-pitch lenslet array; we considered the effects of a 0.2-mm displacement on the range. Second, the incident beam is usually displaced from the optical axis in ocular aberrometry to avoid the corneal reflection. This causes the least-displaced Hartmann-Shack spot to lie away from the pupil center, which can introduce errors if not taken into account.¹⁰ To simulate this effect we modeled a 1.2-mm error in the deviation of the beam from the optical axis. Finally, when the SNR is low, some spots may not be adequately detected. We randomly excluded 10 of the 177 spots in the pupil from the analysis and measured the range achievable (the same 10 spots were excluded for each analysis).

Following the computational modeling just described, we employed a real adaptive optics system to offer a more “real-world” test of each algorithm. The system consisted of a Hartmann-Shack wavefront sensor with identical parameters to that modeled here, in combination with a Mirao 52-d deformable mirror from Imagine Eyes. A model eye was used, consisting of a +20 D lens positioned 3 mm in front of an adjustable aperture (the pupil). A piece of white paper acted as the retina.

To sample the large number of wavefront shapes achievable with the deformable mirror, we systematically applied random voltage signals to the 52 actuators. The random commands were generated on the range –1 to +1 (the maximum signal possible to send to the mirror), scaled by some scaling factor. Scaling factors from 0.05 to 0.50 were used, in increments of 0.05. This means that in the highest amplitude condition it was possible for neighboring actuators to have a difference in voltage of 1.0. For each scaling factor, 100 random sets were generated and applied to the mirror actuators, and the resulting Hartmann-Shack image was recorded. We therefore recorded 1000 separate images resulting from the mirror undergoing various shapes and amplitudes of deformation. Each image was processed in the same way as for the preceding modeling (a higher threshold was used due to the presence of noise in the real images) and analyzed with each of the algorithms. In the case of the Zernike algorithm, it is critical to know the location of the input beam; to measure this we placed a white card in the pupil plane and estimated the centre of the ∼1.2-mm-diam spot that was produced on the wavefront sensor CCD.

Since in this real system case we did not have definitive knowledge of spot-lenslet correspondence, we relied on the internal error checking of each algorithm to reject spots that were probably not correctly assigned. This is a possible source of error, although each of the algorithms used rejects spots in a similar and generally sensible way (primarily by comparing predicted positions to found positions and rejecting those that are greater than a certain number of pixels away). Relying on this internal error checking, we arbitrarily defined failure as 5 or more spots being missed of the ∼280 spots across the 7.5-mm pupil. This criterion is slightly more lax than the single-spot criterion used in the preceding computational modeling, but it seemed necessary to account for occasional poorly sensed spots.

4.

Results and Discussion

The results for the computational modeling are graphed in Fig. 2 under the various conditions already described for defocus, horizontal-vertical astigmatism, vertical coma, and spherical aberration. Data for the various aberration types were stacked to aid visualization. All units are normalized Zernike coefficients in micrometers, in accordance with reporting standards for vision science,¹⁶ and were fit over a 6.0-mm pupil. Coefficients are all negative in sign (e.g., the defocus term represents a myopic eye).

Fig. 2

Range for each Hartmann-Shack sorting algorithm for various types of aberration, expressed in units of (negative) normalized Zernike coefficients for a 6.0-mm pupil. Astigmatism is horizontal-vertical and coma is vertical. In the “decentered pupil” case, there was a 0.2-mm error in assessment of the pupil center; in the “decentered beam” case there was an uncompensated 1.2-mm decentration of the input beam to the eye; in the “spots missing” case there were 10 randomly chosen spots missing from across the pupil.

For the well-centered “normal” configuration, the ranges of each iterative algorithm were greatly improved over the conventional algorithm and were generally comparable to each other. As predicted, the B-spline algorithm showed the greatest range. Up to 16.6 D of defocus and 32.9 D of astigmatism could be accurately measured with this algorithm and our chosen lenslet array. Compared to the B-spline algorithm, the spiral algorithm showed identical range for defocus and coma, but 15 to 17% less range for astigmatism and spherical aberration. The Zernike algorithm showed equivalent range for coma and spherical aberration, but 16 to 18% less range for defocus and astigmatism, again compared to the B-spline algorithm.

The robustness of each algorithm to common sources of error is also explored in Fig. 2. We can see that none of the algorithms suffered a significant decrease in range when random spots were missing in the pupil; despite their serial nature, there were no downstream errors for poorly sensed spots. This is by virtue of the large number of spots considered when making each prediction for the B-spline and Zernike algorithms, and by the use of the “dummy” grid for the spiral algorithm. In fact, the B-spline and spiral algorithms were also highly robust in the presence of the other error types considered; range decrease was no more than 6 to 10% across all of the error conditions for these two algorithms. The Zernike algorithm was not robust against decentrations in the estimated position of either the pupil or the input beam, making it advisable when using this method to use a calibrated scale when decentering the input beam,¹⁰ and possibly a separate pupil monitor if maximal range is desired.

It is of interest to know whether the range limits shown here could be improved on with some other as yet undiscovered sorting algorithm, or whether the obtained limits were a result of physical spot overlap (i.e., limited by image processing difficulties rather than the sorting algorithm). To test this, we repeated the range measurements by using the predicted spot locations from geometric optics rather than relying on the MATLAB centroiding functions. The range for the B-spline algorithm increased by 22 to 42%, depending on the aberration type, confirming that for the isolated aberrations considered here this algorithm surpassed the practical limit set by the overlapping of spots. We can therefore confidently state that any future algorithm would not provide a greater range. For reference, the Zernike algorithm showed minimal improvement in range under the same conditions, and the spiral algorithm showed an improvement only for defocus. This means that as far as a pure sorting task is concerned the B-spline algorithm is superior to the other two algorithms, but that this advantage is diminished under practical conditions due to spot overlap.

It is also important to consider the range results in the context of the maximum levels of each type of aberration that the algorithm must be able to measure. For ocular aberrometry it is desirable to be able to measure as many eyes in the population as possible, both on- and off-axis. For astigmatism, the ranges given here are equivalent to ∼30 D for the improved algorithms (the conversion factor between Zernike coefficient and ophthalmic prescription is ∼0.77 for defocus and ∼1.09 for astigmatism at a pupil size of 6.0 mm). This is far in excess of any case of astigmatism ever reported to our knowledge, even in eyes with severe corneal distortion¹⁷ or far away from the optical axis.¹⁸ The defocus results for the improved algorithms correspond to ∼13.3 to 16.6 D of refractive error. Even allowing for a ∼2-D variation over the central 40 deg of the visual field,¹⁸ this is still sufficient to account for >99.9% of the normal human eye population. Since we have already shown that the improved algorithms operate at or close to the hard limit set by the physical overlap of spots, a larger range would require a lenslet array with a shorter focal length (or some other hardware alteration as discussed in Sec. 1). The range for defocus for each algorithm can be adjusted as a linear function of the focal length of the lenslet array used. For example, if a lenslet array with half the focal length were used instead, a roughly equivalent spot pattern would result only when twice as much defocus is present. This is because the wavefront slope, and therefore the spot displacement, for pure defocus varies linearly in x and y across the pupil.

Spherical aberration is very rarely greater than ∼0.4 μm over a 6.0-mm pupil in the normal population,¹⁹ and changes little off-axis,¹⁸ although it may reach ∼1.0 μm in eyes suffering from corneal pathology.¹⁷ This is well within range of the 2.4 to 2.9 μm shown for the improved algorithms here. Coma is very rarely greater than ∼0.5 μm in the normal population¹⁹ and can change by as much as ∼0.6 μm over the central 40 deg of the visual field,¹⁸ which is similarly well within the improved ranges seen here of ∼4.6 to 4.8 μm. However, keratoconic eyes have been reported to feature up to ∼3.6 μm of vertical coma,¹⁷ which may give any of the algorithms considered here difficulty when combined with other aberration types, due to potential spot overlap. If it were desirable to accommodate a keratoconic population of eyes, a lenslet array with a shorter focal length may be desirable.

The preceding results quantify the performance of each algorithm for particular aberration types in isolation. We next used a real adaptive optics system with a model eye to simulate a worse-case scenario—randomized, irregular wavefronts composed of many different types of aberration. The amplitude of the random wavefronts was varied systematically by altering the allowed amplitude of the commands applied to the deformable mirror, with 100 different wavefronts generated for each of the 10 amplitude conditions. Figure 3 shows an example Hartmann-Shack image that was generated using the maximum amplitude, and illustrates the attempts of each algorithm to analyze it. Closed symbols denote spot positions located by the image processing routine. Circles indicate spots that a given algorithm judged were unambiguously assigned to a lenslet. In this particular example, the spiral algorithm did not make any errors, the B-spline made a small number of errors, and the conventional and Zernike algorithms made a large number of errors.

Fig. 3

Example wavefront from the set of 1000 randomly generated wavefronts. This wavefront was generated with the maximum amplitude considered here. Large circles indicate spots that were correctly located according to the internal error checking of each algorithm. We can see that in this specific example, the spiral algorithm was able to correctly sort the spots; the other algorithms show varying amounts of error.

Figure 4 plots the average rate of failure for each algorithm as a function of the maximum allowed actuator amplitude. The failure criterion was met if 5 or more spots (of the ∼280 in the pupil) were missed. The error bars indicate standard error of the mean.

Fig. 4

Average failure rate for the randomly generated wavefronts. Error bars indicate standard error of the mean. We can see that the spiral and B-spline algorithms were most robust to the irregular aberrations that were generated; the Zernike algorithm was significantly less so, and the conventional algorithm was significantly worse again.

The spiral and B-spline algorithms achieved similar results over much of the amplitude range tested; the spiral algorithm was superior at the highest end of the scale. This is probably because the spiral algorithm is far more local in nature and so is somewhat better suited to handle large wavefront deviations that are not well correlated with deviations in other parts of the pupil (the wavefront measured in Fig. 3 is a good example of this). There is a large difference evident between the Zernike algorithm and the spiral and B-spline algorithms. This is likely because the randomly generated mirror shapes would not often have been well represented by the Zernike expansion. Note that this algorithm still performed significantly better than the conventional algorithm, and that the high-amplitude, irregular aberrations produced here are unlikely to be seen in eyes without significant disease affecting the ocular media.

Finally, we report on the processing time for each algorithm. Absolute processing times, of course, depend highly on the computer hardware used. However, even the relative processing times given here should be taken as a rough approximation, since the MATLAB code for each algorithm has not necessarily been optimized for speed. With that in mind, with 10 trials for the well-centered case on an Intel Core i7 2.66 GHz processor with 2 GB RAM and using a 17×17 lenslet array (i.e., detector size of 6.8 mm), we obtained processing times of 2.3, 4.2, 28, and 41 ms for the conventional, spiral, Zernike, and B-spline algorithms, respectively. This excluded components of each algorithm that involved the initialization of variables, which would optimally be performed only once in a session. The spiral algorithm was faster than the other improved algorithms, as predicted due to the maximally simplistic calculations involved in predicting each spot location.

However, the magnitude of temporal variations in ocular aberrations for a given retinal point are generally small and well within the range of the conventional algorithm. A fast frame rate could therefore be achieved by employing a single iteration of any improved algorithm and feeding these results into the conventional algorithm to guide its search. This approach offers a clear benefit to frame rate because of the unique ability of the conventional algorithm to take advantage of parallel processing and so achieve even greater speed than reported here.

5.

Conclusion

The range of aberrometry was similar for the improved algorithms considered here, despite the differences in processing complexity. These algorithms approached the physical limits to the achievable range imposed by the overlap of Hartmann-Shack spots. The results presented here show that each algorithm is more than sufficient to accurately characterize the vast majority of human eye aberration, both on- and off-axis, with our chosen lenslet array parameters. The B-spline algorithm produced a high range and was highly robust to error, but was an order of magnitude slower than simpler algorithms. The Zernike algorithm was somewhat faster, but care must be taken to accurately estimate decentration of the pupil and of the input beam, and performance is noticeably reduced for wavefronts that are not well approximated by the Zernike expansion. The spiral algorithm was highly robust to error and to irregular aberration types. It was also 7 to 10 times faster than the other improved algorithms, affording the potential for real-time analysis. Real-time analysis for ocular aberrometry could alternatively be achieved by employing only a single initial iteration of any improved algorithm to inform subsequent iterations of the conventional algorithm, since the vast majority of human eye aberration is static.