A note on speed comparison of wavefront reconstruction approaches

The two main reconstruction geometries are namely Hudgin and Fried. In the Hudgin geometry, gradients are the first differences between neighbouring actuators\cite{poyneer2003advanced}. In the Fried geometry, gradients are the average of the two first-differences on a square. Many adaptive optics systems use vector-matrix-multiply (VMM) reconstructors to convert gradient measurements to wavefront phase estimates. The computation time of the reconstruction using the VMM method scales as $O(n^2)$\cite{PoyneerFastFFTreconWFR}.

The reconstruction of the wavefront by means of the FFT was proposed by Freischlad and Koliopoulos\cite{freischlad1985wavefront} for square apertures on the Hudgin geometry. In a further paper\cite{FreischladFFTreconWFR} the authors derived methods for additional geometries, including the Fried geometry, which uses one Shack-Hartmann (SH) sensor. Freischlad also considered the case of small circular apertures\cite{freischladwfrfft} and the boundary problem was identified.


Computational speed comparison
The computational speed of Hudgin-FT and Fried-FT are limited only by the FFT\cite{poyneer2003advanced}. The extra processing to solve the boundary problem is of a lower order of growth computationally. Therefore FFT implementations have computational costs that scale as $O(n \log n)$. However, the implementation of Fried-FT requires potentially 2 times as much total computation as the Hudgin-FT. For the $64 \times 64$ grid, the FFT can be calculated on currently available systems in around 1 ms\cite{poyneer2003advanced}. If $N$ is a power of two the spatial filter, operations can be implemented with FFT's very efficiently. The computational requirements then scale as $O( N^2 log_2 N )$ rather than as $O( N^4 )$ in the direct vector-matrix multiplication approach. The modified Hudgin takes half as much computation as the Fried geometry model\cite{poyneer2003advanced}.


Comparison of FFT WFR and Zernike reconstruction speed

A very interesting paper appeared in the Journal of Refract Surgery\cite{dai2006comparison}. A comparison between Fourier and Zernike reconstructions was performed. In the paper\cite{dai2006comparison}, noise-free random wavefronts were simulated with up to the 15th order of Zernike polynomials. Fourier full reconstruction was more accurate than Zernike reconstruction from the 6th to the 10th orders for low-to-moderate noise levels. Fourier reconstruction was found to be approximately \textbf{100 times faster than Zernike reconstruction}. For Zernike reconstruction, however, the optimal number of orders must be chosen manually. The optimal Zernike order for Zernike reconstruction is lower for smaller pupils than larger pupils. The paper\cite{dai2006comparison} concludes that the FFT WFR is faster and more accurate than Zernike reconstruction, makes optimal use of slope information, and better represents ocular aberrations of highly aberrated eyes.


Noise propagation

Analysis and simulation show that for apertures just smaller than the square reconstruction grid (DFT case), the noise propagations of the FT methods are favourable. For the Hudgin geometry, the noise propagator grows with $O(\ln n)$. For the Fried geometry, the noise propagator is best-fit by a curve that is quadratic in the number of actuators, or $O(\ln^2 n)$. For fixed power-of-two sized grids (required to obtain the speed of the FFT for all aperture sizes) the noise propagator becomes worse when the aperture was much smaller than the grid\cite{poyneer2003advanced}.


Shack-Hartmann sensor gain

The Shack-Hartmann WFS produces a measurement which deviates from the exact wavefront slope. The exact shape of this curve depends on number of pixels used per sub-aperture and the centroid computation method (see \cite{hardyAObook} section 5.3.1 for a representative set of response curves). The most important feature of the response curve is that even within the linear response range, the gain of the sensor is not unity\cite{poyneer2003advanced}. This gain is important in the open-loop: in a closed loop this problem in mitigated by the overall control loop gain, which can be adjusted instead.

References:

\begin{thebibliography}{1}  \bibitem{poyneer2003advanced} L.A. Poyneer. \newblock {Advanced techniques for Fourier transform wavefront reconstruction}. \newblock In {\em Proceedings of SPIE}, volume 4839, page 1023, 2003.  \bibitem{PoyneerFastFFTreconWFR} Lisa~A. Poyneer, Donald~T. Gavel, and James~M. Brase. \newblock Fast wave-front reconstruction in large adaptive optics systems with   use of the fourier transform. \newblock {\em J. Opt. Soc. Am. A}, 19(10):2100--2111, Oct 2002.  \bibitem{freischlad1985wavefront} K.~Freischlad and C.L. Koliopoulos. \newblock {Wavefront reconstruction from noisy slope or difference data using   the discrete Fourier transform}. \newblock 551:74--80, 1985.  \bibitem{FreischladFFTreconWFR} Klaus~R. Freischlad and Chris~L. Koliopoulos. \newblock Modal estimation of a wave front from difference measurements using   the discrete fourier transform. \newblock {\em J. Opt. Soc. Am. A}, 3(11):1852--1861, Nov 1986.  \bibitem{freischladwfrfft} Klaus~R. Freischlad. \newblock Wave-front integration from difference data. \newblock {\em Interferometry: Techniques and Analysis}, 1755(1):212--218,   1993.  \bibitem{dai2006comparison} G.~Dai. \newblock {Comparison of wavefront reconstructions with Zernike polynomials and   Fourier transforms.} \newblock {\em Journal of refractive surgery}, 22(9):943--948, 2006.  \bibitem{hardyAObook} John~W. Hardy. \newblock {\em {Adaptive optics for astronomical telescopes}}. \newblock Oxford University Press, USA, 1998.  \end{thebibliography}

Reconstruction of the wavefront: Boundary problems on a circular aperture

The ``boundary problem'' leads to large uncorrectable errors: if the slopes in the aperture are simply zero-padded, large errors occur across the aperture. The example of the errors can be seen below, where the reconstruction of the wavefront using FFT WFR modified Hudgin method is shown with the circular aperture and central obscuration.

Reconstruction of the wavefront using FFT WFR modified Hudgin method is shown with the circular aperture and central obscuration (outer is 0.9, and inner is 0.1).

In the case of astronomical telescopes, the gradients are typically available on
a circular aperture. The measurement data cannot be zero padded as it leads to large errors on boundaries. Such errors do not decrease with system size\cite{PoyneerFastFFTreconWFR}.

Because of the spatial periodicity and closed-path-loop conditions (see below), zero-padding the gradient measurements is incorrect. Doing so violates both conditions in general\cite{PoyneerFastFFTreconWFR}. These inconsistencies manifest themselves in errors that span the aperture. The errors do not become less significant as the aperture size increases: unlike the squareaperture case, the amplitude of the error remains large and spans the circular aperture.



Assumptions in the wavefront reconstruction that are not always true
There are two key assumptions in the wavefront reconstruction that must be satisfied for the proper work of the algorithm\cite{PoyneerFastFFTreconWFR}:

1. the gradients are spatially periodic (necessary for use of the DFT method, and it must be maintained for a set of gradient measurements). Check if the sum of every row or column in the $N\times N$ gradient signal equals zero.
2. any closed path of gradients must sum to zero (based on the modeling of the gradients as first differences).
   

Hence the two conditions must be satisfied for correct reconstruction:

1. all loops (under Hudgin or Fried geometry) must sum to zero;
2. both slope signals must be spatially periodic (for DFT)
   

This spatial periodicity assumption does not require that the sensed wave front be inherently periodic\cite{PoyneerFastFFTreconWFR}. Like VMM reconstructors work while estimating only the phase values inside the aperture, FT methods accurately reconstruct on grids of aperture size.

The boundary problem is solved by using specific methods to make the gradient
sets consistent\cite{poyneer2003advanced}. These methods involve setting the values of specific gradients outside the aperture in a way to guarantee correct reconstruction inside the aperture. By using these methods FT reconstructors accurately reconstruct all sensed modes inside the aperture.


Solutions to the boundary problem: boundary method
The first method is the boundary method: it estimates the gradients that cross the boundary of the aperture\cite{PoyneerFastFFTreconWFR}. It follows directly from the development
above of inside, boundary, and outside gradients. Only the inside gradients are known from the measurement. The outside gradients can all be set to zero leaving the boundary gradients undetermined. A loop continuity equation can be written for each of the two smallest loops that involve a boundary gradient.

Boundary method. Setting each closed loop across the aperture edge to zero results in an equation relating the unknown boundary gradients to the measured inside gradients and the zeroed outside gradients [picture from the paper \cite{PoyneerFastFFTreconWFR}].

All of these loop continuity equations involving the boundary gradients combine to form a linear system:

$\mathbf{M}u = c$

Here the matrix $M$ is specific for the selected geometry, the $u$ is the vector of all boundary gradients, and the $c$ is a vector containing sums of measured gradients (it has a fixed combination of gradients, but the value of these gradients depends on the actual measurement\cite{PoyneerFastFFTreconWFR}). In the case of noised measurements, the solution can be found using the pseudoinverse of the matrix $M$.


Solutions to the boundary problem: Extention method
The second method is the extension method: it extends out the gradients from inside the aperture. The extension method extends the wavefront shape to
outside the aperture. It does this by a preserving loop continuity\cite{PoyneerFastFFTreconWFR}. The x gradients are extended up and down out of the aperture, while the y gradients are extended to the left and the right.

Extension method, shown for $N=6$. The values of the gradients closest to the aperture edge are repeated outside the aperture [picture from the paper \cite{PoyneerFastFFTreconWFR}].

The extension method produces a completely consistent set of gradients. These provide perfect reconstruction of the phase when there is no noise, except for the piston. If there is noise, the same procedure is done, though loop continuity will not hold on loops involving the seam gradients, just as the boundary gradients in the boundary method were the best, but not an exact, solution when there was noise\cite{PoyneerFastFFTreconWFR}.

References:

\begin{thebibliography}{1}

\bibitem{PoyneerFastFFTreconWFR} Lisa~A. Poyneer, Donald~T. Gavel, and James~M. Brase. \newblock Fast wave-front reconstruction in large adaptive optics systems with use of the fourier transform. \newblock {\em J. Opt. Soc. Am. A}, 19(10):2100--2111, Oct 2002.

\bibitem{poyneer2003advanced} L.A. Poyneer. \newblock {Advanced techniques for Fourier transform wavefront reconstruction}. \newblock In {\em Proceedings of SPIE}, volume 4839, page 1023, 2003.

\end{thebibliography}

Modified Hudgin geometry for Wavefront reconstruction using Fast Fourier Transform

One of the attractive abilities of the FFT WFR is that the filtering construct provides flexibility: reconstruction is accomplished by filtering in the frequency domain and one can modify this filter with negligible computational overhead. It is easy to incorporate filtering options into the
reconstruction filters like noise reduction, modal removal, misalignment, or DM geometry compensation.

Misalignment of the WFS data and the DM geometry.

For example, the misalignment can be compensated: the WFS grid and the DM actuators may be misaligned by shifts along x or y dimensions. If the amount is known, shift slope estimate by a fraction of an actuator spacing such as $(exp[-i 2\pi (\Delta_x k + \Delta_y l)/ N])\Phi$.


Modified Hudgin geometry
For Shack-Hartmann, best the reconstructor is modified Hudgin, for which the
slopes are: $s_x[m,n] = \phi[m+1n,n+0.5] - \phi[m,n+0.5]$.
The reconstruction of the wavefront is performed with the following procedure:
$\begin{matrix} \hat{\Phi}[k,l] & = & \left\{ \begin{matrix} 0 & \mbox{if } k,l = 0 \\ S_x[k,l]\cdot H_x[k,l] + S_y[k,l]\cdot H_y[k,l] & \mbox{otherwise } \end{matrix} \right.\end{matrix}$

The the spatial filters are different: we ahve to shift each slope signal half a
sample along the orthogonal direction:
$ H_x[k,l] = \frac{ (\exp\left[ - \frac{2\pi i \cdot k }{N} \right] -1) \exp(-i\pi
l / N) }{ 4\left( \sin^2\frac{\pi k}{N} + \sin^2\frac{\pi l}{N} \right)} \\ H_y[k,l] = \frac{ (\exp\left[ - \frac{2\pi i \cdot l }{N} \right] -1) \exp(-i\pi
k / N) }{ 4\left( \sin^2\frac{\pi k}{N} + \sin^2\frac{\pi l}{N} \right)} $

Taking the inverse transform $\mathcal{F}^{-1}\hat{\Phi}[k,l] $ produces the
estimate of the wavefront $\hat{\phi}[m, n]$. Such a geometry estimations (see Fig.~\ref{fig:ao2004_3670_modified_Hudgin}) are of
high quality, and it does not suffer from global or local waffle like Fried
geometry.

Modified Hudgin geometry.

That was validated\cite{poyneer2003experimental} in on-sky testing at
Palomar Observatory. In the experiments, the authors of
\cite{poyneer2003experimental} tried out several options for geometries and
filtering, and the modified Hudgin performed best. It is noteworthy that the
regular Hudgin geometry suffered from misalignment-like errors. Finally, the
Fried geometry had excessive local waffle.

Moreover, the modified Hudgin takes half as much computation as the Fried geometry model.

Limitations and Disadvantages of the FTR

Disadvantages
First, if the aperture size in sub-apertures is not a power-of-2, that can cause
performance problems: extensive padding to get to a power-of-2 leads to
increased noise. The FTR requires the square or the pseudo-hex DM geometry. The
non-integer ratio of sub-apertures size or actuator spacing
requires correct re-sampling of estimate.

Advantages
The FTR is fast enough for the ExAO systems and large simulation codes. It provides adaptability with filtering - one can compensate the misalignment and
other errors. The modified Hudgin method does not suffer significantly from the global or local waffle.


The problems in implementation of modified Hudgin
The modified Hudgin still requires the suppression of the piston mode:

if ((rownum == 0)&(colnum == 0))
wavefront(rownum+1, colnum+1) = 0;
else

The original formula for the spatial filters works queer. The additional shift must be like this: $\exp(-i2\pi l / N)$ instead of simply $\exp(-i\pi l / N)$.

Moreover, the code for the different shifts gives slightly different results. The results with the spatial filter:

$$H_x[k,l] = \frac{ (\exp\left[ - \frac{2\pi i \cdot k }{N} \right] -1) \exp(-i2\pi l / N) }{ 4\left( \sin^2\frac{\pi k}{N} + \sin^2\frac{\pi l}{N} \right)} $$

and implemented as a MATLAB code:

H_row = (exp(-(i*2*rownum*pi)/N_row) -1)*(exp(-(i*2*colnum*pi)/N_row)); %% spatial filters for X axis that include the complex conjugate of exponentials
H_column = (exp(-(i*2*colnum*pi)/N_col) -1)*(exp(-(i*2*rownum*pi)/N_col)); %% spatial filters for Y axis that include the complex conjugate of exponentials

gives the reconstruction, shown in Fig.~\ref{fig:wfr_modified_Hudgin_row-col}a. Contrary, the spatial filter of:

$$H_x[k,l] = \frac{ (\exp\left[ - \frac{2\pi i \cdot l }{N} \right] -1) \exp(-i2\pi
k / N) }{ 4\left( \sin^2\frac{\pi k}{N} + \sin^2\frac{\pi l}{N} \right)} $$

with code:

H_row = (exp(-(i*2*rownum*pi)/N_row) -1)*(exp(-(i*2*colnum*pi)/N_row)); %% spatial filters for X axis that include the complex conjugate of exponentials
H_column = (exp(-(i*2*colnum*pi)/N_col) -1)*(exp(-(i*2*rownum*pi)/N_col)); %% spatial filters for Y axis that include the complex conjugate of exponentials

gives better result in the reconstruction, as seen in Fig.~\ref{fig:wfr_modified_Hudgin_row-col}b.

The results of reconstruction for the modified Hudgin algorithm for the row spatial filter: a) row-column relationships with shifts, and b) column-row relationships.

The results of the modified Hudgin geometry are indeed better compared with the standard Hudgin geometry (see the reconstruction with conventional Hudgin FFT WFR (below)

The results of reconstruction with conventional FFT Hudgin.