Reconstruction of Wavefront (no DM) from subaperture slopes

Top  Previous  Next

The goal here is to develop a reconstructor mapping that directly converts sdata into a reconstructed wavefront opd vector dopd : dopd,rec = R sdata . The displacements dopd  are defined with respect to a flat, 0-tilt reference wavefront.  The following sections explain alternate approaches for which code modules exist in the AO tools.

Zonal reconstructor for wavefront

To perform a zonal WF reconstruction using the LightLike AO tools, the approach is to sequentially execute the functions aogeom, lsfptos, and aorecon. aogeom is used to establish the geometry of subapertures and their corner "actuator" positions, as introduced previously.  In the present context, the significance of the so-called "actuator" positions is simply that they define mesh points at which the WF phase will be computed.  We regret the somewhat confusing dual-use terminology.  

In the present reconstruction calculation, misregistration is not allowed:  the "actuators" in aogeommust be placed at the subaperture corners, or the subsequent calculation using lsfptos makes no sense.  Also, there is no significance to slave actuators in this situation: all the actuators should be designated as masters in preparation for using the lsfptos function.

Using the geometry information, the function lsfptos carries out the following calculation.  Consider one subaperture, as illustrated in Figure 5.  The labels ci denote the subaperture corners.  The dopd,idenote the WF displacements at the corners, with respect to a 0-tilt plane wave across the entire complex of subapertures.  One can approximate the average surface slopes of this WF patch by simple bilinear interpolation, namely

       sx = 0.5 * [(dopd,2 -dopd,1)/dx + (dopd,4 - dopd,3)/dx]

       sy = 0.5 * [(dopd,3 -dopd,1)/dy + (dopd,4 - dopd,2)/dy]


Figure 5:  A single subaperture in a Fried geometry

For the whole set of subapertures, one can use these equations and adjacency information to create the matrix relation

      s = S'c dopd,c  ,  where matrix S'c has dims (2nsub x nopd,c).

We use the subscript "opd,c" to emphasize the fact that nopd,c is the number of subaperture corner points.  Note that in the present procedure, the opd mesh has the same spacing as the subaperture mesh (though offset to the corners).  In some later procedures involving DM influence functions, we will introduce opd meshes that are denser than the subaperture mesh.  The Matlab routine lsfptosprovided with the AO tools computes S'c , assuming that dopd,c is in meters, and s is in radians of angle.

CAUTION: the routine lsfptos actually omits the (1/dx,dy) factors, which must be appended manually to obtain the correctly-scaled slope angles in radians of angle.

Now, if we have WFS slope data, sdata, then we can obtain a "reconstructed" WF, dopd,c;rec , by solving the equation sdata = S'c dopd,c  in a least-squares sense.  A least-squares (or some analogous) approach is necessary since S'c has more rows than columns (there are more equations than unknowns).  The AO tools provide several elaborations on the classical least-squares solution, so we will use the term "pseudo-inverse" to generically designate the solution:

       sdata = S'c dopd,c    ==>   dopd,c;rec = (S'c)# sdata = R'c sdata

We use the superscript notation (...)#  to indicate the "pseudo-inverse".  More generally, we will use the notation R to represent a reconstructor matrix (a matrix that transforms a slope data vector into wavefront opd or into actuator commands, as the case may be).  Note that the reconstructor (S'c)# remains fixed, for a given WFS geometry, as the incident wavefront changes in a dynamic system.  The AO tools provide the function aorecon to compute the pseudo-inverse of any specified matrix.

In sum, to perform a zonal WF reconstruction, we require the tools functions aogeom (define subaperture geometry), lsfptos (compute S'c matrix), and aorecon (compute pseudo-inverse (S'c)#).  Full details of the available options, and examples of Matlab terminal sessions, are given in later sections devoted to lsfptos and aorecon.  Background information regarding one key option of aorecon, namely the singular-value decomposition, is discussed next.


Least-squares solution and singular-value decomposition

Given any overdetermined matrix system (more rows than columns) written as

      s = S d,

the classical least-squares solution, derived in numerous textbooks, can be shown to be equivalent to solving the matrix equation

       ST s  =  STS d,   where T=transpose, and (STS)  is a square matrix, presumed to be of full rank.

Assuming (STS)  is not exactly singular, the classical solution can now be written in the matrix form

      drec =  (STS)-1 ST s  =  Rclas s  ,  where (.)-1 designates a true inverse.

We use the symbol R in general to denote a "reconstructor" matrix that maps a given s to the "reconstructed" drec, and Rclas stands for the classical least-squares solution  (STS)-1 ST.

Although the square matrix  (STS)  is usually not exactly singular, so that a numerical inverse does exist, it turns out in practice that  (STS)  often has a very high condition number.  This is equivalent to saying that the matrix is in a certain sense "close" to singular, or that the solution Rclas can produce results with poor numerical stability.  Poor stability in this context means that small changes in the data s (due, e.g., to measurement noise) produce large changes in drec.  In practice, this can introduce meaningless solution jumps and stability problems with the control system.  Modern numerical analysis has developed a method of solving the least-squares problem that allows one to identify and eliminate "modes" which are associated with the numerical instability.  This method is based on a matrix factorization of rectangular matrices called the singular-value decomposition (SVD).  The AO tools function aorecon allows the user to first compute the SVD, and then to compute the reconstructor matrix R after specifying which modes to remove.  This requires executing aorecon twice, with two different sets of arguments.  If no modes are removed, the solution reduces to the classical Rclas.  The later section on aorecon illustrates the Matlab syntax details and usage options.  

The SVD option of aorecon can be used regardless of the precise context in which aorecon is used.  Note that removing the near-singular modes is often critical to the construction of a stable least-squares solution.


(1)  G.E. Forsythe, M.A. Malcolm, and C.B. Moler, Computer Methods for Mathematical Computations, chs. 3 and 9, Prentice-Hall, 1977.

(2)  G. Strang, Linear Algebra and its Applications, 2nd ed., Academic Press, 1980.

(3)  Golub and VanLoan, Matrix Computations, Chs. 2 and 6, Johns Hopkins U. Press, 1983.

Modal reconstructor for wavefront

To perform a modal WF reconstruction using the LightLike AO tools, the approach is to sequentially execute the functions aogeom, slope2mod, and mod2opd. aogeom is used to establish the geometry of the subapertures, as introduced previously.  Even though no DM is present, the corner "actuator" positions in aogeom may be used later to define the points at which the WF phase will be reconstructed;  however, the modal reconstructor can directly yield the reconstructed WF at other points as well.

The previous section explained "zonal" reconstruction:  in that procedure, samples of wavefront opd were defined at the corners of every subaperture "zone", and these samples were considered independent variables on an equal footing, to be solved for via a least-squares fit.  A strong point of the zonal procedure is that it works easily with any WFS geometry.  In particular, the algorithm for computing the S'c influence matrix is the same whether or not there is a central obscuration, or an irregular boundary.  The only thing that is needed is the identification of each subaperture with the corner points that surround it.

An alternate reconstruction procedure, known as "modal reconstruction", begins by expressing the wavefront opd as a superposition of some convenient set of mathematical basis functions ("modes").  In this case, the expansion coefficients become the unknown elements to be reconstructed, after truncating the sum at some finite order.  The elements to be solved for do not directly represent local ("zonal") displacements, but rather the modal amplitudes.  

A feature of modal reconstruction is that it inherently contains an extra degree of freedom relative to the zonal reconstructor, namely the number of modes for which the user wants to solve.  Modal reconstruction has several purposes:

1.There is important theoretical work on AO correction that is based on modal decomposition of the wavefront.  Thus, using a modal reconstructor with data allows one to make close contact with this theory.

2.One may wish to implement an AO system that only corrects a selected set of modes.  This may be desired for various reasons:  (a) complexity or cost limitations, (b) use of more than one corrector element with different properties for different modes, (c) use of post-processing to correct higher-order aberrations, (d) experimental situations (e.g., thermal blooming) where one wishes to ignore certain modes because of special control instabilities that can arise.

3.The subtleties of numerical conditioning (solution stability) may play out differently in the inversions associated with modal and zonal expansions.

A disadvantage of modal reconstruction is that it is more difficult to treat general WFS geometries.  There are several related complications, but the basic problem is that one wants to use basis functions that are orthogonal on the support of the optical pupil.  Therefore, different geometries require algorithm development with different basis sets.  The modal routines presently available in the AO tools allow the modeling of unobstructed circular and annular pupils, via the Zernike circle polynomials and Zernike annular polynomials, respectively.

Consider the following expansion of the wavefront opd:

      dopd(r) = Σ{k=1...nmod}  ak Zk(r)

where the Zk are the basis functions, ak are the expansion coefficients, nmod modes have been retained, and r is an arbitrary position.  The AO tools implement a procedure in which the Zk basis is the Zernike circle or annular polynomial set.  For a discrete set of sample points {ri}, we express the preceding equation in matrix form as

      dopd = Z amod

where  Zi,k = Zk(ri) , amod = [a1,...,anmod ]T,  and dopd = [dopd(r1),...,dopd(rNopd)]T .

The goal of the modal reconstruction is to estimate amod , and finally dopd , from sdata.  The AO tools implement a procedure developed by Gavrielides.  The key concepts are as follows:

(1)  In continuum space, the wavefront slope is the gradient of dopd(r):

      s(r) = grad{dopd(r)} = Σ{k=1...nmod}  ak grad{Zk(r)}

(2)  Therefore, if we had a set of two-vector functions that were orthogonal to the gradient of the Zernikes, we could immediately apply a projection integral to the preceding equation to solve analytically for any ak.

(3) Gavrielides derived a set of vector polynomials that satisfy (2).

(4)  By discretizing the projection integral on the subaperture mesh, we obtain the matrix formulas

      amod, rec  =  V sdata ,  and therefore

      dopd, rec  =  Z amod, rec =  Z V sdata  =  Rmod sdata .

The notation  Rmod =  ZV  defines the final reconstructor matrix that converts a slope data vector into a wavefront opd.

After aogeom is used to define the geometry of the subapertures, the AO tools provide two further functions that implement the modal reconstruction described above:  function slope2mod computes the V matrix for a specified WFS geometry, and function mod2opd computes the Z matrix for a specified opd mesh.  Syntax details of slope2mod and mod2opd are given in a later section.

Notice that this reconstruction method does not require a least-squares solution to obtain the wavefront.  However, when we later extend the procedure to the computation of DM actuator commands that will correct the perturbed wavefront, we will see that a least-squares solution does enter.


A. Gavrielides, "Vector polynomials orthogonal to the gradient of Zernike polynomials", Optics Letters, Vol. 7, No. 11, November 1982, pp. 526-528.

Zonal 2 (Southwell) reconstructor for wavefront

An alternate method of zonal WF reconstruction is due to Southwell.  To perform the Southwell reconstruction using the LightLike AO tools, the approach is to sequentially execute the functions aogeom, ...

This section is under construction.


Southwell, " "