© ESO, 2011

1. Introduction

Astronomical observation using ground-based telescopes is significantly degraded by diffraction-index fluctuations caused by atmospheric turbulence. This turbulence arises from local temperature and density inhomogeneities and results in a time- and space-variant point spread function (PSF). Often the PSF is assumed to be invariant within a short time-period and a small region of space, called an isoplanatic patch. The coherence time and the size of the isoplanatic patch depend on the strength of the turbulence that is usually quantified by Fried’s parameter r_o (Fried 1978) ranging between 10−20 cm for visible wavelengths at astronomical telescope sites. The coherence time for atmospheric turbulence is effectively frozen for images with exposure times shorter than 5−15 ms. Longer exposures effectively perform a time average, and thereby irretrievably wipe out high frequency information, making them bandlimited to angular frequencies smaller than r_o/λ, where λ is the wavelength of the observed light. In contrast, short-exposures encapsulate information up to the diffraction-limited upper frequency bound (which is theoretically given by the ratio D/λ where D denotes the diameter of the telescope’s primary mirror). Figure 1 depicts this issue for the simulated image of a single star and shows the radial averaged modular transfer function (MTF) for diffraction-limited, long- and short-exposure imaging.

Fig. 1 Imaging a single star through atmospheric turbulence: while the angular resolution of a long exposure image (top left) is limited by the strength of the atmospheric turbulence (commonly quantified by Fried’s parameter ro), a short exposure image (bottom left) encapsulates information up to the diffraction-limited upper frequency bound which is proportional to the diameter D of the telescope’s primary mirror. The right panel shows the radial averaged modular transfer function (MTF) for diffraction-limited, short and long exposure imaging. While the long exposure MTF falls to nearly zero at ro/λ, the average short exposure MTF reaches intermediate levels up to the diffraction limited upper frequency bound. The short exposure MTF was averaged over 5000 trials. The simulation was performed for D/ro = 13.3.

The information carried by short exposures was first exploited by Labeyrie (1970), who proposed the averaging of the power spectra of a sequence of short exposures to retain diffraction-limited amplitude information. Shortly thereafter, Knox & Thompson (1974) extended Labeyrie’s idea by suggesting a method for the recovery of the phase information, which is not preserved by Labeyrie’s so-called stellar speckle interferometric method. These early works revolutionized ground-based astronomical observation with large telescopes and have since led to a number of improved signal-processing methods (Lohmann et al. 1983; Stelzer & Ruder 2007) widely referred to as speckle imaging techniques.

An alternative approach was proposed in the seminal work of Ayers & Dainty (1988), who presented a blind deconvolution algorithm (BD) for the problem of atmospherically degraded imaging. BD recovers object information from a blurry and noisy observation without any additional measurement of the distortion. The BD of a single observation is a severely ill-posed problem: there are an infinite number of solutions, and small perturbations of the data result in large deviations in the estimate of the object. The ill-posedness can be alleviated to some degree by confining the set of solutions to physically plausible ones by introducing additional constraints or prior knowledge. Another possibility is to use multiple images or to exploit the partial information about wavefront distortion obtained from wavefront-sensor data, as used in adaptive-optics based myopic deconvolution algorithms.

Since the work of Ayers & Dainty (1988) BD has grown to be a valuable tool in astronomical imaging and has been subject of numerous publications. Today a plethora of algorithms exist that primarily differ in: (i) the data used; (ii) the a-priori knowledge incorporated while deblurring; and (iii) the algorithmic approaches for estimating the object and its blur. For a good overview of BD in the domain of astronomical imaging we refer the reader to Kundur & Hatzinakos (1996); Molina et al. (2001); Pantin et al. (2007).

Recently, electron-multiplying CCD cameras have enabled capturing short-time exposures with negligible noise (Mackay et al. 2001). This in turn has led to a new method: lucky imaging, which can to some degree overcome atmospherically-induced resolution limitations of ground-based telescopes (Law et al. 2006; Oscoz et al. 2008; Hormuth et al. 2008). The lucky imaging idea is based on the work of Fried (1978) (who computed the probability of getting a lucky frame, i.e., an image recorded at a time instant of exceptionally good seeing). This idea proposes to collect only the “best” frames available in a recorded sequence. These “best” frames are subsequently combined to obtain a final image of the object. Usually, out of a thousand images, only a few are selected for the final reconstruction and most of the observed frames are discarded.

This “wastage” can be avoided, and one can indeed use all the frames to obtain an improved reconstruction as we will see in Sect. 5.

Methods for multiframe blind deconvolution (MFBD) aim to recover the image of a fixed underlying object given a sequence of noisy, blurry observations. Each observation has a different and unknown blur, which makes the deconvolution task hard.

Previous approaches to MFBD process all observed frames simultaneously. Doing so limits the total number of frames that can be processed. We show how the computational burden can be greatly reduced by presenting online blind deconvolution (OBD), our online algorithm that processes the input sequence one frame at a time. Each new frame helps to gradually improve the image reconstruction. This simplistic approach is not only natural, but also has several advantages over non-online methods, e.g., lower resource requirements, highly competitive image restoration (Harmeling et al. 2009), low to moderate dependence on regularization or a priori information, and easy extension to super-resolution¹ and saturation-correction.

This paper combines preliminary work (Harmeling et al. 2009, 2010) in the context of astronomical imaging. In particular, the contributions of this paper are as follows:

• (a)

  we show how to incorporate super-resolution whilesimultaneously performing blind deconvolution;

• (b)

  we tackle saturation, a nuisance familiar to anyone who works with astronomical images;

• (c)

  we derive our MFBD algorithm in the framework of stochastic gradient-descent; and

• (d)

  we present results with images taken in a simple astronomer setup, where one does not have access to sophisticated equipment (e.g., adaptive optics), and computational resources might be limited.

Before describing further details, let us put our work into perspective by briefly surveying related work.

2. Related work

MFBD.

A multitude of multi-frame (or multiple-image) deblurring papers discuss the non-blind deconvolution setup, where, in addition to the image sequence the sequence of blur kernels must be known as well. We do not summarize such methods here because ours is a blind deconvolution method. Amongst multiple frame blind approaches, the method of Schulz (1993) is perhaps the earliest. Schulz used penalized likelihood maximization based on a generalized expectation maximization (GEM) framework. Closely related is Li et al. (2004), who also used a GEM framework, but focused on choosing a good objective function and regularizer for optimization. In contrast to our work, both Schulz (1993) and Li et al. (2004) presented batch algorithms that are computationally prohibitive, which greatly limits the number of frames they can simultaneously process.

Sheppard et al. (1998) discussed the MFBD problem and presented a procedure that also processes all frames at the same time. They did, however, mention the possibility of incremental processing of frames, but gave an example only for the non-blind setup. Their blind-deconvolution algorithm was based on conjugate-gradients, for which they had to reparametrize (e.g., x → z²) the variables to enforce nonnegativity. This reparametrization has a long history in image deconvolution (Biraud 1969), but numerically, the ensuing nonlinearity can be damaging as it destroys the convexity of sub-problems.

More recently, Matson et al. (2008) also used the same nonlinear (x → z²) reparametrization for solving MFBD with a parallel implementation of conjugate-gradients. Another approach is that of Zhang et al. (2009), who incorporated a low-pass filter into the MFBD process for suppressing noise, but again at the expense of convexity.

Further MFBD work includes: Anconelli et al. (2006) who considered methods for the reduction of boundary effects; Zhulina (2006) who discussed the Ayers-Dainty algorithm; and Löfdahl (2002) who permitted additional linear inequality constraints. We refer the reader to Matson et al. (2008) for even more references – including those to early works – and a nice summary of blind deconvolution for astronomy. Unlike our algorithm, all the above mentioned blind deconvolution methods are batch procedures; moreover none of them performs either super-resolution or saturation correction.

Super-resolution.

Numerous papers address the standard super-resolution problem. For good surveys we refer the reader to Park et al. (2003); Farsiu et al. (2004). However, most of these works are based on the assumption that the blur is known, and only a few deal with the harder case of blind super-resolution.

The work most closely related to ours is Šroubek et al. (2007), who propose a unifying framework that simultaneously performs blind deconvolution and super-resolution. In Šroubek et al. (2007, 2008) the authors show how a high-resolution image can be obtained from multiple blurry and noise corrupted low-resolution frames. However, their model assumes a priori knowledge about both the image and the blur, and Šroubek et al. (2008) themselves note that their method suffers from numerical instabilities for super-resolution factors larger than 2.5. In contrast, our approach exploits the abundance of available data, which for moderate noise levels does not require imposing any image or blur prior (except nonnegativity), leading to an overall simpler algorithm. Moreover, our method is computationally more efficient, since it is online.

3. The OBD algorithm

3.1. Problem formulation

For simplicity of exposition, our description will focus on one-dimensional images and point spread functions (PSFs). In Appendix A we cover the generalization to two-dimensions.

Let each observed (blurry and noisy) frame be denoted by y_t, the “true” unknown image by x, and each unknown PSF by f_t. Then, we use the observation model (1)where f_t ∗ x represents convolution (circular or non-circular), and n_t denotes measurement noise. Further, on physical grounds we assume both the image x and the PSF f_t to be nonnegative.

3.2. Algorithm

First consider the case where given the next observation y_t and the current image estimate x_t, we wish to compute the PSF f_t. Assuming the noise n_t in Eq. (1) to be Gaussian distributed with zero mean and incorporating nonnegativity, the PSF f_t can be determined by solving a nonnegative least-squares (NNLS) problem². For a given observation frame y_t and a current estimate x_t, we define the loss (2)For a frame sequence y₁,y₂,...,y_T, we aim to minimize the overall loss by computing the image x that solves (3)Problem (3) is not easy, because it is non-convex and its optimal solution requires computing both x as well as the PSFs f₁,...,f_T. Nevertheless, given our formulation, several methods could potentially be used for minimizing L_T(x). For example, an ordinary gradient-projection scheme would be (4)where P₊ denotes projection onto the nonnegative orthant; x_t denotes the current image estimate; and α_t is an appropriate step-size. However, when the number of frames T is large, such an approach rapidly becomes computationally impractical. Hence we turn to a simpler method that processes the input one frame at a time.

3.3. Stochastic gradient descent

A simple and often effective method for minimizing the overall loss in Eq. (3) is stochastic gradient descent (SGD). This method does not process all the frames simultaneously, but at step t it picks (at random) some frame y and updates the current image estimate x_t as (5)where P₊ and α_t are as before; computing ∇ℓ(y;x_t) requires solving Eq. (2). By processing only one frame at a time, SGD leads to huge computational savings. However, there are two main difficulties: update rule (5) converges slowly; and more importantly, it is sensitive to the choice of the step-size α_t; a popular choice is α_t = β/(t₀ + t), where the constants t₀ and β must be tuned empirically.

We propose a practical modification to the step-size computation, wherein we instead use the scaled-gradient version (6)where S_t is a positive-definite matrix. Also update rule (6) can be shown to converge³ under appropriate restrictions on α_t and S_t (Kushner & Yin 2003; Bottou 1998). In general, the matrix S_t is chosen to approximate the inverse of the Hessian of L_T(x^∗) for an optimal x^∗, thereby yielding quasi-Newton versions of SGD. But a more straightforward choice is given by the diagonal matrix (7)where the Diag operator maps a vector x to a diagonal matrix with elements of x along its diagonal. Also note that the division in (7) is element-wise, F_t is the matrix representation of the PSF f_t (see Appendix A), and ε > 0 is a positive constant which ensures that S_t remains positive definite and bounded (both requirements are crucial for convergence of the method). The choice (7) can be motivated with the help of auxiliary functions (e.g., as in Harmeling et al. 2009).

Remark: We note in passing that if one were to use α_t = 1, and set ε = 0, then although convergence is no longer guaranteed, iteration (6) takes a particularly simple form, namely, (8)where  ⊙  denotes the Hadamard (elementwise) product of two vectors – this update may be viewed as an online version of the familiar ISRA (see Daube-Witherspoon & Muehllehner 1986).

Note that for (7) the matrix F corresponds to the PSF f computed via the NNLS problem (2) with y and x = x_t. We call the method based on iteration (6) online blind deconvolution (OBD) and provide pseudo-code as Algorithm 1. We further note that by assuming photon shot noise (Poisson-distributed) in Eq. (1) instead of additive noise, we can also design a Richardson-Lucy type iteration for solving Eq. (3).

Algorithm 1: Online blind deconvolution (OBD).

4. Extending OBD

4.1. Super-resolution

In the OBD setup an entire sequence of frames is at our disposal. Can we exploit this sequence to improve the image reconstruction beyond mere blind deconvolution? The answer is “yes”. With a small increase in computational costs we can augment the basic algorithm and perform super-resolution. For long-exposures that often lose higher-frequency structure (finer details) of the image due to averaging, such increased resolution is particularly desirable.

To incorporate super-resolution into our framework we introduce the resizing matrix (11)where I_n is the n × n identity matrix, 1_n is an n dimensional column vector of ones, and  ⊗  denotes the Kronecker product. The matrix  transforms a vector v of length m into a vector of length n. The sum of v’s entries is preserved (formally verified by applying identity (11) twice). This is a favorable property for images, as the number of photons observed should not depend on the resolution. Note that even if the sizes n and m are not multiples of each other,  will interpolate appropriately. Hence, the super-resolution factor, i.e., the ratio n/m, is not restricted to be integral. Note that for m ≥ n, i.e. downscaling, the matrix operation corresponds to integrating neighboring pixels weighted by their overlap with the new pixel grid. Similarly for m < n, i.e. upscaling, the operation will take the nearest neighbor, if n is divisible by m, or a weighted linear combination of closeby pixels.

To avoid double indexing let n = l_y be the length of y. For super-resolution by a factor of s we choose x and f large enough such that the vector f ∗ x has length sn. Then we replace the loss ℓ(y_t;x) by (cf. Eq. (2)) (12)For this loss, a derivation similar to that for (7) yields the diagonal matrix (13)where F_t corresponds to f_t obtained by solving (12).

4.2. Overexposed pixels

For astronomical images a common problem is saturation of pixels due to overexposure, i.e., some pixels receive so many photons that they exceed the peak intensity permitted by the hardware. This saturation can be particularly confounding if both bright and faint stars are present in the same image, especially when some stars are orders of magnitude brighter. Overexposed pixels impede not only deblurring but also super-resolution and applications such as estimation of star magnitudes.

However, since we have an entire sequence of observed frames, tackling overexposed pictures is feasible. Here the atmospheric blurring proves to be helpful, since it creates non-overexposed margins around a bright star whose center pixels are overexposed. Our method is able to fit these margins and can approximate the true star magnitude. Our approach essentially consists of identifying saturated pixels and excluding them from the computation of the objective function. This approach might seem to be overly simple, but its success is deeply tied to the availability of multiple frames. Specifically, since each frame can have different pixels attaining saturation (different frames are aligned differently), we have to check at each iteration which pixels in the current image are saturated. To ignore these pixels we define a diagonal weighting matrix (per frame) with entries, (14)along its diagonal. Hereby, we assume the value of a saturated pixel to be ρ_max (e.g. in the case of 16 bit images, ρ_max = 65   535). We can modify the updates to ignore saturated pixels by replacing the Euclidean norm with a weighted norm . We replace the loss ℓ(y_t;x) by (15)For this loss, following a derivation similar to (7) yields the diagonal matrix (16)where as before F_t corresponds to f_t obtained by solving (12).

Remark.

One might ask whether we can recover pixels in x that are saturated in most of the frames? The answer is yes, and can be understood as follows. The photons corresponding to such a pixel in x are spread by the PSF across a whole set of pixels in each observed frame. Thus, if not all these pixels are always saturated, the true value for the corresponding pixel in x can be recovered.

5. Results on simulated data

To investigate how our algorithms performs on atmospherically degraded short-exposure images, we first experiment in a controlled setting with simulated data.

Following Harding et al. (1999), we generate a sequence of 200 PSFs with Kolmogorov random phase screens at a specified ratio D/r_o of the telescope diameter to the atmospheric seeing parameter (Fried parameter) equal to 13.3. The strength of the turbulence is chosen to create images recorded by a 26-inch telescope through atmospheric turbulence of a coherence length of approximately r_o = 5 cm.

Figure 2 shows the original object, one out of the 200 PSFs, and the noise-free short exposure image obtained by convolving the shown PSF with the object. The object is a rendered model of the satellite OCNR5 used by Sheppard et al. (1998) and was chosen because of its high dynamic range and its great level of detail.

Fig. 2 Simulation: from left to right: original object image of OCNR5, typical PSF, blurred image.

Before corrupting the images with noise, we add a constant background b to the blurred image f_t ∗ x. To simulate photon noise we scale the pixel values (ranging between 0 and 255) of each short exposure to varying large numbers of photons, i.e. λ(f_t ∗ x + b) and sample a new image z from the corresponding Poisson distribution, i.e. (17)For differing λ we can hereby simulate differing amounts of photon shot noise. After scaling down by 1/λ, we add white Gaussian noise with zero mean and a variance σ² equal to two percent of the maximal image intensity of the whole sequence to model the readout noise common to CCD cameras, To quantify the amount of image noise we define the following signal-to-noise ratio, (20)where x denotes the true satellite image, y_t the noise-corrupted atmospherically degraded observation, and f_t the PSF, respectively. Var(x) denotes the variance of the pixel values of x. For an entire sequence y₁, y₂, ..., y₂₀₀ we average over the computed SNRs of all 200 frames, . Table 1 shows the computed SNR for different parameter settings that we use in our experiments. Note that we use the SNR only to quantify the amount of noise in the simulated data. To measure the quality of the reconstruction we use relative error (explained below).

Figure 3 shows typical frames for different SNRs, each 256  ×  256 pixels in size, and the reconstructed object images of our basic algorithm after having processed all 200 frames within one sequence. The restored images shown are cropped to the size of the observations. As initial estimates for the PSFs we chose constant images of size 60  ×  60 pixels, and as the initial estimate of the object, an average over the first twenty observed frames embedded in a 315  ×  315 array of zeros.

As expected, the quality of the reconstruction suffers as the SNR decreases, which is also reflected quantitatively in Fig. 4, where we plot the relative error  of the reconstructed image  as a function of observed frames and the corresponding SNR.

Evidently, for high SNRs the reconstruction error decreases the more observations have been processed and saturates to a certain value dependent on the SNR. The error is higher the lower the SNR of the available observations. The error does not decrease strictly monotonically from frame to frame, but more in a (long-term) stochastic gradient manner. As expected, for lower SNRs, the unregularized reconstruction process can even diverge. In this noisy regime, additional prior knowledge about the object is necessary and regularization in the restoration process is inevitable.

Figure 5 illustrates that enforcing smoothness by employing Tikhonov regularization on the gradients of the reconstructed image (i.e. a prior term η ∥ ∇x ∥² is added to the loss in (2)) is capable of suppressing noise amplification and stabilizing the deconvolution process even for low SNRs. As expected, when the regularization parameter η is too small, the reconstruction error still diverges (red dotted curve); similarly, when it is too large, the error is increased due to oversmoothing (blue dashed curve). A reasonable choice of the regularization parameter may be obtained by setting it proportional to the noise variance. The color framed image stamps show the reconstruction results for different values of the regularization parameter.

Fig. 3 Simulation: typical observed frames (top row) and reconstructed images (bottom row) after having processed a sequence of 200 blurred frames for different SNRs.

To study the influence of the initialization and the order of frames within one sequence, we reversed and randomly permuted the processing order of the input frames. Figure 6 shows restored object images and the corresponding error curves for a fixed SNR of 18.9 dB, respectively. As can be seen, the error evolution of the deconvolution process is almost independent of the particular ordering of the input frames. All curves converge to a similar value with small variance, and visually, only little (if at all) difference is discernible.

To numerically appraise the quality of our results, we did a quantitative comparison with various state-of-the-art reconstruction methods. Figure 7 shows the visually best observed frame, a reconstruction with AviStack (Theusner 2009), a popular Lucky Imaging software. AviStack partitions the images into small image patches of variable sizes, evaluates the quality of all observed frames for all image patches and then aligns and stacks those image patches, that fulfill a certain quality threshold. For the final reconstruction only the best percent of observed frames was taken. Next to it, a Knox-Thompson reconstruction is shown, which was obtained using Speckle1, a reconstruction software by Stelzer (2009). For the reconstruction, 300 Knox-Thompson and 100 triple correlation phase pairs were used. Finally, the rightmost image shows the result of our basic algorithm without any additional regularisation. In all cases no further post-processing was performed.

Fig. 4 Simulation: evaluation of relative reconstruction error for different SNRs.

Fig. 5 Simulation: evaluation of the relative reconstruction error for different values of the regularization constant η at a fixed SNR of 4.6 dB.

Fig. 6 Simulation: evaluation of the relative reconstruction error and final reconstructed images after having processed 200 frames in chronological, reverse and various random orders at a fixed SNR of 18.9 dB.

Fig. 7 Simulation: final reconstructed images after having processed 200 frames at a fixed SNR of 18.9 dB with Avistack, Knox-Thompson and our proposed method. The relative reconstruction error is overlayed in white.

For a single isoplanatic patch the reconstruction with AviStack is not substantially better than the visually best observed frame, which is also reflected in the relative error overlayed in white. In comparison, both the Knox-Thompson reconstruction and the result by the basic algorithm of our proposed method reveal much greater detail and higher spatial resolution. Subjectively, our result is comparable in quality and resolution to the Knox-Thompson reconstruction, which is quantitatively confirmed by the negligible difference in the reconstruction error. Regarding runtime, the C implementation of Stelzer (2009) takes about 15 min (when invoked carefully by an expert user) for the entire reconstruction on a single core of an Intel(R) Core(TM) i5 processor with 2.67 GHz. Our Matlab implementation that is however not optimized for speed and logs large quantities of intermediate results, takes about thrice as long. A Python implementation using PyCUDA (Klöckner et al. 2009) for GPU enabled computation of the discrete Fourier transform (see Eqs. (A.2) and (A.3)) achieves a runtime of less than 10 min on a low-cost NVIDIA(R) GeForce(TM) GT 430.

Our final set of experiments with simulated data evaluates our algorithm’s super-resolution abilities. We generated three sequences of atmospherically blurred, differently downsampled and noisy observations at a fixed SNR of 18.9 dB. Panel A of Fig. 8 shows typical input images of these sequences together with their corresponding downsampling factors. On each of these three simulations we ran our algorithm with various super-resolution factors. The results are shown in Panel B of Fig. 8. The relative errors overlayed in white are computed by linearly interpolating the reconstructed images to the size of the ground truth image. The numbers suggest that incorporating super-resolution does improve the results of the deconvolution beyond mere interpolation, which validates the merits of our approach.

6. Results on astronomical data

We now present results of our algorithm on a variety of actual astronomical data. Some of the images were taken with an off-the-shelf 12-inch f/10 MEADE LX200 ACF Schmidt-Cassegrain telescope, some with the 24-inch f/8 Hypergraph Cassegrain telescope “Ganymed” of the Capella Observatory located on Mount Skinikas in Crete, Greece. The data consists of short-exposure imagery of star constellations, the lunar Copernicus crater, as well as long-exposure deep-sky images. We compare our results against state-of-the-art methods used by both amateur and professional astronomers, and show that our method yields competitive if not superior results in all case studies.

6.1. Binary star

The first dataset is an image sequence of the binary star system Epsilon Lyrae 2 of the constellation Lyra with an angular separation of 2.3′′ and a relative magnitude of 1.08. As we know precisely what to expect, our results on this dataset serve as an additional proof of concept. The sequence consists of 300 frames, each 132  ×  112 pixels in size, taken with a 24-inch Cassegrain telescope at the Capella Observatory and an Imaging Source DFK 31BU03 CCD camera; the image scale was 0.06′′ per pixel. The seeing was estimated to FWHM  ≈  0.8′′, corresponding to a Fried parameter of r_o ≈ 20   cm at a wavelength of λ = 500   nm.

Figure 10 shows in three columns the first four and the last two frames of the processed sequence and illustrates schematically how our method works. Each row shows from left to right the observed image y_t, the corresponding PSF f_t estimated by our algorithm and the current estimate x_t of the true image we want to recover. The PSF is chosen to be of size 30  ×  30 pixels.

The image x is initialized by the first observed frame y₁. Then f₂ is estimated from the second observed frame y₂ and the current estimate of x. After that we improve the estimate of x by means of (10) and proceed with the next observed frame. Figure 10 shows nicely that already after 40 frames we obtain a good reconstruction.

Fig. 8 Simulation: final reconstructed images (panel B)) after having processed 200 frames for differently downsampled input images (panel A)) and various super-resolution factors at a fixed SNR of 18.9 dB. The downsampling and super-resolution factor is abbreviated with DF and SR respectively. The displayed number corresponds to the relative reconstruction error.

Figure 9 shows an enlarged version of the result of our algorithm after 300 iterations along with the estimated PSFs for each color channel. Note how blurry the observed image y_t is (left), while our estimate of x is almost free of any degradation (middle). Furthermore, we see that both stars have almost identical diffraction patterns which strongly resemble the estimated PSFs (shown on the right for each color channel). This finding justifies our assumption about a constant PSF for the whole image. From the final reconstructed image we determined a separation of 2.28′′ and a magnitude ratio of 1.08, which is in excellent accordance with the literature.

6.2. Copernicus crater

To evaluate our algorithm on an extended celestial object, we applied it to a sequence of short exposures of the Copernicus crater, a prominent lunar crater located in eastern Oceanus Procellarum. The original recording was taken with a 14-inch f/10 Celestron C14 and a DMK 31 AF03 CCD camera from Imaging Source at a frame rate of 30 fps near Frankfurt, Germany (courtesy Mario Weigand). It consists of 2350 frames in total, where each frame is 1024  ×  768 pixels in size. To meet our assumption of a constant PSF, we processed only a small image patch of 70  ×  70 pixels, which corresponds to a angular size of 0.92′′. In this field of view the PSF is assumed to be constant, which is a valid assumption for the seeing conditions at the time of recording.

The top row of Fig. 11 shows the selected region of the central peak in the Copernicus crater and typical observed frames. The image patches were aligned on a pixel scale before processing to reduce computational costs⁴. For reconstruction all 2350 observed frames were taken into account.

The bottom row of Fig. 11 shows a comparison of different reconstruction methods. Panel (a) of Fig. 11 shows the visually best observed frame, Panel (b) a reconstruction with AviStack (Theusner 2009), for which the best ten frames were taken into account. In Panel (c) a Knox-Thompson reconstruction is shown, which was done with Stelzer (2009) using 300 Knox-Thompson and 100 triple correlation phase pairs. Finally, Panel (d) shows the result of our basic algorithm and Panel (e) the result two times super-resolved. In all cases no further post-processing was performed.

As before, within a single isoplanatic patch the result of AviStack seems to be not considerably better than the visually best observed frame. In contrast, the Knox-Thompson reconstruction reveals greater detail and higher spatial resolution. Subjectively, our result is comparable if not superior in quality and resolution to the Knox-Thompson reconstruction. The two times super-resolved reconstruction seems to reveal even more detail.

6.3. Orion Trapezium

In this experiment, we used a 12-inch f/10 Meade LX200 ACF Schmidt-Cassegrain telescope and a AVT PIKE F-032B uncooled CCD camera to record a short video (191 frames acquired at 120 fps) of the Trapezium in the constellation Orion. The exposure time of the individual frames was sufficiently short to “freeze” the atmospheric turbulence and thus retain the high-frequency information which is present in the atmospheric PSF – see Fig. 12 for sample frames.

The Orion Trapezium is formed by four stars ranging in brightness from magnitude 5 to magnitude 8, with angular separations around 10′′ to 20′′. Here it should be mentioned, that our assumption of a constant PSF throughout the field of view is strongly violated. However, by resorting to early stopping in this case, we avoid overfitting the PSF. The first row of Fig. 13 shows from left to right (a) an enlarged unprocessed frame; (b) the deconvolution results obtained by the basic algorithm; (c) the result using the proposed method to handle saturation; and (d) the results if we additionally apply the proposed method for four times super-resolution. The bottom row shows a closeup of the brightest star within the Trapezium. Panel (e) of Fig. 13 shows the star profiles obtained by slicing as indicated by the colored lines in the image stamps (a)−(d).

Fig. 9 Binary star system Epsilon Lyrae 2: typical observed image y300 (left), reconstruction x after 300 iterations (middle), estimated PSFs f300 for each color channel. Note the subtle differences in the PSFs due to wavelength dependent diffraction. Hence, the color channels are not perfectly aligned.

Fig. 10 Binary star system Epsilon Lyrae 2: schematic illustration of the temporal evolution. From left to right: observed image yt (left), estimate of the corresponding PSF ft and reconstruction xt after t timesteps.

An important application in astronomy is the measurement of the brightness of stars and other celestial objects (photometry). To this end, a linear sensor response is required (for our purposes, the used CCD sensor may be assumed linear). The intensity counts can then be translated into stellar magnitudes. Clearly, this is not directly possible for stars that saturate the CCD (i.e., where so many photons are recorded that the capacities of the pixels are exceeded). However, we can use the proposed method for deconvolution with saturation correction and reconstruct the photon counts (image intensities) that would have been recorded had the pixels not been saturated; then we convert these into ratios between star intensities, i.e. differences between stellar magnitudes. For the difference between two star magnitudes, we use the formula m₁ − m₂ = −2.5log ₁₀p₁/p₂ where p₁ and p₂ are the pixel values of two stars in the reconstructed image. We do this for all Trapezium stars relative to the brightest star C and obtain encouraging results (see Table 2).

Fig. 11 Copernicus Crater: top panel: full frame with extracted image patch marked by white square (left) and example sequence of 12 observed frames. Bottom panel: comparison of results of different reconstruction algorithms (from left to right): visually best frame, AviStack, Knox-Thompson, our approach and our approach with two times super-resolution. All image results are shown without any post-processing. This figure is best viewed on screen, rather than in print.

Fig. 12 Orion Trapezium Cluster: example sequence of observed frames, y1, ..., y10.

Fig. 13 Orion Trapezium Cluster (from left to right): a) the first observed frame; b) x191 for basic algorithm; c) x191 for saturation corrected; and d) x191 for saturation corrected and four times super-resolved. Top row shows the overall trapezium; bottom row shows the brightest star enlarged. Panel e) shows the stellar profiles at the positions indicated by the coloured lines in plots a) − d).

Fig. 14 Globular cluster M13: (left) example observed frame, (right) result of saturation corrected, two times super-resolved multi-frame blind deconvolution; (top) overview, (bottom) closeup. For better display the images have been automatically gamma corrected.

6.4. Globular cluster M13

M13 is a globular cluster in the constellation Hercules, around 25 000 light years away, with an apparent size of around 20′. It contains several 100 000 stars, the brightest of which has an apparent magnitude of 12. Such faint stars cannot be imaged using our equipment for short exposures; however, long exposures with budget equipment typically incur tracking errors, caused by telescope mounts that do not perfectly compensate for the rotation of the earth. In our case, the tracking errors induced a significant motion blur in the images, which we attempted to remove using the same algorithm that we used above on short exposures. All raw images were recorded using a 12-inch f/10 MEADE LX200 ACF Schmidt-Cassegrain telescope and a Canon EOS 5D digital single lens reflex (DSLR) camera. The whole sequence consists of 26 images with an exposure time of 60 s each. The top row of Fig. 14 displays a long exposure with motion blur (left panel and the twice super-resolved result of our algorithm (right) applied to 26 motion degraded frames. In the bottom row we clearly see details in our reconstructed image (right) which where hidden in the recorded frames (left). However, note that in the bottom right panel there appear also some JPEG-like artifacts which might suggest that 26 frames were not enough for two times super-resolution.

7. Conclusions and future work

In this paper, we proposed a simple, efficient, and effective multiframe blind deconvolution algorithm. This algorithm restores an underlying static image from a stream of degraded and noisy observations by processing the observations in an online fashion. For moderate signal-to-noise ratios our algorithm does not depend on any prior knowledge other than nonnegativity of the PSFs and the images. Thus, in a sense our reconstruction is unbiased since no specific image model is enforced. Moreover, our formulation exploits the availability of multiple frames to incorporate super-resolution and saturation-correction.

We showed results on both simulated and real world astronomical data to verify and demonstrate the performance of our algorithm. We experimented with not only short-exposure images where the degradation is caused by atmospheric turbulence, but also with long exposures that suffer from saturation and additional blur arising from mechanical inaccuracies in the telescope mount. Our method yields results superior to or at worst comparable to existing frequently used reconstruction methods.

Future work includes further building on the simplicity of our method to improve it to work in real-time. This goal might be achievable by exploiting fast graphical processing unit (GPU) based computing. First attempts already yielded promising results in terms of speedup (see Sect. 5). Beyond computing improvements, two other important aspects are: (i) to explore the spatio-temporal properties of the speckle pattern; and (ii) to incorporate and investigate additional regularization within the reconstruction process. The most challenging subject of future investigation is to extend our method to space-varying PSFs.

Acknowledgments

The authors thank Hanns Ruder, Josef Pöpsel and Stefan Binnewies for fruitful discussions and for their generosity regarding observation time at the Capella Observatory. Many thanks also to Karl-Ludwig Barth (IAS) and Mario Weigand. We finally do thank the reviewer for a thorough review and valuable comments which greatly helped to improve the quality of this article. A Matlab demo with graphical user interface is available at http://pixel.kyb.tuebingen.mpg.de/obd. Corresponding author Michael Hirsch can be reached by e-mail at michael.hirsch@tuebingen.mpg.de.

References

Appendix A: Implementation details

Although in Sects. 3 and 4 we only considered vectors, one-dimensional convolutions, and vector-norms, all results naturally generalize to two-dimensional images. However, efficiently implementing the resulting algorithms for two-dimensional images requires some care and handling of technical details.

A.1. Convolution as matrix-vector multiplication

We introduced f ∗ x as the convolution, which could be either circular or non-circular. Due to linearity and commutativity, we can also use matrix-vector notation to write (A.1)The matrices F and X are given by Matrix W is the discrete Fourier transform matrix, i.e., Wx is the Fourier transform of x. The diagonal matrix Diag(v) has vector v along its diagonal, while I_x, I_f, and I_y are zero-padding matrices which ensure that I_xx, I_ff, and I_yy have the same length. Different choices of the matrices lead to different margin condition of the convolution.

For two-dimensional images and PSFs we have to consider two-dimensional Fourier transforms, which can be written as left- and right-multiplications with W, and represented as a single matrix-vector multiplication using Kronecker products and the vectorization operator vec(x), which stacks columns of the two-dimensional image x into a one-dimensional vector in lexicographical order; formally, (A.4)which follows from the identity (Horn & Johnson 1991) (A.5)The zero-padding operations for two-dimensional images can be written in a similar way.

A.2. Resizing matrices

The resizing matrix can be implemented efficiently using sparse matrices⁵. Resizing two-dimensional images can also be implemented by left- and right-multiplications: let x be an m × n image, then is an image of size p × q. Using Eq. (A.5) we can write this operation as the matrix-vector product (A.6)

All Tables

All Figures

Fig. 1 Imaging a single star through atmospheric turbulence: while the angular resolution of a long exposure image (top left) is limited by the strength of the atmospheric turbulence (commonly quantified by Fried’s parameter ro), a short exposure image (bottom left) encapsulates information up to the diffraction-limited upper frequency bound which is proportional to the diameter D of the telescope’s primary mirror. The right panel shows the radial averaged modular transfer function (MTF) for diffraction-limited, short and long exposure imaging. While the long exposure MTF falls to nearly zero at ro/λ, the average short exposure MTF reaches intermediate levels up to the diffraction limited upper frequency bound. The short exposure MTF was averaged over 5000 trials. The simulation was performed for D/ro = 13.3.

In the text

	Fig. 2 Simulation: from left to right: original object image of OCNR5, typical PSF, blurred image.
In the text

	Fig. 3 Simulation: typical observed frames (top row) and reconstructed images (bottom row) after having processed a sequence of 200 blurred frames for different SNRs.
In the text

	Fig. 4 Simulation: evaluation of relative reconstruction error for different SNRs.
In the text

	Fig. 5 Simulation: evaluation of the relative reconstruction error for different values of the regularization constant η at a fixed SNR of 4.6 dB.
In the text

	Fig. 6 Simulation: evaluation of the relative reconstruction error and final reconstructed images after having processed 200 frames in chronological, reverse and various random orders at a fixed SNR of 18.9 dB.
In the text

	Fig. 7 Simulation: final reconstructed images after having processed 200 frames at a fixed SNR of 18.9 dB with Avistack, Knox-Thompson and our proposed method. The relative reconstruction error is overlayed in white.
In the text

	Fig. 8 Simulation: final reconstructed images (panel B)) after having processed 200 frames for differently downsampled input images (panel A)) and various super-resolution factors at a fixed SNR of 18.9 dB. The downsampling and super-resolution factor is abbreviated with DF and SR respectively. The displayed number corresponds to the relative reconstruction error.
In the text

	Fig. 9 Binary star system Epsilon Lyrae 2: typical observed image y300 (left), reconstruction x after 300 iterations (middle), estimated PSFs f300 for each color channel. Note the subtle differences in the PSFs due to wavelength dependent diffraction. Hence, the color channels are not perfectly aligned.
In the text

	Fig. 10 Binary star system Epsilon Lyrae 2: schematic illustration of the temporal evolution. From left to right: observed image yt (left), estimate of the corresponding PSF ft and reconstruction xt after t timesteps.
In the text

Fig. 11 Copernicus Crater: top panel: full frame with extracted image patch marked by white square (left) and example sequence of 12 observed frames. Bottom panel: comparison of results of different reconstruction algorithms (from left to right): visually best frame, AviStack, Knox-Thompson, our approach and our approach with two times super-resolution. All image results are shown without any post-processing. This figure is best viewed on screen, rather than in print.

In the text

	Fig. 12 Orion Trapezium Cluster: example sequence of observed frames, y1, ..., y10.
In the text

	Fig. 13 Orion Trapezium Cluster (from left to right): a) the first observed frame; b) x191 for basic algorithm; c) x191 for saturation corrected; and d) x191 for saturation corrected and four times super-resolved. Top row shows the overall trapezium; bottom row shows the brightest star enlarged. Panel e) shows the stellar profiles at the positions indicated by the coloured lines in plots a) − d).
In the text

	Fig. 14 Globular cluster M13: (left) example observed frame, (right) result of saturation corrected, two times super-resolved multi-frame blind deconvolution; (top) overview, (bottom) closeup. For better display the images have been automatically gamma corrected.
In the text