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