Issue 
A&A
Volume 531, July 2011



Article Number  A9  
Number of page(s)  11  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/200913955  
Published online  31 May 2011 
Online multiframe blind deconvolution with superresolution and saturation correction
Max Planck Institute for Intelligent Systems,
Spemannstraße 38,
72076
Tübingen,
Germany
email: firstname.lastname@tuebingen.mpg.de
Received: 23 December 2009
Accepted: 21 February 2011
Astronomical images taken by groundbased telescopes suffer degradation due to atmospheric turbulence. This degradation can be tackled by costly hardwarebased approaches such as adaptive optics, or by sophisticated softwarebased methods such as lucky imaging, speckle imaging, or multiframe deconvolution. Softwarebased methods process a sequence of images to reconstruct a deblurred highquality image. However, existing approaches are limited in one or several aspects: (i) they process all images in batch mode, which for thousands of images is prohibitive; (ii) they do not reconstruct a superresolved image, even though an image sequence often contains enough information; (iii) they are unable to deal with saturated pixels; and (iv) they are usually nonblind, i.e., they assume the blur kernels to be known. In this paper we present a new method for multiframe deconvolution called online blind deconvolution (OBD) that overcomes all these limitations simultaneously. Encouraging results on simulated and real astronomical images demonstrate that OBD yields deblurred images of comparable and often better quality than existing approaches.
Key words: methods: numerical / techniques: image processing / techniques: miscellaneous / atmospheric effects
© ESO, 2011
1. Introduction
Astronomical observation using groundbased telescopes is significantly degraded by diffractionindex fluctuations caused by atmospheric turbulence. This turbulence arises from local temperature and density inhomogeneities and results in a time and spacevariant point spread function (PSF). Often the PSF is assumed to be invariant within a short timeperiod 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, shortexposures encapsulate information up to the diffractionlimited 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 diffractionlimited, long and shortexposure 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 r_{o}), a short exposure image (bottom left) encapsulates information up to the diffractionlimited 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 diffractionlimited, short and long exposure imaging. While the long exposure MTF falls to nearly zero at r_{o}/λ, 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/r_{o} = 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 diffractionlimited 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 socalled stellar speckle interferometric method. These early works revolutionized groundbased astronomical observation with large telescopes and have since led to a number of improved signalprocessing 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 illposed 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 illposedness 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 wavefrontsensor data, as used in adaptiveoptics 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 apriori 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, electronmultiplying CCD cameras have enabled capturing shorttime 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 atmosphericallyinduced resolution limitations of groundbased 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 nononline 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 superresolution^{1} and saturationcorrection.
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 superresolution 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 gradientdescent; 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 multiframe (or multipleimage) deblurring papers discuss the nonblind 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 nonblind setup. Their blinddeconvolution algorithm was based on conjugategradients, for which they had to reparametrize (e.g., x → z^{2}) 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 subproblems.
More recently, Matson et al. (2008) also used the same nonlinear (x → z^{2}) reparametrization for solving MFBD with a parallel implementation of conjugategradients. Another approach is that of Zhang et al. (2009), who incorporated a lowpass 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 AyersDainty 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 superresolution or saturation correction.
Superresolution.
Numerous papers address the standard superresolution 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 superresolution.
The work most closely related to ours is Šroubek et al. (2007), who propose a unifying framework that simultaneously performs blind deconvolution and superresolution. In Šroubek et al. (2007, 2008) the authors show how a highresolution image can be obtained from multiple blurry and noise corrupted lowresolution 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 superresolution 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 onedimensional images and point spread functions (PSFs). In Appendix A we cover the generalization to twodimensions.
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 noncircular), 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 leastsquares (NNLS) problem^{2}. For a given observation frame y_{t} and a current estimate x_{t}, we define the loss (2)For a frame sequence y_{1},y_{2},...,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 nonconvex and its optimal solution requires computing both x as well as the PSFs f_{1},...,f_{T}. Nevertheless, given our formulation, several methods could potentially be used for minimizing L_{T}(x). For example, an ordinary gradientprojection scheme would be (4)where P_{+} denotes projection onto the nonnegative orthant; x_{t} denotes the current image estimate; and α_{t} is an appropriate stepsize. 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 stepsize α_{t}; a popular choice is α_{t} = β/(t_{0} + t), where the constants t_{0} and β must be tuned empirically.
We propose a practical modification to the stepsize computation, wherein we instead use the scaledgradient version (6)where S_{t} is a positivedefinite matrix. Also update rule (6) can be shown to converge^{3} 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 quasiNewton 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 elementwise, 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 DaubeWitherspoon & 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 pseudocode as Algorithm 1. We further note that by assuming photon shot noise (Poissondistributed) in Eq. (1) instead of additive noise, we can also design a RichardsonLucy type iteration for solving Eq. (3).
Algorithm 1: Online blind deconvolution (OBD).
4. Extending OBD
4.1. Superresolution
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 superresolution. For longexposures that often lose higherfrequency structure (finer details) of the image due to averaging, such increased resolution is particularly desirable.
To incorporate superresolution 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 superresolution 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 superresolution 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 superresolution 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 nonoverexposed 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 shortexposure 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 26inch 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 noisefree 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 σ^{2} 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 signaltonoise ratio, (20)where x denotes the true satellite image, y_{t} the noisecorrupted 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_{1}, y_{2}, ..., y_{200} 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).
SNRs for different parameter settings of λ and σ^{2}.
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 (longterm) 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 ∥^{2} 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 stateoftheart 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 KnoxThompson reconstruction is shown, which was obtained using Speckle1, a reconstruction software by Stelzer (2009). For the reconstruction, 300 KnoxThompson 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 postprocessing 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, KnoxThompson 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 KnoxThompson 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 KnoxThompson 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 lowcost NVIDIA(R) GeForce(TM) GT 430.
Our final set of experiments with simulated data evaluates our algorithm’s superresolution 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 superresolution 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 superresolution 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 offtheshelf 12inch f/10 MEADE LX200 ACF SchmidtCassegrain telescope, some with the 24inch f/8 Hypergraph Cassegrain telescope “Ganymed” of the Capella Observatory located on Mount Skinikas in Crete, Greece. The data consists of shortexposure imagery of star constellations, the lunar Copernicus crater, as well as longexposure deepsky images. We compare our results against stateoftheart 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 24inch 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_{1}. Then f_{2} is estimated from the second observed frame y_{2} 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 superresolution factors at a fixed SNR of 18.9 dB. The downsampling and superresolution 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 14inch 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^{4}. 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 KnoxThompson reconstruction is shown, which was done with Stelzer (2009) using 300 KnoxThompson and 100 triple correlation phase pairs. Finally, Panel (d) shows the result of our basic algorithm and Panel (e) the result two times superresolved. In all cases no further postprocessing 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 KnoxThompson reconstruction reveals greater detail and higher spatial resolution. Subjectively, our result is comparable if not superior in quality and resolution to the KnoxThompson reconstruction. The two times superresolved reconstruction seems to reveal even more detail.
6.3. Orion Trapezium
In this experiment, we used a 12inch f/10 Meade LX200 ACF SchmidtCassegrain telescope and a AVT PIKE F032B 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 highfrequency 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 superresolution. 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 y_{300} (left), reconstruction x after 300 iterations (middle), estimated PSFs f_{300} 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 y_{t} (left), estimate of the corresponding PSF f_{t} and reconstruction x_{t} 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_{1} − m_{2} = −2.5log _{10}p_{1}/p_{2} where p_{1} and p_{2} 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, KnoxThompson, our approach and our approach with two times superresolution. All image results are shown without any postprocessing. This figure is best viewed on screen, rather than in print. 
Fig. 12 Orion Trapezium Cluster: example sequence of observed frames, y_{1}, ..., y_{10}. 
Fig. 13 Orion Trapezium Cluster (from left to right): a) the first observed frame; b) x_{191} for basic algorithm; c) x_{191} for saturation corrected; and d) x_{191} for saturation corrected and four times superresolved. 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 superresolved multiframe blind deconvolution; (top) overview, (bottom) closeup. For better display the images have been automatically gamma corrected. 
True star magnitudes (note that stars A and B have variable magnitudes), true differences to star C, and estimated difference values estimated after deconvolution without and with saturation correction.
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 12inch f/10 MEADE LX200 ACF SchmidtCassegrain 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 superresolved 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 JPEGlike artifacts which might suggest that 26 frames were not enough for two times superresolution.
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 signaltonoise 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 superresolution and saturationcorrection.
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 shortexposure 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 realtime. 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 spatiotemporal 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 spacevarying PSFs.
This NNLS problem may be solved by various methods; we used the LBFGSB algorithm (Byrd et al. 1995).
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 KarlLudwig 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 email at michael.hirsch@tuebingen.mpg.de.
References
 Anconelli, B., Bertero, M., Boccacci, P., Carbillet, M., & Lanteri, H. 2006, A&A, 448, 1217 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ayers, G. R., & Dainty, J. C. 1988, Opt. Lett., 13, 547 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Biraud, Y. 1969, A&A, 1, 124 [NASA ADS] [Google Scholar]
 Bottou, L. 1998, in Online Learning and Neural Networks, ed. D. Saad (Cambridge University Press), 9 [Google Scholar]
 Byrd, R., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM J. Scientific Comp., 16, 1190 [Google Scholar]
 DaubeWitherspoon, M. E., & Muehllehner, G. 1986, IEEE Tran. Medical Imaging, 5, 61 [Google Scholar]
 Farsiu, S., Robinson, D., Elad, M., & Milanfar, P. 2004, International Journal of Imaging Systems and Technology, 14, 47 [Google Scholar]
 Fried, D. L. 1978, J. Opt. Soc. Amer., 86 [Google Scholar]
 Harding, C. M., Johnston, R. A., & Lane, R. G. 1999, Appl. Opt., 38, 2161 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Harmeling, S., Hirsch, M., Sra, S., & Schölkopf, B. 2009, in Proceedings of the IEEE Conference on Computational Photography [Google Scholar]
 Harmeling, S., Sra, S., Hirsch, M., & Schölkopf, B. 2010, in Image Processing (ICIP), 17th IEEE Inter. Conf., 3313 [Google Scholar]
 Hormuth, F., Hippler, S., Brandner, W., Wagner, K., & Henning, T. 2008, SPIE Conf. Ser., 7014 [Google Scholar]
 Horn, R. A., & Johnson, C. R. 1991, Topics in Matrix Analysis (Cambridge: Cambridge University Press) [Google Scholar]
 Klöckner, A., Pinto, N., Lee, Y., et al. 2009, Performance Computing, 2009, 11 [Google Scholar]
 Knox, K. T., & Thompson, B. J. 1974, ApJ, 193, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Kundur, D., & Hatzinakos, D. 1996, IEEE Signal Processing Mag., 13, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Kushner, H. J., & Yin, G. G. 2003, Stochastic Approximation and Recursive Algorithms and Applications, 2nd edn., Applications of Mathematics (SpringerVerlag) [Google Scholar]
 Labeyrie, A. 1970, A&A, 6, 85 [NASA ADS] [Google Scholar]
 Law, N. M., Mackay, C. D., & Baldwin, J. E. 2006, A&A, 446, 739 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, B., Cao, Z., Sang, N., & Zhang, T. 2004, Electron. Lett., 1478 [Google Scholar]
 Löfdahl, M. G. 2002, Proc. SPIE, 4792, 146 [Google Scholar]
 Lohmann, A. W., & Weigelt, G., Wirnitzer, B. 1983, Appl. Opt., 22, 4028 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Mackay, C. D., Tubbs, R. N., Bell, R., et al. 2001, in SPIE Conf. Ser. 4306, ed. M. M. Blouke, J. Canosa, & N. Sampat, 289 [Google Scholar]
 Matson, C. L., Borelli, K., Jefferies, S., et al. 2008, Appl. Opt., 48 [Google Scholar]
 Molina, R., Núnez, J., Cortijo, F. J., & Mateos, J. 2001, IEEE Signal Proc. Mag. [Google Scholar]
 Oscoz, A., Rebolo, R., López, R., et al. 2008, SPIE Conf. Ser., 7014 [Google Scholar]
 Pantin, E., Starck, J. L., & Murtagh, F. 2007, Deconvolution and Blind Deconvolution in Astronomy: Theory and Applications (CRC Press) [Google Scholar]
 Park, S., Park, M., & Kang, M. 2003, IEEE Signal Processing Magazine [Google Scholar]
 Schulz, T. J. 1993, J. Opt. Soc. Amer., 10, 1064 [NASA ADS] [CrossRef] [Google Scholar]
 Sheppard, D. G., Hunt, B. R., & Marcellin, M. W. 1998, J. Opt. Soc. Amer. A, 15, 978 [NASA ADS] [CrossRef] [Google Scholar]
 Šroubek, F., Cristobál, G., & Flusser, J. 2007, IEEE Tran. Imag. Proc. [Google Scholar]
 Šroubek, F., Cristobál, G., & Flusser, J. 2008, J. Phys.: Conf. Ser., 124, 012048 (8pp) [NASA ADS] [CrossRef] [Google Scholar]
 Stelzer, C. 2009, Speckle1, http://www.tat.physik.unituebingen.de/stelzer/, version number 0.1.2 [Google Scholar]
 Stelzer, C., & Ruder, H. 2007, A&A, 475, 771 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Theusner, M. 2009, AviStack, http://www.avistack.de/index_en.html, version number 1.80 [Google Scholar]
 Zhang, J., Zhang, Q., & He, G. 2009, Appl. Opt., 48, 2350 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Zhulina, Y. V. 2006, Appl. Opt., 45, 7342 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
Appendix A: Implementation details
Although in Sects. 3 and 4 we only considered vectors, onedimensional convolutions, and vectornorms, all results naturally generalize to twodimensional images. However, efficiently implementing the resulting algorithms for twodimensional images requires some care and handling of technical details.
A.1. Convolution as matrixvector multiplication
We introduced f ∗ x as the convolution, which could be either circular or noncircular. Due to linearity and commutativity, we can also use matrixvector 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 zeropadding matrices which ensure that I_{x}x, I_{f}f, and I_{y}y have the same length. Different choices of the matrices lead to different margin condition of the convolution.
For twodimensional images and PSFs we have to consider twodimensional Fourier transforms, which can be written as left and rightmultiplications with W, and represented as a single matrixvector multiplication using Kronecker products and the vectorization operator vec(x), which stacks columns of the twodimensional image x into a onedimensional vector in lexicographical order; formally, (A.4)which follows from the identity (Horn & Johnson 1991) (A.5)The zeropadding operations for twodimensional images can be written in a similar way.
A.2. Resizing matrices
The resizing matrix can be implemented efficiently using sparse matrices^{5}. Resizing twodimensional images can also be implemented by left and rightmultiplications: 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 matrixvector product (A.6)
All Tables
True star magnitudes (note that stars A and B have variable magnitudes), true differences to star C, and estimated difference values estimated after deconvolution without and with saturation correction.
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 r_{o}), a short exposure image (bottom left) encapsulates information up to the diffractionlimited 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 diffractionlimited, short and long exposure imaging. While the long exposure MTF falls to nearly zero at r_{o}/λ, 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/r_{o} = 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, KnoxThompson 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 superresolution factors at a fixed SNR of 18.9 dB. The downsampling and superresolution 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 y_{300} (left), reconstruction x after 300 iterations (middle), estimated PSFs f_{300} 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 y_{t} (left), estimate of the corresponding PSF f_{t} and reconstruction x_{t} 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, KnoxThompson, our approach and our approach with two times superresolution. All image results are shown without any postprocessing. This figure is best viewed on screen, rather than in print. 

In the text 
Fig. 12 Orion Trapezium Cluster: example sequence of observed frames, y_{1}, ..., y_{10}. 

In the text 
Fig. 13 Orion Trapezium Cluster (from left to right): a) the first observed frame; b) x_{191} for basic algorithm; c) x_{191} for saturation corrected; and d) x_{191} for saturation corrected and four times superresolved. 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 superresolved multiframe blind deconvolution; (top) overview, (bottom) closeup. For better display the images have been automatically gamma corrected. 

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