## Abstract

Scene-based wavefront sensing currently uses the periodic-correlation algorithm based on fast Fourier transforms. However, when the object scene contains features at the field edges, the performance of the algorithm is poor due to the periodicity of fast Fourier transforms, called wraparound effect. In this paper, we propose an algorithm based on the gradient cross-correlation. Both simulation and experiment results show its dramatic effectiveness against the wraparound effect, and a considerable improvement is obtained in image resolution with closed loop adaptive optics correction.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Adaptive Optics (AO) is a technique that works successfully in many areas, including large aperture astronomical telescopes, optical microscopy and laser communications [1–4]. A key component of an AO system is the wavefront sensor. Due to the advantages of simple geometric setup, high processing speed and low cost, Shack-Hartmann Wavefront Sensor (S-H WFS) becomes a popular wavefront sensing instrument in AO system [5, 6]. It consists of a lenslet array and a camera at the focal plane of the lenslet array. The lenslet array can splits the incoming distorted wavefront into many sub-wavefronts for each lenslet, and produces a subimage array of the sensing object on the camera. Normally, all of the sub-wavefront gradients over the corresponding lenslets are measured, and the whole wavefront is reconstructed from those measurements. When the object is a point source, a spot array is formed on the camera. The shift between the centroid position of each spot and a predetermined reference position, which is measured by centroid algorithm, decide the corresponding sub-wavefront gradient. When the object of interest is an extended scene (the extended scene discussed in this paper is within the isoplanatic area), the lenslets form small images of the observed scene and each of these sub-images will be shifted by a small amount just like a point source is. In this case, the centroid algorithm would be invalid, and the procedure of measuring the shifts must be substituted by cross-correlation analysis [7, 8].

There are many situations where extended source AO is desired, such as the observation of the sun during the day [8], imaging conducted along horizontal or slant paths through the atmosphere [9], and earth-observation satellites [10]. Poyneer found that the periodic-correlation algorithm, which employs fast Fourier transforms (FFT), works well in these applications in terms of both performance with noise and computational simplicity [9]. In this approach, each sub-image is treated as a single period of an infinite periodic signal, and the shift between two sub-images is determined from the location of the cross-correlation peak of the two sub-images.

Recently, many improvements have been made for periodic-correlation approach to increase its precision. Sidick et al. proposed an iterative approach, the adaptive periodic-correlation (APC) algorithm [11], which can attain very low errors by using iterations. Townson proposed an algorithm based on parameters optimization using centroid algorithm for sub-pixel shift determination [12]. Despite all those improvements, the periodic-correlation approach still has a serious problem: the wraparound effect. Because of the periodicity of the sub-images, as one end of the sub-image moves away, it wraps around from the other side. For images containing features at the edge, which are common whether in the horizontal observation or in the observation of the Earth from space, the non-common image content at the edges leads to very poor performance, even when the contrast of the images is high.

In this paper, we proposed a gradient cross-correlation algorithm to solve the wraparound effect: calculate the gradient of the sub-images of WFS at first, and then compute the shifts with cross-correlation using gradient images. The gradient images not only provides the same shift information, but also can minimize or even eliminate the wraparound effect when we use them to compute cross-correlation with FFT. For some specific scenes, both digital simulation and data from real experiment showed that the use of gradient cross-correlation algorithm can effectively remove the wraparound effect.

## 2. Periodic-correlation algorithm

#### 2.1 Description of the periodic-correlation algorithm

Consider two sub-images formed by two distinct sub-apertures of WFS (They are both N × N pixels and N is a power-of-2). One is the reference sub-image, $r(x,y)$, which will be used to compare with all the other sub-apertures. The other is the test sub-image, $s(x,y)$.

Assume$s\left(x,y\right)=r\left(x-{x}_{0},y-{y}_{0}\right)$, where the x_{0} and y_{0} are not necessarily integer amounts. The shifts x_{0} and y_{0} can be obtained by finding the peak location of the cross-correlation of the two sub-images. The cross-correlation of the two sub-images is:

^{−1}denotes the inverse FFT. According to the Fourier shift propertyHence the cross-correlation is given by:

_{0}is obtained by parabolic interpolation:

#### 2.2 Wraparound effect

Since the cross-correlation of two finite images is computed using FFT, the finite images have to be periodically extended with periods which divide N pixels. Calculating the correlation of the two sub-images via FFT is exactly a kind of filtering operation, and the reference sub-image is treated as a filter template. When the template sub-image moves relatively on the periodic extension test sub-image, the part of periodic test sub-image, containing periodic extension sections at the edges, is also used for correlation computation. However, this causes serious errors for shifts estimations, especially when the sub-images have features at the edges.

Consider a simplified model: the object scene is a rectangle with bottom linked to the edge of the field, as shown in Fig. 1(a). The rectangle is dark corresponding to intensity 0 and the background is white corresponding to intensity 1. Note that in this paper we are interested in certain kinds of scenes, whose contrast between features and the background is high, and the intensity of the background is relatively uniform. This simplified model is the common case, especially in daytime horizontal atmosphere adaptive optics, whose sources are mountains or buildings. It is selected to simply and clearly explain the wraparound effect and show the good performance of the gradient cross-correlation algorithm for this case.

To illustrate the wraparound effect, we consider a specific scenario: there is only shift along y-axis between the test sub-image and the reference sub-image. As shown in Fig. 1(b) and 1(c), they are parts of the object scene image shown in Fig. 1(a). They are both N × N pixels and N is a power-of-2. The rectangle is longer in the test sub-image than reference sub-image, which represents test sub-image shifts downward relatively to the reference sub-image in the object scene. According to periodic-correlation algorithm, the reference sub-image would move through the 3 × 3 periodic extension test sub-image to estimate the shift between two images, from the coincidence location with the center test sub-image of the extension image as origin. The dynamic range is [-N/2, N/2-1] pixel for x axis direction and [-N/2, N/2-1] pixel for y axis direction. The cross-correlation between the reference sub-image and the corresponding section of the periodic extension test sub-image is calculated once the reference sub-image moves a pixel, where the correlation value is exactly the dot product summation of the two images. The shift corresponding to the maximal cross-correlation value is generally chosen as the estimation of the actual shift between two sub-images.

However, for object scenes contains feature at edges such as the Fig. 1(a), the shift corresponding to the maximal cross-correlation value cannot necessarily be chosen as the estimation of the actual shift anymore. As shown in Fig. 1(d), the black image is the 3 × 3 periodic extension test sub-image. The red and green images indicate two different locations the reference sub-image shift to. The red one coincides with the central test sub-image, and the green one represents the coincidence of the rectangle tops between the central test sub-image and the reference sub-image. On these two locations, the sections of the periodic extension test sub-image participating in cross-correlation calculation are shown in Fig. 1(e) and 1(f), corresponding to green reference sub-image location and red reference sub-image location respectively. It can be seen from Fig. 1(e) that the section contains a small piece of the upper extension test sub-image as the red reference sub-image moved away from the central test sub-image. Apparently, due to the small piece, the cross-correlation results of these two pairs of images are equal and equal to the maximal cross-correlation value, even though the shift amount, which the reference sub-image need for moving to the green section location, is corresponding to the actual shift between reference and test sub-image, not the red section. In this condition, the cross-correlation value can no longer correctly guide the shift estimation, and hence the periodic-correlation algorithm cannot apply.

## 3. Gradient cross-correlation algorithm

Gradient information is originally used in image registration to estimate image motion [13,14]. The gradient of an image not only provides the same shift information as the original image, but also emphasizes the image features. The x-axis direction gradient images ${s}_{G}{}_{x}\left(x,y\right)$ and ${r}_{G}{}_{x}\left(x,y\right)$ of the images $s\left(x,y\right)$and $r\left(x,y\right)$ are defined as Eqs. (7) and (8).

The gradient cross-correlation ${C}_{G}{}_{x}\left(x,y\right)$is then calculated by replacing $s\left(x,y\right)$and $r\left(x,y\right)$ with ${s}_{G}{}_{x}\left(x,y\right)$and ${r}_{G}{}_{x}\left(x,y\right)$ in Eqs. (1) and (2), where ${R}_{G}{}_{x}\left(u,v\right)$ and ${S}_{G}{}_{x}\left(u,v\right)$ are the Fourier transforms of ${r}_{G}{}_{x}\left(x,y\right)$ and ${s}_{G}{}_{x}\left(x,y\right)$ respectively, “*” denotes a complex conjugate and F^{−1} denotes the inverse FFT.

According to the feature of the model object, we employ the y-axis gradient of the sub-images in Figs. 1(b) and 1(c) to calculate the correlation as shown in Figs. 2(b) and 2(c) respectively. The red and green images in Fig. 2(d) have the same location as Fig. 1(d). Since the gradients of the sub-images concentrate the features within the fields, the section in Fig. 2(e) of the periodic extension y-axis gradient test sub-image do not bring features of the upper test sub-image to correlation calculation. It can be seen that the correlation value between the y-axis gradient reference sub-image and the red section of the y-axis gradient periodic extension test sub-image is zero, and the maximal cross-correlation value would be obtained only when the y-axis gradient reference sub-image shift to the green section (where the peak intensity pixels of two images coincide). Therefore, using gradients of sub-images to estimate shifts can definitely obtain high-accuracy results, with only the shifts information of the images delivered. So the wraparound effect is totally avoided.

## 4. Numerical simulation

#### 4.1 Simulation of the shifts between two sub-images

The performance of the gradient cross-correlation algorithm was evaluated using two images. These two images simulate two sub-images of two different sub-apertures of SH-WFS. As shown in Fig. 3, they were generated in the following way: the object scene image shown in Fig. 1(a) is used as the object. The reference sub-image is produced by convolving the object with an ideal spot image generated by Fourier optics. Assuming that there is only tip/tilt phase across each lenslet, the test sub-image is then produced by convolving a known distorted spot image, which is exactly the normalized point spread function of a tip/tilt wavefront. The y-axis gradient sub-images and original sub-images are used for shifts estimation with periodic-correlation algorithm respectively. It can be find easily in the lower-left corner of the Fig. 3 that the top(the blue circle) of the gradient correlation function is exclusive and more noticeable than the one(the blue circle) of the original correlation function.

The original and gradient correlation functions of the object scene shown in Fig. 4(a) are also calculated, in which the background of the object is fiber beam considering the background in practice is not always uniform. The Fig. 4(f) is obtained by employing thresholding processing for the y-axis gradient reference sub-image — the pixels with value less than 40% of the brightest light intensity are set as zero, to reduce the effect of the uneven background in correlation calculation. And it can be obviously found that the characteristics (the red circles) of the sub-images in Figs. 4(f) and 4(g) are more noticeable and therefore the top(the blue circle) of the gradient correlation function in Fig. 4(i) is also more exclusive by comparing with Fig. 4(e).

We processed sets of shifts data of the test sub-image by varying the (peak-to-valley) PV values of the tilt wavefront it convolved with, and excellent performance is obtained using the gradient cross-correlation algorithm. In Fig. 5, the estimate shifts and estimate errors (estimate shift minus theoretical shift) of two algorithms for object without and with fiber beam background are shown. As we can see, in both Figs. 5(a) and 5(b) the estimate shifts calculated by the gradient cross-correlation algorithm are very close to the theoretical lines, with their estimate errors fluctuating around 0 in Figs. 5(c) and 5(d), while the periodic-correlation algorithm is almost invalid for shift estimation. For the object without background, the shifts estimated by periodic-correlation algorithm are all exactly close to 0 whatever the theoretical shifts are, which confirms the prediction. However, due to the participation of the features of the background in the correlation function calculation, the estimate shift differences for the object with background between two algorithms is not as big as that for the object without background.

#### 4.2 Simulation of wavefront reconstruction of S-H WFS

The performance of the gradient cross-correlation algorithm was also evaluated by wavefront reconstruction. The necessary S-H WFS extended scene was generated in the following way: a known Zernike wavefront, such as defocus, is used as the distorted wavefront at the pupil plane of the lenslet array for simplicity. Each lenslet then images the sub-wavefront in front of it as a distorted spot image. Then, a distorted spot image array, as shown in Fig. 6(a), is formed, and convolved with an object scene image as shown in Fig. 6(b). This produced an extended-scene image array similar to the one obtainable with extended-scene SH-WFS, as shown in Fig. 6(c). Each sub-aperture image occupies 28 × 28 pixels and has a field of view (FOV) of 14.9 × 14.9″. The central 16 × 16 pixels are chosen for periodic correlation function calculation because in this paper we mainly consider the situation when the atmospheric turbulence is not very strong.

As shown in Fig. 7, the added distorted wavefront is a defocus wavefront and its PV value is 6λ (λ = 785nm) corresponding to the coefficient of the 4^{th} Zernike mode is 3. Shifts of all the sub-apertures in Fig. 6(c) relative to the reference sub-aperture, which is marked by pink pentagram, are calculated using the periodic-correlation algorithm and the gradient cross-correlation algorithm respectively, and Zernike modal wavefront reconstruction is used. We set threshold for the reference gradient sub-image to reduce the effect of the uneven background in shifts estimation: the pixels with values less than 40% the brightest light intensity are set as zero. As we can see, the reconstructed distorted wavefront and coefficients in Fig. 8(a), which employing the periodic-correlation algorithm, mainly contains not only the defocus but also an error mode (the fifth mode, astigmatism) due to the wraparound effect. Besides, those modes in both figures reconstructed with small coefficients derive from the reconstruction error. The coefficient error of the defocus and the error mode are 1.94 and −1.70 respectively. However, the reconstructed distorted wavefront and coefficients shown in Fig. 8(b) employing the gradient cross-correlation algorithm mainly contain only the defocus. The coefficient of this mode is 3.02, much closer to the added value. This shows the effectiveness of the gradient cross-correlation algorithm for the elimination of the wraparound effect.

## 5. Experimental validation of the gradient cross-correlation algorithm

We perform a series of experiments in laboratory to verify the theoretical analysis, which can be divided into two parts: detection and correction. The optical layout of the experiment is shown in Fig. 9. A xenon lamp with fiber, used as the white light source, is zoomed and illuminates the object through a pair of lens. The object is a piece of glass with a piece of black rectangular paper on it, like the object scene shown in Fig. 1(a). Then the light is transmitted and falls on the liquid-crystal-on-silicon(LCOS) device, which acts as a phase modulator and wavefront corrector. The modulated light is reflected again by a mirror and split into two beams by using a beam splitter (BS), with one beam going to the WFS and the other captured by using an EM CCD camera. The two aperture stops S1 and S2 restrict the width of the beam and the field of each sub-image captured by the S-H WFS respectively. Due to polarization and diffraction of the LCOS, a 700nm-900nm filter F and a polarizer P are placed between L1 and L2. The Lens L3 is used as double pass with input at left part and output at right part, and these two parts are very close to the center of L3. The included angle of the two light beams was within 5°.

The sub-image array collected by the EM CCD of the S-H WFS is shown in Fig. 10. Similar to the simulation in section 4.2, each sub-aperture had 28 × 28 pixels and a FOV of 14.9 × 14.9″. The central 16 × 16 pixels sub-image is used to compute correlation. We chose the 16 × 16 pixels sub-image of the sub-aperture located at center as the reference sub-image, and computed the correlation between it and all other sub-images with and without gradient operation. As we can see, due to the illumination of the xenon lamp with fiber, the bright background of the object scene is not uniform but its contrast is still high. We set threshold for the reference gradient sub-image to reduce the effect of the uneven background in shifts estimation, the pixels with values less than 40% the brightest light intensity are set as zero.

#### 5.1 Experiment of wavefront reconstruction

In the wavefront reconstruction experiment, the LCOS generates a set of Zernike defocus wavefronts, and their PV values are 2λ, 4λ, 6λ and 8λ respectively. λ (785nm) is the central operating wavelength of the LCOS. The performance of the gradient cross-correlation algorithm was evaluated by comparing the reconstructed distorted wavefronts and the residual wavefronts with the periodic correlation algorithm. Same as the simulation, the added wavefront is defocus and the coefficient of the fourth mode is 3 as shown in Fig. 7, which represents that the defocus PV value is 6λ. In Fig. 11(a), the reconstructed wavefront employing periodic-correlation algorithm and the Zernike coefficients are shown. As we can see, the reconstructed wavefront mainly contains defocus and the fifth error mode (astigmatism), and employing gradient sub-images for shifts estimation could obviously reduce the coefficient of the error mode, which highly agrees with the simulation result. Besides, it also can be found that the coefficients of other modes become bigger than them in simulation. This is inevitable in experiment because of the poor quality of the edges of the light beams.

The residual wavefronts could better illustrate how the gradient cross-correlation algorithm works. We can see in Fig. 12 and Fig. 13 that the residual errors employing the gradient cross-correlation algorithm are smaller than those employing the periodic-correlation algorithm. The residual errors were decreased by 28.1% averagely. This agrees with the theoretical analysis and simulation result that employing the gradient cross-correlation algorithm can reduce the wraparound effect errors.

#### 5.2 Experiment of wavefront correction

In the correction experiment, we put a phase plate between L4 and the BS as the aberration which needs to be corrected. The phase plate is a dioptric glass, and it introduced 1.65λ (PV value) defocus into the optical system. The performance of the gradient cross-correlation algorithm was evaluated by comparing the resolution of the images corrected with the periodic-correlation algorithm and the gradient cross-correlation algorithm. The images of object without aberration and with aberration but without correction are shown in Figs. 14(a) and 14(e) respectively. It can be seen that for same aberration, the two algorithms give different results as shown in Figs. 14(b) and 14(f). After correction according to the detected aberrations, the images Figs. 14(c) and 14(g) are all clearer than before while the image corrected using the gradient cross-correlation algorithm is clearer than the other. With almost same small RMS after correction as shown in Figs. 14(d) and 14(h), every fiber beam can be separated from the others distinctly in Fig. 14(g), while in Fig. 14(c) the specifics of the fiber beams and the outline of the object can be vaguely observed. This contradiction shows that the detected result of the periodic-correlation algorithm is inaccurate and emphasizes the effectiveness of the gradient cross-correlation algorithm against the wraparound effect further.

The computational time of the two different algorithms were estimated in a personal computer, with Intel Core i3-4170 CPU (clock speed 3.70GHz) and Windows 7 operating system. The computation time of once wavefront reconstruction was 0.022s of periodic correlation algorithm and 0.124s of gradient cross-correlation algorithm using Matlab R2014a. The gradient correlation algorithm is more accurate but slower. However, due to scene-based wavefront sensor generally collects huge amount of data, hardware acceleration is usually needed for its data processing [15]. FPGA (Field Programmable Gate Array) technology is now an increasingly powerful technology, which can provide amazing processing ability by high speed and parallel processing. Therefore, the gradient correlation algorithm can be accelerated by employing it to a FPGA-based Shack-Hartmann wavefront sensor to satisfy the demand of real-time detection.

## 6. Conclusion

In this paper, a shift estimation algorithm based on the gradient cross-correlation for S-H WFS was presented, well solving the wraparound effect existing in the periodic-correlation algorithm. Agreed well with the theoretical analysis, the gradient cross-correlation algorithm enabled the shift estimation with high precision, and eliminated the error reconstruction mode caused by the wraparound effect in simulation. Agreed well with the simulation results, the residual errors of the reconstructed wavefronts decreased by 28.1% averagely by employing the gradient cross-correlation algorithm in experiment. Close loop AO correction also showed considerable improvement in image resolution by employing the gradient cross-correlation algorithm. Due to its excellent performance against the wraparound effect, it can be a significant wavefront sensing algorithm broadening the range of application of scene-based AO system.

## Funding

National Natural Science Foundation of China (NSFC) (61205021, 61475152, 61405194).

## Acknowledgments

This work is supported by State Key Laboratory of Applied Optics, Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences.

## References and links

**1. **P. Wizinowich, D. S. Acton, C. Shelton, P. Stomski, J. Gathright, K. Ho, W. Lupton, K. Tsubota, O. Lai, C. Max, J. Brase, J. An, K. Avicola, S. Olivier, D. Gavel, B. Macintosh, A. Ghez, and J. Larkin, “First light adaptive optics images from the Keck II telescope: a new era of high angular resolution imagery,” Publ. Astron. Soc. Pac. **112**(769), 315–319 (2000). [CrossRef]

**2. **R. W. Wilson and C. R. Jenkins, “Adaptive optics for astronomy theoretical performance and limitations,” Mon. Not. R. Astron. Soc. **278**(1), 39–61 (1996). [CrossRef]

**3. **M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light Sci. Appl. **3**(4), e165 (2014). [CrossRef]

**4. **R. K. Tyson, “Adaptive optics and ground-to-space laser communications,” Appl. Opt. **35**(19), 3640–3646 (1996). [CrossRef] [PubMed]

**5. **B. C. Platt and R. Shack, “History and principles of Shack-Hartmann wavefront sensing,” J. Refract. Surg. **17**(5), S573–S577 (2001). [PubMed]

**6. **S. Thomas, T. Fusco, A. Tokovinin, M. Nicolle, V. Michau, and G. Rousset, “Comparison of centroid computation algorithms in a Shack-Hartmann sensor,” Mon. Not. R. Astron. Soc. **371**(1), 323–336 (2006). [CrossRef]

**7. **T. R. Rimmele, “Solar adaptive optics,” in *Adaptive Optical Systems Technology*, P. L. Wizinowich, ed., Proc. SPIE 4007**,** 218–231(2000).

**8. **T. R. Rimmele and R. R. Radick, “Solar Adaptive Optics at the National Solar Observatory,” in *Adaptive Optical System Technologies*, D. Bonaccini and R. K. Tyson, eds., Proc. SPIE 3353**,** 72–81 (1998).

**9. **L. A. Poyneer, “Scene-based Shack-Hartmann wave-front sensing: analysis and simulation,” Appl. Opt. **42**(29), 5807–5815 (2003). [CrossRef] [PubMed]

**10. **A. Montmerle Bonnefois, T. Fusco, S. Meimon, L. Mugnier, J. F. Sauvage, C. Engel, C. Escolle, M. Ferrari, E. Hugot, A. Liotard, M. Bernot, M. Carlavan, F. Falzon, T. Bret-Dibat, and D. Laubier, “Comparative theoretical and experimental study of a Shack-Hartmann and a Phase Diversity Sensor, for high-precision wavefront sensing dedicated to Space Active Optics,” Proc. SPIE **10563**, 105634B (2014).

**11. **E. Sidick, J. J. Green, R. M. Morgan, C. M. Ohara, and D. C. Redding, “Adaptive cross-correlation algorithm for extended scene Shack-Hartmann wavefront sensing,” Opt. Lett. **33**(3), 213–215 (2008). [CrossRef] [PubMed]

**12. **M. J. Townson, A. Kellerer, and C. D. Saunter, “Improved shift estimates on extended Shack–Hartmann wavefront sensor images,” Mon. Not. R. Astron. Soc. **452**(4), 4022–4028 (2015). [CrossRef]

**13. **D. I. Barnea and H. F. Silverman, “A class of algorithms for fast digital image registration,” IEEE Trans. Comput. **C-21**(2), 179–186 (1972). [CrossRef]

**14. **J. F. Andrus, C. W. Campbell, and R. R. Jayroe, “Digital image registration algorithm using boundaly maps,” IEEE Trans. Comput. **C-24**(9), 935–940 (1975). [CrossRef]

**15. **X. Peng, M. Li, and C. Rao, “Architecture Design of FPGA-Based Wavefront Processor for Correlating Shack-Hartmann Sensor,” Proc. SPIE **7156**, 71561B (2008). [CrossRef]