Highly accurate wavefront reconstruction algorithms over broad spatial-frequency bandwidth

Seung-Whan Bahk

doi:10.1364/OE.19.018997

1. Introduction

Frequency-domain wavefront reconstruction methods are as old as the very early wavefront reconstructors [1, 2]. Freischlad placed this subject on solid ground [3]. The rectangular map constraint of the conventional Fourier method has been removed in an iterative Gerchberg-type algorithm dealing with an arbitrary boundary shape [4]. A series of recent papers by Poyneer discusses improvements on handling boundary conditions and applications in extreme adaptive optics [5, 6]. Similar principles have been applied in shearing interferometers [7]. More serious attention has been paid to the accuracy of the reconstruction methods in Refs. [8–10]. The works of Campos and Yaroslavsky presented a one-dimensional solution based on band-limited integration technique in frequency domain. The two-dimensional (2-D) extension of the same method was not discussed. Complementary to their works, Bahk has introduced a full two-dimensional wavefront reconstructor based on the band-limited derivative calculation [11]. Both approaches emphasize the frequency response of the reconstructed signals. The frequency response of wavefront reconstruction has been discussed earlier in the analysis of lateral shearing interferometry [12]. Frequency response characteristic of a reconstructor is important in focal spot diagnostic for high power lasers where the focal spot is indirectly characterized using wavefront information reconstructed from Shack-Hartmann slopes data [13]. The error in the high-frequency components of the reconstructed wavefront results in an error in the estimate of encircled energy.

This paper develops a set of encompassing mathematical tools for wavefront reconstruction problems, where many additional benefits naturally arise interconnecting the results of previous works. The benefits are exemplified by the development of two new wavefront reconstructors which exhibit high accuracy over broad spatial-frequency bandwidth. It also provides a unified analytic expression for the noise propagation coefficients of several well-known traditional wavefront reconstructors as well as the ones developed here.

This paper consists of four parts. Section 2 introduces the mathematical tools and symbols regarding band-limited derivative operations, which are needed for the analyses in the subsequent sections. In Sec. 3, a way to improve the accuracy of the finite difference method is discussed in connection with wavefront reconstruction. The Simpson rule is adopted for developing a new spatial-domain iterative reconstruction algorithm. The exact details of the algorithm and its frequency domain property are described. In Sec. 4, the band-limited reconstruction algorithm developed for rectangular geometry [11] is extended to hexagonal geometry, which greatly enhances the flexibility of this scheme. Finally, in Sec. 5 the noise propagation curve is analytically derived and compared with numerical simulations.

2. Band-limited derivative

The main results of band-limited derivative techniques in the context of wavefront reconstruction were summarized in [11]. The full derivation of the results will be presented here for the sake of completeness. Additional new notations are introduced which simplify the expressions in Sec. 4.

The motivation for band-limited derivatives, especially for discrete samples, lies in the fact that it provides analytical tool for converting back and forth between slopes measurements and wavefront signal. We start by asking what the exact interpolation formula is for derivatives in discrete samples. According to sampling theorem [14], a band-limited signal can be exactly reconstructed at any point by convolving a sinc-function with discrete samples.

φ (x) = \sum_{n = - \infty}^{\infty} φ (n Δ x) sinc [2 W (x - n Δ x)],

where 2W = 1/Δx and sinc(x) = sin(π x)/(π x).

Using the result

\frac{d}{d x} sinc (x) = - π j_{1} (π x),

where j ₁ is a spherical Bessel function of the first kind, the derivative of function φ(x) is

\frac{d φ (x)}{d x} = - \frac{1}{Δ x} \sum_{n = - \infty}^{\infty} φ (n Δ x) π j_{1} [2 π W (x - n Δ x)] .

Using the fact that j ₁ is an odd function, the derivative interpolation expression at discrete points is

\frac{d φ}{d x} | (x = m Δ x) \equiv φ_{x} (m) = \frac{1}{Δ x} \sum_{n = 1}^{\infty} π j_{1} (n π) {φ [(m + n) Δ x] - φ [(m - n) Δ x]},

where the spherical Bessel function evaluated at integer multiples of π is equivalent to

π j_{1} (n π) = \frac{{(- 1)}^{n + 1}}{n}, n = 1, 2, \dots

The summation of the left hand side of Eq. (4) for all sample points can be shown to be equal to zero by taking advantage of the right hand side expression and using the periodicity condition of discrete samples:

\sum_{m = 0}^{N - 1} \frac{d φ}{d x} | (x = m Δ x) = 0 .

Equation (4) is easier to handle in frequency domain. Discrete-Fourier-transform (DFT) of Eq. (4) is

{\tilde{φ}}_{x} (k) = \frac{2 i}{Δ x} [\sum_{n = 1}^{\infty} \frac{{(- 1)}^{n + 1}}{n} sin (\frac{2 π n}{N} k)] \tilde{φ} (k),

where the tilde notation means DFT of the symbol beneath it.

The coefficient of $\tilde{φ}$ (k) on the right hand side is a Fourier series which represents a sawtooth wave. Thus Eq. (7) is further simplified as

{\tilde{φ}}_{x} (k) = i \frac{2 π}{Δ x} 𝖲 (k) \tilde{φ} (k),

where 𝖲(k) (sawtooth wave) is defined as

𝖲 (k) = {\begin{array}{l} for even N & {\begin{array}{l} k / N, & k = 0, \dots, N / 2 - 1 \\ 0, & k = N / 2 \\ k / N - 1, & k = N / 2 + 1, \dots, N - 1 \end{array} \\ for odd N & {\begin{array}{l} k / N, & k = 0, \dots, (N - 1) / 2 \\ k / N - 1, & k = (N + 1) / 2, \dots, N - 1 . \end{array} \end{array}

Equation (8) provides a convenient way of calculating exact derivatives from band-limited signals.

For the sampling points of a derivative signal offset by a half sampling space from the sampling points of the original signal, a slightly different form is derived going through similar steps as in Eqs. (1)–(4),

\frac{d φ}{d x} | (x = m Δ x + \frac{1}{2} Δ x) \equiv φ_{x, \frac{1}{2}} (m) = \frac{1}{Δ x} \sum_{n = 1}^{\infty} π j_{1} (n π - \frac{1}{2} π) {φ [(m + n) Δ x] - φ [(m - n + 1) Δ x]},

where the spherical-Bessel coefficients can be replaced again with an integer expression

π j_{1} (n π - \frac{1}{2} π) = \frac{4}{π} \frac{{(- 1)}^{n + 1}}{{(2 n - 1)}^{2}}, n = 1, 2, \dots

Employing similar Fourier series analysis that leads to Eq. (8), the frequency-domain expression of Eq. (10) is reduced to

{\tilde{φ}}_{x, \frac{1}{2}} (k) = i \frac{2 π}{Δ x} exp (i \frac{π}{N} k) 𝖳 (k) \tilde{φ} (k)

where 𝖳(k) (triangular wave) is defined as

𝖳 (k) = {\begin{array}{l} k / N, & k = 0, \dots, ⌊ \frac{N - 1}{2} ⌋ \\ 1 - k / N, & k = ⌊ \frac{N - 1}{2} ⌋ + 1, \dots, N - 1 . \end{array}

N can be either an odd or even integer as the floor notation in the above equation will apply to either type.

We also need an interpolation formula for creating a signal shifted by half sample spacing for the Fried geometry,

φ (x = m Δ x + \frac{1}{2} Δ x) \equiv φ_{\frac{1}{2}} (m) = \sum_{n = 1}^{\infty} sinc (n - \frac{1}{2}) {φ [(m + n) Δ x] - φ [(m - n + 1) Δ x]},

The DFT of Eq. (14) is

{\tilde{φ}}_{\frac{1}{2}} (k) = exp (i \frac{π}{N} k) 𝖱 (k) \tilde{φ} (k)

where 𝖱(k) (rectangular wave) is

𝖱 (k) = {\begin{array}{l} for even N & {\begin{array}{l} 1, & k = 0, \dots, N / 2 - 1 \\ 0, & k = N / 2 \\ - 1, & k = N / 2 + 1, \dots, N - 1 \end{array} \\ for odd N & {\begin{array}{l} 1, & k = 0, \dots, (N - 1) / 2 \\ - 1, & k = (N + 1) / 2, \dots, N - 1 . \end{array} \end{array}

Thus the partial derivative in x-direction for Fried geometry in frequency domain is

{\tilde{φ}}_{x, \frac{1}{2}, \frac{1}{2}} (p, q) = i \frac{2 π}{Δ x} exp [i \frac{π}{N} (p + q)] 𝖳 (q) 𝖱 (p) \tilde{φ} (p, q)

Equation (17) has an additional degree of freedom (index p for y-direction) because the reconstructed sample point in the Fried geometry has to be first shifted in y-direction by half-sample size before applying the half-sample shifted derivative operation in x-direction. For consistency, we can verify that the sequential operations of half-pixel shift and derivative operation using 𝖱(k) and 𝖲(k) produce the same result as single operation of 𝖳(k), that is, 𝖳(k) = 𝖱(k) 𝖲(k). This relation, however, does not hold for the value at N/2 for even N, where the left hand side is 0.5 whereas the right hand side is 0. To remove this paradox for even number of samples, we choose to use 𝖲(N/2) = 0.5 and 𝖱(N/2) = 1 or 𝖲(N/2) = −0.5 and 𝖱(N/2) = −1. A similar choice was made in [9] in defining band-limited integration operators to prevent information loss at the boundary. The redefinition of 𝖲 and 𝖱 at the mid-point is implied hereafter. Using the new definition, the half-pixel operator used in the right hand side of Eq. (15) can be alternatively expressed as

exp (\pm \frac{π i}{N} k) 𝖱 (k) = exp [\pm π i 𝖲 (k)] .

𝖲 can be considered as a discrete angular frequency vector circularly shifted by ⌊N/2⌋ or ⌊N/2⌋ – 1 using the new definition at the mid-point. For an odd number of samples,

\frac{2 π}{Δ x} 𝖲 (p) = ↺ {\frac{2 π}{L} [- \frac{N - 1}{2} + p], \frac{N - 1}{2}},

where ↺ {A, N′} means circular shift of vector A to the left by N′. L = NΔx and p = 0,...,N – 1.

For an even number of samples with 𝖲(N/2) = −0.5,

\frac{2 π}{Δ x} 𝖲 (p) = ↺ {\frac{2 π}{L} (- \frac{N}{2} + p), \frac{N}{2}},

and with 𝖲(N/2) = 0.5,

\frac{2 π}{Δ x} 𝖲 (p) = ↺ {\frac{2 π}{L} [- (\frac{N}{2} - 1) + p], \frac{N}{2} - 1} .

We can collectively define the right hand side of Eqs. (19)-(21) as $\bar{k_{x}}$ where the bar over k_x implies a circular shift of the discrete angular frequency vector centered around zero by ⌊N/2⌋ or ⌊N/2⌋ – 1 as prescribed above. Using the new notation and Eq. (18) as well as the relation 𝖳(k) = 𝖱(k)𝖲(k), Eqs (8), (12), and (15) can be alternatively expressed as

{\tilde{φ}}_{x} = i \bar{k_{x}} \tilde{φ},

{\tilde{φ}}_{x, \frac{1}{2}} = i \bar{k_{x}} exp (i \frac{Δ x}{2} \bar{k_{x}}) \tilde{φ},

{\tilde{φ}}_{\frac{1}{2}} = exp (i \frac{Δ x}{2} \bar{k_{x}}) \tilde{φ} .

The $\bar{k_{x}}$ notation allows one to see the formal similarity with continuous variable derivative expressions. The continuous variable counterpart for Eq. (22), for example, is $FT (\frac{d φ}{d x}) = {i k}_{x} FT (φ)$ (where FT denotes Fourier transform and k_x here is a continuous spatial frequency variable. This connection allows one to derive other derivative formulas for discrete band-limited samples by taking advantage of continuous space results.

In many practical situations the band-limited calculations may not produce exact results depending on the nature of signals. The magnitudes of the Fourier coefficients of a linear function, for example, decrease as 1/(spatial frequency) whereas Eq. (22) indicates that the coefficients of the derivative is multiplied by spatial frequency vector. Therefore the highest spatial frequency coefficient does not vanish even if N approaches ∞. The linear terms thus are not band-limited and need to be treated separately using a polynomial fit.

Equations (22)–(24) form the basis of the following analysis.

3. Simpson reconstructor

The analysis in Sec. 2 suggests that the accurate derivative calculation at discrete samples requires the superposition sum of the whole available samples. This can be done more conveniently in the frequency domain, which results in a band-limited reconstructor with unity frequency response [11]. On the other hand, it is still worth while to investigate an improved finite difference scheme for purely spatial domain operation. In finite difference methods involving only a few points, a high degree of accuracy can be obtained by distributing the finite difference over both the measured derivative samples (i.e. slopes) and the integrated samples. Denoting the wavefront estimate as $\hat{φ}$ and the measured slope as S, one can start from a general finite difference expression such as

\sum_{j} a (j) \hat{φ} (i + j) = \sum_{k} b (k) S (i + k),

for one-dimensional problems. a and b are coefficients belonging to specific finite difference scheme. For example, Southwell [15] showed a reconstructor based on

\hat{φ} (i + 1) - \hat{φ} (i) = \frac{Δ x}{2} [S (i) + S (i + 1)] .

The frequency response of Southwell reconstructor is low at high spatial frequency [11]. We enhance the frequency response using the Simpson rule which is

\hat{φ} (i + 1) - \hat{φ} (i - 1) = \frac{Δ x}{3} [S (i - 1) + 4 S (i) + S (i + 1)] .

An iterative wavefront reconstructor based on this scheme will be developed in the next section. Other finite schemes which involve more than three samples (slopes) or have asymmetric forms were not considered at this point as the iteration scheme might become unnecessarily complicated.

3.1. Simpson iterator

Casting the local 1-D equation (27) into a least squares form in 2-D, we obtain an error metric (ɛ) as

\begin{array}{l} ɛ = \sum_{i, j} {\frac{1}{2 Δ x} [\hat{φ} (i, j + 1) - \hat{φ} (i, j - 1)] - \frac{1}{6} [S_{x} (i, j + 1) + 4 S_{x} (i, j) + S_{x} (i, j - 1)]}^{2} + \\ {\frac{1}{2 Δ y} [\hat{φ} (i + 1, j) - \hat{φ} (i - 1, j)] - \frac{1}{6} [S_{y} (i + 1, j) + 4 S_{y} (i, j) + S_{y} (i - 1, j)]}^{2} . \end{array}

Δx and Δy are moved around to the $\hat{φ}$ side so that squared terms are in units of slopes. This provides equal weight to the differences in x and y directions on the assumption that the magnitude of slopes is comparable in either direction.

The condition ∂ɛ/∂ $\hat{φ}$ (i, j) = 0 leads to an equation that can be used for an iterative algorithm. At the stationary point, the error is not decreasing anymore which indicates that the phase estimate $\hat{φ}$ has reached a solution. It is assumed that the boundary of phase and slope points can be an arbitrary shape. The differentiation of the error metric results in four groups, which are indicated by different colors in Fig. 1. Each group can be used in the equation only if all its elements exist. This strategy is realized by using g-parameters as shown in the following:

SIMPSON:

\begin{matrix} g_{L} [\hat{φ} (i, j) - \hat{φ} (i, j - 2)] + g_{R} [\hat{φ} (i, j) - \hat{φ} (i, j + 2)] \\ + g_{U} {(\frac{Δ x}{Δ y})}^{2} [\hat{φ} (i, j) - \hat{φ} (i - 2, j)] + g_{D} {(\frac{Δ x}{Δ y})}^{2} [\hat{φ} (i, j) - \hat{φ} (i + 2, j)] \\ = g_{L} [S_{x} (i, j - 2) + 4 S_{x} (i, j - 1) + S_{x} (i, j)] \frac{Δ x}{3} \\ - g_{R} [S_{x} (i, j) + 4 S_{x} (i, j + 1) + S_{x} (i, j + 2)] \frac{Δ x}{3} \\ + g_{U} {(\frac{Δ x}{Δ y})}^{2} [S_{y} (i - 2, j) + 4 S_{y} (i - 1, j) + S_{y} (i, j)] \frac{Δ y}{3} \\ - g_{D} {(\frac{Δ x}{Δ y})}^{2} [S_{y} (i, j) + 4 S_{y} (i + 1, j) + S_{y} (i + 2, j)] \frac{Δ y}{3} \\ \equiv Δ S (i, j) \end{matrix}

Fig. 1 Simpson geometry.

Download Full Size | PDF

g_L, g_R, g_U, g_D are flags with values 0 or 1 where L, R, U, D indicate left, right, up, and down directions, respectively. They are zeros if the slopes quantity in the parenthesis next to them are incalculable or ones otherwise. g_L at the point (i, j), for example, does not vanish only if all the slopes measurements exist at the additional points (i, j – 1) and (i, j – 2). The scope of each flag is graphically indicated in Fig. 1. For a comparison, the iterative equation for the Southwell reconstructor is written here using the same format.

SOUTHWELL:

\begin{matrix} g_{L} [\hat{φ} (i, j) - \hat{φ} (i, j - 1)] + g_{R} [\hat{φ} (i, j) - \hat{φ} (i, j + 1)] \\ + g_{U} {(\frac{Δ x}{Δ y})}^{2} [\hat{φ} (i, j) - \hat{φ} (i - 1, j)] + g_{D} {(\frac{Δ x}{Δ y})}^{2} [\hat{φ} (i, j) - \hat{φ} (i + 1, j)] \\ = g_{L} [S_{x} (i, j - 1) + S_{x} (i, j)] \frac{Δ x}{2} - g_{R} [S_{x} (i, j) + S_{x} (i, j + 1)] \frac{Δ x}{2} \\ + g_{U} {(\frac{Δ x}{Δ y})}^{2} [S_{y} (i - 1, j) + S_{y} (i, j)] \frac{Δ y}{2} - g_{D} {(\frac{Δ x}{Δ y})}^{2} [S_{y} (i, j) + S_{y} (i + 1, j)] \frac{Δ y}{2} \end{matrix}

The successive-over-relaxation technique in [15] can be applied to the Simpson iterative reconstructor:

{\hat{φ}}^{(m + 1)} (i, j) = {\hat{φ}}^{(m)} (i, j) + ω [\bar{{\hat{φ}}^{(m)}} (i, j) - {\hat{φ}}^{(m)} (i, j)],

where

\bar{{\hat{φ}}^{(m)}} (i, j) = [\bar{{\hat{φ}}_{0}^{(m)}} (i, j) + Δ S (i, j)] / (g_{L} + g_{R} + {g'}_{U} + {g'}_{D}),

\bar{{\hat{φ}}_{0}^{(m)}} (i, j) = g_{L} {\hat{φ}}^{(m)} (i, j - 2) + g_{R} {\hat{φ}}^{(m)} (i, j + 2) + {g'}_{U} {\hat{φ}}^{(m)} (i - 2, j) + {g'}_{D} {\hat{φ}}^{(m)} (i + 2, j) .

ω is the over-relaxation parameter. g′ is

{(\frac{Δ x}{Δ y})}^{2} \cdot g

.

3.2. Frequency response and regularization

The frequency response of the Simpson reconstructor will be calculated using a frequency-domain least-squares method. The sum of the squared error in spatial domain in Eq. (28) is equivalent to the sum of the squared error of the Fourier transformed component by the Parseval theorem:

ɛ = \frac{1}{N^{2}} Σ {{| D_{x} \tilde{\hat{φ}} - A_{x} {\tilde{S}}_{x} |}^{2} + {| D_{y} \tilde{\hat{φ}} - A_{y} {\tilde{S}}_{y} |}^{2}},

where

D_{x} = \frac{1}{2 Δ x} [exp (i \bar{k_{x}} Δ x) - exp (- i \bar{k_{x}} Δ x)] \equiv D_{x, Simpson}

and

A_{x} = \frac{1}{6} [exp (i \bar{k_{x}} Δ x) + 4 + exp (- i \bar{k_{x}} Δ x)] \equiv A_{x, Simpson} .

The solution for $\tilde{\hat{φ}}$ in Eq. (34) is

\tilde{\hat{φ}} = \frac{(D_{x}^{*} A_{x} {\tilde{S}}_{x} + D_{y}^{*} A_{y} {\tilde{S}}_{y})}{{| D_{x} |}^{2} + {| D_{y} |}^{2}} .

As ${\tilde{S}}_{x} = i \bar{k_{x}} \tilde{φ}$ and ${\tilde{S}}_{y} = i \bar{k_{y}} \tilde{φ}$ , the frequency response H defined as the ratio of the reconstructed wavefront amplitude to the true wavefront amplitude associated with the measured slopes at a given spatial frequency point is

H = \frac{\tilde{\hat{φ}}}{\tilde{φ}} = \frac{(D_{x}^{*} A_{x} D_{x, 0} + D_{y}^{*} A_{y} D_{y, 0})}{{| D_{x} |}^{2} + {| D_{y} |}^{2}},

where new notations D_x _,0, D_y _,0 were introduced in place of

i \bar{k_{x}}

and

i \bar{k_{y}}

, respectively.

Applying the Simpson derivative and average operators [Eqs. (35), (36)], we obtain

H_{Simpson} = \frac{1}{2} \frac{ω_{x} sin ω_{x} (2 + cos ω_{x}) + ω_{y} sin ω_{y} (2 + cos ω_{y})}{{sin}^{2} ω_{x} + {sin}^{2} ω_{y}},

where

\bar{k_{x}} Δ x \equiv ω_{x} = 2 π f_{x}

and

\bar{k_{y}} Δ y \equiv ω_{y} = 2 π f_{y}

. f_x and f_y are normalized frequencies ranging from −0.5 to 0.5. It is assumed Δx = Δy. The frequency response of H _Simpson has eight singularities on the four corners and side centers. Except for the region near the poles, the frequency response is nearly unity everywhere, which proves better accuracy of the Simpson rule than the conventional reconstructors (Refer Fig. 1 of [11]) over all spectrum.

The singularities can be removed by introducing the Phillips regularization term [16] to the right hand side of Eq. (28):

\begin{array}{l} ɛ_{reg} = λ \sum_{i, j} {\frac{1}{{(2 Δ x)}^{2}} {[\hat{φ} (i, j + 1) - 2 \hat{φ} (i, j) + \hat{φ} (i, j - 1)]}^{2} \\ + \frac{1}{{(2 Δ y)}^{2}} {[\hat{φ} (i + 1, j) - 2 \hat{φ} (i, j) + \hat{φ} (i - 1, j)]}^{2}} . \end{array}

In frequency domain, this transforms into

ɛ_{reg} = λ \frac{1}{N^{2}} \sum_{i, j} {{| D_{x, reg} \tilde{\hat{φ}} |}^{2} + {| D_{y, reg} \tilde{\hat{φ}} |}^{2}},

where

D_{x, reg} = \frac{2}{Δ x} {sin}^{2} (\frac{\bar{k_{x}} Δ x}{2})

and similarly for D _y,reg.

The denominator of the Simpson frequency response will have an additional term of λ (|D_x, _reg|² + |D_y, _reg|²) that removes the singularity. The regularized frequency response is

H_{Simpson} (λ) = \frac{1}{3} \frac{ω_{x} sin ω_{x} (2 + cos ω_{x}) + ω_{y} sin ω_{y} (2 + cos ω_{y})}{{sin}^{2} ω_{x} + {sin}^{2} ω_{y} + 4 λ ({sin}^{4} \frac{1}{2} ω_{x} + {sin}^{4} \frac{1}{2} ω_{y})} .

The 1-D frequency response (ω_y = 0) with the regularization term has a second maximum near high spatial frequency for sufficiently small λ (<0.08) [Fig. 2(b)]. The free parameter λ can be fixed to a value such that the second maximum is equal to one. Numerically determined value of λ for such condition is 0.07489. This choice of λ gives only 3% error in wavefront amplitude over 80% of the frequency range. Another choice can be λ = 0.07026 which balances the local maximum and minimum around 1. The second option reduces the maximum deviation below 2.2% within 85% of the spectral range. Figure 2(a) shows three-dimensional view of the frequency response of the Simpson-rule reconstructor with λ = 0.07489. The 1-D response is shown in Fig. 2(b). The solid line is calculated from analytic expression [Eq. (43)] whereas the circles are from numerical simulations. The numerical simulation consists of steps of generating slopes from sinusoid wavefronts at a given spatial frequency using analytic derivatives and of reconstructing the wavefront and comparing the ratio between the original and the reconstructed wavefront amplitude at that frequency. The reconstruction algorithm used in the simulation will be explained in detail in the following section. The result shows good agreement with the analytic curve.

Fig. 2 Frequency response of Simpson reconstructor with λ =0.07489 (a) 3-D view of the frequency response. (b) cross-section along f_x axis. Solid line is calculated from the analytic expression and circles are from simiulations

Download Full Size | PDF

3.3. Iterative algorithm with regularization terms

The frequency-domain analysis does not provide a detailed picture how the successive over-relaxation method can be applied in spatial domain iteration especially around the boundary of slopes map. Resolving the stationary condition with the regularization term introduces additional terms on the left hand side of Eq. (29). These are fully written out using g flags.

\begin{array}{l} {(2 Δ x)}^{2} \frac{\partial ɛ_{reg}}{\partial \hat{φ} (i, j)} = λ g_{L} [\hat{φ} (i, j) - 2 \hat{φ} (i, j - 1) + \hat{φ} (i, j - 2)] \\ + λ g_{R} [\hat{φ} (i, j + 2) - 2 \hat{φ} (i, j + 1) + \hat{φ} (i, j)] \\ + λ {g'}_{U} [\hat{φ} (i, j) - 2 \hat{φ} (i - 1, j) + \hat{φ} (i - 2, j)] + λ {g'}_{D} [\hat{φ} (i + 2, j) - 2 \hat{φ} (i + 1, j) + \hat{φ} (i, j)] \\ - 2 λ g_{L R} [\hat{φ} (i, j + 1) - 2 \hat{φ} (i, j) + \hat{φ} (i, j, - 1)] - 2 λ {g'}_{U D} [\hat{φ} (i + 1, j) - 2 \hat{φ} (i, j) + \hat{φ} (i - 1, j)] . \end{array}

g_LR or g_UD is one only if there exist both points to the right and left or up and down of the point (i, j), respectively and zero otherwise. The iteration formula (32) will be modified to

\bar{\hat{φ}} (i, j) = [\bar{\hat{φ_{0}}} (i, j) + Δ S (i, j) + λ Δ \hat{φ} (i, j)] / [(1 + λ) (g_{L} + g_{R} + {g'}_{U} + {g'}_{D}) + 4 λ (g_{L R} + {g'}_{U D})],

where

\begin{array}{l} Δ \hat{φ} \equiv g_{L} [2 \hat{φ} (i, j - 1) - \hat{φ} (i, j - 2)] + g_{R} [2 \hat{φ} (i, j + 1) - \hat{φ} (i, j + 2)] \\ + {g'}_{U} [2 \hat{φ} (i - 1, j) - \hat{φ} (i - 2, j)] + {g'}_{D} [2 \hat{φ} (i + 1, j) - \hat{φ} (i + 2, j)] \\ + 2 g_{L R} [\hat{φ} (i, j + 1) + \hat{φ} (i, j - 1)] + 2 {g'}_{U D} [\hat{φ} (i + 1, j) + \hat{φ} (i - 1, j)] . \end{array}

4. Hexagonal band-limited reconstructor

Iterative Fourier-domain band-limited reconstruction method as proposed in [11] provides unity frequency response over all spatial bandwidth as long as the signals are band-limited. The band-limited reconstructor was designed for the Southwell geometry; the reconstruction point and the corresponding measured slopes are all at the same point on a rectangular grid. It was shown that band-limited derivative operators are also available for Hudgin and Fried geometries [11]. Table 1 summarizes the three operators belonging to each geometry.

Table 1. Summary of Band-Limited Derivative Operators (D _x,0)

View Table | View all tables in this article

Here we present band-limited reconstructors for hexagonal arrays. Hexagonal geometry may be well suited for adaptive optic systems for large telescopes with hexagonal mirror arrays (e.g., James Webb) [17]. Large deformable mirrors used in laser fusion facilities (National Ignition Facility) [18] also have hexagonal actuator patterns. The number density of lenslets is slightly higher in hexagonal geometry than square. Figure 3 shows two possible hexagonal arrays. In Fig. 3(a) the unit hexagon is lying on its facet (prostrate hexagon), whereas in Fig. 3(b) the unit hexagon is standing on the apex (standing hexagon). The circles indicate the measured slopes points and the ×’s are reconstruction points. The reconstruction points are set on a rectangular grid because rectangular arrays are easier to visualize in most computer programs.

Fig. 3 (a) Prostrate hexagon array; (b) standing hexagon array.

Download Full Size | PDF

As will become clear in the discussion of the algorithm, we need procedures for performing DFTs on a hexagonal array. This can be achieved by grouping the slopes into two interleaving groups, which are indicated by red (group 1) and black (group 2) circles in Fig. 3. Starting from Fig. 3(a) geometry, it can be shown that the DFT of the slopes on the rectangular array (on × points) can be expressed as the linear sum of smaller DFTs of the sub-group 1 and 2 (see Appendix).

{\tilde{S}}_{x} = [\begin{array}{l} {\tilde{S}}_{1, x} + T_{hex} \circ {\tilde{S}}_{2, x} \\ {\tilde{S}}_{1, x} - T_{hex} \circ {\tilde{S}}_{2, x} \end{array}] \equiv [\begin{matrix} {\tilde{S}}_{x, I} \\ {\tilde{S}}_{x, II} \end{matrix}]

T _hex is a matrix whose size is M by N (i.e. the size of either array 1 or array 2) and “○” denotes entrywise matrix multiplication. The pth row and qth column element of T _hex is

T_{hex} (p, q) = exp [- \frac{π i}{M} p - π i 𝖲 (q)] .

The combined total DFT array is therefore a vertical concatenation of the two matrices. On the other hand, the resulting total DFT matrix for Fig. 3(b) geometry is a horizontal concatenation,

\begin{array}{l} {\tilde{S}}_{x} & = & [{\tilde{S}}_{1, x} + T_{hex} \circ {\tilde{S}}_{2, x} {\tilde{S}}_{1, x} - T_{hex} \circ {\tilde{S}}_{2, x}] \\ \equiv & [{\tilde{S}}_{x, I} {\tilde{S}}_{x, II}] \end{array}

T_{hex} (p, q) = exp [- π i 𝖲 (p) - \frac{π i}{N} q] .

The same combination rule applies to y-slope measurements.

The above decomposition technique can be inverted such that each subgroup of the hexagonal array can also be expressed as the linear sum of the block I and II of the rectangular array. This inversion is only used for wavefront points in the algorithm, which is

\tilde{φ_{1}} = \frac{1}{2} (\tilde{φ_{I}} + \tilde{φ_{II}}),

\tilde{φ_{2}} = \frac{1}{2} T_{hex}^{*} \circ (\tilde{φ_{I}} - \tilde{φ_{II}}) .

Using the basic results obtained in Sec. 2 and the DFT procedures for the hexagonal arrays in this section, the band-limited reconstruction algorithm for hexagonal slopes arrays can be implemented as shown in the flowchart in Fig. 4.

Fig. 4 Flowchart of band-limited reconstruction for a hexagonal geometry. ℱ is the band-limited filter function [Eq. (53)]. Ω₁ and Ω₂ are the regions where the slopes data exist. The subscripts ‘m’ in Step 2 and 7 denote the measured slopes. IDFT stands for inverse discrete Fourier transform.

Download Full Size | PDF

Step 1 consists of fitting the slopes over low order polynomials, e.g. third order, which will significantly reduce non-band-limited components of the wavefront. If the regions of interest are disconnected, the fitting has to be performed per each region. Owing to the sum requirement [Eq. (6)], a column and a row are appended to the edge of the measured slope matrices (group 1 and 2 separately), which will satisfy the zero sum conditions in x and y directions.

Step 2 initializes the slopes with measured values. Steps 3-8 form a closed loop required for extrapolating slopes outside the non-rectangular region. The iteration is not required if the region is rectangular.

Slopes in group 1 and 2 are separately Fourier-transformed using Eqs. (47)–(50) in Step 3. In Step 4, wavefront matrices corresponding to each block (I or II) are reconstructed in Fourier domain using the band-limited filter function which is

ℱ [{\tilde{S}}_{x} (p, q), {\tilde{S}}_{y} (p, q)] = \frac{D_{x}^{*} (p, q) {\tilde{S}}_{x} (p, q) + D_{y}^{*} (p, q) {\tilde{S}}_{y} (p, q)}{{| D_{x} (p, q) |}^{2} + {| D_{y} (p, q) |}^{2}},

where the band-limited derivative operators D_x and D_y are defined as

D_{x} (p, q) = \frac{2 π i}{Δ x} 𝖲 (q),

D_{y} (p, q) = \frac{2 π i}{(2 Δ y)} 𝖲 (p),

for a prostrate hexagon array and

D_{x} (p, q) = \frac{2 π i}{(2 Δ x)} 𝖲 (q),

D_{y} (p, q) = \frac{2 π i}{Δ y} 𝖲 (p),

for a standing hexagon array.

Step 5 creates wavefront group 1 and 2 by using Eqs. (51)–(52). In Step 3 and 5, the correct T _hex must be used according to its geometry. In Step 6, derivative operators are applied to these temporary wavefront matrices to obtain slopes in group 1 and 2, respectively. These new slopes are different from the measured slopes. We leave the values external to the boundary untouched while restoring the internal values to the original measured slopes. The difference between the measured slope and the calculated slopes decreases over the course of iterations. Step 8 determines whether this difference is within tolerance. Once the convergence criterion is met, the wavefront matrices generated in Step 4 ( $\tilde{φ_{I}}$ , $\tilde{φ_{II}}$ ) are combined to form a single matrix either by vertical or horizonal concatenation depending on the hexagon geometry and inverse Fourier-transformed to the spatial domain to produce the final result in Step 9. Small terms in the imagingary part of the solution can be neglected.

The band-limited algorithms as shown in Fig. 2 of [11] and Fig. 4 can be used together with a non-band-limited filter function, which enables a convenient way to switch between different algorithms. The reconstruction algorithms proposed here are not limited to a specific boundary shape.

5. Error propagation

The wavefront reconstructors have traditionally been characterized with so called error propagation curve. This indicates the sensitivity of the noise in the reconstructed phase to the noise in the slopes measurements. Early numerical and theoretical works show that this sensitivity is a logarithmic function of the number of measurement points [1, 2, 15]. The noise propagation coefficient will be calculated using discrete samples and frequency domain filter functions. We limit the scope to the rectangular area.

Let σ_φ be the root-mean-square of the reconstructed phase φ. According to the Wiener-Khintchin theorem,

σ_{φ}^{2} = \frac{1}{N_{t}} Σ 〈 Δ φ^{2} 〉 = \frac{1}{N_{t} N^{2}} Σ 〈 {| \tilde{Δ φ} |}^{2} 〉,

where 〈·〉 denotes ensemble average of the quantity inside.

According to linear stochastic system theory, the power spectrum of input and output signal is related by the absolute square of the linear system function. In the case of wavefront reconstruction dictated by the following linear response

\tilde{φ} = \frac{(D_{x}^{*} A_{x} {\tilde{S}}_{x} + D_{y}^{*} A_{y} {\tilde{S}}_{y})}{{| D_{x} |}^{2} + {| D_{y} |}^{2} + λ ({| D_{x, reg} |}^{2} + {| D_{y, reg} |}^{2})},

the corresponding stochastic response in power spectrum is

〈 {| \tilde{Δ φ} |}^{2} 〉 = \frac{({| D_{x} A_{x} \tilde{Δ S_{x} |}}^{2} + {| D_{y} A_{y} \tilde{Δ S_{y}} |}^{2})}{{[{| D_{x} |}^{2} + {| D_{y} |}^{2} + λ ({| D_{x, reg} |}^{2} + {| D_{y, reg} |}^{2})]}^{2}} .

Assuming that $\tilde{Δ S_{x}}$ and $\tilde{Δ S_{y}}$ are uncorrelated white noise with variance of $\tilde{σ_{S}^{2}}$ for each, and since $σ_{S}^{2} = \frac{1}{N_{t}} {\tilde{σ_{S}}}^{2}$ , the noise propagation coefficient η is

η = \frac{σ_{φ}^{2}}{h^{2} σ_{S}^{2}} = \frac{1}{L^{2}} \sum_{k_{x}, k_{y}} \frac{({| D_{x} A_{x} |}^{2} + {| D_{y} A_{y} |}^{2})}{{[{| D_{x} |}^{2} + | {D_{y} |}^{2} + λ ({| D_{x, reg} |}^{2} + {| D_{y, reg} |}^{2})]}^{2}},

where

σ_{φ}^{2} = \frac{1}{N_{t}} 〈 {| \tilde{Δ φ} |}^{2} 〉

, h = Δx = Δy and L is the aperture size. This result is equivalent to Noll’s [19] in the case of band-limited operators.

Table 2 summarizes finite-difference derivative/averaging operators for four geometries to be used with Eq. (61).

Table 2. Summary of the Frequency-Domain Operators for the Associated Finite Difference Schemes

View Table | View all tables in this article

The right hand side of Eq. (61) is inversely proportional to |D ₀ _,x|² + |D ₀ _,y|² for band-limited reconstruction and difficult to visualize in linear scale. We define a ‘noise response function’ S_N with the inverse power dependence removed as follows:

S_{N} = \frac{({| D_{0, x} |}^{2} + {| D_{0, y} |}^{2}) ({| D_{x} A_{x} |}^{2} + {| D_{y} A_{y} |}^{2})}{{[{| D_{x} |}^{2} + {| D_{y} |}^{2} + λ ({| D_{x, reg} |}^{2} + {| D_{y, reg} |}^{2})]}^{2}} .

It can be shown by the Cauchy-Schwartz inequality that the noise response is always larger than or equal to absolute squared of the frequency response,

{| H |}^{2} \leq S_{N} .

The inequality (63) shows that the error propagation is intimately related with the frequency response of a reconstructor. The lower bound of η is

\sum_{k_{x}, k_{y}} \frac{{| H |}^{2}}{(k_{x}^{2} + k_{y}^{2})} \leq \sum_{k_{x}, k_{y}} \frac{S_{N}}{(k_{x}^{2} + k_{y}^{2})} = η

From this one can expect that the Southwell reconstructor will have the lowest lower bound and the Fried reconstructor the highest. It agrees with the result of Zou [20].

The analytic expression for η can be calculated and fit to a logarithmic curve, although the logarithm dependence is only approximate except for the band-limited reconstructors. The result is summarized in the second column of Table 3. Singularity points were excluded in the summation over spatial-frequency space. The third column shows the simulated η obtained by running actual reconstructors with zero slopes input with Gaussian noise. The Fried reconstructor was not tested. Two hundred realizations were performed at each N, where N ² is the number of points. N was varied from 10 to 100 by 10. The logarithm fit over the averaged η is shown in the column. The multiplicative coefficients roughly agree with the analytic ones up to the second decimal point, but the additive constants from simulation are always estimated higher than the calculated ones. The offset is about 0.2771 on average. The discrepancy appears to come from the apparent inconsistency in assuming white noise in slopes power spectrum and the use of band-limited derivative formalism. For example, the reconstructed wavefront from white noise spectrum always has some amount of low-order polynomial terms, which cannot be represented by Eq. (59). The constant offset 0.2771 therefore can be considered as the ratio of energy conversion from white noise to non-band-limited signals.

Table 3. Summary of Noise Propagation

View Table | View all tables in this article

The legacy formulas of noise propagation for each reconstructor are also shown in the fourth column of Table 3, quoted from the three authors’ original publications [1, 2, 15]. The quoted Southwell η is estimated only from the graph in the original paper since no explicit formula was given. Noll’s calculation essentially corresponds to band-limited case. Considering the fact that there is some ambiguity in the determination of the constant offset, at least the multiplicative coefficient of the Fried formula comes close to our analytic result; whereas there is about a factor-of-2 difference in the Southwell and Noll’s expression compared with ours. On the other hand, Hudgin’s formula does not agree with our results. Fried’s formulas is based on comparatively large number of N (≤ 39) compared with Southwell and Hudgin’s calculations (N ≤ 20).

6. Conclusions

We have presented detailed derivations of band-limited derivative operators in frequency domain. These are important tools for characterizing and improving frequency response of wavefront reconstructors over broad bandwidth. Two new wavefront reconstructors were proposed utilizing these tools. The reconstructors were designed to be accurate up to high spatial frequency. The first one is based on the Simpson integration rule. The bandwidth of the frequency response of this reconstructor, after being regularized, is excellent up to 85% of the spatial frequency range. A successive-over-relaxation iterative solver was presented in detail where the outermost samples are elegantly handled using g flags. The frequency response behavior of the iterative solver agrees well with the predicted frequency response curve. The second reconstructor is an extension of the band-limited reconstruction algorithm previously developed; the measurement points are on a hexagonal array instead of a rectangular array. A Fourier-domain iterative algorithm was proposed for two types of hexagonal arrays. As was pointed before in [11], the reconstruction process must be preconditioned with low-order polynomial fit. The Simpson-rule based algorithm works purely in spatial domain; therefore it is computationally less complex than band-limited algorithms, whereas the latter provides flexibility against any geometry change. Fourier-domain algorithms have a potential of boosting reduction speed with the help of digital signal processors.

The new wavefront reconstructors are compared with the traditional reconstructors in terms of noise propagation properties through a generalized noise propagation expression. The analytically calculated noise propagation coefficients are consistent with the numerical fit deduced from our own simulations. However we did not find universal agreement with the published results.

The broad-bandwidth wavefront reconstructors developed here are used in wavefront reduction software for characterizing focal spot of OMEGA EP laser beams [13]. The importance of band-limited reconstructor was well illustrated in [21] for a closed-loop wavefront shaping application. One may also find applications in the study of metrology and atmospheric turbulence [22].

Appendix

Equation (47) will be derived for the prostrate hexagon geometry. Supposing the x-slopes are measured on the M by N rectangular grid instead of the hexagonal grid and denoting them with primes, the DFT of S′_x

{\tilde{S}'}_{x} (p, q) = \sum_{n = 0}^{N - 1} \sum_{m = 0}^{M - 1} {S'}_{x} (m, n) exp (- \frac{2 π i}{M} m p - \frac{2 π i}{N} n q)

is made of two smaller DFTs as follows. Grouping the odd rows and even rows separately,

\begin{array}{l} {\tilde{S}'}_{x} (p, q) = \sum_{n = 0}^{N - 1} \sum_{m' = 0}^{M' - 1} {S'}_{x} (2 m', n) exp (- \frac{2 π i}{M'} m' p - \frac{2 π i}{N} n q) \\ + \sum_{n = 0}^{N} \sum_{m ″ = 0}^{M' - 1} {S'}_{x} (2 m ″ + 1, n) exp (- \frac{2 π i}{M'} m ″ p - \frac{2 π i}{N} n q - \frac{π i}{M'} p), \end{array}

where the odd rows are grouped into the first term and the even rows into the second term. M is assumed even and M = 2M′. Note that M in the main text around Eqs. (47), (48) is same as M′ here. The Fourier domain indices are p and q which range from 0 to M – 1 and from 0 to N – 1, respectively.

Denoting the first half range of index p as p′ (p′ = 0, ...,M′ − 1), the above equation becomes

{\tilde{S}'}_{x} (p = p', q) = \tilde{{S'}_{x, 1}} (p', q) + exp (- \frac{π i}{M'} p') \tilde{{S'}_{x, 2}} (p', q)

for the first half rows and

{\tilde{S}'}_{x} (p = p' + M', q) = \tilde{{S'}_{x, 1}} (p', q) - exp (- \frac{π i}{M'} p') \tilde{{S'}_{x, 2}} (p', q)

for the other half rows. We used subscript 1 for odd rows and 2 for even rows.

Slopes on the hexagonal grid are $\tilde{S_{x, 1}} (p, q) = \tilde{{S'}_{x, 1}} (p, q)$ and $\tilde{S_{x, 2}} (p, q) = exp [π i 𝖲 (q)] \tilde{{S'}_{x, 2}} (p, q)$ using half-sample shift operator 𝖲. The shift operator and the extra phase factor constitute T _hex as defined in the main text [Eq. (48)]. Therefore Eq. (67), (68) can be expressed as a block matrix form as in Eq. (47). The first and second block will be denoted as $\tilde{S_{x, I}}$ and $\tilde{S_{x, II}}$ , respectively. The DFT of the group 1 or 2 slopes is to be interpreted as operating on matrices of size M′ by N.

The relation between the slopes group 1 and 2 on the hexagonal grid and the first and second blocks on the rectangular grid can be inverted. In other words, the slope group 1 or 2 can be expressed as the linear sum of the first and the second blocks in the Fourier domain.

\tilde{S_{x, 1}} = \frac{1}{2} (\tilde{S_{x, I}} + \tilde{S_{x, II}}),

\tilde{S_{x, 2}} = \frac{1}{2} T_{hex}^{*} \circ (\tilde{S_{x, I}} - \tilde{S_{x, II}}) .

The block-decomposition technique and the inversion formula can be applied to y-slopes and wavefront points. The results for the standing hexagon geometry can be similarly derived.

Acknowledgments

This work was supported by the U.S. Department of Energy Office of Inertial Confinement Fusion under Cooperative Agreement No. DE-FC52-08NA28302, the University of Rochester, and the New York State Energy Research and Development Authority. The support of DOE does not constitute an endorsement by DOE of the views expressed in this article.

References and links

1. D. L. Fried, “Least-squares fitting a wavefront distortion estimate to an array of phasedifference measurements,” J. Opt. Soc. Am. 67, 370–375 (1977). [CrossRef]

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

3. K. Freischlad and C. L. Koliopoulos, “Modal estimation of a wavefront from difference measurements using the discrete Fourier transform,” J. Opt. Soc. Am. A 3, 1852–1861 (1986). [CrossRef]

4. F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” Appl. Opt. 30, 1325–1327 (1991). [CrossRef] [PubMed]

5. 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. A 19, 2100–2111 (2002). [CrossRef]

6. L. A. Poyneer and B. Macintosh, “Spatially filtered wave-front sensor for high-order adaptive optics,” J. Opt. Soc. Am. A 21, 810–819 (2004). [CrossRef]

7. A. Dubra, C. Paterson, and C. Dainty, “Wave-front reconstruction from shear phase maps by use of the discrete Fourier transform,” Appl. Opt. 43, 1108–1113 (2004). [CrossRef] [PubMed]

8. C. Elster and I. Weingärtner, “High-accuracy reconstruction of a function f(x) when only $\frac{d}{d x} f (x)$ or $\frac{d^{2}}{d x^{2}} f (x)$ is known at discrete measurement points,” in X-Ray Mirrors, Crystals, and Multilayers II, A. K. Freund, A. T. Macrander, T. Ishikawa, and J. L. Wood, eds., Proc. SPIE 4782, 12–20 (2002).

9. J. Campos, L. P. Yaroslavsky, A. Moreno, and M. J. Yzuel, “Integration in the Fourier domain for restoration of a function from its slope: comparison of four methods,” Opt. Lett. 27, 1986–1988 (2002). [CrossRef]

10. L. P. Yaroslavsky, A. Moreno, and J. Campos, “Frequency responses and resolving power of numerical integration of sampled data,” Opt. Express 13, 2892–2905 (2005). [CrossRef] [PubMed]

11. S.-W. Bahk, “Band-limited wavefront reconstruction with unity frequency response from Shack-Hartmann slopes measurements,” Opt. Lett. 33, 1321–1323 (2008). [CrossRef] [PubMed]

12. M. Servin, D. Malacara, and J. L. Marroquin, “Wave-front recovery from two orthogonal sheared interferograms,” Appl. Opt. 35, 4343–4348 (1996). [CrossRef] [PubMed]

13. J. Bromage, S.-W. Bahk, D. Irwin, J. Kwiatkowski, A. Pruyne, M. Millecchia, M. Moore, and J. D. Zuegel, “A focal-spot diagnostic for on-shot characterization of high-energy petawatt lasers,” Opt. Express 16, 16561–16572 (2008). [PubMed]

14. C. E. Shannon, “Communication in the presence of noise,” Proc. Inst. Radio Eng. 37, 10–21 (1949).

15. W. H. Southwell, “Wavefront estimation from wavefront slope measurements,” J. Opt. Soc. Am. 70, 998–1009 (1980). [CrossRef]

16. B. L. Phillips, “A technique for the numerical solution of certain integral equations of the first kind,” J. ACM 9, 84–97 (1962). [CrossRef]

17. See http://www.jwst.nasa.gov/.

18. S. E. Winters, J. H. Chung, and S. A. Velinsky, “Modeling and control of a deformable mirror,” J. Dyn. Sys. Meas. Control 124, 297–302 (2002). [CrossRef]

19. R. J. Noll, “Phase estimates from slope-type wave-front sensors,” J. Opt. Soc. Am. 68, 139–140 (1978). [CrossRef]

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

21. S.-W. Bahk, E. Fess, B. E. Kruschwitz, and J. D. Zuegel, “A high-resolution, adaptive beam-shaping system for high-power lasers,” Opt. Express 18, 9151–9163 (2010). [CrossRef] [PubMed]

22. T. W. Nicholls, G. D. Boreman, and J. C. Dainty, “Use of a Shack-Hartmann wave-front sensor to measure deviations from a Kolmogorov phase spectrum,” Opt. Lett. 20, 2460–1462 (1995). [CrossRef] [PubMed]

Finite-difference scheme	D_x	A_x
Southwell [15]	$\frac{1}{Δ x} [exp (i \bar{k_{x}} Δ x) - 1]$	$\frac{1}{2} [exp (i \bar{k_{x}} Δ x) + 1]$
Hudgin [2]	$\frac{1}{Δ x} [exp (i \bar{k_{x}} Δ x) - 1]$	1
Fried [1]	$\frac{1}{2 Δ x} [exp (i \bar{k_{x}} Δ x) - 1] [exp (i \bar{k_{y}} Δ y) + 1]$	1
Simpson	$\frac{1}{2 Δ x} [exp (i \bar{k_{x}} Δ x) - exp (- i \bar{k_{x}} Δ x)]$	$\frac{1}{6} exp (i \bar{k_{x}} Δ x) + 4 + exp (- i \bar{k_{x}} Δ x)]$

Numerical scheme	Calculated η	Simulated η	Quoted η
Southwell difference	−0.1211 + 0.1591lnN	0.1356 + 0.1638lnN	−0.1237 + 0.3164lnN
Hudgin difference	0.0485 + 0.1592lnN	0.3456 + 0.1528lnN	0.561 + 0.103lnN
Fried difference	−0.0865 + 0.3202lnN	—	0.6558 + 0.3206lnN
Band-limited	−0.0285 + 0.1594lnN	0.2449 + 0.1615lnN	0.1072 + 0.318lnN
Simpson difference	0.0314 + 0.1610lnN	0.3125 + 0.1647lnN	—

Finite-difference scheme	D_x	A_x
Southwell [15]	$\frac{1}{Δ x} [exp (i \bar{k_{x}} Δ x) - 1]$	$\frac{1}{2} [exp (i \bar{k_{x}} Δ x) + 1]$
Hudgin [2]	$\frac{1}{Δ x} [exp (i \bar{k_{x}} Δ x) - 1]$	1
Fried [1]	$\frac{1}{2 Δ x} [exp (i \bar{k_{x}} Δ x) - 1] [exp (i \bar{k_{y}} Δ y) + 1]$	1
Simpson	$\frac{1}{2 Δ x} [exp (i \bar{k_{x}} Δ x) - exp (- i \bar{k_{x}} Δ x)]$	$\frac{1}{6} exp (i \bar{k_{x}} Δ x) + 4 + exp (- i \bar{k_{x}} Δ x)]$

Numerical scheme	Calculated η	Simulated η	Quoted η
Southwell difference	−0.1211 + 0.1591lnN	0.1356 + 0.1638lnN	−0.1237 + 0.3164lnN
Hudgin difference	0.0485 + 0.1592lnN	0.3456 + 0.1528lnN	0.561 + 0.103lnN
Fried difference	−0.0865 + 0.3202lnN	—	0.6558 + 0.3206lnN
Band-limited	−0.0285 + 0.1594lnN	0.2449 + 0.1615lnN	0.1072 + 0.318lnN
Simpson difference	0.0314 + 0.1610lnN	0.3125 + 0.1647lnN	—

Highly accurate wavefront reconstruction algorithms over broad spatial-frequency bandwidth

Abstract

1. Introduction

2. Band-limited derivative

3. Simpson reconstructor

3.1. Simpson iterator

3.2. Frequency response and regularization

3.3. Iterative algorithm with regularization terms

4. Hexagonal band-limited reconstructor

5. Error propagation

6. Conclusions

Appendix

Acknowledgments

References and links

Cited By

Figures (4)

Tables (3)

Equations (70)

Optics Express

Geometry	D _x,0
Southwell	$i \bar{k_{x}}$
Hudgin	$i \bar{k_{x}} exp (i \bar{k_{x}} \frac{Δ x}{2})$
Fried	$i \bar{k_{x}} exp (i \bar{k_{x}} \frac{Δ x}{2} + i \bar{k_{y}} \frac{Δ y}{2})$