The Matched Filter Centroiding (MFC) algorithm measures the centroid location by maximizing the cross correlation of the spot with an assumed reference spot (Leroux & Dainty, 2010). Interpolation is performed on the resultant cross correlation matrix to locate the centroid with sub-pixel accuracy (Poyneer, 2003).

#### **3. Wavefront reconstruction from SHWS slope measurements**

#### **3.1 Vector Matrix Multiply (VMM) method**

The wavefront phase can be obtained by assuming the Southwell's sampling geometry (as shown in Fig. 2) which defines the relation between the wavefront phase and the local slope measurements. Other sampling geometries (Fried, 1977; Hudgin, 1977) can also be adopted instead of Southwell's. From the grid geometry shown in the Fig. 2, it can be shown that (Southwell, 1980) for a SHWS with a pitch of '*h*' and *N* × *N* subapertures,

$$\frac{\mathcal{S}\_{i+1,j}^{x} + \mathcal{S}\_{ij}^{x}}{2} = \frac{\phi\_{i+1,j} - \phi\_{ij}}{h}, \qquad i = 1, 2, \dots \\ N - 1 \qquad \text{and} \qquad j = 1, 2, \dots \\ N$$

$$\frac{\mathcal{S}\_{i+1,j}^{y} + \mathcal{S}\_{ij}^{y}}{2} = \frac{\phi\_{i,j+1} - \phi\_{ij}}{h}, \qquad i = 1, 2, \dots \\ N \qquad \text{and} \qquad j = 1, 2, \dots \\ N - 1 \tag{8}$$

The application of least square curve fitting model, on the SHWS slope data, reduces the problem into a matrix form, *DS* = *Aφ*, where *D* and *A* are sparse matrices of sizes 2*N*2 × 2*N*2 and 2*N*2 × *N*2 respectively. *S* is a vector containing the measured slope values, '*x*' slopes followed by '*y*' slopes. The wavefront phase can therefore be evaluated using the following expression,

$$\phi = (A^{\dagger}A)^{-1}A^{\dagger}DS\tag{9}$$

To solve for *φ* using Matlab, we can apply the in-built command, LSQR on the linear system of equations, *Aφ* = *B*, where *B*(= *DS*) is a vector of size 2*N*2 × 1.

**3.3.1 Temporally evolving phase screens**

(*Pl*

*Pl*

phase screens (*Pl*

evolving phase screens.

are insignificant

Atmospheric statistics are determined by two important parameters, turbulence strength

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 173

velocity profile, *vw*(*h*). These metrics are functions of altitude from the surface of the Earth, '*h*'. Hence, it can be viewed as if the atmosphere is made up of infinite number of turbulence layers. In most situations, it is enough to consider a finite number of discrete layers1 to closely depict the atmosphere. Hence, we limit ourself to a few layers while simulating atmospheric wavefronts. FFT techniques are generally used for the generation of large phase screens in quick time. The limitations like inefficient representation of low frequency turbulence characteristics in the FFT method of phase screen generation can be minimized by the addition

Random wavefronts with turbulence parameters corresponding to different layers can be generated and superposed to simulate a wavefront which closely represents the effect of the full column of the atmosphere. This method of simple superposition of wavefronts is well suited to obtain static and temporally uncorrelated wavefronts. To simulate temporally evolving wavefronts, frozen in turbulence approximation can be used on individual layers for simplicity. A very large wavefront (*Xl* of size *M* × *M* square pixels corresponding to layer *l*) with turbulence strength defined for a single layer may be simulated as a first step. A small portion of this very large wavefront can then be chosen as the initial phase screen

1(*t* = 0), '*t*' represents time) in the temporal evolution of this particular layer. Subsequent

portion on the very large wavefront *Xl* in a definite direction and well defined velocity read from a pre-assumed wind velocity profile (as shown in Fig. 3). The number of such phase

*i*(*t*) ∀ *i* is *N* × *N* square pixels, then the maximum value 'n' can take is equal to (*M* − *N* + 1). The time interval between two temporally adjacent phase screens is defined as: Δ*T* = *dl* / *vl*, where *vl* is the layer velocity and *dl* is the distance moved on the phase screen in a time Δ*T*. In the simulation of temporally evolving phase screens, Δ*T* is kept a constant for all the layers and *dl* (representative of the number of pixels moved within Δ*T* in accordance with the wind velocity of the layer '*l*') is varied for different layers. The phase screens representing the evolution of atmospheric turbulence are finally obtained by superposing different evolving layers. The correlation of the phase screen *P*(*t*) at time *t* = 0 with the phase screen at any other time, *t* = *t*� reduces with increasing *t*� as shown in Fig. 4 for the simulated temporally

The simulation of temporally evolving phase screens described above has a serious problem when a finite number of pixels are used for their representation. The structure of the wind profile would not always allow to move the phase screen by integral number of pixels on the large wavefront, *Xl* in time intervals of Δ*T*. 2 Hence the question of moving the phase screen

is through bilinear interpolation. This method is although simple and well established, it does not completely retain the phase statistics of the simulated wavefronts. Therefore, methods like random mid point displacement (interpolates at the mid point of the pixel) and statistical

by sub-pixel value arises. The most commonly used mathematical tool for moving *Pl*

1 A single atmospheric layer is defined as the depth over which the variations of *C*2

2 *dl* is always not an integral multiple of pixel values

screens generated is limited by the number of pixels on the large wavefront, *Xl*

*i*(*t*), *i*=2,3,...n) at later times are formed by translating the initially selected

*N*(*h*) and vertical wind

. If the size of

*i*(*t*) on *Xl*

*N*(*h*) and *vw*(*h*) profiles

determined by the profile of the refractive index structure constant, *C*2

of low frequency subharmonics (Johansson & Gavel, 1994; Sedmak, 2004).

Fig. 2. Southwell slope geometry (for case of 5 × 5 SHWS) which defines the relation between the wavefront phase and the local slope measurements. The dots represent points where the phase is being estimated and the horizontal and vertical lines over the dots represent the '*x*' and '*y*' slopes measured by the SHWS. The dotted lines show the separation between the lenslets.

#### **3.2 Fast Fourier Transform (FFT) based reconstructor**

This is a modal reconstruction process where complex exponentials are used as basis functions and the wavefront can be reconstructed from its local discrete slopes by a simple multiplicative filtering operation in the spatial frequency domain (Freischlad & Koliopoulos, 1986). Taking the case of Hudgin geometry where the first differences, of the phase values, are the measured slope values,

$$\begin{aligned} \mathcal{S}\_{i,j}^{\chi} &= \phi\_{i+1,j} - \phi\_{i,j} \\ \mathcal{S}\_{i,j}^{\mathcal{Y}} &= \phi\_{i,j+1} - \phi\_{i,j} \end{aligned} \tag{10}$$

Applying the Fourier transform and using the shift property of the discrete Fourier transform on Eq. 10 gives,

$$\begin{aligned} \mathbf{S}\_{k,l}^x &= \Phi\_{k,l} \left[ \exp\left(\frac{\hat{\imath}2\pi k}{N}\right) - 1 \right] \\ \mathbf{S}\_{k,l}^y &= \Phi\_{k,l} \left[ \exp\left(\frac{\hat{\imath}2\pi l}{N}\right) - 1 \right] \end{aligned} \tag{11}$$

where **S** and Φ represent the Fourier transforms of vectors containing slopes and phase values respectively. Φ can now be solved and an inverse Fourier transform performed to arrive at the wavefront phase profile, *φ*. The suitability of different geometries for FFT technique is studied in comparison with the VMM method (Correia et al., 2008).

#### **3.3 Monte Carlo simulations**

Monte Carlo simulations are used for testing adaptive optics systems. It involves the generation of atmosphere like turbulence phase screens, simulation of the SHWS spot pattern and the retrieval of the wavefront shape using centroiding and wavefront reconstruction algorithms. The performance of the SHWS and the wavefront sensing algorithm is quantified by calculating the correlation coefficient between the simulated (initial) phase screen and the reconstructed wavefront.

#### **3.3.1 Temporally evolving phase screens**

6 Adaptive Optics

Fig. 2. Southwell slope geometry (for case of 5 × 5 SHWS) which defines the relation between the wavefront phase and the local slope measurements. The dots represent points where the phase is being estimated and the horizontal and vertical lines over the dots represent the '*x*' and '*y*' slopes measured by the SHWS. The dotted lines show the separation

*Sx*

*Sy*

**S***x k*,*l* = Φ*k*,*l*

**S***y k*,*l* = Φ*k*,*l*

in comparison with the VMM method (Correia et al., 2008).

This is a modal reconstruction process where complex exponentials are used as basis functions and the wavefront can be reconstructed from its local discrete slopes by a simple multiplicative filtering operation in the spatial frequency domain (Freischlad & Koliopoulos, 1986). Taking the case of Hudgin geometry where the first differences, of the phase values, are the measured

*i*,*j* = *φi*+1,*j* − *φi*,*j*

Applying the Fourier transform and using the shift property of the discrete Fourier transform

where **S** and Φ represent the Fourier transforms of vectors containing slopes and phase values respectively. Φ can now be solved and an inverse Fourier transform performed to arrive at the wavefront phase profile, *φ*. The suitability of different geometries for FFT technique is studied

Monte Carlo simulations are used for testing adaptive optics systems. It involves the generation of atmosphere like turbulence phase screens, simulation of the SHWS spot pattern and the retrieval of the wavefront shape using centroiding and wavefront reconstruction algorithms. The performance of the SHWS and the wavefront sensing algorithm is quantified by calculating the correlation coefficient between the simulated (initial) phase screen and the

 ˆ*i*2*πk N*

> ˆ*i*2*πl N*

 − 1

 − 1

 *exp*

> *exp*

*i*,*j* = *φi*,*j*+1 − *φi*,*j* (10)

(11)

between the lenslets.

slope values,

on Eq. 10 gives,

**3.3 Monte Carlo simulations**

reconstructed wavefront.

**3.2 Fast Fourier Transform (FFT) based reconstructor**

Atmospheric statistics are determined by two important parameters, turbulence strength determined by the profile of the refractive index structure constant, *C*2 *N*(*h*) and vertical wind velocity profile, *vw*(*h*). These metrics are functions of altitude from the surface of the Earth, '*h*'. Hence, it can be viewed as if the atmosphere is made up of infinite number of turbulence layers. In most situations, it is enough to consider a finite number of discrete layers1 to closely depict the atmosphere. Hence, we limit ourself to a few layers while simulating atmospheric wavefronts. FFT techniques are generally used for the generation of large phase screens in quick time. The limitations like inefficient representation of low frequency turbulence characteristics in the FFT method of phase screen generation can be minimized by the addition of low frequency subharmonics (Johansson & Gavel, 1994; Sedmak, 2004).

Random wavefronts with turbulence parameters corresponding to different layers can be generated and superposed to simulate a wavefront which closely represents the effect of the full column of the atmosphere. This method of simple superposition of wavefronts is well suited to obtain static and temporally uncorrelated wavefronts. To simulate temporally evolving wavefronts, frozen in turbulence approximation can be used on individual layers for simplicity. A very large wavefront (*Xl* of size *M* × *M* square pixels corresponding to layer *l*) with turbulence strength defined for a single layer may be simulated as a first step. A small portion of this very large wavefront can then be chosen as the initial phase screen (*Pl* 1(*t* = 0), '*t*' represents time) in the temporal evolution of this particular layer. Subsequent phase screens (*Pl i*(*t*), *i*=2,3,...n) at later times are formed by translating the initially selected portion on the very large wavefront *Xl* in a definite direction and well defined velocity read from a pre-assumed wind velocity profile (as shown in Fig. 3). The number of such phase screens generated is limited by the number of pixels on the large wavefront, *Xl* . If the size of *Pl i*(*t*) ∀ *i* is *N* × *N* square pixels, then the maximum value 'n' can take is equal to (*M* − *N* + 1). The time interval between two temporally adjacent phase screens is defined as: Δ*T* = *dl* / *vl*, where *vl* is the layer velocity and *dl* is the distance moved on the phase screen in a time Δ*T*. In the simulation of temporally evolving phase screens, Δ*T* is kept a constant for all the layers and *dl* (representative of the number of pixels moved within Δ*T* in accordance with the wind velocity of the layer '*l*') is varied for different layers. The phase screens representing the evolution of atmospheric turbulence are finally obtained by superposing different evolving layers. The correlation of the phase screen *P*(*t*) at time *t* = 0 with the phase screen at any other time, *t* = *t*� reduces with increasing *t*� as shown in Fig. 4 for the simulated temporally evolving phase screens.

The simulation of temporally evolving phase screens described above has a serious problem when a finite number of pixels are used for their representation. The structure of the wind profile would not always allow to move the phase screen by integral number of pixels on the large wavefront, *Xl* in time intervals of Δ*T*. 2 Hence the question of moving the phase screen by sub-pixel value arises. The most commonly used mathematical tool for moving *Pl i*(*t*) on *Xl* is through bilinear interpolation. This method is although simple and well established, it does not completely retain the phase statistics of the simulated wavefronts. Therefore, methods like random mid point displacement (interpolates at the mid point of the pixel) and statistical

1 A single atmospheric layer is defined as the depth over which the variations of *C*2 *N*(*h*) and *vw*(*h*) profiles are insignificant

2 *dl* is always not an integral multiple of pixel values

in precisely estimating the wavefront shape depends dominantly on the number of detector pixels corresponding to a single subaperture (*NP*), the detector readout noise and its quantum efficiency (Hardy, 1998). As '*NS*' increases, the wavefront sensing accuracy increases at the cost of increased detector readout time and additional computational effort. Increasing the subaperture size will enhance the number of collected photons and in turn reduce the wavefront sensor error. The subaperture size cannot be increased beyond the Fried parameter, '*r*0' so that the angular resolution is not limited by the spatial scale of turbulence, but by the subaperture size. The focal length of the microlens array determines the maximum detectable tilt range of each subaperture. To increase the dynamic range of the SH sensor, detectors with large '*NP*' could be selected, but at the cost of reduced wavefront sampling frequency. There can be a significant degree of error due to ill-positioning of the imaging device at the focal plane of the microlenses. The experimental limits to the performance are light intensity levels,

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 175

This section describes a working model of the low cost, simple laboratory SHWS at the AO Laboratory, Indian Institute of Astrophysics. A CCD camera (model Pulnix TM-1325CL, readout noise = 50 *e*−, pitch = 6.45*μ*m) placed at the focal plane of an array of microlenses (from Flexible Optical B.V., pitch = 200*μ*m, *f* = 4 cm) makes our wavefront sensor. A continuous membrane deformable mirror (Multi-DM from Boston Micromachines, 140 actuators, pitch = 450*μ*m, maximum stroke = 5.5*μ*m) was used to generate distorted wavefronts. He-Ne laser (15mW, *λ*=632.8*nm*) was used as the light source. Also, 4-f geometry was used to remove high frequency noise. The number of pixels per subaperture is 31 in our case and for sensing the whole DM aperture, we used 25 × 25 microlenses. Low resolution higher order Zernike polynomials (matrix of size: 12 × 12) are simulated and corresponding voltage values to be addressed on the DM actuators are calculated. The maximum voltage value that is applied to a single DM actuator is 40V, in order to avoid crosstalk between the

*n* = 10 and azimuthal index, *m* = 4). The maximum gray scale in the image corresponds to 40V and minimum gray scale is 0V. The CCD image of the SHWS spot pattern that was

5b. The wavefronts reconstructed using CoG and WCoG methods are shown in Figs. 5c and 5d respectively. Readout noise does not significantly affect wavefront reconstruction accuracy in our wavefront sensing experiment. This was tested by performing multiple trials of reconstructing the same wavefront distortion introduced on the DM. We find that the percentage error in the wavefront reconstruction accuracy due to readout noise is ∼0.24%. The existence of a good light level makes CoG a better option over WCoG (Roopashree et al.,

This section discusses the problems associated with the centroiding techniques described in the earlier sections and brings out advanced methods to deal with those issues in the case of astronomical sensing. The following subsections discuss the algorithms: improved iteratively

*m* with radial index,

4 on the DM is shown in Fig.

Fig. 5a shows the simulation of a low resolution Zernike polynomial (*Zn*

captured by addressing the simulated Zernike polynomial, *Z*10

exposure time scales and operating speeds.

adjacent subapertures of the SHWS.

**6. Advanced centroiding algorithms**

**5. Working model**

2011).

Fig. 3. Vertical wind profile used for simulation of multilayered evolving turbulence

Fig. 4. Simulation of temporal evolution of turbulence based on a seven layer wind model - Drop in the correlation coefficient with time

interpolation (Wu et al., 2009) have been devised. Statistical interpolation is preferred over random mid point displacement method since it can interpolate to any fraction of pixel pitch. A comparison of the performance of statistical interpolation with other interpolation techniques in simulating evolving phase screens is reported by Roopashree et al. (2010).

#### **3.3.2 Metric for comparison**

Wavefront reconstruction accuracy is quantified through the correlation coefficient (**CC**) which is used as a metric for comparing the simulated phase screen and the reconstructed one,

$$\mathbf{CC} = \frac{E(XY) - E(X).E(Y)}{\sigma\_X \sigma\_Y} \tag{12}$$

where *E* represents the expectation value, *X* and *Y* are vectors containing the pixel values of the phase profiles under comparison and *σX*, *σY* represents the standard deviations.

#### **4. Optimizing the SHWS design parameters**

The accuracy of wavefront sensing depends largely on the design parameters of the SHWS, which are highly interdependent and they include, the number of subapertures (*NS*), subaperture size (d) and the focal length of the microlenses (*f*). The ability of the SHWS in precisely estimating the wavefront shape depends dominantly on the number of detector pixels corresponding to a single subaperture (*NP*), the detector readout noise and its quantum efficiency (Hardy, 1998). As '*NS*' increases, the wavefront sensing accuracy increases at the cost of increased detector readout time and additional computational effort. Increasing the subaperture size will enhance the number of collected photons and in turn reduce the wavefront sensor error. The subaperture size cannot be increased beyond the Fried parameter, '*r*0' so that the angular resolution is not limited by the spatial scale of turbulence, but by the subaperture size. The focal length of the microlens array determines the maximum detectable tilt range of each subaperture. To increase the dynamic range of the SH sensor, detectors with large '*NP*' could be selected, but at the cost of reduced wavefront sampling frequency. There can be a significant degree of error due to ill-positioning of the imaging device at the focal plane of the microlenses. The experimental limits to the performance are light intensity levels, exposure time scales and operating speeds.

### **5. Working model**

8 Adaptive Optics

**0 0.5 1 1.5 2 2.5**

**Altitude (km)**

**0 0.1 0.2 0.3 0.4 0.5 −0.2**

Fig. 4. Simulation of temporal evolution of turbulence based on a seven layer wind model -

interpolation (Wu et al., 2009) have been devised. Statistical interpolation is preferred over random mid point displacement method since it can interpolate to any fraction of pixel pitch. A comparison of the performance of statistical interpolation with other interpolation techniques in simulating evolving phase screens is reported by Roopashree et al. (2010).

Wavefront reconstruction accuracy is quantified through the correlation coefficient (**CC**) which is used as a metric for comparing the simulated phase screen and the reconstructed

> **CC** = *E*(*XY*) − *E*(*X*).*E*(*Y*) *σXσY*

where *E* represents the expectation value, *X* and *Y* are vectors containing the pixel values of

The accuracy of wavefront sensing depends largely on the design parameters of the SHWS, which are highly interdependent and they include, the number of subapertures (*NS*), subaperture size (d) and the focal length of the microlenses (*f*). The ability of the SHWS

the phase profiles under comparison and *σX*, *σY* represents the standard deviations.

**Time (sec)**

Fig. 3. Vertical wind profile used for simulation of multilayered evolving turbulence

**5**

**0 0.2 0.4 0.6 0.8 1 1.2**

**Temporal correlation coefficient**

Drop in the correlation coefficient with time

**4. Optimizing the SHWS design parameters**

**3.3.2 Metric for comparison**

one,

**10**

**15**

**Wind speed (m/s)**

**20**

**25**

**x 104**

(12)

This section describes a working model of the low cost, simple laboratory SHWS at the AO Laboratory, Indian Institute of Astrophysics. A CCD camera (model Pulnix TM-1325CL, readout noise = 50 *e*−, pitch = 6.45*μ*m) placed at the focal plane of an array of microlenses (from Flexible Optical B.V., pitch = 200*μ*m, *f* = 4 cm) makes our wavefront sensor. A continuous membrane deformable mirror (Multi-DM from Boston Micromachines, 140 actuators, pitch = 450*μ*m, maximum stroke = 5.5*μ*m) was used to generate distorted wavefronts. He-Ne laser (15mW, *λ*=632.8*nm*) was used as the light source. Also, 4-f geometry was used to remove high frequency noise. The number of pixels per subaperture is 31 in our case and for sensing the whole DM aperture, we used 25 × 25 microlenses. Low resolution higher order Zernike polynomials (matrix of size: 12 × 12) are simulated and corresponding voltage values to be addressed on the DM actuators are calculated. The maximum voltage value that is applied to a single DM actuator is 40V, in order to avoid crosstalk between the adjacent subapertures of the SHWS.

Fig. 5a shows the simulation of a low resolution Zernike polynomial (*Zn m* with radial index, *n* = 10 and azimuthal index, *m* = 4). The maximum gray scale in the image corresponds to 40V and minimum gray scale is 0V. The CCD image of the SHWS spot pattern that was captured by addressing the simulated Zernike polynomial, *Z*10 4 on the DM is shown in Fig. 5b. The wavefronts reconstructed using CoG and WCoG methods are shown in Figs. 5c and 5d respectively. Readout noise does not significantly affect wavefront reconstruction accuracy in our wavefront sensing experiment. This was tested by performing multiple trials of reconstructing the same wavefront distortion introduced on the DM. We find that the percentage error in the wavefront reconstruction accuracy due to readout noise is ∼0.24%. The existence of a good light level makes CoG a better option over WCoG (Roopashree et al., 2011).

#### **6. Advanced centroiding algorithms**

This section discusses the problems associated with the centroiding techniques described in the earlier sections and brings out advanced methods to deal with those issues in the case of astronomical sensing. The following subsections discuss the algorithms: improved iteratively

Fig. 6. Laser guide star spot pattern as seen at the focal plane of the lenslet array

**60 80 100 120**

**Altitude**

1.2025 in bimodal profile arises due to the unit normalization of

; single modal profile

�� ; bimodal profile (14)

**0 0.2 0.4 0.6 0.8 1**

− (h−94.5)2 2.352

**Vertical Density of**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 177

Fig. 7. Single mode and bimodal vertical line profiles of the column of Sodium layer in the

altitude of 90 km. It was also experimentally observed that there exists two peaks instead of one and hence it is called as bimodal vertical sodium density profile (Drummond et al., 2004). The line profile can be approximated with a single mode (Gaussian distribution) or modeled as if it is bimodal (linear combination of two Gaussian functions) in nature as shown in Fig 7.

Here, "h" represents the altitude from sea level and (*h*) is the unit normalized density as a

density profile when h ∈ [65 115]. The single mode is centered on 90 km with FWHM of 10 km. The bimodal profile has two peaks centered on 84 km and 94.5 km with FWHM of 8.24 km and 2.35 km respectively. The effect of temporally varying sodium layer profile as shown in Fig. 8 can be included in the simulations of elongated laser guide star spots. To generate randomly varying single modal sodium layer profile, as a first iteration, pseudo random number picked from a normal distribution is added to the mean height of the profile (*h*=90 km) in Eq. 14. In the subsequent iterations, '*i*' pseudo random numbers are added to the earlier iteration value,

**Sodium Layer (Normalized)**

**60 80 100 120**

**Altitude**

**0 0.2 0.4 0.6 0.8 1**

mesosphere. Altitude in km units

(*h*) =

⎧ ⎨ ⎩ exp �

1 1.2025 �

function of "h". A factor of 1

− (h−90)2 102

exp �

�

− (h−84)2 8.242

�

+ exp �

**Vertical Density of**

**Sodium Layer (Normalized)**

Fig. 5. (a) Simulated low resolution Zernike polynomial, *Z*10 4 addressed on the DM (b) Spot pattern on the SHWS: *I* = 0.28*μ*W at the location of the photodetector, SLM in off state (c) Reconstructed phase using CoG, *c* = 0.733 (d) using WCoG, *c* = 0.718.

weighted center of gravity, thresholded Zernike reconstructor based centroiding and Gaussian pattern matching algorithm.

#### **6.1 Improved iteratively weighted center of gravity**

The advent of large telescopes made the sensing process more challenging by posing multiple problems due to the use of the LGS. To make it even more complicated, the field of view can only be widened by making use of more than one LGS. The finite thickness of the sodium layer makes it difficult to generate spots with an artificial source that is point-like. The situation becomes worse while dealing with large aperture telescopes (≥10 m). The observed spots are elongated with the elongation length being a function of the distance of the spot from the laser launch point and the elongation axis in the direction of the line joining the subaperture center and the launch point. The spot elongation (*�*) can be estimated by,

$$
\epsilon \approx \frac{L\delta H}{H^2 + H.\delta H} \tag{13}
$$

where *L* is the physical distance from the subaperture center to the launch point of LGS, *δH* is the sodium layer thickness and *H* is the altitude of the sodium layer. The result from the simulation of an elongated spot pattern (no atmosphere) with the launch point matching with the center of the telescope aperture is shown in Fig 6. Another problem along with the non uniform elongation is the temporal variability of the vertical sodium density profile (Davis et al., 2006). Hence, an advanced image processing tool is necessary for accurate detection of the centroid position in the case of LGS based Shack Hartmann sensor (Schreiber et al., 2009; Thomas et al., 2008).

Any iterative process like the IWCoG carries along with it problems like - saturation of performance, non uniform convergence and slow convergence, which can be confirmed from a detailed Monte Carlo analysis (Vyas et al., 2010b). The density of sodium peaks nearly at an 10 Adaptive Optics

**(a) (b)**

**(c) (d)**

pattern on the SHWS: *I* = 0.28*μ*W at the location of the photodetector, SLM in off state (c)

weighted center of gravity, thresholded Zernike reconstructor based centroiding and Gaussian

The advent of large telescopes made the sensing process more challenging by posing multiple problems due to the use of the LGS. To make it even more complicated, the field of view can only be widened by making use of more than one LGS. The finite thickness of the sodium layer makes it difficult to generate spots with an artificial source that is point-like. The situation becomes worse while dealing with large aperture telescopes (≥10 m). The observed spots are elongated with the elongation length being a function of the distance of the spot from the laser launch point and the elongation axis in the direction of the line joining the subaperture center

*�* ≈ *LδH*

where *L* is the physical distance from the subaperture center to the launch point of LGS, *δH* is the sodium layer thickness and *H* is the altitude of the sodium layer. The result from the simulation of an elongated spot pattern (no atmosphere) with the launch point matching with the center of the telescope aperture is shown in Fig 6. Another problem along with the non uniform elongation is the temporal variability of the vertical sodium density profile (Davis et al., 2006). Hence, an advanced image processing tool is necessary for accurate detection of the centroid position in the case of LGS based Shack Hartmann sensor

Any iterative process like the IWCoG carries along with it problems like - saturation of performance, non uniform convergence and slow convergence, which can be confirmed from a detailed Monte Carlo analysis (Vyas et al., 2010b). The density of sodium peaks nearly at an

4 addressed on the DM (b) Spot

*H*2 + *H*.*δH* (13)

Fig. 5. (a) Simulated low resolution Zernike polynomial, *Z*10

**6.1 Improved iteratively weighted center of gravity**

(Schreiber et al., 2009; Thomas et al., 2008).

pattern matching algorithm.

Reconstructed phase using CoG, *c* = 0.733 (d) using WCoG, *c* = 0.718.

and the launch point. The spot elongation (*�*) can be estimated by,

Fig. 6. Laser guide star spot pattern as seen at the focal plane of the lenslet array

Fig. 7. Single mode and bimodal vertical line profiles of the column of Sodium layer in the mesosphere. Altitude in km units

altitude of 90 km. It was also experimentally observed that there exists two peaks instead of one and hence it is called as bimodal vertical sodium density profile (Drummond et al., 2004). The line profile can be approximated with a single mode (Gaussian distribution) or modeled as if it is bimodal (linear combination of two Gaussian functions) in nature as shown in Fig 7.

$$\varrho(h) = \begin{cases} \exp\left[-\frac{(\mathbf{h} - 90)^2}{10^2}\right] & \text{; single modal profile} \\ \frac{1}{1.2025} \left[ \exp\left(-\frac{(\mathbf{h} - 84)^2}{8.24^2}\right) + \exp\left(-\frac{(\mathbf{h} - 94.5)^2}{2.35^2}\right) \right] & \text{; bimodal profile} \end{cases} \tag{14}$$

Here, "h" represents the altitude from sea level and (*h*) is the unit normalized density as a function of "h". A factor of 1 1.2025 in bimodal profile arises due to the unit normalization of density profile when h ∈ [65 115]. The single mode is centered on 90 km with FWHM of 10 km. The bimodal profile has two peaks centered on 84 km and 94.5 km with FWHM of 8.24 km and 2.35 km respectively. The effect of temporally varying sodium layer profile as shown in Fig. 8 can be included in the simulations of elongated laser guide star spots. To generate randomly varying single modal sodium layer profile, as a first iteration, pseudo random number picked from a normal distribution is added to the mean height of the profile (*h*=90 km) in Eq. 14. In the subsequent iterations, '*i*' pseudo random numbers are added to the earlier iteration value,

**0 5 10 15 20**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 179

**SNR = 0.2 SNR = 0.8 SNR = 1.4 SNR = 2.0**

**Iteration Number**

where (*xAct*, *yAct*) is the actual centroid position of the spot that is although unknown in real

We expect that increasing the number of iterations would reduce the CEE which does not happen in the case of iterative algorithms in general. The error tends to saturate after a finite number of iterations. In the IWCoG algorithm, in most cases, the CEE saturates within a few (3-10) iterations as shown in Fig 10. This error saturation effect can prove to become a serious problem since the CEE cannot be reduced further at the cost of increased computational time

The saturation of CEE in IWCoG algorithm can be circumvented by implementing a simple and effective technique called the Iterative Addition of Random Numbers (IARN). After the *nth* iteration, the new weighting function for the (*n* + 1)*th* iteration is calculated such that it is

*R*, *yn*

*c* + *Nxrn*

variance, which assumes different values after each iteration. *Nx* and *Ny* are normalization constants (crudely they define the magnitude of randomness introduced; in our case *Nx* = *Ny* = *N*). The addition of these pseudo random numbers would not allow the saturation of CEE. This process although creates other problems like irregular convergence or no convergence as shown in Fig. 11, the minimum CEE that can be reached reduces by a large amount. This

It can also be observed from Fig. 10 that at different SNR, the iteration number with minimum CEE is different. Due to the non uniform convergence as seen in Fig. 11, the CEE after the last iteration may not be the minimum CEE among all the iterations. For SNR of 0.2, the minimum occurs at n=10; and for SNR=2, minimum CEE occurs at n=5. This suggests that the optimum iteration number is a function of SNR. Hence, we need a technique to identify the iteration

*R*)=(*xn*

*R*) instead of (*xn*

*y* are the generated pseudo random numbers with zero mean and unit

*c* + *Nyrn*

*x* , *yn*

*c* , *yn*

*c* ) such that,

*y* ) (16)

Fig. 10. Saturation and non uniform convergence of IWCoG algorithm at different SNR

**0.2 0.4 0.6 0.8 1 1.2**

and effort, as a result leading to redundant calculations.

(*xn R* , *yn*

method is more significant at low SNR as depicted in Fig. 12.

centered on a modified centroid location, (*xn*

situations, is predefined while performing Monte Carlo simulations.

**Centroid Estimation**

**6.1.1 Error saturation**

where, *rn*

*x* and *rn*

**6.1.2 Non uniform convergence**

with minimum error.

 **Error**

Fig. 8. The fluctuations in the altitude of the sodium layer

Fig. 9. Spot patterns (using single mode vertical sodium profile) with different SNRs (a)∞ (b) 1 (c) 2 (d) 3

'*hi*−1'. In the case of bimodal profile, the peaks of the two Gaussian functions are modified by adding pseudo random numbers, similar to the case of single modal profile.

Instead of simulating the entire adaptive optics loop, we concentrate on simulating the case of centroiding the noisy spot pattern corresponding to a single subaperture. A single subaperture spot pattern is simulated by convolving the vertical sodium layer density profile (length proportional to the elongation length defined in Eqn. 13 with atmosphere like phase screen (Schreiber et al., 2009). The simulated phase screen includes the effect of the round trip that the laser has undertaken and the phase screens can be easily simulated based on a Fourier technique (Harding et al., 1999). The phase screens simulated are based on spatial and temporal superposition of multiple layers of the atmosphere (Roopashree et al., 2010).

Once the spot pattern is simulated, photon noise and readout noise are added. Poisson noise is added to the simulated spot pattern by replacing intensity values at individual pixels by pseudo Poisson random numbers picked up from a Poisson distribution generated with mean equal to the intensity value at those pixels. Readout noise is added by generating *M*2 (as many as the number of pixels in the spot) pseudo Gaussian random numbers (zero mean and standard deviation defined by SNR) and adding to individual intensity values in the spot pattern. The resultant spot pattern images for different SNR are shown in Fig 9.

The Centroid Estimation Error (CEE) can be used as a performance metric in determining the accuracy of centroiding in the Monte Carlo simulations. It is defined as,

$$CEE = \sqrt{(\mathbf{x}\_{Act} - \mathbf{x}\_s)^2 + (y\_{Act} - y\_s)^2} \tag{15}$$

Fig. 10. Saturation and non uniform convergence of IWCoG algorithm at different SNR

where (*xAct*, *yAct*) is the actual centroid position of the spot that is although unknown in real situations, is predefined while performing Monte Carlo simulations.

#### **6.1.1 Error saturation**

12 Adaptive Optics

**0 50 100 150 200**

**Trail Number**

Fig. 9. Spot patterns (using single mode vertical sodium profile) with different SNRs (a)∞ (b)

'*hi*−1'. In the case of bimodal profile, the peaks of the two Gaussian functions are modified by

Instead of simulating the entire adaptive optics loop, we concentrate on simulating the case of centroiding the noisy spot pattern corresponding to a single subaperture. A single subaperture spot pattern is simulated by convolving the vertical sodium layer density profile (length proportional to the elongation length defined in Eqn. 13 with atmosphere like phase screen (Schreiber et al., 2009). The simulated phase screen includes the effect of the round trip that the laser has undertaken and the phase screens can be easily simulated based on a Fourier technique (Harding et al., 1999). The phase screens simulated are based on spatial and temporal superposition of multiple layers of the atmosphere (Roopashree et al., 2010).

Once the spot pattern is simulated, photon noise and readout noise are added. Poisson noise is added to the simulated spot pattern by replacing intensity values at individual pixels by pseudo Poisson random numbers picked up from a Poisson distribution generated with mean equal to the intensity value at those pixels. Readout noise is added by generating *M*2 (as many as the number of pixels in the spot) pseudo Gaussian random numbers (zero mean and standard deviation defined by SNR) and adding to individual intensity values in the spot

The Centroid Estimation Error (CEE) can be used as a performance metric in determining the

(*xAct* − *xs*)2 + (*yAct* − *ys*)2 (15)

pattern. The resultant spot pattern images for different SNR are shown in Fig 9.

accuracy of centroiding in the Monte Carlo simulations. It is defined as,

*CEE* =

adding pseudo random numbers, similar to the case of single modal profile.

**89**

Fig. 8. The fluctuations in the altitude of the sodium layer

**90**

**91**

**89.5**

**90.5**

**Altitude in km**

1 (c) 2 (d) 3

We expect that increasing the number of iterations would reduce the CEE which does not happen in the case of iterative algorithms in general. The error tends to saturate after a finite number of iterations. In the IWCoG algorithm, in most cases, the CEE saturates within a few (3-10) iterations as shown in Fig 10. This error saturation effect can prove to become a serious problem since the CEE cannot be reduced further at the cost of increased computational time and effort, as a result leading to redundant calculations.

The saturation of CEE in IWCoG algorithm can be circumvented by implementing a simple and effective technique called the Iterative Addition of Random Numbers (IARN). After the *nth* iteration, the new weighting function for the (*n* + 1)*th* iteration is calculated such that it is centered on a modified centroid location, (*xn R*, *yn R*) instead of (*xn c* , *yn c* ) such that,

$$(\mathbf{x}\_R^n \; \prime \; y\_R^n) = (\mathbf{x}\_c^n + \mathbf{N}\_\mathbf{x} r\_\mathbf{x}^n \; \prime \; y\_c^n + \mathbf{N}\_\mathbf{y} r\_\mathbf{y}^n) \tag{16}$$

where, *rn x* and *rn y* are the generated pseudo random numbers with zero mean and unit variance, which assumes different values after each iteration. *Nx* and *Ny* are normalization constants (crudely they define the magnitude of randomness introduced; in our case *Nx* = *Ny* = *N*). The addition of these pseudo random numbers would not allow the saturation of CEE. This process although creates other problems like irregular convergence or no convergence as shown in Fig. 11, the minimum CEE that can be reached reduces by a large amount. This method is more significant at low SNR as depicted in Fig. 12.

#### **6.1.2 Non uniform convergence**

It can also be observed from Fig. 10 that at different SNR, the iteration number with minimum CEE is different. Due to the non uniform convergence as seen in Fig. 11, the CEE after the last iteration may not be the minimum CEE among all the iterations. For SNR of 0.2, the minimum occurs at n=10; and for SNR=2, minimum CEE occurs at n=5. This suggests that the optimum iteration number is a function of SNR. Hence, we need a technique to identify the iteration with minimum error.

**0 20 40 60 80 100**

**5 10 15 20 25 30**

**After IARN Before IARN**

**Trial Number**

Fig. 14. Improved CEE after application of modified IARN based IWCoG algorithm

The maximum number of iterations (*nmax*) used although increases the accuracy in some cases, the number of iterations it needs to go below a predefined error threshold or minimum error sets a limit on the maximum speed of computation. IWCoG algorithm with *nmax >* 1 is the slowest among all the other centroiding techniques discussed here. The need to compute large number of iterations for cases of low SNR causes additional burden on the processor. An improved speed of convergence reduces the number of iterations to be performed and hence

The time taken to compute the centroid position of a single subaperture in the lenslet array as a function of the number of iterations is shown in Fig 15. The Monte Carlo simulations were

performed on a 1.4GHz Intel(R) Core(TM)2 Solo CPU with 2GB RAM.

**Error Propagation in IARN Calculated Correlation Coefficients**

**for iteration selection**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 181

**Iteration Number**

Fig. 13. Error propagation in the case of a single noisy LGS spot pattern. There is a negative correlation between the CEE and the correlation coefficient calculated by comparing the noisy spot pattern with a ideal spot pattern positioned at the new centroid estimated in a particular iteration of IWCoG algorithm, i.e. dips in CEE correspond to peaks in the

**0**

**0**

**0.05**

**0.1**

**0.15**

**0.2**

**1**

**2**

**0.5**

**1.5**

**Centroid Estimation Error**

correlation coefficient

**Minimum Centroid**

**6.1.3 Speed of convergence**

improves speed.

**Estimation Error**

 **from 100 iterations**

Fig. 11. The effect of error saturation can be overcome by using the IARN based IWCoG algorithm. This procedure largely reduces the CEE.

Fig. 12. IARN based IWCoG algorithm. Monte Carlo simulations were performed for 1000 different spot patterns at different SNR and the CEE is taken as a minimum among 1000 iterations.

A simple technique based on iterative computation of correlation can be used to identify the iteration with minimum error. After the *nth* iteration, an ideal spot pattern (simulated using PSF with atmospheric turbulence excluded and with zero noise) is calculated with its center at (*xn c* , *yn c* ) and compared with the actual noisy spot pattern image using the correlation coefficient as the metric for comparison. There exists a strong negative correlation between this metric and the CEE as shown in Fig 13. At this point we make a hypothesis that the iteration number with minimum error corresponds to the one with maximum correlation coefficient. However, there exists a finite decorrelation between this metric and CEE which cannot be always neglected. To thoroughly validate the obtained result, the centroid location with minimum error is estimated multiple number of times based on this hypothesis.

This modified version of IARN based IWCoG method is a consistent performer. Simulations show that there is a significant improvement in the obtained CEE after the application of this modified IARN method on IWCoG as against simple IWCoG (Fig 14).

Fig. 13. Error propagation in the case of a single noisy LGS spot pattern. There is a negative correlation between the CEE and the correlation coefficient calculated by comparing the noisy spot pattern with a ideal spot pattern positioned at the new centroid estimated in a particular iteration of IWCoG algorithm, i.e. dips in CEE correspond to peaks in the correlation coefficient

Fig. 14. Improved CEE after application of modified IARN based IWCoG algorithm

#### **6.1.3 Speed of convergence**

14 Adaptive Optics

**Effect of Error Saturation After using IARN**

> **Before IARN After IARN**

**0 20 40 60 80 100**

**Number of Iterations**

**1 2 3 4 5**

Fig. 12. IARN based IWCoG algorithm. Monte Carlo simulations were performed for 1000 different spot patterns at different SNR and the CEE is taken as a minimum among 1000

A simple technique based on iterative computation of correlation can be used to identify the iteration with minimum error. After the *nth* iteration, an ideal spot pattern (simulated using PSF with atmospheric turbulence excluded and with zero noise) is calculated with its

coefficient as the metric for comparison. There exists a strong negative correlation between this metric and the CEE as shown in Fig 13. At this point we make a hypothesis that the iteration number with minimum error corresponds to the one with maximum correlation coefficient. However, there exists a finite decorrelation between this metric and CEE which cannot be always neglected. To thoroughly validate the obtained result, the centroid location

This modified version of IARN based IWCoG method is a consistent performer. Simulations show that there is a significant improvement in the obtained CEE after the application of this

with minimum error is estimated multiple number of times based on this hypothesis.

modified IARN method on IWCoG as against simple IWCoG (Fig 14).

*c* ) and compared with the actual noisy spot pattern image using the correlation

**Signal to Noise Ratio (RON)**

Fig. 11. The effect of error saturation can be overcome by using the IARN based IWCoG

**0**

**0 0.1 0.2 0.3 0.4**

**Centroid Estimation Error**

iterations.

center at (*xn*

*c* , *yn*

algorithm. This procedure largely reduces the CEE.

**0.2**

**0.4**

**Centroid Estimation Error**

**0.6**

**0.8**

The maximum number of iterations (*nmax*) used although increases the accuracy in some cases, the number of iterations it needs to go below a predefined error threshold or minimum error sets a limit on the maximum speed of computation. IWCoG algorithm with *nmax >* 1 is the slowest among all the other centroiding techniques discussed here. The need to compute large number of iterations for cases of low SNR causes additional burden on the processor. An improved speed of convergence reduces the number of iterations to be performed and hence improves speed.

The time taken to compute the centroid position of a single subaperture in the lenslet array as a function of the number of iterations is shown in Fig 15. The Monte Carlo simulations were performed on a 1.4GHz Intel(R) Core(TM)2 Solo CPU with 2GB RAM.

**5 10 15 20**

**Signal to Noise Ratio**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 183

**0 0.2 0.4 0.6 0.8 1**

**Signal to Noise Ratio**

centroid position and large *N* introduces large variations. Simulation results shown here use N=1, i.e. the centroid is moved within a pixel. The disadvantage of the correlation based identification of iteration with minimum error is that it is time-consuming. Also, the dependence of this minimum error iteration detection technique has to be investigated further

In this section, a centroiding technique based on thresholded Zernike reconstructor is presented. This technique relies on removing the high spatial frequency noise present in the spots by reconstructing them via the calculation of lower order complex Zernike moments. The reconstructed spots are thresholded to remove unwanted noise features and the position

of the centroid is then estimated using the weighted centroiding algorithm.

Fig. 18. The consistency of IARN method at very low SNR makes it a more reliable

Fig. 17. A comparison of IARN based IWCoG with other algorithms. It can be observed that IARN is most appropriate at SNR ≤ 5. At high SNR, simple IWCoG performs better. A proper choice of "*N*" may allow IARN based IWCoG method to dominate at higher SNR also.

**0.2**

**0**

**0.2**

for its performance at different SNR.

**6.2 Thresholded Zernike reconstructor**

**0.4**

**Centroid Estimation**

centroiding technique

 **Error**

**0.6**

**0.8**

**0.4**

**Centroid Estimation Error**

**0.6**

**0.8**

**CoG IWC IWCoG IARN**

**CoG IWC IWCoG IARN**

Fig. 15. Time taken for computations. The graph looks like a staircase. This is because of the discrete measurements made by the computer using MATLAB R2008a.

Fig. 16. Increased speed of convergence at low iteration number

A simple and efficient technique can improve the speed significantly. In the first (initial) iteration, the weighting function is centered on the point of maximum intensity (*XI*−*Max* , *YI*−*Max*) in the spot pattern image. The improved speed of convergence at low iteration number in this faster IWCoG as against regular IWCoG method is shown in Fig 16 (Monte Carlo simulations performed using 1000 spot pattern images for SNR = 1). Although it seems as if only the first few (three) iterations converge faster using this modified IWCoG, for nearly 10% of the spot images the convergence improves by atleast five iterations.

A comparison of the improved IWCoG algorithm based on IARN is shown in Figs 17 and 18 for high and low SNR. It can be observed that the performance of improved IWCoG based on IARN is better at low SNR (≤ 5).

The novelty of this improved IARN based IWCoG algorithm is its consistency and constancy of CEE irrespective of the SNR. In addition, it is highly accurate and fast centroiding technique. The normalization constant *N* plays an important role and is a subject of further numerical understanding. Small value of *N* introduces small variations in the estimated 16 Adaptive Optics

**Data Points Linear Fit**

discrete measurements made by the computer using MATLAB R2008a.

**0 200 400<PHONE_NUMBER>**

**Modified IWCoG IWCoG**

**Number of Iterations**

**0 5 10 15**

**Iteration Number**

A simple and efficient technique can improve the speed significantly. In the first (initial) iteration, the weighting function is centered on the point of maximum intensity (*XI*−*Max* , *YI*−*Max*) in the spot pattern image. The improved speed of convergence at low iteration number in this faster IWCoG as against regular IWCoG method is shown in Fig 16 (Monte Carlo simulations performed using 1000 spot pattern images for SNR = 1). Although it seems as if only the first few (three) iterations converge faster using this modified IWCoG,

A comparison of the improved IWCoG algorithm based on IARN is shown in Figs 17 and 18 for high and low SNR. It can be observed that the performance of improved IWCoG based on

The novelty of this improved IARN based IWCoG algorithm is its consistency and constancy of CEE irrespective of the SNR. In addition, it is highly accurate and fast centroiding technique. The normalization constant *N* plays an important role and is a subject of further numerical understanding. Small value of *N* introduces small variations in the estimated

for nearly 10% of the spot images the convergence improves by atleast five iterations.

Fig. 15. Time taken for computations. The graph looks like a staircase. This is because of the

**0**

**0**

Fig. 16. Increased speed of convergence at low iteration number

**0.2**

**0.4**

**Centroid Estimation Error**

IARN is better at low SNR (≤ 5).

**0.6**

**0.8**

**0.02 0.04 0.06 0.08 0.1 0.12**

**Time (sec)**

Fig. 17. A comparison of IARN based IWCoG with other algorithms. It can be observed that IARN is most appropriate at SNR ≤ 5. At high SNR, simple IWCoG performs better. A proper choice of "*N*" may allow IARN based IWCoG method to dominate at higher SNR also.

Fig. 18. The consistency of IARN method at very low SNR makes it a more reliable centroiding technique

centroid position and large *N* introduces large variations. Simulation results shown here use N=1, i.e. the centroid is moved within a pixel. The disadvantage of the correlation based identification of iteration with minimum error is that it is time-consuming. Also, the dependence of this minimum error iteration detection technique has to be investigated further for its performance at different SNR.

#### **6.2 Thresholded Zernike reconstructor**

In this section, a centroiding technique based on thresholded Zernike reconstructor is presented. This technique relies on removing the high spatial frequency noise present in the spots by reconstructing them via the calculation of lower order complex Zernike moments. The reconstructed spots are thresholded to remove unwanted noise features and the position of the centroid is then estimated using the weighted centroiding algorithm.

Fig. 20. The first row contains spots (32 × 32 pixels) with SNR = 1,2,3,4. The second row contains spots reconstructed using Zernike moments with *n* ≤ 24. The third row contains the

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 185

Fig. 21. Left most spot pattern is the noisy spot of size 30 × 30 pixels with SNR = 1, the reconstructed spot is shown on its right using n ≤ 29 and thresholded spots (65, 70, 75%) are

where '*a*' and '*b*' are the major and minor radii of the ellipse. The elongation of the LGS spot has little effect on the centroiding estimation error after application of TZR as shown in Fig. 23. When compared against smaller spots or streaks, the performance of the algorithm is better in

The effect of spot orientation on centroid estimation accuracy is shown in Fig. 24. The estimate for the centroid location is best when the spots are oriented along the direction of the pixels (at 0 orientation angle). The performance of the algorithm can be improved by using a polar coordinate detector where the major axis of each rectangular pixel is aligned with the axis of

The temporal fluctuations in the vertical sodium layer density profile is simulated and shown in Fig. 25. The centroiding estimation error in the case of fluctuating profile is shown in Fig. 25. It can be observed that the fluctuations in the sodium profile nearly doubles the centroid estimation error. It is hence important to continuously monitor the fluctuations and

thresholded spots with intensity greater than 60% of maximum intensity.

shown in the same order.

the case of larger ones.

elongation (Thomas et al., 2008).

Fig. 19. Spot corresponding to a single subaperture with SNR = 1,2,3,4,5 from left to right.

The parameters of interest include the spot elongation, spot orientation (depends on the position of the laser launch point with respect to the subaperture), the amount of shift in the spot, spot size in non-elongated direction and the SNR. The sample spots at different SNR are shown in Fig. 19.

The proposed method of centroiding is based on image reconstruction using Zernike polynomials. A detailed description of the mathematical procedure for the calculation of Zernike moments for individual sub aperture images is given in Vyas et al. (2010c). The major parameters of interest include number of Zernike moments used for reconstruction and the thresholding percentage. Once the spot is reconstructed using Zernike polynomials, it is ready to be thresholded. A threshold percentage is chosen in the first place. Generally preferred between 50 and 90%; anything lower or higher would either retain noise features after thresholding or remove pixels that are very essential for centroiding purposes. Individual pixels on single subaperture spots are treated as 'object pixels' if the value of the pixel is greater than the threshold. Imposing an intensity threshold on the spot pattern allows us to remove most of the high spatial frequency events (or noise features). A sample spot image at different SNR that is reconstructed and threholded is shown in Fig 20. Also, an improper choice of the thresholding percentage can lead to erroneous features in the reconstructed spot images as shown in Fig. 21.

The Zernike reconstructed and thresholded spots are then subjected to weighted centroiding. In the case of NGS, a Gaussian weighting function is assumed and for the case of LGS, the weighting function used is a standard reference of the spot pattern without noise, by assuming the most probabilistic simple sodium layer profile. The thresholded Zernike reconstructor based centroiding has been numerically tested for the case of NGS and LGS based sensing (Vyas et al., 2010a;c).

Based on the Monte Carlo simulations, we could see a significant improvement in the centroid estimation accuracy when thresholded Zernike reconstructor is applied in conjugation with the weighted centroiding technique. A comparison of this technique with simple center of gravity (CoG), weighted center of gravity (WCoG), iterative addition of random numbers (IARN) based iteratively weighted center of gravity (Baker & Moallem, 2007; Vyas et al., 2010b) (IWCoG) is shown in Fig. 22. It can be observed that TZR algorithm performs equally in comparison with other algorithms at SNR between 1 and 2. At a SNR less than 1, IWCoG performs better than TZR. When TZR is combined with WCoG instead of CoG, it performs better than any other algorithm at low SNR.

The elongation (e) of the LGS spot pattern is assumed to be equal to the elongation of the equivalent closest ellipse.

$$e = \frac{a - b}{a} \tag{17}$$

18 Adaptive Optics

Fig. 19. Spot corresponding to a single subaperture with SNR = 1,2,3,4,5 from left to right.

are shown in Fig. 19.

images as shown in Fig. 21.

(Vyas et al., 2010a;c).

equivalent closest ellipse.

better than any other algorithm at low SNR.

The parameters of interest include the spot elongation, spot orientation (depends on the position of the laser launch point with respect to the subaperture), the amount of shift in the spot, spot size in non-elongated direction and the SNR. The sample spots at different SNR

The proposed method of centroiding is based on image reconstruction using Zernike polynomials. A detailed description of the mathematical procedure for the calculation of Zernike moments for individual sub aperture images is given in Vyas et al. (2010c). The major parameters of interest include number of Zernike moments used for reconstruction and the thresholding percentage. Once the spot is reconstructed using Zernike polynomials, it is ready to be thresholded. A threshold percentage is chosen in the first place. Generally preferred between 50 and 90%; anything lower or higher would either retain noise features after thresholding or remove pixels that are very essential for centroiding purposes. Individual pixels on single subaperture spots are treated as 'object pixels' if the value of the pixel is greater than the threshold. Imposing an intensity threshold on the spot pattern allows us to remove most of the high spatial frequency events (or noise features). A sample spot image at different SNR that is reconstructed and threholded is shown in Fig 20. Also, an improper choice of the thresholding percentage can lead to erroneous features in the reconstructed spot

The Zernike reconstructed and thresholded spots are then subjected to weighted centroiding. In the case of NGS, a Gaussian weighting function is assumed and for the case of LGS, the weighting function used is a standard reference of the spot pattern without noise, by assuming the most probabilistic simple sodium layer profile. The thresholded Zernike reconstructor based centroiding has been numerically tested for the case of NGS and LGS based sensing

Based on the Monte Carlo simulations, we could see a significant improvement in the centroid estimation accuracy when thresholded Zernike reconstructor is applied in conjugation with the weighted centroiding technique. A comparison of this technique with simple center of gravity (CoG), weighted center of gravity (WCoG), iterative addition of random numbers (IARN) based iteratively weighted center of gravity (Baker & Moallem, 2007; Vyas et al., 2010b) (IWCoG) is shown in Fig. 22. It can be observed that TZR algorithm performs equally in comparison with other algorithms at SNR between 1 and 2. At a SNR less than 1, IWCoG performs better than TZR. When TZR is combined with WCoG instead of CoG, it performs

The elongation (e) of the LGS spot pattern is assumed to be equal to the elongation of the

*e* = *a* − *b a*

(17)

Fig. 20. The first row contains spots (32 × 32 pixels) with SNR = 1,2,3,4. The second row contains spots reconstructed using Zernike moments with *n* ≤ 24. The third row contains the thresholded spots with intensity greater than 60% of maximum intensity.

Fig. 21. Left most spot pattern is the noisy spot of size 30 × 30 pixels with SNR = 1, the reconstructed spot is shown on its right using n ≤ 29 and thresholded spots (65, 70, 75%) are shown in the same order.

where '*a*' and '*b*' are the major and minor radii of the ellipse. The elongation of the LGS spot has little effect on the centroiding estimation error after application of TZR as shown in Fig. 23. When compared against smaller spots or streaks, the performance of the algorithm is better in the case of larger ones.

The effect of spot orientation on centroid estimation accuracy is shown in Fig. 24. The estimate for the centroid location is best when the spots are oriented along the direction of the pixels (at 0 orientation angle). The performance of the algorithm can be improved by using a polar coordinate detector where the major axis of each rectangular pixel is aligned with the axis of elongation (Thomas et al., 2008).

The temporal fluctuations in the vertical sodium layer density profile is simulated and shown in Fig. 25. The centroiding estimation error in the case of fluctuating profile is shown in Fig. 25. It can be observed that the fluctuations in the sodium profile nearly doubles the centroid estimation error. It is hence important to continuously monitor the fluctuations and

**0 2 4 6 8 10**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 187

**Signal to Noise Ratio**

**with sodium layer height fluctuations without fluctuations in sodium layer profile mean centroid estimation error; with fluctuations mean centroid estimation error; with no fluctuations**

**0 20 40 60 80 100 120 140 160 180 200**

**Trail Number**

Fig. 25. A comparison of the performance of centroiding accuracy with and without sodium layer fluctuations in the case of TZR method. Since the weighting function was not modified

fluctuations exist. Hence, it becomes important to track the profile fluctuations continuously.

The number of pixels used for a single subaperture can also influence centroiding estimation while using TZR. This is because image reconstruction in the case of smaller images becomes easier with fewer Zernike moments. The effect of changing image size is shown in Fig. 27. At a low SNR, reconstructing a 24 × 24 spot pattern using n ≤ 29 leads to inefficient representation of all the required spot features, whereas representation of a 12 × 12 spot with the same number of Zernike modes is more efficient. Hence, in this case (n ≤ 29 and 60% thresholding) at low SNR, smaller spots must be preferred. On the other hand, at a high SNR, fewer Zernike modes are sufficient to represent all the features of the spot. Centroiding accuracy is good for

with changing sodium layer profile, the mean centroiding accuracy reduces when

Fig. 24. Spot orientation affects the centroiding estimation accuracy. It is best when the

streaks are oriented in the direction of the sides of the detector pixels.

**0**

**1**

**2**

**1.25 1.5 1.75**

**Centroiding estimation error**

**0.25 0.5 0.75**

**0**

larger spot pattern at high SNR.

**0.1**

**0.2**

**0.3**

**Centroid estimation error**

**0.4**

**0.5**

Fig. 22. A comparison of different centroiding algorithms with the TZR method in conjugation with WCoG. It can be seen that at low SNR, TZR performs much better than any other centroiding technique. The performance is also stable. Although when applied alone without WCoG, TZR is not better than IWCoG.

Fig. 23. TZR method performs well at different elongation lengths.

make corrections to the weighting function used in the centroiding of the spots. Numerical observations suggest that by doing so, the centroiding estimation error reduces by little above 10%, but does not reach the value that is obtained with zero fluctuations.

Using a very few Zernike moments for reconstruction weakens the spot of its details and having too many of them include the higher order effects which can be identified with noise. As the image size increases, it becomes important to include more and more Zernike terms to closely reconstruct the features of the spot. Percentage thresholding has significant effect on centroid estimation as can be seen in Fig. 26. Centroid estimation error dropped with reducing thresholding percentage up to 60%. Reducing the thresholding percentage further is risky and would retain noise features which may lead to large errors in the case of a spot of size 30 × 30 pixels with 29 Zernike modes.

20 Adaptive Optics

**TZR IWCoG**

**TZR+WCoG CoG**

**IARN based IWCoG**

**0.0067 0.2260 0.3181**

**0 1 2 3 4**

conjugation with WCoG. It can be seen that at low SNR, TZR performs much better than any other centroiding technique. The performance is also stable. Although when applied alone

**0 2 4 6 8 10**

**Signal to Noise Ratio**

make corrections to the weighting function used in the centroiding of the spots. Numerical observations suggest that by doing so, the centroiding estimation error reduces by little above

Using a very few Zernike moments for reconstruction weakens the spot of its details and having too many of them include the higher order effects which can be identified with noise. As the image size increases, it becomes important to include more and more Zernike terms to closely reconstruct the features of the spot. Percentage thresholding has significant effect on centroid estimation as can be seen in Fig. 26. Centroid estimation error dropped with reducing thresholding percentage up to 60%. Reducing the thresholding percentage further is risky and would retain noise features which may lead to large errors in the case of a spot of size 30 × 30

Fig. 22. A comparison of different centroiding algorithms with the TZR method in

**Signal to Noise Ratio**

**0 0.5 1 1.5 2 2.5 3 3.5 4**

without WCoG, TZR is not better than IWCoG.

**2**

**0**

Fig. 23. TZR method performs well at different elongation lengths.

10%, but does not reach the value that is obtained with zero fluctuations.

**0.5**

**1**

**Centroiding estimation error**

pixels with 29 Zernike modes.

**1.5**

**Centroiding estimation error**

Fig. 24. Spot orientation affects the centroiding estimation accuracy. It is best when the streaks are oriented in the direction of the sides of the detector pixels.

Fig. 25. A comparison of the performance of centroiding accuracy with and without sodium layer fluctuations in the case of TZR method. Since the weighting function was not modified with changing sodium layer profile, the mean centroiding accuracy reduces when fluctuations exist. Hence, it becomes important to track the profile fluctuations continuously.

The number of pixels used for a single subaperture can also influence centroiding estimation while using TZR. This is because image reconstruction in the case of smaller images becomes easier with fewer Zernike moments. The effect of changing image size is shown in Fig. 27. At a low SNR, reconstructing a 24 × 24 spot pattern using n ≤ 29 leads to inefficient representation of all the required spot features, whereas representation of a 12 × 12 spot with the same number of Zernike modes is more efficient. Hence, in this case (n ≤ 29 and 60% thresholding) at low SNR, smaller spots must be preferred. On the other hand, at a high SNR, fewer Zernike modes are sufficient to represent all the features of the spot. Centroiding accuracy is good for larger spot pattern at high SNR.

**6.3 Gaussian pattern matching**

identification.

**6.3.1 Feature recognition**

**6.3.2 Shape identification**

**6.3.3 Profile identification**

after the denoising procedure is implemented.

At high noise conditions, the thresholded Zernike reconstructed images lead to multiple noisy features along with the actual spot pattern as depicted in Figs. 20 and 21. This is due to the fact that at high noise level conditions there can be large scale features (sometimes the features can be scaled to dimensions comparable to that of the spot) which might not be removed even

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 189

To use this erroneous spot for accurate centroiding, we take the advantage of the fact that in the case of NGS, the spot finally formed at the focal plane of a lens must assume an airy pattern which can be approximated to a Gaussian like structure. The pattern matching algorithm is implementing in three steps: feature recognition, shape identification and profile

In this step, the features on the spot pattern image are counted and identified. Feature recognition can be performed by using many existing pattern recognition algorithms. In this algorithm, we used a simple Hough peak identification method to detect the features and number them in the order of peak height. To eliminate small scale features which may arise due to unavoidable scintillations, we impose threshold conditions on the size of the features.

Most features don't have a circular shape. The circularity or the extent of the feature being circular is measured for each of the features. This can be measured in many ways. In the proposed algorithm, the local centroid for the feature is calculated and the distance from the centroid at which intensity becomes zero is measured. This distance is called the distance parameter. The distance parameter is estimated at different angles (0-360) from the feature centroid position. We define circularity, C of a single feature as the inverse of the variance of the distance parameter computed over different angles from the centroid position. For an ideal circular feature, the circularity is infinity since the variance is zero. A lower cutoff for

In the previous step, the shape of the spot was used for selective elimination. In this step, the intensity profile is used to choose the actual spot feature. The number of features are recounted and identified. The fall off in the intensity from the centroid of individual features is measured and compared with a standard Gaussian shape. In this process, the features that do not follow a Gaussian like structure are eliminated. This technique is more suitable to use

It can be shown that the wavefront reconstruction accuracy is not a constant, but fluctuates about a mean value irrespective of the sensing geometry at place (see Figs. 28 and 29), by applying Monte Carlo simulations on VMM and FFT methods. These inconsistencies will have a serious effect on the stability of imaging and maintenance of a good Strehl

this parameter is chosen to eliminate features that are not close to a circular shape.

along with the CoG and IWC methods (Vyas et al., 2010a).

**7. Improved consistency: Dither based SHWS**

Fig. 26. The effect of centroiding at different thresholding percentage on spots reconstructed with n ≤ 29.

Fig. 27. Spots reconstructed using n ≤ 29 and 60% thresholding. Spot orientation was 300; with spot size along non-elongated direction being 3 pixels uniformly for all image sizes.

#### **6.3 Gaussian pattern matching**

22 Adaptive Optics

 **60% 65% 70% 75% 80% 85%**

**1 2 3 4 5 6 7 8 9**

Fig. 26. The effect of centroiding at different thresholding percentage on spots reconstructed

**Signal to Noise Ratio**

**0 5 10 15 20**

**Signal to Noise Ratio**

Fig. 27. Spots reconstructed using n ≤ 29 and 60% thresholding. Spot orientation was 300; with spot size along non-elongated direction being 3 pixels uniformly for all image sizes.

**0.5**

**0.2**

**0.4**

**0.6**

**0.8**

**Centroid estimation error**

**1.2**

**1.4**

**1.6**

**1**

**1.5**

**Centroid estimation error**

with n ≤ 29.

**2.5**

**3.5**

**1**

**2**

**3**

At high noise conditions, the thresholded Zernike reconstructed images lead to multiple noisy features along with the actual spot pattern as depicted in Figs. 20 and 21. This is due to the fact that at high noise level conditions there can be large scale features (sometimes the features can be scaled to dimensions comparable to that of the spot) which might not be removed even after the denoising procedure is implemented.

To use this erroneous spot for accurate centroiding, we take the advantage of the fact that in the case of NGS, the spot finally formed at the focal plane of a lens must assume an airy pattern which can be approximated to a Gaussian like structure. The pattern matching algorithm is implementing in three steps: feature recognition, shape identification and profile identification.

#### **6.3.1 Feature recognition**

In this step, the features on the spot pattern image are counted and identified. Feature recognition can be performed by using many existing pattern recognition algorithms. In this algorithm, we used a simple Hough peak identification method to detect the features and number them in the order of peak height. To eliminate small scale features which may arise due to unavoidable scintillations, we impose threshold conditions on the size of the features.

#### **6.3.2 Shape identification**

Most features don't have a circular shape. The circularity or the extent of the feature being circular is measured for each of the features. This can be measured in many ways. In the proposed algorithm, the local centroid for the feature is calculated and the distance from the centroid at which intensity becomes zero is measured. This distance is called the distance parameter. The distance parameter is estimated at different angles (0-360) from the feature centroid position. We define circularity, C of a single feature as the inverse of the variance of the distance parameter computed over different angles from the centroid position. For an ideal circular feature, the circularity is infinity since the variance is zero. A lower cutoff for this parameter is chosen to eliminate features that are not close to a circular shape.

#### **6.3.3 Profile identification**

In the previous step, the shape of the spot was used for selective elimination. In this step, the intensity profile is used to choose the actual spot feature. The number of features are recounted and identified. The fall off in the intensity from the centroid of individual features is measured and compared with a standard Gaussian shape. In this process, the features that do not follow a Gaussian like structure are eliminated. This technique is more suitable to use along with the CoG and IWC methods (Vyas et al., 2010a).

### **7. Improved consistency: Dither based SHWS**

It can be shown that the wavefront reconstruction accuracy is not a constant, but fluctuates about a mean value irrespective of the sensing geometry at place (see Figs. 28 and 29), by applying Monte Carlo simulations on VMM and FFT methods. These inconsistencies will have a serious effect on the stability of imaging and maintenance of a good Strehl

**0 0.1 0.2 0.3 0.4 0.5**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 191

**Realization 1 Realization 2**

**Time (sec)**

Fig. 29. Inconsistency of the wavefront reconstruction accuracy in the case of vector matrix multiply wavefront reconstruction technique when applied to evolving atmosphere like

Probabilistically the point of the best wavefront reconstruction occurs close to the center of the wavefront as shown in Fig. 34 (computed for the case of a single phase screen). Hence it is enough to apply a dither close to the center of the wavefront. The choice of the spatial range over which the dither signal may be applied varies with the number of Shack Hartmann

The dither should be applied at a frequency three times the frequency of adaptive optics wavefront correction. This gives enough time to check the strehl ratio and make suitable corrections to the dither applied. The variance of wavefront reconstruction accuracy as a function of the time interval between the application of two dither signals is shown in Fig. 35. It can be seen that by applying dither within shorter intervals of time gives smaller variance

This dither based sensor can be realized using a Liquid Crystal based Spatial Light Modulator (LC-SLM) by projecting a digital diffractive optical lenslet array on it (Vyas et al., 2009a;

subapertures, the Fried parameter and the degree of freedom in the lenslet array.

in the wavefront reconstruction accuracy and hence is more consistent.

**0.5**

Fig. 30. Schematic demonstration of the dither based sensor

**0.55**

**0.6**

**Wavefront reconstruction accuracy**

turbulence phase screens

**0.65**

**0.7**

Fig. 28. Inconsistency of the wavefront reconstruction accuracy in the case of FFT technique when applied to evolving atmosphere like turbulence phase screens

ratio. In this section, a possible solution to this problem is addressed. The wavefront reconstruction accuracy in SHWS depends strongly on the way wavefront distortion points match with the points of phase estimation (position of the lenslets in the SHWS). A small dither signal which acts like a translating operator on the wavefront sensor with respect to the phase screen can be used to improve the wavefront reconstruction accuracy. For temporally evolving turbulence, the dither signal applied on the sensor can improve the consistency of the wavefront reconstruction accuracy.

For phase screens simulated using the Kolmogorov model, the wavefront reconstruction accuracy is calculated by applying dither such that the sensor shifts in all directions and by different magnitudes. Probabilistically the point of the best wavefront reconstruction is found to be near the center of the phase screen. We show through numerical simulations that the consistency and the accuracy of wavefront reconstruction can be significantly improved using this technique. In real-time systems, the dither signal to be applied can be obtained from the wavefront sensor data of the immediate past within the wavefront decorrelation time. The practicality of building such a sensor is also discussed.

A dither signal is applied on the wavefront sensor (ie. the lenslet array; the detector is placed at a fixed position) so that it is displaced to a new location such that the wavefront reconstruction accuracy is maximized (Fig. 30). This idea is evaluated through numerical simulations. Considering the case of a telescope with an effective diameter of 1m at a site with *r*0 ∼ 10cm and a Shack Hartmann sensor with 100 subapertures Monte Carlo simulations were performed. The wavefront reconstruction accuracy depends on the way in which the wavefront distortion points match with the centers of subapertures. The center of the lenslet array was displaced to different discrete locations with respect to the center of the incoming wavefront within the distance of a single subaperture. The point of the best wavefront reconstruction is the point to which the center of the lenslet array has to be moved to. It can be seen that the application of dither signal not only increases the accuracy of reconstruction, but also significantly improves the consistency in the case of Fourier and VMM reconstruction procedures (Figs. 31 and 32). In a real situation, the position of best wavefront reconstruction must be obtained by looking at the strehl ratio. The position of the maximum intensity point varies with time as shown in Fig. 33.

24 Adaptive Optics

**Realization 1 Realization 2**

**0 0.1 0.2 0.3 0.4 0.5**

**Time (sec)**

Fig. 28. Inconsistency of the wavefront reconstruction accuracy in the case of FFT technique

ratio. In this section, a possible solution to this problem is addressed. The wavefront reconstruction accuracy in SHWS depends strongly on the way wavefront distortion points match with the points of phase estimation (position of the lenslets in the SHWS). A small dither signal which acts like a translating operator on the wavefront sensor with respect to the phase screen can be used to improve the wavefront reconstruction accuracy. For temporally evolving turbulence, the dither signal applied on the sensor can improve the consistency of

For phase screens simulated using the Kolmogorov model, the wavefront reconstruction accuracy is calculated by applying dither such that the sensor shifts in all directions and by different magnitudes. Probabilistically the point of the best wavefront reconstruction is found to be near the center of the phase screen. We show through numerical simulations that the consistency and the accuracy of wavefront reconstruction can be significantly improved using this technique. In real-time systems, the dither signal to be applied can be obtained from the wavefront sensor data of the immediate past within the wavefront decorrelation time. The

A dither signal is applied on the wavefront sensor (ie. the lenslet array; the detector is placed at a fixed position) so that it is displaced to a new location such that the wavefront reconstruction accuracy is maximized (Fig. 30). This idea is evaluated through numerical simulations. Considering the case of a telescope with an effective diameter of 1m at a site with *r*0 ∼ 10cm and a Shack Hartmann sensor with 100 subapertures Monte Carlo simulations were performed. The wavefront reconstruction accuracy depends on the way in which the wavefront distortion points match with the centers of subapertures. The center of the lenslet array was displaced to different discrete locations with respect to the center of the incoming wavefront within the distance of a single subaperture. The point of the best wavefront reconstruction is the point to which the center of the lenslet array has to be moved to. It can be seen that the application of dither signal not only increases the accuracy of reconstruction, but also significantly improves the consistency in the case of Fourier and VMM reconstruction procedures (Figs. 31 and 32). In a real situation, the position of best wavefront reconstruction must be obtained by looking at the strehl ratio. The position of the maximum intensity point

**0.4**

when applied to evolving atmosphere like turbulence phase screens

**0.5**

**0.6**

**Wavefront Reconstruction Accuracy**

the wavefront reconstruction accuracy.

varies with time as shown in Fig. 33.

practicality of building such a sensor is also discussed.

**0.7**

**0.8**

Fig. 29. Inconsistency of the wavefront reconstruction accuracy in the case of vector matrix multiply wavefront reconstruction technique when applied to evolving atmosphere like turbulence phase screens

Fig. 30. Schematic demonstration of the dither based sensor

Probabilistically the point of the best wavefront reconstruction occurs close to the center of the wavefront as shown in Fig. 34 (computed for the case of a single phase screen). Hence it is enough to apply a dither close to the center of the wavefront. The choice of the spatial range over which the dither signal may be applied varies with the number of Shack Hartmann subapertures, the Fried parameter and the degree of freedom in the lenslet array.

The dither should be applied at a frequency three times the frequency of adaptive optics wavefront correction. This gives enough time to check the strehl ratio and make suitable corrections to the dither applied. The variance of wavefront reconstruction accuracy as a function of the time interval between the application of two dither signals is shown in Fig. 35.

It can be seen that by applying dither within shorter intervals of time gives smaller variance in the wavefront reconstruction accuracy and hence is more consistent.

This dither based sensor can be realized using a Liquid Crystal based Spatial Light Modulator (LC-SLM) by projecting a digital diffractive optical lenslet array on it (Vyas et al., 2009a;

Fig. 34. Probability case of VMM

**8. References**

temporally evolving phase screens (100 in number each)

**Variance of**

**wavefront reconstruction accuracy**

**12 x 10−3**

**0 0.2 0.4 0.6 0.8 1 4**

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 193

Fig. 35. Dependence of variance of wavefront reconstruction on the time interval between application of two dither signals on the wavefront sensor. This plot is a result of 20 sets of

Zhao et al., 2006). The application of the dither signal is much simpler here since the physical movement of the sensor is not needed. The speed and physical dimensions of the SLM limits its application to this method in different situations. It has been recently shown that the application of multiple dither sensors improves the wavefront reconstruction accuracy and

Baker, K. L. & Moallem, M. M. (2007). Iteratively weighted centroiding for shack- hartmann

Birch, G. C., Descour, M. R. & Tkaczyk, T. S. (2010). Hyperspectral shack–hartmann test, *Appl.*

Cha, J. W., Ballesta, J. & So, P. T. C. (2010). Shack-hartmann wavefront-sensor-based adaptive

Correia, C., Kulcsár, C., Conan, J.-M. & Raynaud, H.-F. (2008). Hartmann modelling in the

discrete spatial frequency domain application to real-time reconstruction in adaptive

its consistency in the case of large telescope AO systems (Vyas et al., 2011).

URL: *http://www.opticsexpress.org/abstract.cfm?URI=oe-15-8-5147*

wave-front sensors, *Opt. Express* 15(8): 5147–5159.

URL: *http://ao.osa.org/abstract.cfm?URI=ao-49-28-5399*

URL: *http://dx.doi.org/doi/10.1117/1.3475954*

optics system for multiphoton microscopy, 15(4): 046022.

*Opt.* 49(28): 5399–5406.

**Time interval for application of dither**

Fig. 31. Inconsistency of the wavefront reconstruction accuracy in the case of vector matrix multiply wavefront reconstruction technique when applied to evolving atmosphere like turbulence phase screens

Fig. 32. Inconsistency of the wavefront reconstruction accuracy in the case of vector matrix multiply wavefront reconstruction technique when applied to evolving atmosphere like turbulence phase screens

Fig. 33. The propagation of the position coordinate with the maximum wavefront reconstruction accuracy. The analysis performed at time intervals of 0.005 sec

Fig. 34. Probability case of VMM

26 Adaptive Optics

**Before dither based sensor After dither based sensor**

> **Before dither After dither**

**x coordinate y coordinate**

**0 0.1 0.2 0.3 0.4 0.5**

**0 0.1 0.2 0.3 0.4 0.5**

**Time (sec)**

**0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0**

**Time (sec)**

Fig. 33. The propagation of the position coordinate with the maximum wavefront reconstruction accuracy. The analysis performed at time intervals of 0.005 sec

Fig. 32. Inconsistency of the wavefront reconstruction accuracy in the case of vector matrix multiply wavefront reconstruction technique when applied to evolving atmosphere like

**Time (sec)**

Fig. 31. Inconsistency of the wavefront reconstruction accuracy in the case of vector matrix multiply wavefront reconstruction technique when applied to evolving atmosphere like

**0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1**

**0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9**

**Position coordinate**

**Wavefront reconstruction accuracy**

**Wavefront reconstruction accuracy**

turbulence phase screens

turbulence phase screens

Fig. 35. Dependence of variance of wavefront reconstruction on the time interval between application of two dither signals on the wavefront sensor. This plot is a result of 20 sets of temporally evolving phase screens (100 in number each)

Zhao et al., 2006). The application of the dither signal is much simpler here since the physical movement of the sensor is not needed. The speed and physical dimensions of the SLM limits its application to this method in different situations. It has been recently shown that the application of multiple dither sensors improves the wavefront reconstruction accuracy and its consistency in the case of large telescope AO systems (Vyas et al., 2011).

#### **8. References**

Baker, K. L. & Moallem, M. M. (2007). Iteratively weighted centroiding for shack- hartmann wave-front sensors, *Opt. Express* 15(8): 5147–5159.

URL: *http://www.opticsexpress.org/abstract.cfm?URI=oe-15-8-5147*

Birch, G. C., Descour, M. R. & Tkaczyk, T. S. (2010). Hyperspectral shack–hartmann test, *Appl. Opt.* 49(28): 5399–5406.

URL: *http://ao.osa.org/abstract.cfm?URI=ao-49-28-5399*

López-Quesada, C., Andilla, J. & Martín-Badosa, E. (2009). Correction of aberration

Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor 195

Poyneer, L. A. (2003). Scene-based shack-hartmann wave-front sensing: Analysis and

Primmerman, C. A., Murphy, D. V., Page, D. A., Zollars, B. G. & Barclay, H. T. (1991).

Roopashree, M. B., Vyas, A. & Prasad, B. R. (2010). Multilayered temporally evolving phase

Roopashree, M. B., Vyas, A. & Prasad, B. R. (2011). Experimental evaluation of centroiding

Schreiber, L., Foppiani, I., Robert, C., Diolaiti, E., Conan, J. M. & Lombini, M. (2009). Laser

Shack, R. V. & Platt, B. C. (1971). Production and use of a lenticular hartmann screen, *J. Opt.*

Silbaugh, E. E., Welsh, B. M. & Roggemann, M. C. (1996). Characterization of atmospheric

Southwell, W. H. (1980). Wave-front estimation from wave-front slope measurements, *J. Opt.*

Thomas, S. J., Adkins, S., Gavel, D., Fusco, T. & Michau, V. (2008). Study of optimal wavefront

Thompson, L. A. & Castle, R. M. (1992). Experimental demonstration of a rayleigh-scattered

Vyas, A., Roopashree, M. B. & Prasad, B. R. (2009a). Digital long focal length lenslet array

Vyas, A., Roopashree, M. B. & Prasad, B. R. (2009b). Performance of centroiding algorithms

screens based on statistical interpolation, Vol. 7736, SPIE, p. 77363Z.

algorithms at different light intensity and noise levels, AIP.

URL: *http://ao.osa.org/abstract.cfm?URI=ao-50-20-3519*

URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-13-12-2453*

URL: *http://dx.doi.org/10.1111/j.1365-2966.2008.13110.x*

pp. 366–369. DOI:10.1109/ARTCom.2009.30.

URL: *http://www.opticsinfobase.org/abstract.cfm?URI=josa-70-8-998*

URL: *http://ao.osa.org/abstract.cfm?URI=ao-48-6-1084*

URL: *http://ao.osa.org/abstract.cfm?URI=ao-42-29-5807*

URL: *http://link.aip.org/link/?PSI/7736/77363Z/1*

simulation, *Appl. Opt.* 42(29): 5807–5815.

48(6): 1084–1090.

353: 141–143.

50(20): 3519–3528.

13(12): 2453–2460.

*Soc. Am.* 20(61): 656–660.

*Soc. Am.* 70(8): 998–1006.

*Society* 387(1): 173–187.

17: 1485–1487.

*photonics*.

in holographic optical tweezers using a shack-hartmann sensor, *Appl. Opt.*

Compensation of atmospheric optical distortion using a synthetic beacon, *Nature*

guide stars for extremely large telescopes: efficient shack-hartmann wavefront sensor design using weighted center-of-gravity algorithm, *Mon. Not. R. Astron. Soc.* . Sedmak, G. (2004). Implementation of fast-fourier-transform-based simulations of extra-large atmospheric phase and scintillation screens, *Appl. Opt.* 43: 4527–4538. Sergeyev, A. V. & Roggemann, M. C. (2011). Monitoring the statistics of turbulence:

Fried parameter estimation from the wavefront sensor measurements, *Appl. Opt.*

turbulence phase statistics using wave-front slope measurements, *J. Opt. Soc. Am. A*

sensing with elongated laser guide stars, *Monthly Notices of the Royal Astronomical*

laser guide star at 351 nm, *Monthly Notices of the Royal Astronomical Society*

using spatial light modulator, *Proceedings of the international conference on optics and*

at low light level conditions in adaptive optics, *Advances in Recent Technologies in Communication and Computing, 2009. ARTCom '09. International Conference on*,

optics, 7015(1): 701551.

URL: *http://dx.doi.org/doi/10.1117/12.788455*

Davis, D. S., Hickson, P., Herriot, G. & She, C.-Y. (2006). Temporal variability of the telluric sodium layer, *Opt. Lett.* 31(22): 3369–3371.

URL: *http://ol.osa.org/abstract.cfm?URI=ol-31-22-3369*

URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-3-11-1852*

	- URL: *http://ao.osa.org/abstract.cfm?URI=ao-38-11-2161*

URL: *http://www.opticsinfobase.org/abstract.cfm?URI=josa-67-3-375*

López-Quesada, C., Andilla, J. & Martín-Badosa, E. (2009). Correction of aberration in holographic optical tweezers using a shack-hartmann sensor, *Appl. Opt.* 48(6): 1084–1090.

URL: *http://ao.osa.org/abstract.cfm?URI=ao-48-6-1084*

28 Adaptive Optics

Davis, D. S., Hickson, P., Herriot, G. & She, C.-Y. (2006). Temporal variability of the telluric

Drummond, J., Telle, J., Denman, C., Hillman, P., Spinhirne, J. & Christou, J. (2004).

Freischlad, K. R. & Koliopoulos, C. L. (1986). Modal estimation of a wave front from

Fried, D. L. (1977). Least-square fitting a wave-front distortion estimate to an array of

Fried, D. L. (1982). Anisoplanatism in adaptive optics, *Monthly Notices of the Royal Astronomical*

Fugate, R. Q., Fried, D. L., Ameer, G. A., Boeke, B. R., Browne, S. L., Roberts, P. H., Ruane,

distortion using scattered light from a laser guide-star, *Nature* 353: 144–146. Fusco, T., Nicolle, M., Rousset, G., Michau, V., Beuzit, J.-L. & Mouillet, D. (2004). Optimization of a shack-hartmann-based wavefront sensor for xao systems, 5490(1): 1155–1166.

Gilles, L. & Ellerbroek, B. (2006). Shack-hartmann wavefront sensing with elongated sodium laser beacons: centroiding versus matched filtering, *Appl. Opt.* 45: 6568–6576. Harding, C. M., Johnston, R. A. & Lane, R. G. (1999). Fast simulation of a kolmogorov phase

Hartmann, J. (1900). Bemerkungen uber den bau und die justirung von spektrographen, *Z.*

Hudgin, R. H. (1977). Wave-front reconstruction for compensated imaging, *J. Opt. Soc. Am.*

Johansson, E. M. & Gavel, D. T. (1994). Simulation of stellar speckle imaging, 2200(1): 372–383.

Kellerer, A. N. & Kellerer, A. M. (2011). Error propagation: a comparison of shack–hartmann

Leroux, C. & Dainty, C. (2010). Estimation of centroid positions with a matched-filter algorithm: relevance for aberrometry of the eye, *Opt. Express* 18(2): 1197–1206.

URL: *http://www.opticsinfobase.org/abstract.cfm?URI=josa-67-3-375*

URL: *http://www.opticsexpress.org/abstract.cfm?URI=oe-18-2-1197*

and curvature sensors, *J. Opt. Soc. Am. A* 28(5): 801–807. URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-28-5-801*

Hardy, J. W. (1998). *Adaptive Optics for Astronomical Telescopes*, Oxford University Press. Hart, M. L., Jacobsen, A., Angel, J. R. P., Wittman, D., Dekany, R., Mccarthy, D. & Beletic,

R. E., Tyler, G. A. & Wopat, L. M. (1991). Measurement of atmospheric wavefront

J. (1995). Adaptive optics experiments using sodium laser guide stars, *Astrophys J*

Photometry of a sodium laser guide star from the starfire optical range. ii. compensating the pump beam, *Publications of the Astronomical Society of the Pacific*

difference measurements using the discrete fourier transform, *J. Opt. Soc. Am. A*

optics, 7015(1): 701551.

116: 952–964.

3(11): 1852–1861.

*Society* 72: 52–61.

439: 455–473.

67(3): 375–378.

*Instrumentenkd.* 20: 47–58.

URL: *http://dx.doi.org/doi/10.1117/12.788455*

sodium layer, *Opt. Lett.* 31(22): 3369–3371.

URL: *http://dx.doi.org/doi/10.1117/12.549562*

URL: *http://dx.doi.org/doi/10.1117/12.177254*

URL: *http://ao.osa.org/abstract.cfm?URI=ao-38-11-2161*

screen, *Appl. Opt.* 38(11): 2161–2170.

URL: *http://ol.osa.org/abstract.cfm?URI=ol-31-22-3369*

URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-3-11-1852*

phase-difference measurements, *J. Opt. Soc. Am.* 67(3): 370–375. URL: *http://www.opticsinfobase.org/abstract.cfm?URI=josa-67-3-370* Poyneer, L. A. (2003). Scene-based shack-hartmann wave-front sensing: Analysis and simulation, *Appl. Opt.* 42(29): 5807–5815.

URL: *http://ao.osa.org/abstract.cfm?URI=ao-42-29-5807*

URL: *http://ao.osa.org/abstract.cfm?URI=ao-50-20-3519*

URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-13-12-2453*

Southwell, W. H. (1980). Wave-front estimation from wave-front slope measurements, *J. Opt. Soc. Am.* 70(8): 998–1006.

URL: *http://www.opticsinfobase.org/abstract.cfm?URI=josa-70-8-998*

Thomas, S. J., Adkins, S., Gavel, D., Fusco, T. & Michau, V. (2008). Study of optimal wavefront sensing with elongated laser guide stars, *Monthly Notices of the Royal Astronomical Society* 387(1): 173–187.

URL: *http://dx.doi.org/10.1111/j.1365-2966.2008.13110.x*

**9**

 *China*

**Measurement Error of**

**Shack-Hartmann Wavefront Sensor**

Chaohong Li, Hao Xian, Wenhan Jiang and Changhui Rao *1Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu,*

*2School of Ophthalmology and Optometry, Wenzhou Medical School,*

A Shack-Hartmann sensor is one of the most important and popular wavefront sensors used in an adaptive optics system to measure the aberrations caused by either atmospheric turbulence, laser transmission, or the living eye [1-7]. Its design was based on an aperture array that was developed in 1900 by Johannes Franz Hartmann as a means to trace individual rays of light through the optical system of a large telescope, thereby testing the quality of the image.[8] In the late 1960s Roland Shack and Platt modified the Hartmann screen by replacing the apertures in an opaque screen by an array of lenslets [9-10]. The terminology as proposed by Shack and Platt was "Hartmann-screen". The fundamental principle seems to be documented

The schematic of a Shack-Hartman wavefront sensor is shown in Figure 1. It consists of an array of lenses (called lenslets, see Figure 1) of the same focal length. Each is focused onto a photon sensor (typically a CCD array or quad-cell). The local tilt of the wavefront across each lens can then be calculated from the position of the focal spot on the sensor. Any phase aberration can be approximated to a set of discrete tilts. By sampling an array of lenslets, all of these tilts can be measured and the whole wavefront can be approximated. Since only tilts are measured, the Shack-Hartmann can not measure the discontinuous steps of wavefront.

even before Huygens by the Jesuit philosopher, Christopher Scheiner [11].

Fig. 1. Schematic of a Shack-Hartmann wavefront sensor

**1. Introduction**

URL: *http://ao.osa.org/abstract.cfm?URI=ao-45-1-90*

Zou, W., Thompson, K. P. & Rolland, J. P. (2008). Differential shack-hartmann curvature sensor: local principal curvature measurements, *J. Opt. Soc. Am. A* 25(9): 2331–2337. URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-25-9-2331*

## **Measurement Error of Shack-Hartmann Wavefront Sensor**

Chaohong Li, Hao Xian, Wenhan Jiang and Changhui Rao *1Institute of Optics and Electronics, Chinese Academy of Sciences, Chengdu, 2School of Ophthalmology and Optometry, Wenzhou Medical School, China*

## **1. Introduction**

30 Adaptive Optics

196 Topics in Adaptive Optics

Vyas, A., Roopashree, M. B. & Prasad, B. R. (2010a). Centroid detection by gaussian pattern

Vyas, A., Roopashree, M. B. & Prasad, B. R. (2010b). Improved iteratively weighted

Vyas, A., Roopashree, M. B. & Prasad, B. R. (2010c). Noise reduction in the centroiding of laser guide star spot pattern using thresholded zernike reconstructor, 7736(1): 77364E.

Vyas, A., Roopashree, M. & Prasad, B. R. (2011). Multi-dither shack hartmann sensor for large

Wei, X. & Thibos, L. (2010). Design and validation of a scanning shack hartmann aberrometer

Wu, D., Wu, S.-Z., Niu, L.-G., Chen, Q.-D., Wang, R., Song, J.-F., Fang, H.-H. & Sun, H.-B. (2010). High numerical aperture microlens arrays of close packing, 97(3): 031109.

Wu, H.-L., Yan, H.-X., Li, X.-Y. & Li, S.-S. (2009). Statistical interpolation method of turbulent

Yin, X., Li, X., Zhao, L. & Fang, Z. (2009). Adaptive thresholding and dynamic windowing

Zhao, L., Bai, N., Li, X., Ong, L. S., Fang, Z. P. & Asundi, A. K. (2006). Efficient implementation

Zou, W., Thompson, K. P. & Rolland, J. P. (2008). Differential shack-hartmann curvature sensor: local principal curvature measurements, *J. Opt. Soc. Am. A* 25(9): 2331–2337.

method for automatic centroid detection of digital shack-hartmann wavefront sensor,

of a spatial light modulator as a diffractive optical microlens array in a digital

URL: *http://www.opticsinfobase.org/abstract.cfm?URI=AO-2011-ATuA4*

URL: *http://www.opticsexpress.org/abstract.cfm?URI=oe-18-2-1134*

URL: *http://www.opticsexpress.org/abstract.cfm?URI=oe-17-17-14649*

Published By Foundation of Computer Science.

URL: *http://dx.doi.org/doi/10.1117/12.856640*

URL: *http://dx.doi.org/doi/10.1063/1.3464979*

*Appl. Opt.* 48(32): 6088–6098.

phase screen, *Opt. Express* 17(17): 14649–14664.

URL: *http://ao.osa.org/abstract.cfm?URI=ao-48-32-6088*

URL: *http://ao.osa.org/abstract.cfm?URI=ao-45-1-90*

shack-hartmann wavefront sensor, *Appl. Opt.* 45(1): 90–94.

URL: *http://josaa.osa.org/abstract.cfm?URI=josaa-25-9-2331*

*and Applications*, Optical Society of America, p. ATuA4.

pp. 758806–1–11.

matching in adaptive optics, *International Journal of Computer Applications* 1(25): 30–35.

centroiding for accurate spot detection in laser guide star based shack hartmann sensor, *Atmospheric and Oceanic propagation of electromagnetic waves IV*, SPIE,

telescopes: A numerical performance evaluation, *Adaptive Optics: Methods, Analysis*

for measurements of the eye over a wide field of view, *Opt. Express* 18(2): 1134–1143.

A Shack-Hartmann sensor is one of the most important and popular wavefront sensors used in an adaptive optics system to measure the aberrations caused by either atmospheric turbulence, laser transmission, or the living eye [1-7]. Its design was based on an aperture array that was developed in 1900 by Johannes Franz Hartmann as a means to trace individual rays of light through the optical system of a large telescope, thereby testing the quality of the image.[8] In the late 1960s Roland Shack and Platt modified the Hartmann screen by replacing the apertures in an opaque screen by an array of lenslets [9-10]. The terminology as proposed by Shack and Platt was "Hartmann-screen". The fundamental principle seems to be documented even before Huygens by the Jesuit philosopher, Christopher Scheiner [11].

The schematic of a Shack-Hartman wavefront sensor is shown in Figure 1. It consists of an array of lenses (called lenslets, see Figure 1) of the same focal length. Each is focused onto a photon sensor (typically a CCD array or quad-cell). The local tilt of the wavefront across each lens can then be calculated from the position of the focal spot on the sensor. Any phase aberration can be approximated to a set of discrete tilts. By sampling an array of lenslets, all of these tilts can be measured and the whole wavefront can be approximated. Since only tilts are measured, the Shack-Hartmann can not measure the discontinuous steps of wavefront.

Fig. 1. Schematic of a Shack-Hartmann wavefront sensor

Measurement Error of Shack-Hartmann Wavefront Sensor 199

1 2 1 2 1 1 1 1

*nm nm nm nm*

(3)

*x I y I*

*I I*

2

(4)

(5)

*cr* is the variance of the centroid fluctuation in one

(6)

*cr* and 2

(7)

*cs* are both

2 *s cs s G V*

2 2 2 1 1 2 2

*N L L L*

*cs* is the variance of the centroid fluctuation in one direction ( *X* or *Y*), induced by the

defined in pixel2 units. *Gs* is the equivalent Gauss width of the signal spot and is defined in

where η is the positive constant. η is 0.353 when the diffraction aperture is square and 0.431

Based on Eq. (3), the centroid position in the *X* direction can be expressed by Eq. (7). The

1 1 *c cs cb*

where *c x* is the calculated centroid position of the signal in the *X* direction, *cs x* is the real centroid position of the signal in the *X* direction, and *cb x* is the centroid position induced by the total noise, except for the signal in the *X* direction, where the total noise largely includes the readout noise of the detector and the heterogeneous light noise. *sbr S B n m*, , / *n m* , *<Snm>* is the collective average of the signal intensity in the (n,m)th pixel, and *<Bnm>* is the collective average of the total noise intensity in the (n,m)th

*l D*

1

1 12

*r cr r*

direction ( *X* or *Y*), induced by the readout noise of the detector. 2

detailed derivation process of this expression is shown in appendix 1.1.

*V*

*s f G*

*sbr xx x sbr sbr*

*nm nm*

,

*m n m n*

1 1 1 1

The formulas which evaluate the centroid error associated with the signal's photon shot noise and the readout noise of the detector, respectively, have been derived by Cao *et al.* When the photon shot noise of the signal is considered alone, the centroid fluctuation error was obtained by introducing the Gauss width of the signal. When the readout noise of the detector is considered alone, the centroid fluctuation error was also obtained, and its results are shown in Eq. (4) and Eq. (5), respectively [14]. *Nr* is the rms error induced by the fluctuation of the readout noise in each pixel (with units of photoelectron counts). *Vr* is the sum of the readout noise's photoelectron count among all of the pixels within the

1 2 1 2

*L L L L*

*m n m n i i L L L L*

*x y*

corresponding subaperture.

pixel units by the expression:

pixel (with units of ADU).

photon shot noise of signal itself, and 2

when the diffraction aperture is circular.

2

Tyler and Fried have obtained the theory expression, which evaluates the angular position error when a quadrant detector is used in the SHWFS [12]. The formula they obtained, based on circular aperture diffraction, is shown in Eq. (1)

$$
\sigma = \frac{3\pi}{16} \frac{1}{SNR} \frac{\lambda}{D} \tag{1}
$$

where *SNR* is defined as the ratio of the signal's photoelectron counts to the noise's fluctuation intensity within the detection area, *λ* is the wavelength, and *D* is the diameter of the aperture. Their analysis did not discuss the size of the incoming light spot on the detection area in detail. The formula was obtained based on a quadrant detector alone. Generally, the theory expression obtained by Tyler and Fried is not suitable for describing the angular position error when the scale of the discrete detector arrays is greater than 2×2 pixels.

Hardy has described formulas that can be used to evaluate the angular position error [13], under the conditions that the photon shot noise of signal is dominant. Although his formulas discussed the size of the diffraction-limited spot on the discrete detector arrays, it is reliable only under the approximation condition that *l/f>>λ/D* or *l/f<<λ/D* is satisfied, where *l* is the length of a pixel, *f* is the focal length, and *Vs* is the count of signal photoelectrons. Eq. (2) shows the formulas based on square aperture diffraction.

$$\sigma = \begin{cases} 0.277 \frac{\lambda}{D} / \, V\_s^{1/2} & \text{(when } 1 / f << \lambda / D\text{)}\\ 0.500 \frac{\lambda}{D} / \, V\_s^{1/2} & \text{(when } 1 / f >> \lambda / D\text{)} \end{cases} \tag{2}$$

Cao *et al*. have also analyzed the measurement error of a SHWFS. Their work emphasized the discrete sampling error of a CCD and obtained a formula which is used to describe the centroid position error induced by the readout noise of the CCD and the photon shot noise of the signal [14]. Their research results are only the approximation of some cases discussed in this article. Jiang *et al*. have partially analyzed the measurement error of SHWFS and performed their method by setting a fixed threshold to suppress the impact of the random noise [15].

In this chapter, the wavefront error of a Shack-Hartmann wavefront sensor was analyzed in detail based on the research results of angular position error and wavefront error [16-17]. and the formula used to evaluate the wavefront error was derived, it concerns with the signal to noise ratio, number of photons and reconstruction matrix also.

#### **2. The angular position error caused by random noise**

The wavefront to be measured is segmented into many subwavefronts by lenslet arrays, and the light spots at the focal plane of the subapertures are detected by the CCD. Particularly, the analysis is based on the notion that the wavefront is essentially flat over each subaperture and *r0>D* (*r0* is the coherent length of incoming wavefront). The centroid position can be calculated by Eq. (3) [18]. The detection area of the subaperture is *L1*×*L2* pixels, and *xnm* and *ynm* are the (n,m)th pixel's *X* coordinate and *Y* coordinates, respectively. *Inm* is the total intensity in the (n,m)th pixel, including the signal photons and all other noises.

Tyler and Fried have obtained the theory expression, which evaluates the angular position error when a quadrant detector is used in the SHWFS [12]. The formula they obtained, based

> 3 1 16 *SNR D*

where *SNR* is defined as the ratio of the signal's photoelectron counts to the noise's fluctuation intensity within the detection area, *λ* is the wavelength, and *D* is the diameter of the aperture. Their analysis did not discuss the size of the incoming light spot on the detection area in detail. The formula was obtained based on a quadrant detector alone. Generally, the theory expression obtained by Tyler and Fried is not suitable for describing the angular position error when the scale of the discrete detector arrays is greater than 2×2

Hardy has described formulas that can be used to evaluate the angular position error [13], under the conditions that the photon shot noise of signal is dominant. Although his formulas discussed the size of the diffraction-limited spot on the discrete detector arrays, it is reliable only under the approximation condition that *l/f>>λ/D* or *l/f<<λ/D* is satisfied, where *l* is the length of a pixel, *f* is the focal length, and *Vs* is the count of signal

photoelectrons. Eq. (2) shows the formulas based on square aperture diffraction.

signal to noise ratio, number of photons and reconstruction matrix also.

**2. The angular position error caused by random noise**

1/2

*s*

*s*

0.277 / ( / / )

*V when l f D D*

(2)

*V when l f D D*

0.500 / ( / / )

1/2

Cao *et al*. have also analyzed the measurement error of a SHWFS. Their work emphasized the discrete sampling error of a CCD and obtained a formula which is used to describe the centroid position error induced by the readout noise of the CCD and the photon shot noise of the signal [14]. Their research results are only the approximation of some cases discussed in this article. Jiang *et al*. have partially analyzed the measurement error of SHWFS and performed their method by setting a fixed threshold to suppress the impact of the random

In this chapter, the wavefront error of a Shack-Hartmann wavefront sensor was analyzed in detail based on the research results of angular position error and wavefront error [16-17]. and the formula used to evaluate the wavefront error was derived, it concerns with the

The wavefront to be measured is segmented into many subwavefronts by lenslet arrays, and the light spots at the focal plane of the subapertures are detected by the CCD. Particularly, the analysis is based on the notion that the wavefront is essentially flat over each subaperture and *r0>D* (*r0* is the coherent length of incoming wavefront). The centroid position can be calculated by Eq. (3) [18]. The detection area of the subaperture is *L1*×*L2* pixels, and *xnm* and *ynm* are the (n,m)th pixel's *X* coordinate and *Y* coordinates, respectively. *Inm* is the total intensity in the (n,m)th pixel, including the signal photons and all other noises.

(1)

on circular aperture diffraction, is shown in Eq. (1)

pixels.

noise [15].

$$\mathbf{x}\_{i} = \frac{\sum\_{m=1}^{L\_{1}} \sum\_{n=1}^{L\_{2}} \mathbf{x}\_{nm} I\_{nm}}{\sum\_{m=1}^{L\_{1}} \sum\_{n=1}^{L\_{2}} I\_{nm}}, \mathbf{y}\_{i} = \frac{\sum\_{m=1}^{L\_{1}} \sum\_{n=1}^{L\_{2}} \mathbf{y}\_{nm} I\_{nm}}{\sum\_{m=1}^{L\_{1}} \sum\_{n=1}^{L\_{2}} I\_{nm}} \tag{3}$$

The formulas which evaluate the centroid error associated with the signal's photon shot noise and the readout noise of the detector, respectively, have been derived by Cao *et al.* When the photon shot noise of the signal is considered alone, the centroid fluctuation error was obtained by introducing the Gauss width of the signal. When the readout noise of the detector is considered alone, the centroid fluctuation error was also obtained, and its results are shown in Eq. (4) and Eq. (5), respectively [14]. *Nr* is the rms error induced by the fluctuation of the readout noise in each pixel (with units of photoelectron counts). *Vr* is the sum of the readout noise's photoelectron count among all of the pixels within the corresponding subaperture.

$$\left| \left. \phi\_{\rm cs} \right\|^{2} = \frac{\left. \mathcal{G}\_{s} \right|^{2}}{V\_{s}} \tag{4}$$

$$\left|\phi\_{cr}\right|^2 = \frac{N\_r^{-2}}{{V\_r}^2} L\_1 L\_2 \frac{L\_1^2 - 1}{12} \tag{5}$$

2 *cs* is the variance of the centroid fluctuation in one direction ( *X* or *Y*), induced by the photon shot noise of signal itself, and 2 *cr* is the variance of the centroid fluctuation in one direction ( *X* or *Y*), induced by the readout noise of the detector. 2 *cr* and 2 *cs* are both defined in pixel2 units. *Gs* is the equivalent Gauss width of the signal spot and is defined in pixel units by the expression:

$$\mathbf{G}\_s = \eta \frac{f}{l} \frac{\lambda}{D} \tag{6}$$

where η is the positive constant. η is 0.353 when the diffraction aperture is square and 0.431 when the diffraction aperture is circular.

Based on Eq. (3), the centroid position in the *X* direction can be expressed by Eq. (7). The detailed derivation process of this expression is shown in appendix 1.1.

$$\mathbf{x}\_c = \frac{sbr}{1+sbr}\mathbf{x}\_{cs} + \frac{1}{1+sbr}\mathbf{x}\_{cb} \tag{7}$$

where *c x* is the calculated centroid position of the signal in the *X* direction, *cs x* is the real centroid position of the signal in the *X* direction, and *cb x* is the centroid position induced by the total noise, except for the signal in the *X* direction, where the total noise largely includes the readout noise of the detector and the heterogeneous light noise. *sbr S B n m*, , / *n m* , *<Snm>* is the collective average of the signal intensity in the (n,m)th pixel, and *<Bnm>* is the collective average of the total noise intensity in the (n,m)th pixel (with units of ADU).

Measurement Error of Shack-Hartmann Wavefront Sensor 201

The intensity distribution of the signal's light spot at the focal plane of the subaperture can be calculated by circular or square aperture diffraction approximations. On the other hand, the Gauss distribution can also be used to approximately describe the intensity distribution

different approximation conditions. The detailed derivation process is shown in appendix

*Dl f* can be calculated by Eq. (11) .

1 /(1 ( ) ( ))

1 /[1 exp( )] 2

2 2 01 11

*Jr Jr*

1 0

*k*

22 2 2 2

*f f f*

When the direct current part of the noise (including the heterogeneous light noise and the readout noise of the CCD) is subtracted, it can be considered as white noise, and 0 *Bnm* . Then, the standard deviation of the angular position error in the *X* direction caused by noise

> ( , ,, ) [ ( 1) /12] 1 *Dl f LL L l C snr f*

2 1/2 2 2 1/2 2 2 1/2 ( , ,, ) 2

/ 1 [( 1) /(12 )] ( ) /

*L N l f L L*

2 1/2

where *SNR* has the same definition as in Eq. (1). 2*L Nb* is the sum of the rms error of total noise among all of the pixels within the detection area, and expresses the total intensity of

*b b s*

/( ) [ ( 1) /12] max( ) /( )

*N LN V l L L*

*b bs D l f nm b s*

*LN V D SNR D*

1/2 expresses the photon shot noise induced only by the incoming signal.

*C S LN V f*

2 1/2

, *n* is a positive constant,

*nlD f*

sin ( 1) (2 ) [ /( )] 2 (2 1) (2 1)!

2 2 1 1 1 2

*k k*

(12)

( )

*circular aperture diffraction approximation x x x kk*

( )

*square aperture diffraction approximation*

*tribution approximation*)

*Dl f* are described in Eq. (11) with

(11)

(13)

of the light spot. The analytic expressions of *C*( , ,, )

( , ,, )

*C*

can be described by Eq. (12):

and θ is the diffraction angle.

*Dl f*

(

*x r x*

2

When 1 2 *LLL* , the Eq. (12) can also be expressed by Eq. (13)

2 2 1/2

1

noise fluctuation. *Vs*

*SNR D*

where 1 1 1

 *Gauss dis*

/ 2 4 4 sin( ) sin( ) , 2 2 *D D l Dl Dl*

12 1

2.1~2.3. The value of *C*( , ,, )

Based on the error transition principles, the rms error of centroid measurement induced by random noise in the *X* direction can be written in Eq. (8) as:

$$\phi\_c = \text{I}(\frac{sbr}{1+sbr})^2 \phi\_{cs}{}^2 + (\frac{1}{1+sbr})^2 \phi\_{cb}{}^2 \text{I}^{1/2} \tag{8}$$

where *cs* is the rms error of centroid measurement in the *X* direction induced by the signal's photon shot noise and*cb* is the rms error of centroid measurement in the *X* direction induced by all the other noise. The signal is mostly comprised of heterogeneous light and readout noise.

If there were no heterogeneous light and readout noise in the detection area, the signal's photon shot noise should be the unique noise resource which affects centroid measurement. Based on Eq. (4) and Eq. (6), when the discrete sampling error of the detector is ignored, the rms error of angular position in the *X* direction caused by the photon shot noise of the signal can be written as:

$$
\sigma\_1 = (\frac{G\_s}{V\_s})^{1/2} \frac{l}{f} = \eta \frac{\lambda}{D} V\_s^{-1/2} \tag{9}
$$

When the photon shot noise of the signal is small compared with the readout noise and the heterogeneous light noise, the heterogeneous light noise and the readout noise become the primary noise, which affects centroid calculation. When the heterogeneous light noise can be considered as a uniform noise, like the readout noise of the detector, it exists in each pixel and it has the same fluctuation characteristics among the pixels in the detection area. So, the noise in one pixel (including the heterogeneous light noise and the readout noise of the CCD) can be summed and described by *Nb* . *Nb* is defined as the rms error of the heterogeneous light noise and the readout noise photoelectron count in one pixel, and it has the same fluctuation characteristics as the readout noise of the detector. *Nb* has units of ADU. Subsequently, the rms error of centroid measurement in the *X* direction caused by the heterogeneous light noise and the readout noise of the CCD can be written as:

$$\begin{aligned} \phi\_c &= \frac{1}{1 + sbr} \phi\_{cb} \\ &= \left(1 + \sum S\_{nm} \,/\sum B\_{nm}\right)^{-1} \left[L\_1 L\_2 \left(L\_1 \, ^2 - 1\right) / 12\right]^{1/2} N\_b / V\_b \\ &= \left(\sum B\_{nm} \,/\, N\_b + \left(\Box\_{\{\lambda, D, l, f\}} \cdot snr\right)^{-1} \left[L\_1 L\_2 \left(L\_1 \, ^2 - 1\right) / 12\right]^{1/2}\right) \end{aligned} \tag{10}$$

where max( ) / *nm b snr S N* , *V B b nm* , and max( ) *Snm* is the signal's peak intensity. *snr* is defined as the ratio of the signal's peak intensity to the rms error induced by the background noise, *Bnm* is the average intensity of noise in the *(n, m)*th pixel, which includes the heterogeneous light noise and the readout noise of the CCD. *C*( , ,, ) *Dl f* is the light spot constant which is defined as the ratio of the total signal intensity to the signal's peak intensity in the subaperture, and its value can be measured or calculated exactly by ( , ,, ) *D l* /max( ) *C SS f nm nm* .

Based on the error transition principles, the rms error of centroid measurement induced by

2 2 1 2 2 1/2 [( ) ( ) ] 1 1 *c cs cb*

*cs* is the rms error of centroid measurement in the *X* direction induced by the

1/2 1/2

  1 2 1/2

1 2 1/2

(10)

*Dl f* is the light spot

12 1

*nm nm b b*

*S B LL L N V*

(1 / ) [ ( 1) /12] /

( / ) [ ( 1) /12]

where max( ) / *nm b snr S N* , *V B b nm* , and max( ) *Snm* is the signal's peak intensity. *snr* is defined as the ratio of the signal's peak intensity to the rms error induced by the background noise, *Bnm* is the average intensity of noise in the *(n, m)*th pixel, which includes

constant which is defined as the ratio of the total signal intensity to the signal's peak intensity in the subaperture, and its value can be measured or calculated exactly

*B N C snr L L L*

( , ,, ) 12 1

*cb* is the rms error of centroid measurement in the *X*

(9)

(8)

*sbr sbr*

direction induced by all the other noise. The signal is mostly comprised of heterogeneous

If there were no heterogeneous light and readout noise in the detection area, the signal's photon shot noise should be the unique noise resource which affects centroid measurement. Based on Eq. (4) and Eq. (6), when the discrete sampling error of the detector is ignored, the rms error of angular position in the *X* direction caused by the photon shot noise of the signal

2

*s*

heterogeneous light noise and the readout noise of the CCD can be written as:

*nm b D l f*

the heterogeneous light noise and the readout noise of the CCD. *C*( , ,, )

1 1

*c cb*

by ( , ,, ) *D l* /max( ) *C SS*

*f nm nm* .

*sbr*

1 ( ) *s s*

*G l V V fD*

When the photon shot noise of the signal is small compared with the readout noise and the heterogeneous light noise, the heterogeneous light noise and the readout noise become the primary noise, which affects centroid calculation. When the heterogeneous light noise can be considered as a uniform noise, like the readout noise of the detector, it exists in each pixel and it has the same fluctuation characteristics among the pixels in the detection area. So, the noise in one pixel (including the heterogeneous light noise and the readout noise of the CCD) can be summed and described by *Nb* . *Nb* is defined as the rms error of the heterogeneous light noise and the readout noise photoelectron count in one pixel, and it has the same fluctuation characteristics as the readout noise of the detector. *Nb* has units of ADU. Subsequently, the rms error of centroid measurement in the *X* direction caused by the

random noise in the *X* direction can be written in Eq. (8) as:

where

can be written as:

light and readout noise.

signal's photon shot noise and

*sbr*

The intensity distribution of the signal's light spot at the focal plane of the subaperture can be calculated by circular or square aperture diffraction approximations. On the other hand, the Gauss distribution can also be used to approximately describe the intensity distribution of the light spot. The analytic expressions of *C*( , ,, ) *Dl f* are described in Eq. (11) with different approximation conditions. The detailed derivation process is shown in appendix 2.1~2.3. The value of *C*( , ,, ) *Dl f* can be calculated by Eq. (11) .

$$\mathbf{C}\_{\left(\lambda,D,l,f\right)} = \begin{cases} 1 / \left(1 - f\_0^{\;2} \left(r\_1\right) - f\_1^{\;1} \left(r\_1\right)\right) \\ \text{(circular\;\;aperture\;\;diffraction\;\;approximation)} \\ \left[\frac{\pi}{2} / \left(-\frac{\sin^2 x\_1}{x\_1} + \sum\_{k=0}^{+\infty} \frac{(-1)^k \times \left(2x\_1\right)^{2k+1}}{(2k+1)\times \left(2k+1\right)!}\right) \right] \\ \text{(square\;\;aperture\;\;diffraction\;\;approximation)} \\ 1 / \left[1 - \exp\left(-\frac{n^2 l^2 D^2}{2\pi f^2 \lambda^2}\right)\right] \end{cases} \tag{11}$$

( *Gauss dis tribution approximation*)

where 1 1 1 / 2 4 4 sin( ) sin( ) , 2 2 *D D l Dl Dl x r x f f f* , *n* is a positive constant, and θis the diffraction angle.