Implementation and benchmarking of a crosstalk-free method for wavefront Zernike coefficients reconstruction using Shack-Hartmann sensor data

R. S. Biesheuvel; A. J. E. M. Janssen; P. Pozzi; S. F. Pereira

doi:10.1364/OSAC.1.000581

1. Introduction

Determination of wavefronts is an important field in optics. In adaptive optics, for example, the wavefront is constantly monitored in order to feed the servo system that corrects distortions on it due to the optical system self or due to the environment. Adaptive optics started in the field of astronomy, but nowadays it has been extended, among others, to the fields of microscopy and retinal imaging. Apart from that, it is worth mentioning the application in accurate characterization of complex optical systems for aberrations.

In many applications involving wavefront reconstruction, Shack-Hartmann sensors are frequently used. In these sensors, an array of microlenses focuses the incoming beam into a camera, resulting in a grid of spots that are displaced from their center position (flat wave reference) because of aberrations contained in the wavefront. The displacement of each spot is then determined and used to calculate the local slope of the wavefront across each lens of the array. Once the slopes are known, an often-used polynomial basis to expand the phase of the scalar field (wavefront) is the set of Zernike circle polynomials. The wavefront Zernike coefficients can be estimated with for example, the least-squares method, to be called Method A henceforth. However, already since the 1980s, it is known that there is a fundamental problem in the combination of the least-squares method and Zernike circle polynomials: this method presents crosstalk between coefficients. In 2014, Janssen introduced an analytic result in the form of a new relation between the wavefront derivative data and the wavefront Zernike coefficients [1]. This relation theoretically solves the crosstalk problem that is present in the least-squares method. In particular, the wavefront coefficients, once estimated, do not change anymore when the number of Zernike circle polynomials involved in the fit is increased further.

In this paper, we show an implementation of Janssen’s analytic result in the reconstruction of several wavefronts that have been generated using a spatial light modulator. In particular, we show that this method, to be called Method B henceforth, does not suffer from cross-talk between the coefficients.

The paper is organized as follows. In Section 2 we present theoretical aspects of the problem, and in Section 3, the experimental setup and methods are explained. The actual results are presented and discussed in Section 4, and the conclusions from these measurements are presented in Section 5.

2. Theory

In this section, Zernike circle polynomials are discussed and used as a basis to describe the wavefront deviation with respect to a suitably defined reference sphere of an optical system. This wavefront deviation or wavefront aberration function will be loosely called ’wavefront aberration’ or ’aberration’ in what follows. We start this section by showing the definition of Zernike circle polynomials and a description of the data that is obtained with the Shack-Hartmann sensor is discussed. Next, we discuss two reconstruction methods: the well-known least-squares method (Method A) and the recently method introduced by Janssen [1] (Method B).

2.1. Zernike circle polynomials as a basis for wavefront expansion

In systems where aberrations are desired to be known, such as microscopes and telescopes, a circular aperture is usually present. In order to describe these aberrations, Zernike circle polynomials are commonly used. In this paper, the American National Standards Institute (ANSI) definition for Zernike circle polynomials is considered, which is commonly used to describe wavefront reconstruction from Shack-Hartmann data [2]. As defined in [3], the polynomials in polar coordinates (ρ, θ) are given by

z_{n}^{m} (ρ, θ) = N_{n}^{m} R_{n}^{| m |} (ρ) Θ_{m} (θ),

where

N_{n}^{m} = \sqrt{(2 - δ_{m 0}) (n + 1)},

R_{n}^{| m |} (ρ) = \sum_{s = 0}^{\frac{n - | m |}{2}} \frac{{(- 1)}^{s} (n - s)!}{s! (\frac{n - | m |}{2} - s)! (\frac{n + | m |}{2} - s)!} ρ^{n - 2 s},

Θ_{m} (θ) = {\begin{array}{l} \cos (m θ), & if m \geq 0, \\ - \sin (m θ), & if m < 0, \end{array}

such that

z_{n}^{m} (ρ, θ)

is a real-valued, orthonormal (on the unit disc) expression for Zernike circle polynomials, n = 0, 1, 2, … is the degree of the Zernike polynomial, and m = 0, ±1, ±2, … the azimuthal order satisfying:

n \geq 0,

n - | m | is even,

| m | \leq n,

ρ ≤ 1 being the radius on the unit disc, and $δ_{n n^{'}}$ the Kronecker delta function. Here, the n^th “order” polynomial is referred to as the n^th degree polynomial.

As mentioned above, these circle polynomials are orthonormal on the unit disc, i.e.,

\frac{1}{π} \int_{0}^{1} \int_{0}^{2 π} z_{n}^{m} (ρ, θ) z_{n^{'}}^{m^{'}} (ρ, θ) ρ d θ d ρ = δ_{n n^{'}} δ_{m m^{'}} .

With this definition, any real-valued wavefront function (defined on the unit disc) can be described as a linear combination of the Zernike circle polynomials given by

W (ρ, θ) = \sum_{m = - \infty}^{\infty} \sum_{n \in η_{m}} a_{n}^{m} z_{n}^{m} (ρ, θ),

where η_m is the set of allowed n values dependent on m, namely η_m = {|m|, |m| + 2, |m| + 4, …}, when m ≠ 0, and η_m = {2, 4, …, ∞} when m = 0 and

z_{0}^{0} = 1

. This set ensures that the constraints set in n and m are all met. Finally,

a_{n}^{m}

is the real-valued Zernike coefficient.

2.2. Complex-valued definition of the Zernike circle polynomials

Following the complex-valued definition of the Zernike circle polynomials as used in [1] in polar coordinates (ρ, θ) on the unit disc, we have that

Z_{n}^{m} (ρ, θ) = R_{n}^{| m |} (ρ) e^{i m θ} .

The radial polynomial

R_{n}^{| m |}

is given as in Eq. (1) above, and can also be given as

R_{n}^{| m |} (ρ) = ρ^{| m |} P_{\frac{n - | m |}{2}}^{(0, | m |)} (2 ρ^{2} - 1),

where

P_{k}^{(α, β)} (x)

is the Jacobi polynomial of degree k, which is orthogonal with respect to the weight (1 − x)^α(1 + x)^β on the interval [−1, 1]. Note that the factor ρ^|^m| is missing in [1], Eq. (6), but it is merely a typo and has no consequences in the further developments in that reference. We set

Z_{n}^{m} = 0

for all values of n and m where n − |m| is odd or negative. There is the normalization condition

\int_{0}^{1} \int_{0}^{2 π} Z_{n}^{m} (ρ, θ) {(Z_{n^{'}}^{m^{'}} (ρ, θ))}^{*} ρ d θ d ρ = \frac{π}{n + 1} δ_{n n^{'}} δ_{m m^{'}} .

This orthogonality means that, like the real-valued Zernike polynomial, any sufficiently smooth (complex) wavefront can be described by

W (ρ, θ) = \sum_{m = - \infty}^{\infty} \sum_{n \in η_{m}} α_{n}^{m} Z_{n}^{m} (ρ, θ),

where

α_{n}^{m}

are generally complex-valued coefficients corresponding to the complex Zernike polynomial and η_n is as in Eq. (9). Expanding the complex exponential in Eq. (10) leads to the following conversion between this definition of complex Zernike cirlce polynomials and the ANSI standard

N_{n}^{m} Z_{n}^{m} (ρ, θ) = {\begin{array}{l} z_{n}^{| m |} (ρ, θ) + i z_{n}^{- | m |} (ρ, θ), & if m > 0, \\ z_{n}^{| m |} (ρ, θ) - i z_{n}^{- | m |} (ρ, θ), & if m < 0, \\ z_{n}^{m}, & if m = 0, \end{array}

where

N_{n}^{m}

is defined in Eq. (2).

From the definition of the complex Zernike polynomial in Eq. (10) we can also see that the complex conjugate ${(Z_{n}^{| m |})}^{*}$ is equal to $Z_{n}^{- | m |}$ . This observation, together with the relations between the real and complex Zernike circle polynomials leads to the expression

z_{n}^{m} = {\begin{array}{l} \frac{N_{n}^{m}}{2} (Z_{n}^{| m |} + Z_{n}^{- | m |}) = N_{n}^{m} ℜ (Z_{n}^{| m |}), & if m > 0 \\ \frac{N_{n}^{m}}{2 i} (Z_{n}^{| m |} - Z_{n}^{- | m |}) = N_{n}^{m} J (Z_{n}^{| m |}), & if m < 0 \\ N_{n}^{m} Z_{n}^{| m |}, & if m = 0 . \end{array}

When using complex polynomials to describe a real wavefront, the coefficients are usually complex valued. In order to transform them back to the coefficients of the real-valued Zernike $z_{n}^{m}$ , the following relations can be used

a_{n}^{m} = {\begin{array}{l} \frac{1}{N_{n}^{m}} ℜ (α_{n}^{| m |} + α_{n}^{- | m |}), & if m > 0 \\ \frac{- 1}{N_{n}^{m}} J (α_{n}^{| m |} - α_{n}^{- | m |}), & if m < 0 \\ \frac{1}{N_{n}^{m}} ℜ (α_{n}^{m}), & if m = 0 . \end{array}

Here, the Zernike circle polynomials are evaluated on the computer, and data is saved in vectors and matrices. To make it easy to loop over all polynomials, a single index is introduced. In the notation by [3], this single index j = 0, 1, … is given by

j = \frac{n (n + 2) + m}{2},

where n is the Zernike degree and m the azimuthal order. From now on, we will write

Z_{j} = z_{n}^{m}

. The reverse can also be done, that is finding n and m from j as follows

n = ⌈ \frac{- 3 + \sqrt{9 + 8 j}}{2} ⌉,

m = 2 j - n (n + 2),

where ⌈x⌉ denotes the “ceiling” function, that is the smallest integer greater or equal than x.

2.3. Data from the Shack-Hartmann sensor

It is well known that aberrations in an optical system degrade its imaging quality. The phase aberration is defined by a function φ_W = arg(W). A typical way of measuring this function is by using a Shack-Hartmann sensor. A Shack-Hartmann sensor consists of a camera chip and lenslet array. The lenslets are placed at the focal distance from the camera chip such that an incoming plane wave will be focused as many spots on the camera. If the incoming wave is aberrated, the position of each spot on the sensor will changed.

It can be derived from the theory how much a spot is displaced due to an aberration. Because this displacement is linear, the reverse problem can be solved as well. In other words, one can retrieve the aberration given the spot displacement.

In Ref. [4], it has proven that for a varying wavefront over the subaperture, the slope needs to be averaged over that sub-aperture. The expressions for the average slopes in terms of displacement of the spot on the camera becomes

{\begin{array}{l} \frac{1}{A_{Σ}} \int_{Σ} \frac{\partial W}{\partial x} d x d y & = r \frac{Δ x}{f} \\ \frac{1}{A_{Σ}} \int_{Σ} \frac{\partial W}{\partial y} d x d y & = r \frac{Δ y}{f} \end{array}

where Δx, Δy are the shift in the x– and y– positions of the spots, r is the radius of the incoming beam on the Shack-Hartmann sensor, and f is the focal length of the sensor lens array, Σ is the illuminated sub-aperture domain, with surface area A_Σ. Here Σ and A_Σ change when the sub-aperture is only partially illuminated (i.e., at the edge of the beam). This averaging of the slope is of great matter in recovering the wavefront.

2.3.1. Method A

The least-squares (LSQ) fit is based on the real Zernike circle polynomials, and uses the fact that the p coefficients are not dependent on x and y. The wavefront in x, y coordinates (with $\sqrt{x^{2} + y^{2}} \leq 1$ 1) can be described as

W (x, y) = \sum_{m = - \infty}^{\infty} \sum_{n \in η_{m}} a_{n}^{m} z_{n}^{m} (x, y) .

The measured quantity, however, is not W but $\frac{\partial W}{\partial x}$ and $\frac{\partial W}{\partial y}$ . Taking the partial derivatives to x and y results in the over determined system of

{\begin{cases} \frac{\partial W}{\partial x} = \sum_{m = - \infty}^{\infty} \sum_{n \in η_{m}} a_{n}^{m} \frac{\partial z_{n}^{m}}{\partial x} \\ \frac{\partial W}{\partial y} = \sum_{m = - \infty}^{\infty} \sum_{n \in η_{m}} a_{n}^{m} \frac{\partial z_{n}^{m}}{\partial y} . \end{cases}

This system can be solved for a finite set of polynomials. Using Eq. (20), given the displacements Δx and Δy, one can create a vector s containing the slopes as such

s = {[{\bar{\frac{\partial W}{\partial x}} |}_{1} {\bar{\frac{\partial W}{\partial x}} |}_{2} \dots {\bar{\frac{\partial W}{\partial x}} |}_{n_{spots}} {\bar{\frac{\partial W}{\partial y}} |}_{1} {\bar{\frac{\partial W}{\partial y}} |}_{2} \dots {\bar{\frac{\partial W}{\partial y}} |}_{n_{spots}}]}^{T} .

The partial derivatives of the Zernikes in the x- and y-direction can also be put in a matrix, called the geometry matrix. We recall the convention that

Z_{j} = z_{n}^{m}

with j, m, n related as in Eqs. (17)–(19). The geometry matrix G can be built up

{\bar{\frac{\partial Z_{j}}{\partial x}} |}_{n}

, the average gradient of Zernike mode j at the position of subaperture n. This averaging is done due to the fact that the spot displacement measured with the Shack-Hartmann sensor is proportional to the average slope of the wavefront, as expressed in Eq. (20). It should be noted that the positions over which the averaging is done is normalized to the unit disc. These windows are the same defined in Eq. (20). The matrix will have a size of (2n_spot × J), where J is the maximum index for the Zernike modes used to have a good approximation of the true wavefront. The expression for G becomes

G = {[\begin{matrix} {\bar{\frac{\partial Z_{1}}{\partial x}} |}_{1} & {\bar{\frac{\partial Z_{1}}{\partial x}} |}_{2} & \dots & {\bar{\frac{\partial Z_{1}}{\partial x}} |}_{n_{spot}} & {\bar{\frac{\partial Z_{1}}{\partial y}} |}_{1} & {\bar{\frac{\partial Z_{1}}{\partial y}} |}_{2} & \dots & {\bar{\frac{\partial Z_{1}}{\partial y}} |}_{n_{spot}} \\ {\bar{\frac{\partial Z_{2}}{\partial x}} |}_{1} & {\bar{\frac{\partial Z_{2}}{\partial x}} |}_{2} & \dots & {\bar{\frac{\partial Z_{2}}{\partial x}} |}_{n_{spot}} & {\bar{\frac{\partial Z_{2}}{\partial y}} |}_{1} & {\bar{\frac{\partial Z_{2}}{\partial y}} |}_{2} & \dots & {\bar{\frac{\partial Z_{2}}{\partial y}} |}_{n_{spot}} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ {\bar{\frac{\partial Z_{J}}{\partial x}} |}_{1} & {\bar{\frac{\partial Z_{J}}{\partial x}} |}_{2} & \dots & {\bar{\frac{\partial Z_{J}}{\partial x}} |}_{n_{spot}} & {\bar{\frac{\partial Z_{J}}{\partial y}} |}_{1} & {\bar{\frac{\partial Z_{J}}{\partial y}} |}_{2} & \dots & {\bar{\frac{\partial Z_{J}}{\partial y}} |}_{n_{spot}} \end{matrix}]}^{T} .

The system of Eq. (22) can then be written as

s \approx G \cdot a,

where a is the vector containing the Zernike coefficients. The least-squares estimation of

a_{n}^{m}

becomes

a \approx G^{+} \cdot s,

where G⁺ the generalized inverse of the geometry matrix. This is an approximation as G only contains the information of a finite number of Zernike modes, and their contribution is averaged over the lenslet array.

2.3.2. Method B

Method B relies on an analytical relation found between the local derivatives of the wavefront and Zernike circle polynomials. This is in contrast with method A, where there is a link between the local derivatives of the wavefront and the derivatives of the Zernike circle polynomials.

In Ref. [1], it has found that the LMS complex coefficients are given as (see Appendix)

{\hat{α}}_{n}^{m} = C_{n}^{m} φ_{n}^{m} - C_{n + 2}^{m} φ_{n + 2}^{m},

where

C_{n}^{m} = \frac{1 + δ_{n | m |}}{2 n}

φ_{n}^{m} = \frac{1}{2} {(β_{+})}_{n - 1}^{m + 1} + \frac{1}{2} {(β_{-})}_{n - 1}^{m - 1},

and where

δ_{n n^{'}}

is the Kronecker delta equal to 1 if n = n′ and 0 otherwise. Furthermore β₊ and β₋ are the Zernike coefficients of

\frac{\partial W}{\partial x} \pm i \frac{\partial W}{\partial y}

so that

\frac{\partial W}{\partial x} \pm i \frac{\partial W}{\partial y} = \sum_{m = - \infty}^{\infty} \sum_{n \in η_{m}} {(β_{\pm})}_{n}^{m} Z_{n}^{m} .

Note that the LMS coefficients ${\hat{α}}_{n}^{m}$ are analytically related to four β coefficients, namely ${(β_{+})}_{n + 1}^{m + 1}$ , ${(β_{+})}_{n - 1}^{m + 1}$ , ${(β_{-})}_{n - 1}^{m - 1}$ , and ${(β_{-})}_{n + 1}^{m - 1}$ . As a consequence of orthogonality of the $Z_{n}^{m}$ in the right-hand side of Eq. (30), the LMS coefficients, once estimated, do not change anymore when the number of Zernike terms used in a finitized version of Eq. (30) is increased further.

The fact that ${\hat{α}}_{n}^{m}$ is analytically related to β coefficients is desirable, because β-coefficients can directly be estimated (in a least-squares sense) from measurable quantities. This fit is made in the same way as the least-squares method. This means that also the complex Zernike circle polynomials need to be averaged over the lenslets. The fit to get the β coefficients, however, is done with a different basis than in the least-squares method to find the a-coefficients. The effects of this are discussed in the following section.

As a note for this method, when n = |m|, there will be non-existent combinations of n and m in Eq. (29). In that case the value of β will be set to 0. For instance, ${\hat{α}}_{1}^{1}$ is among others dependent on ${(β_{+})}_{0}^{2}$ , which goes against the constraint |m| ≤ n. To go from complex coefficient $α_{n}^{m}$ to the real coefficient $a_{n}^{m}$ , the relations in Eq. (16) can be used.

In [5], a seemingly different approach, based on vector polynomials, is used to express the wavefront Zernike coefficients in terms of wavefront slope data. However, we show in the appendix that this method is essentially equivalent to Janssen’s algorithm in [1].

2.3.3. Wavefront reconstruction using Method A and Method B

The main difference between Method A and Method B for finding the coefficients is in how the fitting is done. Both methods use a least-squares fit using a geometry matrix, but the matrix elements are constructed differently. In Method A, the geometry matrix elements are evaluations of the average gradient of the real-valued Zernike circle polynomials over certain subdisks, while for Method B, the average of the complex-valued Zernike circle polynomials over the same windows [see Eq. (30)].

The gradients of the Zernike circle polynomials are known not to be orthogonal. This can cause problems called crosstalk when fitting the coefficients, especially when there are more aberrations present in the system than are being fit.

If a is an M-dimensional vector containing the coefficients of the aberrations present in the system, the slopes on the Shack-Hartmann sensor can be determined as s = Ga, where s is an 2n_spot long vector containing the x- and y-displacement on the Shack-Hartmann sensor and G an 2n_spot × J geometry matrix defined in Eq. (24).

The notation $Z_{j} = z_{n}^{m}$ is used since j, n, m are related [see Eq. (17)–(19)]. Z₀ is not included in the matrix since $Z_{0} = z_{0}^{0} = 1$ , and its partial derivatives are equal to 0. Note that the first column contains the x- and y-derivatives of the first Zernike polynomial evaluated in all n_spot points. When a least-squares estimation of the coefficients $\hat{a}$ (where the hat means to indicate that it is an estimated parameter) is made using less Zernike circle polynomials, up to Zernike polynomial J < M, crosstalk will occur. The estimator $\hat{a}$ can be expressed as

\begin{array}{l} \hat{a} \approx G_{l}^{+} s \\ \hat{a} \approx G_{l}^{+} G a, \end{array}

where G_l is the geometry matrix containing the columns of the first J Zernike circle polynomials. The estimator will estimate the lower-order values of the coefficients with influence of the higher-order values, because the matrix

G_{l}^{+} G

will not be an identity matrix.

When estimating the coefficient $a_{n}^{m}$ , a higher-order aberration $a_{n^{'}}^{m^{'}}$ will influence the estimation if it is not accounted for in G_l (i.e., the single index of $a_{n^{'}}^{m^{'}} j > J$ ) and if

{(n, m, n^{'}, m^{'}) \in ℤ | n \in η_{m}, m = m^{'} or m = m^{'} \pm 2, n^{'} > n, n \geq m^{'} + 2, n^{'} \in η_{m^{'}}, m^{'} \neq 0},

or if

{(n, m, n^{'}, m^{'}) \in ℤ | n \in η_{m}, m = 0 or m = 2, n^{'} > n, n \geq 2, n^{'} \in η_{m^{'}}, m^{'} = 0},

where in both cases η_m and

η_{m^{'}}

η_m′ are the sets of allowed values for n and n′ dependent on m and m′ such that Eqs. (5)–(7) are all met.

Because of this dependence of an expansion coefficient on the maximum degree of the system of equations, it is expected that methods such as the least-squares method will incorrectly estimate the coefficients when there are higher-order aberrations present that are not accounted for in the geometry matrix G_l. For Method B, the geometry matrix contains the Zernike circle polynomials themselves, and therefore it is not expected to present any crosstalk. This is experimentally verified and shown in this paper.

3. Experimental setup and methods

In Fig. 1, a schematic view of the experimental setup is shown. The phase of an uniform, collimated laser beam (HeNe laser) is modified by a spatial light modulator (SLM, Holoeye PLUTO-2-VIS-056) in such a way that known aberrations are added to the beam. In order to remove the zeroth- order reflected light from the SLM, a phase ramp is added to all SLM patterns so that the unaberrated spot is blocked by the iris after being focused by lens L3. The lenses L1 and L2 and L3 and L4 form two 4 f systems that projects the phase imposed on the SLM onto the Shack Hartmann sensor (SHS). For each set of aberrations added in the SLM we obtain experimentally the raw data of the SHS, i.e., the image of the x-y displaced spots detected by the camera of the SHS. This data is used to find the average slopes of the wavefront as explained in Section 2.3.

Fig. 1 Scheme of the experimental setup. FR is a Faraday isolator, WP is a halfwave plate, SLM is a spatial light modulator, SHS is the Shack Hartmann wavefront sensor. The zeroth order light of the SLM is blocked by the iris.

Download Full Size | PDF

Initially, the system is tested to identify its intrinsic aberrations (from lenses, SLM, etc). After calibration of the system response, theses aberrations are canceled by the SLM (flattening). Afterwards, additional and controlled aberrations are applied to the SLM with specific strength of a few Zernike coefficients. The raw measurements made with the SHS are further processed and used to compare the reconstruction methods A and B. In Algorithm 1 we describe all steps of the experiments.

Algorithm 1. Complete measurement and comparison of Shack-Hartmann phase retrieval algorithms

View Table | View all tables in this article

3.1. Addressing a phase pattern to the SLM

3.1.1. Canceling the intrinsic aberrations

In order to be sure that only the aberration added on the SLM is measured, first all initial aberrations are eliminated. These aberrations can include alignment errors and the surface of the SLM itself, which might not be completely flat [6]. Removing the initial aberrations has been done by using the Shack-Hartmann sensor and retrieving the wavefront using the least-squares method (Method A).

3.1.2. Calibration of the SLM response

The way a phase pattern is addressed to the SLM is as follows. First, an aperture on the SLM is defined. In the current research this is a circular aperture, but the same method is valid if an annular aperture is used. All pixels i within this aperture are used. A “Zernike matrix” can be set up, such that for each pixel i within the aperture the value of all necessary Zernike circle polynomials can be computed. In matrix form this would be

Z = [\begin{matrix} Z_{1} |_{1} & Z_{2} |_{1} & \dots & Z_{J} |_{1} \\ Z_{1} |_{2} & Z_{2} |_{2} & \dots & Z_{J} |_{2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ Z_{1} |_{I} & Z_{2} |_{I} & \dots & Z_{J} |_{I} \end{matrix}],

where J is the total number of Zernike circle polynomials evaluated, and I is the total number of pixels within the aperture. Note that here only the real valued polynomials were added in the wavefront, since this is compatible with the experiment. Due to the cyclic nature of the phase pattern and the limits of the SLM, the phase difference assigned to the SLM should be between 0 and 2π. This can be done using the modulo (or mod) operation: a mod n = a − n ⌊a/n⌋.

If p is a vector containing the values of the individual pixels of the SLM, it can be constructed from the vector a containing the coefficients for the to-be-added aberration by

p = (Z a) \mod 2 π .

Therefore, in the case of updating the SLM pattern we use

p_{n e w} = (p_{o l d} - γ Z \hat{a}) \mod 2 π,

where γ is the gamma function of the SLM, calibrated for the used wavelength.

After the correction is done and the wavefront from the SLM is flattened, the phase pattern of the “flat” phase is saved. To this phase pattern, a phase ramp is added to separate the 0^th and 1^st order of the SLM together with the desired aberration, in the same fashion as the correction is added. After this, it is necessary to check if the maximum phase change over 4 pixels is not being exceeded. Such a check is necessary to avoid aliasing of the SLM-phase.

In the current research, this aliasing constraint is simplified to the constraint that the difference between two neighboring pixels should not exceed 0.5π. This is evaluated by letting p_i,_j be the value of the pixel located at position i, j on the SLM electrode matrix. Then, first two matrices are constructed:

\begin{array}{l} Δ p_{x} = p_{i, j + 1} - p_{i, j}, \\ Δ p_{y} = p_{i + 1, j} - p_{i, j} . \end{array}

Afterward, the element-wise minimum is taken between Δp and 2π − |Δp| for both x and y, in order to account for the modulated phase. If any of the values of this piecewise minimum is above 0.5π, it is said to break the aliasing constraint.

3.2. Constructing the matrix G

Recalling from the theory described in the previous section, in order to recover the coefficients describing the aberrated wavefront, s and G have to be constructed. Here, s relies on the two Hartmannograms, one of the flat wavefront and one of the aberrated wavefront from the SLM, and the radius r_SH of beam hitting the Shack-Hartmann sensor. Two G matrices will then be constructed for the respective Methods A and B.

3.2.1. Initial wavefront reconstruction pupil

To construct G, it is necessary to average the gradients of the Zernike circle polynomials for method A or the complex valued Zernike circle polynomials for method B. The window over which the polynomial has to be averaged can be seen as a scaled version of the lenslet, scaled so that all illuminated lenslets fit the unit disc. In order to compute the Zernike circle polynomials in these windows, the center position c and the radius r_SH of the beam on the Shack-Hartmann sensor have to be known. The window size is estimated by the average distance between the nearest neighbor spots on the Shack-Hartmann sensor. The center position c follows the average position of all the midpoints of the spots on the Shack-Hartmann sensor. The radius is estimated by calculating the length between the center and the furthest spot from the center.

3.2.2. Wavefront pupil boundaries optimization (with RMS optimization)

After this first estimation of the center position and radius is made, a better estimation can be found by optimization. For this optimization, an error has been defined that can be minimized. In this research, a root-mean-square, or RMS type error is chosen. Using the known added aberration and the measured aberration coefficients, an RMS error can be defined. To this purpose, a new Zernike matrix similar to Eq. (34) is constructed, this time with N points on a grid within the unit disc. The reference phase p_ref and the recovered phase p_rec can be constructed as

p = Z a,

where for p_ref, the reference vector a_ref is used, while for p_rec the estimated coefficient vector

\hat{a}

is used. The RMS error is then defined

ϵ = {‖ \frac{p_{ref} - p_{rec}}{N} ‖}_{2},

where ||x||₂ is the Euclidean vector norm of vector x. This RMS error is then minimized for center position and radius of the beam on the Hartmannogram. A limited memory bound Broyden-Fletcher-Goldfarb-Shanno (L-BFGS-B) minimization algorithm is applied to find the center position and radius [7]. The termination conditions for this optimization are:

\frac{f^{k} - f^{k + 1}}{\max {| f^{k} |, | f^{k + 1} |, 1}} \leq 10^{- 5},

\max {| proj (g_{i}) | i = 1, \dots, n} \leq 10^{- 5},

k \geq 10^{3},

where f^k is the value of the RMS error of the k^th iteration of the minimization algorithm, and proj(g_i) is i^th component of the projected gradient where n projections are made. If any of these statements were true, the optimization was terminated.

After looking at the RMS error landscapes, it was found that not every minimum found was a global minimum. If this was the case, the global minimum coordinates were estimated using the RMS landscape graphs, and a brute force optimization was run around those coordinates. This brute force optimization calculates the RMS error value in a grid of points. From the coordinates with the lowest RMS error, a new downhill simplex minimization algorithm is started. This way the global minimum was attempted to be found, and the optimal center and radius positions were determined.

These optimized parameters can then be used to find the RMS value fitting any degree of Zernike circle polynomials, and can also be used to determine the error landscape by calculating the RMS values when the center and radius differ slightly from the optimal value.

3.2.3. Application to the Zernike coefficients determination

Following the rule of thumb defined in [8], the amount of Zernike circle polynomials that can be fit given k spots on the Shack-Hartmann sensor is ^k/3. In the current research, there were about 144 spots on the Shack-Hartmann sensor, so, in this way, the optimum amount of polynomials was determined to be 45, i.e., Zernike expansion with a maximum degree of eight. In Figs. 5 and 6 other fitting orders are shown for some Zernike coefficients, and as one can see, this fitting order is good for all cases.

Fig. 2 Reference wavefronts (left column) for the first four experiments of Table 1. Reconstructed wavefronts from experimental data using Method A (middle column) and Method B (right column). The RMS errors under each row have been obtained by comparing the reconstructed (middle and right column) with reference wavefront (left column). For the reconstructions, the first 8 degrees in the Zernike expansion has been used.

Download Full Size | PDF

Fig. 3 Reference wavefronts (left column) for the second four experiments of Table 1. Reconstructed wavefronts from experimental data using Method A (middle column) and Method B (right column). The RMS errors under each row have been obtained by comparing the reconstructed (middle and right column) with reference wavefront (left column). For the reconstructions, the first 8 degrees in the Zernike expansion has been used.

Download Full Size | PDF

Fig. 4 Reference wavefronts (left column) for the last four experiments of Table 1. Reconstructed wavefronts from experimental data using Method A (middle column) and Method B (right column). The RMS errors under each row have been obtained by comparing the reconstructed (middle and right column) with reference wavefront (left column). For the reconstructions, the first 8 degrees in the Zernike expansion has been used.

Download Full Size | PDF

Fig. 5 Convergence of coefficients for the sub_zerns_1 experiment.

Download Full Size | PDF

Fig. 6 Convergence of coefficients for the sub_zerns_3 experiment.

Download Full Size | PDF

For method A, the amount of polynomials fit is equal to the amount of retrieved coefficients, while for Method B, if a fit is made with up to maximum degree of eight, the first seven degrees can be retrieved.

It should be noted that this optimization can take long due to the fact that the geometry matrix needs to be calculated in every iteration, as the values in the matrix depend on the center position and radius.

4. Results

The two methods have been implemented for three different test cases: (a) single Zernike aberrations, (b) a combination of three aberrations, assumed not to show crosstalk, and (c) four cases of aberrations where crosstalk is present. The aberrations are listed in Table 1. We remind that the values of the Zernike coefficients presented in the table are the nominal values that are added to the SLM, in addition to the intrinsic aberrations of the system. These nominal values have been used to calculate the reference wavefront in order to compare it with the estimated wavefront from experimental data using Method A and B. In Figs. 5 and 6 one can see the nominal values (dotted red line) and the values that are estimated using Method A (triangles) and Method B (dots) for several maximum fitted Zernike degree when some of the added Zernike coefficients (Experiment sub_zerns_1 and sub_zers_3).

Table 1. Coefficients Used in Specific Zernike Experiments, Rounded Off to 3 Significant Numbers

View Table | View all tables in this article

The reconstructions using the recovered coefficients are shown in Figs. 2–4. For these reconstructions, the first eight degrees in the Zernike expansion were used. The left column shows the reference wavefront that is constructed by using the nominal aberrations that have been added on the SLM, while the middle and right column show the wavefront reconstruction starting from the Shack Hartmann data and using our own implementation of the least-squares method (Method A) and Method B, respectively, as explained in the previous two sections. The RMS error, defined in Eq. (39) is the error between the reconstructed wavefronts using Method A and B (middle and right columns) and the reference (left column) is indicated below each row in the figures. From these results, it can be seen that the Method B compares very well with Method A and the input reference wavefront.

Based on the errors in the reconstruction found in Figs. 2, 3, and 4, one can see that both methods provide comparable accuracy in reconstructing the wavefront. However, it is known (and also observed here) that Method A presents crosstalk of coefficients when less Zernike powers are fit than there are aberrations present in the system. In the following, we show experimentally that this is not the case for Method B.

Using the same center position and radius, different amount of Zernike degrees can be fit in order to see the convergence behavior of gathered coefficients. Due to the fact that there are only 1, 2, or 3 Zernike modes present in the specific Zernike experiment, the convergence behavior of these experiments can be visualized.

Based on analysis of the aberrations listed in Table 1, the single Zernike experiments 5_1 and 6_4 and 3 random Zernike experiments 3_zerns_1 and 3_zerns_3, one can see in Figs. 2 and 3 that the coefficients seem to be measured independently from each other, and there is a little difference between the two methods. However, looking at the subsequent Zernike experiments sub_zerns_1 and sub_zerns_3, a large difference can be seen between the two reconstruction methods. The initial guesses of Method A over- or underestimate the presence of the aberration when not enough powers are fit. From Fig. 5, it can be seen that defocus is overestimated by more than 50% until 4 degrees are fit. Until a maximum degree of four, the spherical aberration coefficient $a_{4}^{0}$ is overestimated. All values seem to be within normal range when fitting a maximum degree of eight. The same can be seen in Fig. 6, where the coefficient $a_{3}^{3}$ is overestimated at maximum degree of three and four. Method B does not present this over-estimation. This behavior is also observed in the sets sub_zerns_2 and sub_zerns_4 (not shown in the figures). We also note that in the case of sub_zerns_3, both methods present some overestimation for large fitted degree. We don’t know the exact reason for that, but we observe that for this experiment, the RMS error between both methods and the reference wavefront was higher than for the case of sub_zerns_1, sub_zerns_2 and sub_zerns_4 (see Fig. 4).

5. Conclusion and outlook

In conclusion, we have implemented a new wavefront reconstruction method (Method B) for Shack-Hartmann sensors based on Zernike expansion of derivatives of the Zernike circle polynomials that was introduced in [1]. We have shown with experiments that Method B is advantageous as compared to other known methods such as a least-squares method (Method A) due to the lack of cross-talk between the coefficients. We have also compared Method B to Method A for reconstructing the wavefront using a Shack-Hartmann sensor for some sets of aberrations.

Based on the error values of the reconstruction with optimal amount of Zernike coefficients, the quality of the fit is in general similar for Method A and Method B. When less than the optimal Zernike powers were fit, it was seen that Method B estimates the coefficients more accurately. For the single Zernike aberration case, it was shown that only Method A shows cross-coupling of higher order aberrations, while Method B does not.

As to computational load, we may point out that these are of comparable order, with Method B somewhat more demanding due to the complex arithmetics and the alignment needed for compatibility with the ANSI-format of the Shack-Hartmann sensor.

Finally, we would like to point out that there are also alternative expansions to treat this problem, namely the use of eigenfunctions of the Laplacian with Neumann boundary conditions [9,10] that do have orthogonal gradients. It would be interesting as future work to compare this method with Method A and Method B considered in the present paper.

Appendix: Relation with orthogonal vector polynomials

In the course of our investigations, we became aware of a seemingly different approach to obtain the wavefront aberration coefficients from wavefront derivative data. We shall describe and relate this approach to our method in the framework and notations of [1], so that the Zernike circle polynomials are unnormalized and have exponential azimuthal dependence facilitating mathematical developments. The reader will have no particular problems in reformulating the main results of this appendix in terms of the ANSI-style circle polynomials using Eqs. (27)–(29). Thus, we have for integer n and m such that n − |m| is even and non-negative:

Z_{n}^{m} (ν, μ) \equiv Z_{n}^{m} (ρ, θ) = R_{n}^{| m |} (ρ) e^{i m θ},

with real ν, µ such that ν² + µ² ≤ 1 and

ν + i μ = ρ e^{i θ}; ν = ρ \cos θ; μ = ρ \sin θ,

and the radial polynomials

R_{n}^{| m |} (ρ)

given by Eq. (11). We now sketch the approach in [5], where the notations and conventions differ from our present analysis. The approach [5] uses the notion of vector polynomials

{\underline{G}}_{n}^{m} = {\underline{G}}_{n}^{m} (ρ, θ) = (G_{n, 1}^{m} (ρ, θ), G_{n, 2}^{m} (ρ, θ)) \in ℂ^{2}

for integer n, m such that n − |m| is even and non-negative with n ≠ 0, that satisfy

\frac{1}{π} \int_{0}^{1} \int_{0}^{2 π} \nabla Z_{n}^{m} (ρ, θ) . {\underline{G}}_{n^{'}}^{m^{'} *} (ρ, θ) ρ d ρ d θ = δ_{m, m^{'}} δ_{n, n^{'}}

for integer n, m, n′, m′ such that n − |m| and n′ − |m′| are even and non-negative while n ≠ 0 ≠ n′. In Eq. (46) we denote for (z, w) ∈ ℂ² and (g, h) ∈ ℂ²

(z . w) \cdot (g, h) = z g + w h,

the ∗ denotes complex conjugation, and

(\nabla Z_{n}^{m}) (ν, μ) = (\frac{\partial Z_{n}^{m}}{\partial ν} (ν, μ), \frac{\partial Z_{n}^{m}}{\partial μ} (ν, μ)) \in ℂ^{2} .

With

α_{n}^{m}

the Zernike coefficients of W (to be found) according to Eq. (13), where, due to orthogonality [see Eq. (12)],

α_{n}^{m} = \frac{n + 1}{π} \int_{0}^{1} \int_{0}^{2 π} W (ρ, θ) Z_{n}^{m *} (ρ, θ) ρ d ρ d θ,

we have

\nabla W = \sum_{n, m} α_{n}^{m} \nabla Z_{n}^{m},

and

α_{n}^{m} = \frac{1}{π} \int_{0}^{1} \int_{0}^{2 π} \nabla W (ρ, θ) {\underline{G}}_{n}^{m *} (ρ, θ) ρ d ρ d θ .

In [5], Section 2.1, the condition Eq. (46) is elaborated using Green’s theorem, where it should be noted that in [5] the attention is restricted to the cases n = n′, m = m′ in Eq. (46). It is shown that the problem of satisfying Eq. (46), with n = n′, m = m′ is solved when we can find

{\underline{G}}_{n}^{m}

such that on the pupil

\nabla \cdot {\underline{G}}_{n}^{m} = \frac{\partial G_{n, 1}^{m}}{\partial ν} + \frac{\partial G_{n, 2}^{m}}{\partial μ} = - (n + 1) Z_{n}^{m},

while the boundary condition

G_{n, 1}^{m} (ρ = 1, θ) \cos θ + G_{n, 2}^{m} (ρ = 1, θ) \sin θ = 0

should be satisfied at the rim ρ² = ν² + µ² = 1 of the pupil. Next, in [5], Section 2.2.1, the

{\underline{G}}_{n}^{m}

are required to be irrotational, meaning that the

\nabla \times {\underline{G}}_{n}^{m} = 0

, because of considerations of minimal noise propagation. This condition of irrotationality is satisfied when there is a scalar function

U_{n}^{m}

on the pupil such that

{\underline{G}}_{n}^{m} = \nabla U_{n}^{m} .

Combining Eq. (54) and Eq. (52), we see that we want

U_{n}^{m}

to satisfy

\nabla^{2} U_{n}^{m} = Δ U_{n}^{m} = - (n + 1) Z_{n}^{m}

while the boundary condition Eq. (53)

\frac{\partial U_{n}^{m}}{\partial ν} (ρ = 1, θ) \cos θ + \frac{\partial U_{n}^{m}}{\partial μ} (ρ = 1, θ) \sin θ = 0

should hold on the rim ρ = 1 of the pupil. In [5] Section 2.2.2–2.2.3, the

U_{n}^{m}

are found by writing the condition Eq. (55) in polar coordinates, with separate consideration of the cases m = 0 and m ≠ 0, and explicitly using the series representation of radial polynomials. This yields the two components of

G_{n}^{m}

in trigonometric polynomial form, for which it can be shown that the boundary condition Eq. (53) is satisfied as well. In [1], Section 4, it is shown that Eq. (55) has a solution

U_{n}^{m} = - [\frac{Z_{n + 2}^{m}}{4 (n + 2)} - \frac{(n + 1) Z_{n}^{m}}{2 n (n + 2)} + \frac{Z_{n - 2}^{m}}{4 n}] .

We shall verify below that, when n = |m| + 2, |m| + 4, …, this

U_{n}^{m}

also satisfies the boundary condition Eq. (56), and that a concise formula for

{\underline{G}}_{n}^{m} = \nabla U_{n}^{m}

in terms of Zernike circle polynomials results. For the case that n = |m|, we shall show that

U_{n}^{m} = - [\frac{Z_{n + 2}^{m}}{4 (n + 2)} - \frac{(3 n + 4) Z_{n}^{m}}{4 n (n + 2)}], n = | m |,

satisfies Eqs. (55) and (56), and we find a concise formula in terms of the Zernike circle polynomials for

{\underline{G}}_{| m |}^{m}

as well. This then can be used to show that

α_{n}^{m}

of Eq. (51) actually coincide with the LMS estimator found in [1], Section 3. The trigonometric/polynomial solution form of

{\underline{G}}_{n}^{m}

found in [5] involves the coefficients in the series representation of

R_{n}^{| m |} (ρ)

in Eq. (58), and these become awkward to use when the degree becomes large (n should be limited to ≤ 44 when using double precision). This problem is virtually absent when the representation of the

{\underline{G}}_{n}^{m}

in terms of the Zernike circle polynomials is used since there are nowadays several methods for reliably computing (the radial parts of the) Zernike circle polynomials of arbitrary large degree n and azimuthal order m, [11], [12].

We shall now show for m = 1, 2, … and n = m + 2, m + 4, …. that $U_{n}^{m}$ of 57 satisfies Eq. (56). Recalling the convention in [1] that any $Z_{n^{'}}^{m^{'}}$ with |m′| > n′ is set to 0, we have

(\frac{\partial}{\partial ν} \pm i \frac{\partial}{\partial μ}) Z_{n^{'}}^{m^{'}} (ν, μ) = 2 \sum_{l = 0}^{\frac{n^{'} - | m^{'} |}{2}} (n^{'} - 2 l) Z_{n^{'} - 1 - 2 l}^{m^{'} \pm 1},

see [1], Eq. (13). Using Eq. (59) with ν + iµ = e^iθ, noting that

Z_{n - 1 - 2 l}^{m \pm 1} (1, θ) = e^{i (m \pm 1) θ}, l = 0, 1, \dots \frac{1}{2} (n - | m | - 1),

while (as m > 0)

Z_{m - 1}^{m - 1} (1, θ) = e^{i (m - 1) θ}, Z_{m - 1}^{m + 1} (1, θ) = 0,

we have

(\frac{\partial}{\partial ν} + i \frac{\partial}{\partial μ}) Z_{n}^{m} = 2 \sum_{l = 0}^{\frac{n - m}{2}} (n - 2 l) e^{i (m + 1) θ} - 2 m e^{i (m + 1) θ},

(\frac{\partial}{\partial ν} - i \frac{\partial}{\partial μ}) Z_{n}^{m} = 2 \sum_{l = 0}^{\frac{n - m}{2}} (n - 2 l) e^{i (m - 1) θ} .

Adding and subtracting Eqs. (62) and (63) from one another then gives

\frac{\partial Z_{n}^{m}}{\partial ν} = 2 e^{i m θ} \cos θ \sum_{l = 0}^{\frac{n - m}{2}} (n - 2 l) - m e^{i (m + 1) θ},

\frac{\partial Z_{n}^{m}}{\partial μ} = 2 e^{i m θ} \sin θ \sum_{l = 0}^{\frac{n - m}{2}} (n - 2 l) - \frac{1}{i} m e^{i (m + 1) θ} .

We also observe that

\sum_{l = 0}^{\frac{n - m}{2}} (n - 2 l) = \frac{n + m}{2} (\frac{n + m}{2} + 1) .

Then, from Eqs. (57), (64) and (66), we get

\begin{matrix} \frac{\partial U_{n}^{m}}{\partial ν} (ρ = 1, θ) = \\ = - [\frac{1}{4 (n + 2)} \frac{\partial Z_{n}^{m}}{\partial ν} - \frac{n + 1}{2 n (n + 2)} \frac{\partial Z_{n}^{m}}{\partial ν} + \frac{1}{4 n} \frac{\partial Z_{n - 2}^{m}}{\partial ν}] \\ = - 2 D_{n}^{m} e^{i m θ} \cos θ + m E_{n} e^{i (m + 1) θ}, \end{matrix}

where

D_{n}^{m} = \frac{\frac{n + 2 + m}{2} (\frac{n + 2 - m}{2} + 1)}{4 (n + 2)} - \frac{(n + 1) \frac{n + m}{2} (\frac{n - m}{2} + 1)}{2 n (n + 2)} + \frac{\frac{n - 2 + m}{2} (\frac{n - 2 - m}{2} + 1)}{4 n} = 0,

and

E_{n} = \frac{1}{4 (n + 2)} - \frac{n + 1}{2 n (n + 2)} + \frac{1}{4 n} = 0 .

Hence,

(\partial U_{n}^{m} / \partial ν) (ρ, θ) = 0

for ρ = 1 and, similarly

(\partial U_{n}^{m} / \partial μ) (ρ, θ) = 0

for ρ = 1, and this implies that Eq. (56) holds. This handles the case that m = 1, 2, … and n = m + 2, m + 4, …. The case the that m = −1, −2, … and n = |m| + 2, |m| + 4, … follows from the case already handled by using that

Z_{n^{'}}^{m^{'}} (ν, μ) = {(Z_{n^{'}}^{- m^{'}} (ν, μ))}^{*} .

The case that m = 0 and n = 2, 4, … can be proved in the same way as the case m = 1, 2, … and n = m + 2, m + 4, …, where now in Eqs. (64) and (65) the terms involving eⁱ⁽^m⁺¹⁾^θ disappear.

We shall next show that for n = |m| + 2, |m| + 4, …

\frac{\partial U_{n}^{m}}{\partial ν} = - \frac{1}{4} [Z_{n + 1}^{m + 1} + Z_{n + 1}^{m - 1} - Z_{n - 1}^{m + 1} - Z_{n - 1}^{m - 1}],

\frac{\partial U_{n}^{m}}{\partial μ} = - \frac{1}{4 i} [Z_{n + 1}^{m + 1} - Z_{n + 1}^{m - 1} - Z_{n - 1}^{m + 1} + Z_{n - 1}^{m - 1}] .

Then Eqs. (71) and (72) yield the announced concise expression in terms of the Zernike circle polynomials for

{\underline{G}}_{n}^{m} = (\frac{\partial U_{n}^{m}}{\partial ν}, \frac{\partial U_{n}^{m}}{\partial μ}) .

Using Eqs. (59) and (57), we have that

\begin{matrix} \frac{\partial U_{n}^{m}}{\partial ν} = - [\frac{1}{4 (n + 2)} \sum_{l = 0}^{\frac{n + 2 - | m |}{2}} (n + 2 - 2 l) (Z_{n + 1 - 2 l}^{m + 1} + Z_{n + 1 - 2 l}^{m - 1}) \\ - \frac{n + 1}{2 n (n + 2)} \sum_{l = 0}^{\frac{n - | m |}{2}} (n - 2 l) (Z_{n - 1 - 2 l}^{m + 1} + Z_{n - 1 - 2 l}^{m - 1}) \\ + \frac{1}{4 n} \sum_{l = 0}^{\frac{n - 2 - | m |}{2}} (n - 2 - 2 l) (Z_{n - 3 - 2 l}^{m + 1} + Z_{n - 3 - 2 l}^{m - 1})] \end{matrix}

We now observe that the three series in the last member of Eq. (74) have the same terms, except that the second series omits the term l = 0 from the first series and the third series omits the terms with l = 0, 1 from the first series. Using Eq. (69), we see that the terms in the first series with l = 2, 3, …(n + 2 − |m|)/2 are canceled, and we get

\begin{matrix} \frac{\partial U_{n}^{m}}{\partial ν} = - [\frac{1}{4 (n + 2)} ((n + 2) (Z_{n + 1}^{m + 1} + Z_{n + 1}^{m - 1} + n (Z_{n - 1}^{m + 1} + Z_{n - 1}^{m - 1})) \\ - \frac{n + 1}{2 n (n + 2)} n (Z_{n - 1}^{m + 1} + Z_{n - 1}^{m - 1})] \\ = - [\frac{1}{4} (Z_{n + 1}^{m + 1} + Z_{n - 1}^{m + 1}) - \frac{1}{4} (Z_{n - 1}^{m + 1} + Z_{n - 1}^{m - 1})], \end{matrix}

and this is Eq. (71). In a similar fashion, we get Eq. (72).

We now consider the case that n = m = 1, 2, …. Then $Z_{n - 2}^{m} = 0$ and $Δ Z_{n}^{m} = 0$ , and so it follows from Eq. (64) that $Δ U = - (m + 1) Z_{m}^{m}$ holds for any U of the form

U = - [\frac{Z_{m + 2}^{m}}{4 (m + 2)} - C Z_{m}^{m}] .

We shall determine C such that Eq. (56) holds for n = m. Now

Z_{m}^{m} (ρ, θ) = ρ^{m} e^{i m θ}, Z_{m + 2}^{m} (ρ, θ) = ((m + 2) ρ^{m + 2} - (m + 1) ρ^{m}) e^{i m θ},

and so the condition Eq. (56) for n = m and θ = 0 yields

C = \frac{3 m + 4}{4 m (m + 2)} .

It can be shown that with this value of C the U of Eq. (76) satisfies Eq. (56) also for θ, 0. This handles the case that n = m = 1, 2, ….. The case n = m = 0 is non-existent, and the case withn = −m = 1, 2, … follows from the case already handled by complex conjugation. We now also compute for this U, using Eq. (59), and m > 0

\begin{matrix} \frac{\partial U}{\partial ν} = - [\frac{1}{4 (m + 2)} \frac{\partial Z_{m + 2}^{m}}{\partial ν} - \frac{3 m + 4}{4 m (m + 2)} \frac{\partial Z_{m}^{m}}{\partial ν}] \\ = - [\frac{1}{4} Z_{m + 1}^{m + 1} + \frac{1}{4} Z_{m + 1}^{m - 1} + \frac{m}{4 (m + 2)} Z_{m - 1}^{m - 1} - \frac{3 m + 4}{4 m (m + 2)} m Z_{m - 1}^{m - 1}] \\ = - [\frac{1}{4} Z_{m + 1}^{m + 1} + \frac{1}{4} Z_{m + 1}^{m - 1} - \frac{1}{2} Z_{m - 1}^{m - 1}], \end{matrix}

and similarly

\frac{\partial U}{\partial μ} = - \frac{1}{i} [\frac{1}{4} Z_{m + 1}^{m + 1} - \frac{1}{4} Z_{m + 1}^{m - 1} + \frac{1}{2} Z_{m - 1}^{m - 1}] .

The results Eqs. (79) and (80) continue to hold when m < 0 and n = |m|. Hence, also in this case we get a concise result for

{\underline{G}}_{| m |}^{m}

as in Eqs. (71)–(73).

We finally show that the $α_{n}^{m}$ obtained in [1], Section 3 coincide with the $α_{n}^{m}$ of Eq. (51). Observe that the vector polynomials ${\underline{G}}_{m}^{m}$ in Eq. (51) were derived in [5] from the conditions Eq. (46) and Eq. (53) under the assumption that they have vanishing curls, whereas the LMS estimator found in [1], section 3 has been derived under the condition that a natural mean-square error functional involving expansion coefficients is minimized. This LMS estimate uses the expansion coefficients ${(β_{\pm})}_{n}^{m}$ in

\frac{\partial W}{\partial ν} \pm i \frac{\partial W}{\partial μ} = \sum_{n, m} β_{\pm}_{n}^{m} Z_{n}^{m} .

That is, we have for m ≠ 0

α_{n}^{m} = C_{n}^{m} φ_{n}^{m} - C_{n + 2}^{m} φ_{n + 2}^{m}, n = | m |, | m | + 2, \dots,

with

φ_{n}^{m} = \frac{1}{2} {(β_{+})}_{n + 1}^{m + 1} + {(β_{-})}_{n - 1}^{m - 1},

C_{| m |}^{m} = \frac{1}{| m |}, C_{n}^{m} = \frac{1}{2 n}, n = | m | + 2, | m | + 4, \dots .

(when m = 0 we only consider n = 2, 4, …). In the case that n = |m| + 2, |m| + 4, …, we see that Eq. (82) and (83) gives

α_{n}^{m} = \frac{1}{4 n} {(β_{+})}_{n - 1}^{m + 1} + \frac{1}{4 n} {(β_{-})}_{n - 1}^{m - 1} - \frac{1}{4 (n + 2)} {(β_{+})}_{n + 1}^{m + 1} - \frac{1}{4 (n + 2)} {(β_{-})}_{n + 1}^{m - 1} .

Now we have from Eq. (81)

\frac{\partial W}{\partial ν} = \sum_{n, m} \frac{1}{2} ({(β_{+})}_{n}^{m} + {(β_{-})}_{n}^{m}) Z_{n}^{m}

\frac{\partial W}{\partial μ} = \sum_{n, m} \frac{1}{2 i} ({(β_{+})}_{n}^{m} - {(β_{-})}_{n}^{m}) Z_{n}^{m} .

Therefore, from Eq. (71) and (72) and the orthogonality/normalization of the

Z_{n}^{m}

, see Eq. (12), we see that the integral expression at the right-hand side of Eq. (51) becomes

\begin{matrix} \frac{1}{π} \int_{0}^{1} \int_{0}^{2 π} (\frac{\partial W}{\partial v} \frac{\partial U_{n}^{m *}}{\partial v} + \frac{\partial W}{\partial μ} \frac{\partial U_{n}^{m *}}{\partial μ}) ρ d ρ d θ \\ = - \frac{1}{8 π} \int_{0}^{1} \int_{0}^{2 π} \sum_{n^{'}, m^{'}} ({(β_{+})}_{n^{'}}^{m^{'}} + {(β_{-})}_{n^{'}}^{m^{'}}) Z_{n^{'}}^{m^{'}} \\ \times {[Z_{n + 1}^{m + 1} + Z_{n + 1}^{m - 1} - Z_{n - 1}^{m + 1} - Z_{n + 1}^{m - 1}]}^{*} ρ d ρ d θ \\ = - \frac{1}{8 π} \int_{0}^{1} \int_{0}^{2 π} \sum_{n^{'}, m^{'}} ({(β_{+})}_{n^{'}}^{m^{'}} - {(β_{-})}_{n^{'}}^{m^{'}}) Z_{n^{'}}^{m^{'}} \\ \times {[Z_{n + 1}^{m + 1} - Z_{n + 1}^{m - 1} - Z_{n - 1}^{m + 1} + Z_{n + 1}^{m - 1}]}^{*} ρ d ρ d θ \\ = - \frac{1}{8 (n + 2)} ({(β_{+})}_{n + 1}^{m + 1} + {(β_{-})}_{n + 1}^{m + 1} + {(β_{+})}_{n + 1}^{m - 1} + {(β_{-})}_{n + 1}^{m - 1}) \\ + \frac{1}{8 n} ({(β_{+})}_{n - 1}^{m + 1} + {(β_{-})}_{n - 1}^{m + 1} + {(β_{+})}_{n - 1}^{m - 1} + {(β_{-})}_{n - 1}^{m - 1}) \\ - \frac{1}{8 (n + 2)} ({(β_{+})}_{n + 1}^{m + 1} - {(β_{-})}_{n + 1}^{m + 1} - {(β_{+})}_{n + 1}^{m - 1} + {(β_{-})}_{n + 1}^{m - 1}) \\ + \frac{1}{8 n} ({(β_{+})}_{n - 1}^{m + 1} - {(β_{-})}_{n - 1}^{m + 1} - {(β_{+})}_{n - 1}^{m - 1} + {(β_{-})}_{n - 1}^{m - 1}), \end{matrix}

and this coincides with the right-hand side of Eq. (86) when the various cancellations are noted.

In a similar fashion the case n = |m| > 0 can be handled using Eqs. (79) and (80), and this yields that the $α_{| m |}^{m}$ from Eqs. (82) and (83) coincide with the right-hand side integral expression in Eq. (51).

Funding

European Research Council (ERC), FP7/2007-2013, 339681.

Acknowledgments

The authors would like to acknowledge Prof. J. Braat for useful discussions and comments on this article. We also thank Dr. S. van Haver for discussions.

References

1. A. J. E. M. Janssen, “Zernike expansion of derivatives and laplacians of the zernike circle polynomials,” J. Opt. Soc. Am. A 31, 1604–1613 (2014). [CrossRef]

2. The choice for ANSI convention is motivated here by how the Shack-Hartmann sensor deals with the circle polynomials. We would like to point out that using the Born and Wolf convention (that uses exponential instead of cosine/sine functions), renders the results in [1] more transparent and concise.

3. L. N. Thibos, R. A. Applegate, J. T. Schwiegerling, and R. Webb, “Standards for reporting the optical aberrations of eyes,” J. Refract. Surg. 18, S652–S660 (2002). [PubMed]

4. G. M. Dai, Wavefront Optics for Vision Correction, vol. 179 (SPIE Press, 2008). [CrossRef]

5. V. N. Mahajan and E. Acosta, “Vector polynomials for direct analysis of circular wavefront slope data,” J. Opt. Soc. Am. 34, 1908–1913 (2017). [CrossRef]

6. N. Matsumoto, T. Ando, T. Inoue, Y. Ohtake, N. Fukuchi, and T. Hara, “Generation of high-quality higher-order Laguerre-Gaussian beams using liquid-crystal-on-silicon spatial light modulators,” J. Opt. Soc. Am. A 25, 1642–1651 (2008). [CrossRef]

7. R. Fletcher, Practical Methods of Optimization (John Wiley & Sons, 1987).

8. O. Soloviev and G. Vdovin, “Estimation of the total error of modal wavefront reconstruction with zernike polynomials and hartmann-shack test,” in 5th International Workshop on Adaptive Optics for Industry and Medicine, vol. 6018 (International Society for Optics and Photonics, 2006), p. 60181D.

9. H. Shengyang, X. Fengjie, L. Changhai, and J. Zongfu, “Eigenfunctions of Laplacian for phase estimation from wavefront gradient or curvature sensing,” Opt. Commun. 284, 2781–2783 (2011). [CrossRef]

10. H. Shengyang, N. Yu, X. Fengjie, and J. Zongfu, “Modal wavefront reconstruction with Zernike polynomials and eigenfunctions of the Laplacian,” Opt. Commun. 288, 7–12 (2013). [CrossRef]

11. A. J. E. M. Janssen and P. Dirksen, “Computing zernike polynomials of arbitrary degree using the discrete fourier transform,” J. Eur. Soc. Rapid Publ. 2, 07012 (2007). [CrossRef]

12. B. H. Shakibaei and R. Paramesran, “Recursive formula to compute Zernike radial polynomials,” Opt. Lett. 14, 2487–2489 (2013). [CrossRef]

Experiment name	Values of the applied Zernike coefficients (Units of waves)
5_1	$a_{5}^{1}$ 0.750
5_5	$a_{5}^{5}$ 1.250
6_2	$a_{6}^{2}$ 0.500
6_4	$a_{6}^{4}$ 1.000
6_6	$a_{6}^{6}$ 2.000
3_zerns_1	$a_{2}^{0}$ 2.000	$a_{4}^{4}$ 1.500	$a_{6}^{- 2}$ 0.500
3_zerns_2	$a_{2}^{0}$ 2.000	$a_{3}^{- 1}$ 3.000	$a_{5}^{5}$ 1.500
3_zerns_3	$a_{2}^{- 2}$ 2.000	$a_{4}^{0}$ 1.500	$a_{6}^{- 4}$ 0.750
sub_zerns_1	$a_{2}^{0}$ 4.000	$a_{4}^{0}$ 1.500	$a_{6}^{0}$ 0.750
sub_zerns_2	$a_{1}^{1}$ 4.000	$a_{3}^{1}$ 1.500	$a_{5}^{1}$ 0.750
sub_zerns_3	$a_{3}^{3}$ 2.500	$a_{5}^{3}$ 0.750
sub_zerns_4	$a_{2}^{2}$ 2.500	$a_{4}^{2}$ 1.000	$a_{6}^{1}$ 0.500

Experiment name	Values of the applied Zernike coefficients (Units of waves)
5_1	$a_{5}^{1}$ 0.750
5_5	$a_{5}^{5}$ 1.250
6_2	$a_{6}^{2}$ 0.500
6_4	$a_{6}^{4}$ 1.000
6_6	$a_{6}^{6}$ 2.000
3_zerns_1	$a_{2}^{0}$ 2.000	$a_{4}^{4}$ 1.500	$a_{6}^{- 2}$ 0.500
3_zerns_2	$a_{2}^{0}$ 2.000	$a_{3}^{- 1}$ 3.000	$a_{5}^{5}$ 1.500
3_zerns_3	$a_{2}^{- 2}$ 2.000	$a_{4}^{0}$ 1.500	$a_{6}^{- 4}$ 0.750
sub_zerns_1	$a_{2}^{0}$ 4.000	$a_{4}^{0}$ 1.500	$a_{6}^{0}$ 0.750
sub_zerns_2	$a_{1}^{1}$ 4.000	$a_{3}^{1}$ 1.500	$a_{5}^{1}$ 0.750
sub_zerns_3	$a_{3}^{3}$ 2.500	$a_{5}^{3}$ 0.750
sub_zerns_4	$a_{2}^{2}$ 2.500	$a_{4}^{2}$ 1.000	$a_{6}^{1}$ 0.500

Implementation and benchmarking of a crosstalk-free method for wavefront Zernike coefficients reconstruction using Shack-Hartmann sensor data

Abstract

1. Introduction

2. Theory

2.1. Zernike circle polynomials as a basis for wavefront expansion

2.2. Complex-valued definition of the Zernike circle polynomials

2.3. Data from the Shack-Hartmann sensor

2.3.1. Method A

2.3.2. Method B

2.3.3. Wavefront reconstruction using Method A and Method B

3. Experimental setup and methods

3.1. Addressing a phase pattern to the SLM

3.1.1. Canceling the intrinsic aberrations

3.1.2. Calibration of the SLM response

3.2. Constructing the matrix G

3.2.1. Initial wavefront reconstruction pupil

3.2.2. Wavefront pupil boundaries optimization (with RMS optimization)

3.2.3. Application to the Zernike coefficients determination

4. Results

5. Conclusion and outlook

Appendix: Relation with orthogonal vector polynomials

Funding

Acknowledgments

References

Cited By

Figures (6)

Tables (2)

Equations (88)

OSA Continuum