J.-L. Starck
CEA, DSM/DAPNIA, CE-Saclay, F-91191 Gif-sur-Yvette Cedex, France
F. Murtagh
Space Telescope--European Coordinating Facility, European Southern
Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany
Affiliated to Astrophysics Division, Space Science Department,
European Space Agency
A. Bijaoui
Observatoire de la Côte d'Azur, B.P. 229, F-06304 Nice Cedex 4, France
Extensive literature exists on the wavelet transform and its
application (Chui 1992; Daubechies 1992; Meyer 1989). A discrete
wavelet transform approach can be obtained from multiresolution
analysis (Mallat 1989).
Multiresolution analysis results
from the embedded subsets generated by interpolations at
different scales.
A function is projected at each step j onto the subset
. This
projection is defined by the scalar product
of
with the scaling
function
which is dilated and translated:
is a scaling function which has the property
Equation permits
the set
to be computed directly from
.
If we start from the set
, we compute all the sets
, with j>0, without directly computing any other scalar
product:
At each step, the number of scalar products is divided by 2;
the signal is smoothed and information is lost. The remaining
information can be restored using the complementary subspace of
in
.
This subspace can be generated by a suitable wavelet function
with translation and dilation:
We compute the scalar products with
With this analysis, we have built the first part of a filter bank. In
order to restore the original data, Mallat uses the properties of
orthogonal wavelets, but the theory has been generalized to a
large class of filters (Cohen, Daubechies, & Feauveau 1992).
Filters and
, the conjugates of h and g, have been introduced
(Daubechies 1992), and the restoration is performed with
In order to get an exact restoration, two conditions are required for the conjugate filters: (1) the De-aliasing condition:
and (2) exact restoration:
Many sets of filters have been proposed, especially for coding. It has been shown (Daubechies 1992) that the choice of these filters must be guided by the regularity of the scaling and the wavelet functions. The complexity is proportional to N. The algorithm provides a pyramid of N elements. The 2D algorithm is based on separable variables leading to x and y directions being prioritized. This implies a non-isotropic analysis, which often runs counter to physical patterns.
Feauveau (1990) introduced the quiconx analysis according to
Adelson (Adelson, Simoncelli, & Hingorani 1987).
This analysis is not dyadic and
allows an image decomposition with a resolution factor equal to .
By this method, we have one wavelet image at each scale, and not three
as with the previous method.
If has a cut-off frequency, the discrete filter
used to compute the coefficients at a given resolution from the
previous one (Starck & Bijaoui 1994), is
.
Therefore, between two resolutions,
the cut-off frequency is divided by 2, allowing us to reduce the number
of samples.
The information lost in the filtering is given by the discrete filter
defined by
.
The signal is reconstructed with the two filters
,
:
The described algorithm is easily adapted to image processing. As opposed
to classical multiresolution analysis, it is a decomposition
with isotropic scaling and wavelet functions.
The frequency band is also reduced by a factor 2 at each step.
Applying the sampling theorem, we can build a pyramid of
elements.
For image analysis the number of elements is
. The
increase in storage is small.
The `` à trous'' algorithm is a powerful algorithm for the following reasons: (1) the computational requirement is reasonable, (2) the algorithm is easy to program, (3) in two dimensions, the transform is practically isotropic, (4) we can use compact scaling functions, (5) the reconstruction algorithm is trivial, (6) the transform is known at each pixel, allowing detection without any error, and without interpolation, (7) we can follow the evolution of the transform from one scale to the next, and (8) invariance under translation is completely verified.
Details of the algorithm are given in Bijaoui, Starck, & Murtagh (1994)
and Shensa (1992).
The wavelet transform of an image by this algorithm
produces, at each scale j, a set which
we will call a wavelet plane throughout the following discussion. This has
the same number of pixels as the image. The original image
can be expressed as the sum over wavelet planes and the
smoothed array
The pyramidal algorithm (Starck 1993) is derived from
the `` à trous'' algorithm.
At each scale, is obtained by computing the difference
between
and
(where F denotes a linear filtering).
But
is equal to
(U being the undersampling
or decimating operator),
and not equal to
(we first undersample
by keeping one
pixel
out of two). This transform produces data similar to those obtained
by the Burt & Adelson (1983) pyramidal Laplacian algorithm, but
it is done with one wavelet. Burt's algorithm needs three wavelets and is
therefore
not isotropic (Bijaoui 1991). The number of elements is
,
as in the case of the transform using the FFT, but this method
has the advantage to allow all computation to be carried out in direct space.
An iterative algorithm is necessary for an exact restoration.
The 2-dimensional extension of Mallat's algorithm leads to a wavelet transform with three wavelet functions (three wavelet coefficient sub-images at each scale). An isotropic wavelet seems more appropriate, specially in astronomical imaging where objects are often isotropic (e.g., stars).
Mallat's and Feauveau's methods provide a remarkable framework to code a signal, and especially an image, with a pyramidal set of values. But, contrary to the continuous wavelet transform, these analyses are not covariant under translation. At a given scale, we derive a decimated number of wavelet coefficients. We cannot restore the intermediate values without using the approximation at this scale and the wavelet coefficients at smaller scales. Since the multiresolution analysis is based on scaling functions without cut-off frequency, the application of the Shannon interpolation theorem is not possible. The interpolation of the wavelet coefficients can only be done after reconstruction and shift. This is unimportant for a signal coding that does not modify the data, but is not the same if we want to analyze or restore an image.
If the image I we want to analyze is the convolution product of
an object O by a point spread function (PSF) (), we
have
where are the wavelet coefficients of z, and a is the scale parameter.
We deduce that
We can directly analyze the object from the wavelet coefficients of the image. But due to decimation effects in Mallat's, Feauveau's, and pyramidal methods, this equation fails. The wavelet transform using the FFT (which decimates using Shannon's theorem) and the `` à trous'' algorithm (which does not decimate) are the only ones which respect the scale separation property.
By definition, the wavelet coefficient mean is null. Every time we have a positive structure at a scale, we have negative values surrounding it. These negative values often create artifacts during the restoration process, or complicate the analysis. For instance, if we threshold small values (noise, insignificant structures, etc.) in the wavelet transform, and if we reconstruct the image at the full resolution, the structure's flux will be modified. Furthermore, if an object is very high, the negative values will be important, too, and will lead to false structure detections.
We often have bright point objects in astronomical imaging (stars, cosmic ray hits, etc.) The convolution of a Dirac function by the wavelet function is equal to the wavelet function, so at each scale, and at each point source, we will have the wavelet. Cosmic rays can pollute all the scales of the wavelet transform.
There is no ideal wavelet transform algorithm---the selection will depend on the application. The two last problems cannot be solved by the wavelet transform, and they lead to other multiresolution tools which we now present.
The search for new multiresolution tools was motivated by problems related to the wavelet transform. We would prefer that a point structure (represented in one pixel in the image) is present only at the first scale, and that a positive structure in the image not create negative values in the multiresolution space. We will see how such an algorithm can be found, using morphological filters such as the median filter.
By modifying the ``à trous'' algorithm, we easily obtain the desired transform. The algorithm becomes:
The reconstruction is carried out by a simple addition of all the scales:
Such an algorithm has several advantages: (1) the transform can be carried out with integer values, (2) energy in structures at each scale is real, and is not modified by the treatment, (3) structure contours are better respected, (4) there are no negative values around positive structures, and (5) the algorithm can easily be modified to work at intermediate scales (in order to have a non-dyadic analysis). We just multiply the size s at the step 4 by a coefficient different from 2 (if we want half resolution, we multiply by 1.5).
However this algorithm has a big disadvantage: computation time. The mask increases by two at each scale, and we have to sort a number of values which increases considerably. To solve this problem, we introduce decimation.
In this transform, the mask size is constant, and there is no time computation problem anymore. However, the reconstruction is not exact, and an iterative procedure has to be introduced during the transformation or at the reconstruction. If we choose the iterative transform, we have:
Images generally contain noise (Gaussian or Poisson, usually)
and hence the wavelet coefficients
are noisy, too. In most applications, it is necessary to know whether a
coefficient is due to signal or to noise.
The wavelet transform yields a set of resolution-related views of the
input image. A wavelet image plane at level j has coefficients
given by .
If we obtain the distribution of the coefficient
for each plane,
based on the noise, we can introduce a statistical significance test for this
coefficient.
Given stationary Gaussian noise, it suffices to compare to
, where
is the standard deviation of wavelet plane j,
and k is often chosen as 3.
If
is small it is not significant, and could be due to noise. If
is large, it is significant.
If the noise in the data I is Poisson, the transform
acts as if the data arose from the
Gaussian white noise model with unit standard deviation (Anscombe 1948).
Taking the wavelet transform of
,
will
be significant if
is above a given threshold.
(The superscript on the wavelet coefficients indicates the image to which
the wavelet transform was applied.)
The appropriate value of
in the succession of wavelet planes is assessed
from the standard deviation of the noise
in the original image
and from study of the noise in the wavelet space. Details of this study
can be found in Starck & Murtagh (1994).
We will say that a multiresolution support (Starck, Bijaoui, & Murtagh 1994)
of an image describes in a logical or boolean way whether
an image I contains information at a
given scale j and at a given position .
If
(or true), then I contains information at
scale j and at the position
.
M depends on several parameters:
(1) the input image,
(2) the
algorithm used for the multiresolution decomposition,
(3) the noise, and
(4) all additional constraints we want the support to satisfy.
Such a support results from the data, the treatment (noise estimation, etc.), and from knowledge on our part of the objects contained in the data (size of objects, linearity, etc.). In the most general case, a priori information is not available to us.
The multiresolution support of an image is computed in several steps:
The multiresolution support indicates where the information is, and allows all information we have about the image (data, noise, catalog information, etc.) to be integrated. It can then be used to visualize or to further process the data. The fact that we can add constraints and laws to the support makes this a very convenient data structure for introducing information into the treatment.
The visualization of the wavelet transform allows us to have better knowledge of the structure of the image. The `` à trous'' algorithm is particularly efficient, because each scale has the same size as the original image, and can be visualized, processed, and compared very easily. See Starck (1993) and Bijaoui, Starck, & Murtagh (1994) which use various ways to visualize a wavelet transform.
In general in astronomy, it is not the image which is of interest but rather
the objects in the image! The multiresolution support provides all of the
information needed to demarcate and then label the objects.
Varying background is automatically catered for by the multiresolution
transform (equation or its equivalent in the case of a pyramidal
transform). This is a very real advantage, since alternative
approaches must devote considerable processing attention to adaptive or other
methods in order to allow for variable backgrounds or superimposed objects.
Multiresolution provides much information on the scale of the objects. Thus foreground objects (which may not be of interest, e.g., cosmic ray hits) as well as background objects (e.g., faint galaxies) can be found---assuming an appropriate wavelet (or transform kernel) and processing chain---at different resolution scales. It is not difficult to determine all objects, in this way, and then to retain only those which are of immediate interest. Examples of the use of the multiresolution approach to object detection can be found in Murtagh, Zeilinger, Starck, & Bijaoui (1994).
Objective comparison of two images is often necessary, for instance when we want to evaluate the restoration quality. Very few quantitative parameters can be extracted. The correlation between the original image and the restored one provides a classical criterion. Another way to compare two pictures is to determine the mean-square error. These criteria, however, are not sufficient. They give no information about the resulting resolution. A comprehensive criterion must take into account the resolution. We can compute for each dyadic scale the correlation coefficient and the quadratic error between the wavelet transforms of the original and the restored images. This allows us to compare, for each resolution, the restoration quality (Starck & Bijaoui 1994).
An easy way to filter an image is to keep only the significant wavelet coefficients and to reconstruct an image (Starck & Bijaoui 1994; Donoho 1992). We have shown (Starck, Bijaoui, & Murtagh 1994) that, with few iterations and by using the multiresolution support, we can have excellent results.
The object-image relation is often given by:
O is the observed object, I is the obtained image, P is the point spread function (PSF) of the imaging system, and N is an additive noise. We want to determine O knowing I and P, the main difficulties being the existence of: (1) a cut-off frequency of the PSF, and (2) the additive noise (e.g., Cornwell 1988). The multiresolution approach (Starck & Murtagh 1994) allows the noise to be controlled during the deconvolution process, and leads to a regularization of this ill-posed problem. The Poisson noise case and the use of a multiresolution support framework have been treated in (Murtagh, Starck, & Bijaoui 1994).
In interferometric imaging, measurements are carried out in the
Fourier space but the plane is not completely covered. The image
(``dirty map''), is obtained by a simple inverse Fourier transform
of the data, and the PSF (``dirty beam''),
by an inverse Fourier transform of the
plane coverage.
The presence of secondary lobes in the dirty beam creates very important
artifacts in the dirty map and a deconvolution is necessary. By applying
the CLEAN method (Högbom 1974) at each scale of the wavelet transform using
the FFT, we can localize significant structures.
An iterative reconstruction algorithm allows solutions to be found which
satisfy the positivity constraint, and the fidelity to measurements constraint
(i.e., at each measured
, we require that
the solution O satisfies
.
More details can be found in Starck & Bijaoui (1994)
and Starck, Bijaoui, Lopez, & Perrier (1994).
Several studies have used the orthogonal wavelet transform to compress astronomical data (Press 1991; White 1994). Using the multiresolution support should improve the quality of the decompressed image. Results presented in Starck, Murtagh, & Louys (1994) show how good results can be.
The wavelet transform, and more generally the multiresolution transform, provides a powerful and versatile framework for astronomical image processing. A data structure is created which allows such objectives as the following to be handled: visualization, object detection, filtering, deconvolution, and compression. Studies in all of these areas have been cited here, and in many cases the results obtained are considerably better than traditional approaches.
Anscombe, F. J. 1948, Biometrika, 15, 246
Bijaoui, A. 1991, Ondelettes et Paquets d'Ondes, Inria, Rocquencourt
Bijaoui, A., Starck, J. L., & Murtagh, F. 1994, Traitement du Signal, 3, 11
Burt, P. J., & Adelson A. E. 1983, IEEE Trans. on Communications, 31, 532
Cornwell, T. J. 1988, Proc. NATO Advanced Study Institute on Diffraction-Limited Imaging with Very Large Telescopes, Cargèse, 273
Chui, C. H. 1992, Wavelet Analysis and its Application (New York, Academic Press)
Cohen, A., Daubechies, I., & Feauveau, J. C. 1992, Comm. Pur. Appl. Math., 45, 485
Daubechies, I. 1992, Ten Lectures on Wavelets, Philadelphia.
Donoho, D. L. 1992, Proc. Progress in Wavelet Analysis and Applications (Toulouse, Ed. Frontières)
Feauveau, J. C. 1990, Thesis, University Paris Sud.
Högbom, J. A. 1974, A&AS, 15, 417
Mallat, S. 1989, IEEE Trans. on Pattern Analysis and Machine Intelligence, 11, 574
Meyer, Y. 1989, Wavelets, ed. J. M. Combes et al., (Berlin, Springer Verlag), 21
Murtagh, F., Starck, J. L., & Bijaoui, A. 1994, A&A, submitted
Murtagh, F., Zeilinger, W., , Starck, J. L., & Bijaoui, A. 1994,
Press, W. H. 1991, in Astronomical Data Analysis Software and Systems I, ASP Conf. Ser., Vol. 25, eds. D.M. Worrall, C. Biemesderfer, & J. Barnes (San Francisco, ASP), p. 25
Shensa, M. J. 1992, Proc. IEEE Transactions on Signal Processing, 40, 2464
Starck, J. L. 1993, The Wavelet Transform, MIDAS Manual (Garching, ESO)
Starck, J. L., & Bijaoui, A. 1994, Signal Processing, 35, 195
Starck, J. L., Bijaoui, A., Lopez, A., & Perrier, C. 1994, A&A, 283, 349
Starck, J. L., & Bijaoui, A. 1994, J. Opt. Soc. Am. A, 11, No. 4
Starck, J. L., & Murtagh, F. 1994, A&A, 288, 343
Starck, J. L., Bijaoui, A., & Murtagh, F., 1994, CVIP: Graphical Models and Image Processing, submitted
Starck, J. L., Murtagh, F., & Louys, M., 1994,
Van Cittert, P. H. 1931, Physik, Z., 69, 298
White, R. L. 1993, in Space and Earth Science Data Compression Workshop, ed. James C. Tilton, NASA Conference Publication 3183, p. 117