© ESO, 2011

1. Introduction

Speckle imaging techniques have been evolving since the fundamental idea was presented almost forty years ago (Labeyrie 1970; Labeyrie et al. 2006; Knox & Thompson 1974; Weigelt 1977). Critically, the early work demonstrated how moments of the Fourier transform of an ensemble of short exposures contain information out to the diffraction limit. Many variations on this theme have been implemented, but in all cases, an ensemble of short-exposure images is collected and then post-processed to restore the object. In this presentation, we will develop object spectrum phasor reconstruction methods. We concentrate on phasor-based phase reconstruction algorithms in order to avoid problems associated with branch points in the object spectrums (Dente 2000; Ayers et al. 1988). In each case, by minimizing a phasor-based error function, we develop iterative improvement algorithms that converge to the phasor array that minimizes the error function; we call this the least-squares-best (LSB) two-dimensional phasor array for each object spectrum. In addition, the optimum LSB algorithms automatically preserve the unit modulus and Hermitian symmetry of each updated phasor array on each iteration. These phasor reconstruction algorithms offer improvements to prior work, even when viewed against the backdrop of considerable prior work on phasor-based reconstructions.

Most prior work approximately falls into two classes. In the first class, researchers took least-squares or weighted-least-squares methods that reconstructed phase from phase differences and then converted the methods to phasor-based algorithms (Ayers et al. 1988; Meng et al. 1990; Matson 1991). Specifically, phase sums or differences were converted to phasor multiply or divide operations. These can be practical, intuitive methods for making phase reconstruction algorithms tolerant to branch points. However, the methods remain suspect since they have not been developed by minimizing a phasor-based error function.

In the second class, researchers started with a phasor-based error function. However, instead of directly solving for the minimum of the error, they performed a gradient-guided search in the high-dimensional space spanned by the object spectrum phase elements (Negrete-Regagnon 1996). Done this way, there are myriad local minima en route to the true minimum where all the gradients of the error function are zero. Searching for the minimum error in this way can be a poor choice of algorithm that can easily stagnate (Negrete-Regagnon 1996).

Juxtaposed to these methods, our algorithms directly solve for the phasor array that minimizes a phasor-based error function. Once we choose an error function, the remaining steps are systematically developed by setting the error function variation, with respect to each phasor element, to zero. The resulting coupled nonlinear equations for the minimum-error phasor array are then solved iteratively. There are no issues with local minima or other search complexities in the high-dimensional space of phasor elements. Once we have selected the error function, our converged solution for the phasor array represents a true minimum for this error function; it cannot be improved. The only element of subjectivity, the only element of doubt, is the original choice for the phasor-based error function.

In our simulations, we will implement three speckle methods on a simple object at low photon-per-frame light levels. Then, we will apply the methods to a complex extended object, once again varying the photons per frame. We assume that the only aberrations are those introduced by atmospheric turbulence setting the ratio of the telescope diameter, D, to the Fried Parameter, r₀, equal to ten (Fried 1966). Additionally, we assume that the focal-plane detector array is photon-noise limited, the illumination is narrow-band (essentially monochromatic) and the atmosphere is static during each data frame. In all cases, we apply each of the methods to aberrated image ensembles of two hundred short-exposure data frames. Finally, we will describe the strengths and weaknesses of each of the three methods applied to both simple and extended objects.

2. Reconstructing the object spectrum modulus

We label the arrays of short-exposure data as (1)We will assume that the short exposures have all been centered with gross displacements removed. The index l = 1,2,...   R labels the frame, while the next two indices, (i,j), label the pixel in the array. We take a two-dimensional finite-Fourier-transform (FFT) on the pixel indices as (2)The array indices, (m,n), label the elements of the two-dimensional FFT. The spatial frequency indices are rearranged so that the zero spatial frequency is centered at the point (m = 0,   n = 0). We then form the unbiased average (Dente 2000) (3)in which is the modulus of the object spectrum (FFT), and is the speckle transfer function; this can be calculated with a reference star or even analytically. A standard Wiener-Helstrom deconvolution filter then yields an estimate of the modulus-squared of the object spectrum. In the following comparisons, we use the same modulus for each of the speckle methods; therefore, they only differ in the phasor reconstruction methods.

3. Reconstructing the object spectrum phase

The object spectrum, essentially the FFT of the object, is given by the modulus times the phase factor as (4)in which is the object spectrum phase array and is defined as the object spectrum phasor array. In what follows, we will develop and compare methods that use a Knox-Thompson or bispectrum data base to generate phasor estimates. In each case, prior to applying an inverse FFT to the reconstructed object spectrum, we apply the diffraction-limited optical transfer function for the telescope; this avoids signal-to-noise problems at the higher spatial frequencies.

In the Knox-Thompson method, called KT, we start by forming the unbiased cross-spectrums from an ensemble of R short exposures: (5)Here, we have picked one unit offset in each direction in the FFT; other offsets can easily be analyzed, but the unit offset choice gives the best accuracy on the high-spatial-frequency phase. The final terms remove the expectation of the photon noise. The ensemble averaging over many realizations of the turbulence-induced phase screens makes the cross-spectrum transfer functions real-valued (Knox & Thompson 1974). Therefore, the turbulence-induced phase errors are averaged away, and the phase of the cross-spectrum can be directly related to phase differences in the object spectrum (Knox & Thompson 1974; Dente 2000). The remaining phases of these two cross-spectrums lead to two sets of equations for the phasor array of the object spectrum given by: (6)For the least-squares-best KT method, we solve for the phasor array by minimizing the least-squares error function defined as (7)This error function follows naturally from the cross-spectrum equations, Eq. (6). However, note that the complex conjugate phasor factor has been replaced with the inverse. For a unit modulus complex number, a phasor, the complex conjugate operation and the inverse are identical, so that . By using the equivalent inverse form in Eq. (7), we develop an error function that only depends on the phasor array elements and carries no dependence on the conjugate phasors. With this form, when we set the variation of the error with respect to each phasor element to zero, we obtain N² equations that automatically incorporate the phasor constraint of unit modulus. If we had instead used the complex conjugate form, we would have had to consider independent variations of the error functions with respect to each phasor and conjugate phasor element. These 2   N² equations, along with the N² unit modulus constraints, would eventually lead to the same final equations for the least-squares-best phasor array, Eq. (10), below. By incorporating the P^-1 = P^∗ form directly into the error function, we can derive the correct minimum error equations much more directly. This same method is repeated in all of the error functions.

Equation (7) is a subjective element in the procedure. Certainly other error functions could be generated, but it seems difficult to develop a more natural choice. To proceed, we set the variation of the error with respect to each equal to zero. The resulting equations that locate the error function minimum, can be put in the form (8a)Each of these derivatives, for each array element, can be rearranged to N² equations (8b)in which (9)Equation (8b) is solved by (10)We see that each element of the optimized phasor array is determined by the four nearest-neighbor array elements, with each term weighted by a connecting phasor and the squared modulus of a cross-spectrum. The influence of nearest neighbors is similar to that found in phase reconstruction algorithms for phase difference data sets, but in this case, automatic weighting of each nearest neighbor occurs (Hudgin 1977). These coupled nonlinear equations, Eq. (10), must be simultaneously solved, yielding a least-squares-best phasor array, , in which the unit modulus of the phasor is automatically preserved on each iteration. In addition, the algorithm preserves the Hermitian conjugation property of the object spectrum phasor array, so that .

An iterative improvement solution to Eq. (10) works quite well. The updated phasor estimates for all array elements, , are calculated by substituting the previous phasor estimates, , into the right-hand side of Eq. (10). The iterative convergence can be considerably improved if we relax the phasor towards its updated value as (11)with the new phasor estimate averaged with the old phasor estimate and then renormalized; then replaces the old phasor array. We typically use a weighting factor ω = .25 ↔ .5. We initialize the iterations by setting all and then stop the iterations when the array elements are stabilized and unchanging; these converged phasor elements have forced the variations of the error function to zero and no further improvements are possible. This least-squares-best algorithm is graphically illustrated in Fig. 1a.

Fig. 1 a) Graphical form of methods KT and TCA. b) Graphical form of method TCB.

In the two bispectrum, or Triple-Correlation methods, referred to as TCA and TCB, we start by forming several sets of the short-exposure ensemble averages defined by (12)in which range over all spatial frequencies in the object spectrum, while (m′,n′) are fixed array elements. Each value of the pair (m′,n′) defines a bispectrum plane. The “noise terms” refer to the photon noise terms that remove the expectation of the photon noise (Ayers et al. 1988).

In TCA, we choose two bispectrum planes defined by m′ = 1,n′ = 0 and m′ = 0,n′ = 1. This choice is similar to the unit shift used in Knox-Thompson. Once again, the object spectrum phasor array is taken to minimize the least-squares phasor error defined by (13)Again, we use the phasor condition, , to obtain an error that only depends on the phasor array. While this choice of error function is subjective, it is difficult to make a more natural choice. Here again, the variation, with respect to , leads to a nonlinear system of equations that can be put in the form of Eq. (8) and then solved by Eq. (10). In this case, we find the complex factor on the right-hand side as (14)These optimum phasor equations, locating the error function minimum, can then be iteratively solved to convergence, yielding the least-squares-best phasor array, . Previous estimates generate updated estimates according to Eqs. (10) and (11). The iterations are stopped when the phasor updates are converged. This weighted least-squares-best algorithm is graphically illustrated in Fig. 1a.

Our third implementation, method TCB, is based on four closely-spaced bispectrum planes with (m′ = 1,n′ = 0), (m′ = 0,n′ = 1), (m′ = 2,n′ = 0), and (m′ = 0,n′ = 2). Here, the least-squares phasor error function depends on four bispectrums as (15)Here again, the least-squares-best values for the phasor array results from taking the derivative of this error with respect to and setting the result to zero. This leads to an iterative procedure that minimizes this error. The final weighted least-squares-best algorithm is graphically illustrated in Fig. 1b.

We can consider generalizations to even more bispectrum planes, by defining an error function per plane as (16)We then sum over any set of bispectrum planes, , in order to set a total error. In this notation, the error function for the TCA method is , while the error function for the TCB method is . We explicitly show in the second line that each phasor error term is multiplied by a weighting, the squared modulus of the bispectrum.

4. Reconstruction results

In the simulations, we will assume that the only aberrations are those introduced by atmospheric turbulence, setting the ratio of the telescope diameter, D, to the Fried Parameter, r₀, equal to ten (Fried 1966). We could also generalize the simulations to include static and dynamic telescope aberrations. Additionally, we assume that the focal-plane detector array is photon-noise limited, so that we neglect dark-current, as well as read noise. Also, we only process the data on (64 × 64) pixel elements. We fix the final (F/#) of the telescope so that the diffraction-limited point-spread-function is approximately two to three pixels in diameter. In all cases, we used the same speckle-derived modulus, while each phasor reconstructor used the algorithms described for KT, TCA and TCB in the last section. In each case, we applied the diffraction-limited optical transfer function for the telescope to the object spectrum reconstruction prior to applying an inverse FFT.

In the examples that follow, convergence to the final phasor array required approximately two thousand iterations of the N² coupled nonlinear equations given by Eq. (10). We always continue the self-consistent iterations until the phasor array elements are converged and unchanging; this point corresponds to the true minimum of the error function. Since these minimum error phasor equations only couple nearest neighbors (or next nearest neighbors for TCB), it takes many iterations for all of the phasor array elements to be self-consistently determined so that each phasor element is properly influenced by every other phasor element. The myriad paths that connect any two points in the two-dimensional spatial-frequency space are all automatically included during the iterations, while each of these paths automatically carries the proper weighting (Typically, the image reconstructions for method TCA take less than one minute of computation time on a personal computer.).

Fig. 2 a) Asterism, single frame diffraction aberrated image. b) Asterism, long term aberrated image. c) Knox-Thompson method KT. ΔP = 1.2e-2. d) Bispectrum method TCA. ΔP = 1.3e-2. e) Bispectrum method TCB. ΔP = 3.8e-2.

First, we implemented the methods on a simple four-star asterism at low photon-per-frame light levels. For each case, we used an ensemble of two hundred short-exposure frames with five hundred photons per frame. Figures 2a, b show a single frame diffraction-limited asterism image with photon noise, as well as the long-term aberrated image. The reconstruction results are shown in Figs. 2c–e. In this case, all of the methods offered acceptable reconstructions.

Fig. 3 a) Extended object, single frame diffraction limited image. b) Extended object, long term aberrated limited image. c) Knox-Thompson method KT. ΔP = 1.6e-3. d) Bispectrum method TCA. ΔP = 0.9e-3. e) Bispectrum method TCB. ΔP = 3.2e-3.

Next, we reconstructed a more complicated and extended object at higher light levels. We used an ensemble of two hundred short-exposure frames, but increased the photons per frame to ten thousand. Figures 3a, b show a single-frame, diffraction-limited image with photon noise, as well as the long-term aberrated image. The reconstruction results are shown in Figs. 3c–e. For this extended object, the differences between the three speckle reconstruction approaches are more pronounced, with KT and the closely-related TCA giving the best results and TCB the worst.

5. Discussion and conclusion

The phasor-based error functions, Eqs. (7), (13), (15) and (16) lead to a new set of iterative phasor reconstruction algorithms that are least-squares-best. Therefore, the final converged phasor array, , is located at the minimum of each error function. In addition, the algorithms automatically preserve the unit modulus, as well as the Hermitian symmetry property of the object spectrum phasor array, so that and on every iteration. The examples demonstrated reconstructions on 200 short-exposure frames in the presence of significant turbulence and photon noise. For the simple object case shown in Figs. 2a–e, all of the methods, KT, TCA, and TCB work well. For the extended object shown in Figs. 3a–e, our method TCB, using four bispectrum planes, gave poorer reconstructions than either KT or TCA.

To further quantify the reconstruction errors, we calculated a phasor-based error function that measured the average error between the raw data phase differences and the phasor reconstructions. For each method, a mean squared phasor error was developed from the cross-spectrums as (17)The cross-spectrum powers normalize the terms so that this gives an estimate of the squared error per phasor array element. We allowed for the possibility of constant phase errors, α and β, between the raw data and the reconstructions and always minimized the error function with respect to these two parameters; this removes the penalty for slight tilts in the reconstructed object spectrum. The same error function, using the array products and , was also evaluated for methods TCA and TCB. In all cases, as the number of photons per frame increases, these error functions decrease and approach zero. The figure captions for the reconstructions list this phasor error. By this error measure, TCB generated the largest errors in these examples.

These preceding observations, indicating inferior image reconstructions for TCB, especially when the object is an extended object, demand an explanation. First, we emphasize that these observations only refer to the phasor reconstruction algorithms presented here. Other reconstruction algorithms, as developed over the last twenty-five years, do not share this property. Those alternative algorithms do seem to benefit from using a larger number of bispectrum planes (Lohmann et al. 1983; Hofmann & Weigelt 1995). Currently, we do not fully understand why our minimum least-squares reconstructions do not improve, and even degrade upon increasing to four rather than two bispectrum planes. One possible explanation follows: since the nonlinear least-squares equations, Eqs. (10) and (14) for method TCA, always couple a phasor element to its four nearest neighbors, the multiple iterations en route to convergence ultimately link a single phasor element anywhere in the array to all other phasor elements. This iterative convergence to self-consistency could circumvent the need for more than two bispectrum planes. Alternatively, our extended object example showing poorer performance for TCB could be an isolated effect that is peculiar to our choice of extended object, as well as our limit to just four bispectrum planes.

Although the phasor-based error functions, Eqs. (7), (13), (15) and (16), certainly have a strong intuitive appeal, it is always possible that better error functions could be derived. Meanwhile, the phasor-based reconstruction algorithms that result from minimizing these error functions remain least-squares-best and rapidly converging, while preserving the unit modulus, as well as the Hermitian symmetry of the object spectrum phasor array on every iteration; these are powerful features in their favor.

References

All Figures

	Fig. 1 a) Graphical form of methods KT and TCA. b) Graphical form of method TCB.
In the text

	Fig. 2 a) Asterism, single frame diffraction aberrated image. b) Asterism, long term aberrated image. c) Knox-Thompson method KT. ΔP = 1.2e-2. d) Bispectrum method TCA. ΔP = 1.3e-2. e) Bispectrum method TCB. ΔP = 3.8e-2.
In the text

	Fig. 3 a) Extended object, single frame diffraction limited image. b) Extended object, long term aberrated limited image. c) Knox-Thompson method KT. ΔP = 1.6e-3. d) Bispectrum method TCA. ΔP = 0.9e-3. e) Bispectrum method TCB. ΔP = 3.2e-3.
In the text