\documentclass[11pt,twoside]{article}  % Leave intact
\usepackage{graphicx}
\usepackage{tabularx}
\usepackage{amsmath,amssymb,bm,xspace}
%\usepackage{theorem,cite,citesort}

\usepackage{adassconf}

\newcommand{\Otwonexmath}[1]{{\ensuremath{#1}\xspace}}
\newcommand{\Otwonebmath}[1]{\Otwonexmath{\bm{#1}}}

\newcommand{\Otwonebx}{\Otwonebmath{x}}
\newcommand{\Otwonebhx}{\Otwonebmath{\widehat{x}}}
\newcommand{\Otwonebn}{\Otwonebmath{n}}
\newcommand{\Otwoneby}{\Otwonebmath{y}}
\newcommand{\Otwonebz}{\Otwonebmath{z}}
\newcommand{\Otwonehf}{\widehat{f}}
\newcommand{\Otwonehy}{\widehat{y}}
\newcommand{\Otwonebhz}{\widehat{\Otwonebz}}
\newcommand{\Otwonebp}{\Otwonebmath{p}}
\newcommand{\Otwonebtheta}{\Otwonebmath{\theta}}
\newcommand{\Otwonebhtheta}{\widehat{\Otwonebtheta}}
\newcommand{\OtwonesP}{\mathcal{P}}
\newcommand{\Otwoneetal}{{\em et al}}

\begin{document}   % Leave intact

%-----------------------------------------------------------------------
%			    Paper ID Code
%-----------------------------------------------------------------------
% Enter the proper paper identification code.  The ID code for your
% paper is the session number associated with your presentation as
% published in the official conference proceedings.  You can           
% find this number locating your abstract in the printed proceedings
% that you received at the meeting or on-line at the conference web
% site; the ID code is the letter/number sequence proceeding the title 
% of your presentation. 
%
% This will not appear in your paper; however, it allows different
% papers in the proceedings to cross-reference each other.  Note that
% you should only have one \paperID, and it should not include a
% trailing period.
%
% EXAMPLE: \paperID{O4-1}
% EXAMPLE: \paperID{P7-7}
%

\paperID{O2-1}
%%%% ID=O2-1

%-----------------------------------------------------------------------
%		            Paper Title 
%-----------------------------------------------------------------------
% Enter the title of the paper.
%
% EXAMPLE: \title{A Breakthrough in Astronomical Software Development}
% 
% If your title is so long as to fill the page header when you print it,
% then please supply a short form as a \titlemark.
%
% EXAMPLE: 
%  \title{Rapid Development for Distributed Computing, with Implications
%         for the Virtual Observatory}
%  \titlemark{Rapid Development for Distributed Computing}
%

\title{Wavelet-Based Superresolution in Astronomy}
%\titlemark{ }

%-----------------------------------------------------------------------
%		          Authors of Paper
%-----------------------------------------------------------------------
% Enter the authors followed by their affiliations.  The \author and
% \affil commands may appear multiple times as necessary (see example
% below).  List each author by giving the first name or initials first
% followed by the last name.  Authors with the same affiliations
% should grouped together. 
%
% EXAMPLE: \author{Raymond Plante, Doug Roberts, 
%                  R.\ M.\ Crutcher\altaffilmark{1}}
%          \affil{National Center for Supercomputing Applications, 
%                 University of Illinois Urbana-Champaign, Urbana, IL
%                 61801}
%          \author{Tom Troland}
%          \affil{University of Kentucky}
%
%          \altaffiltext{1}{Astronomy Department, UIUC}
%
% In this example, the first three authors, "Plante", "Roberts", and
% "Crutcher" are affiliated with "NCSA".  "Crutcher" has an alternate 
% affiliation with the "Astronomy Department".  The fourth author,
% "Troland", is affiliated with "University of Kentucky"

\author{Rebecca Willett, Robert Nowak}
\affil{Rice University Department of Electrical and Computer Engineering,
Houston, TX and University of Wisconsin Department of Electrical and Computer
Engineering, Madison, WI} \author{Ian Jermyn, Josiane Zerubia}
\affil{INRIA, Project Ariana, Sophia-Antipolis, France}

%-----------------------------------------------------------------------
%			 Contact Information
%-----------------------------------------------------------------------
% This information will not appear in the paper but will be used by
% the editors in case you need to be contacted concerning your
% submission.  Enter your name as the contact along with your email
% address.
% 
% EXAMPLE:  \contact{Dennis Crabtree}
%           \email{crabtree@cfht.hawaii.edu}
%

\contact{Rebecca Willett}
\email{willett@rice.edu}

%-----------------------------------------------------------------------
%		      Author Index Specification
%-----------------------------------------------------------------------
% Specify how each author name should appear in the author index.  The 
% \paindex{ } should be used to indicate the primary author, and the
% \aindex for all other co-authors.  You MUST use the following
% syntax: 
%
% SYNTAX:  \aindex{Lastname, F. M.}
% 
% where F is the first initial and M is the second initial (if
% used).  This guarantees that authors that appear in multiple papers
% will appear only once in the author index.  
%
% EXAMPLE: \paindex{Crabtree, D.}
%          \aindex{Manset, N.}        
%          \aindex{Veillet, C.}        
%
% NOTE: this information is also used to build the author list that
% appears in the table of contents.  Authors will be listed in the order
% of the \paindex and \aindex commmands.
%

\paindex{Willett, R. M.}
\aindex{Jermyn, I.}
\aindex{Nowak, R. D.}
\aindex{Zerubia, J.}    

%-----------------------------------------------------------------------
%		      Author list for page header	
%-----------------------------------------------------------------------
% Please supply a list of author last names for the page header. in
% one of these formats:
%
% EXAMPLES:
% \authormark{Lastname}
% \authormark{Lastname1 \& Lastname2}
% \authormark{Lastname1, Lastname2, ... \& LastnameN}
% \authormark{Lastname et al.}
%
% Use the "et al." form in the case of seven or more authors, or if
% the preferred form is too long to fit in the header.

\authormark{Willett, Jermyn, Nowak, \& Zerubia}

%-----------------------------------------------------------------------
%			Subject Index keywords
%-----------------------------------------------------------------------
% Enter a comma separated list of up to 6 keywords describing your
% paper.  These will NOT be printed as part of your paper; however,
% they will be used to generate the subject index for the proceedings.
% There is no standard list; however, you can consult the indices
% for past proceedings (http://adass.org/adass/proceedings/).
%
% EXAMPLE:  \keywords{visualization, astronomy: radio, parallel
%                     computing, AIPS++, Galactic Center}
%
% In this example, the author noticed that "radio astronomy" appeared
% in the ADASS VII Index as "astronomy" being the major keyword and
% "radio" as the minor keyword.  The colon is used to introduce another
% level into the index.

\keywords{multi: resolution, reconstruction, registration, inverse problems}

%-----------------------------------------------------------------------
%			       Abstract
%-----------------------------------------------------------------------
% Type abstract in the space below.  Consult the User Guide and Latex
% Information file for a list of supported macros (e.g. for typesetting 
% special symbols). Do not leave a blank line between \begin{abstract} 
% and the start of your text.

\begin{abstract}        
High-resolution astronomical images can be reconstructed from several blurred
and noisy low-resolution images using a computational process known as
superresolution reconstruction. Superresolution reconstruction is closely
related to image deconvolution, except that the low-resolution images are not
registered and their relative translations and rotations must be estimated in
the process. The novelty of our approach to the superresolution problem is the
use of wavelets and related multiresolution methods within an
expectation-maximization reconstruction process to improve the accuracy and
visual quality of the reconstructed image. Simulations demonstrate the
effectiveness of the proposed method, including its ability to distinguish
between tightly grouped stars with a small set of observations.
\end{abstract}

\section{Introduction}
The physical resolution of astronomical imaging devices such as space telescopes
is limited by system parameters such as lens aperture and CCD array properties
and by physical effects such as the atmosphere and the
interstellar/intergalactic medium. However, such systems typically do not take a
single snapshot of a celestial body, but rather collect a series of images. Due
to mechanical vibrations of the instrument and movement of the satellite
relative to the body being images, the images collected are all slightly
different observations of the same scene. Superresolution image reconstruction
refers to the process of reconstructing a new image with a higher resolution
using this collection of low resolution, shifted, rotated, and often noisy
observations. This allows users to see image detail and structures which are
difficult if not impossible to detect in the raw data.

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 log-likelihood 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 non-quadratic 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. Expectation-Maximization 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 well-suited 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 ({\em 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 low-complexity,
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 E-step
based on the fast Fourier transform (FFT) and a DWT-based M-step, resulting in
an efficient iterative process requiring $O(N \log N)$ operations per iteration,
where $N$ is the number of pixels in the superresolution image.


\section{Problem Formulation}
In the proposed method, each observation, $y_{k}$, is modeled as a shifted,
rotated, blurred, downsampled, and noisy version of the superresolution $\Otwonebx$.
The shift and rotation is caused by the movement of the instrument, and the blur
is caused by the point spread function (PSF) of the instrument optics and the
integration done by the CCD array. The downsampling (subsampling) operator
models the change in resolution between the observations and the desired
superresolution image. If the noise can be modeled as additive white Gaussian
noise, then we have the observation model
$$y_{k} = DBS_{k}\Otwonebx + n_{k}, \;\; k = 1, \ldots, n$$
where $D$ is the downsampling operator, $B$ is the blurring operator, $S_{k}$ is
the shift and rotation operator for the $k^{th}$ observation, and $n_{k}$ is
additive white Gaussian noise with variance $\sigma^{2}$. By collecting the
series of observations into one array, $\Otwoneby$, the noise observations into
another array, $\Otwonebn$, and letting $H$ be a block diagonal matrix composed of the
$n$ matrixes $DBS_{k}$ for $k = 1, \ldots, n$, then we have the model
\begin{equation}
\Otwoneby = H\Otwonebx+\Otwonebn.
\label{O2-1:eq:inverse}
\end{equation}
Similarly, if the noise is better modeled a Poisson, then we have the model
\begin{equation}
\Otwoneby \sim \mbox{Poisson}(H\Otwonebx).
\label{O2-1:eq:poisson}
\end{equation}

From the formulation above, it is clear that superresolution image
reconstruction is a type of inverse problem in which the operator to be
inverted, $H$, 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 low-resolution observations to one another. Using these
estimates, we reconstruct an initial superresolution image estimate $\Otwonebhx$. This
estimate is used in the third step, where we re-estimate the shift and rotation
parameters by registering each of the low resolution observations to the initial
superresolution estimate. Finally, we use a wavelet-based EM algorithm to solve
for $\Otwonebhx$ using the registration parameter estimates. We begin by describing a
wavelet-based 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.

\section{Registration of the Observations}

The first step in the proposed method is to register the observed low-resolution
images to one another using a Taylor series expansion. This was proposed by
Irani and Peleg, 1991. In particular, let $f_{1}$ and $f_{2}$ be the continuous
images underlying the sampled images $y_{1}$ and $y_{2}$, respectively. If
$f_{2}$ is equal to a shifted, rotated version of $f_{1}$, then we have the
relation
$$f_{2}(t_{x},t_{y}) = f_{1}(t_{x}\cos r - t_{y} \sin r + s_{x}, t_{y}\cos r + t_{x}\sin r + s_{y}).$$ 
where $(s_{x}, s_{y})$ is the shift and $r$ is the rotation.
A first order Taylor series approximation of $f_{2}$ is then
$$\Otwonehf_{2}(t_{x},t_{y}) = f_{1}(t_{x},t_{y}) + \left( s_{x} - t_{x}r -
t_{y}r^{2}/2 \right) \frac{\partial f_{1}}{\partial t_{x}} + \left( s_{y} -
t_{y}r - t_{x}r^{2}/2 \right) \frac{\partial f_{1}}{\partial t_{y}}. $$ Let
$\Otwonehy_{2}$ be a sampled version of $\Otwonehf_{2}$; then $y_{1}$ and $y_{2}$ can be
registered by finding the $s_{x}$, $s_{y}$, and $r$ which minimize
$\|y_{1}-\Otwonehy_{2}\|^{2}_{2}$, where $\|x\|^{2}_{2} = \sum_{i} x_{i}^{2}$. This
minimization is calculated with an iterative procedure which ensures that the
motion being estimated at each iteration is small enough for the Taylor series
approximation to be accurate; see (Irani \& Peleg, 1991) for details. This
method was applied to the earth image data in Figure \ref{O2-1:fig:EarthResults}. The
sixteen true shifts and rotations are displayed in black in Figure
\ref{O2-1:fig:registration}, and the results of this registration method are
displayed in hollow circles and bars.

\begin{figure}[ht]
\centerline{
\begin{tabular}{ccc}
\includegraphics[width=2.2in]{O2-1_shiftEst.eps} &&
\includegraphics[width=2.2in]{O2-1_rotationEst.eps} \\
(a) & & (b)
\end{tabular} 
}
\caption{Image registration results. (a) True shifts (black), initial shift
estimates (hollow), and final shift estimates (gray). (b) True rotations
(black), initial rotation estimates (hollow), and final rotation estimates
(gray).} 
\label{O2-1:fig:registration}
\end{figure}

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 \ref{O2-1:sec:em}. 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, $y_{2}$, to another low
resolution estimate, $y_{1}$, we instead register it to $DBS_{1}\Otwonebhx$. The
results of this refinement are displayed in gray in Figure
\ref{O2-1:fig:registration}. 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.

\section{Multiscale Expectation-Maximization}
\label{O2-1:sec:em}
Maximization is facilitated within the EM framework through the introduction of
a particular ``unobservable'' or ``missing'' data space. The key idea in the EM
algorithm is that the indirect (inverse) problem can be broken into two
subproblems; one which involves computing the expectation of unobservable data
(as though no blurring or downsampling took place) and one which entails
estimating the underlying image from this expectation. By carefully defining the
unobservable data for the superresolution problem, we derive an EM algorithm
which consists of linear filtering in the E-step and image denoising in the
M-step.

The Gaussian observation model in (\ref{O2-1:eq:inverse}) can be written with respect
to the DWT coefficients $\Otwonebtheta$, where $\Otwonebx = W\Otwonebtheta$ and $W$ denotes the
inverse DWT operator (Mallat, 1998): $$\Otwoneby = H W \Otwonebtheta + \Otwonebn.$$ Clearly, if we
had $\Otwoneby = W\Otwonebtheta + \Otwonebn$ ({\em i.e.} if no subsampling or blurring had
occurred), we would have a pure image denoising problem with white noise, for
which wavelet-based denoising techniques are very fast and nearly optimal
(Mallat, 1998).
Next note that the noise in the observation model can be decomposed into two
different Gaussian noises (one of which is non-white): $$\Otwonebn = \alpha H \Otwonebn_{1}
+ \Otwonebn_{2}$$
where $\alpha$ is a positive parameter, and $\Otwonebn_{1}$ and $\Otwonebn_{2}$ are
independent zero-mean Gaussian noises with covariances $\Sigma_{1} = I$ and
$\Sigma_{2} = \sigma^{2}I-\alpha^{2}HH^{T}$, respectively. Using $\Otwonebn_{1}$ and
$\Otwonebn_{2}$, we can rewrite the Gaussian observation model as $$\Otwoneby =
H\underbrace{(W \Otwonebtheta + \alpha \Otwonebn_{1})}_{\Otwonebz} +\Otwonebn_{2}.$$
This observation is the key to our approach since it suggests treating $\Otwonebz$ as
missing data and using the EM algorithm to estimate $\Otwonebtheta$.

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
$\Otwonebx$. Those photons emitted from location $m$ are detected (in the detection
space) at position $n$ with transition probability $H_{m,n}$, where $H$ is the
superresolution operator from (\ref{O2-1:eq:poisson}). In such cases, the measured
data are distributed according to 
\begin{equation}
\Otwoneby_{n} \sim \mbox{Poisson}\left( \sum_{m} H_{m,n} x_{m}\right).
\label{O2-1:eq:Poiss}
\end{equation}
In this formulation of the EM algorithm, the missing data is defined as $\Otwonebz =
\{z_{m,n}\}$, where $z_{m,n}$ denotes the number of photons emitted from $m$ and
detected at $n$. The complete data are Poisson distributed according to
$$z_{m,n} \sim \mbox{Poisson}(H_{m,n} \Otwonebx_{m} ).$$ Hence the observed data $\Otwoneby$
in (\ref{O2-1:eq:Poiss}) are given by $\Otwoneby_{n} = \sum_{m} z_{m,n}$. Additionally,
were we able to observe $\Otwonebz = \{z_{m,n}\}$, the direct emission data for each
location $m$ is given by sums of the form $\sum_{n} z_{m,n} \sim
\mbox{Poisson}(\Otwonebx_{m})$. Therefore, if $\Otwonebz$ were known, we could avoid the
inverse problem altogether and simply deal with the issue of estimating a
Poisson intensity given direct observations.

From these formulations of the problem for Gaussian and Poisson data, the EM
algorithm produces a sequence of estimates $\{ \Otwonebx^{(t)}, t = 0, 1, 2, \ldots
\}$ by alternately applying two steps: \begin{description}
\item[E-Step:] Updates the estimate of the missing data using the relation:
$$\Otwonebhz^{(t)} = E \left[\Otwonebz | \Otwoneby, \Otwonebhtheta^{(t)}\right].$$
In the Gaussian case this can be reduced to a Landweber iteration (Landweber, 1951):
$$\Otwonebhz^{(t)} = \Otwonebhx^{(t)} + \frac{\alpha^{2}}{\sigma^{2}}H^{T}\left(\Otwoneby-H\Otwonebhx^{(t)}\right).$$
Here, computing $\Otwonebhz^{(t)}$ simply involves applications of the operator $H$.
Recall that $H$ consists of shifting, blurring, and rotation (which can be
computed rapidly with the 2D FFT) and downsampling (which can be computed
rapidly in the spatial domain). Thus the complexity of each E-Step is 
$O(N \log N)$.

In the Poisson case, this step is reduced to 
$$z^{(t)}_{n,m} = 
   \frac{\Otwoneby_{n}\Otwonebhx_{m}^{(t)}H_{n,m}}{\sum_{\ell}\Otwonebhx_{\ell}^{(t)}H_{n,\ell}},
$$
which can also be computed in $O(N \log N)$ operations.

\item[M-Step:] Updates the estimate of the superresolution image $\Otwonebx$. In the
Gaussian case, this constitutes updating the wavelet coefficient vector
$\Otwonebtheta$ according to
$$\Otwonebhtheta^{(t+1)} = \arg \min_{\Otwonebtheta} \left\{ \frac{\| W \Otwonebtheta -\Otwonebhz^{(t)} 
   \|^{2}_{2}}{2 \alpha^{2}} + \mbox{pen}(\Otwonebtheta) \right\}$$
and setting $\Otwonebhx^{(t+1)} = W\Otwonebhtheta^{(t+1)}$.
This optimization can be preformed using any wavelet-based denoising procedure.
For example, under an {\em i.i.d.} Laplacian prior, $\mbox{pen}(\Otwonebtheta) = -\log
p(\Otwonebtheta) \propto \tau \|\Otwonebtheta\|_{1}$ (where $\|\Otwonebtheta\|_{1} = \sum_{i}
|\Otwonebtheta_{i}|$ denotes the $l_{1}$ norm), $\Otwonebhtheta^{(t+1)}$ is obtained by
applying a {\em soft-threshold} function to the wavelet coefficients of
$\Otwonebhz^{(t)}$. For the simulations presented in this paper, we applied a similar
denoising method described in (Figueiredo \& Nowak, 2001), which requires $O(N)$
operations.

In the Poisson case, the M-Step is equivalent to photon-limited 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 Haar-based maximum
penalized likelihood estimation method we developed in Willett \& Nowak, 2003;
the MPLE function employed here is

\begin{equation}
L(\Otwonebx) \ \equiv \ \log p(\Otwoneby\,|\,\sum_{m}z^{(t)}_{m,n})\, - \,  \frac{1}{2} \, \mbox{pen}(\sum_{m}z^{(t)}_{m,n})\log n,
\label{O2-1:rple} 
\end{equation}

where $p(\Otwoneby\,|\,\Otwonebx)$ is the Poisson likelihood of observing photon counts
$\Otwoneby$ given intensities $\sum_{m}z^{(t)}_{m,n}$ and
$\mbox{pen}(\sum_{m}z^{(t)}_{m,n})$ is the number of non-zero coefficients in
the Haar-based multiscale representation of $\sum_{m}z^{(t)}_{m,n}$ and $n$ is
the total number of photons in the observation vector $\Otwoneby$. 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 bottom-up tree pruning algorithm which requires
$O(N \log N)$ operations.

\end{description}


In both the Gaussian and Poisson noise cases, the proposed method has the
advantage that the M-step is a denoising procedure, and the multiscale methods
employed here are both near-minimax 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 $\widehat{\Otwonebx}^{(t)}$, 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.
 
\section{Simulation Results}

We reconstructed two superresolution images to demonstrate the practical
effectiveness of the proposed algorithms. A sample low-resolution satellite
image of earth is displayed in Figure \ref{O2-1:fig:EarthResults}(a). Sixteen such
$64 \times 64$ images were generated using an original $256 \times 256$ image
(Figure \ref{O2-1:fig:EarthResults}(c), with values ranging from 0 to 255), which was
shifted and rotated by a random, subpixel amount, distorted by a $ 4 \times 4$
uniform blur, and contaminated with additive white Gaussian noise with variance
$\sigma^{2} = 1/2$. The superresolution image in Figure
\ref{O2-1:fig:EarthResults}(b) was reconstructed using the wavelet-based EM algorithm
described above for $\alpha^{2} = \sigma^{2} = 1/2$. In the simulation, the
estimate is initialized with a least-squares 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.

\begin{figure}[ht]
\centerline{
\begin{tabular}{cc}
\multicolumn{2}{c}{
\includegraphics[width=2.2in]{O2-1_EarthObs.ps}
}\\
\multicolumn{2}{c}{(a)}\\
\includegraphics[width=2.2in]{O2-1_EarthABE.ps} &
\includegraphics[width=2.2in]{O2-1_EarthTrue.ps} \\
(b)  & (c)
\end{tabular} 
}
\caption{Superresolution results for earth image. (a) One of $16$ observation
images ($64 \times 64$), contaminated with Gaussian noise, $\sigma^{2} = 0.5$
and a $4 \times 4$ uniform blur. (b) $256 \times 256$ result. (c) True high
resolution image.} 
\label{O2-1:fig:EarthResults}
\end{figure}

In our second simulation, we studied the effect of the proposed method on an
image of stars. The original image is displayed in Figure
\ref{O2-1:fig:StarsResults}(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 \ref{O2-1:fig:StarsResults}(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
\ref{O2-1:fig:StarsResults}(b). The superresolution image in Figure
\ref{O2-1:fig:StarsResults}(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 \ref{O2-1:fig:StarsResults}(e) displayed the output of our multiscale
method for Poisson noise. In the case, the mean photon count was $527$. As in
the Gaussian case, the improvement in resolution is significant.

\begin{figure}[ht]
\centerline{
\begin{tabular}{cc}
\multicolumn{2}{c}{
\includegraphics[width=2.2in]{O2-1_StarsObs.ps}
}\\
\multicolumn{2}{c}{(a)}\\
\includegraphics[width=2.2in]{O2-1_StarsABE.ps} &
\includegraphics[width=2.2in]{O2-1_StarsTrue.ps} \\
(b)  & (c)\\
\includegraphics[width=2.2in]{O2-1_StarsABE02.ps} &
\includegraphics[width=2.2in]{O2-1_StarsHaarP.ps} \\
(d)  & (e)
\end{tabular} 
}
\caption{Superresolution results for stars image. (a) One of $16$ observation
images ($64 \times 64$), contaminated with Gaussian noise, $\sigma^{2} = 0.5$
and a $4 \times 4$ uniform blur. (b) $256 \times 256$ result. (c) True high
resolution image. (d) $256 \times 256$ result obtained using only two $64 \times
64$ observation images. (e) $256 \times 256$ result obtained using $16$ $64
\times 64$ observation images contaminated with Poisson noise; mean photon count
per observation image = 527.} 
\label{O2-1:fig:StarsResults}
\end{figure}

\section{Conclusions and Future Work}

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 wavelet-based superresolution has the potential
to make a significant impact on astronomical imaging.





%-----------------------------------------------------------------------
%			      References
%-----------------------------------------------------------------------
% List your references below within the reference environment
% (i.e. between the \begin{references} and \end{references} tags).
% Each new reference should begin with a \reference command which sets
% up the proper indentation.  Observe the following order when listing
% bibliographical information for each reference:  author name(s),
% publication year, journal name, volume, and page number for
% articles.  Note that many journal names are available as macros; see
% the User Guide listing "macro-ized" journals.   
%
% EXAMPLE:  \reference Hagiwara, K., \& Zeppenfeld, D.\  1986, 
%                Nucl.Phys., 274, 1
%           \reference H\'enon, M.\  1961, Ann.d'Ap., 24, 369
%           \reference King, I.\ R.\  1966, \aj, 71, 276
%           \reference King, I.\ R.\  1975, in Dynamics of Stellar 
%                Systems, ed.\ A.\ Hayli (Dordrecht: Reidel), 99
%           \reference Tody, D.\  1998, \adassvii, 146
%           \reference Zacharias, N.\ \& Zacharias, M.\ 2003,
%                \adassxii, \paperref{P7.6}
% 
% Note the following tricks used in the example above:
%
%   o  \& is used to format an ampersand symbol (&).
%   o  \'e puts an accent agu over the letter e.  See the User Guide
%      and the sample files for details on formatting special
%      characters.  
%   o  "\ " after a period prevents LaTeX from interpreting the period 
%      as an end of a sentence.
%   o  \aj is a macro that expands to "Astron. J."  See the User Guide
%      for a full list of journal macros
%   o  \adassvii is a macro that expands to the full title, editor,
%      and publishing information for the ADASS VII conference
%      proceedings.  Such macros are defined for ADASS conferences I
%      through XI.
%   o  When referencing a paper in the current volume, use the
%      \adassxii and \paperref macros.  The argument to \paperref is
%      the paper ID code for the paper you are referencing.  See the 
%      note in the "Paper ID Code" section above for details on how to 
%      determine the paper ID code for the paper you reference.  
%
\begin{references}

\reference Donoho, D., 1999,
Ann. Statist., 27, 859

\reference Figueiredo, M. \& Nowak, R., 2002,
IEEE Transactions on Image Processing, 29, 1033

\reference  Figueiredo, M. \& Nowak, R., 2001,
IEEE Transactions on       Image Processing, 10, 1322

\reference Green, P.\ J., 1990,
 IEEE Trans. Med. Imaging, 9, 84
 % “Bayesian reconstruction from emission tomography data using a modified EM algorithm,”
 
\reference Hardie, R.\ C., Barnard, K.\ J., \& Armstrong, E.\ E., 1997,
IEEE Transactions on Image Processing, 6, 1621

\reference Hebert, T.\ \& Leahy, R., 1989,
IEEE Trans. Med. Imaging, 8, 194
%“A generalized EM algorithm for 3-d Bayesian reconstruction from Poisson data using Gibbs priors,”

\reference Irani, M. \& Peleg, S.\ 1991,
Computer Vis. Graph. Image Process.: Graph. Models Image Process, 53, 231
% "Improving resolution by image registration,"

\reference Kolaczyk, E., 1999,
J. Amer. Statist. Assoc., 94, 920

\reference Kolaczyk, E.\ \& Nowak, R., 2003,
to appear in the Annals of Statistics, ``Multiscale Likelihood Analysis and Complexity Penalized Estimation'',
\\ \verb+http://www.ece.wisc.edu/~nowak/pubs.html+

\reference Krim, H., \& Schick, I. C., 1999,
IEEE Trans. on Information Theory, 45
%“Minmax description length for signal denoising and optimal representation,”

\reference Landweber, L., 1951,
American Journal of Mathematics, 73, 615

 \reference Liu, J. \& Moulin, P., 2001,
 IEEE Transactions on Image Processing, 10, 841
%“Complexity-regularized image denoising,”

\reference Mallat, S., 1998,
A Wavelet Tour of Signal Processing

\reference Nowak, R., \& Kolaczyk, E., 2000,
IEEE Transactions on Information Theory, 46, 1811

\reference Nguyen, N., Milanfar, P., \& Golub, G., 2001,
IEEE Transactions on Image Processing, 10, 573

\reference Saito, N, 1994,
Wavelets in Geophysics, ed.\ Foufoula-Georgiou and Kumar
% “Simultaneous noise suppression and signal compression using a library of orthonormal bases and the minimum description length criterion,”


\reference Salakhutdinov, R.\ \& Roweis, S., 2003,
International Conference on Machine Learning, 664

\reference Schultz, R.\ R.\ \& Stevenson, R.\ L., 1996,
IEEE Transactions on Image Processing, 5, 996

\reference Starck, J., Murtagh, F., \& Bijaoui, A., 1998,
Image Processing and Data Analysis: The Multiscale Approach

\reference Timmermann, K. \& Nowak, R, 1999,
IEEE Transactions on Information Theory, 45, 846

\reference Willett, R.\ \& Nowak, R., 2003,
IEEE Transactions of Medical Imaging, 22, 332



\end{references}

% Do not place any material after the references section

\end{document}  % Leave intact
