Nonlinear spline wavefront reconstruction through moment-based Shack-Hartmann sensor measurements

M. Viegers; E. Brunner; O. Soloviev; C. C. de Visser; M. Verhaegen

doi:10.1364/OE.25.011514

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:

\begin{matrix} I_{x} (n) \approx ϕ_{x}^{'} (x_{n}), & n = 1, \dots, N, \end{matrix}

where 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) ∈ ℝ², and x_n defines the location of the 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:

p (b (x)) = {\begin{array}{l} \sum_{| κ | = d} c_{κ}^{t} B_{κ}^{d} (b (x)), & x \in t \\ 0, & x \notin t, \end{array}

where

c_{κ}^{t} \in ℝ

are called the B-coefficients, which define the unique local function on the triangle t, and

B_{κ}^{d} (b (x)) \in ℝ

are the Bernstein basis polynomials for multi-index κ = (κ₀, κ₁, κ₂) ∈ ℕ³. For a given degree d, the possible number of combinations of κ₀, κ₁ and κ₂ that satisfy |κ| = d, determines the total number of basis polynomials, which is given by

\hat{d} = \frac{(d + 2)!}{2 d!}

.

The wavefront is approximated on a triangulation $T$ by a global simplex B-spline polynomial $s_{d}^{r} (b (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 $T$ consisting of J triangles is written in vector notation as

\begin{matrix} ϕ (x) \approx s_{d}^{r} (b (x)) = B^{d} (b (x)) c, & x \in T \end{matrix},

where

B^{d} (x) \in ℝ^{1 \times J \hat{d}}

and

c \in ℝ^{J \hat{d} \times 1}

respectively are the vectors containing the individual basis polynomials

B_{κ}^{d} (b (x))

and the corresponding B-coefficients

c_{κ}^{t}

for all triangles t in

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

Ac = 0

and is solved in the least squares sense. The constraint matrix

A : = {[H^{⊤} h^{⊤}]}^{⊤} \in ℝ^{(E V + 1) \times J \hat{d}}

is introduced in [3] as a combination of the smoothness matrix

H \in ℝ^{E V \times J \hat{d}}

, by which continuity of order r is imposed, and of the anchor vector

h \in ℝ^{1 \times J \hat{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 $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.

Fig. 1 Type-II triangulations (black lines) for first moment (left) and second moment (right) measurements of a 3 × 3 SH array (blue lines) for the SABRE and the moment-based SABRE.

Download Full Size | PDF

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_{x^{2}} : = \frac{\int_{ℝ^{2}} (I (u) - I_{0} (u)) u^{2} d u}{\int_{ℝ^{2}} I (u) d u} = \frac{1}{4 π^{2}} \frac{\int_{ℝ^{2}} P^{2} (x) ϕ_{x}^{' 2} (x) d x}{\int_{ℝ^{2}} P^{2} (x) d x},

where I(u) and I₀(u) define the intensity at a certain location in the image plane u = (u, v) ∈ ℝ², 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_n, is replaced by an integral over the subaperture:

\begin{matrix} I_{x} (n) = c_{1} \int_{P_{n} (x)} ϕ_{x}^{'} (x) d x, & n = 1, 2, \dots, N, \end{matrix}

for the example of direction 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 π} \frac{1}{\int_{P_{n} (x)} 1 d 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:

\begin{matrix} I_{x^{2}} (n) = c_{2} \int_{P_{n} (x)} ϕ_{x}^{' 2} (x) d x - I_{x} {(n)}^{2}, & n = 1, 2, \dots, N, \end{matrix}

for the example of the direction x, with the constant term defined as

c_{2} : = \frac{1}{4 π^{2}} \frac{1}{\int_{P_{n} (x)} 1 d 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:

\begin{matrix} I_{x y} (n) = c_{2} \int_{P_{n} (x)} ϕ_{x}^{'} (x) ϕ_{y}^{'} (x) d x - I_{x} (n) I_{y} (n), & n = 1, 2, \dots, N . \end{matrix}

Following the same reasoning as given previously for the second moment equations from Eq. (7), the 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 ∈ ℝ² of the B-form polynomial p(b(x)) introduced in (2) is given by

D_{e}^{1} p (b (x)) = \frac{d!}{(d - 1)!} B_{t}^{d - 1} (b (x)) P^{d - 1, d} (a_{e}) c^{t},

where

P^{d - 1, d} (a_{e}) \in ℝ^{\hat{d - 1} \times \hat{d}}

is the Casteljau matrix of degree d to d – 1, which is expressed in terms of the barycentric directional coordinate a_e := b(v) – b(w) ∈ ℝ³ of unit vector e = v − w ∈ ℝ² with respect to the triangle t [13]. Basis polynomial vector

B_{t}^{d - 1} (b (x)) \in ℝ^{1 \times \hat{d - 1}}

for triangle t contains the basis polynomials of reduced degree d − 1. Here

\hat{d - 1}

defines the number of B-coefficients on a simplex for a spline of degree d − 1 according to the definition of

\hat{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

\int_{t} B_{κ}^{d} (b (x)) d x = \frac{A_{t}}{\hat{d}},

which depends only on the polynomial degree d and the area A_t of simplex 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_{κ}^{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^t. 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.

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

\int_{t} B_{κ}^{d_{1}} (b (x)) B_{γ}^{d_{2}} (b (x)) d x = \frac{d_{1}! d_{2}!}{(d_{1} + d_{2})!} \frac{(γ + κ)!}{γ! κ!} \frac{A_{t}}{\hat{d_{1} + d_{2}}},

where shorthand notation κ! = κ₀!κ₁!κ₂! is used for the factorial of the multi-index κ, and

\hat{d_{1} + d_{2}}

defines the number of B-coefficients on a simplex for a spline of degree d₁ + d₂ according to the definition of

\hat{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:

\begin{matrix} ϕ_{n} (x) \approx B_{n}^{d} (b (x)) c_{n}, & x \in T_{n}, \end{matrix}

with local B-form matrix

B_{n}^{d} (b (x)) \in ℝ^{1 \times J_{n} \hat{d}}

and local coefficient vector

c_{n} \in ℝ^{J_{n} \hat{d} \times 1}

, where

T_{n}

is a sub-triangulation [5] of

T

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

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

\int_{t} ϕ_{x}^{'} (x) d x = \int_{t} \frac{d!}{(d - 1)!} B_{t}^{d - 1} (b (x)) P^{d - 1, d} (a_{x}) c^{t} d x

= \frac{d!}{(d - 1)!} \frac{A_{t}}{\hat{d - 1}} 𝟙 P^{d - 1, d} (a_{x}) c^{t}

= \frac{2 A_{t}}{d + 1} 𝟙 P^{d - 1, d} (a_{x}) c^{t},

where

𝟙 \in ℝ^{1 \times \hat{d - 1}}

is a row vector containing only elements equal to one. The directional coordinate of unit vector e_x in the derivative direction 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

I_{x} (n) = \frac{1}{2 π} \frac{2}{J_{n} (d + 1)} \sum_{j = 1}^{J_{n}} 𝟙 P_{j}^{d - 1, d} (a_{x}) c^{t_{j}},

with

P_{j}^{d - 1, d} (a_{x})

the de Casteljau matrix for the directional coordinate a_x with respect to the triangle t_j. It can be noted that A_t disappears in the equations, since the integral in (6) is divided by the total area of the subaperture

\int_{P_{n} (x)} 1 d x = J_{n} A_{t}

in order to obtain the average slope.

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

\begin{array}{l} I_{x} (n) = \frac{1}{2 π} \frac{2}{J_{n} (d + 1)} I_{1, n} (a_{x}) c_{n} + η_{x} (x), \\ I_{y} (n) = \frac{1}{2 π} \frac{2}{J_{n} (d + 1)} I_{1, n} (a_{y}) c_{n} + η_{y} (x), \end{array}

with

I_{1, n} (a_{e}) : = [𝟙 P_{1}^{d - 1, d} (a_{e}) \dots 𝟙 P_{J_{n}}^{d - 1, d} (a_{e})] \in ℝ^{1 \times J_{n} \hat{d}}

a vector containing the sum of each column of the de Casteljau matrices corresponding to the triangles t_j, j = 1, 2,…J_n in subaperture 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):

\begin{matrix} \int_{t} ϕ_{x}^{' 2} (x) d x = \int_{t} {(\frac{d!}{(d - 1)!} B^{d - 1} (b (x)) P^{d - 1, d} (a_{x}) c^{t})}^{2} d x \\ = {(\frac{d!}{(d - 1)!})}^{2} c^{t ⊤} P^{d - 1, d} {(a_{x})}^{⊤} (\int_{t} B^{d - 1} {(b (x))}^{⊤} B^{d - 1} (b (x)) d x) P^{d - 1, d} (a_{x}) c^{t} \\ = {(\frac{d!}{(d - 1)!})}^{2} c^{t ⊤} P^{d - 1, d} {(a_{x})}^{⊤} \frac{A_{t}}{2 (\hat{d - 1})} B_{I_{γ, κ}}^{d - 1} P^{d - 1, d} (a_{x}) c^{t} \\ = \frac{d A_{t}}{(2 d - 1)} c^{t ⊤} P^{d - 1, d} {(a_{x})}^{⊤} B_{I_{γ, κ}}^{d - 1} P^{d - 1, d} (a_{x}) c^{t} . \end{matrix}

Using Eq. (11),

B_{I_{γ, κ}}^{d - 1} \in ℝ^{\hat{d - 1} \times \hat{d - 1}}

is defined as a matrix containing the scalars

\frac{(d - 1)! (d - 1)!}{(2 (d - 1))!} \frac{(γ + κ)!}{γ! κ!}

of each inner product

\int_{t} B_{κ}^{d - 1} (b (x)) B_{γ}^{d - 1} (b (x)) d x

, with |κ| = |γ| = 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):

I_{x^{2}} (n) = \frac{1}{4 π^{2}} (\frac{d}{J_{n} (2 d - 1)} \sum_{j = 1}^{J_{n}} c^{t ⊤} P_{j}^{d - 1, d} {(a_{x})}^{⊤} B_{I_{γ, κ}}^{d - 1} P_{j}^{d - 1, d} (a_{x}) c^{t} - {(\frac{2}{J_{n} (d + 1)} \sum_{j = 1}^{J_{n}} 𝟙 P_{j}^{d - 1, d} (a_{x}) c^{t})}^{2}) .

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

\begin{array}{l} I_{x^{2}} (n) = \frac{1}{4 π^{2}} c_{n}^{⊤} (\frac{d}{J_{n} (2 d - 1)} I_{2, n} (a_{x, x}) - \frac{4}{J_{n}^{2} {(d + 1)}^{2}} I_{1, n}^{⊤} (a_{x}) I_{1, n} (a_{x})) c_{n} + η_{x^{2}} (x), \\ I_{y^{2}} (n) = \frac{1}{4 π^{2}} c_{n}^{⊤} (\frac{d}{J_{n} (2 d - 1)} I_{2, n} (a_{y, y}) - \frac{4}{J_{n}^{2} {(d + 1)}^{2}} I_{1, n}^{⊤} (a_{y}) I_{1, n} (a_{y})) c_{n} + η_{y^{2}} (x), \\ I_{x y} (n) = \frac{1}{4 π^{2}} c_{n}^{⊤} (\frac{d}{J_{n} (2 d - 1)} I_{2, n} (a_{x, y}) - \frac{4}{J_{n}^{2} {(d + 1)}^{2}} I_{1, n}^{⊤} (a_{x}) I_{1, n} (a_{y})) c_{n} + η_{x y} (x), \end{array}

where

I_{2, n} (a_{e_{1}, e_{2}}) \in ℝ^{J_{n} \hat{d} \times J_{n} \hat{d}}

is a block diagonal matrix with a total of J_n blocks of the form

P_{j}^{d - 1, d} {(a_{e_{1}})}^{⊤} B_{I_{γ, κ}}^{d - 1} P_{j}^{d - 1, d} (a_{e_{2}})

corresponding to triangles t_j, j = 1, 2,…, J_n, in subaperture n. The Cartesian unit vectors e₁, e₂ are defined according to the considered moment. The terms

η_{x^{2}} (x)

,

η_{y^{2}} (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

b_{n} : = [I_{x} (n), I_{y} (n), I_{x^{2}} (n), I_{y^{2}} (n), I_{x y} (n)] \in ℝ^{5 \times 1}

is associated, consisting of the moment measurements obtained from the SH pattern.

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 $c_{n} \in ℝ^{J_{n} \hat{d}}$ for the first moments (m = 1, 2) or the second and mixed moments (m = 3, 4, 5) respectively:

\begin{matrix} {[r_{n} (c_{n})]}_{m} : = b_{n, m} - (l_{n, m}^{T} c_{n} + c_{n}^{T} Q_{n, m} c_{n}), & m = 1, \dots, 5, \end{matrix}

where

l_{n, m}^{T} c_{n}

is given by the sensor model from Eq. (17) for m = 1, 2 and is zero for m = 3, 4, 5, and

c_{n}^{T} Q_{n, m} 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:

\min_{c \in ℝ^{J \hat{d}}} f (c) = \sum_{n = 1}^{N} f_{n} (c_{n}) = \frac{1}{2} \sum_{n = 1}^{N} {‖ r_{n} (c_{n}) ‖}_{2}^{2}

\begin{matrix} s . t . & Ac = 0, \end{matrix}

where

c \in ℝ^{J \hat{d}}

denotes the global coefficient vector and

A \in ℝ^{E V + 1 \times J \hat{d}}

the global constraint matrix, which define the B-spline wavefront model introduced in Eqs. (3) and (4).

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_d = 4 pixels, a focal distance of f = 20 mm and a wavelength of λ = 633 nm. In the simulations an N_i × N_j = 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 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_p = 128. The focal spot in the n-th subaperture is then obtained with the zero-shifted two-dimensional discrete Fourier transform:

i_{n} = {| \frac{e^{i k (u^{2} + v^{2}) / (2 f)}}{i λ f} ℱ (P_{n} e^{- i ϕ_{n}}) δ^{2} |}^{2} \in ℝ^{M_{CCD} \times M_{CCD}},

with the wave number k = 2π/λ, the pupil function

P_{n} \in ℝ^{M_{CCD} \times M_{CCD}}

containing ones and zeros, with the corresponding wavefront segment

ϕ_{n} \in ℝ^{M_{CCD} \times M_{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_p = 320.

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:

R M S E_{r e l} = \frac{∥ ϕ - \hat{ϕ} ∥}{∥ ϕ ∥},

which is defined by the norm of the difference between the piston removed input phase ϕ and estimated phase

\hat{ϕ}

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

Fig. 2 Performance of the SABRE and SABRE-M on a 10 × 10 SH array, using a second degree spline, for different order Zernike input aberrations. The best fit of the aberration shows the maximal performance that can be obtained with a second degree spline.

Download Full Size | PDF

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₀ = 10D 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₀ = [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.

Fig. 3 Turbulent phase screens from Zernike modes of different polynomial order n and of different aberration strength D/r₀, and corresponding SH patterns. (a) and (b): weak aberration with RMS ≈ 3 rad; (c) and (d): strong aberration with RMS ≈ 20 rad.

Download Full Size | PDF

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.

Fig. 4 Performance of the SABRE and SABRE-M on a 10 × 10 SH array, using respectively a spline model of second and third degree. Stochastic results are shown, for input aberrations with Kolmogorov statistics created using Zernike modes up to different orders. The best fit of the aberration shows the maximal performance that can be obtained with a second, respectively third, degree spline.

Download Full Size | PDF

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.

Fig. 5 Cross-section of Fig. 4, showing the performance of SABRE and SABRE-M on a 10 × 10 SH array, for Kolmogorov models including 25 orders of Zernike modes, at different aberration strengths D/r₀. Both methods show to be nearly independent of aberration strength.

Download Full Size | PDF

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₀ = 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.

Fig. 6 Performance of the SABRE and SABRE-M for a 10 × 10, 15 × 15 and 20 × 20 SH array. Results are obtained for phase screens including up to 25 orders of Zernike modes with aberration strength D/r₀ = 40.

Download Full Size | PDF

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₀ = 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.

Fig. 7 Performance of the SABRE and SABRE-M on a 10 × 10 SH array at different SNR levels. The dashed line shows the performance for the noiseless case. Results are obtained for phase screens including Zernike modes of polynomial order up to 25 with aberration strength D/r₀ = 40.

Download Full Size | PDF

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]

Nonlinear spline wavefront reconstruction through moment-based Shack-Hartmann sensor measurements

Abstract

1. Introduction

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

2.1. SABRE: a spline-based wavefront reconstruction method

2.2. A second moment technique

3. Wavefront reconstruction with the SABRE-Moments

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

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

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

3.4. The global problem formulation and nonlinear least squares estimator

4. Simulations with the SABRE-M

4.1. Proof of concept

4.2. Power of a higher degree approximation model

4.3. Analysis on the number of subapertures

4.4. Analysis on measurement noise

5. Conclusions

Funding

Acknowledgments

References and links

Cited By

Figures (7)

Equations (26)

Optics Express