Comparison of bispectrum, multiframe blind deconvolution and hybrid bispectrum-multiframe blind deconvolution image reconstruction techniques for anisoplanatic, long horizontal-path imaging

Abstract. The potential benefits of real-time, or near-real-time, image processing hardware to correct for turbulence-induced image defects for long-range surveillance and weapons targeting are sufficient to motivate significant resource commitment to their development. Quantitative comparisons between potential candidates are necessary to decide on a preferred processing algorithm. We begin by comparing the mean-square-error (MSE) performance of speckle imaging (SI) methods and multiframe blind deconvolution (MFBD), applied to long-path horizontal imaging of a static scene under anisoplanatic seeing conditions. Both methods are used to reconstruct a scene from three sets of 1000 simulated images featuring low, moderate, and severe turbulence-induced aberrations. The comparison shows that SI techniques can reduce the MSE up to 47%, using 15 input frames under daytime conditions. The MFBD method provides up to 40% improvement in MSE under the same conditions. The performance comparison is repeated under three diminishing light conditions, 30, 15, 8 photons per pixel on average, where improvements of up to 39% can be achieved using SI methods with 25 input frames, and up to 38% for the MFBD method using 150 input frames. The MFBD estimator is applied to three sets of field data and representative results presented. Finally, the performance of a hybrid bispectrum-MFBD estimator that uses a rapid bispectrum estimate as the starting point for the MFBD image reconstruction algorithm is examined.


Introduction
There is interest in the development of human-portable surveillance systems capable of observing a wide field of view over long horizontal or near-horizontal paths. In addition, remotely piloted surveillance systems must transmit images to remote users over wireless channels that have limited bandwidth. As a result, currently fielded surveillance systems compress images prior to transmission. 1 An image processing system that could reconstruct images corrupted by turbulence prior to compression and transmission would allow such systems to make more efficient use of available bandwidth. Thus, effort must be directed toward the examination and selection of computationally efficient and robust image reconstruction techniques. Under all but the most benign conditions, temperature inhomogeneities in the atmosphere result in turbulence and variations in the index of refraction along the imaging path. Even small variations in the index of refraction cause changes in the optical path length that result in phase aberrations at the aperture. 2 As the optical path length increases, or the turbulence strength increases, the aberrations become stronger, and the isoplanatic angle decreases. Terrestrial surveillance scenarios often involve imaging scenes that subtend angles much larger than the isoplanatic angle, a situation referred to as anisoplanatic imaging. Under these conditions, increasing the size of the aperture will not improve the quality of the image. Unless the seeing conditions are very favorable, anisoplanatism will dominate most practical horizontal imaging situations. The resulting distortion will limit the performance of any optical system operating in such a turbulent atmosphere, frequently causing the observed scene to be blurred beyond usefulness. 3 Many strategies have been proposed to estimate the true scene from turbulence-corrupted images, not all of which are explored in this work. Each represents a compromise between how quickly the technique produces an estimate of the scene, the accuracy of the estimate, and the robustness of the technique. Quick and accurate convergence under a wide range of atmospheric seeing conditions, regardless of the scene illumination, while remaining transportable by a single operator are desirable characteristics. Power and portability constraints may eliminate traditional adaptive optics approaches that have been applied to similar problems in the imaging of space objects 4 using currently available technology. An additional limitation to adaptive optic systems is that they only correct well over a few multiples of the isoplanatic angle, while common horizontal surveillance scenarios observe scenes that extend over many multiples of the isoplanatic angle.
Ayers and Dainty, 5 among others, and Lohmann et al. 6 determined that near diffraction-limited phase information could be recovered from multiple-speckled images through the use of Knox-Thompson and bispectrum techniques. 7 The recovered phase and intensity information are used to reconstruct nearly diffraction-limited images of astronomical objects. Carrano 8 proposed the use of speckle imaging (SI) for horizontal imaging and EM Photonics, (Newark, Delaware) has developed an embedded system using the technique to reconstruct turbulent images over horizontal paths in near real time. 9 Bos and Roggemann have extended this technique for use over wide fields of view and also explored the performance limits of SI over a range of turbulence levels. Both Carrano and Bos ignore anisoplanatism in their processing but nevertheless produce reconstructions with reduced error. 3 Other blind deconvolution techniques have also been applied to reconstructing turbulence-corrupted images. Richardson 10 and Lucy 11 pioneered the application of an iterative blind deconvolution technique to a single image degraded by atmospheric turbulence. Their method relied on the enforcement of positivity, and finite support constraints on the object estimate. Schulz 12 extended that method to include multiple input images and developed a penalized-maximum-likelihood algorithm to avoid the trivial solution that incorrectly concludes that the optical system's point spread function (PSF) is a Dirac delta function and that the most likely object estimate is the observed turbulent image. Little research has been done to compare the performance characteristics of the multi-frame blind deconvolution (MFBD) and SI techniques under anisoplanatic conditions; that is the subject of this work.
This article compares the performance of SI techniques and MFBD at estimating object intensities from simulated images exhibiting extreme levels of anisoplanatism, using the number of input frames and scene illumination levels as the operational parameters. To enable this comparison, the authors developed a horizontal imaging simulation model and used it to create three sets of 1000 simulated, turbulent images based on the image processing Lena image, common in the literature. In order to explore the technique's performance under multiple illumination conditions, each set of simulated turbulent images was used to generate a set of speckle images with a mean photon count of 30, 15, and 7 photons per pixel. The effectiveness of SI and MFBD is assessed by calculating the mean-square-error (MSE) between the resulting object estimate and the diffraction limited image.
We find that the both SI-and MFBD-reconstructed objects show significant improvement in MSE compared with the average MSE over the entire data set. SI reduces the MSE 46%, 42%, and 47% on average for low, moderate, and severe turbulence cases respectively, using 15 input frames, under daytime illumination. Similarly, the MFBD method provides 40%, 25%, and 36% improvements in MSE on average under the same conditions. Under lowlight conditions (fewer than 30 photons per pixel), improvements of 39%, 29%, and 27% are available using SI with 25 input frames used in the reconstruction, and 38, 34 and 33, respectively, for the MFBD method using 150 input frames.
The remainder of this article is organized as follows. In Sec. 2, we discuss horizontal imaging as it applies to the reconstruction of turbulent images under anisoplanatic conditions and develop foundations of the bispectrum and the objective functions to be applied to the optimization routine. In Sec. 3, the object recovery methods used for SI and MFBD are described for both Gaussian and Poisson noise models. In Sec. 4, the comparative results of processing for SI and MFBD are presented under all three turbulence cases and, in the case of the Poisson noise model, for three turbulence conditions and three values of mean photon counts. Section 5 presents the results of MFBD performance using field data followed by the results of a processing strategy that uses the bispectrum estimate to jumpstart the MFBD estimator in Sec. 6. Summary conclusions are offered in Sec. 7.

Horizontal Imaging
We begin our discussion by examining some assumptions regarding the horizontal-imaging problem. The atmospheric coherence radius, r 0 , is commonly used to define the effective aperture radius of an imaging system. Increasing the size of the aperture beyond r 0 will not improve the resolution of the system, when imaging in the presence of anisoplanatism. For spherical wave propagation, the atmospheric coherence radius, r 0 , is defined by Frieds 13 and Roggemann and Welsh 2 as (1) where λ is the mean wavelength, the wave number k ¼ 2π∕λ.
The limits of the integral are from z ¼ 0 at the pupil plane to the scene at z ¼ L. For purposes of this study, it is assumed that in Eq. (1) that the refractive-index structure constant C 2 n ðzÞ does not vary with location, for paths that are horizontal or at shallow angles from the aperture. 14 This assumption simplifies the expression of the atmospheric coherence radius given in Eq. (1) to It is further assumed that the strength of the turbulence is such that scintillation effects can be neglected. Evaluating these expressions for the simulated propagation path L ¼ 1 km and index of refraction structure C 2 n for the lowturbulence condition C 2 n ¼ 2.25 × 10 −14 m ð−2∕3Þ , and r 0 ¼ 3.33 cm. Simulation conditions for the moderate condition C 2 n ¼ 3.75 × 10 −14 m ð−2∕3Þ yields r 0 ¼ 2.45 cm. For the severe turbulence condition, C 2 n ¼ 5.25 × 10 −14 m ð−2∕3Þ and r 0 ¼ 2.01 cm. A complete description of the simulator used to create these data sets is available in Ref. 15.
Using Fried's 16 definition for the isoplanatic angle, θ o , and evaluating for the conditions under consideration, we find the isoplanatic angle to be between 10.8 and 6.63 μrad.
For the conditions simulated here, a single pixel in the simulated imaging system captures 2.79 μrad. Expressing the θ 0 values for the low, medium, and severe turbulence conditions of the simulation, we see that the isoplanatic patch covers 4, 3, and 2 pixels in the simulated imaging system and that anisoplanatic effects are present in the simulated images. These results are summarized in Table 1, and samples of the simulated images are shown in Fig. 1.

Speckle Imaging
Atmospheric turbulence can also be described by a correlation time constant τ c . The phase perturbation at every multiple of τ c can be considered independent of the previous and next interval. One strategy to mitigate the effects of atmospheric turbulence is the use of short exposure images, where the detector is exposed to light for a period less than τ c to freeze the atmospheric turbulence in a single, independent realization of the phase perturbation. The term "speckle imaging" arises from astronomical imaging, where short exposure images of stellar objects exhibit high-frequency intensity fluctuations, or "speckles." Noise-free image formation in the spatial domain is the convolution of the atmospheric PSF and the object intensity distribution, wherex is a two-dimensional (2-D) coordinate system in the image plane, andx o is a 2-D coordinate system in the object plane. A set of images that have been corrupted by additive Gaussian noise can be described.
where n k ðxÞ represents an additive noise term characterized by an independent, identically distributed Gaussian random variable with zero mean and variance σ 2 . Taking the Fourier transform of both sides of Eq. (4) allows a multiplication in the Fourier spatial-frequency domain. The product of the complex object spectrum OðfÞ with the atmospheric optical transfer function (OTF) H k ðfÞ, results in the received complex image spectrum D k ðfÞ.
Two independent steps are required to recover the amplitude and phase information associated with the received complex object spectrum. First, inverse pseudo-Wiener filtering the ensemble average power spectral density (PSD) yields an estimate of the power spectrum of the object in the scene.
where H LE ðfÞ represents the long exposure OTF of the imaging system and turbulent atmosphere and h•i represents an ensemble average.
For our processing, the theoretical long-exposure atmospheric OTF of the simulated imaging system is applied to estimate the amplitude of the complex object spectrum. 13,17 The additive term α is a constant adjusted by the user that compensates for the fact that at high signal-to-noise ratios the PSD of the OTF may go to zero.
The amplitude of the PSD is insufficient to form an image. A separate process is required to recover the phase information in the speckle images. The object phase information is recovered by computing the bispectrum, 7 defined as where D Ã ðfÞ represents the complex conjugate of DðfÞ, and Δf represents a constant spatial frequency offset. It can be shown 2 that the object phase information can be recovered from the bispectrum by substituting the image spectrum DðfÞ defined in Eq. (6) into the the definition  in Eq. (9). The estimate of the object phase can be determined at any arbitrary spatial frequencyf 0 ¼f þ Δf in terms of the bispectrum estimate, assuming the phases at f and Δf are known.
Because images are real valued, the phase off ¼ 0 is zero. The adjacent frequency points ð0; ΔfÞ, ð0; −ΔfÞ, ðΔf; 0Þ, ð−Δf; 0Þ can also arbitrarily be set to zero without loss of image quality. This assumption will result in loss of image registration information, and reconstructed images are often not centered within the image frame. Further processing can recenter the reconstructed images.

MFBD for the Gaussian Noise Model
MFBD is another common image reconstruction method 5 that uses an iterative process to recover object amplitude and phase information. The term iðxÞ is the noise-free image characterized in Eq. (4) as the convolution of a deterministic object and the imaging system's incoherent PSF. The incoherent PSF hðxÞ is the modulus squared of the coherent PSF.
where the coherent PSF is given by where A is a binary aperture function whose value is zero outside the aperture,x p is a 2-D coordinate system in the aperture plane, the • operator represents the inner product, and The phase term θðx p Þ represents the accumulated phase perturbation in the aperture of the imaging system decomposed to a linear sum of orthonormal basis functions, ϕ k ðx p Þ and α k are the weighting coefficients for each term and image. Zernike polynomials are a convenient set of basis functions commonly used in describing classical aberrations in optical systems.
Using a Gaussian noise model, each point x in each image d k ðxÞ is a random variable with a Gaussian probability density function (PDF). The PDF of d k ðxÞ is parameterized by the object intensities oðxÞ and the vector of aberration weighting coefficientsα k . 4 14) and the likelihood of the complete data set consisting of all the pixel intensities in all the corrupted images is given by Taking the natural log of Eq. (15) makes the analysis more tractable. The resulting summation, neglecting a constant term, yields the log-likelihood function The processing objective is to maximize this cost function, parameterized in terms of the true object o x and the vector of Zernike coefficients,α k that expand the phase terms θðx p Þ. This is a complex nonconvex optimization problem that may not yield a global minimum. Further discussion is provided in Sec. 3.1. To maximize this cost function, the optimization routine needs at least an estimate of the gradient, ∇½iðx;α k Þ, and Hessian, ∇ 2 ½iðx;α k Þ, of the likelihood function. Although an analytic form of the Hessian is not required, optimization is much more efficient if an analytic form of the gradient is provided. The expression of the gradient is derived in Appendix A.

Poisson Noise Model MFBD
Not all images are taken in full daylight. At low light levels, photon noise may dominate image frames. This is often characterized by a grainy quality of the images due to the random arrival times of discrete photons. Photon noise in images is described by modeling the number of photons detected in an image frame at each pixel as a Poisson random variable with a mean photon count rate λ, which is proportional to average pixel intensity. For this simulation, the number of photons detected at each detector pixel is assumed to be an independent, Poisson distributed random variable with a mean rate given by a noiseless diffraction-limited image iðxÞ. The random nature of the PSF is neglected. The probability of detecting d k ðxÞ photons at a specific pixel location is given by The distribution over the entire set of pixel locations d k is given by 4 As before, taking the natural log yields a modified loglikelihood function L poisson ½iðx;αÞ ¼ − where the last term is a constant and can be neglected. The expression for the gradient of the Poisson noise model log-likelihood function is derived in Appendix B.

Methods
Our simulations assume that the propagation occurs over horizontally homogeneous conditions with both the object and the imaging system immersed in a turbulent atmosphere. Furthermore, we assume that the height above ground does not vary significantly between imaging system and a static scene. In addition, C 2 n is assumed to be constant over the propagation path. 14 We assume that the simulated data has effectively frozen the turbulence at the moment the turbulent image is created. Reconstruction processing begins by selecting images from the complete data set in groups of N f ¼ 2, Prior to applying the simulated turbulent images to the reconstruction algorithms, they must be recentered, as tilt is not estimated in the MFBD algorithm. This was accomplished by using a correlation filter to compare each image in the stack to an ensemble average image and then shifting the turbulent image to recenter it. When the optimization routine exits, the recovered object is compared with the diffractionlimited image and the MSE is averaged over all pixels after both images are normalized to a maximum value of 255, as follows where fðxÞ is the normalized diffraction-limited image,fðxÞ is the current-normalized estimate of the object, and N 2 is the total number of pixels in the image. The image stack is incremented to begin processing the next group of K turbulent images and the process repeated. When the entire data set has been processed, the average of the vector of MSEs for images processed N f at a time is calculated. To facilitate direct comparison between the reconstruction methods, exactly the same simulated turbulent images were applied to each estimator in every trial for both methods, and the MSE was calculated from the recovered images.

Multiframe Blind Deconvolution
In order to reduce edge effects in the reconstructed object, each image in the data set was preprocessed to pad the recentered image by replicating the edges of the image outward and then adding a border of zeros. The images as applied to the estimator were padded with 8 replicated pixels followed by 5 zero pixels at the margins of each image, bringing the total size of the image to 256 × 256 pixels.
We jointly process all images and all Zernike coefficients, thus for a data set of K, N × N images, using J Zernike polynomial terms, there will be N 2 þ J × K parameters that must be jointly estimated. To make the optimization tractable and reduce run times, the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) 15 method is used to optimize the cost functions in Eqs. (16) and (19) to find the object and aberration coefficients most likely to have produced the simulated images in the data set. One of the characteristics of line-search methods like L-BFGS is the need for an estimate of the Hessian ∇ 2 fðxÞ. The limited memory form of the BFGS does not require an explicit expression for ∇ 2 fðxÞ. It estimates the value of the Hessian matrix by maintaining the last few updates of fðxÞ and ∇fðxÞ. L-BFGS is a quasi-Newtonian, "hill-climbing" technique that begins with an initial guess at a solution forx 0 and then proceeds along a line in the direction pointed to by the gradient of the objective function evaluated at each location. For MFBD processing, the spatial average of all the frames used in the trial was used as starting point for the estimator. In Sec. 6, we propose a hybrid bispectrum-MFBD (B-MFBD) method, in which a bispectrum object estimate is used to provide the starting point for the MFBD estimator. The initial guess Zernike coefficients applied to the estimator for both techniques are random Gaussian numbers with a mean of 0.5 and unity variance. Based on previous work, 18 for the image reconstruction processing in this article, the number of iterations was limited to 25 and the number of Zernike terms was fixed at 30. Figure 2 shows the plots of MSE for the Gaussian case MFBD, and bispectrum estimators, as the number of images, N f are varied for low, moderate, and severe turbulence conditions. For comparison, the average MSE for the each data set is also provided. We see that at N f ¼ 2 and thereafter both estimators can be expected to perform better than the average MSE for all three simulated turbulence conditions. For the low turbulence case, C 2 N ¼ 2.25 × 10 −14 m ð−2∕3Þ , marginal improvement in MSE declines at N f ¼ 12, reaching a maximum improvement of ∼40% over the average MSE over the entire data set. For the moderate turbulence case, C 2 N ¼ 3.75 × 10 −14 m ð−2∕3Þ , the improvement in MSE available by including additional input frames hits a maximum of ∼25% of full scale at N f ¼ 14. For the severe turbulence case, C 2 N ¼ 5.25 × 10 −14 m ð−2∕3Þ , the improvement in MSE available by including additional input frames hits a maximum of ∼36% of full scale at N f ¼ 14 and neither the MSE nor the standard deviation improves significantly as additional input images are added to the processing stack. However, for low and moderate turbulence cases, if processing time is not of consequence, the MSE and its standard deviation continues to improve as additional images are included or if the number of iterations is increased. 18 A sample of the simulated turbulent image data set is compared with a sample-reconstructed object for MFBD and SI as shown in Figs. 3-5.

Poisson Noise Model, Photon Rate 2 × 10 6 , 30 Photons Per Pixel
Each set of 1000 turbulent images, representing the the three turbulence cases, was used to generate a set of speckle images with a mean photon count per image of 2 × 10 6 . Each set of images was processed using the MFBD and SI methods described above. The MFBD method used the cost function and gradients described in Eqs. (19), (31), and (32). Figure 6 shows the plots of MSE for the Gaussian MFBD and SI reconstructors, as the number of images, N f are varied for (a) Case 1: low turbulence provided. For the low turbulence case, MFBD performance is less than the data set average until 50 input frames are used in each reconstruction. At N f ¼ 50 and thereafter the estimator can be expected to produce an estimate that has a value lower than 2095, the average MSE of the images in the simulated data set. Marginal improvement in MSE continues as additional frames are added to the image stack, reaching a maximum of about 38% improvement over the average MSE across the data set. Bispectrum achieves similar results with far fewer input frames. For the low   turbulence case, MFBD performance approximates that of the bispectrum technique after 125 images are processed in each image stack. For moderate and severe turbulence cases, the use of 150 input frames produces a slightly lower MSE than the SI method. The reduced number of trials for large frame groups is conjectured to be the reason why the MSE plots are not consistently monotonic. Processing times for SI can be a short as 1 min, MFBD can be as short as 30 min when run to convergence. A sample simulated turbulence image is compared with MFBD and SI reconstructed objects in Figs. 7-9. The MFBD sample was processed using 175 frames; the SI method used 200 frames in the reconstruction. For comparison, the average MSE for each turbulence condition is also provided. For the low-turbulence case, MFBD performance is less than the input images until 50 input frames are used in each reconstruction. The MFBD estimator produces a lower MSE than that of the   bispectrum estimator when the number of frames processed N f > 50. MFBD error reduces as additional frames are added to the image stack, reaching a maximum of about 38% improvement over the average MSE across the data set.
A sample simulated turbulent image is compared with reconstructed objects in Figs. 11-13. The MFBD sample was processed using 175 frames, the bispectrum method used 200 frames in the reconstruction. Processing times for bispectrum are significantly shorter than for MFBD.

Poisson Noise Model, Mean Photon
Rate 5 × 10 15 , 8 Photons Per Pixel Figure 14 shows the plots of MSE for the Gaussian MFBD and bispectrum estimators, as the number of images N f are varied. For comparison, the average MSE for the entire data set is also provided. On average, MFBD performance is less than the input images until 25 input frames are used in each reconstruction. At N f ¼ 25 and thereafter, the estimator can be expected to produce an estimate that has a lower MSE than the average MSE of the images in the simulated data   set. One possible reason that the SI estimate performance has diminished is that the estimator was optimized for high light levels. Only five bispectrum subplanes were used to reconstruct the phase estimate, which is sufficient for high-light levels, but the inclusion of higher-frequency offsets may increase the noise in the final estimate, thereby degrading the MSE performance. The MSE produced by the MFBD estimator will on average be less than that of SI techniques when the number of frames processed N f > 100. After N f > 125, MSE is not further reduced by MFBD processing. For the moderate and severe turbulence cases shown in Figs. 14(b) and 14(c), the MFBD estimator produces a lower MSE then SI after N f > 100. Sample simulated turbulent images are compared with reconstructed objects in Figs. 15-17. The MFBD sample was processed using 175 frames; the bispectrum method used 200 frames in the reconstruction.

MFBD Estimator Performance Using Field Data
MFBD estimator performance is also evaluated with field images that were gathered as an ancillary part of a freespace-laser communications experiment. 19 The transmitter   side of the experiment consisted of an 808-nm laser transmitter and a pinwheel target. The receiver end of the experiment was a 300-cm Celestron telescope located 3.046-km downrange. The imaging path extended horizontally over both land and water at ∼250 m above sea level. As part of the communications experiment, simultaneous images and turbulence measurements were taken using a Shack-Hartmann wave front sensor and a point gray CCD camera. The images used in MFBD reconstruction were recorded in the absence of continuous wave front sensor (WFS) measurements because the active laser present in the images recorded for the communications experiment regularly saturated the image sensor. Turbulence measurements were taken for three 27-s intervals pre-and postimage collection. WFS sensor data was not concurrent with the images. Our assumption is that the atmospheric turbulence will not change more rapidly than the data could be recorded. The Fried parameter estimates, r 0 , are shown in Table 2.
Images from the July 1, 2009, data set were processed using the MFBD estimator. A 250 × 250 pixel region of   interest was extracted from the recorded images and padded up to 256 × 256 by replicating the pixels at the margin of the image outward. Fifteen images were processed using a spatially averaged initial estimate. The images were applied to the estimator for 20 iterations, and the result is shown in Figs. 18 and 19.
Images from July 24, 2009, 1900 EDT data set were processed using the MFBD estimator. A 250 × 250 pixel region of interest was extracted from the recorded images and padded up to 256 × 256 by replicating the pixels at the margin of the image outward. R 0 for this data set ranged from 2.3 to 2.4 cm, substantially smaller than data set 1. To provide more information for the estimator to work with, 25 images were processed using 35 Zernike coefficients and 20 iterations of the estimator. A spatially averaged image was used as a starting point for the optimization routine, and the result is shown in Fig. 20.
6 Using Bispectrum as the MFBD Initial Estimate Both bispectrum and MFBD have strengths to recommend them. Bispectrum produces rapid results but requires an estimate of C 2 n and α to get started. MFBD converges more slowly and has the potential to produce a higher-quality image reconstruction. But, MFBD requires an initial estimate     of the object to get started and is subject to trapping in a local minimum of the cost function. The reduced run time of the bispectrum estimator suggests that it may offer a quick way to provide a better starting point for the MFBD estimator. Reconstruction processing in Sec. 3.1 used the spatial average of the input images as a starting point for the L-BFGS optimization routine. A hybrid strategy was attempted to combine the best of both techniques. In the first case, a bispectrum estimate for the object intensities was used to replace the ensemble average of the input images as the starting point for the MFBD estimator. The MSE was calculated as a function of the number of input frames used in the estimate. Setting the number of processing iterations at N p ¼ 40 and limiting the number of input frames to N f ¼ 15 the hybrid method, B-MFBD, reduced the MSE by 16.9%, 18.6%, and 14.9% over MFBD processing using the ensemble average as the starting point when processing the low, moderate, and severe turbulence data sets, respectively. When compared with the SI processing, reductions of 9.4%, 8.4%, and 6.6% in MSE were observed. Reconstruction times for each frame averaged ∼1 min for bispectrum, 10 min for MFBD, and 35 min for the B-MFBD hybrid. These results are summarized in Table 3. The hybrid B-MFBD improves on the performance of the standard MFBD estimator and approximates the performance of the SI estimator. The significant characteristic of the B-MFBD estimator is a 4% to 6% improvement on the variance of the estimate over SI and MFBD applied by themselves. The plots for MSE versus the number of frames are shown in Fig. 21. Gaussian (a) Case 1: low turbulence    unpredictable, this technique is unlikely to see application in an embedded system.

Conclusions
The performances of an SI estimator and an unconstrained optimization-based MFBD estimator were compared in terms of the MSE between the object reconstructed using the method and a diffraction-limited image. Three 1000image data sets of a single image distorted by low, moderate, and severe turbulence that includes anisoplanatic effects were applied to both methods and their MSE performance was evaluated. At normal illumination levels, a wide range of turbulence cases would be well served by either SI or MFBD. Point performance estimates, using a data set of 1000 simulated turbulence-corrupted images, indicate that speckle-imaging techniques reduce the MSE 46%, 42%, and 47% on average for low, moderate, and severe turbulence conditions, respectively, using 15 input frames under daytime conditions and moderate frame rates. Similarly, the MFBD method provides 40, 29, and 36 improvements in MSE on average under the same conditions. But the speckle techniques are significantly faster. The comparison is repeated under decreasing illumination conditions (less than 100 photons per pixel) where improvements of 39, 29, and 27 are available using SI methods and 25 input frames and 38, 34, and 33, respectively, for the MFBD method and 150 input frames, respectively, under the assumption that the phase errors can be characterized as a Gaussian distribution. For all simulated turbulence cases, significant reductions were observed with as few as two input images. For the Poisson case, significant results were achieved with as few as 50 frames for MFBD and 2 frames for speckle. Processing with field data with the MFBD produced qualitatively similar results over a modest range of atmospheric turbulence. The B-MFBD estimator using the SI estimator to provide the MFBD processing algorithm with an initial starting point produces modest MSE performance improvement over SI alone using 15 input frames, 15 Zernike terms and 40 iterations of the MFBD algotithm. The B-MFBD method produced more substantial improvement when compared with MFBD processing that used a spatial average of input frames as the optimizer's starting point when using similar inputs and exit criteria. Qualitative improvements were achieved by letting the optimization run until convergence, but this technique may not be suitable for embedded system implementation because of long, unpredictable processing time.