Issue 
A&A
Volume 529, May 2011



Article Number  A58  
Number of page(s)  7  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/201016335  
Published online  01 April 2011 
New phasor reconstruction for speckle imaging
^{1} GCD Associates,
2100 Alvarado NE, Albuquerque, NM
87110,
USA
^{2}
Boeing LTS, PO
Box 5670, Kirtland
AFB, NM
87185,
USA
email: Michael.Tilton@kirtland.af.mil
Received:
15
December
2010
Accepted:
2
March
2011
We will develop and then compare object spectrum phasor reconstruction results for several speckle imaging approaches. Each phasor reconstruction algorithm results from minimizing a very naturally defined weightedleastsquares error function. Once we pick a phasorbased error function, the remaining steps in our algorithms are 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. In the applications, we will specifically compare and contrast three implementations: 1) KnoxThompson; 2) bispectrum, using only two bispectrum planes; 3) bispectrum, using four bispectrum planes. In each application of the three approaches, we first calculate the modulus of the object spectrum using a WienerHelstrom filter to remove the speckle transfer function. The methods then differ only in their object spectrum phasor reconstructions. In the simulations, we will implement all three methods on a simple object at low photonperframe light levels. Next, we will apply the methods to a complex extended object. Although we develop and minimize error functions for three specific speckle methods, the approach readily generalizes to other cases.
Key words: techniques: image processing / methods: analytical / methods: data analysis / atmospheric effects
© 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 shortexposure images is collected and then postprocessed to restore the object. In this presentation, we will develop object spectrum phasor reconstruction methods. We concentrate on phasorbased 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 phasorbased error function, we develop iterative improvement algorithms that converge to the phasor array that minimizes the error function; we call this the leastsquaresbest (LSB) twodimensional 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 phasorbased reconstructions.
Most prior work approximately falls into two classes. In the first class, researchers took leastsquares or weightedleastsquares methods that reconstructed phase from phase differences and then converted the methods to phasorbased 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 phasorbased error function.
In the second class, researchers started with a phasorbased error function. However, instead of directly solving for the minimum of the error, they performed a gradientguided search in the highdimensional space spanned by the object spectrum phase elements (NegreteRegagnon 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 (NegreteRegagnon 1996).
Juxtaposed to these methods, our algorithms directly solve for the phasor array that minimizes a phasorbased 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 minimumerror phasor array are then solved iteratively. There are no issues with local minima or other search complexities in the highdimensional 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 phasorbased error function.
In our simulations, we will implement three speckle methods on a simple object at low photonperframe 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_{0}, equal to ten (Fried 1966). Additionally, we assume that the focalplane detector array is photonnoise limited, the illumination is narrowband (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 shortexposure 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 shortexposure 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 twodimensional finiteFouriertransform (FFT) on the pixel indices as (2)The array indices, (m,n), label the elements of the twodimensional 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 WienerHelstrom deconvolution filter then yields an estimate of the modulussquared 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 KnoxThompson 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 diffractionlimited optical transfer function for the telescope; this avoids signaltonoise problems at the higher spatial frequencies.
In the KnoxThompson method, called KT, we start by forming the unbiased crossspectrums 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 highspatialfrequency phase. The final terms remove the expectation of the photon noise. The ensemble averaging over many realizations of the turbulenceinduced phase screens makes the crossspectrum transfer functions realvalued (Knox & Thompson 1974). Therefore, the turbulenceinduced phase errors are averaged away, and the phase of the crossspectrum can be directly related to phase differences in the object spectrum (Knox & Thompson 1974; Dente 2000). The remaining phases of these two crossspectrums lead to two sets of equations for the phasor array of the object spectrum given by: (6)For the leastsquaresbest KT method, we solve for the phasor array by minimizing the leastsquares error function defined as (7)This error function follows naturally from the crossspectrum 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^{2} 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^{2} equations, along with the N^{2} unit modulus constraints, would eventually lead to the same final equations for the leastsquaresbest 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^{2} 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 nearestneighbor array elements, with each term weighted by a connecting phasor and the squared modulus of a crossspectrum. 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 leastsquaresbest 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 righthand 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 leastsquaresbest 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 TripleCorrelation methods, referred to as TCA and TCB, we start by forming several sets of the shortexposure 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 KnoxThompson. Once again, the object spectrum phasor array is taken to minimize the leastsquares 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 righthand side as (14)These optimum phasor equations, locating the error function minimum, can then be iteratively solved to convergence, yielding the leastsquaresbest 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 leastsquaresbest algorithm is graphically illustrated in Fig. 1a.
Our third implementation, method TCB, is based on four closelyspaced bispectrum planes with (m′ = 1,n′ = 0), (m′ = 0,n′ = 1), (m′ = 2,n′ = 0), and (m′ = 0,n′ = 2). Here, the leastsquares phasor error function depends on four bispectrums as (15)Here again, the leastsquaresbest 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 leastsquaresbest 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_{0}, equal to ten (Fried 1966). We could also generalize the simulations to include static and dynamic telescope aberrations. Additionally, we assume that the focalplane detector array is photonnoise limited, so that we neglect darkcurrent, 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 diffractionlimited pointspreadfunction is approximately two to three pixels in diameter. In all cases, we used the same specklederived modulus, while each phasor reconstructor used the algorithms described for KT, TCA and TCB in the last section. In each case, we applied the diffractionlimited 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^{2} coupled nonlinear equations given by Eq. (10). We always continue the selfconsistent 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 selfconsistently determined so that each phasor element is properly influenced by every other phasor element. The myriad paths that connect any two points in the twodimensional spatialfrequency 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) KnoxThompson method KT. ΔP = 1.2e2. d) Bispectrum method TCA. ΔP = 1.3e2. e) Bispectrum method TCB. ΔP = 3.8e2. 
First, we implemented the methods on a simple fourstar asterism at low photonperframe light levels. For each case, we used an ensemble of two hundred shortexposure frames with five hundred photons per frame. Figures 2a, b show a single frame diffractionlimited asterism image with photon noise, as well as the longterm 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) KnoxThompson method KT. ΔP = 1.6e3. d) Bispectrum method TCA. ΔP = 0.9e3. e) Bispectrum method TCB. ΔP = 3.2e3. 
Next, we reconstructed a more complicated and extended object at higher light levels. We used an ensemble of two hundred shortexposure frames, but increased the photons per frame to ten thousand. Figures 3a, b show a singleframe, diffractionlimited image with photon noise, as well as the longterm 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 closelyrelated TCA giving the best results and TCB the worst.
5. Discussion and conclusion
The phasorbased error functions, Eqs. (7), (13), (15) and (16) lead to a new set of iterative phasor reconstruction algorithms that are leastsquaresbest. 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 shortexposure 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 phasorbased 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 crossspectrums as (17)The crossspectrum 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 twentyfive 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 leastsquares reconstructions do not improve, and even degrade upon increasing to four rather than two bispectrum planes. One possible explanation follows: since the nonlinear leastsquares 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 selfconsistency 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 phasorbased 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 phasorbased reconstruction algorithms that result from minimizing these error functions remain leastsquaresbest 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
 Ayers, G. R., Northcott, M. J., & Dainty, J. C. 1988, J. Opt. Soc. Am., 5, 963 [Google Scholar]
 Dente, G. C. 2000, Appl. Opt., 39, 1480 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Fried, D. L. 1966, J. Opt. Soc. Am., A56, 1372 [NASA ADS] [CrossRef] [Google Scholar]
 Hofmann, K. H., & Weigelt, G. 1995, A&A, 300, 403 [NASA ADS] [Google Scholar]
 Hudgin, R. H. 1977, J. Opt. Soc. Am., 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Knox, K. T., & Thompson, B. J. 1974, AJ, 193, L45 [Google Scholar]
 Labeyrie, A. 1970, A&A, 6, 85 [NASA ADS] [Google Scholar]
 Labeyrie, A., Lipson, S. G., & Nisenson, P. 2006 (Cambridge University Press) [Google Scholar]
 Lohmann, A. W., Weigelt, G., & Wimitzer, B. 1983, Appl. Opt., 22, 4028 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Matson, C. L. 1991, J. Opt. Soc. Am., A8, 1905 [Google Scholar]
 Meng, J., Aitken, G., Hege, E., & Morgan, J. 1990, J. Opt. Soc. Am., A7, 1243 [Google Scholar]
 NegreteRegagnon, P. 1996, J. Opt. Soc. Am., A13, 1557 [NASA ADS] [CrossRef] [Google Scholar]
 Weigelt, G. 1977, Opt. Comm., 21, 55 [Google Scholar]
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) KnoxThompson method KT. ΔP = 1.2e2. d) Bispectrum method TCA. ΔP = 1.3e2. e) Bispectrum method TCB. ΔP = 3.8e2. 

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

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.