Wavefront reconstruction in adaptive optics systems using nonlinear multivariate splines

doi:10.1364/JOSAA.30.000082

« Show journal navigation

Wavefront reconstruction in adaptive optics systems using nonlinear multivariate splines

Cornelis C. de Visser and Michel Verhaegen »View Author Affiliations

JOSA A, Vol. 30, Issue 1, pp. 82-95 (2013)
http://dx.doi.org/10.1364/JOSAA.30.000082

View Full Text Article

Acrobat PDF (1464 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
INTRODUCTION
PRELIMINARIES ON MULTIVARIATE SIMPLEX B-SPLINES
WAVEFRONT RECONSTRUCTION WITH SIMPLEX B-SPLINES
SIMULATIONS WITH THE SABRE
CONCLUSIONS
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (32)
Cited By
Figures (10)
Metrics

Abstract

This paper presents a new method for zonal wavefront reconstruction (WFR) with application to adaptive optics systems. This new method, indicated as Spline based ABerration REconstruction (SABRE), uses bivariate simplex B-spline basis functions to reconstruct the wavefront using local wavefront slope measurements. The SABRE enables WFR on nonrectangular and partly obscured sensor grids and is not subject to the waffle mode. The performance of SABRE is compared to that of the finite difference (FD) method in numerical experiments using data from a simulated Shack–Hartmann lenslet array. The results show that SABRE offers superior reconstruction accuracy and noise rejection capabilities compared to the FD method.

1. INTRODUCTION

Active control of the phase of the photon wavefront for aberration compensation is essential to many photonics applications in science and engineering. It is the field of adaptive optics (AO) that is concerned with measuring and reshaping the wavefront phase in real-time. One of the most important applications of AO is in the field of astronomy, where AO systems are used in optical telescopes to compensate for atmospheric turbulence induced wavefront aberrations that degrade the quality of scientific observations [1

J. M. Beckers, P. Lena, O. Lai, P. Y. Madec, G. Rousset, M. Séchaud, M. J. Northcott, F. Roddier, J. L. Beuzit, F. Rigaut, and D. G. Sandler, Adaptive Optics in Astronomy (Cambridge University, 1999).

Any AO system can be divided into three parts. The first is the wavefront sensing part in which a wavefront sensor (WFS) measures the slopes or curvature of the wavefront. The second part is the wavefront reconstruction (WFR) part that uses the WFS measurements to reconstruct the global wavefront. The third part is the control part that uses the reconstructed wavefront to generate control laws for the actuators of a deformable mirror.

WFR is a complicated process because current WFS’s cannot measure the wavefront directly but rather its slope or its curvature [2

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

J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]

F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. 27, 1223–1225 (1988). [CrossRef]

–5

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

]. This means that the global wavefront must be estimated from WFS measurements, which are often contaminated with sensor noise and biased due to unmodelled parts in the sensor-phase relationship. The WFR estimation problem, together with the large scale of modern WFS arrays, leads to a computationally expensive process that is challenging to implement in real-time control systems. The three most important aspects of any WFR method are its reconstruction accuracy, its computational efficiency, and its ability to deal with real-life sensor and actuator geometries. There exist many different WFR methods, which can roughly be divided into zonal (local) methods and modal (global) methods.

One of the best known and most widely used zonal WFR methods is the finite difference (FD) method, which comes in a number of different forms [2

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

J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]

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

]. The FD method is characterized by the fact that it assumes that the wavefront is defined on a rectangular grid with linear polynomials interpolating between grid points (phase points). Fueled by the development of a new generation of extremely large optical telescopes, innovations in FD methods are mainly focused on increasing their computational efficiency [6

M. Kissler-Patig, “Overall science goals and top level AO requirements for the E-ELT,” presented at First AO4ELT Conference, Victoria, Canada, B.C., September 25 and 30, 2010.

V. Korkiakoski and C. Vérinaud, “Simulations of the extreme adaptive optics system for EPICS,” Proc. SPIE 7736, 773643 (2010).

]. Examples of such innovations are a computationally efficient sparse matrix inversion method by Ellerbroek [8

B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. 19, 1803–1816 (2002). [CrossRef]

] and Vogel [9

C. R. Vogel, “Sparse matrix methods for wavefront reconstruction, revisited,” Proc. SPIE 5490, 1327–1335 (2004).

] and a multigrid preconditioned conjugate-gradient method by Gilles et al. [10

L. Gilles, C. R. Vogel, and B. L. Ellerbroek, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. 19, 1817–1822 (2002). [CrossRef]

] and Vogel and Yang [11

C. R. Vogel and Q. Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt. 45, 705–715 (2006). [CrossRef]

]. More recently, Rosensteiner presented a cumulative reconstruction method based on line integrals that further reduces computational complexity [12

M. Rosensteiner, “Cumulative reconstructor: fast wavefront reconstruction algorithm for extremely large telescopes,” J. Opt. Soc. Am. 28, 2132–2138 (2011). [CrossRef]

Modal methods based on polynomials in polar coordinates, such as the Zernike polynomials and Karhoenen–Loève (KL) functions, are widely used in AO systems [13

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

,14

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

]. A more recent type of modal methods are the Fourier domain methods, which have been developed for the purpose of improving the computational efficiency of the WFR problem [15

L. A. Poyneer, D. T. Gavel, and J. M. Brase, “Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform,” J. Opt. Soc. Am. 19, 2100–2111 (2002). [CrossRef]

,16

L. A. Poyneer, B. A. Machintosh, and J. P. Veran, “Fourier transform wavefront control with adaptive prediction of the atmosphere,” J. Opt. Soc. Am. 24, 2645–2660 (2007). [CrossRef]

]. More recently, Hampton et al. introduced a modal WFR method based on Haar wavelets that was shown to have a linear computational complexity [17

P. J. Hampton, P. Agathoklis, and C. Bradley, “A new wave-front reconstruction method for adaptive optics systems using wavelets,” IEEE J. Select. Topics Signal Process. 2, 781–792 (2008). [CrossRef]

While current WFR methods are well established and proven in practice, they suffer from a lack of generality. FD methods on the one hand are relatively straightforward to implement on rectangular WFS arrays but are essentially linear while the physical wavefront is certainly not linear. Additionally, it is not trivial to implement FD methods on AO systems with nonrectangular WFS arrays and obstructions in the field of view of the telescope pupil (e.g., spiders). Modal methods based on Zernike and KL polynomials have a limited approximation power and are subject to oscillations on the edges of the domain (i.e., Runge’s phenomenon). Additionally, the polar nature of Zernike and KL polynomial methods makes their implementation on noncircular WFS arrays nontrivial [13

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

]. Fourier methods on the other hand have a transparent implementation on rectangular WFS grids. However, on nonrectangular, nonhomogeneous, or partially obscured WFS grids, their implementation becomes more complex [15

The objective of this paper is to present a new method for WFR from wavefront slope measurements that uses bivariate simplex B-splines (Fig. 1) inside a linear regression framework [18

G. Awanou, M. J. Lai, and P. Wenston, “The multivariate spline method for scattered data fitting and numerical solutions of partial differential equations,” in Wavelets and Splines , G. Chen and M. J. Lai, eds. (Nashboro, 2005), pp. 24–75.

C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).

C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]

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]

–22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

]. This new method, which we call Spline based ABerration REconstruction (SABRE), aims to be a truly general WFR method. In essence, SABRE locally models the wavefront with linear and nonlinear simplex B-spline basis functions on triangular subpartitions of the WFS domain using local WFS measurements. While this paper presents only a least squares (LS) estimator for the simplex spline coefficients (i.e., B-coefficients), it is compatible with any more advanced linear regression parameter estimation technique.

SABRE has five important advantages over other WFR methods. First, SABRE is invariant of WFS geometry in the sense that it can be used without any modification on nonrectangular, nonhomogeneous (misaligned), and partially obscured sensor grids. This is a significant advantage in real life AO setups, where nonrectangular WFS grids with misaligned lenslet images are often encountered. Second, SABRE allows WFR using nonlinear basis functions, resulting in a more accurate approximation of the physical wavefront. Third, SABRE has an inherent noise smoothing capability that makes it highly resilient to sensor noise. Fourth, in contrast to the Fried-geometry based FD method, the SABRE is not subject to the waffle mode [23

M. D. Oliker, “Sensing waffle in the Fried geometry,” Proc. SPIE 3353, 964–971 (1998). [CrossRef]

,24

W. Zou and J. P. Rolland, “Quantifications of error propagation in slope-based wavefront estimations,” J. Opt. Soc. Am. 23, 2629–2638 (2006). [CrossRef]

]. Finally, the local nature of SABRE means that it can be implemented on multicore hardware resulting in a distributed SABRE (D-SABRE) that can significantly increase computational efficiency. This paper will focus primarily on the first four of these advantages, while the D-SABRE will be explored in a forthcoming paper.

This paper is outlined as follows. In Section 2 we provide preliminaries on bivariate simplex B-splines as they are central to the new WFR method. Then, in Section 3 we introduce the SABRE method and present an LS estimator for estimating the SABRE model parameters. Additionally, we show in Section 3 that for fundamental reasons the SABRE is not subject to the waffle mode. In Section 4 the results from a number of numerical experiments utilizing a Fourier optics based Shack–Hartmann (SH) sensor simulation are presented. In the experiments it is shown that SABRE can reconstruct nonlinear wavefronts that are closer to physical reality than wavefronts produced by any FD method. Subsequently, the ability of SABRE to reconstruct wavefronts on nonrectangular domains is demonstrated. In Section 5 we conclude the paper and provide pointers for future research.

2. PRELIMINARIES ON MULTIVARIATE SIMPLEX B-SPLINES

To provide the reader with a basic understanding of the theory behind SABRE a brief introduction into the theory of bivariate simplex B-splines is given. For a more complete and general introduction into the theory of multivariate simplex B-splines we refer to [22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

A. Two-Simplex and Barycentric Coordinates

Let

t

be a two-simplex (triangle) formed by the convex hull of its three nondegenerate vertices

(v_{0}, v_{1}, v_{2}) \in R^{2}

as follows:

t ≔ 〈 [\begin{array}{c} v_{0_{x}} \\ v_{0_{y}} \end{array}], [\begin{array}{c} v_{1_{x}} \\ v_{1_{y}} \end{array}], [\begin{array}{c} v_{2_{x}} \\ v_{2_{y}} \end{array}] 〉 \in R^{2} .

(1)

The basis polynomials of the simplex B-splines are functions in terms of Barycentric coordinates. The Barycentric coordinate system is a local coordinate system that is defined on an individual simplex. If

x \in R^{2}

is a point on the Cartesian plane, then the normalized Barycentric coordinate

b \in R^{3}

x

with respect to the triangle

t

can be determined using the following equations:

[\begin{array}{c} x_{1} \\ x_{2} \end{array}] = [\begin{array}{ccc} v_{0_{x}} & v_{1_{x}} & v_{2_{x}} \\ v_{0_{y}} & v_{1_{y}} & v_{2_{y}} \end{array}] [\begin{array}{c} b_{0} \\ b_{1} \\ b_{2} \end{array}], b_{0} + b_{1} + b_{2} = 1 .

(2)

The condition

b_{0} + b_{1} + b_{2} = 1

is a normalization condition, which ensures that any

x \in R^{2}

has a unique representation

b \in R^{3}

in Barycentric coordinate space. Now define the normalized simplex vertex matrix

V

as follows:

V ≔ [\begin{array}{cc} (v_{1_{x}} - v_{0_{x}}) & (v_{2_{x}} - v_{0_{x}}) \\ (v_{1_{y}} - v_{0_{y}}) & (v_{2_{y}} - v_{0_{y}}) \end{array}] .

(3)

Using the matrix

V

from (3), the Barycentric coordinate

b = (b_{0}, b_{1}, b_{2}) \in R^{3}

of the Cartesian coordinate

x \in R^{2}

is calculated as follows:

[\begin{matrix} b_{1} \\ b_{2} \end{matrix}] = V^{- 1} [\begin{matrix} x_{1} \\ x_{2} \end{matrix}] b_{0} = 1 - b_{1} - b_{2} .

(4)

In the remainder of the paper, we shall use the following shorthand notation for the Cartesian-to-Barycentric coordinate transformation from

x \in R^{2}

b \in R^{3}

on the triangle

t

b (x) ≔ (b_{0}, b_{1}, b_{2}) \in R^{3}, x \in R^{2} .

(5)

B. Triangulations of Simplices

The approximation power of the multivariate simplex B-spline can be increased by combining many simplices into a structure called a triangulation. A triangulation

T

is a special partitioning of a domain into a set of

J

nonoverlapping simplices:

T ≔ ⋃_{i = 1}^{J} t_{i}, t_{i} \cap t_{j} \in {\emptyset, \tilde{t}}, \forall t_{i}, t_{j} \in T

(6)

with the edge simplex

\tilde{t}

either a line (1-simplex), or a vertex (0-simplex) in the case of a two-dimensional triangulation consisting of triangles.

A number of algorithms can be used to create triangulations from a given set of vertices. The most widely known of these is the Delaunay triangulation method, which has a standard implementation in Matlab. The triangulations used in this paper were all created with a different (non-Delaunay) technique based on the grid cell subdivision scheme introduced in [19

C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).

C. Basis Functions of the Simplex B-Splines

The basis polynomials of the simplex B-splines are Bernstein polynomials in terms of Barycentric coordinates. The basis polynomials are derived using the multinomial theorem, which states that any polynomial of total degree

d

can be expanded into a sum of monomials. In

R^{3}

the result of the multinomial theorem is the following:

{(b_{0} + b_{1} + b_{2})}^{d} = \sum_{κ_{0} + κ_{1} + κ_{2} = d} \frac{d!}{κ_{0}! κ_{1}! κ_{2}!} b_{0}^{κ_{0}} b_{1}^{κ_{1}} b_{2}^{κ_{2}}

(7)

with

κ = (κ_{0}, κ_{1}, κ_{2})

a multi-index with the properties

| κ | = κ_{0} + κ_{1} + κ_{2} = d, κ_{0} \geq 0, κ_{1} \geq 0, κ_{2} \geq 0 .

(8)

The Bernstein basis polynomials of the simplex B-splines are defined as the individual monomials in (7), with the additional rule that they are equal to 0 by definition when the evaluation point

x

is outside of the triangle

t

B_{κ}^{d} (b (x)) ≔ {\begin{cases} \frac{d!}{κ_{0}! κ_{1}! κ_{2}!} b_{0}^{κ_{0}} b_{1}^{κ_{1}} b_{2}^{κ_{2}}, & x \in t \\ 0, & x \notin t \end{cases} .

(9)

Any polynomial

p (b (x))

of degree

d

on a simplex

t

can be written as a linear combination of basis polynomials in what is known as the B-form as follows [25

C. de Boor, “B-form basics,” in Geometric Modeling: Algorithms and New Trends (SIAM, 1987).

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

(10)

with

c_{κ}^{t}

the B-coefficients that uniquely determine the polynomial

p (b (x))

on the triangle

t

. The B-coefficients have a special geometric ordering inside their parent simplex [see Fig. 2]. This ordering is called the B-net and is essential for defining continuity between simplices and for enforcing local external constraints on the simplex B-spline function [20

C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]

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]

–22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

Fig. 1. Principle of the multivariate simplex spline: A 5th degree spline function with

C^{1}

continuity defined on four triangles with (left) the four individual spline pieces [

p_{1} (b)

p_{2} (b)

p_{3} (b)

, and

p_{4} (b)

] and (right) the global spline function

p (b)

formed by combining the four spline pieces.

Fig. 2. B-net for a 4th degree spline function on a triangulation consisting of the three triangles

t_{i}

t_{j}

, and

t_{k}

The total number of B-coefficients and basis polynomials is equal to

\hat{d}

, which for the two-dimensional case and a given degree

d

is given by

\hat{d} ≔ \frac{(d + 2)!}{2 d!} .

(11)

D. Vector Formulations of the B-Form

In [20

C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]

] a vector formulation of the B-form from (10) was introduced. With (5) the vector formulation for a B-form polynomial

p (b (x))

in Barycentric

R^{3}

p (b (x)) ≔ {\begin{cases} B^{d} (b (x)) \cdot c^{t}, & x \in t \\ 0, & x \notin t \end{cases}

(12)

with

b (x)

the Barycentric coordinate of the Cartesian

x

according to (5).

The row vector

B^{d} (b (x))

in (12) is a vector that is constructed from individual basis polynomials that are sorted lexicographically according to [26

X. L. Hu, D. F. Han, and M. J. Lai, “Bivariate splines of various degrees for numerical solution of partial differential equations,” SIAM J. Sci. Comput. 29, 1338–1354 (2007). [CrossRef]

B^{d} (b (x)) ≔ [\begin{matrix} B_{d, 0, 0}^{d} (b (x)) & B_{d - 1, 1, 0}^{d} (b (x)) & \dots \end{matrix} \begin{matrix} \dots & B_{0, 1, d - 1}^{d} (b (x)) & B_{0, 0, d}^{d} (b (x)) \end{matrix}] \in R^{1 \times \hat{d}} .

(13)

The column vector

c^{t}

the vector of lexicographically sorted B-coefficients on the triangle

t

c^{t} = {[\begin{matrix} c_{d, 0, 0} & c_{d - 1, 1, 0} & \dots & c_{0, 1, d - 1} & c_{0, 0, d} \end{matrix}]}^{⊤} \in R^{\hat{d} \times 1} .

(14)

For example, for a B-form polynomial

p (b (x))

of degree

d = | κ | = 1

in Barycentric

R^{3}

on the triangle

t

we have

κ \in {(1, 0, 0), (0, 1, 0), (0, 0, 1)}

. In this case the vector formulation of the B-form from (12) is

p (b (x)) = B^{1} (b (x)) \cdot c^{t} = [\begin{matrix} b_{0}^{1} b_{1}^{0} b_{2}^{0} & b_{0}^{0} b_{1}^{1} b_{2}^{0} & b_{0}^{0} b_{1}^{0} b_{2}^{1} \end{matrix}] {[\begin{matrix} c_{1, 0, 0}^{t} & c_{0, 1, 0}^{t} & c_{0, 0, 1}^{t} \end{matrix}]}^{⊤} .

The simplex B-spline function

s_{d}^{r} (b (x))

of degree

d

and continuity order

r

, defined on a triangulation

T_{J}

consisting of

J

triangles, is defined as follows:

s_{d}^{r} (b (x)) ≔ B^{d} \cdot c \in R, x \in T_{J},

(15)

where the continuity order

r

, also denoted by

C^{r}

, signifies that all derivatives up to order

r

of two B-form polynomials defined on two neighboring triangles are equal on the edge between the two triangles. For example,

C^{0}

continuity means that only the values of the B-form polynomials are equal on an edge between two neighboring triangles, while

C^{1}

continuity means that both the first order derivatives and the values of the B-form polynomials match on the edge.

In (15),

B^{d}

is the global vector of basis polynomials that is constructed using (13) as follows:

B^{d} ≔ [\begin{array}{cccc} B_{t_{1}}^{d} (b (x)) & B_{t_{2}}^{d} (b (x)) & \dots & B_{t_{J}}^{d} (b (x)) \end{array}] \in R^{1 \times J \cdot \hat{d}} .

(16)

Note that according to (12) we have

B_{t_{j}}^{d} (b (x)) = 0

for all evaluation locations

x

that are located outside of the triangle

t_{j}

. The result of this is that the global vector of basis polynomials

B^{d}

is a sparse vector.

The global vector of B-coefficients

c

in (15) is constructed as follows:

c ≔ {[\begin{array}{cccc} {c^{t_{1}}}^{⊤} & {c^{t_{2}}}^{⊤} & \dots & {c^{t_{J}}}^{⊤} \end{array}]}^{⊤} \in R^{J \cdot \hat{d} \times 1}

(17)

with each

c^{t_{j}}

a per-simplex vector of lexicographically sorted B-coefficients from (14).

E. Spline Spaces

A spline space is defined as the space of all spline functions

s_{d}^{r}

of a given degree

d

and continuity order

r

on a given triangulation

T

. Such spline spaces have been well studied in the past (see, e.g., [22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

,27

M. J. Lai and L. L. Schumaker, “On the approximation power of bivariate splines,” Adv. Comput. Math. 9, 251–279(1998). [CrossRef]

,28

M. J. Lai, “Some sufficient conditions for convexity of multivariate Bernstein–Bezier polynomials and box spline surfaces,” Studia Scient. Math. Hung. 28, 363–374 (1990). [CrossRef]

]). We use the definition of the spline space from [22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

S_{d}^{r} (T) ≔ {s_{d}^{r} \in C^{r} (T) : {s_{d}^{r} |}_{t} \in P_{d}, \forall t \in T}

(18)

with

P_{d}

the space of polynomials of degree

d

F. Continuity between Simplices

By definition, a spline function is a piecewise defined polynomial function with a predefined continuity order between its polynomial pieces. The continuity order

r

signifies that all

m

th order derivatives, with

0 \leq m \leq r

, of two B-form polynomials defined on two neighboring triangles are equal on the edge between the two triangles. For simplex B-splines, continuity between neighboring triangles is enforced by continuity conditions.

The principle of the continuity conditions in Cartesian

R^{2}

is the following. Let two neighboring triangles

t_{i}

and

t_{j}

, differing by only the vertex

w

, be defined as follows:

t_{i} ≔ 〈 v_{0}, v_{1}, w 〉, t_{j} ≔ 〈 v_{0}, v_{1}, v_{2} 〉 .

(19)

Then

t_{i}

and

t_{j}

meet along the line

\tilde{t}

given by

\tilde{t} ≔ t_{i} \cap t_{j} = 〈 v_{0}, v_{1} 〉 .

(20)

The formulation for the continuity conditions for

C^{r}

continuity between

t_{i}

and

t_{j}

is the following [18

,22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

- c_{(κ_{0}, κ_{1}, m)}^{t_{i}} + \sum_{| γ | = m} c_{(κ_{0}, κ_{1}, 0) + γ}^{t_{j}} B_{γ}^{m} (b (w)) = 0, 0 \leq m \leq r,

(21)

with

γ = (γ_{0}, γ_{1}, γ_{2})

a multi-index independent of

κ

It was shown in [20

C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]

] that the formulation provided in (21) is only accurate for specific B-net orientations, such as that of

t_{i}

and

t_{j}

in Fig. 2. A more general formulation of the continuity conditions requires the utilization of a B-net orientation rule such as that introduced in [19

C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).

For

C^{r}

continuity there are a total of

Q

continuity conditions of the form (21) per edge:

Q ≔ \sum_{m = 0}^{r} (d - m + 1) .

(22)

The continuity conditions for all

E

edges are collected into a single set of linear equations:

H c = 0 .

(23)

In (23), the matrix

H \in R^{E Q \times J \hat{d}}

is the so-called global smoothness matrix. Each row in

H

contains a single continuity condition. The vector

c

is the global vector of B-coefficients from (17). Constructing

H

is not a trivial task; for particular details on its construction we refer the reader to [19

C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).

,20

C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]

,22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

G. Matrix Form of the Directional Derivative

In [21

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]

] a formulation of the directional derivatives of B-form polynomials in matrix form was presented. Because of its crucial role in WFR with splines, this formulation will be repeated here to increase the understanding of the reader.

Before introducing the formulation, however, the concept of the directional coordinate must be explained. First, let

u

be a unit vector in Cartesian

R^{2}

; then we define

a \in R^{3}

as the directional coordinate of

u

with respect to the triangle

t

. The directional coordinate

a

should be seen as the Barycentric representation of the Cartesian unit vector

u \in R^{2}

with respect to a given triangle. That is, if

u = v - w

, with

v

and

w

vectors in Cartesian

R^{2}

, then using the short-hand notation from (5), the directional coordinate is defined as follows:

a ≔ b (v) - b (w) \in R^{3}

(24)

with

b (v)

and

b (w)

the Barycentric coordinates with respect to

t

v

and

w

, respectively.

Using the formulation from [21

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]

], the directional derivative of order

m

in the direction

u

of a B-form polynomial

p (b (x))

on a single triangle

t

can be expressed in terms of the original vector of B-coefficients as follows:

D_{u}^{m} p (b (x)) = \frac{d!}{(d - m)!} B^{d - m} (b (x)) P^{d, d - m} (a) \cdot c^{t}

(25)

with

P^{d, d - m} (a) \in R^{((d - m + 2)! / 2 (d - m)!) \times \hat{d}}

the de Casteljau matrix of degree

d

d - m

from [21

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]

] expressed in terms of the directional coordinate

a

u

with respect to the triangle

t

. In (25),

B^{d - m} (b (x))

is the vector of basis polynomials from (13) for degree

d - m

, and

c^{t}

is the vector of B-coefficients for a single triangle from (14).

It will be shown in the next section that (25) is the key to spline based WFR, because it directly relates the spline function model of the wavefront to its directional derivatives modeling the slopes and curvatures in terms of

c^{t}

For example, in the case of the first order directional derivative of a simplex B-spline

p (b (x))

of degree 1 on a single triangle, (25) can be simplified as follows:

D_{u}^{1} p (b (x)) = B^{0} (b (x)) P^{1, 0} (a) \cdot c^{t},

(26)

which for

b

a \in R^{3}

reduces to

D_{u}^{1} p (b (x)) = {\begin{cases} [\begin{array}{ccc} a_{0} & a_{1} & a_{2} \end{array}] \cdot c^{t}, & x \in t, \\ 0, & x \notin t \end{cases} .

(27)

Before continuing, a full-triangulation formulation of de Casteljau matrix for a triangulation

T_{J}

of the form (6) consisting of

J

triangles needs to be defined. For a single observation, this full-triangulation formulation is a block diagonal matrix with a total of

J

blocks of the form

P^{d, d - m} (a_{u})

on the main diagonal, with

j = 1, 2, \dots, J

P_{u}^{d, d - m} ≔ diag {(P_{j}^{d, d - m} (a_{u}))}_{j = 1}^{J} \in R^{J \frac{(d - m + 2)!}{2 (d - m)!} \times J \hat{d}},

(28)

where

a_{u}

is the directional coordinate of the derivative direction

u

with respect to the triangle

t_{j}

The full triangulation form of the directional derivative of the simplex B-spline for a single observation

b (x)

then follows from the combination of (16), (28), and (17):

D_{u}^{m} s_{d}^{r} (b (x)) = \frac{d!}{(d - m)!} B^{d - m} P_{u}^{d, d - m} \cdot c .

(29)

3. WAVEFRONT RECONSTRUCTION WITH SIMPLEX B-SPLINES

In this section the SABRE method for WFR is introduced and connected to the literature. Additionally, we show that the SABRE is not subject to the waffle mode. Finally, the steps that make up the SABRE algorithm are presented, its computational aspects are discussed, and a tutorial example is given.

A. Wavefront Reconstruction from Slope Measurements

The relationship between the slopes of the wavefront phase and the wavefront phase can be described in the form of the following system of first order partial differential equations (PDEs) [3

J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]

σ_{x} (x, y) = \frac{\partial ϕ (x, y)}{\partial x}, σ_{y} (x, y) = \frac{\partial ϕ (x, y)}{\partial y}

(30)

with

ϕ (x, y)

the unknown wavefront, and with

σ_{x} (x, y)

and

σ_{y} (x, y)

the wavefront slopes at

(x, y)

in the directions

x

and

y

, respectively.

B. Wavefront Reconstruction with the Finite Difference Method

One of the best known zonal methods for the reconstruction of wavefronts from wavefront slopes is the FD method. The FD method comes in a number of different specific forms that depend on the sampling geometry, such as the Fried [2

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

], Hudgin [29

R. H. Hudgin, “Wavefront reconstruction for compensated imaging,” J. Opt. Soc. Am. 67, 375–378(1977). [CrossRef]

], and Southwell [5

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

] geometries. All FD methods have in common that they reduce (30) into a discrete form whose solution is approximated with linear interpolating polynomials [3

J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]

In order to connect the SABRE method with the literature we shall compare it directly with the Fried form of the FD method [2

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

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

]. The Fried FD sensor model assumes that the average slope in a grid cell formed by four phase points is equal to the average difference in phase between the phase points [see Fig. 3]. The Fried sensor model is

σ_{x} (i, j) = [(ϕ (i + 1, j) - ϕ (i, j)) + (ϕ (i + 1, j + 1) - ϕ (i, j + 1))] / (2 h) + n_{x} (i, j) σ_{y} (i, j) = [(ϕ (i, j + 1) - ϕ (i, j)) + (ϕ (i + 1, j + 1) - ϕ (i + 1, j))] / (2 h) + n_{y} (i, j)

(31)

with

h

the step size between phase points, with

n_{x} (i, j)

and

n_{y} (i, j)

residual terms that contain both noise and modeling errors. The indices

i = 1, 2, \dots, M - 1

, and

j = 1, 2, \dots, N - 1

are grid coordinates on the grid domain

Ω (M, N)

, which is defined as follows:

Ω (M, N) \in R^{2}, M, N \in N .

(32)

Fig. 3. Southwell geometry (top left) and Fried (top right) sensor geometries compared with four different SABRE geometries (middle and bottom rows). Black dots are the phase point locations and horizontal and vertical lines are the slope measurements in the

x

and

y

directions, respectively. The open circles are the locations of the slope measurements. Gray lines in the SABRE geometries are the triangle edges, while the shaded area inside the triangulations is the area in which

ϕ (x, y)

is defined.

By enumerating all

M \cdot N

phase points

ϕ (i, j)

into a single row vector

ϕ

, and by using all slope measurements, the following matrix equation can be constructed:

σ = G ϕ + n

(33)

with

σ = {[\begin{matrix} σ_{x}^{⊤} & σ_{y}^{⊤} \end{matrix}]}^{⊤}

the vector of slope measurements that is constructed by appending the vector of all slopes in the

x

direction with the vector of all slopes in the

y

direction. The matrix

G

is the geometry matrix that is constructed such that each row multiplied with

ϕ

results in a single expression of the form (31).

Based on (33), an estimate for the global wavefront

\hat{ϕ}

can now be obtained using any parameter estimation technique (see, e.g., [3

J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]

B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. 19, 1803–1816 (2002). [CrossRef]

C. R. Vogel, “Sparse matrix methods for wavefront reconstruction, revisited,” Proc. SPIE 5490, 1327–1335 (2004).

,11

C. R. Vogel and Q. Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt. 45, 705–715 (2006). [CrossRef]

]). For example, the ordinary LS estimator for

ϕ

{\hat{ϕ}}_{FD} = G^{+} σ

(34)

with

G^{+}

the Moore–Penrose pseudo inverse of

G

G^{⊤} G

is singular in general [1

]. In particular, for the Southwell FD geometry, the nullspace of

G^{⊤} G

has dimension 1 (piston mode), while for the Fried FD geometry, the nullspace has dimension 2 (piston mode and waffle mode) [23

M. D. Oliker, “Sensing waffle in the Fried geometry,” Proc. SPIE 3353, 964–971 (1998). [CrossRef]

,24

W. Zou and J. P. Rolland, “Quantifications of error propagation in slope-based wavefront estimations,” J. Opt. Soc. Am. 23, 2629–2638 (2006). [CrossRef]

]. The matrix

G^{+}

is the FD reconstruction matrix, which is computed only once for a given FD sensor geometry.

C. Wavefront Reconstruction with the SABRE

The main contribution of this paper is the SABRE, which is a new method for WFR that uses the basis functions of bivariate simplex B-splines to locally model the wavefront using local WFS measurements. The SABRE can be seen as a generalization of the FD method to the nonlinear case [see Appendix].

The approach taken with the SABRE method differs from the FD method in the sense that it is assumed that the unknown wavefront

ϕ (x, y)

can be approximated with a parametrized model in the form of a bivariate simplex B-spline of degree

d \geq 1

and continuity order

r

from (15):

ϕ (x, y) \approx s_{d}^{r} (b (x, y)) = B^{d} (b (x, y)) \cdot c, d \geq 1, (x, y) \in t .

(35)

The right-hand side of (35) then is the SABRE model that approximates the wavefront

ϕ (x, y)

Under the assumption that (35) holds, we can then use (25) to rewrite (30) into the following slope sensor model:

σ_{x} (x, y) = \frac{d!}{(d - 1)!} B^{d - 1} (b (x, y)) P^{d, d - 1} (a_{x}) \cdot c^{t} + ν_{x} (x, y), σ_{y} (x, y) = \frac{d!}{(d - 1)!} B^{d - 1} (b (x, y)) P^{d, d - 1} (a_{y}) \cdot c^{t} + ν_{y} (x, y)

(36)

with

d \geq 1

the polynomial degree of the spline, and with

a_{x}

and

a_{y}

the directional coordinates respectively of

u_{x}

and

u_{y}

with respect to the triangle

t

. The terms

ν_{x} (x, y)

and

ν_{y} (x, y)

in (36) are residual terms which contain both the sensor noise as well as the modeling errors.

Use of the SABRE method requires the introduction of a new sensor geometry that is triangular in nature. In Fig. 3 the basic SABRE sensor geometry is compared to that of the Southwell and Fried FD sensor geometries. In Table 1 the properties of the different SABRE geometries are shown. The fundamental difference between the FD and the SABRE methods is that for all FD methods the wavefront phase is defined only at the grid locations

(m, n) \in N^{2}

, while for the SABRE method the phase is defined at any location

(x, y) \in T

. The result is that the SABRE method allows for the decoupling between slope measurement locations and phase point locations.

Table 1. Properties of SABRE Geometries

Sensor Geometry	Triangulation Type	Rectangular Domain	Data at Vertices
Type-I	I	Yes	Yes
Type-IG	I	Yes	No
Type-IB	I	No	Yes
Type-IBG	I	No	No
Type-II	II	Yes	Yes
Type-IIG	II	Yes	No
Type-IIB	II	No	Yes
Type-IIBG	II	No	No

D. Anchor Constraint

Before continuing with the definition of the SABRE, a new type of constraint is introduced in the form of the anchor constraint. Effectively, the anchor constraint predefines the value of the unknown constant of integration (i.e., the piston mode) that arises when solving the first order partial differential equation (PDE) from (30). The anchor constraint is essential for producing a well-conditioned parameter estimation problem for the B-coefficients of the SABRE model.

The anchor constraint is derived as follows. First, we reformulate (29) with

m = 1

in integral form as follows:

\int_{T} D_{u}^{1} s_{d}^{r} (b (x)) d u = \int_{T} d B^{d - 1} P^{d, d - m} c d u, = B^{d} c + k

(37)

with

k

an unknown constant that is proportional to the piston mode. The resulting spline model for the wavefront then is

s_{d}^{r} (b (x)) = B^{d} c + k . 1

(38)

The affine property of the B-coefficients allows us to rewrite (38) as follows:

s_{d}^{r} (b (x)) = B^{d} (c + k \cdot 1)

(39)

with

1 \in R^{J \cdot \hat{d} \times 1}

a row vector containing only 1’s.

The global B-coefficient vector can be partitioned into two parts, the first part containing only the first B-coefficient, and the second part the remainder of B-coefficients:

c = [\begin{matrix} c_{d, 0, 0}^{t_{1}} + k \\ \tilde{c} + k \cdot 1 \end{matrix}]

(40)

with

1 \in R^{(J \cdot \hat{d} - 1) \times 1}

The anchor constraint then is the following:

k = - c_{d, 0, 0}^{t_{1}} .

(41)

Substitution of (41) in (39) then results in the SABRE model with fixed piston mode:

s_{d}^{r} (b (x)) = B^{d} \cdot [\begin{matrix} 0 \\ \tilde{c} - c_{d, 0, 0}^{t_{1}} \cdot 1 \end{matrix}] .

(42)

In the following, we require a vector form of the anchor constraint. If

h = [\begin{matrix} 1 & 0 & \dots & 0 \end{matrix}]

is the anchor vector, the anchor constraint becomes:

h \cdot {[\begin{matrix} (c_{d, 0, 0}^{t_{1}} + k) & {(\tilde{c} + k \cdot 1)}^{⊤} \end{matrix}]}^{⊤} = 0 .

(43)

E. Linear SBS

The linear SABRE problem for a total of

K

slope measurements on a complete triangulation consisting of

J

triangles is derived from (36) as follows:

σ = B^{0} P_{u}^{1, 0} c + n, 0 = A c

(44)

with

σ \in R^{K \times 1}

the vector of measured wavefront slopes,

B^{0} \in R^{K \times J}

the global matrix of basis polynomials from (16) of polynomial degree 0,

P_{u}^{1, 0} \in R^{J \times 3 J}

the de Casteljau matrix from (28), and

c \in R^{3 J \times 1}

the global vector of B-coefficients from (17). The matrix

A

in (44) is a constraint matrix defined as follows:

A ≔ [\begin{array}{c} H \\ h \end{array}] \in R^{(E V + 1) \times J \hat{d}}

(45)

with

H \in R^{(E V) \times J \hat{d}}

the smoothness matrix from (23) and

h \in R^{1 \times J \hat{d}}

the anchor constraint from (43).

For a single triangle

t

, the slope equation from (44) can be simplified using (27):

σ_{x} (x, y) = a_{x_{0}} c_{1, 0, 0} + a_{x_{1}} c_{0, 1, 0} + a_{x_{2}} c_{0, 0, 1} + n_{x} (x, y), σ_{y} (x, y) = a_{y_{0}} c_{1, 0, 0} + a_{y_{1}} c_{0, 1, 0} + a_{y_{2}} c_{0, 0, 1} + n_{y} (x, y)

(46)

with

(a_{x_{0}}, a_{x_{1}}, a_{x_{2}})

the directional coordinates of a vector

u_{x} = (1, 0)

in the

x

direction and

(a_{y_{0}}, a_{y_{1}}, a_{y_{2}})

the directional coordinates of a vector

u_{y} = (0, 1)

in the

y

direction, all with respect to

t

. Note that the SABRE model does not depend on the location

(x, y)

; the wavefront slope is assumed to be constant over the complete triangle.

F. Nonlinear SABRE

The nonlinear SABRE problem follows from extending the per-triangle nonlinear sensor model from (36) into the following full triangulation form for all

K

slope and

L

curvature measurements:

σ = d B^{d - 1} P_{u}^{d, d - 1} c + n, 0 = A c

(47)

with

σ \in R^{K \times 1}

the vector of measured wavefront slopes, and with

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

the global vector of B-coefficients from (17). The matrix of basis functions in (47) is constructed according to (16), while the de Casteljau matrix is constructed using (28). Finally, the constraint matrix

A

in (47) is constructed using (45).

An important implication of nonlinear SABRE is that the SABRE model is a function of the geometric location of the slope measurements [i.e.,

(x, y)

]. The reason for this is that the matrix with B-form regressors

B^{d - 1}

in (47) for

d \geq 2

depends on the (Barycentric) location of the slope measurement

b (x, y)

For example, if we assume a quadratic polynomial model structure for the wavefront, that is,

ϕ (x, y) \approx B^{2} c

, the slope equation in (47) becomes

σ = 2 B^{1} P_{u}^{2, 1} c + n .

(48)

For a single triangle

t

, (48) can be expanded as follows:

σ_{x} (x, y) = 2 b_{0} (a_{x_{0}} c_{2, 0, 0} + a_{x_{1}} c_{1, 1, 0} + a_{x_{2}} c_{1, 0, 1}) + 2 b_{1} (a_{x_{0}} c_{1, 1, 0} + a_{x_{1}} c_{0, 2, 0} + a_{x_{2}} c_{0, 1, 1}) + 2 b_{2} (a_{x_{0}} c_{1, 0, 1} + a_{x_{1}} c_{0, 1, 1} + a_{x_{2}} c_{0, 0, 2}) + n_{x} (x, y), σ_{y} (x, y) = 2 b_{0} (a_{y_{0}} c_{2, 0, 0} + a_{y_{1}} c_{1, 1, 0} + a_{y_{2}} c_{1, 0, 1}) + 2 b_{1} (a_{y_{0}} c_{1, 1, 0} + a_{y_{1}} c_{0, 2, 0} + a_{y_{2}} c_{0, 1, 1}) + 2 b_{2} (a_{y_{0}} c_{1, 0, 1} + a_{y_{1}} c_{0, 1, 1} + a_{y_{2}} c_{0, 0, 2}) + n_{y} (x, y)

(49)

with

b (x, y) = (b_{0}, b_{1}, b_{2})

the Barycentric coordinate of the Cartesian point

(x, y)

with respect to

t

. Clearly, (49) is dependent on the actual location of the slope measurement inside the triangle.

G. Generalized SABRE

As a final generalization, the SABRE method enables the use of the resultant directional derivative instead of the set of directional derivatives in the

x

and

y

directions [see Fig. 4]. In this case, the two-component PDE from (36) can be replaced with a single PDE defined in the general direction

u

as follows:

σ_{u} (x, y) = \frac{d!}{(d - 1)!} B^{d - 1} (b (x, y)) P^{d, d - 1} (a_{u}) \cdot c^{t} + n (x, y)

(50)

with

u

a Cartesian vector in

R^{2}

and with

a_{u}

the corresponding directional coordinate in Barycentric

R^{3}

. The advantage of (50) is that compared to the two-component PDE from (36) it requires only half the number of slope measurements. This translates into lower computational and memory loads during the construction of the reconstruction matrix. The disadvantage of (50), however, is that

u

will vary for time-varying wavefronts, requiring a complete rebuild of the reconstruction matrix at each time step.

Fig. 4. General simplex Type-IB geometry (left) and simplex Type-IG geometry (right). Open circles are the lenslet apertures.

H. Least Squares Estimator for the B-Coefficients

In this paragraph an LS estimator is presented for the B-coefficients of the SABRE model approximating a wavefront using wavefront slope measurements.

The LS estimator requires that the constraint equations in (44) and (47) are eliminated from the system by projection on the nullspace of the constraint matrix

A

from (45). In that case, the WFR problem can be formulated in the form of an unconstrained linear regression problem as follows:

σ = D c + n

(51)

with

n

a vector containing Gaussian distributed white noise and with

D

the system matrix given by

D = d B^{d - 1} P_{u}^{d, d - 1} N_{A}

(52)

with

N_{A}

a basis for

null (A)

. The matrix

D

in (51) is the SABRE equivalent of the FD geometry matrix

G

from (33).

The OLS cost function in terms of the B-coefficients is

J (c) = {(σ - D c)}^{⊤} (σ - D c) .

(53)

The LS estimator for the B-coefficients of the SABRE model is

{\hat{c}}_{LS} = N_{A} {(D^{⊤} D)}^{- 1} D^{⊤} σ, = Q σ

(54)

with

Q = N_{A} {(D^{⊤} D)}^{- 1} D^{⊤}

the SABRE reconstruction matrix, which is computed only once for a given geometry. Matrix

D^{⊤} D

is guaranteed to be of full rank and therefore has an empty null-space. This is in contrast with the FD equivalent

G^{⊤} G

from (33), which in general has a rank deficiency of two and consequently a two-dimensional null-space. As a result the SABRE has a predefined piston mode through the anchor constraint, and is not subject to the waffle mode through the action of the continuity constraints.

Note that the matrix

D^{⊤} D

is invertible, which is a direct result of the use of the anchor constraint from (43) together with the smoothness constraints from (23).

Using (54) the resulting LS SABRE model is

{\hat{ϕ}}_{LS} (x) = B^{d} ({b x \hat{) c}}_{LS} .

(55)

I. SABRE Algorithm

The SABRE algorithm consists of two parts: an initialization part and a real-time part. The SABRE algorithm differs from the FD reconstruction algorithm only during the initialization part in which the SABRE reconstruction matrix is determined. The initialization part of the SABRE algorithm consists of the following four steps:

Input: The inputs to the algorithm are the SABRE model structure in the form of the triangulation type (see Table 1), the spline degree $d$ , and the continuity order $r$ . Additionally, a reference WFS image is supplied, and the user specifies a parameter estimator.
Step 1: The reference centers of the SH lenslet array are determined using the reference WFS image. A triangulation $T_{J}$ of the type specified by the user is created by triangulating the reference centers using a (Delaunay) triangulation algorithm.
Step 2: The system of equations from (47) are formulated. For this, the matrix with B-form regressors $B^{d - 1}$ from (28), and the global de Casteljau matrix $P_{u}^{d, d - 1}$ from (28) are constructed. The constraint matrix $A$ from (45) is assembled using the smoothness matrix $H$ from (23) and the anchor constraint $h$ from (43). A basis for the null space of $A$ is calculated.
Step 3: The system matrix $D$ is constructed using (52), and the problem is stated in the linear regression form from (51).
Step 4: The SABRE reconstruction matrix $Q$ is calculated based on the choice of parameter estimator. In the case of an LS estimator, $Q$ is constructed according to (54).

The real-time part of the SABRE algorithm uses the SABRE reconstruction matrix

Q

in an matrix-vector-multiplication (MVM) operation to reconstruct the wavefront phase from the SH slope measurements.

J. Computational Aspects of the SABRE

An important aspect of any WFR method is its computational complexity. While it is not the main emphasis of this paper, some basic indications of the computational complexity of the SABRE method will be given here together with some pointers to a more efficient D-SABRE method. In the form presented in this paper, the SABRE method has a computational complexity that is of the same order as that of the FD method, i.e.,

O (N^{3})

for the construction of the reconstruction matrix and

O (N^{2})

for the reconstruction itself.

In particular, the computational complexity of the SABRE method is

O ({(rank N_{A})}^{3})

for the construction of the reconstruction matrix, and

O ({(rank N_{A})}^{2})

for the WFR operation, with

N_{A}

a basis for the nullspace of the constraint matrix

A

from (45). The rank of

N_{A}

in turn is a function of the total number of triangles and internal triangle edges in a triangulation and the degree and continuity order of the SABRE model.

However, the local nature of the basis functions of the simplex B-splines enable the definition of a D-SABRE method that solves the global

N \times N

WFR problem by splitting it into

G

parts. Each of the

G

parts is formed by a subtriangulation of the global triangulation, and the WFR problem is solved for each subtriangulation individually using only local sensor data. This would reduce the computational complexity of the reconstruction problem from

O (N^{2})

O (N^{2} / G^{2})

when run on

G

parallel processors, leading to a speedup factor of

O (G^{2})

over the ordinary global SABRE method.

K. Example of the Linear SABRE

In this paragraph we will present a tutorial example of the linear SABRE method. In the example, we use the Simplex Type-I sensor geometry from Fig. 3. We assume that there are two wavefront slope sensor apertures at which the local wavefront slope in the

(x, y)

directions is measured. Furthermore, we assume a linear SABRE model structure of the form (44). The corresponding spline function is defined on a triangulation consisting of two triangles, with a symmetric B-net orientation [see Fig. 5].

Fig. 5. Simplex Type-I sensor geometry (left) and the triangulation and B-net for a linear simplex B-spline function (right).

Before starting with the reconstruction, the various entities in (44) must be declared. We start with the vector of B-coefficients, which in this case consists of six B-coefficients:

c = {[\begin{matrix} c_{1, 0, 0}^{t_{1}} & c_{0, 1, 0}^{t_{1}} & c_{0, 0, 1}^{t_{1}} & c_{1, 0, 0}^{t_{2}} & c_{0, 1, 0}^{t_{2}} & c_{0, 0, 1}^{t_{2}} \end{matrix}]}^{⊤} .

The construction of the de Casteljau matrix

P^{1, 0}

from (44) requires us to first define directional coordinates for the

x

and

y

directions for each of the two triangles. Using (24), we find for the directional coordinates

a_{x}^{t_{1}}

and

a_{y}^{t_{1}}

with respect to

t_{1}

and the directional coordinates

a_{x}^{t_{2}}

and

a_{y}^{t_{2}}

with respect to

t_{2}

a_{x}^{t_{1}} = {[\begin{matrix} 0 & 1 & - 1 \end{matrix}]}^{⊤}, a_{y}^{t_{1}} = {[\begin{matrix} 1 & 0 & - 1 \end{matrix}]}^{⊤}, a_{x}^{t_{2}} = {[\begin{matrix} - 1 & 0 & 1 \end{matrix}]}^{⊤}, a_{y}^{t_{2}} = {[\begin{matrix} 0 & - 1 & 1 \end{matrix}]}^{⊤} .

(56)

Using (56), together with (27), the full-triangulation de Casteljau matrix

P^{1, 0}

P^{1, 0} = [\begin{matrix} 0 & 1 & - 1 & 0 & 0 & 0 \\ - 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & - 1 \\ 0 & 0 & 0 & 0 & - 1 & 1 \end{matrix}] .

We have two continuity conditions of the form (21):

- c_{1, 0, 0}^{t_{1}} + c_{1, 0, 0}^{t_{2}} = 0

and

- c_{0, 1, 0}^{t_{1}} + c_{0, 1, 0}^{t_{2}} = 0

, which lead to the smoothness matrix

H

given by

H = [\begin{matrix} - 1 & 0 & 0 & 1 & 0 & 0 \\ 0 & - 1 & 0 & 0 & 1 & 0 \end{matrix}] .

(57)

The anchor constraint is constructed using (41)

h = [\begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 \end{matrix}] .

(58)

The complete constraint matrix

A

is formed by combining (57) with (58) as shown in (45). We now proceed with reformulating the constrained model structure from (44) into the equivalent constrained form from (51) by use of the following basis for the null space of

A

N_{A} = {[\begin{matrix} 0 & 0 & - 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}]}^{⊤} .

We are now ready to define the reduced regression matrix

D

from (51):

D = P^{1, 0} N_{A} = [\begin{matrix} 1 & 1 & 0 \\ - 1 & 0 & 0 \\ 0 & 0 & - 1 \\ 0 & - 1 & 1 \end{matrix}],

which we can use to define an LS estimator for

c

using (54).

4. SIMULATIONS WITH THE SABRE

In this section the WFR performance of the SABRE method is compared to that of the Southwell and Fried FD methods in two numerical experiments utilizing a Fourier optics based SH WFS simulation.

A. Setup of the Numerical Experiments

For the experiments, simulated turbulence phase screens are used in combination with a Fourier optics based SH WFS simulation to generate simulated wavefront phase slope measurements. The simulated SH WFS sensor effectively performs a local spatial averaging of the wavefront at each lenslet.

A total of 100 simulated phase screens are created using a von Kármán turbulence model [30

J. M. Conan, G. Rousset, and P. Y. Madec, “Wave-front temporal spectra in high-resolution imaging through turbulence,” J. Opt. Soc. Am. 12, 1559–1570 (1995). [CrossRef]

,31

R. Conan, “Mean-square residual error of a wavefront after propagation through atmospheric turbulence and after correction with Zernike polynomials,” J. Opt. Soc. Am. 25, 526–536 (2008). [CrossRef]

]. For the von Kármán turbulence model a turbulence outer scale of 50 m is used with a Fried coherence length of 0.2 m. The diameter of the simulated aperture is 2 m. No continuous inner scale is assumed. The average phase range for the 100 phase screens is 23 rad. The WFR takes place on the grid

X_{r}

which is of the same resolution as the simulated SH lenslet array used in the simulation, i.e.,

17 \times 17

. The phase screen

ϕ

is generated on a

272 \times 272

grid indicated by

X_{wf}

such that each lenslet images a

16 \times 16

segment of

ϕ

The overall quality of the reconstruction is calculated on the high resolution grid

X_{wf}

. For this, the reconstructed wavefronts must be interpolated. For the SABRE method, this interpolation is implicit, as the SABRE model is an analytical function that can be evaluated at any point inside the model domain. For the FD method, two different interpolation schemes are used; a linear and a cubic interpolation scheme.

Three performance metrics are defined to assess the performance of the WFR methods. The first is the estimation error that is defined as follows:

ϵ = {‖ \hat{ϕ} - ϕ ‖}_{2}

(59)

with

\hat{ϕ}

the estimated wavefront that is interpolated on

X_{wf}

. The normalized estimation error

E_{R}

is given by

E_{R} = \frac{ϵ}{{‖ ϕ ‖}_{2}} .

(60)

The third performance metric is the noise gain. The noise gain is defined as the ratio of

ϵ

to the standard deviation of the noise, indicated as

σ_{ν}

, for a flat phase screen (i.e.,

ϕ = 0

K_{ν} = {\frac{ϵ}{σ_{ν}} |}_{ϕ = 0} .

(61)

B. Fourier Optics-Based SH Simulator

In the numerical experiments, a simulated SH lenslet array is used to produce slope measurements, see Fig. 6. The SH lenslet array simulation is based on the Fresnel approximation to the Maxwell equations, and was taken from [32

Y. Dai, F. Li, X. Cheng, Z. Jiang, and S. Gong, “Analysis on Shack–Hartmann wave-front sensor with fourier optics,” Opt. Laser Technol. 39, 1374–1379 (2007). [CrossRef]

]. For the numerical experiments, a

K \times L = 17 \times 17

lenslet array is used for a total of 289 lenslets. Each (circular) lenslet has a diameter of

D_{l} = 150 μm

and a focal distance of

f = 5 \cdot 10^{- 3} m

. The observation wavelength is

λ = 633 nm

. For the nonrectangular WFS numerical experiments, a predefined mask is placed on the lenslet array, resulting in a nonrectangular array with a total of 240 lenslets, see Fig. 7.

Fig. 6. SH array image for a turbulence phase screen realization used in the numerical experiment.

Fig. 7. Four SABRE sensor geometries used in the numerical experiments: A Type II geometry (top left) containing 1024 triangles, Type-IIG geometry (top right) containing 256 triangles, Type-IIB geometry containing 816 triangles (bottom left), and a Type-IIBG geometry containing 208 triangles (bottom right). The open circles are the locations of the SH-lenslets.

A simulated CCD with

4352 \times 4352

pixels is placed in the focal plane of the lenslet array. It is assumed that the pupil plane function of each lenslet acts on a

16 \times 16

pixel segment of the global wavefront

ϕ

. In order to increase the resolution of the focal plane lenslet image, this segment is embedded at the center of a larger

M_{p} \times M_{p} = 256 \times 256

grid containing only zeros.

The focal plane intensity function

I_{k l}

for each of the lenslets is then calculated as follows:

I_{k l} = | F F {A_{k l} \cdot e^{i \cdot ϕ_{k l}}} | \in R^{256 \times 256}

(62)

with

ϕ_{k l} \in R^{256 \times 256}

the wavefront phase segment imaged by the lenslet

(k, l)

A_{k l} (m, n) \in R^{256 \times 256}

the pupil mask of lenslet

(k, l)

, and

F F

the zero-shifted two-dimensional discrete Fourier transform.

The complete SH-array focal intensity function

I (x, y) \in R^{4352 \times 4352}

is then constructed using the lenslet intensity functions from (62) as follows:

I (k \cdot M_{p} + m, l \cdot M_{p} + n) = I_{k l} (m, n), m = 1, 2, \dots, M_{p}, n = 1, 2, \dots, M_{p}

(63)

for each

k = 0, 1, 2, \dots, K - 1

and each

l = 0, 1, 2, \dots, L - 1

The local wavefront slope is obtained by calculating the shift of the center of mass (COM) of each lenslet image with respect to a reference COM. The COM for each lenslet image is calculated as follows:

x_{c_{k l}} = \frac{\sum_{m = 0}^{M_{p} - 1} \sum_{n = 0}^{M_{p} - 1} m \cdot I (k \cdot M_{p} + m, l \cdot M_{p} + n)}{\sum_{m = 0}^{M_{p} - 1} \sum_{n = 0}^{M_{p} - 1} I (k \cdot M_{p} + m, l \cdot M_{p} + n)}, y_{c_{k l}} = \frac{\sum_{m = 0}^{M_{p} - 1} \sum_{n = 0}^{M_{p} - 1} n \cdot I (k \cdot M_{p} + m, l \cdot M_{p} + n)}{\sum_{m = 0}^{M_{p} - 1} \sum_{n = 0}^{M_{p} - 1} I (k \cdot M_{p} + m, l \cdot M_{p} + n)} .

(64)

Finally, the slopes at lenslet

(k, l)

are calculated using the result from (64) as follows:

σ_{x_{k l}} = \frac{(x_{c_{k l}} - x_{0_{k l}})}{f} + w ν σ_{y_{k l}} = \frac{(y_{c_{k l}} - y_{0_{k l}})}{f} + w ν

(65)

with

x_{0_{k l}}

and

y_{0_{k l}}

the calibrated (reference) centroid locations, with

0.01 \cdot ‖ σ ‖ \leq w \leq 10 \cdot ‖ σ ‖

the variable noise magnitude and

ν

a uniformly distributed random number sequence in the interval

[\begin{matrix} - 1 & 1 \end{matrix}]

C. Reconstruction Performance on Rectangular SH Arrays

In the first numerical experiment a quantitative comparison is made between the reconstruction performance of the SABRE and the Fried and Southwell FD methods using wavefront slope measurements obtained with the simulated rectangular

17 \times 17

SH array described in Section 4.B. For the SABRE, a Type-II triangulation consisting of 1024 triangles and a Type-IIG triangulation consisting of 256 triangles are used. Linear, quadratic, and cubic SABRE model structures are used in the experiment. The LS estimator from Section 3.H is used to estimate the B-coefficients of the SABRE models. For the FD methods, both linear and cubic interpolators are used to interpolate from the SH grid to the validation grid

X_{wf}

The reconstruction performance of the SABRE is compared with that of the Fried and Southwell FD methods for 20 different signal-to-noise ratio (SNR) values. For every SNR value, the results from 100 turbulence realizations and 50 noise realizations are averaged. In total, there are 100,000 reconstructions for each reconstruction experiment. The quality of the resulting WF models is assessed on the validation grid

X_{wf}

. In Fig. 8 an example is shown of a quadratic SABRE and a linear Fried FD reconstruction of a single (zero noise) turbulence realization.

Fig. 8. Linear Fried FD reconstruction (bottom left) and a quadratic SABRE reconstruction (bottom right) of a turbulence phase screen.

In Table 2 the realization-averaged normalized error (

{\bar{E}}_{R}

) for the zero noise case is shown for all tested reconstruction methods. These results show that only the quadratic SABRE model with first order continuity has a lower minimum

{\bar{E}}_{R}

than the linear Fried FD model. On the other hand, all SABRE models have a lower minimum

{\bar{E}}_{R}

than the cubic interpolated Fried FD model and both Southwell FD models. In Table 2 it can also be seen that the averaged noise gain

{\bar{K}}_{ν}

for all but one of the SABRE models is lower than that of both Fried FD models. Note that both Southwell FD models have a lower average noise gain than all other models.

Table 2. Results of the Zero-Noise Experiments

SABRE/FD Space	Rectangular Domain	$M in ({\bar{E}}_{R})$	$M ax ({\bar{E}}_{R})$	${\bar{K}}_{ν}$
Southwell (linear)	Yes	0.2127	1.2152	0.0388
Southwell (cubic)	Yes	0.2118	1.2738	0.0407
Fried (linear)	Yes	0.0533	1.5751	0.0514
Fried (cubic)	Yes	0.0630	1.9096	0.0582
SABRE ( $d = 1$ , $r = 0$ , $J = 1024$ )	Yes	0.0535	1.3801	0.0426
SABRE ( $d = 1$ , $r = 0$ , $J = 816$ )	No	0.0523	1.4519	0.0474
SABRE ( $d = 2$ , $r = 1$ , $J = 1024$ )	Yes	0.0525	1.4335	0.0442
SABRE ( $d = 2$ , $r = 1$ , $J = 816$ )	No	0.0510	1.4946	0.0487
SABRE ( $d = 3$ , $r = 1$ , $J = 256$ )	Yes	0.0579	1.6233	0.0498
SABRE ( $d = 3$ , $r = 1$ , $J = 208$ )	No	0.0570	1.6929	0.0547

In Fig. 9 the performance in terms of

{\bar{E}}_{R}

of the various WFR methods is plotted as a function of the SNR. The figure shows that the quadratic SABRE model has a higher performance than the linear and cubic Fried FD models for all SNR values, although the difference is most pronounced at high noise levels. The linear SABRE model has a higher performance than both Fried FD models for SNR

< 15

, but has a somewhat lower performance for SNR

> 15

. The Southwell FD models have a higher performance than the SABRE models for SNR

< 0.3

, but for SNR

> 0.6

all SABRE models outperform the Southwell FD method. Note that the cubic interpolated Fried FD method has a lower performance than the linear interpolated Fried FD method as well as all SABRE models for all noise levels; this is due to the Lagrangian bicubic interpolation operation that effectively smooths the Fried FD phase points with cubic polynomials of full continuity order and consequently lower approximation power than the linear Fried FD and SABRE models. Note also that in Fig. 9 only the results from the linearly interpolated Southwell FD method are presented; it produced virtually identical results to its cubic interpolated counterpart.

Fig. 9. Average normalized residual

{\bar{E}}_{R}

as a function of signal to noise ratio for a number of different linear and nonlinear SABRE models compared with the Southwell and Fried reconstructors using wavefront slope and curvature measurements. A total of 100 turbulence realizations and 50 noise realizations are used per SNR magnitude.

D. Reconstruction Performance on Nonrectangular SH Arrays

In the second numerical experiment the reconstruction performance of the SABRE on nonrectangular sensor arrays is investigated. For this experiment, an octagonal mask with a central obstruction is used on the

17 \times 17

SH array described in Section 4.B [see Fig. 7]. The nonobscured lenslet reference centroids are triangulated with a Type-IIB consisting of 816 triangles and a Type-IIBG triangulation consisting of 208 triangles. The same linear, quadratic, and cubic SABRE model structures are used as in the previous experiment. The LS estimator from Section 3.H is used to estimate the B-coefficients of the SABRE models. It is important to note that the SABRE algorithm for the nonrectangular SH array is exactly equal to that of the rectangular SH array.

In Table 2 the realization-averaged normalized error of the nonrectangular SABRE models for the zero-noise case are presented. In Fig. 10 three plots are shown that compare the realization-averaged normalized estimation error

{\bar{E}}_{R}

of SABRE models of equal degree and continuity order defined on rectangular and nonrectangular SH arrays. The plots show that the performance of the SABRE models on nonrectangular arrays is higher than SABRE models of the same degree defined on rectangular arrays for low noise levels while performance is somewhat lower for high noise levels.

Fig. 10. Investigation of the influence of sensor geometry on the average normalized residual

{\bar{E}}_{R}

for SABRE models of equal degree and continuity order. A total of 100 turbulence realizations and 50 noise realizations are used for each noise setting.

5. CONCLUSIONS

A new zonal method for WFR using multivariate simplex B-spline basis functions in a linear regression framework is introduced. This new method, indicated as the SABRE, approximates solutions to the first order partial differential equations that relate wavefront slope to the wavefront phase using bivariate simplex B-splines. While the current paper only introduces an LS estimator for the B-coefficients of the SABRE model, the SABRE is fully compatible with any more advanced linear regression parameter estimation technique. The SABRE method is aimed at AO applications in which linear WFR methods cannot provide the required accuracy levels at a given sensor noise level, or in AO applications in which wavefront slope or curvature measurements are made on nonrectangular, partially obstructed, and misaligned sensor grids. Furthermore, the SABRE reconstruction matrix is of full rank, and as a consequence the SABRE is not subject to the waffle mode.

The new method is compared to existing FD methods from the literature on both a theoretical and a numerical level. On a theoretical level it is proved that any FD method is a special case of a linear SABRE model with

C^{0}

continuity on a Type-II triangulation. Additionally, we prove that the approximation power of a class of nonlinear SABRE models exceeds the approximation power of the class of linear SABRE models. This class of nonlinear SABRE models will ultimately result in reconstructed wavefronts of higher accuracy than those produced by any linear SABRE, and consequently, any FD model.

The computational complexity of the SABRE method in its basic form is

O (N^{2})

for the reconstruction operation. This is equal to the MVM method used for the FD method. However, the local nature of the simplex B-spline basis polynomials allows the definition of a D-SABRE algorithm that can theoretically obtain a computational complexity of

O (N^{2} / G^{2})

when a total of

G

parallel processors are used. Additionally, the sparse matrix methods from [8

B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. 19, 1803–1816 (2002). [CrossRef]

C. R. Vogel, “Sparse matrix methods for wavefront reconstruction, revisited,” Proc. SPIE 5490, 1327–1335 (2004).

] and the multigrid preconditioned conjugate-gradient methods from [10

L. Gilles, C. R. Vogel, and B. L. Ellerbroek, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. 19, 1817–1822 (2002). [CrossRef]

,11

C. R. Vogel and Q. Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt. 45, 705–715 (2006). [CrossRef]

] can be used with the SABRE to further increase computational efficiency.

Two numerical experiments are conducted using a Fourier optics based simulation of a SH lenslet array. Both experiments use turbulence phase screens generated using a von Kármán turbulence model together with a variable magnitude sensor noise model. In the first numerical experiment, the reconstruction performance of the SABRE method is compared to that of the Fried and Southwell FD methods on rectangular SH lenslet arrays. In this experiment it is shown that the quadratic SABRE model outperforms all tested FD methods for all but the highest noise levels. In the second experiment the reconstruction performance of the SABRE on nonrectangular SH sensor arrays is investigated. It is found that reconstruction performance is not negatively influenced by the particular geometry of the underlying sensor grid.

Appendices

APPENDIX A

In this appendix a theorem is provided that on a theoretical level connects the SABRE method to the FD method. In particular, it is proved that any linear FD model is in fact a special case of a linear SABRE model on a Type-II triangulation. Additionally, it is proved that there exist nonlinear SABRE models that can approximate a wavefront with higher accuracy than any linear FD method.

A. SABRE as a Generalization of the FD Method

In the following sections, we require a definition for the space of linear FD approximators, that is, the space that contains all possible linear FD reconstructions for a given sensor grid. We indicate this space as

F_{1}^{0} (Ω (M, N))

and define it as follows:

Definition 1.

F_{1}^{0} (Ω (M, N))

is defined as the space of all linear FD approximators of continuity order

C^{0}

on the

M \times N

grid

Ω (M, N)

from (32) such that we have for every (estimated) phase point

{\hat{ϕ}}_{FD} (x_{i}, y_{j})

{\hat{ϕ}}_{FD} (x_{i}, y_{j}) \in F_{1}^{0} (Ω (M, N)), 1 \leq i \leq M, 1 \leq j \leq N .

(A1)

The following lemma proves that any linear FD model is in fact a special case of a linear SABRE model defined on a Type-II triangulation.

Lemma 1. Let

S_{1}^{0} (T_{II} (Ω (M, N)))

be the linear spline space of

C^{0}

continuity defined on the Type-II triangulation

T_{II} (Ω (M, N))

with

Ω (M, N)

the FD grid domain from (32). In that case, the space of FD approximators

F_{1}^{0} (Ω (M, N))

from (A1) is a subset of this spline space as follows:

F_{1}^{0} (Ω (M, N)) \subset S_{1}^{0} (T_{II} (Ω (M, N))) .

(A2)

Proof. The proof is based on the existence of an injective mapping of the FD phase points to the B-coefficients of the spline space. First, let

ω

be any grid cell in

Ω (M, N)

consisting of the four vertices

V_{ω} = {v_{i, j}, v_{i + 1, j}, v_{i, j + 1}, v_{i + 1, j + 1}}

. These vertices thus coincide exactly with the location of the FD phase points. Then, let

T_{4} (ω)

be the Type-II triangulation of

V_{ω}

with vertex set

V_{T_{4}} = {V_{ω}, v_{*}}

in which

v_{*}

is an additional vertex located at the center of

ω

. Now let

s_{1}^{0} (x, y)

be a linear bivariate spline function with

C^{0}

continuity defined on

T_{4} (ω)

as follows:

s_{1}^{0} (x, y) = B^{1} (b (x, y)) c_{ω} \in S_{1}^{0} (T_{4} (ω))

(A3)

with

c_{ω}

a vector of B-coefficients that can be decomposed as follows:

c_{ω} = {[\begin{matrix} c_{v}^{⊤} & c_{*}^{⊤} \end{matrix}]}^{⊤} \in R^{12}

(A4)

with

c_{v} \in R^{8}

the B-coefficients located at the vertices in

V_{ω}

, and with

c_{*} \in R^{4}

the B-coefficients located at the center vertex

v_{*}

From (21) it follows that for

C^{0}

continuity, it is required that the values of all B-coefficients located at a vertex in

V_{T_{4}}

are equal. As a result, we find the following for

c_{v}

and

c_{*}

c_{v} = [\begin{matrix} c_{v_{1}} & c_{v_{1}} & c_{v_{2}} & c_{v_{2}} & c_{v_{3}} & c_{v_{3}} & c_{v_{4}} & c_{v_{4}} \end{matrix}], c_{*} = [\begin{matrix} c_{v_{*}} & c_{v_{*}} & c_{v_{*}} & c_{v_{*}} \end{matrix}] .

(A5)

Evaluation of the spline function (A3) at any vertex

v_{m} \in V_{ω}

results in

s_{1}^{0} (v_{m}) = B^{1} (b (v_{m})) \cdot c_{ω} = c_{v_{m}}, 1 \leq m \leq 4

(A6)

with

b (v_{m})

the Barycentric coordinate transformation of

v_{m}

For the theory to hold, the FD phase points must be equal to the SABRE wavefront phase at all vertices in

V_{ω}

. From (A6) it therefore follows that

{\hat{ϕ}}_{FD} (v_{m}) = c_{v_{m}}, 1 \leq m \leq 4

(A7)

with

{\hat{ϕ}}_{FD} (v_{m})

the FD wavefront phase as defined in (A1).

Now let

c_{ϕ}

be the vector of B-coefficients resulting from the substitution of the results from (A7) in

c_{v}

from (A5):

c_{ϕ} = [\begin{matrix} {\hat{ϕ}}_{FD} (v_{1}) & {\hat{ϕ}}_{FD} (v_{1}) & \dots & {\hat{ϕ}}_{FD} (v_{4}) & {\hat{ϕ}}_{FD} (v_{4}) \end{matrix}] .

(A8)

Equation (A8) implies that there exists an injective mapping of

{\hat{ϕ}}_{FD} (V_{ω})

c_{ω}

. Evaluation of the spline function from (A3) at the vertex set

V_{T_{4}}

then results in

s_{1}^{0} (V_{T_{4}}) = B^{1} (b (V_{T_{4}})) \cdot [\begin{matrix} c_{ϕ} \\ c_{*} \end{matrix}] = [\begin{matrix} {\hat{ϕ}}_{FD} (V_{ω}) \\ v_{*} \end{matrix}]

(A9)

with

b (V_{T_{4}})

the transformation to Barycentric coordinates of the vertices in

V_{T_{4}}

Equation (A9) implies that the set of all FD approximators

{\hat{ϕ}}_{FD} (V_{ω})

is a subset of

S_{1}^{0} (T_{4} (ω))

for any real value of the B-coefficients

c_{*}

located at the center vertex

v_{*}

F_{1}^{0} (ω) \subset S_{1}^{0} (T_{4} (ω)) .

(A10)

By induction, the result from (A10) can be extended to the complete domain

Ω (M, N)

, thereby proving the lemma.

The following lemma reproduces theorems from [22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

] that prove that the approximation accuracy of a given wavefront phase screen

ϕ (x, y)

with simplex B-splines depends on the spline degree as well as on the configuration of a triangulation.

Lemma 2 ([22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

]: Theorem 10.2, Theorem 10.6, Theorem 10.13, Theorem 10.21) Let

S_{d}^{r} (T (Ω))

be the spline space of degree

d

and continuity order

r

defined on the triangulation

T (Ω)

of the domain

Ω

. If

ϕ (x, y)

is a differentiable function with at least

d + 1

derivatives such that

ϕ (x, y) \in W_{q}^{m} (Ω)

, with

W_{q}^{d + 1} (Ω)

the Sobolev space of order

m

equipped with the

L_{q}

norm, then we have for the

q

-distance

d {(ϕ (x, y), S_{d}^{r})}_{q}

between

ϕ (x, y)

and

S_{d}^{r}

d {(ϕ (x, y), S_{d}^{r})}_{q} = {\begin{cases} O ({| T |}^{0}) & if d < \frac{3 r + 2}{2}, r > 0 \\ O ({| T |}^{d}) & if \frac{3 r + 2}{2} \leq d \leq 2 r + 1, r > 0 \\ O ({| T |}^{d + 1}) & if d \geq 3 r + 2, r \geq 0 \end{cases}

with

| T |

the longest triangle edge in the triangulation

T

Proof. For a proof of the lemma, we refer the reader to Theorem 10.2, Theorem 10.6, Theorem 10.13, and Theorem 10.21 in [22

M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

The results of Lemma 1 and Lemma 2 lead to the main theorem, which states that there exists an estimator that the SABRE method can reconstruct the wavefront

ϕ (x, y)

with equal or higher accuracy than the FD method.

Theorem 1. If

f \in F_{1}^{0} (Ω (M, N))

is an FD approximator in the space of linear FD approximators from (A1), then there exists a spline function

s_{d}^{r} \in S_{d}^{r} (T_{II} (Ω (M, N)))

defined on a Type-II triangulation of

Ω (M, N)

for which it holds that it can approximate any continuous function

ϕ (x, y)

with higher accuracy than any linear FD approximator:

{‖ ϕ (x, y) - s_{d}^{r} ‖}_{q} < {‖ ϕ (x, y) - f ‖}_{q}

(A11)

for all

d > 1

r = 0

, and for all

d \geq (3 r + 2) / 2

r > 0

Proof. The proof is based on the results from Lemma 1 and Lemma 2. In Lemma 1 it was proved that

F_{1}^{0} (Ω (M, N))

is in fact a subset of

S_{1}^{0} (T_{II} (Ω (M, N)))

, which implies the following:

{‖ ϕ (x, y) - s_{1}^{0} ‖}_{q} \leq {‖ ϕ (x, y) - f ‖}_{q}

(A12)

with

s_{1}^{0} \in S_{1}^{0} (T_{II} (Ω (M, N)))

Additionally, in Lemma 2 it was proved that the distance between

ϕ (x, y)

and

S_{d}^{r} (T)

depends on

| T |

as well as on

d

and

r

. In particular it was proved that

d (ϕ (x, y), S_{d}^{r} (T))_{q} < d (ϕ (x, y), S_{1}^{0} (T))_{q},

(A13)

holds for all

d > 1

r = 0

, and for all

d \geq (3 r + 2) / 2

r > 0

. Combining (A12) with (A13) then proves the theorem.

REFERENCES

1.	J. M. Beckers, P. Lena, O. Lai, P. Y. Madec, G. Rousset, M. Séchaud, M. J. Northcott, F. Roddier, J. L. Beuzit, F. Rigaut, and D. G. Sandler, Adaptive Optics in Astronomy (Cambridge University, 1999).
2.	D. L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. 67, 370–375 (1977). [CrossRef]
3.	J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]
4.	F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. 27, 1223–1225 (1988). [CrossRef]
5.	W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. 70, 998–1006 (1980). [CrossRef]
6.	M. Kissler-Patig, “Overall science goals and top level AO requirements for the E-ELT,” presented at First AO4ELT Conference, Victoria, Canada, B.C., September 25 and 30, 2010.
7.	V. Korkiakoski and C. Vérinaud, “Simulations of the extreme adaptive optics system for EPICS,” Proc. SPIE 7736, 773643 (2010).
8.	B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. 19, 1803–1816 (2002). [CrossRef]
9.	C. R. Vogel, “Sparse matrix methods for wavefront reconstruction, revisited,” Proc. SPIE 5490, 1327–1335 (2004).
10.	L. Gilles, C. R. Vogel, and B. L. Ellerbroek, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. 19, 1817–1822 (2002). [CrossRef]
11.	C. R. Vogel and Q. Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt. 45, 705–715 (2006). [CrossRef]
12.	M. Rosensteiner, “Cumulative reconstructor: fast wavefront reconstruction algorithm for extremely large telescopes,” J. Opt. Soc. Am. 28, 2132–2138 (2011). [CrossRef]
13.	G. M. Dai, “Modal wave-front reconstruction with Zernike polynomials and Karhunen–Loève functions,” J. Opt. Soc. Am. 13, 1218–1225 (1996). [CrossRef]
14.	R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66, 207–211 (1976). [CrossRef]
15.	L. A. Poyneer, D. T. Gavel, and J. M. Brase, “Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform,” J. Opt. Soc. Am. 19, 2100–2111 (2002). [CrossRef]
16.	L. A. Poyneer, B. A. Machintosh, and J. P. Veran, “Fourier transform wavefront control with adaptive prediction of the atmosphere,” J. Opt. Soc. Am. 24, 2645–2660 (2007). [CrossRef]
17.	P. J. Hampton, P. Agathoklis, and C. Bradley, “A new wave-front reconstruction method for adaptive optics systems using wavelets,” IEEE J. Select. Topics Signal Process. 2, 781–792 (2008). [CrossRef]
18.	G. Awanou, M. J. Lai, and P. Wenston, “The multivariate spline method for scattered data fitting and numerical solutions of partial differential equations,” in Wavelets and Splines , G. Chen and M. J. Lai, eds. (Nashboro, 2005), pp. 24–75.
19.	C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).
20.	C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]
21.	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]
22.	M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).
23.	M. D. Oliker, “Sensing waffle in the Fried geometry,” Proc. SPIE 3353, 964–971 (1998). [CrossRef]
24.	W. Zou and J. P. Rolland, “Quantifications of error propagation in slope-based wavefront estimations,” J. Opt. Soc. Am. 23, 2629–2638 (2006). [CrossRef]
25.	C. de Boor, “B-form basics,” in Geometric Modeling: Algorithms and New Trends (SIAM, 1987).
26.	X. L. Hu, D. F. Han, and M. J. Lai, “Bivariate splines of various degrees for numerical solution of partial differential equations,” SIAM J. Sci. Comput. 29, 1338–1354 (2007). [CrossRef]
27.	M. J. Lai and L. L. Schumaker, “On the approximation power of bivariate splines,” Adv. Comput. Math. 9, 251–279(1998). [CrossRef]
28.	M. J. Lai, “Some sufficient conditions for convexity of multivariate Bernstein–Bezier polynomials and box spline surfaces,” Studia Scient. Math. Hung. 28, 363–374 (1990). [CrossRef]
29.	R. H. Hudgin, “Wavefront reconstruction for compensated imaging,” J. Opt. Soc. Am. 67, 375–378(1977). [CrossRef]
30.	J. M. Conan, G. Rousset, and P. Y. Madec, “Wave-front temporal spectra in high-resolution imaging through turbulence,” J. Opt. Soc. Am. 12, 1559–1570 (1995). [CrossRef]
31.	R. Conan, “Mean-square residual error of a wavefront after propagation through atmospheric turbulence and after correction with Zernike polynomials,” J. Opt. Soc. Am. 25, 526–536 (2008). [CrossRef]
32.	Y. Dai, F. Li, X. Cheng, Z. Jiang, and S. Gong, “Analysis on Shack–Hartmann wave-front sensor with fourier optics,” Opt. Laser Technol. 39, 1374–1379 (2007). [CrossRef]

OCIS Codes
(000.3860) General : Mathematical methods in physics
(000.4430) General : Numerical approximation and analysis
(010.1080) Atmospheric and oceanic optics : Active or adaptive optics
(010.7350) Atmospheric and oceanic optics : Wave-front sensing
(350.1260) Other areas of optics : Astronomical optics
(010.1285) Atmospheric and oceanic optics : Atmospheric correction

ToC Category:
Atmospheric and Oceanic Optics

History
Original Manuscript: June 29, 2012
Revised Manuscript: September 27, 2012
Manuscript Accepted: October 20, 2012
Published: December 13, 2012

Citation
Cornelis C. de Visser and Michel Verhaegen, "Wavefront reconstruction in adaptive optics systems using nonlinear multivariate splines," J. Opt. Soc. Am. A 30, 82-95 (2013)
http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-30-1-82

Sort: Author | Year | Journal | Reset

References

J. M. Beckers, P. Lena, O. Lai, P. Y. Madec, G. Rousset, M. Séchaud, M. J. Northcott, F. Roddier, J. L. Beuzit, F. Rigaut, and D. G. Sandler, Adaptive Optics in Astronomy (Cambridge University, 1999).
D. L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. 67, 370–375 (1977). [CrossRef]
J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]
F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. 27, 1223–1225 (1988). [CrossRef]
W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. 70, 998–1006 (1980). [CrossRef]
M. Kissler-Patig, “Overall science goals and top level AO requirements for the E-ELT,” presented at First AO4ELT Conference, Victoria, Canada, B.C., September 25 and 30,2010.
V. Korkiakoski and C. Vérinaud, “Simulations of the extreme adaptive optics system for EPICS,” Proc. SPIE 7736, 773643 (2010).
B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. 19, 1803–1816 (2002). [CrossRef]
C. R. Vogel, “Sparse matrix methods for wavefront reconstruction, revisited,” Proc. SPIE5490, 1327–1335 (2004).
L. Gilles, C. R. Vogel, and B. L. Ellerbroek, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. 19, 1817–1822 (2002). [CrossRef]
C. R. Vogel and Q. Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt. 45, 705–715 (2006). [CrossRef]
M. Rosensteiner, “Cumulative reconstructor: fast wavefront reconstruction algorithm for extremely large telescopes,” J. Opt. Soc. Am. 28, 2132–2138 (2011). [CrossRef]
G. M. Dai, “Modal wave-front reconstruction with Zernike polynomials and Karhunen–Loève functions,” J. Opt. Soc. Am. 13, 1218–1225 (1996). [CrossRef]
R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66, 207–211 (1976). [CrossRef]
L. A. Poyneer, D. T. Gavel, and J. M. Brase, “Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform,” J. Opt. Soc. Am. 19, 2100–2111(2002). [CrossRef]
L. A. Poyneer, B. A. Machintosh, and J. P. Veran, “Fourier transform wavefront control with adaptive prediction of the atmosphere,” J. Opt. Soc. Am. 24, 2645–2660 (2007). [CrossRef]
P. J. Hampton, P. Agathoklis, and C. Bradley, “A new wave-front reconstruction method for adaptive optics systems using wavelets,” IEEE J. Select. Topics Signal Process. 2, 781–792 (2008). [CrossRef]
G. Awanou, M. J. Lai, and P. Wenston, “The multivariate spline method for scattered data fitting and numerical solutions of partial differential equations,” in Wavelets and Splines, G. Chen and M. J. Lai, eds. (Nashboro, 2005), pp. 24–75.
C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).
C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]
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]
M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).
M. D. Oliker, “Sensing waffle in the Fried geometry,” Proc. SPIE 3353, 964–971 (1998). [CrossRef]
W. Zou and J. P. Rolland, “Quantifications of error propagation in slope-based wavefront estimations,” J. Opt. Soc. Am. 23, 2629–2638 (2006). [CrossRef]
C. de Boor, “B-form basics,” in Geometric Modeling: Algorithms and New Trends (SIAM, 1987).
X. L. Hu, D. F. Han, and M. J. Lai, “Bivariate splines of various degrees for numerical solution of partial differential equations,” SIAM J. Sci. Comput. 29, 1338–1354 (2007). [CrossRef]
M. J. Lai and L. L. Schumaker, “On the approximation power of bivariate splines,” Adv. Comput. Math. 9, 251–279(1998). [CrossRef]
M. J. Lai, “Some sufficient conditions for convexity of multivariate Bernstein–Bezier polynomials and box spline surfaces,” Studia Scient. Math. Hung. 28, 363–374 (1990). [CrossRef]
R. H. Hudgin, “Wavefront reconstruction for compensated imaging,” J. Opt. Soc. Am. 67, 375–378(1977). [CrossRef]
J. M. Conan, G. Rousset, and P. Y. Madec, “Wave-front temporal spectra in high-resolution imaging through turbulence,” J. Opt. Soc. Am. 12, 1559–1570 (1995). [CrossRef]
R. Conan, “Mean-square residual error of a wavefront after propagation through atmospheric turbulence and after correction with Zernike polynomials,” J. Opt. Soc. Am. 25, 526–536 (2008). [CrossRef]
Y. Dai, F. Li, X. Cheng, Z. Jiang, and S. Gong, “Analysis on Shack–Hartmann wave-front sensor with fourier optics,” Opt. Laser Technol. 39, 1374–1379 (2007). [CrossRef]

Agathoklis, P.
Awanou, G.
Beckers, J. M.
Beuzit, J. L.
Bradley, C.
Brase, J. M.
Cheng, X.
Chu, Q. P.
Conan, J. M.
Conan, R.
Dai, G. M.
Dai, Y.
de Boor, C.
de Visser, C. C.
Ellerbroek, B. L.
Fried, D. L.
Gavel, D. T.
Gilles, L.
Gong, S.
Hampton, P. J.
Han, D. F.
Herrmann, J.
Hu, X. L.
Hudgin, R. H.
Jiang, Z.
Kissler-Patig, M.
Korkiakoski, V.
Lai, M. J.
Lai, O.
Lena, P.
Li, F.
Machintosh, B. A.
Madec, P. Y.
Mulder, J. A.
Noll, R. J.
Northcott, M. J.
Oliker, M. D.
Poyneer, L. A.
Rigaut, F.
Roddier, F.
Rolland, J. P.
Rosensteiner, M.
Rousset, G.
Sandler, D. G.
Schumaker, L. L.
Séchaud, M.
Southwell, W. H.
Veran, J. P.
Vérinaud, C.
Vogel, C. R.
Wenston, P.
Yang, Q.
Zou, W.

Adv. Comput. Math.

M. J. Lai and L. L. Schumaker, “On the approximation power of bivariate splines,” Adv. Comput. Math. 9, 251–279(1998). [CrossRef]

Appl. Opt.

Automatica

C. C. de Visser, Q. P. Chu, and J. A. Mulder, “A new approach to linear regression with multivariate splines,” Automatica 45, 2903–2909 (2009). [CrossRef]
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]

IEEE J. Select. Topics Signal Process.

P. J. Hampton, P. Agathoklis, and C. Bradley, “A new wave-front reconstruction method for adaptive optics systems using wavelets,” IEEE J. Select. Topics Signal Process. 2, 781–792 (2008). [CrossRef]

J. Opt. Soc. Am.

M. Rosensteiner, “Cumulative reconstructor: fast wavefront reconstruction algorithm for extremely large telescopes,” J. Opt. Soc. Am. 28, 2132–2138 (2011). [CrossRef]
G. M. Dai, “Modal wave-front reconstruction with Zernike polynomials and Karhunen–Loève functions,” J. Opt. Soc. Am. 13, 1218–1225 (1996). [CrossRef]
R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am. 66, 207–211 (1976). [CrossRef]
L. A. Poyneer, D. T. Gavel, and J. M. Brase, “Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform,” J. Opt. Soc. Am. 19, 2100–2111(2002). [CrossRef]
L. A. Poyneer, B. A. Machintosh, and J. P. Veran, “Fourier transform wavefront control with adaptive prediction of the atmosphere,” J. Opt. Soc. Am. 24, 2645–2660 (2007). [CrossRef]
W. H. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. 70, 998–1006 (1980). [CrossRef]
D. L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. 67, 370–375 (1977). [CrossRef]
J. Herrmann, “Least-squares wave front errors of minimum norm,” J. Opt. Soc. Am. 70, 28–35 (1980). [CrossRef]
B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. 19, 1803–1816 (2002). [CrossRef]
L. Gilles, C. R. Vogel, and B. L. Ellerbroek, “Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction,” J. Opt. Soc. Am. 19, 1817–1822 (2002). [CrossRef]
W. Zou and J. P. Rolland, “Quantifications of error propagation in slope-based wavefront estimations,” J. Opt. Soc. Am. 23, 2629–2638 (2006). [CrossRef]
R. H. Hudgin, “Wavefront reconstruction for compensated imaging,” J. Opt. Soc. Am. 67, 375–378(1977). [CrossRef]
J. M. Conan, G. Rousset, and P. Y. Madec, “Wave-front temporal spectra in high-resolution imaging through turbulence,” J. Opt. Soc. Am. 12, 1559–1570 (1995). [CrossRef]
R. Conan, “Mean-square residual error of a wavefront after propagation through atmospheric turbulence and after correction with Zernike polynomials,” J. Opt. Soc. Am. 25, 526–536 (2008). [CrossRef]

Opt. Laser Technol.

Y. Dai, F. Li, X. Cheng, Z. Jiang, and S. Gong, “Analysis on Shack–Hartmann wave-front sensor with fourier optics,” Opt. Laser Technol. 39, 1374–1379 (2007). [CrossRef]

Proc. SPIE

M. D. Oliker, “Sensing waffle in the Fried geometry,” Proc. SPIE 3353, 964–971 (1998). [CrossRef]
V. Korkiakoski and C. Vérinaud, “Simulations of the extreme adaptive optics system for EPICS,” Proc. SPIE 7736, 773643 (2010).

SIAM J. Sci. Comput.

X. L. Hu, D. F. Han, and M. J. Lai, “Bivariate splines of various degrees for numerical solution of partial differential equations,” SIAM J. Sci. Comput. 29, 1338–1354 (2007). [CrossRef]

Studia Scient. Math. Hung.

M. J. Lai, “Some sufficient conditions for convexity of multivariate Bernstein–Bezier polynomials and box spline surfaces,” Studia Scient. Math. Hung. 28, 363–374 (1990). [CrossRef]

Other

C. de Boor, “B-form basics,” in Geometric Modeling: Algorithms and New Trends (SIAM, 1987).
J. M. Beckers, P. Lena, O. Lai, P. Y. Madec, G. Rousset, M. Séchaud, M. J. Northcott, F. Roddier, J. L. Beuzit, F. Rigaut, and D. G. Sandler, Adaptive Optics in Astronomy (Cambridge University, 1999).
C. R. Vogel, “Sparse matrix methods for wavefront reconstruction, revisited,” Proc. SPIE5490, 1327–1335 (2004).
M. Kissler-Patig, “Overall science goals and top level AO requirements for the E-ELT,” presented at First AO4ELT Conference, Victoria, Canada, B.C., September 25 and 30,2010.
G. Awanou, M. J. Lai, and P. Wenston, “The multivariate spline method for scattered data fitting and numerical solutions of partial differential equations,” in Wavelets and Splines, G. Chen and M. J. Lai, eds. (Nashboro, 2005), pp. 24–75.
C. C. de Visser, “Global nonlinear model identification with multivariate splines,” Ph.D. thesis (Delft University of Technology, 2011).
M. J. Lai and L. L. Schumaker, Spline Functions on Triangulations (Cambridge University, 2007).

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper