Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.stsci.edu/stsci/meetings/irw/proceedings/lucyl.ps.gz Äàòà èçìåíåíèÿ: Sat May 21 01:43:06 1994 Äàòà èíäåêñèðîâàíèÿ: Sun Dec 23 02:50:04 2007 Êîäèðîâêà: IBM-866 Ïîèñêîâûå ñëîâà: ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï

The Restoration of HST Images and Spectra II
Space Telescope Science Institute, 1994
R. J. Hanisch and R. L. White, eds.
Image Restorations of High Photometric Quality
L. B. Lucy 1
Space Telescope ‘ European Coordinating Facility, European Southern Observatory,
KarlíSchwarzschildíStrasse 2, Dí85748 Garching bei M# unchen, Germany
Abstract. An image restoration technique is described that achieves high photometric
accuracy for both stars and distributed emission. The technique makes use of observerí
supplied information that some objects in the îeld are point sources and thereby eliminates
the ringing and photometric bias that arise with conventional restoration procedures. A
hierarchy of codes is described based on this twoíchannel decomposition of astronomical
images. The more sophisticated of these codes incorporate simultaneous estimation of the
PSF and are thus especially relevant for groundíbased imaging. Observing programs are
described in which these codes effectively transfer the resolution of HST to groundíbased
images.
1. Introduction
Astronomical images are commonly taken with the intent of making photometric measurements
of objects in the îeld. For preíCOSTAR HST images, the question then arises: should these
measurements be made on the original or on the restored image? When the objects of interest are
stellar, we might intuitively expect that higher accuracy will be achieved using the restored image.
The reason for this is that the sharpening of point sources should allow the photometry to be carried
out with a smaller aperture, thus reducing contamination by light from nearby sources and from the
background.
Unfortunately, two effects have hindered the realization of this expected gain in accuracy.
The îrst arises when the stars to be measured are superposed on an extended background source
# nebular emission or numerous unresolved stars. In this circumstance, nonílinear restoration
techniques # e.g., RichardsoníLucy (RíL) or Maximum Entropy (ME) # give rise to Gibbs
oscillations or #ringing# as they attempt the impossible task of recovering delta functions. Because
of this effect, image restoration does not necessarily increase the rate at which a photometric
measurement of a stellar object converges with increasing aperture.
The second effect is statistical bias and is of particular concern for ME restorations. When
entropy is deîned relative to a uniform prior (or default) image, the restored ME image departs
from uniformity because of real structure in the observed image, but only by the amount required
by the adopted goodnessíofíît criterion. Accordingly, the restored ME image is biased in favor of
uniformity (#atness), and this implies that the highest intensity peaks are underestimated.
For stellar photometry, these problems have led to the suggestion that the restored HST image
be used merely to identify and locate point sources, with the actual photometric measurements
being then made on the original image. Here, however, the possibility is explored of developing
more powerful restoration procedures that eliminate the above difîculties. Some years ago already,
Frieden and Wells (1978) implemented a scheme that effectively solved the îrst of the above
problems. Their scheme comprised the following discrete steps: 1) a model of the slowly varying
background is constructed; 2) this model is subtracted from the original image and the resulting
1 Afîliated to Astrophysics Division, Space Science Department, European Space Agency
79

80 Lucy
residual image is deconvolved; 3) the background model is added to the deconvolved residual image
to obtain the înal restored image. In this paper, a somewhat related procedure is described, leading
to a hierarchy of fully automatic codes of increasing sophistication which yield restorations of high
photometric precision for both stars and distributed emission. The performance of these codes is
tested and illustrated in Paper II (Hook & Lucy 1994).
2. Method
The method adopted might be described as image restoration incorporating designated sources
having speciîed shape and known positions. For astronomical images, the designated sources are
naturally taken to be point sources (stars, QSOs, AGNs), so the speciîed shape is the delta function.
Clearly, an observer can commonly, with high conîdence, identify many objects in an image as
point sources. Accordingly, it makes sense to develop codes that use this extra information and, at
least for these designated point sources, thereby eliminate the ringing associated with failed attempts
to recover delta functions.
In 1íD vector notation, the adopted model for \Psi(¦; j), the intensity distribution on the celestial
sphere, is
\Psi j = / j + / \Lambda
j : (1)
Here / j represents distributed emission and / \Lambda
j the designated point sources. This latter representaí
tion is achieved by deîning / \Lambda
j to be zero everywhere except at the locations (j = j 1 ; j 2 ; : : : ; j N )
of N point sources.
This model of the object can be mapped onto the image plane using the equation of image
formation, and the resulting prediction \Phi i for the intensity distribution in the image plane can
then be compared to the observed distribution ~
\Phi i . The posed image restoration problem is thus to
determine j and \Lambda
j in order to achieve an optimum ît of \Phi i to ~
\Phi i .
An obvious îrst thought is to adopt the RíL procedure and determine j and \Lambda
j with an iterative
scheme that asymptotically yields Maximum Likelihood estimates subject to the appropriate noní
negativity and normalization constraints. But this results in indeterminacy since an individual star
can be represented either in the \Lambda
j vector or by the corresponding singleípixel peak added to the j
vector, or indeed by any linear combination of these extremes.
To eliminate this indeterminacy, we must impose a constraint on j that excludes singleípixel
peaks. In other words, we must impose a resolution limit on j that is somewhat larger than the
limit implied by the numerical discretization. This is achieved by introducing regularization into
the optimization problem for determining j and \Lambda
j . Of the many possible regularization terms
(see, e.g., Titterington 1985), we choose the entropic form
S = \Gamma\Sigma
/ j
oe
ln / j
¼ j
; (2)
where oe = \Sigma/ j and where the #oating default (Horne 1985, Lucy 1994)
¼ j = \SigmaR jm m : (3)
Here the resolution kernel R jm is a bellíshaped function whose width is the required resolution
limit. Structure in j on a smaller scale results in a decrease in the entropy S and so is disfavored.
3. Code I: Known PSF
In image restoration, the PSF P (xj¦) is usually regarded as known. For this standard case, the
ideas of x2 lead to the following optimization problem: the quantities j and \Lambda
j are determined by
maximizing the objective function
Q( j ; \Lambda
j ; ï) = H + ffS + ï(\Sigma\Psi j \Gamma 1) : (4)

Image Restorations of High Photometric Quality 81
Here H denotes the reduced logílikelihood of the observed image and is given by
H = \Sigma ~
\Phi i ln \Phi i ; (5)
with
\Phi i = \SigmaP ij \Psi j : (6)
The regularization parameter ff controls the weight given to the regularization of / j , and the
Lagrange multiplier ï allows overall #ux conservation to be imposed.
4. Construction of Algorithms
The RíL algorithm (Richardson 1972; Lucy 1974) was initially derived from Bayes' theorem, but
in fact its use is really justiîed by the subsequent discovery (Lucy 1974; Shepp and Vardi 1982)
that every iteration increases the likelihood assigned to the observed image. Thus, in the usual case
where no stars are designated point sources, the objective function that is asymptotically maximized
by the RíL algorithm is
Q( j ; ï) = H + ï(\Sigma j \Gamma 1) : (7)
Knowing this, we can rewrite the RíL algorithm as an operation on Q. The result for the
increment / r+1 \Gamma r is
\Delta j = j
/
@Q
@/ j
! \Lambda
; (8)
where the asterisk indicates that ï has been evaluated by requiring that
\Sigma\Delta/ j = 0 : (9)
Eq. (8) suggests an operational procedure for deriving potentially effective algorithms: write down
the function Qwhose optimization deînes the restoration problem, then derive an iterative correction
scheme by applying Eq. (8).
With Q given by Eq. (4) and S by Eq. (2), the resulting algorithm simpliîes to
\Delta/ j = / j
/
C j \Gamma ff
@S
@/ j
!
(10)
and
\Delta/ \Lambda
j = / \Lambda
j C j ; (11)
where C j is the RíL correction factor.
5. Comments
The use and performance of Code I will perhaps be clariîed by the following remarks:
Initialization The positions of the designated point sources (j = j 1 ; j 2 ; : : : ; j N ) must be given
as well as initial amplitudes for j and \Lambda
j `
. These are taken to be constant and such that \Sigma/ j =
\Sigma/ \Lambda
j `
= 1
2 .
Normalization The integrated #ux of the point sources (/ \Lambda
j ) plus that of the distributed emission
(/ j ) is constrained to equal that of the observed image ( ~
\Phi i ). The allocation of #ux between these
two channels (initially 1:1, see above) is adjusted iteratively as Q is maximized.

82 Lucy
Regularization Applies only to the distributed emission. This implies that the converged pointí
source amplitudes / \Lambda
j `
are Maximum Likelihood estimates relative to the highly #exible model of
background emission provided by / j . This latter aspect is to be compared with the mathematically
simple background models used in standard photometric software packages.
SpatiallyíVariant PSFs As implemented, Code I assumes a spatiallyíinvariant PSF but only in
order to use FFTs. As with the RíL algorithm, this assumption is not otherwise necessary.
Point Sources These are modelled as single active pixels in a sea of zeroed pixels. If more accurate
centering is necessary, subípixelation is readily introduced (Lucy and Baade 1989; White 1990).
Such representations of delta functions facilitate the use of discrete FFTs and lead to simple code.
Nevertheless, for highest accuracy, continuous positions should be used.
Floating Default A ME restoration with uniform prior gives the most probable distribution of
photons in the restored image for a speciîed goodnessíofíît to the observed image. Although intuí
itively appealing, this solution, as noted earlier, is photometrically biased as a result of redistributing
photons from high to low intensity regions. Moreover, when the PSF is compact, this redistribution
is physically inconsistent with the independence of different #picture elements.# However, for a
#oating default computed with a compact kernel, this bias is greatly reduced (e.g., Lucy 1994) as is
this inconsistency.
Because the entropic form with #oating default is probably not derivable by combinatorial arí
guments, it should be regarded as a useful mathematical device and not ascribed a more fundamental
status than other regularization functions (Titterington 1985).
NoníNegativity As implemented, Code I maximizes Q subject to the constraints / j Ö 0 for all j
and / \Lambda
j ` Ö 0 for all `. There would be some merit in changing the second constraint to / j +/ \Lambda
j ` Ö 0
for all `. This would then allow / \Lambda
j `
! 0, which, though unphysical, avoids bias when estimating
magnitudes near the noise limit, as happens, for example, in following the fading of a SN. An
algorithm that achieves this has been constructed.
NoníDesignated Sources Given the resolution limit imposed on distributed emission, the question
arises as to how well point sources are represented that are not designated as such. This is investigated
in Fig. 1 where all emission derives from point sources but only one # the brightest # is designated
as such. Despite the resolution limit, the restored background (/) shows sharp peaks at the positions
of several stars. Accordingly, this model of the background could be used to identify and locate
additional point sources and thus form a step in an iterative procedure that maximally decomposes
the crowded îeld into individual stars, whose magnitudes are then derived from the înal solution
for / \Lambda
j `
.
6. Code II: Unknown PSF
In the standard case of known PSF (Code I), knowledge of the PSF is often derived from the
very image to be restored. This is of necessity true for groundíbased images where the seeingí
dominated and thus timeívariable PSF must be determined from suitably isolated stars in the îeld.
For such cases, a natural generalization of Code I is a code that incorporates these #PSF stars# into
the list of designated point sources and treats the PSF as an unknown function to be determined
simultaneously with j and \Lambda
j `
. Clearly, this generalization effectively mandates the assumption
of a spatiallyíinvariant PSF.
The objective function whose maximization deînes the above restoration problem is
Q( j ; \Lambda
j ; ae k ; ï; Ö) = H + ffS + ï(\Sigma\Psi j \Gamma 1) + Ö(\Sigmaae k \Gamma 1) : (12)

Image Restorations of High Photometric Quality 83
Figure 1. Oneídimensional illustration of application of Code I to crowded star îeld.
The brightest star at x = 0 is a designated point source, the others are not. The observed
image ( ~
\Phi) is plotted as a histogram and the predicted intensity distribution in the image
plane (\Phi) as a dashed curve. The continuous curve is the restored and regularized model
(/) of the distributed emission # background stars, in this case. Further individual stars
can be recognized in /.

84 Lucy
The new symbols introduced here are ae k the unknown PSF and a second Lagrangian multiplier Ö to
impose normalization of the PSF. Because the PSF is now assumed to be translationally invariant,
the predicted intensity distribution in the image plane is
\Phi i = \Sigmaae i\Gammaj \Psi j : (13)
As with Code I, an algorithm for maximizing Q is derived following the operational procedure of
x4. Details are omitted.
6.1. Remarks
Code II is not without pitfalls. The îrst is that Q has multiple maxima of approximately equal
height. In fact, if N point sources are designated, there are N + 1 such maxima, with each of the
N spurious maxima corresponding to the entire image ( ~
\Phi) being attributed to the PSFíbroadening
of just one of the point sources. These spurious solutions (m = 1; 2; : : : ; N ) thus have the form:
j = 0 all j, \Lambda
j `
= 0 for ` 6= m, / \Lambda
jm = 1, and ae k = ~
\Phi jm+k . Fortunately, these unwelcome
solutions are readily avoided with sensible initialization (see x5) and would in any case be readily
recognized as spurious.
A more serious problem arises when all the designated point sources are superposed on
distributed emission. The resulting PSF and the allocation of emission between point sources
and background are then determined by and sensitive to the regularization procedure. Limited
experiments suggest that sensible results in this circumstance require strong regularization with a
rather broad resolution kernel in Eq. (3).
Code II might seem to be an example of blind deconvolution (Nisenson et al. 1990) since both
the restored image and the PSF are derived simultaneously from a single observed image. However,
in contrast to blind deconvolution, the designated point sources here play an essential role in making
the PSF determinate. Accordingly, as noted above, this code should be regarded as consolidating
conventional imageíprocessing procedures into a single, automatic code.
7. Code III: Regularized PSF
Even when there are numerous designated point sources, the empirical PSF derived with Code II
will be noisy at low intensity levels. For the seeingídominated PSFs of groundíbased images, this
must be expected in the outer halo that surrounds the approximately Gaussian core (King 1971). An
obvious next step, therefore, is to regularize the PSF in order to damp out such noise #uctuations,
and this can be achieved simultaneously with the estimation of / j , / \Lambda
j `
and ae k if an extra entropic
term is added to Eq. (12). Moreover, if this new term also has a #oating default, its kernel can be
chosen to optimize the regularization by taking advantage of prior information about the PSF. For
example, apart from diffraction spikes, a groundíbased PSF should be circularlyísymmetric. It is
useful then to take the #oating default to be the azimuthallyíaveraged ae k . The înal ae k then departs
from circular symmetry only if the stellar images demand it # perhaps due to tracking errors #
and such distortions are corrected in the restored image.
Thus far only a 1íD experimental version of Code III has been tested.
8. Scientiîc Opportunities
The successful tests of Codes I and II reported in Paper II prompt thoughts as to the scientiîc
programs that could beneît from their application. Some of these might be described as softwareí
driven observing programs.
Host galaxies Conventional deconvolution of images of QSOs does not usefully contribute to
morphological studies of host galaxies because of the ringing centered on the deconvolved quasars.

Image Restorations of High Photometric Quality 85
But with Codes I or II, the quasar can be treated as a designated point source. This not only
eliminates ringing but also provides the restored and regularized image (/ j ) of the host galaxy with
quasar removed. With the alternative technique of PSF îtting and subtraction, the residual image
is not restored and also suffers from the Poisson noise associated with the bright quasar and left
behind on subtraction.
Measurement of H 0 Major programs underway with HST aim to discover Cepheids in relatively
distant galaxies and thus to determine H 0 using the locally calibrated period‘luminosity relation. In
these programs, the HST resolves crowded îelds into stars thus allowing Cepheids to be discovered
and photometered.
An alternative approach is the following: First, a single HST image is used to resolve the îeld
into stars and to determine accurate star positions. Multiple groundíbased images of the îeld are
then decomposed into stars using Codes I or II with the HST star positions as input. Finally, from
the values \Lambda
j `
on the different frames, Cepheids are discovered numerically and their periods found.
When so used, Codes I and II effectively transfer HST's resolving power to groundíbased
telescopes. Besides thus being highly economical of HST time, this technique has the further
advantage that the photometry is in wellíunderstood systems established with and for groundíbased
telescopes.
Imaging of Nebulae Many planetary nebulae and supernovae remnants occur in extremely crowded
star îelds, to the point that scientiîcally useful measurements are hard to make. Again Codes I or
II could be used with star positions from one HST image to obtain starífree groundíbased images
of such nebulosities in diagnosticallyíimportant emission lines.
Supernovae Light Curves Because of variable seeing, conventional îxedíaperture photometric
measurements of SNe are contaminated with variable amounts of light from the host galaxy. But
with Codes I or II the SN can be designated as a point source; its intensity is then determined relative
to a detailed model of the host galaxy with seeing taken into account.
References
Frieden, B. R., & Wells, D. C. 1978, J. Opt. Soc. Am., 68, 93
Hook, R. N., & Lucy, L. B., this volume (Paper II)
Horne, K. 1985, MNRAS, 213, 129
King, I. R. 1971, PASP, 83, 199
Lucy, L. B. 1974, AJ, 79, 745
Lucy, L. B. 1994, A&A, submitted
Lucy, L. B., & Baade, D. 1989, in 1st ESO/STíECF Data Analysis Workshop, P. Grosb#l, R. H.
Warmels, & F. Murtagh, eds., European Southern Observatory, Garching, 219
Nisenson, P., Standley, C., & Gay, D. 1990, in The Restoration of HST Images and Spectra, R. L.
White & R. J. Allen, eds., Space Telescope Science Institute, Baltimore, 103
Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55
Shepp, L. A., & Vardi, Y. 1982, IEEE Trans. Med. Imag., MIí1, 113
Titterington, D. M. 1985, A&A, 144, 381
White, R. L. 1990, in The Restoration of HST Images and Spectra, R. L. White & R. J. Allen, eds.,
Space Telescope Science Institute, Baltimore, 139