Superresolution is a useful technique in a variety of applications (Schultz & Stevenson, 1996; Hardie, Barnard, & Armstrong, 1997), and recently, researchers have begun to investigate the use of wavelets for superresolution image reconstruction (Nguyen, Milanfar, & Golub, 2001). We present a new method for superresolution image reconstruction based on the wavelet transform in the presence of Gaussian noise and on an analogous multiscale approach in the presence of Poisson noise. To construct the superresolution image, we use an approach based on maximum penalized likelihood estimation. The reconstructed image is the argument (in our case the superresolution image) that maximizes the sum of a loglikelihood function and a penalizing function. The penalizing function can be specified by an ad hoc smoothness measure, a Bayesian prior distribution for the image (Hebert & Leahy, 1989; Green, 1990), or a complexity measure (Liu & Moulin, 2001). Smoothness measures include simple quadratic functions that measure the similarity between the intensity values of neighboring pixels, as well as nonquadratic measures that better preserve edges. Similar penalty functions result from Markov Random Field (MRF) priors. Complexity measures are usually associated with an expansion of the intensity image with respect to a set of basis functions (e.g. Fourier or wavelet) and count the number of terms retained in a truncated or pruned series (Saito, 1994; Krim & Schick, 1999); the more terms (basis functions) used to represent the image, the higher the complexity measure. Many algorithms (e.g. ExpectationMaximization algorithms, the Richardson Lucy algorithm, or close relatives) have been developed to compute MPLEs under various observation models and penalization schemes.
Wavelets and multiresolution analysis are especially wellsuited for astronomical image processing because they are adept at providing accurate, sparse representations of images consisting of smooth regions with isolated abrupt changes or singularities ( e.g. stars against a dark sky). Many investigators have considered the use of wavelet representations for image denoising, deblurring, and image reconstruction; for examples, see Mallat, 1998, and Starck, Murtagh, & Bijaoui, 1998. The proposed approach uses an EM algorithm for superresolution image reconstruction based on a penalized likelihood formulated in the wavelet domain. Regularization is achieved by promoting a reconstruction with lowcomplexity, expressed in terms of the wavelet coefficients, taking advantage of the well known sparsity of wavelet representations.
The EM algorithm proposed here extends the work of Nowak & Kolaczyk, 2000, and Figueiredo & Nowak, 2002, which addressed image deconvolution with a method that combines the efficient image representation offered by the discrete wavelet transform (DWT) with the diagonalization of the convolution operator obtained in the Fourier domain. The algorithm alternates between an Estep based on the fast Fourier transform (FFT) and a DWTbased Mstep, resulting in an efficient iterative process requiring operations per iteration, where is the number of pixels in the superresolution image.
From the formulation above, it is clear that superresolution image reconstruction is a type of inverse problem in which the operator to be inverted, , is partially unknown due to the unknown shifts and rotations of the observations. The first step of our approach is to estimate these parameters by registering the lowresolution observations to one another. Using these estimates, we reconstruct an initial superresolution image estimate . This estimate is used in the third step, where we reestimate the shift and rotation parameters by registering each of the low resolution observations to the initial superresolution estimate. Finally, we use a waveletbased EM algorithm to solve for using the registration parameter estimates. We begin by describing a waveletbased method for the Gaussian noise model, and follow that by a discussion of a multiscale technique for Poisson data. Each of these steps is detailed below.
The first step in the proposed method is to register the observed lowresolution
images to one another using a Taylor series expansion. This was proposed by
Irani and Peleg, 1991. In particular, let and be the continuous
images underlying the sampled images and , respectively. If
is equal to a shifted, rotated version of , then we have the
relation

After the registration parameters have been initially estimated using the above method, we use these estimates to calculate an initial superresolution image as described in Section 4.. This initial image estimate is then used to refine the registration parameter estimates. The method is the same as above, but instead of registering a low resolution estimate, , to another low resolution estimate, , we instead register it to . The results of this refinement are displayed in gray in Figure 1. From these plots, it is clear that the Taylor series based approach can produce highly accurate results. However, in low SNR scenarios, where confidence in registration parameter estimates may be low, the estimates can be further refined at each iteration of the proposed EM algorithm, as discussed in the following sections.
Note that the motion model considered here encompasses only shift and rotation movement. When recording still or relatively still objects distant from the imaging device, this model is sufficient. More sophisticated models are an area of open research.
The Gaussian observation model in (1) can be written with respect
to the DWT coefficients
, where
and denotes the
inverse DWT operator (Mallat, 1998):
An analogous formulation is possible for the Poisson noise model. In this case,
photon projections can be described statistically as follows. Photons are
emitted (from the emission space) according to a high resolution intensity
. Those photons emitted from location are detected (in the detection
space) at position with transition probability , where is the
superresolution operator from (2). In such cases, the measured
data are distributed according to
From these formulations of the problem for Gaussian and Poisson data, the EM algorithm produces a sequence of estimates by alternately applying two steps:
In the Poisson case, this step is reduced to
In the Poisson case, the MStep is equivalent to photonlimited image denoising. In practice, this commonly consists of performing maximum likelihood estimation, but these estimates are known to diverge from the true image after several iterations. The denoising can also be accomplished using the Haarbased maximum penalized likelihood estimation method we developed in Willett & Nowak, 2003; the MPLE function employed here is
where is the Poisson likelihood of observing photon counts given intensities and is the number of nonzero coefficients in the Haarbased multiscale representation of and is the total number of photons in the observation vector . In maximizing this function, the resulting reconstruction will be one that has a relatively high likelihood value as well as a relatively low complexity Haar representation. This denoising procedure is a bottomup tree pruning algorithm which requires operations.
In both the Gaussian and Poisson noise cases, the proposed method has the advantage that the Mstep is a denoising procedure, and the multiscale methods employed here are both nearminimax optimal.
In low SNR cases, where confidence in the initial estimates of the shift and rotation parameters may be low, the proposed algorithm can be modified by simply inserting an additional step in which the parameter estimates are updated based on the current estimate of the superresolution image , as described in the previous section. Similarly, if the blurring operator is only partially known, its parameters can also be updated at each iteration of the proposed algorithm. The resulting scheme is not a standard EM algorithm, but it is guaranteed not to decrease the penalized likelihood function.
We reconstructed two superresolution images to demonstrate the practical effectiveness of the proposed algorithms. A sample lowresolution satellite image of earth is displayed in Figure 2(a). Sixteen such images were generated using an original image (Figure 2(c), with values ranging from 0 to 255), which was shifted and rotated by a random, subpixel amount, distorted by a uniform blur, and contaminated with additive white Gaussian noise with variance . The superresolution image in Figure 2(b) was reconstructed using the waveletbased EM algorithm described above for . In the simulation, the estimate is initialized with a leastsquares superresolution estimate of relatively poor quality. While not presented here, experiments have shown that the proposed approach is competitive with the state of the art in superresolution image reconstruction.

In our second simulation, we studied the effect of the proposed method on an image of stars. The original image is displayed in Figure 3(c). The data for this simulation was generated using the same procedure as for the Earth image. Note from the sample observation image in Figure 3(a) that several stars are indistinguishable prior to superresolution image reconstruction, but that after the application of the proposed method these stars are clearly visible in Figure 3(b). The superresolution image in Figure 3(d) was generated using the same procedure but with only two observation images. Clearly the quality of this result is less than the quality achievable with more observations, but significantly more stars are distinguishable in the output than in any one of the observations. This implies that the proposed method may be useful even for small sets of observations. Finally, Figure 3(e) displayed the output of our multiscale method for Poisson noise. In the case, the mean photon count was . As in the Gaussian case, the improvement in resolution is significant.

In this paper we present a multiscale approach for superresolution image reconstruction in the presence of either Gaussian or Poisson noise. This approach uses several shifted, rotated, blurred, noisy observations to construct a new image with higher resolution than any of the observations. Because of the multiscale complexity regularization used in the the EM reconstruction algorithm, our approach is robust to noise. This is demonstrated for several practical examples. Notably, we demonstrated that stars which are irresolvable in any individual observation can be clearly distinguished after superresolution image reconstruction, even when using only two observation images.
Future work includes speeding the convergence of the proposed method using a new technique described in (Salakhutdinov & Roweis, 2003) and a more Bayesian approach in which the shift and rotation parameters are integrated out of the problem formulation. In addition, there are several questions still to be addressed, including: How does the blur radius impact reconstruction quality? How much can resolution be recovered at what accuracy? How can the proposed multiscale methods be extended to optimally process data collected by photon detector cells with spatially varying sensitivities? Despite these open areas of research, however, we feel that waveletbased superresolution has the potential to make a significant impact on astronomical imaging.
Donoho, D., 1999, Ann. Statist., 27, 859
Figueiredo, M. & Nowak, R., 2002, IEEE Transactions on Image Processing, 29, 1033
Figueiredo, M. & Nowak, R., 2001, IEEE Transactions on Image Processing, 10, 1322
Green, P. J., 1990, IEEE Trans. Med. Imaging, 9, 84
Hardie, R. C., Barnard, K. J., & Armstrong, E. E., 1997, IEEE Transactions on Image Processing, 6, 1621
Hebert, T. & Leahy, R., 1989, IEEE Trans. Med. Imaging, 8, 194
Irani, M. & Peleg, S. 1991, Computer Vis. Graph. Image Process.: Graph. Models Image Process, 53, 231
Kolaczyk, E., 1999, J. Amer. Statist. Assoc., 94, 920
Kolaczyk, E. & Nowak, R., 2003,
to appear in the Annals of Statistics, ``Multiscale Likelihood Analysis and Complexity Penalized Estimation'',
http://www.ece.wisc.edu/~nowak/pubs.html
Krim, H., & Schick, I. C., 1999, IEEE Trans. on Information Theory, 45
Landweber, L., 1951, American Journal of Mathematics, 73, 615
Liu, J. & Moulin, P., 2001, IEEE Transactions on Image Processing, 10, 841
Mallat, S., 1998, A Wavelet Tour of Signal Processing
Nowak, R., & Kolaczyk, E., 2000, IEEE Transactions on Information Theory, 46, 1811
Nguyen, N., Milanfar, P., & Golub, G., 2001, IEEE Transactions on Image Processing, 10, 573
Saito, N, 1994, Wavelets in Geophysics, ed. FoufoulaGeorgiou and Kumar
Salakhutdinov, R. & Roweis, S., 2003, International Conference on Machine Learning, 664
Schultz, R. R. & Stevenson, R. L., 1996, IEEE Transactions on Image Processing, 5, 996
Starck, J., Murtagh, F., & Bijaoui, A., 1998, Image Processing and Data Analysis: The Multiscale Approach
Timmermann, K. & Nowak, R, 1999, IEEE Transactions on Information Theory, 45, 846
Willett, R. & Nowak, R., 2003, IEEE Transactions of Medical Imaging, 22, 332