## Abstract

We propose a spline-based aberration reconstruction method through moment measurements (SABRE-M). The method uses first and second moment information from the focal spots of the SH sensor to reconstruct the wavefront with bivariate simplex B-spline basis functions. The proposed method, since it provides higher order local wavefront estimates with quadratic and cubic basis functions can provide the same accuracy for SH arrays with a reduced number of subapertures and, correspondingly, larger lenses which can be beneficial for application in low light conditions. In numerical experiments the performance of SABRE-M is compared to that of the first moment method SABRE for aberrations of different spatial orders and for different sizes of the SH array. The results show that SABRE-M is superior to SABRE, in particular for the higher order aberrations and that SABRE-M can give equal performance as SABRE on a SH grid of halved sampling.

© 2017 Optical Society of America

## 1. Introduction

In the field of adaptive optics (AO), Shack-Hartmann (SH) sensors are commonly used to estimate the wavefront for turbulence induced aberration compensation. The microlens array of the SH sensor creates a focal spot pattern and the first moment of the intensity distribution of each focal spot gives an approximate of the local averaged spatial gradient of the wavefront aberration.

The most used approaches for wavefront reconstruction (WFR) that process the local slope measurements from a SH sensor are the zonal (local) finite difference (FD) method [1] or modal methods that are based on a set of (global) basis functions for the wavefront representation [2]. A recently developed method for WFR was presented by de Visser et al. in [3]. The SABRE (Spline based ABerration REconstruction) method uses bivariate simplex splines [4] to compute an estimate of the wavefront from the local wavefront slope measurements of a SH sensor. This way the linearity in the phase estimation problem is preserved, while the use of modal functions in this zonal reconstruction method, represented by a set of local nonlinear spline basis functions, makes it possible to obtain higher order local approximations of the wavefront.

Simulations have shown that this method, compared to the classical finite difference method [1], is highly resilient to sensor noise and invariant of wavefront sensor geometry. In addition, the local nature of SABRE allows an implementation in a distributed way, that can significantly increase computational efficiency [5]. However, since [3] uses only the first moment information from the SH sensor, the potential of obtaining higher order local wavefront estimates with the nonlinear spline functions is not achieved.

Since each focal spot in a SH pattern is an intensity distribution containing more information than just the average slope, this work investigates if extra information from the sensor can be used in order to exploit this potential. In this paper, an extension of the SABRE is presented, the SABRE-Moments (SABRE-M), that uses additional higher order information from the focal spots of the SH sensor, for more accurate wavefront reconstruction. For this purpose, wavefront sensorless techniques—modal methods that use the complete intensity distribution from the science camera for nonlinear WFR—were considered for an application to each of the individual intensity distributions from the subapertures of the SH microlens array. Recently, in the works of Booth [6], Linhai and Rao [7] and Yang *et al.* [8], the linear relation between the change of the *second moment of intensity* of the image and the averaged squared spatial gradient of the wavefront aberration was proven. Based on this relation, a novel sensor model for the SABRE model is presented in this work that applies the first and second moment information from the SH sensor. This allows modelling of the wavefront with bivariate simplex splines of higher polynomial degree.

The SABRE-M method, that is introduced in this work, has some important benefits. Firstly, this method provides a more accurate wavefront reconstruction, especially for higher order aberrations. Secondly, the method makes it possible to reduce the number of subapertures in the SH grid without losing reconstruction accuracy compared to SABRE for the original grid. As a consequence the increased subaperture size can have beneficial effects on the dynamic range, the sensitivity and the signal to noise ratio. For example in the case of low light level conditions, caused be a faint natural or artificial guide star, e.g. in astronomy [9], ophthalmology [10], or live-cell fluorescence microscopy [11], increasing the subaperture size could be advantageous [12].

The organisation of this paper is as follows. In Section 2, preliminaries on the SABRE and the second-moments technique of [8] are given followed by a motivation for the suggested combination of both techniques. Then, the novel moment-based sensor model is derived and presented in Section 3. In Section 4, the results from numerical experiments are shown and discussed. Finally, Section 5 concludes this paper, along with recommendations for further research.

## 2. A combination of the SABRE and a second moment technique

#### 2.1. SABRE: a spline-based wavefront reconstruction method

The SABRE [3] is a wavefront reconstruction method estimating the phase aberration *ϕ*(**x**) from spatial gradient measurements, which are approximately given by the change of the first moments of the focal spots of a SH array. In the *x* direction, the following relation is applied:

*I*(

_{x}*n*) and

*I*(

_{y}*n*) denote the first moment of the

*n*-th subaperture of the SH sensor, in the directions

*x*and

*y*respectively, using the notation for the pupil plane coordinates

**x**= (

*x, y*) ∈ ℝ

^{2}, and

**x**

*defines the location of the*

_{n}*n*-th subaperture center.

The method uses bivariate simplex splines for the modelling of the wavefront, which are defined in the barycentric coordinate system according to [4]. Let *b*(**x**) be the barycentric coordinates of **x** ∈ *t*, for some given simplex *t*, then, as shown in [3], any polynomial *p*(*b*(**x**)) of degree *d* can be written as a linear combination of basis polynomials:

*t*, and ${B}_{\kappa}^{d}(b(\mathbf{x}))\in \mathbb{R}$ are the Bernstein basis polynomials for multi-index

*κ*= (

*κ*

_{0},

*κ*

_{1},

*κ*

_{2}) ∈ ℕ

^{3}. For a given degree

*d*, the possible number of combinations of

*κ*

_{0},

*κ*

_{1}and

*κ*

_{2}that satisfy |

*κ*| =

*d*, determines the total number of basis polynomials, which is given by $\widehat{d}=\frac{(d+2)!}{2d!}$.

The wavefront is approximated on a triangulation
$\mathcal{T}$ by a *global* simplex B-spline polynomial
${s}_{d}^{r}(b(\mathbf{x}))$, that consists of local polynomials *p*(*b*(**x**)) of degree *d*, and has a predefined continuity order *r* at the edges of adjacent triangles. As shown in [3], the wavefront model defined on a triangulation
$\mathcal{T}$ consisting of *J* triangles is written in vector notation as

*t*in $\mathcal{T}$.

The WFR problem consists of the sensor slope model of Eq. (1), expressed in terms of the spline model from Eq. (3) through the derivative of a B-form polynomial, and a set of global constraints

and is solved in the least squares sense. The constraint matrix $\mathbf{A}:={\left[{\mathbf{H}}^{\top}{\mathbf{h}}^{\top}\right]}^{\top}\in {\mathbb{R}}^{(EV+1)\times J\widehat{d}}$ is introduced in [3] as a combination of the*smoothness matrix*$\mathbf{H}\in {\mathbb{R}}^{EV\times J\widehat{d}}$, by which continuity of order

*r*is imposed, and of the

*anchor vector*$\mathbf{h}\in {\mathbb{R}}^{1\times J\widehat{d}}$, which is used to fix the unknown piston mode.

Different triangulation types were introduced in [3], of which the regular Type-II triangulation is used in this work. To obtain the SABRE model for SH slope measurements, the triangulation $\mathcal{T}$ is constructed on the (reference) centers of the subapertures, i.e. the locations in the SH array for which gradient information is provided through the sensor model in Eq. (1). Figure 1(a) shows an example of a Type-II triangulation on a 3 × 3 SH array.

The SABRE has the approximation power to obtain a very accurate estimate of the wavefront, due to the local modelling with higher degree (*d* > 1) basis polynomials. However, it is limited in two ways by the fact that only first order information of the focal spots is processed: firstly, only gradient information is extracted from the SH patterns, and therefore no use can be made of the potential of higher degree modelling. Secondly, we can only use up to second degree basis polynomials (*d* ≤ 2), in order to guarantee the resulting system equations from Eq. (1) are fully determined given the constraints in Eq. (4). The use of higher degree (*d* ≥ 3) basis polynomials will not result in a unique solution, since the system is under-determined due to the lack of data.

Therefore, this work proposes to extract additional higher order information from the individual SH sensor focal spots, by deriving a novel sensor model. A second-moments technique, presented in the work of Yang *et al.* [8], is considered for this purpose.

#### 2.2. A second moment technique

Commonly, WFR methods for SH sensors are based on the well known relation between the change of the first moments *I _{x}* (

*n*),

*I*(

_{y}*n*) of an intensity distribution and the averaged gradient of the prevailing wavefront. In [8], a new linear relation was proven, stating that the normalized change of the second moments ${I}_{{x}^{2}}$, ${I}_{{y}^{2}}$ of an image is proportional to the averaged squared gradient of the wavefront aberration:

*I*(

**u**) and

*I*

_{0}(

**u**) define the intensity at a certain location in the image plane

**u**= (

*u, v*) ∈ ℝ

^{2}, for the aberrated wavefront

*ϕ*(

**x**) and an unaberrated wavefront respectively, and

*P*(

**x**) denotes the pupil function which is one inside the subaperture and zero elsewhere.

The relation for the SH model in Eq. (1) is based on the principle that a tip and tilt in the wavefront gives a displacement of the focal spot in the image plane, *i.e.* the change of the first moments *I _{x}* (

*n*),

*I*(

_{y}*n*) of intensity. The change of the second moment of intensity gives the change of the averaged width of the intensity distribution. In the novel sensor model for the SABRE, presented in the next chapter, three additional relations are derived for the change of the higher order moments of the SH focal spots: two relations for the second moments ${I}_{{x}^{2}}(n)$, ${I}_{{y}^{2}}(n)$ in the

*x*and

*y*direction, and one relation for a mixed moment

*I*(

_{xy}*n*), showing the cross-correlation between the two second moments.

## 3. Wavefront reconstruction with the SABRE-Moments

In this section we derive the SABRE-Moments WFR problem. The linear phase-moment relationships for all considered moments *I _{x}* (

*n*),

*I*(

_{y}*n*), ${I}_{{x}^{2}}(n)$, ${I}_{{y}^{2}}(n)$ and

*I*(

_{xy}*n*) are presented in Section 3.1. After the introduction of closed form expressions for derivation and integration of B-spline polynomials in Section 3.2, the SABRE-M sensor model in terms of B-coefficients is derived in Section 3.3. The section closes with the LS problem formulation of the SABRE-M method.

#### 3.1. Principle of a moment-based Shack-Hartmann sensor model

The novel moment-based SH sensor model consists of five equations, two for the first moment measurements *I _{x}* (

*n*) and

*I*(

_{y}*n*), two for the second moment measurements ${I}_{{x}^{2}}(n)$ and ${I}_{{y}^{2}}(n)$ and one for the mixed moment measurement

*I*(

_{xy}*n*), defined in the

*n*-th subaperture of the SH sensor.

Based on the analytical expression for the change of the first moments of intensity *I _{x}* (

*n*) and

*I*(

_{y}*n*), which states that the change of the focal spot centroid is proportional to the average of the wavefront spatial gradient, the system of equations (1), which relates the change of the first moment of intensity by approximation to the local spatial gradient in

**x**

*, is replaced by an integral over the subaperture:*

_{n}*x*, where

*P*(

_{n}**x**) is the square pupil function of subaperture

*n*, which is one inside the subaperture and zero elsewhere, and ${c}_{1}:=\frac{1}{2\pi}\frac{1}{{\displaystyle {\int}_{{P}_{n}(x)}1d\mathbf{x}}}$ a constant term including the division by the total power of the light inside the pupil.

The expressions for the change of the second moments
${I}_{{x}^{2}}(n)$ and
${I}_{{y}^{2}}(n)$ are based on the relation proven in the work of Yang *et al.* [8], relating the second moment measurements to the averaged squared spatial gradients:

*x*, with the constant term defined as ${c}_{2}:=\frac{1}{4{\pi}^{2}}\frac{1}{{\displaystyle {\int}_{{P}_{n}(x)}1d\mathbf{x}}}$. A simplification of the equations in [8] was achieved by considering

*central*second moments around the centroid. Therefore the squared centroid is subtracted on the right side in Eq. (7), causing the linear components of the sensor model to fall out.

Finally, the expression for the change of the mixed moment *I _{xy}* (

*n*), can be derived in a similar way as the second moment equations in [8], and is given by the cross-correlation between the second moments in the

*x*and

*y*direction:

*central*mixed moment is considered yielding the subtraction of the first moment product on the right side of Eq. (8).

These five equations form the moment-based SH sensor model which is the foundation for the extension of the SABRE presented in this work. An important difference to the SABRE is that the moment measurements are related to the complete respective subaperture by the integration in the Eqs. (6)–(8). As shown in Fig. 1(a), the standard SABRE triangulation is built on the vertices located at the SH subaperture centers. This geometry is needed in the SABRE using basis polynomials of degree *d* = 2 in order to obtain a well determined system. A triangulation that is more in line with the novel SABRE-M model is selected, such that the vertices coincide with the corners of the subapertures, which allows integration over the complete simplices, as shown in Fig. 1(b). The increased amount of data in the novel sensor model ensures a well determined system on this triangulation grid, even for basis polynomials of degree *d* = 3.

#### 3.2. The directional derivative and integral of a B-form polynomial

The derivation and integration of spline functions has to be discussed for the derivation of the SABRE-M sensor model.

The directional derivative of a B-form polynomial can be expressed in terms of the original vector of B-coefficients [13]. This relation is used for the modelling of the wavefront slopes in the SABRE in [3]. On a simplex *t*, the first order directional derivative in the direction of a Cartesian unit vector **e** ∈ ℝ^{2} of the B-form polynomial *p*(*b*(**x**)) introduced in (2) is given by

*Casteljau matrix*of degree

*d*to

*d*– 1, which is expressed in terms of the barycentric directional coordinate

**a**:=

_{e}*b*(

**v**) –

*b*(

**w**) ∈ ℝ

^{3}of unit vector

**e**=

**v**−

**w**∈ ℝ

^{2}with respect to the triangle

*t*[13]. Basis polynomial vector ${\mathbf{B}}_{t}^{d-1}(b(\mathbf{x}))\in {\mathbb{R}}^{1\times \widehat{d-1}}$ for triangle

*t*contains the basis polynomials of reduced degree

*d*− 1. Here $\widehat{d-1}$ defines the number of B-coefficients on a simplex for a spline of degree

*d*− 1 according to the definition of $\widehat{d}$ as given in Section 2.

In [4] and [14], a formulation of the integrals of B-form polynomials was presented. An explicit expression of integrals of the Bernstein basis polynomials over the complete area of a simplex *t* is given as

*d*and the area

*A*of simplex

_{t}*t*. The integration of B-form

*p*(

*b*(

**x**)) from (2) over simplex

*t*is then computed as the sum of the corresponding B-coefficients ${c}_{\kappa}^{t}$, |

*κ*| =

*d*multiplied by the right hand side of (10). Therefore, like the directional derivative, also integration over a simplex can be expressed in terms of the B-coefficient vector

**c**

*. This useful property of spline functions forms the basis for the modelling of averaged slopes in the novel sensor model for the spline based wavefront reconstruction method.*

^{t}Further, the inner product of two basis polynomials integrated over simplex *t* is given by [14]

*κ*! =

*κ*

_{0}!

*κ*

_{1}!

*κ*

_{2}! is used for the factorial of the multi-index

*κ*, and $\widehat{{d}_{1}+{d}_{2}}$ defines the number of B-coefficients on a simplex for a spline of degree

*d*

_{1}+

*d*

_{2}according to the definition of $\widehat{d}$.

#### 3.3. Derivation of the SABRE-M model in spline coefficients

In this section the equations for the moment-based Shack-Hartmann sensor from Eqs. (6)–(8) are derived in terms of the spline functions, resulting in the system equations for the SABRE-M.

The wavefront is defined locally at each subaperture *n* = 1, 2,*…, ,N* in terms of a spline function:

*J*triangles on one subaperture domain, where the vertices are at the corners of the subaperture as shown in Fig. 1(b).

_{n}The averaged local gradient from Eq. (6) in *x* or *y* direction inside the subaperture *n* is calculated as the sum of the *J _{n}* integrals of the wavefront spatial gradients over each simplex
$t\in {\mathcal{T}}_{n}$ divided by the total area of the subaperture domain.

The integral over a simplex *t* of the derivative in direction *x* of the simplex B-spline wavefront model from Eq. (3), can be obtained using the definition of the derivative of a B-form polynomial from Eq. (9) and the expression of the integral in Eq. (10):

**e**in the derivative direction

_{x}*x*with respect to the triangle

*t*is denoted by

**a**.

_{x}Using (15), and summing over all *J _{n}* triangles in subaperture

*n*, the first moment equation with respect to

*x*for a subaperture

*n*of the sensor model from Eq. (6) can now be rewritten to

*de Casteljau*matrix for the directional coordinate

**a**with respect to the triangle

_{x}*t*. It can be noted that

_{j}*A*disappears in the equations, since the integral in (6) is divided by the total area of the subaperture $\int}_{{P}_{n}(x)}1d\mathbf{x}={J}_{n}{A}_{t$ in order to obtain the average slope.

_{t}Repeating the same for the *y* coordinate and accounting for the measurement noise and modelling errors, the linear part of the moment-based sensor model is obtained, which in terms of the local coefficient vector **c*** _{n}* is given by

*de Casteljau*matrices corresponding to the triangles

*t*= 1, 2,…

_{j}, j*J*in subaperture

_{n}*n*The terms

*η*(

_{x}**x**),

*η*(

_{y}**x**) ∈ ℝ contain sensor noise and modelling errors.

The second moment and the mixed moment equations are derived in a similar way as for the first moments, whereas now the averaged *squared* gradients from Eq. (7) are calculated.

The integral over a simplex *t* of the squared derivative in the *x* or *y* direction of the wavefront model from Eq. (3) is obtained using the definition of the derivative of a spline from Eq. (9) and the definition of the integral of the inner product of basis polynomials over a triangle *t* from Eq. (11):

*κ*| = |

*γ*| =

*d*− 1.

The central second moment equation with respect to *x* for a subaperture *n* of the sensor model from Eq. (7) can now be expressed in terms of the spline functions by summing of Eq. (18) for *J _{n}* triangles and subtracting the squared first moment from Eq. (16):

In an analogous manner, the equations for the second moment with respect to *y* and for the mixed moment from Eq. (8) are written in terms of the subaperture local coefficient vector **c*** _{n}*. The quadratic part of the moment-based sensor model is then obtained as

*J*blocks of the form ${\mathbf{P}}_{j}^{d-1,d}{\left({\mathbf{a}}_{{\mathbf{e}}_{1}}\right)}^{\top}{\mathbf{B}}_{{\mathbf{I}}_{\gamma ,\kappa}}^{d-1}{\mathbf{P}}_{j}^{d-1,d}({\mathbf{a}}_{{\mathbf{e}}_{2}})$ corresponding to triangles

_{n}*t*,

_{j}*j*= 1, 2,…,

*J*, in subaperture

_{n}*n*. The Cartesian unit vectors

**e**

_{1},

**e**

_{2}are defined according to the considered moment. The terms ${\eta}_{{x}^{2}}(\mathbf{x})$, ${\eta}_{{y}^{2}}(\mathbf{x})$,

*η*(

_{xy}**x**) ∈ ℝ contain sensor noise and modelling errors.

#### 3.4. The global problem formulation and nonlinear least squares estimator

In this section the nonlinear least squares estimator is presented for the B-coefficients of the global SABRE-M model, in order to obtain an approximation of the wavefront from the moment measurements.

The parameter estimation problem, like for the SABRE, consists of the minimisation of the error between the measurements and the sensor model in a least squares sense subjected to the continuity constraints.

Since the moments are only related to the splines inside their respective subaperture, a local least squares problem can be defined on each subaperture domain consisting of *J _{n}* triangles. For this to each subaperture

*n*= 1,…,

*N*a measurement vector

The local residual vector *r _{n}* (

**c**

*) ∈ ℝ*

_{n}^{5×1}is defined in the following manner. For each of the five moment measurements

**b**

*,*

_{n}*,*

_{m}*m*= 1,…, 5, a residual value is defined that is strictly linear or quadratic in ${\mathbf{c}}_{n}\in {\mathbb{R}}^{{J}_{n}\widehat{d}}$ for the first moments (

*m*= 1, 2) or the second and mixed moments (

*m*= 3, 4, 5) respectively:

*m*= 1, 2 and is zero for

*m*= 3, 4, 5, and ${\mathbf{c}}_{n}^{\mathrm{T}}{\mathbf{Q}}_{n,m}{\mathbf{c}}_{n}$ is given by the sensor model from Eq. (20) for

*m*= 3, 4, 5 and is zero for

*m*= 1, 2. It should be noted that the sensor model is identical for each of the subapertures due to the local nature of the splines.

The global SABRE-M problem for a total of *N* subapertures is written as the following nonlinear least squares problem:

In order to solve for the B-coefficients in Eq. (23), a Gauss-Newton algorithm is used and the constraints are imposed by KKT conditions [15].

## 4. Simulations with the SABRE-M

In numerical experiments the performance of the second moment (SM) method, the SABRE-M, presented in Section 3, is compared to the first moment (FM) method SABRE, that was described in Section 2.1. The purpose of these experiments is to validate the SM model and to evaluate the improvements that are achieved by the use of additional information from the SH sensor.

For all the experiments described in this section, a Fourier-optics-based SH WFS simulation according to the Fraunhoffer diffraction principle is used in order to generate the first and second moment measurements in each subaperture. Aiming for the application on SH arrays consisting of a smaller number of bigger subapertures, considered are a large pitch between the square subapertures of *M _{p}* = 32 pixels, a width of the diffraction limited PSF of

*M*= 4 pixels, a focal distance of

_{d}*f*= 20 mm and a wavelength of

*λ*= 633 nm. In the simulations an

*N*×

_{i}*N*= 10 × 10 microlens array is used. Each intensity spot of the SH pattern is created by embedding the segment of the wavefront corresponding to the individual lens in the center of a

_{j}*M*

_{CCD}×

*M*

_{CCD}grid containing zeros, where

*M*

_{CCD}is defined by the pitch times the diffraction limit, so here

*M*

_{CCD}=

*M*·

_{d}*M*= 128. The focal spot in the

_{p}*n*-th subaperture is then obtained with the zero-shifted two-dimensional discrete Fourier transform:

*k*= 2

*π/λ*, the pupil function ${P}_{n}\in {\mathbb{R}}^{{M}_{\text{CCD}}\times {M}_{\text{CCD}}}$ containing ones and zeros, with the corresponding wavefront segment ${\varphi}_{n}\in {\mathbb{R}}^{{M}_{\text{CCD}}\times {M}_{\text{CCD}}}$, and the sampling interval

*δ*. The complete SH pattern is obtained by cutting out each focal spot at the subaperture size and placing the images in the array, forming an

*M*×

*M*grid, where

*M*=

*N*·

_{i}*M*= 320.

_{p}In the first part of the experiments (Sections 4.1 and 4.2), besides the·performance through WFR from SH measurements, also the performance of the spline models is shown through a fit of the wavefront. This is showing us the differences between the modelling ability with the used (first moment or first and second moment) information from the sensor and the potential modelling ability of the splines.

In order to evaluate the performance of the wavefront reconstruction, the relative root mean squared error (RMSE) is used as performance metric:

which is defined by the norm of the difference between the piston removed input phase*ϕ*and estimated phase $\widehat{\varphi}$, normalised by the root mean square (RMS) of the input phase.

#### 4.1. Proof of concept

The first part of the experiments is a proof of concept of the new sensor model derived in Section 3. In order to evaluate the gain in reconstruction accuracy that can be achieved through the additional second moment information retrieved from the SH sensor by Eq. (20), the SABRE-M is compared to the SABRE, that only uses an approximation of Eq. (17) as sensor model. In the experiments, a spline model of polynomial degree *d* = 2 and continuity order *r* = 1 on a Type-II triangulation grid is used for both methods. The reconstruction accuracy is tested for aberrations of different spatial orders, for which deterministic Zernike modes up to the fifth order with an amplitude of *A* = 1.9 rad are used as input phase. Because the reconstruction is made on a square aperture, the Zernike polynomials are defined on a circle with the aperture diagonal as diameter.

In Fig. 2, the performance of the SABRE-M is compared to the SABRE, plotted for the different order Zernike aberrations denoted by
${Z}_{n}^{m}$. As a reference, also the approximation error of the best possible fit with the second degree spline model is given, showing the potential in modelling accuracy with the *d* = 2, *r* = 1 splines. A fitting error of zero can be observed for aberrations up to second order Zernike modes, since these orders of modes can be exactly reconstructed by the second degree splines. Increasing the order of the mode in general implies a decreasing fitting ability. An exception is seen for the specific modes
${Z}_{n}^{\pm n}$, which shape is relatively flat in the middle and extremely steep at the corners of the square reconstruction domain, worsening the fitting ability with the spline model. The WFR results show that the SABRE-M, by the use of higher order information from the SH sensor, gives a performance that is significantly superior to that of the SABRE for input phase Z-modes from the third order onwards. The SABRE-M method outperforms the SABRE by 10% to 70%. Further it can be observed that even though the reconstruction accuracy is dependent on the shape of the deterministic Zernike modes, overall the advantage of SABRE-M is increasing, from less pronounced improvements at the third order Z-modes to very pronounced improvements for the fifth order Z-modes.

In addition, it is seen that the SM method approaches the modelling accuracy of the best possible fit at the specific modes ${Z}_{n}^{\pm n}$. This clearly shows us the difference between modelling power of the spline degree and reconstruction power by the SH sensor information. It also gives an example of the limitation of the reconstruction accuracy by the degree of the spline model.

#### 4.2. Power of a higher degree approximation model

The performance of the wavefront reconstruction with simplex splines depends partly on the selection of the spline parameters. Depending on the size of the SH array, the input aberration and the amount of data provided by the method, a different spline degree, continuity order or triangulation type is desirable, as a trade-off exists between smoothness and degrees of freedom of the global spline. The highest spline degree that can be used by the SABRE is *d* = 2, as discussed in Section 2.1. Due to the use of additional second moment information, with the SABRE-M, a *d* = 3 spline model can be used, which contains more degrees of freedom and therefore increases the fitting ability of the wavefront model.

For the evaluation of the performance of the second moment method, stochastic results are used, with *d* = 2, *r* = 1 and *d* = 3, *r* = 1 spline models for the SABRE-M and a *d* = 2, *r* = 1 model for the SABRE, all constructed on a Type-II triangulation grid. Random input phase screens are created with Zernike functions according to a Kolmogorov turbulence statistics model [16]. For this model, a turbulence outer scale of *L*_{0} = 10*D* m and five different turbulence strengths are used, with the severity defined by the ratio of the telescope diameter and the Fried coherence length *D/r*_{0} = [5, 10, 20, 40, 60]. Furthermore, the number of orders of Zernike modes included in the phase screens is varied, using *n* = [2, 5, 10, 15, 20, 25, 30] polynomial orders, corresponding to a total number of Zernike modes of [5, 20, 65, 135, 230, 350, 495]. For each number of included orders of Zernike modes, sets of 30 wavefront realisations for each turbulence strength are used in the experiment, giving a total of 150 realisations per order. In Fig. 3, four examples are shown of the phase screens and their corresponding SH patterns for a low and a high order aberration, each at a weak and a strong aberration strength.

In Fig. 4 the averaged results from these simulations for aberrations of different spatial orders are shown. In comparison to SABRE(*d* = 2, *r* = 1), the SABRE-M (*d* = 2, *r* = 1) using quadratic B-splines shows a significant advantage in reconstruction accuracy for aberrations that include Zernike polynomials of spatial orders in the range of 5 – 25. As seen from the best fit for the *d* = 2 spline, the SABRE-M (*d* = 2, *r* = 1) follows the maximal modelling accuracy much closer than the SABRE (*d* = 2, *r* = 1), due to the additional second moment information. For aberrations including > 15 orders a better accuracy can be achieved by increasing the spline degree to *d* = 3, for which a higher amount of degrees of freedom in the model is demonstrated by the superior best fit of the cubic spline model. The SABRE-M (*d* = 3, *r* = 1) shows an almost constant relative reconstruction error of less than 0.03. For aberrations that include Zernike polynomials of spatial orders up to 25 an improvement compared to SABRE of 65% is obtained. However, for very low order aberrations including < 15 orders, there is not enough second moment information for the high amount of degrees of freedom of the cubic spline model, and SABRE-M (*d* = 2, *r* = 1) should be used.

In the experiments aberrations including orders ≥ 30 (495 modes) are not considered because the results are affected by simulation errors due to the limited number of samples in the CCD grid. Because the second moment method considers the square of the phase gradient, which oscillates twice faster than the gradient used by the first moment method, the simulation for the second moments reaches the limitation imposed by the Nyquist criterion earlier. This can explain the growing reconstruction error when including Zernike modes of order ≥ 30 for the SABRE-M (*d* = 3, *r* = 1), as seen in Fig. 4. Hence, in order to guarantee a fair comparison of both methods through numerical simulations, the results from aberrations of extremely high spatial orders were not considered, and aberrations of Zernike modes of the first 25 spatial orders are used in the following simulations. According to [17] and [18] this number of orders is sufficient for an accurate representation of the turbulent phase.

Figure 5 shows a cross-section of the results from Fig. 4 at the level of 25 included orders of Zernike modes, in which the reconstruction accuracy is plotted for input phase screens of different aberration strengths. The results, which also applies to all the other levels of included Zernike orders, show for the cases where SABRE-M outperforms SABRE, that the suitable SM method is always superior, independent of aberration strength. Also, it is observed that the methods are mostly not dependent on aberration strength, but that the performance is mainly determined by the included modes. This would allow selecting smaller SH arrays of bigger lenses, for reconstruction in the range of orders in which the second moment method performs best.

#### 4.3. Analysis on the number of subapertures

As discussed in the introduction, one of the main goals of the SM method is to achieve better performance with smaller SH arrays of bigger lenses. For this purpose, the performance of the method is tested for different sizes of the SH array.

Figure 6 shows the relative RMSE for a 10 × 10, a 15 × 15 and a 20 × 20 microlens array, in which the results for 30 phase screens of *D/r*_{0} = 40 and including Zernike modes of polynomial order 25 or smaller are seen. The reduced loss of performance of the SM method for a decreasing number of subapertures shows clearly the main advantage of the novel method. By the use of the higher order information from the SH sensor, the method is able to reconstruct at a coarser sampling of the wavefront. The results show that a doubled sampling of the grid is needed for SABRE (at 20 × 20), in order to achieve a comparable accuracy as the SABRE-M at 10 × 10.

#### 4.4. Analysis on measurement noise

In this part of the experiments, the influence of sensor readout noise on the performance is investigated. The noise that affects the SH wavefront sensor measurements is simulated through Gaussian-distributed white noise added to the intensity patterns used for the moment computation. Different signal to noise ratio (SNR) levels are defined by the ratio of the intensity and the noise variance in the decibel scale, where a SNR of 0 dB corresponds to a magnitude of the noise that is equal to that of the signal. The images are preprocessed by applying a threshold, in order to remove the biggest influence of the noise on the measurements.

Figure 7 shows the performance of the methods for 30 wavefront realisations at each SNR level, using phase screens of *D/r*_{0} = 40 and 25 orders of Zernike modes. The results show that the SM method loses performance in contrast to the extremely noise resilient SABRE method. It is seen that for high signal to noise ratio levels ≥ 24 the SABRE-M still clearly outperforms the SABRE. For SNR levels < 24, the advantage compared to the SABRE is strongly reduced. Further research is required to improve the SM measurements in the presence of noise, for this the work of [19] and [20] is suggested.

## 5. Conclusions

A new SH sensor model that uses first and second moment information of the focal spots for wavefront reconstruction with multivariate simplex B-spline basis functions is introduced. This new wavefront reconstruction method, the SABRE-M (Moment) method, can be seen as an extension of the SABRE method which is based on an approximate model of the change in the first moments of the focal spots, commonly referred to as centroids. The SABRE-M sensor model includes next to the exact equations for the two first moment measurements, three additional equations that relate the change of the second moments of the focal spot to the local averaged square gradient of the wavefront.

The SABRE-M is intended for more accurate wavefront reconstruction in particular in the presence of higher order aberrations. First-moment-based methods only give a measure of the averaged slopes in each subaperture. However the second moment measurements allow the sensing of higher order aberrations in the subapertures. Also, whilst the original SABRE method is restricted to the use of linear or quadratic B-spline polynomials because of the limited number of measurements, the SABRE-M can employ cubic polynomials enabling the modelling of higher spatial frequencies in the wavefront.

The two claims of additional information retrieval and increased approximation power are validated in numerical experiments with a Fourier-based simulation of a 10 × 10 SH array.

Considering Zernike modes of the first 5 polynomial orders as aberrations, the reconstruction accuracies obtained for SABRE and SABRE-M are compared, with both models using a quadratic B-spline model. With equal achievable approximation power, the SABRE-M outperforms the first-moment-based SABRE method for all considered modes, with the advantage ranging between 10% and 70%. Superiority due to the use of the second order measurements is most pronounced for Zernike modes of polynomial order 4 or higher.

To analyse the reconstruction performance of SABRE-M for quadratic and cubic B-spline models a Monte Carlo simulation was performed for random aberrations created with Zernike modes according to a Kolmogorov turbulence statistics model. Increasing the number of included polynomial orders, it was observed that the spatial frequencies in the aberrations determine which B-spline model is the optimal choice for the SABRE-M. For aberrations including Zernike modes of 15 or less polynomial orders the quadratic model outperforms the cubic model. The cubic model is superior for aberrations including more than 15 polynomial orders because in this range the SABRE-M truly benefits from the additional approximation power. With the adequate B-spline model the second-moment-based SABRE-M shows an improvement of 55% to 65% compared to the SABRE for the considered high spatial frequency aberrations including more than 5 polynomial orders.

To highlight the benefit of the SABRE-M for possible application in low light scenarios, SH arrays of different sizes where considered. It was shown that the standard first-moment-based SABRE method needs a SH array of 20 × 20 subapertures in order to achieve reconstruction accuracies comparable to the performance achieved with the SABRE-M on a much coarser 10 × 10 array. Including the second moment measurements permits the use of a smaller number of big SH lenses at equal performance, which will naturally increase the signal to noise ratio in the focal spots.

Finally, in numerical experiments the SABRE-M performance has shown to be sensitive to the influence of Gaussian noise. The standard procedure of thresholding did ease the effect, however for decreasing signal to noise ratio levels the advantage to SABRE reduces. Further research is required to improve the second moment measurements in noisy conditions.

We conclude that the SABRE-M is suitable for higher order wavefront reconstruction within the single subaperture domains of a SH array, allowing an application on SH grids with a reduced number of subapertures and an increased subaperture size without the loss of reconstruction accuracy, which reduces the scale of the wavefront reconstruction problem and creates favorable signal to noise ratio conditions.

## Funding

European Research Council (ERC), (339681). Dutch NOVA partnership in the European Strategy Forum on Research Infrastructures (ESFRI). Russian Ministry of Education (“5 in 100”).

## Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013) / ERC grant agreement n° 339681. The work of E. Brunner was supported in part by the Dutch NOVA partnership in the European Strategy Forum on Research Infrastructures (ESFRI). The work of O. Soloviev is partially funded by the program “5 in 100” of the Russian Ministry of Education, and by Flexible Optical B.V.

## References and links

**1. **D. L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. **67**, 370 (1977). [CrossRef]

**2. **W. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. **70**, 998 (1980). [CrossRef]

**3. **C. C. de Visser and M. Verhaegen, “Wavefront reconstruction in adaptive optics systems using nonlinear multivariate splines,” J. Opt. Soc. Am. A **30**, 82 (2013). [CrossRef]

**4. **M. J. Lai and L. L. Schumaker, *Spline Functions on Triangulations* (Camebridge University, 2007). [CrossRef]

**5. **C. C. de Visser, E. Brunner, and M. Verhaegen, “On distributed wavefront reconstruction for large-scale adaptive optics systems,” J. Opt. Soc. Am. A **33**, 817–831 (2016). [CrossRef]

**6. **M. J. Booth, “Wavefront sensorless adaptive optics for large aberrations,” Opt. Lett. **32**, 5 (2007). [CrossRef]

**7. **H. Linhai and C. Rao, “Wavefront sensorless adaptive optics: a general model-based approach,” Opt. Express **19**, 371–379 (2011). [CrossRef] [PubMed]

**8. **H. Yang, O. Soloviev, and M. Verhaegen, “Model-based wavefront sensorless adaptive optics system for large aberrations and extended objects,” Opt. Express **23**, 24587 (2015). [CrossRef] [PubMed]

**9. **K. Tyson and Robert, *Principles of Adaptive Optics* (Academic, Inc, 1991).

**10. **J. F. Castejón-Mochón, N. López-Gil, A. Benito, and P. Artal, “Ocular wave-front aberration statistics in a normal young population,” Vision Res. **42**, 1611–1617 (2002). [CrossRef] [PubMed]

**11. **X. Tao, J. Crest, S. Kotadia, O. Azucena, D. C. Chen, W. Sullivan, and J. Kubby, “Live imaging using adaptive optics with fluorescent protein guide-stars,” Opt. Express **20**, 15969–82 (2012). [CrossRef] [PubMed]

**12. **F. Roddier, *Adaptive Optics in Astronomy* (Camebridge University, 1999). [CrossRef]

**13. **C. C. de Visser, Q. P. Chu, and J. A. Mulder, “Differential constraints for bounded recursive identification with multivariate splines,” Automatica **47**, 2059–2066 (2011). [CrossRef]

**14. **H. J. Tol, C. C. de Visser, and M. Kotsonis, “Model reduction of parabolic PDEs using multivariate splines,” Int. J. Control **7179**, 1–16 (2016). [CrossRef]

**15. **A. Bjorck, *Numerical Methods for Least Squares Problems* (SIAM, 1996). [CrossRef]

**16. **R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. **66**, 207 (1976). [CrossRef]

**17. **G.-m. Dai, “Modal wave-front reconstruction with Zernike polynomials and Karhunen-Loève functions,” J. Opt. Soc. Am. A **13**, 1218 (1996). [CrossRef]

**18. **M. C. Roggemann and T. J. Schulz, “Algorithm to increase the largest aberration that can be reconstructed from Hartmann sensor measurements,” Appl. Opt. **37**, 4321 (1998). [CrossRef]

**19. **C. Leroux and C. Dainty, “Estimation of centroid positions with a matched-filter algorithm: relevance for aberrometry of the eye,” Opt. Express **18**, 1197 (2010). [CrossRef] [PubMed]

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