Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.stsci.edu/stsci/meetings/irw/proceedings/adorfhm.ps.gz Äàòà èçìåíåíèÿ: Sat May 21 00:57:34 1994 Äàòà èíäåêñèðîâàíèÿ: Sun Dec 23 05:18:54 2007 Êîäèðîâêà: IBM-866 Ïîèñêîâûå ñëîâà: ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï ï

The Restoration of HST Images and Spectra II
Space Telescope Science Institute, 1994
R. J. Hanisch and R. L. White, eds.
Towards HST Restoration with a SpaceíVariant PSF, Cosmic Rays and Other
Missing Data
HansíMartin Adorf
Space Telescope ‘ European Coordinating Facility, European Southern Observatory,
KarlíSchwarzschildíStrasse 2, Dí85748 Garching bei M# unchen, Germany
Abstract. The open problem of restoring full Wide Field and Planetary Camera (WF/PC)
image frames is considered. Several novel algorithms or modiîcations to existing algorithms
are described useful for restoring undersampled (multií) frames degraded by a spaceívariant
point spread function (SVíPSF) and by irregular faults. The algorithms comprise a faithful
rotation operator for sufîciently sampled images, an alternative method for treating reguí
lar and/or irregular missing data, and an algorithm for SVíPSF restoration combining the
classical sectioned restoration method of Trussel & Hunt with the generalized #coíaddition#
restoration method of Lucy & Hook. With these algorithms a comprehensive fullíframe
WF/PC restoration now appears to be feasible.
1. Introduction
The restoration of images obtained by the Hubble Space Telescope's îrst Wide Field and Planetary
Camera (WF/PCí1), is still an unsolved problem. This is mainly due to the spatial variability of
its point spread function (PSF), in conjunction with a variety of other factors. The purpose of this
contribution is to investigate the algorithmic requirements for a comprehensive spaceívariant PSF
(SVíPSF) restoration procedure, and to show that, based on the developments of the past three years,
a fullíframe WF/PC restoration has become feasible.
A comprehensive WF/PC restoration process (Fig. 1) has to address the following problems:
1. The pointíspread function is spatially variable and relatively large. Restoration with a SVí
PSF is known to be notoriously difîcult, and generally demanding in terms of CPUípower.
The PSF has to be computed from a parametric model, or somehow estimated from the data
frame(s). An efîcient algorithm must be used in view of the large data and PSF frame formats.
2. The data frames are undersampled. Observations must therefore be carefully designed in
terms of splits, subípixel shifts and rotations. Often the registration parameters must be
estimated postífacto, and different overlapping frames have to be combined.
3. Long exposures suffer from numerous cosmic ray (CR) hits. These have to be detected #
a daunting task in the presence of undersampling # and the missing data have to be treated
somehow.
4. There are other problems, such as detector contaminations dubbed #measles# and #atíîeld
uncertainties.
In the following several algorithms addressing problems (1) to (3) above are presented. Firstly,
the problem of restoring undersampled (multií) frames is discussed. It is argued that the cosmic ray
problem and the undersamplingproblem can be treated on the same footing. The proposed algorithm
for irregular missing data with masks largely follows treatments by Adorf et al. (1993), White
(1993), and recently by Freudling (1993), except for the upí and downísampling aspects. Secondly,
an algorithm for SVíPSF restoration is presented which combines the generalized #coíaddition#
72

SpaceíVariant PSF, Cosmic Rays, & Missing Data 73
Data Predictions
Masks
Ratios
PSFs
Correction
factor
Restoration
Figure 1. Flowíchart of a fullí#edged WF/PC restoration process. The data (upper
left) consists of several small, coarsely sampled patches (#sections#) taken from one or
more WF/PC exposures. Associated with each patch is a separate approximately spaceí
invariant (SI) pointíspread function (center). There is only one large, înely sampled joint
restoration (right) from which several coarsely sampled predictions are computed with
the help of the SIíPSFs. A prediction based algorithm for #onítheí#y# fault detection of
cosmic ray hits and other defects is integrated into the restoration scheme.
algorithm of Lucy & Hook (Lucy 1991a; Lucy & Hook 1992, 1993) with the sectioned restoration
method originally devised by Trussel & Hunt (1978a, 1978b). When designing these algorithms the
following principles have been adhered to (in order of decreasing priority): preserve the integrity of
the observational data; leave the RichardsoníLucy (RL) restoration method (Richardson 1972; Lucy
1974) intact as much as possible; strive for theoretically sound, yet efîcient, algorithms rendering
fullíframe WF/PC restoration practical.
2. Treatment of Undersampled Image Frames
In principle, both WF/PC cameras undersample their PSF. In the case of WF/PCí1 the undersampling
problem is usually obscured by the combined effects of spherical aberration and noise, rendering
many WF/PC frames sufîciently sampled in effect. For WF/PCí2 it is expected that, in view of
its sharp PSF, noníideal pixel response function, and reduced detector noise (Burrows 1994), the
undersampling problem will be much more prominent.
The anticipated problems associated with undersampled image frames have been addressed in
a number of articles (Adorf 1989a, 1989b, 1990). Fosbury (1993) discusses several image sampling
strategies for the proposed future Advanced Camera for HST, and lists recent references concerning
undersampled multiíframes. The major difîculty with any undersampled data frame is that it does
not represent a continuous image. Thus generally, it cannot be shifted, stretched or rotated. 1
In order to overcome the undersampling problem, in the context of HST, it is generally
recommended to obtain two or more data frames of the same îeld (i.e., multiíframes) with different
1 Even when the data are sufîciently sampled,the conventional rotation algorithms implemented in the standard astronomical
image processing packages are generally inadequate,since they are not faithful to the data. An algorithm is called #faithful#
here if it preserves the information content (in particular, the high spatial frequency spectrum) of the data. Every faithful
operator has an inverse operator.

74 Adorf
(e.g., shifted or rotated) sampling patterns, and to subsequently combine them. The generic
procedure for reconstructing an image from a set of undersampled multiíframes consists of carrying
out a joint restoration following the standard LucyíHook recipe, which internally requires algorithms
for shifting, stretching, or rotating the model (= best restoration) in order to map it onto the sampling
patterns of the different observations. It is important that these mappings are carried out faithfully.
Of these operators the rotation operator is the most difîcult one to realize in a faithful manner.
2.1. A Novel Faithful Rotation Algorithm
The rotation algorithm described here strictly preserves the spatial Fourieríspectrum and is based on
the fact that a 2 \Theta 2 rotation matrix R(`), which rotates an image by an angle `, can be decomposed
in several ways, including one consisting of three consecutive 2 \Theta 2 shear matrices
R(`) = S x (\Gamma`=2) S y (arctan(sin(`)) S x (\Gamma`=2) (1)
Here S x (`) and S y (`) designate a shear by an angle ` in the xí and yídirection, respectively.
Since a sufîciently sampled image can be faithfully sheared using efîcient 1D sinc interpolaí
tions (based on fast Fourier tranforms, FFTs), the whole rotation operation can also be carried out
faithfully. The only problem to watch out for are potential Nyquist frequency components, which
can only appear in evenísized frames. Since the three shears, out of which the rotation is composed,
are carried out independently on individual rows (or columns), the faithful rotation algorithm is
particularly suited to a vectoríparallel computer architecture.
The rotation algorithm can be incorporated into the generalized LucyíHook #coíaddition#
scheme, in order to faithfully rotate the current best restoration from model to observational space,
and then to rotate the correction factor back
OE k
(x) = p(x) ffi [R(`)y k
(x)] (2)
c k (x) = R(\Gamma`)[ae k (x) \Lambda p(x)] (3)
Note that neither the WF/PC PSF nor the pixel response function are isotropic.
2.2. Restoration to a Finer #Subsampled# Grid
The standard implementations of the RLíalgorithm include a subsampling option (Lucy 1991b;
White 1993), so the problem of restoring to a îner grid is widely considered to be #solved.#
However, it has not been proven so far that the RLíalgorithm with subsampling is unique or has
desirable properties when unaliasing or reconstructing missing high spatial frequency components.
When the RLímethod with subsampling option is used to restore an undersampled #image#
frame, it will create a restoration that, while consistent with the data, is not smooth, but rather
maximally pixelated (see the Appendix). So the question emerges which method to use in the
case of an undersampled optical PSF, and a detector with noníideal pixel response, as expected for
WF/PCí2.
The following observations show that there is indeed some opportunity for algorithmic develí
opment:
1. The sampling theorem states: It is possible to exactly reconstruct a continuous bandlimited
signal from discrete samples, if the signal is sufîciently densely sampled (and noiseífree).
The sampling theorem, however, does not state: an undersampled continuous signal cannot
be reconstructed from discrete samples.
2. An undersampled signal is data with missing values (e.g. every other data row and column is
missing).
3. One possesses some prior knowledge, to help in îlling in the missing values:
‘ the continuous signal is nonínegative (apart from noise);

SpaceíVariant PSF, Cosmic Rays, & Missing Data 75
‘ the observed signal is strictly bandlimited (apart from noise);
‘ the PSF is known, and since its Fourierítransform (the optical transfer function, OTF),
has a înite support, it always stretches beyond more than one pixel (it is in fact inînitely
large).
4. The interplay between the (linear) bandílimitation and the (nonílinear) nonínegativity coní
straint is not yet completely understood. So there is latitude for theoretical developments and
experimentation.
In order to overcome the pixelation effect of the classical RLíalgorithm, the following modiîed
RLíalgorithm is suggested which operates solely on the îne grid:
1. Create a static mask m(x) which is 1 where data values have been observed and 0 elsewhere.
2. Replace the predicted values OE k (x) by the observed ones OE \Lambda (x) wherever the latter are known,
i.e., where the mask m(x) = 1 and apply other constraints such as nonínegativity (projection
onto convex set, POCS); a dynamically relaxing bandílimit constraint should presumably be
included here.
3. Proceed with a standard RLírestoration.
This algorithm allows some direct control over the extent of high frequency components introduced
by the restoration algorithm in its unaliasing attempt.
2.3. Detection and Treatment of Cosmic Ray Hits
The cosmic ray (CR) problem involves two separate parts: detecting CRíhits and somehow #îlling
in# the missing data.
Several standíalone algorithms have previously been published to detect CRíhits (Murtagh
1992; Murtagh & Adorf 1992; Adorf et al. 1993). However, it is clear that if a restoration is carried
out anyway, the CRícontamination problem is best treated within that restoration process (Freudling
1993; White 1993). For sufîciently sampled images these algorithms work well. However, for
undersampled images they are affected by the pixelation problem described above.
Let us, for a moment, consider how the human visual system would înd cosmic ray hits:
suspicious pixels are compared with their respective neighbourhoods. More precisely, the brain,
using PSFíknowledge, predicts how the observation ought to look, and the eye compares this
prediction with the actual data.
The essence of this detection method can be captured in mathematical terms, using Lucy's
(1974) original notation, as: The prediction OE k (x), computed from the current best restoration
k (x) via convolution with the PSF, represents our current understanding of the observation. Thus
any faults in the observed data OE \Lambda (x) are best detected by comparing the data OE \Lambda (x) with the
prediction OE k (x). We are led to the following spatially adaptive, dynamic ßíoe clipping algorithm
for detecting CRíhits (and potentially other faults):
1. Estimate a map of spatially variable standard deviations oe(x). If the Poissonínoise assumption
holds, one can simply take the square root of the signal counts. Otherwise one may estimate
oe(x) from the spatial variance within a small neighborhood around location x.
2. Flag all those pixels whose data values OE \Lambda (x) are outside OE k (x) \Sigma ßoe(x) thereby producing
a dynamic mask m k (x), which is 1 where the data OE \Lambda (x) is #good# judged via OE k (x), and 0
elsewhere. Here ß denotes a threshold value.
Since the map of standard deviations oe(x) may be corrupted by initially missed CRs, it is best made
dynamic too by updating it regularly after every m th RLíiteration. This double dynamization of the

76 Adorf
variance estimation and CRídetection steps should allow to optimally detect CRíhits even on single
exposures down into the noise sea. The problem to watch out for is dynamical instability.
Several methods are available to treat CRíhits once detected. One might replace the predicted
values by the observed ones wherever the latter are known (projection onto convex set) and proceed
with a quasi standard RLírestoration. In mathematical notation this algorithm reads:
OE k
(x) = p(x) ffi / k
(x) (4)
OE \Lambdak
(x) = m k (x)OE \Lambda
(x) +
i
1 \Gamma m k (x)
j
OE k (x) (5)
r k (x) = OE \Lambdak (x)=OE k (x) (6)
c k (x) = r k (x) \Lambda p(x) (7)
/ k+1
(x) = c k
(x) k
(x) (8)
This method of îlling in missing data values 2 might also be beneîcial in situations when the
information is not completely destroyed, e.g., in the case of saturated pixels, where the recorded
value can kept as a lower bound for a number of iterations, but later on would be left free to be
inferred from the data (White, pers. comm.).
3. Restoring a WF/PC Frame with a SpaceíVariant PSF
Let us now turn to the other main topic of this contribution, namely an algorithm for restoring a
full WF/PC frame (or undersampled multiíframes) with its SVíPSF. The wellíknown #bruteíforce#
SVíPSF restoration method is the #sectioned# algorithm (Trussel & Hunt 1978a, 1978b) which
works as follows:
1. Subdivide the observation into a regular set of M patches (the #sections#) with data OE \Lambda
ï (x).
2. Restore each section under the hypothesis that there is a PSF p ï (x) which is spaceíinvariant
(SI) within patch ï.
3. Sew the M restorations k
ï (x) together.
The spatial density of patches (more precisely patch centre locations) necessary to carry out a
restoration depends on the degree of spatial variability of the PSF. It is clear that the patches must
overlap each other by at least the effective PSF diameter in order to exclude border effects, thereby
limiting the minimum size of the patches.
The main problem with the simple sectioned approach is the danger of discontinuities occurring
in the restoration at the patch boundaries, particularly when the PSF is asymmetric and changes its
morphology, as in the case at least for WF/PCí1.
Here is an elegant way out: if one considers the patches as #independent# overlapping obserí
vations of the same îeld, each with its own spaceíinvariant (SI) PSF, then one can simply use the
LucyíHook #coíaddition# algorithm to construct a single restored image k (x) from the multiple
data patches OE \Lambda
ï (x). A spatially variable, #fuzzy# weighting function w ï (x) can be used in order to
combine the correction factors c k
ï (x) belonging to the different patches ï
c k (x) =
X
ï
w ï (x)c k
ï (x) (9)
2 This simple and efîcient heuristic algorithm can be (and hasbeen) criticized for violating the îrst principle of data analysis
#Don't modify the data.# However, this criticism can be rebutted, because bad pixels never attain the same status as good
ones.

SpaceíVariant PSF, Cosmic Rays, & Missing Data 77
It is obvious that the formula above may be generalized to multiple input data frames.
Presently, it is an open question which weighting function w ï (x) to best use. A #hard# hatí
shaped function may again lead to the problem of potential discontinuities at the section boundaries.
Triangular or cosine bellíshaped weighting functions are better behaved in that respect. One could
possibly #derive# a weighting function by crossícorrelating a series of interpolating PSFs and use
the correlation coefîcient as a guide to the spatial variability of w ï (x). Or one might apply a
minimum meanísquareíerror criterion to optimize the choice of the functional form of the w ï (x).
The proposed SVíPSF algorithm can be generalized in various ways. One could adapt the
spatial density (and corresponding size) of the sections to the local signalítoínoise (S=N ) level,
since it does not make sense to restore low S=N regions very accurately. Another way of saving
CPUítime might be to use a denser patch grid where the PSF changes more rapidly, e.g., at the
frame edges, or in areas of increased astronomical interest (White, pers. comm.).
Of course, the success of the proposed SVíPSF algorithm hinges upon the ability to produce
suitable PSFs. Whether one has to rely on theoretical models or use one of the #blind deconvolution#
PSFíestimation techniques (Lucy 1994) is currently an open question.
A înal remark: The proposed SVíPSF algorithm is ideally suited to a coarseígrained vectorí
parallel computer architecture such as TMC's Connection Machine 5.
4. Summary
Two problems related to WF/PC imagery have been addressed: restoration in the presence of
undersampling and restoration in the presence of a spaceívariant pointíspread function.
It has been argued that the process of RichardsoníLucy (RL) restoration with subsampling is not
fully understood in the presence of undersampling and that further work is required. A POCSíbased
variant of the RLíalgorithm with subsampling has been suggested for further examination. Also an
alternative modelíbased algorithm for detecting irregular random missing data such as cosmic ray
(CR) hits has been proposed. The algorithm is based on a spatially adaptive, dynamic ßíoe clipping,
which should be optimal, at least for singleípixel CRíhits.
A novel, easily parallelizable algorithm for faithful image rotation has been stated. It is based
on a sequence of three successive shears which in turn can efîciently be implemented using fast
Fourier transforms (FFTs). The rotation algorithm is required within a joint restoration of multiple
rotated WF/PC frames.
Finally a novel algorithm has been proposed for restoring fullíframe WF/PC images distorted
by a spaceívariant point spread function (SVíPSF). The algorithm is efîcient and parallelizable. The
algorithm allows several interesting generalizations which promise to either increase performance
or restoration quality or both.
With these algorithms a comprehensive fullíframe SVíPSF restoration of single or multiple
WF/PC frames now appears to be feasible.
Appendix
Some simulations have been carried out to investigate how much pixelation the conventional RLí
method with subísampling option (where the dataítoíprediction ratio is upsampled by pixel replií
cation) imprints onto the restoration of an undersampled distorted image frame.
Two grids were used, a "îne" grid with 64 \Theta 64 pixels and a "coarse" grid with 32 \Theta 32 pixels,
i.e., half the sampling rate of the îne grid. As an image model a single sine wave was used with a
wavelength of 6 pixels on the îne grid, which is still sufîciently sampled on the coarse grid (2/3 of
the coarse grid's Nyquist frequency). A value of 1 was added to create a nonínegative image.
Four different PSFs were used: a delta function (PSF0) and 3 diffraction limited PSFs with
"aperture" radii of 31, 15, and 7 pixels (PSF1 to 3), respectively. PSF1 is approximately twoífold

78 Adorf
(fourífold) undersampled on the îne (coarse) grid. PSF2 is sufîciently sampled on the îne grid,
but about twoífold undersampled on the coarse grid. PSF3 is sufîciently sampled on both grids.
The model was convolved with each of the PSFs, then 2 \Theta 2 blockíaveraged and downísampled
(decimated) to 32 \Theta 32 pixels, generating four simulated data frames (REST0 to 3). No noise was
added. The STíECF's IRAF code for the RLímethod with subísampling option was used to restore
the four simulated "data" frames.
As expected, the restoration REST0 (after some 20 accelerated iterations) clearly showed
the coarseígrid pixel pattern: pairs of two pixels on the îne grid had exactly the same value.
The restoration REST3 on the other hand showed very little pixelation, if any. The intermediate
restorations REST2 and REST1 showed a degree of remnant pixelation on the îne grid which
increased with the degree of undersampling.
Acknowledgments. Thanks are due to Bob Hanisch, STScI, for suggesting the SVíPSF restoraí
tion problem and to Richard Hook, STíECF, who helped to improve the written version of this
contribution. I am particularly indepted to Rick White, STScI, for enlightening discussions and
several valuable suggestions.
References
Adorf, H.íM. 1989a, STíECF Newsletter, 12, 9
Adorf, H.íM. 1989b, in 1st ESO/STíECF Data Analysis Workshop, P. J. Grosb#l, F. Murtagh, & R.
H. Warmels, eds., European Southern Observatory, Garching, 215
Adorf, H.íM. 1990, in Errors, Bias, and Uncertainties in Astronomy, C. Jaschek & F. Murtagh, eds.,
Cambridge University Press, Cambridge, 71
Adorf, H.íM., Hook, R. N., & Lucy, L. B. 1993, in Space Astronomical Telescopes and Instrumení
tation II, SPIE, in press
Burrows, C. 1994, this volume
Fosbury, R. 1993, in The Future of Space Imaging, R. Brown, ed., Space Telescope Science Institute,
Baltimore, 83
Freudling, W. 1993, STíECF Newsletter, 20, 8
Hook, R. N., & Lucy, L. B. 1992, STíECF Newsletter, 17, 10
Hook, R. N., & Lucy, L. B. 1993, STíECF Newsletter, 19, 6
Lucy, L. B. 1974, AJ, 79, 745
Lucy, L. B. 1991a, STíECF Newsletter, 16, 6
Lucy, L. B. 1991b, in The Restoration of HST Images and Spectra, R. L. White & R. J. Allen, eds.,
Space Telescope Science Institute, Baltimore, 80
Lucy, L. B. 1994, this volume
Lucy, L. B., & Hook, R. N. 1992, in Astronomical Data Analysis Software and Systems I, ASP
Conference Series 25, D. M. Worrall, C. Biemesderfer, & J. Barnes, eds., 277
Murtagh, F. D. 1992, in AstronomicalData Analysis Software and Systems I, ASP Conference Series
25, D. M. Worrall, C. Biemesderfer, & J. Barnes, eds., 265
Murtagh, F. D., & Adorf, H.íM. 1992, in Data Analysis in Astronomy IV, V. Di GesÑu, et al., eds.,
Plenum Press, New York, 103
Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55
Trussel, H. J., & Hunt, B. R. 1978a, IEEE Trans. Acoust. Speech Sig. Proc., 26, 157
Trussel, H. J., & Hunt, B. R. 1978b, IEEE Trans. Acoust. Speech Sig. Proc., 26, 608
White, R. L. 1993, in Restoration Newsletter, R. J. Hanisch, ed., Space Telescope Science Institute,
Baltimore, 1, 11