Abstract: Many astronomical studies rely upon the accurate reconstruction of spatially distributed phenomena from photon-limited data These measurements are inherently “noisy” due to low photon counts In addition, the behavior of the underlying photon intensity functions can be very rich and complex, and consequently difficult to model a priori Nonparametric multiscale reconstruction methods overcome these challenges and facilitate characterization of fundamental performance limits In this paper, we review several multiscale approaches to photon-limited image reconstruction, including wavelets combined with variance stabilizing transforms, corrected Haar wavelet transforms, multiplicative multiscale innovations, platelets, and the a trous wavelet transform We discuss the performance of these methods in simulation studies, and detail statistical analyses of their performances 1 Photon-limited astronomy Many imaging modalities involve the detection of light or higher energy photons, and often the random nature of photon emission and detection is the dominant source of noise in imaging systems Such cases are referred to as photon-limited imaging applications, since the relatively small number of detected photons is the factor limiting the signal-to-noise ratio In many cases, the intensity underlying the photonlimited observations may be distorted by the point spread function or other physical properties of the imaging system Using these inherently noisy and distorted observations to perform quantitative inference on the underlying astrophysical phenomenon is a challenging problem affecting many researchers in the statistical and astronomical communities The data collected by these imaging systems are usually assumed to obey a spatial Poisson distribution involving a two-dimensional intensity image that describes the probability of photon emissions at different locations in space The mean and variance of a Poisson process are equal to the intensity The intensity/mean is the “signal” of interest and the variability of the data about the mean can be interpreted as “noise” Thus, as the intensity varies spatially as a function of astrophysical structure and function, so does the signal-to-noise ratio In this sense it could be said that the noise in photon-limited imaging is signal-dependent Many investigators have considered the use of wavelet representations for image denoising, deblurring, and other forms of image reconstruction because of the theoretical near-optimality and practical efficacy of wavelets in a variety of image processing contexts; for examples, see (Mallat 1998; Starck et al 1998; Aldroubi and Unser 1996) The procedure for classic wavelet denoising via hard thresholding is the following: compute the wavelet transform of the noisy image, set wavelet coefficients with magnitude less than some threshold to zero, and compute the inverse wavelet transform Wavelet denoising via soft thresholding is very similar, except each coefficient is either set to zero or shrunk depending upon its magnitude The basic idea behind these methods is that wavelet bases form a very parsimonious representation of many images of interest, so that the bulk of the image’s energy in concentrated in just a few wavelet coefficients, which as a result have very high magnitudes Most noise, however, does not share this property, and has its energy distributed relatively evenly among all the wavelet coefficients Thus, by