## Multiframe blind deconvolution with real data: imagery of the Hubble Space Telescope

Optics Express, Vol. 1, Issue 11, pp. 355-362 (1997)

http://dx.doi.org/10.1364/OE.1.000355

Acrobat PDF (221 KB)

### Abstract

Multiframe blind deconvolution - the process of restoring resolution to blurred imagery when the precise form of the blurs is unknown - is discussed as an estimation-theoretic method for improving the resolving power of ground-based telescopes used for space surveillance. The imaging problem is posed in an estimation-theoretic framework whereby the object’s incoherent scattering function is estimated through the simultaneous identification and correction of the distorting effects of atmospheric turbulence. An iterative method derived via the expectation-maximization (EM) procedure is reviewed, and results obtained from telescope imagery of the Hubble Space Telescope are presented.

© Optical Society of America

## 1. Introduction

^{11. M. C. Roggemann and B. Welsh, Imaging Through Turbulence, CRC Press, Inc. (1996).}. These methods can be broadly classified into two categories: post-processing through signal processing whereby the image distortions are identified and compensated by computer processing of the recorded imagery; and pre-processing through adaptive optics whereby the distortions are identified and compensated with wavefront sensors and deformable mirrors so that undistorted imagery is, in principle, recorded directly. With current technology the cost of adaptive compensation can be prohibitive, particularly if full compensation is desired. Because of this, many systems in use today benefit only from partial compensation or none at all, and post-processing methods such as the one discussed in this paper provide an economical means for restoring a telescope’s lost resolving power.

*r*

_{0}, which provides a measure of the largest spatial frequency that is present in the long-exposure imagery. Whereas a diffraction-limited telescope records spatial frequencies as high as the ratio

*D*/

*λd*

_{o}, where

*D*is the telescope diameter,

*λ*is the nominal wavelength being sensed, and

*d*

_{o}is the distance from the telescope to the object, long-exposure imagery acquired through atmospheric turbulence will only contain frequencies to the turbulence-limited cut-off frequency

*r*

_{0}/

*λd*

_{0}, where

*r*

_{0}can take values under 10cm. For a 2m telescope with

*r*

_{0}= 10cm, then, recovery of diffraction-limited resolution requires one to super-resolve the long-exposure imagery by a factor of 20. Such improvements are rarely, if ever, accomplished with real data.

^{22. A. Labeyrie, “Attainment of diffraction limited resolution in large telescopes by Fourier analyzing speckle patterns in star images,” Astron. Astrophys. , 6, 85 (1970).}, and has since spawned a number of signal-processing methods broadly classified as

*speckle imaging techniques*. A review of these is found in Ref. 1 and the many references within. Many of these, such as Labeyrie’s speckle interferometry, the Knox-Thompson method

^{33. K. T. Knox and B. J. Thompson, “Recovery of images from atmospherically degraded short exposure images,” Astrophys. J. , 193, L45–L48 (1974). [CrossRef] }, and the triple-correlation method

^{44. A. W. Lohmann, G. Weigelt, and B. Wirnitzer, “Speckle masking in astronomy: triple correlation theory and applications,” Appl. Opt. , 22, 4028–4037 (1983). [CrossRef] [PubMed] }, rely on the existence of a data transformation that is invariant to the effects of atmospheric turbulence. After performing this transformation along with any required calibration, a non-trivial image-recovery problem must be solved to obtain an estimate of the underlying object. As an alternative to these methods, the data can be processed to directly estimate and correct for the effects of atmospheric turbulence. Methods based on this approach are broadly classified as

*multiframe blind deconvolution*

^{55. T. J. Schulz, “Multi-frame blind deconvolution of astronomical images”, J. Opt. Soc. Am. A , 10, 1064–1073 (1993). [CrossRef] }techniques. The application of a multiframe blind deconvolution method to real telescope data is the focus of this paper.

## 2. Mathematical model

*u*and

*y*are two-dimensional spatial variables in the telescope-pupil and image planes, respectively;

*A*(·) is the telescope pupil function;

*θ*

_{t}(·) is the time-varying phase distortion caused by turbulence-induced changes in the refractive index of the atmosphere;

*λ*is the nominal wavelength; and

*d*

_{o}is the distance from the object to the telescope. This is the isoplanatic approximation and results in an imaging model that is linear and spatially invariant. In the absence of other effects, then, the time-varying intensity recorded by the telescope’s camera is:

*o*(·) is the object’s incoherent scattering function, after magnification and inversion

^{66. Strictly speaking, because of image inversion and magnification an imaging system is never spatially invariant; however, if one views the input signal as the inverted and magnified object, many imaging systems are then well-modeled as spatially invariant.}. The notation

*i*(·;

*θ*

_{t},

*o*) shows the image intensity’s explicit dependence on the unknown object

*o*(·) and turbulence-induced phase aberrations

*θ*

_{t}(·).

7. D. L. Snyder, A. M. Hammoud, and R. L. White, “Image recovery from data acquired with a charge-coupled-device camera,” J. Opt. Soc. Am. A , **10**, 1014–1023 (1993). [CrossRef] [PubMed]

*n*th detector element in a CCD camera array during the

*k*th exposure is:

*N*

_{k}[

*n*] is a Poisson distributed random variable that models the number of object-induced photocounts recorded by the detector,

*M*

_{k}[

*n*] is a Poisson distributed random variable that models the number of background-induced (both internal and external) photocounts recorded by the detector,

*g*

_{k}[

*n*] is a Gaussian distributed random variable that models the CCD read-out noise, and

*b*[

*n*] is a deterministic bias induced by the electronic circuitry used to acquire and record the data. All of these variables are usually independent from each other and across detector elements.

*t*

_{k}is the time and

*T*the duration of the

*k*th exposure,

*y*

_{n}is the spatial region over which the

*n*th detector element integrates the image intensity, and

*γ*[

*n*] is a spatially-varying scale factor that accounts for detector efficiency. Pixels that record meaningless data are accommodated with this model by setting

*γ*[

*n*] = 0. The phase aberrations are usually constant over a short exposure, and the detector regions are often small in size relative to the fluctuations in the image intensity so that the integrating operation can be modeled by the sampling operation:

*y*

_{n}is the location of the

*n*th detector element and |

*y*| is its integration area. Furthermore, the object and point-spread functions are often approximated on a discrete grid so that:

_{y}, Δ

_{x}, and Δ

_{u}are the image-, object-, and pupil-plane grid sizes, respectively. When the pupil-plane array size is

*N*×

*N*and Δ

_{u}Δ

_{y}/

*λd*

_{o}= 1/

*N*, the computations required to create the point-spread function can be accomplished efficiently by using a fast Fourier transform (FFT) algorithm.

*M*

_{k}[

*n*] is denoted by

*I*

_{b}[

*n*] and assumed to be the same for each exposure. The CCD read-out noise

*g*

_{k}[

*n*] are typically modeled as zero-mean, independent random variables with variance

*σ*

^{2}[

*n*]. Appropriate values for the gain function

*γ*[·], background mean

*I*

_{b}[·], read-noise variance

*σ*

^{2}[·], and bias

*b*[·] are usually determined through a carefully controlled study of the data acquisition system. This typically involves the analysis of data acquired both without illumination (dark frames) and with uniform illumination (flat fields).

## 3. Maximum-likelihood estimation

*N*object- and background-induced photocounts in the

*n*th detector element during the

*k*th exposure is:

*o*[

*m*] =

*T*|

*y*|

*o*(

*x*

_{m}) contains the dependence on the unknown object, and

*h*(·;

*θ*) contains the dependence on the unknown phase aberrations. Furthermore, the probability density for the read-out noise is

_{tk}*d*

_{k}[·]}, the maximum-likelihood estimate of the unknown object and phase aberrations is selected to maximize the likelihood:

*o, θ*) an overly complicated problem. When the read-out noise variance is large (greater than 50 or so), however, the data can be modified according to:

^{88. D. L. Snyder, C. W. Helstrom, A. D. Lanterman, M. Faisal, and R. L. White, “Compensation for readout noise in CCD images,” J. Opt. Soc. Am. A , 12, 272–283 (1985). [CrossRef] }. Accordingly, the modified data are approximately Poisson distributed with the mean value

*I*

_{k}[

*n*] +

*I*

_{b}[

*n*] +

*σ*

^{2}[

*n*]. The probability mass function for these data is then modeled as:

*o, θ*) subject to any physical constraints that are imposed on these parameters. For all problems we have considered, the object intensity

*o*[·] is constrained to be a nonnegative function. Known object support is another constraint that is commonly used; however, implementation of this constraint in an optimization procedure usually results in a trivial reduction of the number of unknown parameters.

### 3.1 Comments

^{99. G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its application”, Opt. Lett. , 13, 547–549 (1988). [CrossRef] [PubMed] }class of iterative algorithms; and ii) constrained numerical optimization algorithms. The methodology described here is one of constrained numerical optimization, and is a straightforward extension of the method proposed by Schulz

^{55. T. J. Schulz, “Multi-frame blind deconvolution of astronomical images”, J. Opt. Soc. Am. A , 10, 1064–1073 (1993). [CrossRef] }with important modifications to accommodate real sensor artifacts such as read-noise, internal and external background, and non-uniform gain. Other methods of constrained optimization have been proposed

^{1010. R. G. Lane, “Blind deconvolution of speckle images,” J. Opt. Soc. Am. A , 9, 1508–1514 (1992). [CrossRef] ,1111. S. M. Jefferies and J. C. Christou, “Restoration of astronomical images by iterative blind deconvolution”, Astrophys. J. , 63, 862–874 (1993). [CrossRef] ,1212. E. Thiebaut and J.-M. Conan, “Strict a priori constraints for maximum-likelihood blind deconvolution,” J. Opt. Soc. Am. A , 12, 485–492 (1995). [CrossRef] }, and some have relied on essentially the same methodology as the one proposed in Ref. 5

5. T. J. Schulz, “Multi-frame blind deconvolution of astronomical images”, J. Opt. Soc. Am. A , **10**, 1064–1073 (1993). [CrossRef]

^{1212. E. Thiebaut and J.-M. Conan, “Strict a priori constraints for maximum-likelihood blind deconvolution,” J. Opt. Soc. Am. A , 12, 485–492 (1995). [CrossRef] }, and the utilization (or lack thereof) of realistic sensor models. Key points regarding the methodology presented in this paper are:

- A practical and faithful model is used for the data collection process. By using the maximum-likelihood method for estimation, this model is in turn used to induce a cost function for a constrained optimization problem.
- The importance of the point-spread function model (Eq. 7) should not be understated. This point was emphasized by Cornwell in his summary of the 1994 workshop on the Restoration of HST Images and Spectra II as he commented on the subject of blind deconvolution
^{1313. T. J. Cornwell, “Where have we been, where are we now, where are we going?”, in The Restoration of HST Images and Spectra II, B. Hanisch and R. L. White, editors, (Space Telescope Science Institute, Baltimore, MD, 1993) pp. 369–372.}:“One always gains by adding more information in the form of known imaging physics. For example, one can ask whether the closure relations are enforced? Alternatively, is the PSF forced to be the Fourier transform of a cross-correlation of a pupil plane with phase errors? If not, then it should be.”

### 3.2 Numerical optimization

5. T. J. Schulz, “Multi-frame blind deconvolution of astronomical images”, J. Opt. Soc. Am. A , **10**, 1064–1073 (1993). [CrossRef]

*algorithm*performance (computation time). The methodology used to define the problem, which is greatly dependent upon the use of sound mathematical and physical models, will determine

*estimator*performance (bias and variance).

5. T. J. Schulz, “Multi-frame blind deconvolution of astronomical images”, J. Opt. Soc. Am. A , **10**, 1064–1073 (1993). [CrossRef]

## 4. Results with real data

**10**, 1064–1073 (1993). [CrossRef]

14. T. J. Schulz, “Movie of processed Hubble Space Telescope imagery,” http://www.ee.mtu.edu/faculty/schulz/mpeg/hst.mpeg

## 5. Summary

^{1515. R. G. Paxman, T. J. Schulz, and J. R. Fienup, “Joint estimation of object and aberrations using phase diversity,” J. Opt. Soc. Am. A , 9, 1072–1085 (1992). [CrossRef] }; clever methods for generating initial parameter estimates; and the development and implementation of improved methods of optimization.

## 6. Acknowledgments

## References and links

1. | M. C. Roggemann and B. Welsh, |

2. | A. Labeyrie, “Attainment of diffraction limited resolution in large telescopes by Fourier analyzing speckle patterns in star images,” Astron. Astrophys. , |

3. | K. T. Knox and B. J. Thompson, “Recovery of images from atmospherically degraded short exposure images,” Astrophys. J. , |

4. | A. W. Lohmann, G. Weigelt, and B. Wirnitzer, “Speckle masking in astronomy: triple correlation theory and applications,” Appl. Opt. , |

5. | T. J. Schulz, “Multi-frame blind deconvolution of astronomical images”, J. Opt. Soc. Am. A , |

6. | Strictly speaking, because of image inversion and magnification an imaging system is never spatially invariant; however, if one views the input signal as the inverted and magnified object, many imaging systems are then well-modeled as spatially invariant. |

7. | D. L. Snyder, A. M. Hammoud, and R. L. White, “Image recovery from data acquired with a charge-coupled-device camera,” J. Opt. Soc. Am. A , |

8. | D. L. Snyder, C. W. Helstrom, A. D. Lanterman, M. Faisal, and R. L. White, “Compensation for readout noise in CCD images,” J. Opt. Soc. Am. A , |

9. | G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its application”, Opt. Lett. , |

10. | R. G. Lane, “Blind deconvolution of speckle images,” J. Opt. Soc. Am. A , |

11. | S. M. Jefferies and J. C. Christou, “Restoration of astronomical images by iterative blind deconvolution”, Astrophys. J. , |

12. | E. Thiebaut and J.-M. Conan, “Strict |

13. | T. J. Cornwell, “Where have we been, where are we now, where are we going?”, in |

14. | T. J. Schulz, “Movie of processed Hubble Space Telescope imagery,” http://www.ee.mtu.edu/faculty/schulz/mpeg/hst.mpeg |

15. | R. G. Paxman, T. J. Schulz, and J. R. Fienup, “Joint estimation of object and aberrations using phase diversity,” J. Opt. Soc. Am. A , |

**OCIS Codes**

(100.1830) Image processing : Deconvolution

(100.3020) Image processing : Image reconstruction-restoration

(100.3190) Image processing : Inverse problems

**ToC Category:**

Focus Issue: Signal collection and recovery

**History**

Original Manuscript: September 17, 1997

Published: November 24, 1997

**Citation**

Timothy Schulz, Bruce Stribling, and Jason Miller, "Multiframe blind deconvolution with real data:
imagery of the Hubble Space Telescope," Opt. Express **1**, 355-362 (1997)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-1-11-355

Sort: Journal | Reset

### References

- M. C. Roggemann and B. Welsh, Imaging Through Turbulence, CRC Press, Inc. (1996).
- A. Labeyrie, "Attainment of diffraction limited resolution in large telescopes by Fourier analyzing speckle patterns in star images," Astron. Astrophys., 6, 85 (1970).
- K. T. Knox and B. J. Thompson, "Recovery of images from atmospherically degraded short exposure images," Astrophys. J., 193, L45-L48 (1974). [CrossRef]
- A. W. Lohmann, G. Weigelt and B. Wirnitzer, "Speckle masking in astronomy: triple correlation theory and applications," Appl. Opt. , 22, 4028-4037 (1983). [CrossRef] [PubMed]
- T. J. Schulz, "Multi-frame blind deconvolution of astronomical images", J. Opt. Soc. Am. A , 10, 1064-1073 (1993). [CrossRef]
- Strictly speaking, because of image inversion and magnification an imaging system is never spatially invariant; however, if one views the input signal as the inverted and magnified object, many imaging systems are then well-modeled as spatially invariant.
- D. L. Snyder, A. M. Hammoud, and R. L. White, "Image recovery from data acquired with a charge-coupled-device camera," J. Opt. Soc. Am. A , 10, 1014-1023 (1993). [CrossRef] [PubMed]
- D. L. Snyder, C. W. Helstrom, A. D. Lanterman, M. Faisal, and R. L. White, "Compensation for readout noise in CCD images," J. Opt. Soc. Am. A , 12, 272-283 (1985). [CrossRef]
- G. R. Ayers and J. C. Dainty, "Iterative blind deconvolution method and its application", Opt. Lett. , 13, 547-549 (1988). [CrossRef] [PubMed]
- R. G. Lane, "Blind deconvolution of speckle images," J. Opt. Soc. Am. A , 9, 1508-1514 (1992). [CrossRef]
- S. M. Jefferies and J. C. Christou, "Restoration of astronomical images by iterative blind deconvolution", Astrophys. J., 63, 862-874 (1993). [CrossRef]
- E. Thiebaut and J.-M. Conan, "Strict a priori constraints for maximum-likelihood blind deconvolution," J. Opt. Soc. Am. A , 12, 485-492 (1995). [CrossRef]
- T. J. Cornwell, "Where have we been, where are we now, where are we going?", in The Restoration of HST Images and Spectra II, B. Hanisch and R. L. White, editors, (Space Telescope Science Institute, Baltimore, MD, 1993) pp. 369-372.
- T. J. Schulz, "Movie of processed Hubble Space Telescope imagery," http://www.ee.mtu.edu/faculty/schulz/mpeg/hst.mpeg
- R. G. Paxman, T. J. Schulz, and J. R. Fienup, "Joint estimation of object and aberrations using phase diversity," J. Opt. Soc. Am. A , 9, 1072-1085 (1992). [CrossRef]

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

OSA is a member of CrossRef.