1.

Introduction

The Shack–Hartmann wavefront sensor (SHWFS) continues to be the most widely used wavefront sensor in astronomy with the most mature technology. It is essentially composed of an array of microlenses and a detector, typically a charge coupled device (CCD) array, situated at the microlenses focal plane. The system’s pupil is reimaged onto the array of microlenses, thereby being split into subpupils and creating an array of images of the observed object at the detector. The relative displacement of these images between the aberrated and plane incident wavefront cases is proportional to the wavefront tilt over each microlens. The estimation of such displacements allows to retrieve the aberrated incident wavefront profile.¹ This sensor’s parameters, such as the pupil’s sampling factor and subpupil’s image field of view (FoV) and resolution, can be selected according to the application requirements of wavefront estimation precision, sensitivity, and dynamic range. This flexibility is in contrast to other wavefront sensing techniques, for example, the pyramid wavefront sensor (P-WFS).² This other technique achieves greater sensitivity than the SHWFS but relies for its operation on the at least partial correction of the wavefront aberration to reduce its dynamic range.³

Multiobject adaptive optics (MOAO) systems span a wide sensed FoV, of the order of arcminutes, and correct only those reduced portions of the FoV (of the order of arcseconds), where the scientific objects of interest are situated.^4–6 Thus, they operate in an open-loop correction configuration, and a requirement of their wavefront sensors is to maintain their sensitivity in the low light level conditions when sensing a larger FoV than in traditional closed-loop systems, to cope with the uncorrected atmospheric turbulence dynamic range. For SHWFSs, this requirement applies directly to the centroiding method employed to estimate the tilt at the subpupil level from the spot’s position in the focal plane image, and therefore justifies the revision of traditional centroiding methods and the proposal of new ones (as suggested in the future work section in Ref. 6).

Hardware speed and processing capabilities continue to improve. The hundreds, even thousands, of GFlops of processing power of modern graphical processing units (GPUs), for example, encourage us to move from the simple centre of gravity (CoG)-based traditional centroiding methods, which work in the subpupil image domain, toward more robust methods that work in transformed domains where optimum estimators could be more easily obtained.

In this paper, we present the characterization and comparative performance results of a centroiding method formulated as an optimized Bayesian estimator in the Fourier domain, for Shack–Hartmann wavefront sensors and point-like guiding sources, which we have called the weighted Fourier phase slope (WFPS) algorithm. Its performance is compared with the thresholded centre of gravity (TCoG) and cross-correlation (CC) centroiding methods through simulations first at the subpupil level, then at the pupil level, and finally with real SHWFS images obtained with a laboratory setup.

2.

Description of the Weighted Fourier Phase Slope Algorithm

The formulation of the WFPS has already been presented and validated as a tilt estimation method at the subpupil level of SHWFSs.^7,8 This algorithm first calculates the horizontal and vertical displacements $C_{k,l}^x$ and $C_{k,l}^y$ of the subpupil image $I_{xy}$ evaluated at the discretized spatial frequencies determined by the integer indices $k$ and $l$ based on the slopes, or derivatives with respect to both spatial frequency axes, of the phase of the Fourier transform of $I_{xy}$, without the need to compute explicitly the phase or its unwrapping^9,10

$V_{E^x}$ and $V_{E^y}$ are the covariance matrices of errors $E^x$ and $E^y$. Their role is to give more weight to the cleanest or least noisy displacement measurements in $C_v^x$ and $C_v^y$. Moreover, they take advantage of the strong correlation of the errors at neighboring spatial frequencies to cancel out noise in an optimized manner. This is equivalent to smoothing the rapid variations in the Fourier phase shape as a function of spatial frequency and has a symmetrizing effect in the image domain. Thus, it can also be viewed as a low pass filter in the cepstrum domain of the $I_{xy}$ subpupil image.

2.1.

Relationship of the Weighted Fourier Phase Slope Algorithm with Cepstrum Filtering Techniques

The complex cepstrum of image $I_{xy}$ is defined as the image ${\widehat{I}}_{xy}$ with the Fourier transform

Eq. (4)

$${\mathcal{F}\{{\widehat{I}}_{xy}\}|}_{{\omega }_x,{\omega }_y}=\log [{\mathcal{F}\{I_{xy}\}|}_{{\omega }_x,{\omega }_y}]=\log |{\mathcal{F}\{I_{xy}\}|}_{{\omega }_x,{\omega }_y}|+j\arg ({\mathcal{F}\{I_{xy}\}|}_{{\omega }_x,{\omega }_y}),$$

where log is the natural logarithm, $|·|$ is the magnitude, and $\arg (·)$ is the unwrapped phase.¹⁰ Low-pass cepstrum filtering consists in convolving the log of the Fourier domain with a broad shape, thereby smoothing the spectral shape. This has a denoising, dereverberating, and symmetrizing effect on the original image domain. On the other hand, taking the phase slope at the origin of spatial frequencies is equivalent to calculating the pure CoG in the image domain [this can be derived from Eq. (1)]. These facts lead us to the relationship between the proposed algorithm and filtering in the cepstrum domain: averaging the phase slope around the origin of spatial frequencies in the Fourier domain implies smoothing the Fourier’s phase or, equivalently, low-pass cepstrum filtering of the image, and subsequently taking the CoG of the filtered image in the original domain.

Panels (a), (b), and (c) in Fig. 1 show simulated focal plane subpupil images obtained with the simulator described in Sec. 4. Measurement noise is dominated by asymmetric optical aberration in panel (a), by the detector noise in panel (b), and by spot’s truncation due to a limited FoV in panel (c). Panels (d), (e), and (f) are the same images after low-pass cepstrum filtering. In the case of panel (f), the original image has been zero padded to a bigger FoV before the cepstrum transform. In all cases, the low-pass filter has been calculated from the WFPS Bayesian weights. The green crosses in Fig. 1 represent perfect $Z$-tilt estimation. The red and blue dots are the results of pure CoG application over the original and filtered images, respectively. A considerable improvement in the $Z$-tilt estimation is achieved by low-pass cepstrum filtering, owing to its capacity to attenuate distortion and noise, and symmetrize the spot. Relating the WFPS algorithm to the cepstrum deconvolution helps us become aware of its effect in the image domain; the proposed algorithm could be viewed as the subset of operations from the cepstrum filtering technique strictly required to estimate the centroid, to minimize its computational cost.

Fig. 1

Simulated focal plane subpupil images with (a) asymmetric aberration, (b) detector noise, and (c) limited FoV; and (d)–(f) their cepstrum low-pass filtered counterparts. Green crosses correspond to perfect $Z$-tilt estimation; red and blue dots are the results of pure CoG application.

3.

Comparison of Weighted Fourier Phase Slope Computational Cost Against Traditional Centroiding Methods

The computational cost of the WFPS algorithm is dominated by the two-dimensional (2-D) DFTs in Eq. (1). These DFTs need not be fully calculated but only at those $N_f$ spatial frequencies considered in Eq. (2). As will be seen in Sec. 4, six spatial frequencies give a good tilt estimation performance; including more spatial frequencies increases the computational cost with no appreciable tilt estimation improvement. Hence, $N_f=6$ will be considered in the present section.

Specifically, Fourier phase slopes will be assumed to be obtained at frequencies $f_{00}$, $f_{01}$, $f_{02}$, $f_{10}$, $f_{11}$, and $f_{20}$, with $f_{kl}$ being the 2-D spatial frequency composed of the horizontal $f_k$ and vertical $f_l$ frequencies. So, if we compute the $N$ by $N$ 2D-DFTs in Eq. (1) by first obtaining the one-dimensional $N$ point DFTs of the rows, only the first three output values will be necessary. When afterward applying $N$ point DFTs to the resulting three columns, only three, two, and one output values are required, summing up the six output values that we are looking for.

When computing these partial 2D-DFTs for a typical subpupil image size between $12\times 12$ and $16\times 16\text{pixels}$ with real data, the highly parallelizable direct method, which computes every output of the 1D-DFT independently, is advantageous from a computational latency point of view over the Cooley–Tukey, prime factors or FFT pruning schemes, which are recursive or iterative in nature.¹⁰ Application of the direct method to the 1D-DFTs means $N$ complex multiplications and $N-1$ complex sums for every computed DFT output value. Provided the DFTs are applied to rows first and then to columns in the 2D-DFT, the multiplications are real by complex for the row DFTs and complex by complex for the column DFTs. Three such 2D-DFTs are computed in Eq. (1). The product of coordinates $x$ and $y$ by phasors in the same equation can be done offline. Operations in Eq. (2) can be disregarded for the present estimation.

With the purpose of assessing the WFPS algorithm’s computational cost against other control algorithms, the weighted centre of gravity (WCoG), TCoG, and CC algorithms have been selected. As for the CC algorithm, we are assuming the conventional method of implementing the correlation operation in the Fourier domain, with an $N_{\mathrm{fft}}$ by $N_{\mathrm{fft}}$ size of the 2D-FFT large enough to avoid spatial aliasing, and that a TCoG is applied to the correlation figure to obtain the sought-after centroid. Thus, we guarantee the robustness and linearity of the CC method.¹¹ The “phase correlation” technique estimates the displacement of the correlation figure in the Fourier domain, thereby considerably reducing the computational cost of the correlation;^12–14 however, it has been reported to be very susceptible to noise¹² and is only adequate for small displacements between live and reference images due to the Fourier phase wrapping around the $-\pi$ to $\pi$ range, which leaves no reliable phase values when the displacement is of the order of half the FoV, the 2D-DFT is the same size as the FoV and in the presence of noise.¹³ In contrast, the “periodic correlation” technique, consisting of obtaining the correlation figure in the Fourier domain with 2D-DFTs of the same size as the live image, has reported to be adequate only when the image quality is good, due to spatial aliasing.¹²

Table 1 shows the results of computational requirement estimations for an example of SHWFS with $80\times 80$ subapertures, $16\times 16$ detector pixels per subaperture, the imposed condition of finishing the complete centroids computation in $100\mu \mathrm{s}$, and for the four centroiding methods. Regarding the computational cost of the CC method, the 2D-FFTs have to be computed in their entirety. $N_{\mathrm{fft}}$ was varied allowing for different $N_{\mathrm{ref}}$ by $N_{\text{ref}}$ sizes of the reference image. The Cooley–Tukey scheme has been assumed for $N_{\mathrm{fft}}=32$ in the 1D-FFTs, and the Prime Factors scheme for $N_{\mathrm{fft}}$ equal to 20 and 24. The resulting figures in Gflops include the total number of real sums and multiplications per second required. It can be seen that the computational cost of the WFPS algorithm is an order of magnitude higher than that of centre of gravity-based algorithms but still within the Tflop limit. The CC algorithm requires yet a further order of magnitude more computational capacity, reaching a value close to 10 Tflops.

Table 1

Computational cost comparison of WCoG, TCoG, WFPS, and CC algorithms, with an example.

4.

Performance Simulations at the Subpupil Level

4.1.

Description of the Simulator

We have first evaluated and understood the behavior of the WFPS algorithm and compared it with other known centroiding methods through numerical simulations at the subpupil level using the MATLAB tool from Mathworks, Inc. The simulation workflow follows the block diagram shown in Fig. 2. A single Shack–Hartmann subaperture is illuminated with a wavefront fulfilling Kolmogorov phase statistics, simulated as in Ref. 15. A point-like source, uniform light intensity across the subaperture, and statistical independence between frames are assumed. A Fraunhofer integral is then applied to each complex light field frame to obtain the corresponding image at the microlens focal plane, and finally the electron multiplying charge coupled device (EMCCD) detector gain and noise model from Ref. 16 is applied to obtain a close-to-real simulated subpupil image sequence, except that it does not consider temporal turbulence dynamics. The incident phase tilt is then estimated by applying different centroiding methods to this image sequence, and the result is compared with the known applied true tilt, thus obtaining a tilt estimation error for each centroiding method as its figure of merit. Two possible “true” tilts have been considered: the Zernike tilt obtained at the Kolmogorov phase simulation, or $Z$-tilt, and the mean of the phase gradient, or $G$-tilt, which is obtained by applying a true CoG to a version of the image sequence with maximum FoV and exempt from detector noise.

Fig. 2

Block diagram of the simulation workflow at the subpupil level in the MATLAB tool.

A parallel training sequence, with 50,000 frames and with exactly the same system parameters as the sequence under study, is simulated to obtain the necessary weights for the WFPS algorithm. As described in Sec. 2, the Fourier phase slopes of the sequence at a set of spatial frequencies are obtained by applying Eq. (1) to each image. The resulting spot displacement measurements are compared with the true applied tilt to obtain the measurement errors, and MAP weights are a function of the measurement error covariance matrices. The number $N_f$ and order of the spatial frequencies involved in the WFPS algorithm, and consequently the number of weights, are selectable and determined in the simulation as a compromise between computation latency and tilt estimation precision.

Simulations in this section assume an adaptive optics system with an open-loop correction scheme, so the TCoG and CC algorithms have been considered as comparison methods for their good linearity response even with wide FoVs.¹¹ WCoG is in principle adequate for closed-loop controlled systems.¹⁷ TCoG’s threshold is a unique value for the whole sequence under study, optimized by minimizing the tilt estimation error. As for the CC algorithm, the reference image is a 2-D Gaussian shape with a two-pixel full width at half maximum (FWHM) and an $8\times 8\text{pixel}$ FoV. The correlation figure is obtained with no interpolation, i.e., with the same resolution as the live and reference images. The centroid is calculated over the correlation figure with a TCoG algorithm for which the threshold is again a unique optimized value for the whole sequence.

The centroiding algorithms give the result of the spot’s displacement in pixel units. However, r.m.s. wave radians are required for the tilt estimation error results. Equation (5) relates displacement in pixels $d$ with the peak-to-valley phase tilt $\mathrm{\Delta }\varphi$ across the subaperture extent at the telescope pupil plane in wave radians¹¹

Eq. (5)

$$\mathrm{\Delta }\varphi =\left[\frac{2\pi }{\mathrm{FWHM}}\right]d.$$

FWHM being the full width at half maximum of the simulated spot at the diffraction limit in pixels. For circular subapertures, peak-to-valley tilt is divided by a factor of four to obtain r.m.s. radians, whereas for a square subaperture the factor is $2\sqrt{3}$ (as derived by the authors).

All the parameters defining the system at subpupil level can be easily changed in the simulation; hence, the sensitivity of the different algorithms as a function of the light level flux and detector noise, their linearity as the FoV of the subpupil image changes, and their robustness against such effects as higher orders of turbulence phase when the Fried parameter decreases, can be easily evaluated. Moreover, optimization of every centroiding method’s parameters, such as MAP weights for the WFPS, is achieved through this type of simulation.

4.2.

Optimum Field of View

The optimum working FoV for a particular sensing condition is a means of assessing the capacity of each algorithm to maintain its sensitivity when an increase in FoV is required to gain in dynamic range. Extended FoVs imply more sensing pixels, more influence of the detector noise, and less sensitivity. Therefore, those algorithms that are more robust against detector noise will be optimum at larger FoVs, thus allowing for larger dynamic ranges of tilt sensing, as required with strong turbulence and in an open-loop wavefront sensing.

Figure 3 shows the dependence of the tilt estimation error in r.m.s. radians at the sensing wavelength with the FoV in pixels by pixels and light flux level in photons per subpupil, for the TCoG algorithm (dashed red lines), the CC algorithm (dash-dotted green lines), and the WFPS algorithm with $4\times 4$ spatial frequencies (solid blue lines), when estimating $Z$-tilt (a) and estimating $G$-tilt (b). $D_{\text{sub}}/r_0$ is 2.5. The spot’s FWHM size is 2 pixels at diffraction limit. Clock induced charge (CIC) is $0.05{\mathrm{e}}^{-}/\mathrm{pix}/\text{frame}$, and readout noise (RON) is $50{\mathrm{e}}^{-}$ rms. Selected EM gain for each light level is high without saturating a 14-bit detector, with a maximum value of 1000. The detector quantum efficiency (QE) is 97% and its sensitivity is $10{\mathrm{e}}^{-}/\mathrm{ADU}$ (electrons per analog to digital unit). These values for the detector parameters are maintained throughout this section. Table 2 shows the optimal FoV values for the different light flux levels and the three algorithms under study. A general conclusion that can be extracted from this simulation is that the WFPS algorithm is the least affected by detector noise among the three evaluated algorithms; it calls for wider FoVs for each light level under study and gives as a result a better tilt estimation at the optimal FoV. The CC algorithm would be in second place in this regard, with a very similar behavior. However, the TCoG is very much affected by detector noise and requires small optimal FoVs. Moreover, the dashed red traces in Fig. 3 tend to increase rapidly with the FoV value for light levels below 100 photons, thereby indicating that the TCoG algorithm is very sensitive to the optimal selection of FoV. WFPS and CC are much more tolerant in this sense in that they allow for a broader range of FoV selection with little or no penalization in tilt estimation error for light levels above 30 photons. For low light levels, such as 30 to 100 photons, the TCoG must work with FoVs in the order of $10\times 10$ to $12\times 12\text{pixels}$, whereas CC and WFPS can work with FoVs of $12\times 12$ to even $16\times 16\text{pixels}$. Similar conclusions are obtained for the $D_{\text{sub}}/r_0=1\mathrm{case}$.⁹

Fig. 3

Tilt estimation error in r.m.s. radians as a function of FoV and light flux level for the TCoG algorithm (dashed red), the CC algorithm (dash-dotted green) and the WFPS algorithm with $4\times 4$ selected spatial frequencies (solid blue), (a) for $Z$-tilt estimation and (b) $G$-tilt estimation. $D/r_0$ is 2.5. CIC is $0.05{\mathrm{e}}^{-}/\mathrm{pix}/\text{frame}$ and RON is $50{\mathrm{e}}^{-}\mathrm{rms}$.

Table 2

Optimal FoV in pixels × pixels as a function of light flux level for the TCoG, CC, and WFPS (with 4×4 frequencies) algorithms, obtained from Fig. 3. D/r0=2.5. CIC=0.05  e−/pix/frame. RON=50  e− rms.

Optimal FoV for Z-tilt/G-tilt estimation (pixels×pixels). D/r0=2.5
#photons (per subpupil)	TCoG	CC	WFPS (4×4 freqs)
10	$10\times 10/10\times 10$	$10\times 10/10\times 10$	$10\times 10/10\times 10$
20	$10\times 10/10\times 10$	$10\times 10/12\times 12$	$12\times 12/12\times 12$
30	$10\times 10/10\times 10$	$12\times 12/12\times 12$	$12\times 12/12\times 12$
50	$12\times 12/12\times 12$	$12\times 12/12\times 12$	$14\times 14/16\times 16$
100	$12\times 12/12\times 12$	$16\times 16/12\times 12$	$16\times 16/16\times 16$
300	$14\times 14/16\times 16$	$16\times 16/16\times 16$	$16\times 16/16\times 16$
10,000	$16\times 16/16\times 16$	$16\times 16/16\times 16$	$16\times 16/16\times 16$

4.3.

Maximum-a-Posteriori Weights

A key aspect of the WFPS algorithm is obtaining the value of the MAP weights at the spatial frequencies considered for the centroid computation. These are determined at the simulations at the subpupil level, in which the true applied tilts for a large sequence of frames are known, so that the error vectors $E^x$ and $E^y$ in Eq. (3) and their covariance matrices $V_{E^x}$ and $V_{E^y}$ in Eq. (2) can be determined. This is the reason of the 50,000 frame training sequence in Fig. 2.

Table 3 shows the MAP weights for the calculation of the horizontal WFPS centroid value, called $W_{k,l}^x$ [the matrix format version of $W_v^x$ in Eq. (2)]. The weights for the vertical centroid value would be the transpose of those listed; that is, $W_{k,l}^y=W_{l,k}^x$. The following parameter values describe the geometry of the system: $D_{\text{sub}}/r_0$ is 2.5; the FWHM of the diffraction spot is 2 pixels; the FoV is $14\times 14\text{pixels}$, and the detector parameters are as in Sec. 4.2. Three different incident light levels have been considered: 30, 100, and 1000, respectively, photons per subpupil. Simulations have been conducted for circular and square subapertures, $D_{\text{sub}}$ being the diameter for the former shape and the side for the latter. Both $Z$-tilt and $G$-tilt have been considered as true tilts. Finally, the spot displacement has been evaluated at the lowermost six spatial frequencies as it has been seen that, for higher frequencies, the weight values become lower than 1% of the total sum of weights.⁹ Table 4 shows the $W_{k,l}^x$ weights for a $D_{\text{sub}}/r_0$ value of unity and an FoV of $10\times 10\text{pixels}$.

Table 3

List of Wk,lx MAP weights for the horizontal WFPS centroid value. Dsub/r0 is 2.5. FWHM of the diffraction spot is 2 pixels. FoV is 14×14  pixels. CIC is 0.05  e−/pixel/frame. RON is 50  e− rms. Values are shown for three incident light levels (30, 100, and 1000 photons), circular, and square subaperture shapes, and Z-tilt and G-tilt estimation, for the case of six selected spatial frequencies.

Table 4

Same as Table 3 for a Dsub/r0 equaling 1 and an FoV of 10×10  pixels.

For the low light level of 30 photons, the weights tend to be similar, irrespective of the subaperture’s shape or the true tilt being estimated. Here, detector noise dominates in the error estimation, and the actual spot’s shape or symmetricity is of much less importance. It is for higher incident light levels that we can see the influence of the spot’s shape on the weights, and in a more evident manner for the $D_{\text{sub}}/r_0=1$ case, for which the spot’s shape in the focal plane is better defined. When estimating $G$-tilt, the weights tend to concentrate at zero frequency, the estimation then approximating to a pure CoG, whereas for the $Z$-tilt estimation, the weights are more dispersed, thereby symmetrizing the spot. Furthermore, weights for a circular subaperture and $G$-tilt estimation tend to be symmetric, whereas both a square subaperture and $Z$-tilt estimation tend to concentrate the weights at zero vertical frequency, especially for $D_{\text{sub}}/r_0=1$ in Table 4.

4.4.

Sensitivity Performance

As can be seen in Sec. 4.2, a centroiding algorithm’s optimal working FoV is a function of the incident light level per subaperture. In practice, it is the available number of pixels per subaperture at the detector that will determine the working FoV in an SHWFS. In this section, we present a comparative behavior of the WFPS centroiding method as a function of incident light for a number of possible FoV values.

Figure 4 shows the plots of tilt estimation error in r.m.s. radians at the sensing wavelength as a function of incident light level per subaperture in photons, with logarithmic scales for both abscissae and ordinates, for three centroiding algorithms: TCoG, plotted in dashed red, CC plotted in dash-dotted green, and WFPS with six spatial frequencies plotted in solid blue. Panel (a) is for $Z$-tilt and (b) is for $G$-tilt estimation. Each panel comprises four subpanels, each for a different FoV value: $10\times 10$, $12\times 12$, $14\times 14$, and $16\times 16\text{pixels}$. $D_{\text{sub}}/r_0$ is 2.5 for a circular subaperture. Sampling follows Nyquist criterion at the diffraction limit. The detector parameters are as in Sec. 4.2.

Fig. 4

(a) $Z$-tilt and (b) $G$-tilt estimation error (r.m.s. radians) versus incident photons per subaperture, for several FoVs and the TCoG, CC, and WFPS centroiding method with six spatial frequencies. Black dotted plots correspond to the WFPS algorithm with 2D-FFT sizes bigger than the FoV. See text for detailed system’s description.

For light levels from below 20 to beyond 200 photons per subaperture, both CC and WFPS algorithms outperform the TCoG method owing to their capacity to neutralize detector noise, and this improvement becomes more evident for increasing FoVs. For example, the decrease in $Z$-tilt and $G$-tilt estimation error at 50 photons is of $\sim 40\%$ and $\sim 33\%$ for a $16\times 16\text{pixel}$ FoV, respectively, and of $\sim 20\%$ and $\sim 12\%$ for a $10\times 10\text{pixel}$ FoV. For high light levels of 1000 to 10000 photons per subaperture, with measurement noise dominated by high-order turbulent phase modes, truncation due to limited FoV, and light shot (intrinsic Poisson) noise, all the algorithms show r.m.s. errors of one twentieth of a wavelength or less. However, certain minor differences between them are worth noting. When estimating $Z$-tilt, which requires symmetrizing the spot, and with a reduced FoV, as in the case of $10\times 10\text{pixels}$, the WFPS method needs to employ bigger 2D-DFTs than the FoV to cope with spots near the edge of the FoV; this is the reason for the black dotted curves, for which $14\times 14$ 2D-FFTs have been used. $G$-tilt estimation, on the other hand, requires considering the asymmetric part of the spot as well. Here, the WFPS method outperforms the CC and TCoG methods for large FoVs because its robustness against the still present detector noise without the need of thresholding.

The straight line of 2.03 r.m.s. radians corresponds to the energy in the tip and tilt modes for a circular aperture with $D/r_0=2.5$ and Kolmogorov turbulence; whereas 0.79 r.m.s. radians is the turbulence level when tip and tilt are corrected, this being a kind of fitting error at the subpupil level.¹⁸ The same sensitivity study has been conducted for the $D_{\text{sub}}/r_0=1$ case, showing the same kind of differences between algorithms, although attenuated. The whole study was also repeated for a square subaperture and the conclusions are identical as for the circular subaperture.⁹

5.

Performance Simulations at Pupil Level

To simulate the performance of the complete SHWFS with the different centroiding algorithms under study, the object-oriented MATLAB toolbox for adaptive optics (OOMAO) has been employed. This is a freely available extension of the MATLAB language consisting of a library of classes oriented toward the numerical modeling of AO systems.¹⁹ The “shack hartmann” object has been for this purpose extended with the WFPS centroiding algorithm; the necessary MAP weights are obtained in the simulation at the subpupil level described in Sec. 4.1. Similarly, the “detector” class, which is embedded in the “shack hartmann” class, has been modified to implement the EMCCD model described in Ref. 16.

Figure 5 shows the simulation workflow corresponding to an open-loop correction configuration. In the upper optical path, an SHWFS receives the light from a Natural Guide Star (NGS) at a wavelength of 550 nm through a 4.2-m telescope, with an 8.4% central obscuration surface. It has an array of $20\times 20$ square lenslets, with an equivalent side size at the telescope aperture of 21 cm, and a spot sampling at diffraction limit that follows the Nyquist criterion. Two values of $r_0$ (21 and 8.4 cm) have been used in the “atmosphere” class to give $D_{\text{sub}}/r_0$ values of unity and 2.5. For the first case, $10\times 10\text{pixels}$ have been assigned at the detector per subaperture, giving an FoV of $\sim {2.7}^{″}$, and $14\times 14\text{pixels}$ have been assigned for the second case, giving an FoV of $\sim {3.8}^{″}$. The QE at the detector has been set to 97%, the EMCCD gain value to 1000, CIC noise to $0.05{\mathrm{e}}^{-}/\text{pixel}/\text{frame}$, and RON to $50\mathrm{rms}{\mathrm{e}}^{-}$. The sensitivity of the detector is the default $1{\mathrm{e}}^{-}/\mathrm{ADU}$. The system sampling rate is 500 frames per second.

Fig. 5

Simulation workflow at the OOMAO.

In the lower optical path, the scientific object under study is located on sky at the same position as the NGS and also at the same wavelength. Its light goes through the telescope to a “zernike” object that acts as a phase corrector, which receives the necessary Zernike coefficients from the “shack hartmann” object. These are calculated in the “shack hartmann” instance in two different ways depending on the incident light level: when the light level is higher than $\sim 50$ photons per subpupil for the TCoG algorithm or $\sim 30$ photons per subpupil for the CC and WFPS algorithms, a zonal linear minimum mean square error method estimates the phases at the subpupils corners, and then a pseudoinversion of the matrix that relates the modes with these phases is applied; for dimmer situations, directly estimating the Zernike modes from the phase slopes, also by matrix pseudoinversion, gives better results. The number of recovered Zernike modes that seems to be best for this lenslet pitch configuration is 120, which corresponds to modes up to the 14th radial order. Phase correction is applied with the same spatial resolution at the pupil as that defined at the “telescope” and “shack hartmann” objects; i.e., the number of lenslets multiplied by the number of detector pixels per lenslet. After this correction, the final Strehl Ratio (SR) can easily be calculated by the Marechal approximation, even before the computation of the point spread function (PSF). The residual phase is integrated during a period of 4 s for this purpose. Finally, an “imager” object integrates the corrected images during 4 s to obtain the PSF over which encircled energy (EE) graphs are obtained.

Panel (a) for a 21 cm $r_0$ and panel (b) for an 8.4 cm $r_0$ in Fig. 6 show the results of the SR achieved as a function of NGS magnitude when employing the TCoG (circle-marked red traces), CC (star-marked green traces), and WFPS with six spatial frequencies (triangle-marked blue traces) centroiding algorithms, all tuned to estimate $Z$-tilt at the subpupil level, and for 120 and 136 recovered Zernike modes (continuous and discontinuous traces, respectively). We can see that the CC and WFPS algorithms show a similar sensitivity behavior, outperforming the TCoG algorithm by increasing the limiting NGS magnitude for a particular SR. This performance improvement occurs for NGS magnitudes from 8.5 to 9 to 11 to 12 and is in accordance with the sensitivity performance at the subpupil level shown in Fig. 4. In the best cases, the NGS limiting magnitude increases by $\sim 0.6$ to $\sim 0.7\mathrm{mag}$. The fraction of encircled PSF energy as a function of the spanned circle diameter size in arcseconds is shown in panels (c) and (d); panel (c) is for a 21 cm $r_0$ and NGS magnitudes of 10.5 (continuous trace) and 11.5 (discontinuous trace), whereas panel (d) is for an 8.4 cm $r_0$ and NGS magnitudes of 9.5 (continuous trace) and 10.5 (discontinuous trace). The centroiding algorithms were tuned to estimate $G$-tilt at the subpupil level and 120 Zernike modes were recovered at the pupil level. EE results are in good agreement with SR results and show the same increase in NGS limiting magnitude of CC/WFPS over TCoG.

Fig. 6

OOMAO simulation results expressed as Strehl ratio as a function of NGS magnitude for (a) $r_0=21\mathrm{cm}$ and (b) $r_0=8.4\mathrm{cm}$, and encircled energy curves as a function of the spanned portion of the PSF in arcsecs for (c) $r_0=21\mathrm{cm}$ and (d) $r_0=8.4\mathrm{cm}$. TCoG results are the red circle-marked traces, CC results are the green star-marked traces and WFPS results are the blue triangle-marked traces.

6.

Laboratory Test

The WFPS algorithm’s performance has been assessed through a laboratory test at the Instituto de Astrofísica de Canarias (IAC), taking advantage of the equalized and diffraction limited field spectrograph experiment (EDiFiSE) project’s AO preoptics. Figure 7 shows the optical setup used. Inside the IAC Atmosphere and Telescope (IACAT) simulator, a variable white light source feeds a pinholed plate through optical fibers; only one has been illuminated, to simulate a single NGS. Atmospheric turbulence is simulated by up to three motorized Kolmogorov phase plates. The 2-cm circular mask that emulates the system aperture has a secondary central obscuration and spider. The characteristics of the William Herschel Telescope have been chosen, with a focal ratio of $f/10.94$. EDiFiSE’s AO subsystem is closed loop configured. However, for the present test, the corrective mirrors from Physik Instrumente (PI) and Alpao are used as passive reflectors, and the SHWFS sees uncorrected turbulence with a large dynamic range. This WFS comprises $10\times 10$ subapertures for a 4.2-m simulated system aperture. A $12\times 12\text{pixel}$ window is assigned to each subaperture at the Andor iXon 860 detector.

Fig. 7

Portion of the EDiFiSE project laboratory optical setup used for the WFPS algorithm assessment.

To assess the performance degradation of the centroiding algorithms as the light level decreases, static turbulent aberrations were simulated; for each one, the light source level was varied, and for each of these levels, SHWFS image sequences of 5000 frames were taken. These sequences were processed offline with the centroiding methods under evaluation: the TCoG, the CC, and the WFPS with the six lowermost spatial frequencies involved and two 2D-FFT sizes, $12\times 12$ and $14\times 14$. The true applied phase was estimated to be at the halfway point in the range between the minimum and maximum centroid coordinate values at each subpupil, computed with all the algorithms over the accumulated frame with the highest light level.

The computation of the WFPS MAP weights and the optimization of the control centroiding methods were aimed to estimate $Z$-tilt, and were performed by simulating the real system’s geometry at the subpupil level, as in Sec. 4. The $r_0$ introduced in the simulation is specified by the phase plates manufacturer for a wavelength of 600 nm, where the Andor detector shows a maximum QE response. Also, a 2.5-pixel diffraction-limited FWHM spot size, obtained from the wavelength value and the optical scale plate at the SHWFS detector, was introduced in the simulation. The detector was characterized and values of $16{\mathrm{e}}^{-}/\mathrm{ADU}$ for the sensitivity, $0.155{\mathrm{e}}^{-}/\text{pixel}/\text{frame}$ for the CIC noise, and 50 rms ${\mathrm{e}}^{-}$ for the RON were obtained and simulated. The EMCCD gain values applied were 800 for up to 1000 incident photons per subpupil and 160 for higher light levels.

Figure 8 shows the estimated error in the phase tilt determination at the SHWFS subpupil in r.m.s. radians as a function of incident light level for (a) $D/r_0\sim 1.5$ and (b) $D/r_0\sim 2$ and the centroiding algorithms under study. The low error values and low dispersion of results at high light level confirm that the estimated true applied phase, though affected by some amount of bias due to spot truncation and high-order aberrations, is adequate for a comparative sensitivity performance assessment among the algorithms under study. Degradation of performance with decreasing light level is in very good coincidence with the simulations at the subpupil level, thus confirming a higher sensitivity of the WFPS over the TCoG algorithm for light levels of 20 to 200 incident photons per subpupil, very similarly to the results for the CC algorithm.

Fig. 8

Estimated error in tilt determination at subpupil level in r.m.s. radians as a function of incident light level for (a) $D/r_0\sim 1.5$ and (b) $D/r_0\sim 2$, for the various centroiding algorithms under study for a real SHWFS at the IAC lab facilities.

7.

Conclusions and Future Work

The WFPS centroiding algorithm has been formulated as a Bayesian estimator of the tilt, at the subpupil of an SHWFS, from the focal plane image Fourier phase slopes. Phase slopes are obtained without phase computation and unwrapping, and the computational cost is estimated as an order of magnitude lower than for the CC algorithm. Its applicability to wavefront sensing with point-like guiding sources has been shown. When large FoVs are required, such as in open-loop sensing and strong turbulence conditions, the WFPS and CC algorithms show a similar sensitivity, which outperforms the TCoG algorithm, for light levels in the range of 20 to 200 incident photons per subpupil and in the presence of detector noise and turbulence high-order perturbations. Sky coverage is thus improved, as the NGS can be 0.6 to 0.7 mag higher.

Work is now under way to test the WFPS algorithm on the sky. For this purpose, simulations should be completed to assess the robustness of the algorithm in changing signal and turbulence strength situations and the necessary MAP weights refreshment period. Also of interest would be the completion of the laboratory test to actually implement a correction open loop and assess the comparative performance of the algorithm in terms of obtained SR, for example. Finally, the algorithm will be implemented in a real-time control platform, and we will look for interested MOAO systems to test it on the sky. Success in this test would encourage us to adapt the algorithm to extended objects sensing, such as in Laser Guide Star and solar AO systems.