1.

Introduction

The goal of this article is to use a parameterized, multiframe blind deconvolution (MFBD) technique to reconstruct an object estimate from a set of simulated anisoplanatic images, and examine the mean square error (MSE) performance of the estimator as the parameters are varied. This article examines estimator performance as the number of frames used in the estimation is varied, and as the number of Zernike polynomial coefficients used to characterize the phase term of the point spread function are varied under the assumption of both Gaussian and Poisson noise distributions.

Every optical system using light that has propagated any appreciable distance through the atmosphere will suffer, to some degree, from turbulence induced phase aberrations. In addition to phase errors at the aperture, light propagating over longer distances or through stronger turbulence, will cause images to suffer from anisoplanatic, and possibly scintillation effects as well. Often the image blur induced by these phase aberrations is the limiting factor in the ability to recognize details of objects in the scene.

Under isoplanatic conditions, the light coming from all points in the scene can be assumed to experience similar turbulence induced aberrations. The isoplanatic angle ${\theta }_0$ is the angular separation between point sources for which the phase changes at the aperture can be considered to be significantly decorrelated. However, in many near-surface surveillance imaging scenarios, it is reasonable to assume that the field of view of the imaging system will subtend an angle wide enough that this assumption will not be valid. In this case, we describe the viewing as anisoplanatic. The longer the optical path length and the stronger the turbulence, the more severe these aberrations become, and the isoplanatic angle decreases. Increasing the size of the aperture will not improve the quality of the image under anisoplanatic conditions. Unless the seeing conditions are very favorable, anisoplanatism will play a role in most practical horizontal imaging situations. Some technique for reducing the effects of anisoplanatism is desired.

A variety of techniques have been devised to correct the turbulence induced phase aberrations in vertical imaging applications. Adaptive optic strategies using wave front sensors to control deformable mirrors have been used in celestial observation systems for many years.¹ Techniques that exploit deliberately imposed, known phase diversity^2–4 have also been used with some success. Paxman and Schulz explored this problem by creating phase diversity across multiple speckle images. This technique uses two simultaneous measurements—an in-focus image and another with a known degree of defocus applied before the other measurement is taken.⁵ This technique is limited to fields of view that do not appreciably exceed the isoplanatic angle existing at the moment the image was captured, require substantial hardware, and divide photons between two detectors. Post-detection processing of wide field of view images captured with short-exposure times is another alternative. Fraser et al. described a technique for point-by-point registration of anisoplanatic speckle images to reduce motion blur and prepare the images for other deconvolution strategies.⁶ Ayers and Dainty 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 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 the most likely object estimate is the observed turbulent image.⁸ Hybrid hardware–software strategies offer the potential to produce on-the-fly estimates of scenes, but require substantial investment in both hardware and software to produce results.⁹ Bos and Roggemann¹⁰ have reported the use of software reconstruction techniques using the bi-spectrum method in nearly real time. The use of these strategies to surveillance imaging is largely unexplored.

This article describes a method of jointly estimating object intensities and imaging system PSFs from simulated anisoplanatic images that have been corrupted by atmospheric turbulence. The image model that forms the foundation of this estimator is that of a linear shift invariant PSF and a deterministic object. It is conjectured that anisoplanatic effects of the turbulent atmosphere are compensated for by the estimator by reconstructing a spatialy averaged PSF. Bos’s¹¹ work using cross-spectrum and bi-spectrum phase reconstructions points to this potential solution. Carrano¹² has also published a work in this area that neglects the anisoplanatic effects. This investigation will be the subject of another article. The method developed here is applied to three sets of images with varying levels of turbulence, and the effectiveness is assessed by calculating the MSE between the resulting recovered object and the diffraction limited image.

We find that the MFBD reconstructed objects show significant improvement in MSE compared to the average MSE between all the images in a data set and the associated diffraction limited image. The improvement in MSE was 40% for the low turbulence case, 25% for moderate turbulence, and 36% for severe turbulence case. We also provide an estimate of the optimum number of images and Zernike coefficients to use in the future work with MFBD reconstructions.

The remainder of this article is organized as follows. In Sec. 2, we discuss the horizontal imaging problem and briefly describe the simulation that produced the data sets used in the study. In Sec. 3, the object recovery methods for the Gaussian case is described followed by the Poisson case. In Sec. 4, the results of both the Gaussian and Poisson case reconstructions are presented. Finally, some thoughts on processing and conclusions regarding the technique are provided in Sec. 5.

2.

Background

We now describe the MFBD algorithm for the Gaussian and Poisson noise models. In MFBD, the input is a set of measured noisy and turbulence corrupted images. In a stack of $K$ turbulence corrupted, but measurement noise-free images, the $k$’th image can be described as the convolution of an unchanging object in space convolved with the PSF of the optical system $s(\overrightarrow{x})$. Mathematically, this can be expressed as¹³

The PSF is the modulus squared of the coherent impulse response ${|h_k(\overrightarrow{x})|}^2$, which is in turn the inverse Fourier transform of the generalized pupil function. Mathematically, these relationships are given by

In nonblind deconvolution problems, the data collected, $d_k(\overrightarrow{x})$, is used with a known PSF $s_k(\overrightarrow{x})$ to determine $f(\overrightarrow{x})$. In blind deconvolution, we are given $d_k(\overrightarrow{x})$ and use that information to estimate both the object, $f(\overrightarrow{x})$ and the PSF $s_k(\overrightarrow{x})$ jointly. There is no closed form solution to the problem of jointly estimating an object and the aberration parameters for each image frame. Hence, an iterative approach is needed to find the object pixel intensities and Zernike coefficients that are most likely to have resulted in the simulated data for each image. In the next section, we describe two such approaches, one based on a Gaussian noise model and other based on a Poisson noise model.

2.1.

Data Set

It is common to simulate the effects of the turbulent atmosphere by placing layers of uniform turbulence between the object and the imaging system. The data set consisting of 1000 simulated turbulent images used in this article was created using an image common in the literature. Five Kolmogorov phase screens were generated. The image was propagated over a distance of 1000 m. Light from each object pixel was projected through the phase screens, in turn, at 200-m separations using a geometric optics approach to account for the effects of anisoplanatism. Phase errors accumulating from each screen are combined at the pupil to create a turbulence degraded PSF. Each of the PSFs is then scaled by the object pixel intensities to create a turbulence corrupted image for low, moderate, and severe turbulence conditions. Parameters for the simulated imaging system include a 10-cm aperture with a $358\times 358\text{pixel}$ detector and a 0.7-mm pixel pitch. A fuller description of the simulator used to create this data set is available in Ref. 11.

For the conditions simulated here, a single pixel in the simulated imaging system captures 2.79 μrad. Expressing the ${\theta }_0$ values for the low, medium, and severe turbulence conditions of the simulation, we see in Table 1 that the isoplanatic patch covers 4, 3, and 2 pixels in the simulated imaging system (Fig. 1).

Table 1

Atmospheric simulation turbulence conditions.

Atmospheric turbulence parameters
Severity condition	$C_N^2$ (${\mathrm{m}}^{(-2/3)}$)	$r_0$ spherical case (cm)	${\theta }_0$ (μrad)	${\theta }_0$ (pixels)
Low	$2.25\times {10}^{-14}$	3.33	10.8	4
Moderate	$3.75\times {10}^{-14}$	2.45	7.75	3
High	$5.25\times {10}^{-14}$	2.01	6.63	2

Fig. 1

Horizontal imaging simulator output. Single image representative samples of the horizontal image simulator output: (a) diffraction-limited image, (b) image for $C_N^2=2.25\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$, (c) image for $3.75\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$, and (d) image for $5.25\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$.

3.

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 and $C_n^2$ is a constant over the propagation path.¹⁵ We assume that the simulated data has effectively frozen the turbulence at the moment the turbulent image is created. Prior to applying the simulated turbulent images to the reconstruction algorithm, 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 the ensemble average image and then shifting the turbulent image to recenter it. In order to reduce the aliasing associated with the finite support, each frame of the data set was preprocessed to pad the centered image by replicating the edges of the image outward and then adding a border of 0. The abrupt transitions artificially introduced by the padding process can result in high spatial frequency components that are sometimes mitigated by the application of spatial filters. Using 15 frames in the reconstruction, the image stack was padded and then a Tukey¹⁶ tapered filtered was applied to the image. Both the tapered and untapered images were applied to the estimator. The elapsed processing time and MSE of the reconstructed object were determined with the estimator limited to 20 iterations. The amount of padding for subsequent processing was determined by examining the effect on the processing time and the MSE as the amount of padding was varied. All subsequent processing was accomplished by padding each recentered turbulent image but without tapering it. The images are applied to the estimator with eight replicated pixels followed by five 0 pixels at the margins of each image, bringing the total size of the image to $256\times 256\text{pixels}$ as seen in Fig. 2. These results are summarized in Table 2.

Fig. 2

Padding example: unpadded image (a), and padded image (b).

Table 2

Selection of image padding.

3.1.

L-BFGS Optimization

The cost functions in Eqs. (10) and (15) are parameterized by the object pixel intensities and aberration coefficients, and are applied to a nonlinear optimization MATLAB routine to find the object and aberration coefficients most likely to have produced the images that were simulated in the data set. The intensities at each pixel location in each image are vectorized. The vectorized initial guesses for each image’s Zernike polynomial coefficients are appended to the end of the vector of image intensities formatted as shown in Table 3. We are jointly processing all images and all Zernike coefficients, thus for a data set of $K$, $N\times N$ images, using $J$ Zernike polynomial terms, there will be $N^2+J\times K$ parameters that must be jointly estimated. The optimization routine will return a vector of the reconstructed object’s intensities followed by the estimate of the Zernike coefficients for each frame of the input stack as shown in Table 4. Optimization over such a large parameter space is impractical using conventional optimization techniques. To make the optimization tractable, we use the L-BFGS method to process the images. L-BFGS is a quasi-Newtonian “hill-climbing” technique that begins with an initial guess at a solution for ${\overrightarrow{x}}_0$ and then proceeds along a line in the direction pointed to by the gradient of the objective function evaluated at each pixel location. One of the drawbacks to searching along the gradient is the need for the Hessian ${\nabla }^{2}f(\overrightarrow{x})$ to prevent the estimate from hopping back and forth across the sides of the valley. The limited memory form of the BFGS does not require an explicit expression for ${\nabla }^2$. It estimates the value of the Hessian matrix by maintaining the last few updates of $f(\overrightarrow{x})$ and $\nabla f(\overrightarrow{x})$. A Quasi-Newtonian line search optimization can quickly converge to a local minimum for cost functions, but there is no guarantee that the minimum is a global minimum. In processing the initial object estimate, the average of all the frames used in the trial was applied to the estimator.

Table 3

Input object and Zernike coefficients. The current stack of images being processed is spatially averaged and stripped into a vector with the intensities at the beginning and the initial guess at each image’s Zernike coeffiecients at the end.

Estimator input parameters
${\overline{o}}_1$	${\overline{o}}_1$	…	${\overline{o}}_{N^2}$	${\overline{\alpha }}_{1,1}$	${\widehat{\alpha }}_{1,2}$	…	${\widehat{\alpha }}_{1,J}$	…	${\widehat{\alpha }}_{K,1}$	${\widehat{\alpha }}_{K,2}$	…	${\widehat{\alpha }}_{K,J}$

Table 4

Vectorized reconstructed object and Zernike coefficients. The estimate is returned as a vector with the estimated object pixel intensities at the beginning and the estimate of each input image’s Zernike coeffiecients at the end.

Output object estimate and Zernike coefficient vector
${\widehat{o}}_1$	${\widehat{o}}_1$	…	${\widehat{o}}_{N^2}$	${\widehat{\alpha }}_{1,1}$	${\widehat{\alpha }}_{1,2}$	…	${\widehat{\alpha }}_{1,J}$	…	${\widehat{\alpha }}_{K,1}$	${\widehat{\alpha }}_{K,2}$	…	${\widehat{\alpha }}_{K,J}$

It is necessary to provide the estimator with a stopping criterion. Using the low condition data set with the number of frames set to 5, 15, and 35, the log-likelihood function value was monitored during reconstruction. The results are shown in Fig. 3. Regardless of the number of frames used in the reconstruction, assuming one call to the likelihood function per iteration, the reconstruction is essentially complete after 10 iterations. For the image reconstruction processing in this article, the number of iterations was limited to 25.

Fig. 3

Gaussian MSE versus number of function calls.

3.3.

Number of Zernike Terms Needed in the Optimization Process

Recovering a common object estimate from the stack of degraded images is computationally intense regardless of the method used. Using more Zernike polynomial terms require more variables to be estimated, and longer processing will result. Figure 4(a) shows the processing time required for a fixed number of frames when the number of Zernike polynomial terms is varied. Figure 4(b) shows how the processing time varies for a fixed number of Zernike coefficients as the number of input images is varied. Of greater impact on processing time is the number of images used to recover the object. For a set of $J$, $N\times N$ images, the number of variables increases as $J\times N^2$. Previous work indicated that 15 images and 35 Zernike terms would provide a good estimate of the object.¹⁷ Further exploration over a larger data set yielded similar results and additional insight into the estimator’s performance. With the number of images set to 50 in order to reduce the influence of the number of images on the outcome, the number of Zernike coefficients was varied from 10 to 100 terms for all three turbulence conditions as shown in Fig. 5. For all three turbulence cases of $C_n^2$, additional terms beyond 60 do not significantly improve the MSE.

Fig. 4

Processing time versus number of Zernike terms (a), and number of frames used in reconstruction (b).

Fig. 5

Mean square error (MSE) versus number of Zernike coefficients. With the number of frames held constant at 50, the number of Zernike coefficients are varied from 10 to 100 for (a) low condition $C_N^2=2.25\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$, (b) moderate condition $C_N^2=3.75\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$, and (c) severe condition $C_N^2=5.25\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$.

4.

Results

4.1.

Gaussian Noise Model

4.1.1.

Case 1 low condition ${\mathsf{C}}_{\mathsf{N}}^{\mathsf{2}}\mathsf{=}\mathsf{2.25}\mathsf{\times }{\mathsf{10}}^{\mathsf{-}\mathsf{14}}{\mathsf{m}}^{(\mathsf{-}\mathsf{2}/\mathsf{3})}$

Examining Fig. 6, we see that at $N=2$ and thereafter the estimator can be expected to perform better than the average MSE for the simulated image which was 673. Marginal improvement in MSE declines at $N=12$, reaching a maximum improvement of approximately 40% over the average MSE over the entire data set. However, if processing time is not of consequence, the MSE and its standard deviation continues to improve as additional images are added. Examining the results of the Zernike term evaluation shown in Fig. 5, 60 Zernike coefficients were used to characterize the PSFs and the results are compared to the estimator’s performance using 30 Zernike terms in Fig. 6. The use of additional Zernike terms does not add as much processing time as using additional frames, but each reconstruction will take longer as additional Zernike terms are used to characterize the PSFs. The incremental improvement in MSE is not worth the additional time consumed. The diffraction limited image is compared to a sample of the simulated turbulent image data set images and a sample reconstructed object in Fig. 7.

Fig. 6

MSE versus number of frames. Case 1 low turbulence $C_N^2=2.25\times {10}^{-14}$ with 30 and 60 Zernike terms.

Fig. 7

Case 1 sample images. Compares the (a) diffraction-limited image with (b) the single sample image and (c) a sample reconstructed object.

4.1.2.

Case 2 moderate condition ${\mathsf{C}}_{\mathsf{N}}^{\mathsf{2}}\mathsf{=}\mathsf{3.75}\mathsf{\times }{\mathsf{10}}^{\mathsf{-}\mathsf{14}}{\mathsf{m}}^{(\mathsf{-}\mathsf{2}/\mathsf{3})}$

Figure 8 shows (a) the diffraction limited image, (b) a sample recorded image, and (c) a sample reconstructed object. We see in Fig. 9 that the MFBD estimator will consistently perform on average better than the average image’s error as soon as the size of the processing window reaches two frames. At $N=2$ and thereafter the estimator can be expected to perform better than the average MSE of the simulated image. The improvement in MSE available by including additional input frames hits a maximum of approximately 25% of full scale at $N=14$. Neither the MSE nor the standard deviation improves significantly as additional input images are added to the stack. As a consequence of the results of the Zernike term sweeps discussed above, the estimator was run using 60 Zernike coefficients to characterize the PSFs. The results are compared in Fig. 9. The use of additional Zernike terms does not incur as large a computational penalty as that associated with adding additional frames but each reconstruction will take longer. The incremental improvement in MSE is not worth the additional processing time.

Fig. 8

Case 2 sample images. Compares the (a) diffraction-limited image with (b) the single sample image and (c) a sample reconstructed object.

Fig. 9

MSE versus number of frames. Case 2 moderate turbulence $C_N^2=3.75\times {10}^{-14}$ with 30 and 60 Zernike terms.

4.1.3.

Case 3 severe condition ${\mathsf{C}}_{\mathsf{N}}^{\mathsf{2}}\mathsf{=}\mathsf{5.25}\mathsf{\times }{\mathsf{10}}^{\mathsf{-}\mathsf{14}}{\mathsf{m}}^{(\mathsf{-}\mathsf{2}/\mathsf{3})}$

Again we see in Fig. 10 that the estimator requires at least two input frames to reliably produce an estimate of the object that has a lower MSE than 1165, the average simulated image’s MSE. At $N=2$ and thereafter the estimator can be expected to perform better than the average MSE of the recorded image. The improvement in MSE available by including additional input frames hits a maximum of approximately 36% of full scale at $N=14$ and neither the MSE nor the standard deviation improves significantly as additional input images are added to the processing stack. Figure 11 shows (a) the diffraction limited image, (b) a sample recorded image, and (c) a sample reconstructed object. Based on the results of the Zernike term sweeps discussed above, the estimator was run using 60 Zernike coefficients to characterize the PSFs. The results are compared in Fig. 10. The use of additional Zernike terms does not incur as large a computational penalty as that associated with adding additional frames but each reconstruction will take longer. As shown, the incremental improvement in MSE is not worth the additional processing time.

Fig. 10

MSE versus number of frames. Case 3 severe turbulence $C_N^2=5.25\times {10}^{-14}$ with 30 and 60 Zernike terms.

Fig. 11

Case 3 sample images. Compares the (a) diffraction-limited image with (b) the single sample image and (c) a sample reconstructed object.

4.2.

Poisson Noise Model Mean Photon Rate $2\times {10}^6$

Each set of 1000 turbulent images representing the three turbulence cases was used to generate a set of speckle images with a mean photon count per image of $2\times {10}^6$. Each set of images was processed using the MFBD methods described above using the cost function and gradient described in Eqs. (15)–(17).

4.2.1.

Case 1 low condition ${\mathsf{C}}_{\mathsf{N}}^{\mathsf{2}}\mathsf{=}\mathsf{2.25}\mathsf{\times }{\mathsf{10}}^{\mathsf{-}\mathsf{14}}{\mathsf{m}}^{(\mathsf{-}\mathsf{2}/\mathsf{3})}$

Examining Fig. 12, we see that on average, MFBD performance is less than the input images until 50 input frames are used in each reconstruction. At $N=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. However, if processing time is not of consequence, the MSE and its standard deviation continues to improve as additional images are added to the stack of images presented to the estimator, so additional improvement in the quality of the image is available. The diffraction limited image is compared to a sample simulated turbulent image and a sample reconstructed object as shown in Fig. 13.

Fig. 12

MSE versus number of frames. Poisson Case 1 Low turbulence $C_N^2=2.25\times {10}^{-14}$, mean photon rate $2\times {10}^6$.

Fig. 13

Case 1 sample images, mean photon rate $2\times {10}^6$, $C_N^2=2.25\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$. Compares the (a) diffraction-limited image with (b) a single sample image and (c) a sample reconstructed object.

4.2.2.

Case 2 moderate condition ${\mathsf{C}}_{\mathsf{N}}^{\mathsf{2}}\mathsf{=}\mathsf{3.75}\mathsf{\times }{\mathsf{10}}^{\mathsf{-}\mathsf{14}}{\mathsf{m}}^{(\mathsf{-}\mathsf{2}/\mathsf{3})}$

Figure 14 shows (a) the diffraction limited image, (b) a sample recorded image, and (c) a sample reconstructed object using 175 frames to estimate the object. Again we see in Fig. 15 that the MFBD estimator will not perform on average any better than 2285, the average simulated turbulent image error, until the number of images processed reaches 50 frames. At $N=50$ and thereafter the estimator can be expected to perform better than the average MSE of the simulated image reaching a maximum of about 34% improvement. The marginal improvement in MSE available by including additional input frames begins to decline at about $N=175$ and neither the MSE nor the standard deviation improves significantly as additional input images are processed.

Fig. 14

Case 2 sample images, mean photon rate $2\times {10}^6$, $C_N^2=3.75\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$. Compares the (a) diffraction-limited image with (b) a single sample image and (c) a sample reconstructed object.

Fig. 15

MSE versus number of frames. Poisson case 2 moderate turbulence $C_N^2=3.75\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$, mean photon rate $2\times {10}^6$.

4.2.3.

Case 3 severe condition ${\mathsf{C}}_{\mathsf{N}}^{\mathsf{2}}\mathsf{=}\mathsf{5.25}\mathsf{\times }{\mathsf{10}}^{\mathsf{-}\mathsf{14}}{\mathsf{m}}^{(\mathsf{-}\mathsf{2}/\mathsf{3})}$

Again we see in Fig. 16 that the MFBD estimator will not perform on average any better than the average image error until the number of images offered to the estimator reaches 50 frames. At $N=50$ and thereafter the estimator can be expected to perform better than the average MSE (2456) of the simulated image reaching a maximum of about 33%. The marginal improvement in MSE available by including additional input frames begins to decline at about $N=175$ and neither the MSE nor the standard deviation seems to improve significantly from there as additional input images are added to the stack. Figure 17 shows (a) the diffraction limited image, (b) a sample recorded image, and (c) a sample reconstructed object.

Fig. 16

MSE versus number of frames. Poisson case 3 severe turbulence $C_N^2=5.25\times {10}^{-14}$, mean photon rate $2\times {10}^6$.

Fig. 17

Case 3 sample images, mean photon rate $2\times {10}^6$, $C_N^2=5.25\times {10}^{-14}{\mathrm{m}}^{(-2/3)}$. Compares the (a) diffraction-limited image with (b) a single sample image and (c) a sample reconstructed object.

5.

Conclusions

The performance of an unconstrained optimization-based MFBD estimator was evaluated in terms of the MSE between the reconstructed object and a diffraction limited image. Three 1000-image data sets of a single image distorted by low, moderate, and severe turbulence were generated using a horizontal imaging simulator that includes anisoplanatic effects. The data sets were then applied to the estimator and its MSE performance evaluated. If a hardware implementation was to be produced with a fixed, or limited set of operator options, a wide variety of turbulence cases would be well served by a selection of 14 images and 30 polynomial terms for use with the estimator. Point performance estimates, using a data set of 1000 simulated turbulence corrupted images, indicate that the algorithm is capable of producing 40%, 25%, and 36% improvements in MSE for low, moderate, and severe-anisoplanitic turbulence cases, 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, but 175 frames would be a reasonable place to design a system that would be able to cope with a variety of atmospheric turbulence and light levels. For further research, it may be possible to speed up the reconstruction by providing a better initial guess at the object. Simulated annealing techniques could be used to perturb the estimate away from a local minimum and may prove to be an effective answer to local minimum trapping.