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

Wavefront sensor is a device that measures the optical wavefront aberration. The Shack Hartmann wavefront sensor (SHWS), named after Johannes Franz Hartmann and Roland Shack, is one of the most often used optical wavefront sensor. It is made up of an array of microlenses (all having the same focal length and aperture size) and a detector placed at the focal plane of these microlenses. Johannes Franz Hartmann developed a device that consisted of an opaque screen with multiple holes to test the quality of imaging in large telescope systems and fine-tune the telescope focus (Hartmann, 1900). In the process of developing an adaptive optics system to improve the resolution of satellite images taken from the earth, Roland Shack came up with a feasible model of the sensor by using an array of tiny lenses instead of holes (Shack & Platt, 1971). The technological advancements in the field of optical fabrication allowed the industry to make lenses as small as 100 μm using materials like fused silica, ZnS, ZnSe, Si, Ge, etc. Individual microlenses are also called subapertures and the SHWS spots formed at the focal plane where the detector is placed are referred to as “spots” in this chapter.


Introduction
Wavefront sensor is a device that measures the optical wavefront aberration.The Shack Hartmann wavefront sensor (SHWS), named after Johannes Franz Hartmann and Roland Shack, is one of the most often used optical wavefront sensor.It is made up of an array of microlenses (all having the same focal length and aperture size) and a detector placed at the focal plane of these microlenses.Johannes Franz Hartmann developed a device that consisted of an opaque screen with multiple holes to test the quality of imaging in large telescope systems and fine-tune the telescope focus (Hartmann, 1900).In the process of developing an adaptive optics system to improve the resolution of satellite images taken from the earth, Roland Shack came up with a feasible model of the sensor by using an array of tiny lenses instead of holes (Shack & Platt, 1971).The technological advancements in the field of optical fabrication allowed the industry to make lenses as small as 100 µm using materials like fused silica, ZnS, ZnSe, Si, Ge, etc. Individual microlenses are also called subapertures and the SHWS spots formed at the focal plane where the detector is placed are referred to as "spots" in this chapter.
SHWS is widely used in diverse wavefront sensing applications.It is very commonly used in astronomical adaptive optics systems (Gilles & Ellerbroek, 2006), lens testing (Birch et al., 2010), ophthalmology (Wei & Thibos, 2010) and microscopy (Cha et al., 2010).It is used in the correction of errors due to non-flatness of spatial light modulators in holographic optical tweezer applications (López-Quesada et al., 2009).A few modifications over the simple SHWS were also suggested in the literature.A hexagonal arrangement of the lenslet array can increase the sensitivity and dynamic range of the detector (Wu et al., 2010).A differential SHWS was proposed which measures the wavefront slope differentials (Zou et al., 2008).In order to make a dynamic microlens array with greater dynamic range, digital SHWS were developed with the help of Liquid Crystal-Spatial Light Modulators (LC-SLMs) (Zhao et al., 2006).It is also possible to characterize atmospheric turbulence using the SHWS via the structure function measurements, C 2 N profile measurements, wind velocity profile estimation and measurement of deviations from the theoretical model of Kolmogorov turbulence (Sergeyev & Roggemann, 2011;Silbaugh et al., 1996).Also, using the SHWS is advantageous over the curvature sensor when large numbers of sensing elements (subapertures) are to be used (Kellerer & Kellerer, 2011) and hence is the best option in large telescope adaptive optics.
The simplest model of a SHWS is shown in Fig. 1.When a plane wavefront is incident on the SHWS, it produces equidistant focal spots at the detector plane.Any wavefront distortion Fig. 1.Description of SHWS: Phase screen (left) following the Kolmogorov turbulence model of the atmosphere, simulated using the Fast Fourier Transform (FFT) technique.Simulated spot pattern with no turbulence assuming a SHWS with 10 × 10 subapertures is shown in the center and the distorted spot pattern image simulated by taking the turbulence into account isshownontheright.
introduced in the incident wavefront will displace the spots from their original locations.The distance moved by individual spots contains the information of the local wavefront slopes of the incoming distorted wavefront.The retrieval of the wavefront shape in a SHWS is a two step process.In the first step, the positions of the spots corresponding to individual subapertures is determined through spot centroiding techniques and the local slopes, (β x , β y ) of the wavefront are subsequently determined by calculating the shift in the spots from a reference location (calculated from image of the spot pattern captured when no distortion is present).As a second step, the shape of the wavefront is reconstructed from the measured local slope values.
The simplest method of locating the spot is to identify the pixel with maximum intensity.In this technique of peak identification, the finite number of detector pixels per subaperture limits the accuracy of locating the position of the spot.In order to accurately determine the location of the spots, we need to calculate the location of the spot in subpixel units or increase the resolution of imaging.Dealing with high resolution images means more readout time and greater computational effort and hence should be avoided in applications requiring high-speed operation.The slopes are related to the centroid locations in the following manner, where R =(x r , y r ) is the reference coordinate and P =(x p , y p ) is the position within a single subaperture to which the spot is moved after distortion of the wavefront, ∆x = x p − x r and ∆y = y p − y r represent the magnitude of shifts along 'x' and 'y' directions and ' f 'isthefocal length of the microlenses.
Since most of the wavefront reconstruction error arises from inaccurate centroiding, it is a very important step.Centroiding in the case of astronomical adaptive optics is even more challenging because the science objects under observation are weak light sources and hence need longer exposure times (several seconds to a few minutes).This even makes them inappropriate for wavefront sensing, which needs to be done at a much faster rate.For this reason, another star (called Natural Guide Star or NGS) in the close vicinity of the science target is used as a reference source.In the absence of a NGS, an artificial star (called Laser Guide Star or LGS) is generated by shining a high power laser into the atmosphere, and is used as a reference source.
LGS has become an integral part of adaptive optics systems for better sky coverage (Fugate et al., 1991;Primmerman et al., 1991).
LGS can be of two types -Rayleigh or sodium beacon.Rayleigh beacon is formed by the light scattered from molecules at lower altitudes, ranging from 10 km to 15 km depending on the power of the laser and the site conditions (Thompson & Castle, 1992).This low altitude reference source would lead to under sampling and hence the problem of focal anisoplanatism arises (Fried, 1982).A high power 589 nm laser can be used to excite the meteoric origin sodium atom layers in the mesosphere which are present at a mean altitude of 92 km and with a mean thickness of ∼10 km (Hart et al., 1995).The de-excitation of the atoms from the upper states (3p 3 2 and 3p 1 2 ) to the lower state (3s 1

2
)via spontaneous decay results in resonant backscattering and hence an artificial star.Increasing the power of the laser enhances the chances of stimulated emission against spontaneous decay and thereby reducing the intensity of the desired backscattered photons, eventually leading to population inversion and medium saturation.To avoid this problem, an optimum laser power is used.The limit on the power of LGS also limits the number of available photons for wavefront sensing.The number of available photons, p N decides the centroiding accuracy and the minimum exposure time which in turn controls the adaptive optics servo bandwidth (Hardy, 1998).
The other noise concerns include background noise due to Rayleigh scattering of laser light and under sampling of the spot due to servo bandwidth constraints.In the case of laser guide star (LGS) based sensing, the elongation of the spots in large telescopes and the variability in sodium density profile cause further errors.Readout noise of the detectors in addition to the photon noise may also seriously degrade centroiding accuracy in a few cases.The following section discusses a few basic centroiding algorithms discussed in the literature.

Center of gravity
The method of center of gravity (CoG) calculates the centroid location as the weighted mean of the position coordinates, (x, y) and the weight being the spot intensity as a function of position coordinates, I(x, y)=I(X, Y).Here, X, Y are used to denote discreteness.The centroid, (x s , y s ) of a single subaperture spot pattern, (I) is evaluated using, where i, j are row and column indices running from 1 to M,w h e r eM × M is the size of a single subaperture matrix.This is the simplest of all the centroiding techniques and is best suited to situations where the light intensity levels are sufficiently high and the signal to noise ratio (SNR) is good enough.

169
Advanced Methods for Improving the Efficiency of a Shack Hartmann Wavefront Sensor www.intechopen.com

Weighted center of gravity
If the shape of the spot pattern is known beforehand, then this method takes advantage of this additional information for accurate determination of the location of the centroid (Fusco et al., 2004).The mathematical form that is assumed for the shape of the spot is called the weighting function and is multiplied with the intensity function before applying the CoG algorithm as previously discussed.The estimated centroid location becomes: The weighting function W(x, y) is generally assumed to be a Gaussian function when a natural guide star (NGS) is used as the reference in wavefront sensing: The spots cannot be assumed to be Gaussian in a digital SHWS due to the domination of diffractive noise and instability of the light source.Hence, an adaptive thresholding and dynamic windowing method is used in this case for accurate centroid detection (Yin et al., 2009).Also, in the case of LGS adaptive optics system, the spot cannot be assumed to be Gaussian for large telescope systems due to the problem of spot elongation.In this case, the weighting function, W(x, y) can be simulated by assuming ideal conditions (no turbulence) (Vyas et al., 2010c).
This algorithm is best suited in the closed loop adaptive optics systems where the shift in the spots over consecutive temporal measurements is small.It is not suitable for large shift in the spots and hence inappropriate to open loop systems and large phase errors (Vyas et al., 2009b).

Intensity Weighted Centroiding
Intensity Weighted Centroiding (IWC) is similar to WCoG with a difference that the weighting function, W ij is the intensity distribution of the spot pattern, I ij .Hence, in IWC, the estimated centroid position becomes, In comparison to the CoG method, this algorithm performs a better job under low light level conditions and low background and readout noise.

Iteratively weighted center of gravity
The problem of inaccurate centroiding in the case of large shift in the spots can be overcame by using the iteratively weighted center of gravity (IWCoG) algorithm where the centroid location is computed iteratively (Baker & Moallem, 2007).The weighting function is modified 170 Topics in Adaptive Optics www.intechopen.comafter each iteration and is centered around the new centroid location, (x n s , y n s )identifiedbythe n th iteration.The centroid location in the n th (n ≥ 2) iteration is defined as, where The first (initial) iteration uses the weighting function defined in Eq. 4. The width of the Gaussian can also be modified after each iteration, but in general has a very little effect on the accuracy of centroid estimation if the same optimal width is used in all the iterations.Any iterative process carries along with it problems like -saturation of performance, non uniform convergence and speed.These issues are answered by a hybrid algorithm (Vyas et al., 2010b), a combination of the IWCoG and the correlation technique, which is discussed in the coming sections.

Matched Filter algorithm
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).

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'andN × N subapertures, The application of least square curve fitting model, on the SHWS slope data, reduces the problem into a matrix form, DS = Aφ,whereD and A are sparse matrices of sizes 2N 2 × 2N 2 and 2N 2 × N 2 respectively.S is a vector containing the measured slope values, 'x's l o p e s followed by 'y' slopes.The wavefront phase can therefore be evaluated using the following expression, To solve for φ using Matlab, we can apply the in-built command, LSQR on the linear system of equations, Aφ = B,whereB(= DS) is a vector of size 2N

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, Applying the Fourier transform and using the shift property of the discrete Fourier transform on Eq. 10 gives, 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).

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.

172
Topics in Adaptive Optics www.intechopen.com

Temporally evolving phase screens
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, v w (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 layers 1 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 (X l 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 (P l 1 (t = 0),'t' represents time) in the temporal evolution of this particular layer.Subsequent phase screens (P l i (t), i=2,3,...n) at later times are formed by translating the initially selected portionontheverylargewavefrontX l 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, X l .I ft h es i z eo f P l 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 = d l / v l ,w h e r ev l is the layer velocity and d l 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 d l (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, X l 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 P l i (t) on X l 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  (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).

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, 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.

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 (N S ), 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 (N P ), the detector readout noise and its quantum efficiency (Hardy, 1998).As 'N S ' 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 'N P ' 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.

Working model
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 − ,p i t c h=6 .4 5 µ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.8nm)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 (Z n 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).

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 weighted center of gravity, thresholded Zernike reconstructor based centroiding and Gaussian pattern matching algorithm.

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, 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  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 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, 'h i−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 where (x Act , y Act ) is the actual centroid position of the spot that is although unknown in real situations, is predefined while performing Monte Carlo simulations.

Error saturation
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 n th iteration, the new weighting function for the (n + 1) th iteration is calculated such that it is centered on a modified centroid location, (x n R , y n R ) instead of (x n c , y n c ) such that, where, r n x and r n y are the generated pseudo random numbers with zero mean and unit variance, which assumes different values after each iteration.N x and N y are normalization constants (crudely they define the magnitude of randomness introduced; in our case N x = N y = 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.

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.A simple technique based on iterative computation of correlation can be used to identify the iteration with minimum error.After the n th iteration, an ideal spot pattern (simulated using PSF with atmospheric turbulence excluded and with zero noise) is calculated with its center at (x n c , y n 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).The maximum number of iterations (n max ) 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 n max > 1i s 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.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.

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.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.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).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.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.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.

Gaussian pattern matching
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.

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 scintillations, we impose threshold conditions on the size of the features.

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.

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).

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. 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.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;  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).

Fig. 2 .
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.

175AdvancedFig
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.

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

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

Fig. 10 .
Fig. 10.Saturation and non uniform convergence of IWCoG algorithm at different SNR Fig. 11.The effect of error saturation can be overcome by using the IARN based IWCoG algorithm.This procedure largely reduces the CEE.

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. 16 .Fig. 18 .
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.

183Advanced
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.

Fig. 21 .
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.
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

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

Fig. 24 .
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.
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 189 Advanced Methods for Improving the Efficiency of aShack Hartmann Wavefront Sensor   www.intechopen.com

Fig. 29 .
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. 32 .Fig. 33 .
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