1. Introduction

Within the last ten years adaptive optics (AO) correction of astronomical images has become a standard technique on many observatories. For stellar observations the standard wavefront sensing method uses the Shack-Hartmann sensor. The pupil is reimaged onto a lenslet array, thereby splitting the pupil into subpupils. For each lenslet, a subimage of the reference object is produced. If the reference object is a point source, the position of the center of intensity in the subimage defines a vectorial shift that is proportional to the average wavefront gradient over the corresponding lenslet. These vectorial shifts serve as sensor signals for controlling the actuators on a deformable mirror counteracting the wavefront error.

In case the reference object of interest is extended (e.g. the Sun or the Moon), the procedure of measuring the center of intensity must be substituted by cross correlation analysis [1, 2]. The shift of each subimage is estimated relative to a reference subimage for which the shift is estimated relative to the temporally previous subimage. Related to astronomical observations, the main problem is photon starvation resulting in noisy subimage shifts. The effects of noise is well understood for stellar AO [3, 4], and it has also been studied for correlation based wavefront sensing, e.g., in [5, 6].

In this paper we take starting point in a least squares estimate of the shift, expressed in spatial frequency space. The estimate is based on the phase of the cross correlation spectrum. An analytical model for the photon and readout noise induced RMS error to the shift estimate is then derived. In the next section, analytical predictions of how the RMS error depends on object contrast, photon counts and readout noise are verified by simulations, as are the dependencies on sampling format, zero-padding and field of view (FoV). In the last section, the reduction of computational cost for this method compared to the conventional method, based on the localization of the cross correlation peak [6], is discussed.

2. Analytical model

2.1. Shift estimate and noise performance

Given the two images f(r) and g(r), where r = (x,y), and assuming g(r) = f(r - r ₀), a least squares estimate of the shift r ₀ can be obtained by a minimization with respect to r ₀ of the squared error ε ² given by

where F(f) and G(f) are the Fourier transforms of f(r) and g(r) respectively. The squared error ε ² can be minimized setting the derivatives with respect to x ₀ and y ₀ to zero. This leads to

where ℜ extracts the real part and ^* denotes the complex conjugate. Given that φ_C (f) - 2π f·- r ₀ is small, where φ_C (f) is the phase of the cross correlation spectrum C(f) = |C(f) | exp[iφ_C (f)] = F(f)G ^*(f), leads to the following estimate of the shifts x̂ ₀ and ŷ ₀, for details see [7],

Given a sampled version of f(r) and g(r), for which r = (x,y) = Δ(p,q), and using the discrete Fourier transform (DFT), for which f = (f_x , f_y ) = Δ _f (u,v), the shift estimate takes the form

The sampling format m satisfies ΔΔ _f = 1/m, where m × m is the number of pixels (sampling points) in each image. Assuming the two images f(r) and g(r) to be identical but mutually shifted and subjected to additive noise we have F(f) = F ₀(f) + N ₁(f), G(f) = G ₀(f) + N ₂(f) and |F ₀(f)| = |G ₀(f)|. For large signal to noise ratios |F ₀(f)|/|N(f)|, or equivalently |N ₁(f)|/|F ₀(f)| ≪ 1 and |N ₂(f)|/|F ₀(f)| ≪ 1, and a small correlation phase φ_c (f) this leads to

Here φ_{c 0} (f) is the ideal correlation phase and φ_N (f) is the noise contribution. φ_N1 (f) and φ_N2 (f) are the random noise phases and φ_{F 0} (f) and φ_G0 (f) are the phases of the noise-free spectra. Eq. (5) shows that, as expected, a large signal to noise ratio leads to a small noise contribution to the correlation phase and hence also to the shift estimate in Eq. (4). It is therefore important that contributions from spectral components where the signal to noise ratio is small are rejected completely from the shift estimate, as pointed out in references [8, 9].

Assuming 〈|N ₁(f)|²〉 = 〈|N ₂(f)|²〉 = 〈|N(f)|²〉, where 〈〉 stands for ensemble average, and using Eq. (5), it is seen that the sampled version of φ_N (f) satisfies

Replacing φ_C (u,v) in the b-coefficients, last line in Eq. (4), by Φ_N (u,v), the noise contributions δ_x and δ_y to the shift estimates can be calculated. Assuming no correlation of φ_N (u,v) for different (u,v), we get after some calculations

It should be mentioned that since the images are real, the cross correlation spectrum is Her-mitian, i.e., C(u,v) = C ^*(-u, -v), and only half of the points in the spectrum need to be included in Eq. (4). In relation to astronomical AO, where the speed of the calculations is of crucial importance (~ 1 kHz sampling rate), the simplest possible estimate based on only two frequency values (u,v) = (u ₁,0) and (u,v) = (0,v ₁) simplifies the shift estimate in Eq. (4) to

This estimate is independent of the amplitude of the correlation spectrum |C(u,v)|, cf. Eq. (4). It is straightforward to show that this will be the case for all estimates based on only two linearly independent frequency components of the correlation spectrum. As pointed out by [8, 9] the estimate based on u ₁ = v ₁ = 1 has the advantage of being the least sensitive to aliasing. Furthermore, it is an inherent feature of the optical transfer function, that it acts as a spatial low-pass filter leaving the low frequency components in the image with the highest contrast and hence most reliable. Taking u ₁ = v ₁ = 1 it is seen that the condition for avoiding unwrapping the phase φ_C leads to |x ₀| ≤ mΔ/2 and |y ₀| ≤ mΔ/2, or in other words, the shift should be less than half of the field of view. Invoking more than two frequency components as given in Eq. (4) yields the possibility of suppressing contributions from unreliable frequency components having a low |C(u,v)| but still a satisfactory signal to noise ratio. Use of Eqs. (6) and (7) gives the error contributions 〈δ_x 〉^1/2 and 〈δ_y 〉^1/2 to the two point estimate in Eq. (8):

2.2. Performance in the presence of photon and readout noise

Letting f(r) and g(r) be measured in units of photon-counts/arcsec² on the sky, we have

where N _tot = Δ²∑∑f ₀(p,q) is the total photon counts in the noise free image f ₀(p,q) and |V(u,v)| = |F(u,v)|/|F(0,0)| is the visibility, or contrast, of the frequency component (u,v). Assuming the noise contribution to fluctuate independently in different pixels we get

where the index i stands for either 1 or 2. Δ⁴〈$n_i^2$(p,q)〉 = Δ² f ₀(p,q) is valid for photon (Poisson) noise and Δ⁴ 〈$n_i^2$(p, q)〉 = 〈$n_{\mathrm{R}}^2$ (p,q)〉 is valid for readout noise, where 〈$n_{\mathrm{R}}^2$(p, q)〉^1/2 is the readout noise in e-. Introducing these relations in Eq. (11), and assuming uncorrelated noise sources, leads to

Use of Eqs. (9), (10) and (12) gives the result below for the simple two point estimate in Eq. (8):

This is in correspondence with the estimated error given in [3, 4] for a quadrant detector, m = 2, and assuming Nyquist sampling. Several useful conclusions can be drawn from the above expression; (i) the standard error is inversely proportional to the contrast, (ii) neglecting the readout noise, the standard error is proportional to the inverse of the square root of the total photon counts in the subimage, (iii) increasing u ₁ and v ₁ will only lead to improved performance if | V(u ₁, 0) |u ₁ and |V(0, v ₁) |v ₁ are also increased. A few interesting cases can also be considered (neglecting the readout noise).

Case 1: Sampling format. The sampling interval Δ is changed but the field of view (FoV) is fixed. The format m → km, where k < 1 corresponds to less dense sampling and k > 1 corresponds to more dense sampling. The signal N _tot is unchanged, δ → Δ/k and the constant FoV implies that δ _f is unchanged. The continuous image is unchanged and so is the visibility |V(u,v)|. According to Eq. (13) this leads to 〈δ ²〉 = 〈${\delta }_x^2$〉 + 〈${\delta }_y^2$〉 unchanged as long as aliasing can be neglected. If the contribution from the readout noise is included it is seen that the performance degrades with increased sampling.

Case 2: Zero-padding. The FoV, or image, is zero-padded but Δ is fixed. The format m→km, where k≥1. N _tot and Δ are unchanged and Δ _f → Δ _f /k. In case the same physical frequency components are used, that is u→ku and v → kv, this leads to unchanged visibility |V(ku, kv)|. Eq. (13 ) gives that the noise contribution 〈δ ²〉 is unchanged. However, zero-padding leads to the lowest frequency values being diminished, probably leading to higher contrast and a better performance for u ₁=v ₁ = 1 for the zero-padded image.

Case 3: Field of view. The effective FoV is enlarged by increasing the number of pixels which implies that the format m → km. The signal N _tot → k²N_tot, Δ is unchanged and Δ _f → Δ _f /k. The following controlled case may be considered here: The field is enlarged by enlarging the image feature, f(r) → f(r/k), which means that the Fourier transform F(f) → kF(k f). This means that the visibility |V(u, v)| is unchanged and that 〈δ ²〉 is unchanged. However, enlarging the FoV and leaving the overall spectrum unchanged will lead to improved performance for the frequency components showing unchanged visibility. Furthermore the lowest frequency component u ₁ = v ₁ = 1 may show higher contrast further improving the performance.

3. Simulations of the noise performance

3.1. Procedure

The performance of the algorithm was evaluated using two identical, but mutually shifted, images. To estimate a realistic shift between the images, the RMS jitter (in radians) caused by the tip and tilt of the wavefront over a lenslet, given in [10], was used

where the diameter D of the aperture is approximated to the side of a lenslet in the Shack-Hartmann sensor and r ₀ is the Fried parameter (projected onto the lenslet array). We assume that the side of the lenslet is matched to the Fried parameter, D = r ₀. Given that the lenslet image is properly sampled in the focal plane, i.e. Nyquist sampling, the angular image extent of a pixel is

These relationships give that the image RMS jitter is ≈ 1.2 pixels in image space. The reference image, f ₀(p,q), and test image, g ₀(p,q), both scaled to an average intensity level, were shifted such a distance with respect to each other. The fractional pixel shift was obtained by downsampling images from a larger format. Photon and readout noise were then added to the images according to

where ε _P and ε _R are integer random variables described by a zero mean normal probability density function with σ² _εP = f ₀(p, q) and ${\mathrm{\sigma }}_{\epsilon \mathrm{R}}^2$ = 〈$n_{\mathrm{R}}^2$〉 respectively. The signal f ₀(p, q) was assumed to be large, such that the Poisson probability distribution could be approximated with a discrete normal probability distribution. The signal g(p, q) was obtained in the same way.

The DFT’s of f(p, q) and g(p, q) were then calculated and the correlation spectrum is given by C(u,v) =F(u,v)G ^*(u,v). From the correlation spectrum phase, the shift could then be estimated, either by the simple expression given by Eq. (8) or by the amplitude-weighted estimate given by Eq. (4). Also, the conventional method of localization of the cross correlation peak by parabolic interpolation in c(p, q) = ℱ^-1{C(u,v)}, see [6], was used as reference.

As reference image, an ideal negative Gaussian function on a bias level was used to produce a crater-like image, f(r) =b ₁-b ₂ · exp[-r ²/2σ ²]. The parameters b ₁ and b ₂ could be tuned to set the signal and contrast level of the image. The width of the function was set to σ = Δm/6, a sixth of the field of view, in these simulations. Examples of a single noise realisation of the images and the correlation spectrum of this feature are given in Fig. 1. In the same way a more irregular object, cut out from a real image, was used in Fig. 2. In both cases, the average signal per pixel was set to 600 e- and the readout noise was 8 e- per pixel. It is seen in the figures that only the low frequency components exhibit a gradient in the correlation spectrum phase.

 

Fig. 1. Example of images and correlation spectrum for the case of m = 8 and a negative Gaussian function. (a) Reference image f(p,q) with noise (original noise-free high resolution image below). (b) Shifted image g(p,q) with noise (original noise-free high resolution image below). (c) Logarithm of normalized amplitude of correlation spectrum log(|C(u,v)|/|C(0,0)|). (d) Phase of correlation spectrum (|ϕ_C (u,v) in radians.

Download Full Size | PDF

 

Fig. 2. Example of images and correlation spectrum for the case of m = 8 and a more irregular image. See previous figure for explanations.

Download Full Size | PDF

A number of noise realisations could then be used to obtain a measure of the standard error caused by noise. The simulations were also used to investigate bias errors in the estimates, by changing the angle of the separation between the images.

For the case of the ideal Gaussian function, see Fig. 1, the estimates obtained for eight different shift directions, 500 runs for each direction, are given in Fig. 3(a)–(c). In Fig. 3(a), the shift was estimated with Eq. (8), using only the frequency components ϕ_C (1,0) and ϕ_C (0,1). The shift in Fig. 3(b) was estimated with Eq. (4), setting all C(u, v) to zero except for the lowest frequency components C(1, - 1), C(1,0), C(1,1) and C(0,1). In in Fig. 3(c) the shift was estimated using the conventional method of localization of the cross correlation peak by parabolic interpolation, see [6] for details.

In the same way, this was reproduced for the irregular image in Fig. 2. The estimates based on Eq. (8) and ϕ_C (1,0) and ϕ_C (0,1), are shown in Fig. 3(d). The estimates were improved if Eq. (4) was applied, using only the four low frequency components C(1, -1), C(1,0), C(1, 1) and C(0,1). This is seen in Fig. 3(e). Again the shifts given by the the conventional method of localization of the cross correlation peak by parabolic interpolation, are shown in Fig. 3(f).

It is seen that by including more frequency components in the estimate, the precision was improved. Related to Fig. 3(a), the noise contribution 〈${\delta }_x^2$〉^1/2/Δ in the simulation was 0.13 pixels. The parameters used in this simulation were m = 8, N _tot = 38400 e-, |V(1,0)| = 0.053 and 〈$n_{\text{R}}^2$〉^1/2 = 8 e- and the analytical prediction in Eq. (13) leads to the same result. The noise contribution in Fig. 3(b) was 0.10 pixels. Using the conventional shift algorithm, relying on subpixel localization of the cross correlation peak by fitting a quadratic function to the peak and neighbour values, the noise contribution in Fig. 3(c) turned out to be 0.09 pixels.

 

Fig. 3. Estimated shifts from simulations. Each point represents a noise realisation for a separation given by the colors in the legend. Upper row, (a)–(c), represent estimates for the Gaussian function in Fig. 1 and lower row, (d)–(f), represent estimates for the irregular image used in Fig. 2. Leftmost estimates, (a) and (d), are based on Eq. (8) with u ₁ = v ₁ = 1. Middle estimates, (b) and (e), are based on Eq. (4) when setting all C(u,v) to zero but C(1,-1), C(1,0), C(1,1) and C(0,1). Rightmost plots, (c) and (f), give the estimates that are based on parabolic interpolation of the cross correlation peak.

Download Full Size | PDF

3.2. Dependence on contrast, signal level and readout noise

By altering the contrast |V(u,v)|, or the signal level N _tot, or the readout noise 〈$n_{\mathrm{R}}^2$〉^1/2, in the input images, the effect of this on the estimated standard error was compared with the analytical prediction given by Eq. (13). The results are shown in Fig. 4. The values of the parameters used in the simulations were |V(1,0)| = 0.053, N _tot/m ² = 600 e-and 〈$n_{\mathrm{R}}^2$〉^1/2 =0 e-, unless variable, and the image format was m = 8. No shift between the images, with the ideal negative Gaussian object, was introduced in this simulation and 1000 noise realisations were used to estimate each data point in the plots.

 

Fig. 4. Effects of varying parameters on the standard error. Lines - analytical predictions from Eq. (13) using the parameters specified in Section 3.2. Stars - simulated estimation using Eq. (8) with u ₁ = 1. (a) Standard error as function of contrast. (b) Standard error as function of signal. (c) Standard error as function of readout noise.

Download Full Size | PDF

3.3. Dependence on sampling format, zero-padding and field of view

The three cases described in the discussion following Eq. (13) were evaluated by simulations as well. The same images were used for these simulations, the ideal Gaussian image, starting with a sampling format of m = 8. No shift was introduced between the images, and the shift estimate was obtained using Eq. (8). The parameters used were |V(1,0)| = 0.053, N _tot = 38400 e- (unless variable) and 〈$n_{\mathrm{R}}^2$〉^1/2 = 0 e-. 1000 noise realisations were used to estimate each data point.

 

Fig. 5. Lines - analytical predictions from Eq. (13) using the parameters specified in Section 3.3. Stars - simulated estimation using Eq. (8). (a) Standard error, in fraction of the field of view, as a function of sampling points. (b) Standard error, in pixels, as function of zero-padded format. (c) Standard error, in pixels, as function of size of the field of view (stretched object).

Download Full Size | PDF

For case one, i.e. increased sampling within a constant field of view, N _tot was kept constant, whereas the sampling format m was altered to km (the downsampling from the higher format was changed). The estimated shifts were obtained using Eq. (8) with u ₁=v ₁ = 1. The estimated standard error 〈${\delta }_x^2$〉^1/2/Δm, given as a fraction of the field of view, is shown in Fig. 5(a).

For case two, the images were zero-padded before doing the DFT’s, hence extending the format by a factor of k. The estimated shifts were obtained using Eq. (8) with u ₁=v ₁=k. The estimated standard error 〈${\delta }_x^2$〉^1/2/δ, in pixels, is shown in Fig. 5(b).

For case three the same procedure as for case one was used, with the only difference being that the signal was scaled to k ² N _tot. The estimated standard error 〈${\delta }_x^2$〉^1/2/Δ, in pixels, is given in Fig. 5(c). The simulated cases are in agreement with the discussion following Eq. (13). There are no tendencies of variations in 〈${\delta }_x^2$〉^1/2.

4. Gain in computational time

The main reason for investigating this algorithm is the desire for a fast algorithm in order to reduce the servo-lag of the closed loop AO system. Compared to the conventional method of localizing the cross correlation peak by parabolic interpolation, there is a reduction of the computational time when using the method described by Eq. (8).

For the conventional method of fitting a parabolic function to the cross correlation peak there are a number of computations to be carried out. For each shift estimate, or for each subaperture, there is an FFT of the test image, for which the computational time scales as t _F m ²log₂ m ². This is then followed by a multiplication F(u,v)G ^*(u,v) that scales as tMm ². The inverse FFT of the complex cross correlation spectrum is then calculated to obtain the cross correlation function, and the computational time for this scales as 2t _F m ²log₂ m ², where the factor 2 comes from the fact that the correlation spectrum is complex. t_M should be smaller than 2t_F since it only includes a complex multiplication whereas 2t_F also includes a complex addition. The peak value is then to be found and this process scales as t _S m ². The overhead time, t _OH1, consumed by the parabolic interpolation is constant. The total time for the method is given by

As m → ∞ it is seen that this method scales as O(m ² log₂ m ²).

If Eq. (8) is used instead, only the two frequency components G(u ₁,0) and G(0,v ₁) have to be calculated. The computational time for this scales as 2t _M m ². The overhead computational time for the correlation phase and shift estimate is independent of m and is given by t _OH2. The total computational time for the algorithm using Eq. (8) is

and thus scales as O(m ²).

Comparing the computational time of the two algorithms one can see that there is a gain if there is a possibility to use Eq. (8) instead of the conventional method. The gain t _cc/t_φ scales as O(log₂ m ²). As m grows large, it is arguable whether Eq. (8) will give a reasonable estimate. More frequency components could be included to improve the precision, if they show reliable noise performance.

The computational times of the two different algorithms were estimated using a Matlab function written for each routine. The computational speed was optimized as far as possible for both routines. The format m was varied and 10000 runs were carried out for each format to estimate t _cc and t_φ . Matlab 7.0 ran with real time priority under Windows XP on a 1.4 GHz Intel Pentium M processor. Given this non real time operating system, 10 repetitions were carried out to estimate the variation in computational speed.

 

Fig. 6. (a) Time divided by m ² log₂ m ² for conventional method as a function of format. (b) Time divided by m ² for correlation phase method as a function of format. (c) Fraction of the times for the two methods as a function of format.

Download Full Size | PDF

The FFT library, FFTW [11], used in Matlab implies that the computational time for an FFT is not strictly proportional to m ² log₂ m ². The efficiency of the FFTW-library varies with frame format as well as with the processor that the system is running on, see benchmarking in [11].

The measured computational time for a single shift estimate using the traditional cross correlation procedure,t _cc, is normalized with m ² log₂ m ² and plotted in Fig. 6(a) versus frame format. As expected from Eq. (17), the normalized t _cc tends to be constant for large m.

The measured computational time for a single shift estimate using Eq. (8), t_φ , is normalized with m ² and plotted in Fig. 6(b) versus frame format. As expected from Eq. (18), the normalized t_φ tends to be constant for large m. The absolute time for the two methods can be extracted from Fig. 6(a) and Fig. 6(b) respectively, but it is of course computer dependent.

Finally, the fraction between the times of the two methods as a function of frame format is given in Fig. 6(c). The predicted O(log² m ²) is not clearly seen due to the fact that for these small formats there is a considerable overhead time for the methods as seen in Fig. 6(a) and Fig. 6(b). For the format of m = 8 the gain corresponded to approximately seven.

5. Conclusion

A fast algorithm for subpixel shift determination of two images was presented in this paper. Based on only two frequency components in the cross correlation spectrum, the precision of the shift estimate was comparable to the conventional method of localization of the cross correlation peak by parabolic interpolation, for guide objects with a low frequency content. The analytical expression of the precision of the method, depending on noise and image contrast, agreed well with the simulations. This was also valid for the dependence on sampling format, zero-padding and field of view. A reduction of the computational cost was demonstrated when using only two frequency components. This is at the expense of a marginal decrease of performance which could be overcome invoking more frequency components.

The accuracy and precision of this method, similar to all cross correlation methods, will be dependent on the guide object. The remedy to the bias errors in the estimates is to use the algorithm in closed loop AO systems. Furthermore, noncyclic boundaries of the images might necessitate windowing of the images, e.g. with a Hann window, prior to calculation of the cross correlation spectrum.