Iteratively Weighted Centroiding for Shack-Hartmann Wave-Front Sensors

: Several techniques have been used with Shack-Hartmann wavefront sensors to determine the local wave-front gradient across each lenslet. In this article we introduce an iterative weighted technique which is specifically targeted for open-loop applications such as aberrometers and metrology. In this article the iterative centroiding technique is compared to existing techniques such as center-of-mass with thresholding, weighted center-of-gravity, matched filter and cross-correlation. Under conditions of low signal-to-noise ratio, the iterative weighted centroiding algorithm is demonstrated to produce a lower variance in the reconstructed phase than existing techniques. The iteratively weighted algorithm was also compared in closed-loop and demonstrated to have the lowest error variance along with the weighted center-of-gravity, however, the iteratively weighted algorithm removes the bulk of the aberration in roughly half the iterations than the weighted center-of-gravity algorithm. This iterative weighted algorithm is also well suited to applications such as guiding on telescopes. useful discussions with Wim de Vries, Scot Olivier, Steven Jones and Lisa Poyneer. This work was performed under the auspices of the U.S. Department of Energy by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48.


Introduction
Shack-Hartmann wave-front sensors are the most prolific wave-fronts sensors currently in use in closed-loop adaptive optics systems. 1 Their popularity is also growing in a number of metrology applications and in the case of aberrometers, used to measure aberrations in the eye, they represent the foundation for the majority of the commercially available instruments. Other metrology applications for which Shack-Hartmann wave-front sensors are also utilized include; optical metrology 2 , diagnosing aberrations present in lasers 3 , wafer inspection and in determining density and magnetic field profiles in plasmas. 4,5 In these applications, the Shack-Hartmann wave-front sensor is not being used in a closed-loop adaptive optics system but rather as a stand alone measurement of the wave-front. The aberrations measured by these systems can be large and noisy. Often times in aberrometers a large number of CCD pixels are used for every subaperture in the lenslet array, which can lead to extremely noisy measurements due to read noise, dark current and scattered light. Each of these noise terms wreaks havoc with conventional centre-of-mass centroiding techniques. In this article current algorithms used to determine the wave-front gradients across the Shack-Hartmann lenslet array are reviewed. Techniques such as noise filtering and, in the case of matched filtering and correlation sensing, utilizing simulated references to improve the performance of these algorithms is examined. An iterative weighted algorithm is then presented which, as is shown in later sections, is particularly well suited for open-loop measurements of aberrated beams. A comparison of the iteratively weighted algorithm with existing algorithms is then presented both in open-loop, as well as in closed-loop. The results are then summarized in the last section of the article.

Existing algorithms
The most basic and widely used technique for determining the gradients from the spot pattern produced by a Shack-Hartmann lenslet array is the centre-of-mass (COM) centroiding algorithm. This technique is expressed mathematically in equations 2 a-c and 3 below with the weighting function, w(x,y), set to unity. The x and y centroids are simply the sum of the intensity pattern multiplied by x and y, respectively, divided by the sum of the intensity pattern. The problems that arise with this technique are due to the read noise, dark current and background light associated with each of the pixels on the wave-front sensor camera. As is readily apparent from Eq. 2 below, the pixels far away from the centre are multiplied by a larger number. As a consequence, the noise contribution far from the centre can quickly corrupt or overwhelm the actual signal when a large number of detector pixels, greater than the minimum four pixels, are used.
One technique that improves upon the centre-of-mass algorithm when noise is present is to utilize thresholding (T-COM) prior to the application of the COM algorithm. 6 This thresholding process identifies the noise level, subtracts the noise level from the subaperture array and sets all of the resulting negative pixels to zero. This algorithm can be expressed by Eq. 2 a-c and Eq. 3 below with the intensity, I(x,y), in these equations replaced by the intensity minus the noise level. The weighting function in these equations is replaced by a binary mask with the location of all positive pixels(I -noise level) set to 1 and all negative pixels set to 0.
Another variation that has recently been proposed is to implement the centroiding algorithm with Gaussian weighting to overcome the limitations due to noise. 7,8 This technique has been called weighted centre-of-gravity (WCOG). In this technique, the authors use a Gaussian weighting function, which is fully characterized by its full-width-at-half-maximum (FWHM), to determine the local gradients in each of the subapertures. This would then correspond to equations 1, 2 a-c and 3 below, however, the Gaussian weighting function of Eq. 1 would not have the ability change its centroid location, x c =0 and y c =0. By choosing the FWHM of the Gaussian weighting function to be less than the size of the subaperture, the noise contribution from those pixels far from the reference spot is reduced significantly. The application of Gaussian weighting was demonstrated by these authors to greatly reduce the effects of noise on determining the local gradients.
A separate technique to determine the wave-front gradients was developed by the solar physics community which came to be known as correlation sensing. 9 The solar physics community utilizes Shack-Hartmann wave-front sensors in closed-loop adaptive optics systems to look at an extended object, the sun. For this application researchers needed a technique that would allow the gradients in each of the subapertures to be determined but they could not utilize the conventional centre-of-mass technique. Their solution was to use crosscorrelation to determine the local gradients in each of the subapertures. Cross-correlation is performed between a given subaperture and all of the remaining subapertures to determine the local displacement of the spots in each of the subapertues, i.e. the displacement of the reference image that provides the maximum correlation. This technique has also been applied to the post-analysis of open-loop intracavity laser and telescope adaptive optics systems. 10 Analogous to this technique is matched filtering, which is very similar to correlation sensing. In this case a convolution in the spatial domain is performed and the peak of the resultant signal is used to determine the gradients in a given subaperture. 11 For both of these techniques the subpixel shift is determined by fitting the correlation peak using parabolic interpolation. 12,13 In this process a three-by-three subarray is extracted from the correlation array with the correlation peak located at the central pixel. The x and y subpixel shifts are given by the x and y derivative(second order) of the correlation peak divided by its x and y Laplacian at the central pixel, respectively. For a symmetric object, these techniques should perform equivalently.
A separate technique that has recently been developed examines the similarity between Shack-Hartmann sensing and carrier frequency interferometry to determine the phase. 14,15,16 In this case the spot pattern formed by the Hartmann sensor is treated like an interference pattern and Fourier transformed. The resulting spots to the right and left of the central "dc" spot have encoded the x gradients and the spots above and below the central "dc" spot contain the y gradient information. By shifting one of the x spots to the centre and filtering out all the information outside of a circle placed about this spot, an inverse Fourier transform returns the wrapped x gradient. Repeating this process with one of the y spots returns the wrapped y gradient. By performing an unwrapping procedure, the gradients at each of the subapertures are determined. In a regime where the aberrations are small, this approach has been further simplified. 17 An iterative windowing technique has also been developed to reduce the noise in situations where a large number of pixels is contained within a given subaperture. 18 This windowing technique estimates the centroid by conventional centre-of-mass techniques and then decreases the window size by one pixel, catered about the estimated centroid location. The authors continue this process until the window size is approximately the theoretical size of the first lobe of the diffraction pattern. This technique is equivalent to Eq. 2 a-c and Eq. 3 below with the weighting function equal to one within the area of the window and zero outside of the window. The window is then iteratively decreased in size and kept centred on the previous centre-of-mass estimate.

Iterative weighted algorithm
The basis of this algorithm finds its roots in the weak lensing astrophysics community. In this case the desire is to determine the ellipticity correlations between galaxies in an effort to map out the dark matter in the universe. One such weak lensing data pipeline uses a Gaussian weighting function which is characterized by its FWHM and by its centroid location, both of which can be iteratively adjusted. 19,20 From this weighted distribution, the moments, up to second order, of the intensities in each of the subapertures can be found. The x and y centroids are simply the sum of the intensity pattern multiplied by the weighting function and x and y, respectively, divided by the sum of the intensity pattern multiplied by the weighting function. These can then be iteratively applied to the centre-of-mass of the original weighting Gaussian function to refine the location of the centroid of the Gaussian and minimize the noise. In addition, using the higher order moments of the intensity pattern provides an estimation of the FWHM of the Gaussian which can also be iterated upon along with the centroid of the Gaussian. This then forms the basis for an iterative scheme which sets the centroid and FWHM of the weighting function in the presence of noise. In the limit that the ability to iteratively change the location of the centroid and its FWHM is removed, this method reduces to the weighted centre-of-gravity (WCOG) algorithm described above. By allowing the iterative process to occur, however, the algorithm is able to more accurately locate the centroids in cases of large aberrations where the Hartmann spots are not close to the reference spots and noise is present. In addition, this technique could be used to more quickly bring the AO system to the null operating point when operating in closed-loop as discussed below.
The weighting function, w(x, y), is expressed by a Gaussian distribution centred at (x c ,y c ) with a full-width-at-half-maximum equal to 2σSQRT(2ln2). The weighting function is given by with σ representing the root-mean-square (RMS) deviation of position from the mean. This definition of the weighting function results in a normalization of the integral over all space of unity. For the simulations, this integral is approximated by a sum of the weighting function over the given subaperture. In this expression x, x c , y, y c and σ represent real numbers. The various position moments, up to second order, are expressed as where S w , S x , S y , S xx and S yy represent the different moments of the lenslet spot pattern, I(x,y), on the CCD camera with the weighting function. The average location of the Gaussian, or its centroid, is simply where x c and y c represent the x and y centroid locations. For a Gaussian spot, the RMS deviation of position from the mean corresponds to the weighted (x-x c ) 2 + (y-y c ) 2 This set of equations form the basis for the iteratively weighted centroiding algorithm in which the centroid location and the FWHM of the weighting function can be adaptively changed.
As pointed out by the referee, this technique could also be utilized using a recentred average of the Hartmann spots as the weighting function rather than a Gaussian or Sinc 2 function. This technique is also applicable to extended scenes as in the case of correlation sensing. For this choice of the weighting function the FWHM would not be changed such that Eqs. 2 a-c and 3 would be iterated upon with the weighting function consisting of a normalized average of the recentred Hartmann images/spots. For the simulations presented in this article, however, a Gaussian weighting function was used and Eqs. 1-4 were utilized for the iterative weighted algorithm. It is also the case that σ can be fixed and not iterated upon with a Gaussian or Sinc 2 weighting function as well.
To investigate the convergence of the iteratively weighted algorithm, numerical wave optics simulations were performed. The simulations consisted of a lenslet array with 16 lenslets across the diameter of a circular pupil. The phase aberration used to test the iteratively weighted algorithm was based on a scaled version of Zernike coefficients measured in a human eye. 21 The Hartmann spots were then formed by Fourier transforming the electric field at each of the lenslets in the array and multiplying by the complex conjugate to obtain the intensities. Both photon noise and read noise were then introduced to the focal plane containing the Hartmann spots. The iteratively weighted algorithm was used to estimate the slopes which were then reconstructed. The variance between the original and reconstructed phases was then calculated. The variance in the reconstructed wavefront, as a function of the number of iterations was then used to test the convergence of the iteratively weighted algorithm. Figure 1 shows the residual variance as a function of iterations through Eqs. 1-4. This figure shows that the algorithm converges within about four steps for this initial aberration(peak-to-valley ~ 15 radians). Based on this result, all of the simulations in this article which use the iteratively weighted algorithm were performed with four iterations. The numerical simulations can also be utilized to test the width of the Gaussian weighting function after the four iterations as a function of the initial guess. Figure 2 shows the width of the Gaussian weighting function, σ, after four iterations, solid red line, as a function of the initial starting value, σ init , for the Gaussian width. The spread in starting values for σ ranged over a factor of close to 5, from 1.26 to 6 pixels. The resultant values for the weighting functions, σ, after four iterations varied by less than seven percent. The simulated Hartmann spots had a full-width-at-half-maximum of 3.8 pixels which corresponds to a Gaussian σ s of 1.6 pixels. Figure 2 also shows the reconstructed wavefront variance, solid black line, using four iterations of the iterative weighted algorithm as a function of the initial value for the weighting function, σ init . The weighting function σ that resulted in the lowest reconstructed wavefront variance was σ = 1.62 pixels. This illustrates for this particular phase profile that over a broad range of initial starting values of the weighting function, σ init , both larger and smaller, that after four iterations the algorithm converges to a σ value that is within several percent of the Hartmann spot full-width-at-half-maximum. With this particular phase profile, an initial starting value, σ init = 2.25 pixels, that was slightly larger than the spot size, σ s = 1.6 pixels, resulted in the lowest variance.
where N ph is the number of photons, N W is the FWHM of the Gaussian weighting function, N T is the Full Width at Half Maximum (FWHM) of the seeing-limited Point Spread Function (PSF) in pixels, N D is the FWHM of the diffraction limited point spread function and σ detector is the standard deviation of the detector noise. The iteratively weighted algorithm could also be applied to the functions of field acquisition and guiding on astronomical telescopes. The weighting function in Eq. 1 could also be different in each of the subapertures to include effects such as laser guide star parallax for the next generation of large diameter astronomical telescopes. A more general expression for the Gaussian weighting function would then be expressed as w(x,y)=(2πσ x σ y ) - 1

EXP{-[(xx c )cosθ-(y-y c )sinθ] 2 /(2σ x 2 )-[(x-x c )sinθ+(y-y c )cosθ] 2 /(2σ y 2 )}
with σ x and σ y representing the root-mean-square (RMS) deviation of position from the mean along the x and y axes, respectively. The angle θ represents the rotation of the Gaussian from the x-axis. For the results presented in this article, a symmetric weighting function was assumed such that θ = 0 degrees and σ x = σ y = σ. In this case, the more general weighting function reduces to the form given in Eq. 1 above.

Noise filtering and reference generation
Nearly all of the centroiding techniques, except for the conventional centre-of-mass algorithm, perform better after filtering when the spots formed by the lenslet array are larger than a single pixel. In this case, the noise, which varies at the spatial scale of an individual pixel, can be reduced by filtering the signal prior to analysis by the above mentioned gradient determining algorithms without affecting the signal itself. This is not the case when the spot size from the lenslet array is chosen to be a single pixel in width such that its spatial variation is comparable with the background noise. In this case filtering will affect both the noise and the signal you are trying to detect. For this article, where the spot size was larger than a single pixel, several filtering techniques were investigated for their utility in reducing the noise and improving the reconstructed wave-front variance. The techniques examined included filtering with wavelet transforms, Fourier transforms and spatial convolution with Gaussians. In general the wavelet transform resulted in the lowest residual variance followed by Fourier transforms and then filtering using a spatial convolution with a Gaussian. Other authors have used filters previously including triangular and square convolutions in real space 22 and box car averaging in Fourier space, 14,22 which was required for this particular algorithm.
The wavelet transform that was evaluated for noise filtering was specifically the D trous algorithm, which is for use with discreet data. 23 The image to be filtered in this case is the intensity spot pattern recorded on the wave-front sensor due to the lenslet array. The different scales of the image are formed by taking convolutions of the image with a scaling function that maintains its discreet shape but doubles in spatial scale between successive wavelet scales or convolutions. A B 3 -spine was chosen as the scaling function such that the coefficients of the convolution mask are given by the one-dimensional array, (1/16, 1/4, 3/8, 1/4, 1/16), for a vector and the two-dimensional array, for images, where ⊗ represents the Kronecker product. 23 From the above description, it is readily apparent that the method of removing noise using spatial convolution with a Gaussian is quite similar to the implementation of the D trous wavelet transform. Fourier transforms were also investigated as a means of removing the high spatial frequencies associated with the detector noise. In this case the image was Fourier transformed and a circular binary mask applied to block out the high spatial frequencies before inverse Fourier transforming back. This is similar to the filtering method used in carrier frequency interferometry. 16 The spatial filtering in this article was applied to the focal plane containing all of the Hartmann spots, as opposed to each of the subapertures separately. This focal plane was larger than the pupil and so no zero padding was utilized as the distance between the outer Hartmann spots and the edge of the focal plane was greater than the extent of the filters.
The three filtering schemes were tested on the reconstruction of a phase profile with a circular aperture. The read noise on the detector was set at 20 e-root-mean-square and the number of photons per subaperture was varied between 50 and 1000. The results, shown in Fig. 1 below, demonstrate the utility of filtering the signal before determining the gradients when the signal-to-noise-ratio decreases as the unfiltered signal, the solid green line, had the highest residual variance. The D trous wavelet transform, dashed black line, had slightly better performance, particularly at the lowest SNR, than the Fourier transform filter, solid blue line, which performed slightly better than the method utilizing convolution with Gaussians, dashed red line. In addition to filtering the noise, centroid locating techniques such as correlation sensing and matched-filter also need a reference image. This reference can be taken from one of the subapertures, filtered or unfiltered, or it can be generated analytically without noise. This can lead to better performance for high noise scenarios and also enables the techniques to accurately reconstruct the tip/tilt component of the phase aberration, a feature which is compromised when comparing to one of the subapertures. If one of the Hartmann spots from the lenslet array is used as the reference, then that spot can still be compared to a simulated spot to recentre the reference before that spot is compared with all of the remaining spots in the lenslet array to allow determination of the tip/tilt information. Figure 4 below shows the residual variance as a function of photons/subap for the cases where the reference is simulated(dashed black line), taken from the noise filtered centre subaperture(solid blue line) and taken from the unfiltered centre subaperture(dashed red line). The intensity profile on the CCD camera was filtered with the Fourier transform filter and the correlation sensing method was used to determine the centroids. The tip/tilt was also removed from the gradients before and after reconstruction to place all three reference generation methods on a comparable footing. As Fig. 4 shows, the simulated reference performed better than the technique of using one of the subapertures at lower SNRs. The overall variance is lower in Fig. 4, relative to Fig.  3, due to the subtraction of the tip/tilt aberration. Fig. 4. Residual variance as a function of noise for the three different reference generation methods.

Comparison of algorithms
A comparison of the algorithms discussed above was performed. This comparison involved starting from an initial phase aberration, expressed as Zernike coefficients, at a circular aperture. At this plane a lenslet array was simulated such that there were a total of sixteen lenslets across the diameter of the circular aperture and the electric field across each of the lenslets was parsed into 32 by 32 phase elements. In addition, each lenslet represented 16 by 16 pixels on the CCD wave-front camera. The Hartmann spots were then formed by Fourier transforming the electric field at each of the lenslets in the array and multiplying by the complex conjugate to obtain the intensities. Both photon noise and read noise were then introduced to the focal plane containing the Hartmann spots. In all of the simulations the read noise was fixed at 20 e-rms and the number of photons-per-pixel was varied. The residual variance between the original phase and the reconstructed phase was then determined as a function of the number of photons-per-pixel. The different algorithms were then used to identify the displacement of the Hartmann spots caused by the phase aberration and a Multigrid reconstructor 24 was used to reconstruct the phase based on the local gradients determined by the different algorithms. The initial phase aberration used to test the various centroiding algorithms was based on a scaled version of Zernike coefficients measured in a human eye. 21 Through filtering the noise and utilizing a simulated reference spot, the cross-correlation, matched filtering, centre-of-mass, windowed and centre-of-mass threshold algorithms detailed above were optimized before comparing them to the weighted centre-ofgravity and the weighted iterative algorithm expressed in Eqs. 1-4. The results from the various algorithms are shown below in Fig. 5. At the higher photonsper-subaperture end of the graph, all of the algorithms except the weighted centre-of-gravity (WCOG) algorithm, solid black line, converge as expected. This algorithm is higher than the others simply because the aberration is large such that the spots are displaced from the nominal centre of the CCD pixels representing the lenlset array. To encompass the displaced spot, the WCOG algorithm must have a large FWHM Gaussian weighting function which makes it much more like the conventional centre-of-mass algorithm and thus very susceptible to noise. As shown below, when the WCOG algorithm is operated in closed-loop the spot displacement is brought to zero allowing a much smaller FWHM to be used. This enables the WCOG algorithm to achieve the best performance along with the iteratively weighted algorithm. As the number of photons-per-subaperture is decreased, the measurements become noisier and the remaining algorithms begin to diverge in their performance. After the WCOG algorithm, the threshold centre-of-mass (T-COG) algorithm, dashed blue line, and the windowed algorithm, solid red line, have the next largest variance. This is followed by the matched filter, solid green line, and cross-correlation algorithms, solid grey line. These two algorithms have virtually identical performance except at the noisiest point in the graph where the matched filter algorithm has a slightly higher variance than the cross-correlation algorithm. This is expected as the Hartmann spots are symmetric and both techniques use parabolic interpolation to achieve subpixel accuracy. The best performance is achieved by the iteratively weighted algorithm, dashed black line, presented in this article. The performance improvement over the existing algorithms is particularly high under conditions of low signalto-noise ratio as can be seen in Fig. 5. As the number of photons decreases, the slope of the variance vs. photon number, N ph , for the iterative weighted algorithm closely matches the (1/N ph ) 2 , dashed grey line, behaviour expected form Eq. 5b above. It is also of interest to determine how well these algorithms operate in closed-loop since nearly all adaptive optics(AO) systems operate in this manner and AO systems represent a large fraction of the Shack-Hartmann wave-front sensors in operation. The same phase aberration discussed above was used to test the closed-loop operation of the algorithms. The phase was reconstructed as detailed above, but a gain multiplier of 0.8 was applied to it before being subtracted from the original phase. A comparison of the iteratively weighted, WCOG, cross-correlation and centre-of-mass algorithms was undertaken and the results are shown below in Fig. 6. Figure 6a represents the lower signal-to-noise-ratio of the two figures at 105 photons-per-subaperture and 20 e-rms noise. The centre-of-mass, dashed grey line, exhibits the worst performance followed by an approximately equivalent performance by the crosscorrelation algorithm, solid grey line, the threshold centre-of-mass algorithm, solid blue line, and the windowing algorithm, solid red line. The iteratively weighted, dashed black line, and WCOG algorithms, solid black line, converge to the same nominal level. The iterative functionality of the iterative weighted algorithm was turned off after the first three iterations to reduce the computational load. By incorporating the iterative capability in the first few times through the loop, however, the rate at which the aberrations are removed is accelerated considerably over the WCOG algorithm as seen in Fig. 6a and 6b. Figure 6b represents a higher signal-to-noise ratio such that there is less of a disparity between the cross-correlation algorithm and the iteratively weighted and WCOG algorithms. Again the iteratively weighted algorithm removes the bulk of the phase aberrations in roughly half the iterations through the loop, at gain=0.8, than the WCOG algorithm.  Figure 4a illustrates the results with 105 photons/sub-aperture and Fig. 4b represents the case of 500 photons/sub-aperture.

Summary
In this article an iteratively weighted centroiding algorithm was developed for wave-front sensors that measure the gradient of the wave-front, such as Shack-Hartmann and shearing interferometers. This algorithm was specifically designed to mitigate the deleterious effects of noise and is ideally suited for use in open-loop applications where the Hartmann spots are far from their unaberrated positions. A comparison between the iteratively weighted centroiding algorithm was performed with other commonly used centroiding algorithms including; centerof-mass, center-of-mass threshold, cross-correlation, matched filter, windowing and weighted center-of-gravity. In order to provide the best comparison between the algorithms, methods such as filtering the noise on the simulated CCD wave-front camera and using simulated reference spots were investigated to improve the performance of algorithms such as center-ofmass threshold, cross-correlation and matched-filter. In the case of filtering the noise on the simulated CCD wave-front camera, Fourier filtering, filtering with D trous wavelet transform and filtering by Gaussian convolution were investigated. The best performance was achieved with the D trous wavelet transform, however, all three methods of filtering out-performed the unfiltered case. The filtering was performed for the case where the full-width-at-halfmaximum of the Hartmann spots on the CCD wave-front camera was several pixels, 3.8, in width. The matched-filter and cross-correlation algorithms were also compared using simulated reference spots and using a spot from one of the lenslets as a reference, both noise filtered and unfiltered. Because the method of using one of the spots from a lenslet as a reference is more susceptible to tilt errors, the tilt was removed for these comparisons. The results indicated that using a simulated spot as the reference resulted in a lower error variance in the reconstructed wave-front than using a spot from one of the lenslets as a reference.
The performance of all of the algorithms compared was significantly better than the traditional center-of-mass algorithm in both open and closed-loop as expected. The open-loop comparison showed that the iteratively weighted algorithm achieved the lowest residual variance of all of the algorithms tested, as shown in Figs.5, under conditions of low signal-tonoise ratio. This was the case even when the optimizing techniques discussed above were used on the other algorithms. The iteratively weighted centroiding algorithm was also investigated in closed-loop operation, the most common mode of operation for adaptive optics systems. In this case the iteratively weighted and weighted center-of-gravity centroiding algorithms gave comparable levels of performance and resulted in a lower variance than the other algorithms. The iteratively weighted centroiding algorithm did, however, remove the bulk of the phase aberrations in roughly half the iterations than the WCOG algorithm.