A. Llebaria, P. Lamy, P. Malburet
Laboratoire d'Astronomie Spatiale du CNRS, BP 8, F-13376 Marseille
Cedex 12
The accurate detection and subsequent removal of point-like objects in astronomical images is a frequent problem. We mention, as examples, the case of stars pervading an extended objects, or cosmic ray impacts on CCD images. The problem is often complicated by various distortions affecting the point-like objects, such as, saturation, variable Point Spread Function (PSF) in the field-of-view, and irregular images which preclude the use of fitting methods. We propose a new approach to the problem based on the properties of surface curvature which is related to peaks, pits and flats in images and therefore contain the key information to discriminate point-like objects.
The intensity of an image can be viewed as a two-dimensional surface in a three-dimensional space with the two spatial coordinates forming two of the three dimensions and the intensity axis forming the third dimension. Surface curvature will describe how much this surface is bending in a given direction at a particular point; it involves the first and second partial derivatives with respect to the spatial coordinates and thus involves the spatial neighborhood of the point. It is readily understandable that surface curvature will turn out be a powerful means of discriminating objects in images.
Practical experiments with the surface curvature of astronomical images
have quickly revealed the superiority of dealing with the log of the
intensity, , rather than with the intensity itself. The
justifications of this transformation are easily understood because, (1)
the range of curvatures will be narrower on the
images
than on the intensity images, insuring a better detection of faint objects, and
(2) point-like objects, such as stars may be viewed, at least in a first
approximation, as having a Gaussian shape. The logarithmic transformation
will generate similar paraboloids (i.e., having the same parameters)
for which the maximum of curvature takes place at their tops. Practically,
the logarithmic transformation of the star images will exhibit a coarse
paraboloidal shape which can be easily localized from the maximum of
curvature.
A detailed analysis of surface curvature can be found in Peet and Sahota
(1985). This analysis shows that at each regular point of a surface we can
easily determine and
, the local minimum and maximum curvatures,
respectively.
and
are invariant under rigid motions of the surface.
Two other curvatures are often defined; the Gaussian or total curvature
, and the mean curvature
; but it turns out that the
principal curvatures
and
are best suited to the detailed
characterization required by our problem.
As curvature is a geometric property of continuous surfaces, we must apply a
continuous model to the discrete array of 2-D data in order to deduce
local derivatives, and from this the local curvatures. Even if more complex models
can be used, a simplified facet model of second order (Halarick 1984), is good
enough when a small set of points is involved in the estimation.
For the 33 patch window the practical calculation of
and
is
extremely simple. The surface
is approximated by a quadric
and the coefficients are determined using the least squares fit to the real
surface on a
neighborhood of the pixel under consideration.
It is quite convenient to discuss the surface curvature in the ()
space as illustrated in Figure 1. Using the above convention,
,
and only the corresponding half plane is therefore allowed. Also, according
to the above definitions, if a curve on the surface is concave up, the
surface curvature in that direction is positive and vice-versa. The allowed
half-plane of Figure 1 (left side) is divided in three regions by the
and
axis. The two axis themselves correspond to parabolic points, with ridges on
the
axis (
), and valleys on the
axis (
). The second quadrant
, is the region of the saddle points.
The elliptic points are located in the remaining two half-quadrants,
distinguished by their concavity: pits for
, and peaks for
. Obviously, a plane is at the origin
.
Figure:
(Left side).Surface curvature visualized in the plane of the two principal
curvatures and
.
(Right side).A practical
,
diagram for an astronomical image. The
regions of interest are related to specific parts or curvatures of real
objects.
Original PostScript figure (53 kB)
Let us now consider how the above theoretical considerations apply to real
astronomical images. Figure 1 (right side) is a sketch of the ()
diagram based on a CCD image which exhibits an extended source; the tail of
comet P/Halley, a rich star field, and various defects and artifacts. The
features of this (
) diagram are sufficiently general to localize the
various domains of interest and to establish the appropriate criteria for the
discrimination.
The main feature of the diagram is a ``butterfly''
pattern composed of two wings A and B and a central core C. Both wings
pertain to stars, wing A to their top portions and wing B to their bottom
portions. These latter parts are fairly asymmetric saddle
points leaning to valleys, and explain the closeness of the wing B to
the axis. In fact, the orientation of the wing, as given for instance
by the slope
of their symmetry axis, is a very good estimate of
the anisotropy of the surface curvature. For instance, elongated
(trailed) star images will induce a rotation of wing A. The central core C
obviously correspond to the flat part of the image; essentially the
extended object and the background. Finally, two additional features are
conspicuous in the diagram, the blobs E and F respectively localized on
the
and
axis. They correspond to elongated structures in the
original image which involve both ridges (top parts) and valleys (bottom
parts). These structures are typically bad columns in the CCD and saturated
star images which spill off along the column direction.
It is now a very simple matter to establish criteria for discriminating objects
in the original image. These criteria are written in terms on conditions on
and
or on a combination of the two. From the initial image we
define a mask (a binary image) defining the subset of
pixels in the original image which belongs to the ``flat'' regions. The
complementary subset includes all invalid points from stars, saturations, etc.
The purpose of the next and final phase is to replace patches of invalid pixels in the original intensity image by a local interpolation, the so-called pyramidal interpolation, using a non-linear multiresolution method which resembles the Haar wavelet analysis and synthesis.
The procedure starts by shrinking the original pixels image to a
image. Each pixel in the latter corresponds to an area
of
pixels in the former. Their value is a mean of valid pixels
in the
related area. If 3 or 4 pixels in this area are invalid,
the mean is classed as invalid. After this first pass, we obtain an
image of means and a corresponding binary image or mask of valid or invalid
means. Repeating this procedure builds another series of images
(means and mask) of lower resolution. The image of means is called the analysis
image.
The procedure is repeated l times until either the extinction of patches of
invalid pixels is obtained, or or
are too small
(that is equal to or
less than 4 pixels on a side). If some final pixels are invalid, they are
replaced by values estimated from a global adjustment of a quadratic surface
over the remaining valid pixels. After this, all pixels are valid at
the lowest resolution level and the synthesis process can start.
The synthesis is done by (1) expanding the synthesis image built in the
previous l-1st step by , (2) substitute invalid pixels in the
analysis image at this step l with the corresponding pixels from the above
expanded image, and take this image (with substituted pixels) as the new
synthesis image at this step.
The procedure stops when the initial resolution is recovered. The method
reconstitutes, in each patch of invalid pixels, the low resolution estimate
achieved with valid pixel neighbors.
This procedure needs a minimal percentage of valid pixels ( to
obtain significant results. When the number of valid pixels is very low,
the threshold to obtain a significant mean from a
area in the
analysis step can be reduced to one valid pixel.
We have presented a new method for the discrimination of point-like objects using surface curvature. The characterization by the curvature offers a very general and flexible method which can be adapted to any specific problem without any necessary a priori knowledge of the shape of the point-like objects. A few trials are basically what is needed to determine the threshold for discrimination. Another potential application is the removal of cosmic rays events from single exposure CCD images (when the exposure has been duplicated, the preferred approach is obviously the comparison of the two images). These impact events appear as very sharp peaks which may be unambiguously discriminated by their curvature.
Peet, F. G., & Sahorta, T. S. 1985, IEEE Trans. Pattern Anal. Machine Intelligence, Vol. PAMI-7, 734
Halarick 1984, IEEE Trans. Pattern Anal. Machine Intelligence, Vol. PAMI-6, 58