Multiresolution in Astronomical Image Processing ... - Semantic Scholar
Multiresolution in Astronomical Image Processing: A General Framework Fionn Murtagh (1), Jean-Luc Starck (2) and Albert Bijaoui (3) (1) Space Telescope – European Coordinating Facility, ESO, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany. (Affiliated to Astrophysics Division, Space Science Department, European Space Agency.) (2) CEA, DSM/DAPNIA, CE-Saclay, F-91191 Gif-sur-Yvette Cedex, France. (3) Observatoire de la Cote ˆ d’Azur, B.P. 229, F-06304 Nice Cedex 4, France.
Abstract Multiresolution transforms, including a wavelet transform, are applied to image visualization, image restoration, filtering and compression, and object detection. Variance stabilization is used, when appropriate, to cater for common astronomical image noise models. We discuss validation of such methods in the case of astronomical image processing. A range of examples illustrate the effectiveness of this approach in handling point source and extended astronomical objects.
1 Introduction We describe a general framework for image analysis based on the following methodological components: multiresolution; noise separation; object support image; and Minkowski operators. Image analysis applications, most of which will be touched on below, include: restoration; noise filtering; compression; artifact filtering; object detection; and object description. General characteristics of astronomical images include: the presence of noise, to be addressed in the next paragraph; point sources (stars or star-like objects); extended objects (resolved galaxies, or nebulae; also planets and other solar system objects); faint objects and extended structures, or superimposed objects. Image calibration is usual. Detector faults in individual pixels, groups of contiguous pixels, and alignments of pixels, must be catered for. Noise is often taken as additive Poisson (related to arrival of photons) and/or Gaussian. Commonly-used electronic CCD (charge-coupled device) detectors have a range of Poisson noise components, together with Gaussian read-out noise (Snyder et al., 1993). Digitized photographic images are found by Tekalp and Pavlovi´c (1991) to be also additive Poisson and Gaussian (and subject to nonlinear distortions, which we will not further discuss here).
1
2 Multiresolution An image, I0 , is successively transformed as follows:
f (I0) = I1 ; f (I1) = I2 ; f (I2) = I3 ; : : : Letting wj
=
Ij ?1 ? Ij , we can write I0 =
Xw p
j
+
wr
(1)
j =1
The images Ij are versions of the original image at p increasing scales, or decreasing resolution levels; and the image wr is a residual image. Each Ij and wj above is indexed by all pixels (x; y ). An example of such a transform is provided by the discrete wavelet transform known as the a` trous (“with holes”) algorithm (Holdschneider et al., 1989). An interlaced convolution is carried out with a B3 spline scaling function (see definitions in Starck and Murtagh, 1994): at iteration i, i > 0, the convolution makes use of a pixel distance of 2i?1 between a central pixel and adjacent ones. This scaling function is fairly similar in shape to a Gaussian; its shape is, in particular, isotropic; it is compact; it is calculated in direct space, precluding periodization artifacts; reconstruction is trivial (equation (1)); and other properties are mentioned in Starck et al. (1995). We have found this transform to be powerful, when dealing with astronomical imagery. For example, Murtagh et al. (1993) use it very effectively to illustrate faint object features (cometary coma) surrounding a brighter object (the cometary nucleus): see Figs. 1 and 2. A disadvantage of this wavelet function is related to the negative side-lobes, which can give rise to artifacts around bright point sources. An alternative multiresolution transform which we have used extensively on astronomical images is the pyramidal median transform. Consider a square kernel of dimensions n n. One possible set of dyadic resolution scales results from taking n = 2s + 1 for s = 1; 2; 4; : : : A multiresolution median transform is carried out by taking the median transform of I0, to produce I1; then increasing s to provide the next scale, take the median transform of I1 to yield I2; and so on. Multiresolution transform coefficients are formed as in equation (1). The maximum number of resolution levels is set by the user (and is inherently constrained by the dimensions of I0 ). A computational improvement in this algorithm is to note that the images I1 , I2 , : : :, contain increasing redundancy. Decimation is therefore used, i.e. only one pixel out of every 2 2 pixels is retained. The kernel remains the same. This gives rise to a pyramidal sequence of images. A problem manifests itself, however, in that the reconstruction formula defined by equation (1) is no longer valid. B-spline interpolation can be used to recreate an image at greater resolution, starting from the image of lowest resolution (i.e. of smallest dimensions). A further improvement is to iteratively improve the approximation, in going from Ij +1 to Ij (see Starck et al., 1995a).
2
Some advantages and disadvantages of wavelet and other multiresolution transforms are discussed in Starck et al., 1995a. We have found that sensitivity to certain aspects of astronomical objects, notably their isotropy, seems to be aided by the wavelet transform, especially in image restoration. The multiresolution median transform is not a wavelet transform – it is not linear (and a given scale is not guaranteed to be of zero mean) – but it allows rings around point sources produced by a wavelet transform, and which are responsible for artifacts in restoration, to be suppressed.
3 Noise Suppression We consider astronomical images to have additive noise, and the noise associated with adjacent pixels is uncorrelated. In principle noise is Poisson (due to photon arrival events), and commonly-used CCD detectors or digitized photographic plates contribute an additional Gaussian term. The Gaussian case may be an acceptable approximation of Poisson noise, and may also be an acceptable approximation of result of complex detector characteristics and calibration operations (e.g. flat-fielding, to compensate for uneveness in pixel response; geometric corrections; saturation leading to nonlinear detector response; and so on). For this reason, we have considered to date the cases of (i) Gaussian noise; (ii) Poisson noise; and (iii) Poisson plus Gaussian noise. The multiresolution planes, wj in equation (1), provide valuable high-level views of the original image, and we have already noted that they may offer a “microscope” for the study of faint structure. Therefore, we seek significant pixel values in all wj images. Reciprocally, pixel values which are not statistically significant constitute noise at that level. In the Gaussian case, wj (x; y ) is compared to kj at level j . The value of k can be set by the user, and is often chosen as 3. The value of j is defined globally, under the assumption of stationarity. If j wj j kj then wj is significant, and if j wj j is less than this threshold then wj is not significant. The planes of the a` trous transform are defined as a weighted sum of coefficients wj (x; y ) at previous scales. Hence if I0 (x; y ) is Gaussian, for all (x; y ), then wj (x; y ) is also Gaussian, for all j and (x; y ). This is not the case for the pyramidal median transform. However, with this method, we have found noise suppression based on the assumption of Gaussian noise at all levels to be heuristically acceptable and successful. We now look at the cases when the noise associated with I0 is not Gaussian. If it is Poisson, and if the mean pixel value is sufficiently large, then the Anscombe (1948) transformation provides us with the image I0 to work with:
r
3 I0 = 2 I0 + 8
I0 is Gaussian with unit standard deviation.
(2)
In the case of additive Poisson and Gaussian noise, this transformation was generalized in Murtagh et al. (1995) as follows:
r
3 I0 = 2 I0 + 2 + 2 ? 8
3
(3)
where is the gain, or number of (Poisson-distributed) electrons per data number; and and are the mean and standard deviation (in electron units) of the Gaussian noise component. Murtagh et al. (1995) discuss the conditions under which these variance stabilizing transforms are valid (broadly speaking, in regimes of sufficient signal). Noise filtering is carried out as follows. Values of wj which are insignificant, hence noisy, are set to 0, and the image is reconstructed using equation (1). This filtering is seen to be adaptive, taking account of local image conditions. Starck et al. (1995) discuss relationships with work in statistics and elsewhere which uses such thresholding in wavelet space to suppress noise. A satisfactory filtering implies that the error between the image I00 obtained by reconstructing the thresholded coefficients, and the original image I0 , i.e. E = I00 ? I0, contains only noise. Such may not be quite the case in practice (due to approximate noise assumptions, to begin with). We can improve the filtering by iterating: 1. Determine the error, E , between original image and filtered image. 2. Filter E in the same manner. Then E En are the insignificant residuals.
Es ? En, where Es is the filtered result, and
=
3. Add Es to the current estimate of the filtered image, I00 . Define a new error, and reiterate. About a half-dozen iterations usually suffice. Alternatively, we can iterate until convergence of the criterion:
j ((E
?
(n 1)
)
? (E
(n)
= (E (n)) j
))
where n is the iteration number, and is the standard deviation. Noise suppression may be effectively incorporated into image compression. The following example uses the pyramidal median transform, with noise suppression based on a standard deviation of 3.84. The original image is shown in Fig. 3, and the reconstructed image (without noise) in Fig. 4. The latter, when compressed, occupies 3.9% of the storage space required by the 16-bit version of the original. Further details of this compression method may be found in Starck et al. (1995c).
4 Support The quality of operations such as filtering or restoration on astronomical images is not measured with global criteria (e.g. square error, or robust equivalents of this). Instead, local criteria are considered of greatest importance. Prime among these are astrometric and photometric invariance. Generally, given that astronomical objects do not have privileged directions of elongation, the former should be fine if the image operations are based on isotropic convolution kernels. Photometric accuracy needs to be considered in more detail. Filtering, for instance, may very well affect object pixel values at some resolution level. Thus we “protect” such pixels related to an object in the image. A 4
strategy such as one of the following may then be used: a pixel, found as significant as some level, will be protected at all levels; or in an iterative algorithm (we have seen an iterative filtering algorithm above, and will look at iterative restoration algorithms below), a pixel found once to be significant at some level will continue to be protected in subsequent iterations. “Protection” of coefficients wj , i.e. they cannot be set to zero, is implemented with a mask image corresponding to level j . This is called the support image. Visualization of the support may use the image
X p
j
2
M (j )
j =1
where M (j ) is the support at scale j .
5 Application to Image Restoration The techniques introduced above – multiresolution image transforms, noise removal, and the multiresolution support – may be introduced into iterative image restoration. The observed image (“data”), I , is related to the deconvolved image (“image”), O, given additive noise, N , through I = OP +N (4)
where is a convolution, and P is the point spread function (PSF). The Richardson-Lucy (RL) iterative solution for O involves the coupled iterative equations
O(n+1)
? O [I=I P ] (n)
(n)
I (n)
? P O
(n)
(5)
where P is the transpose of the PSF, and the iteration is indicated by n. The residual at iteration n is
R(n) = I ? P O(n) R(n) may be multiresolution-transformed.
(6)
It may therefore, at each iteration n, be noisefiltered. This is found to considerably speed up convergence. It also provides a more smoothed image – the noise peaks are not fitted with PSFs, since they are (if all works well!) removed. This approach to image restoration is a regularized one, where the regularization is local and adaptive. Although we have discussed RL restoration, the noise suppression method is equally relevant for other restoration methods – van Cittert, fixed step gradient, etc. A simulated globular cluster (see short description in Hanisch, 1993), corresponding to the pre-refurbished Hubble Space Telescope Wide Field/Planetary Camera (WF/PC-1) image, will be used as an example. The number of point sources used in the simulation was 467, their centers being located at sub-pixel locations. Fig. 5 (upper left) shows the image following addition of noise and PSF-blurring. Fig. 5 (upper right) shows the result using multiresolution-based regularization and the generalized Anscombe (Poisson and Gaussian noise) formula. The regularized result 5
based on Gaussian noise was visually much worse than the regularized result based on Poisson noise, and the latter is shown in Fig. 5 (lower left). Fig. 5 (lower right) shows the Richardson-Lucy result with noise ignored. In these figures, the contours used were 4, 12, 20 and 28 (and note that the central regions of stars have been purposely omitted). The number of iterations used in all cases was 100. For regularization, a 3 noise threshold was applied, and 4 levels were used in the a` trous wavelet transform. Initialization of the support image also took the known object central positions into account (the practical usefulness of taking central positions, only, was limited by the spread of the objects). The “truth” locations of the stars were used to carry out aperture photometry on these results. A concentric circular aperture of 3 pixels in diameter was used. The results are shown in Fig. 6. The “Noise ignored” case corresponds to the unregularized RL restoration. A robust regression (least median of squares linear regression) was also used on the results shown in Fig. 6 to confirm that the Gaussian/Poisson noise model provided photometric results of the best quality. Further discussion of these results can be found in Murtagh et al. (1995). An investigation of object detectability was also carried out, and results are summarized in Fig. 7 (models used: “P+G”: Poisson plus Gaussian; “G”: Gaussian; “P”: Poisson; “Unreg”: unregularized). The clear distinction between the “P+G” and “P” cases, on the one hand, and the “G” and “Unreg” cases, on the other, can be seen. The results of “P+G” and “P” are relatively very similar: “P+G” is found to have a higher hit rate, at low detection thresholds, leading to both a higher number of correct objects, and a higher false alarm rate. Further details are provided in Murtagh et al. (1995).
6 Application to Object Detection The multiresolution support, discussed in section 4, provides an excellent first phase of processing for astronomical object detection and characterization. The use of the Minkowski closing operator (erosion with a suitable structuring element, followed by dilation) may allow spurious “objects” to be removed: very linear detector artifacts such as bleeding or charge-overflow effects, or cosmic ray hits on CCD images. An example of this is shown in Figs. 8, 9 and 10. NGC 4636 was discussed by Kissler et al. (1993), and characterized as a rich globular cluster system in a normal elliptical galaxy. An extracted 512 512 subimage from a European Southern Observatory New Technology Telescope (NTT) image is shown in Fig. 8. A multiresolution support of this image was obtained. This multiresolution support image with values 0, 2, 4, and 8 represents, in one image, the 4 resolution levels examined. The result of booleanizing this support image – any pixel with value greater than 2 was assigned a 1-value – and then carrying out 2 openings, is shown in Fig. 10. The object “islands” are labeled, and a report file produced, with Gaussian profile fits and other features. This information can be used to discriminate between objects of interest.
6
7 Conclusions A powerful theoretical framework has been overviewed, based on wavelet and other multiresolution image transforms. A principal objective has been to separate signal from noise. A data structure which is a by-product of this processing is the multi-level support image. The unique characteristics of astronomical image data imply that assessment of quality of results is carried out on the basis of photometric, astrometric and morphological faithfulness to the “real” data.
Acknowledgements Figs. 1 and 2 draw on work with H. Bohnhardt. ¨ Figs. 3 and 4 draw on work with B. Pirenne. Figs. 8, 9 and 10 draw on work with W. Zeilinger.
References 1. F.J. Anscombe, “The transformation of Poisson, binomial and negative-binomial data”, Biometrika 15, 246–254 (1948). 2. B. Hanisch, “WFPC simulation data sets available”, in B. Hanisch, ed., Restoration (1993), 76–77. 3. M. Holdschneider, R. Kronland-Martinet, J. Morlet and Ph. Tchamitchian, “A real-time algorithm for signal analysis with the help of the wavelet transform”, in Wavelets, J.M. Combes et al., eds, Springer-Verlag, Berlin (1989), 286–297. 4. M. Kissler, T. Richtler, E.V. Held, E.K. Grebel, S. Wagner and M. Cappaccioli, “NGC 4636 – a rich globular cluster system in a normal elliptical galaxy”, The Messenger No. 73 (1993), 32–34. 5. F. Murtagh, W.W. Zeilinger, J.-L. Starck and H. Bohnhardt, ¨ “Detection of faint extended structures by multiresolution wavelet analysis”, The Messenger No. 73 (1993), 37–39. 6. F. Murtagh, J.-L. Starck and F. Murtagh, “Image restoration with noise suppression using a multiresolution support”, Astronomy and Astrophysics (1995) in press. 7. D.L. Snyder, A.M. Hammoud and R.L. White, “Image recovery from data acquired with a charge-coupled-device camera”, Journal of the Optical Society of America 10, 1014–1023 (1993). 8. J.L. Starck and F. Murtagh, “Image restoration with noise suppression using the wavelet transform”, Astronomy and Astrophysics 288, 342–348 (1994).
7
9. J.-L. Starck, F. Murtagh and A. Bijaoui, “Multiresolution and astronomical image processing”, in Astronomical Data Analysis Software and Systems III, D. Shaw, H. Payne and J. Hayes, eds., Astron. Soc. Pac., 1995a, in press. (Preprint available at: http://http.hq.eso.org/fmurtagh/wavelets.html). 10. J.-L. Starck, F. Murtagh and A. Bijaoui, “Multiresolution support applied to image filtering and restoration”, 1995b, submitted. (Preprint available at: http://http.hq.eso.org/fmurtagh/wavelets.html). 11. J.-L. Starck, F. Murtagh and M. Louys, “Astronomical image compression using the pyramidal median transform”, in Astronomical Data Analysis Software and Systems III, D. Shaw, H. Payne and J. Hayes, eds., Astron. Soc. Pac., 1995c, in press. (Preprint available at: http://http.hq.eso.org/fmurtagh/wavelets.html). 12. A.M. Tekalp and G. Pavlovi´c, “Restoration of scanned photographic images”, in Digital Image Restoration, A.K. Katsaggelos, ed., New York, Springer-Verlag, 209–238 (1991).
8
List of Figure Captions Figure 1: Comet P/Swift-Tuttle, October-November 1992. Figure 2: Wavelet transform of Fig. 1 image, enhancing coma. Figure 3: Galaxy NGC 2683, Guide Star Scan image, original. Figure 4: Galaxy NGC 2683, Guide Star Scan image, following noise removal, compression, and uncompression. Figure 5: Top left: globular cluster, blurred and noisy image. Top right: regularized RL restoration with Poisson + Gaussian noise. Bottom left: regularized RL restoration with Poisson noise assumed. Bottom right: unregularized RL restoration. Figure 6: Linearity tests, using aperture photometry, for three different noise assumptions and for the unregularized (“noise ignored”) RL restoration. Zero-point lines are shown for reference. Figure 7: Plots of false alarms/non-detections for different noise models. Figure 8: Part of ESO NTT image of NGC 4636. Figure 9: Multiresolution support corresponding to Fig. 8, booleanized. Figure 10: Booleanized multiresolution support, following 2 Minkowski closings.
9
516 5 5
516
Figure 1: Comet P/Swift-Tuttle, October-November 1992. 10
516 5 5
516
Figure 2: Wavelet transform of Figure 1 image, enhancing coma 11
171 2 2
171
Figure 3: Galaxy NGC 2683, Guide Star Scan image, original. 12
171 2 2
171
Figure 4: Galaxy NGC 2683, Guide Star Scan image, following noise removal, compression, and uncompression. 13
Figure 5: Top left: globular cluster, blurred and noisy image. Top right: regularized RL restoration with Poisson + Gaussian noise. Bottom left: regularized RL restoration with Poisson noise assumed. Bottom right: unregularized RL restoration. 14
4 0
1
2
3
2 0 -2 -4
-1
0
1
2
0
1
Input magnitudes
2
3
3
4
4
•
• •
2
4
-1
-2
Poisson noise
2
•
•
4
-2
• •
• • • • ••• • • •• ••••• •••• •• ••••• • •• ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• •• • • • • ••• • •• •• •• •••• • ••••• • • • • • • •• • • • • • •• • • • • • • ••• • • •• • • • • •
•
-4
-4 -2
•
Gaussian noise
•• • • •• • •• •• • • • ••• •••••••••••••••••••••• •••••••••••••••••••••• ••••••••••••••••••••••••••••••••••••••• • •• • •• • • •••• ••••••• •••••••••••••••••••••• ••••••• ••• • •••• • ••••• ••• •••• ••••••••• •••• • • • •• • • • • • •• • • ••• • • •• • • • •
••
•
Input magnitudes
•
•
•
• • •• ••• •• • •• •••• ••• ••••••• •••••••• •• ••••• • •••••••••••••••••••••••••••••••••••••••••••••••••••• ••••••••••••••••••••••••••••••••••••••••••••••• ••••••••••••• • •• • • • •• • • • • ••• • •••• •• • ••• • • • • •• •• • •• •• ••• ••• • • • • •• • •• • • • • • •
Input magnitudes
•
-2
0
4
0
-1
•
-2
•
Delta magnitudes
•
Delta magnitudes
•
-2
0
2
• • • • •• • • • •• •• •• •••• •• • • • • • • •• ••••• • •• ••••••••••••••••••• • •••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• • • • •••• •••• • • • •• • • •• • • •• • •••• • • •• • •• •• • ••• • •• • •• • •• • • ••• • • • • •• • • • •
-4
Delta magnitudes
4
Noise ignored
•
-2
Delta magnitudes
15
Figure 6: Linearity tests, using aperture photometry, for three different noise assumptions and for the unregularized (“noise ignored”) RL restoration. Zero-point lines are shown for reference.
Gaussian and Poisson noise
-1
0
1
Input magnitudes
2
3
4
1.2 1.0 0.8 0.6 0.4
G
0.2
P
P+G Unreg
0.0
Number of false alarms / Number of non-detections
1.4
Object detection confusion plots
0
5
10
15
20
25
Detection threshold
Figure 7: Plots of false alarms/non-detections for different noise models.
16
30
1748 950 400
1198
Figure 8: Part of ESO NTT image of NGC 4636. 17
799 1 1
799
Figure 9: Multiresolution support corresponding to Fig. 8, booleanized. 18
799 1 1
799
Figure 10: Booleanized multiresolution support, following 2 Minkowski closings. 19
Abstract Multiresolution transforms, including a wavelet transform, are applied to image visualization, image restoration, filtering and compression, and object detection. Variance stabilization is used, when appropriate, to cater for common astronomical image noise models. We discuss validation of such methods in the case of astronomical image processing. A range of examples illustrate the effectiveness of this approach in handling point source and extended astronomical objects.
1 Introduction We describe a general framework for image analysis based on the following methodological components: multiresolution; noise separation; object support image; and Minkowski operators. Image analysis applications, most of which will be touched on below, include: restoration; noise filtering; compression; artifact filtering; object detection; and object description. General characteristics of astronomical images include: the presence of noise, to be addressed in the next paragraph; point sources (stars or star-like objects); extended objects (resolved galaxies, or nebulae; also planets and other solar system objects); faint objects and extended structures, or superimposed objects. Image calibration is usual. Detector faults in individual pixels, groups of contiguous pixels, and alignments of pixels, must be catered for. Noise is often taken as additive Poisson (related to arrival of photons) and/or Gaussian. Commonly-used electronic CCD (charge-coupled device) detectors have a range of Poisson noise components, together with Gaussian read-out noise (Snyder et al., 1993). Digitized photographic images are found by Tekalp and Pavlovi´c (1991) to be also additive Poisson and Gaussian (and subject to nonlinear distortions, which we will not further discuss here).
1
2 Multiresolution An image, I0 , is successively transformed as follows:
f (I0) = I1 ; f (I1) = I2 ; f (I2) = I3 ; : : : Letting wj
=
Ij ?1 ? Ij , we can write I0 =
Xw p
j
+
wr
(1)
j =1
The images Ij are versions of the original image at p increasing scales, or decreasing resolution levels; and the image wr is a residual image. Each Ij and wj above is indexed by all pixels (x; y ). An example of such a transform is provided by the discrete wavelet transform known as the a` trous (“with holes”) algorithm (Holdschneider et al., 1989). An interlaced convolution is carried out with a B3 spline scaling function (see definitions in Starck and Murtagh, 1994): at iteration i, i > 0, the convolution makes use of a pixel distance of 2i?1 between a central pixel and adjacent ones. This scaling function is fairly similar in shape to a Gaussian; its shape is, in particular, isotropic; it is compact; it is calculated in direct space, precluding periodization artifacts; reconstruction is trivial (equation (1)); and other properties are mentioned in Starck et al. (1995). We have found this transform to be powerful, when dealing with astronomical imagery. For example, Murtagh et al. (1993) use it very effectively to illustrate faint object features (cometary coma) surrounding a brighter object (the cometary nucleus): see Figs. 1 and 2. A disadvantage of this wavelet function is related to the negative side-lobes, which can give rise to artifacts around bright point sources. An alternative multiresolution transform which we have used extensively on astronomical images is the pyramidal median transform. Consider a square kernel of dimensions n n. One possible set of dyadic resolution scales results from taking n = 2s + 1 for s = 1; 2; 4; : : : A multiresolution median transform is carried out by taking the median transform of I0, to produce I1; then increasing s to provide the next scale, take the median transform of I1 to yield I2; and so on. Multiresolution transform coefficients are formed as in equation (1). The maximum number of resolution levels is set by the user (and is inherently constrained by the dimensions of I0 ). A computational improvement in this algorithm is to note that the images I1 , I2 , : : :, contain increasing redundancy. Decimation is therefore used, i.e. only one pixel out of every 2 2 pixels is retained. The kernel remains the same. This gives rise to a pyramidal sequence of images. A problem manifests itself, however, in that the reconstruction formula defined by equation (1) is no longer valid. B-spline interpolation can be used to recreate an image at greater resolution, starting from the image of lowest resolution (i.e. of smallest dimensions). A further improvement is to iteratively improve the approximation, in going from Ij +1 to Ij (see Starck et al., 1995a).
2
Some advantages and disadvantages of wavelet and other multiresolution transforms are discussed in Starck et al., 1995a. We have found that sensitivity to certain aspects of astronomical objects, notably their isotropy, seems to be aided by the wavelet transform, especially in image restoration. The multiresolution median transform is not a wavelet transform – it is not linear (and a given scale is not guaranteed to be of zero mean) – but it allows rings around point sources produced by a wavelet transform, and which are responsible for artifacts in restoration, to be suppressed.
3 Noise Suppression We consider astronomical images to have additive noise, and the noise associated with adjacent pixels is uncorrelated. In principle noise is Poisson (due to photon arrival events), and commonly-used CCD detectors or digitized photographic plates contribute an additional Gaussian term. The Gaussian case may be an acceptable approximation of Poisson noise, and may also be an acceptable approximation of result of complex detector characteristics and calibration operations (e.g. flat-fielding, to compensate for uneveness in pixel response; geometric corrections; saturation leading to nonlinear detector response; and so on). For this reason, we have considered to date the cases of (i) Gaussian noise; (ii) Poisson noise; and (iii) Poisson plus Gaussian noise. The multiresolution planes, wj in equation (1), provide valuable high-level views of the original image, and we have already noted that they may offer a “microscope” for the study of faint structure. Therefore, we seek significant pixel values in all wj images. Reciprocally, pixel values which are not statistically significant constitute noise at that level. In the Gaussian case, wj (x; y ) is compared to kj at level j . The value of k can be set by the user, and is often chosen as 3. The value of j is defined globally, under the assumption of stationarity. If j wj j kj then wj is significant, and if j wj j is less than this threshold then wj is not significant. The planes of the a` trous transform are defined as a weighted sum of coefficients wj (x; y ) at previous scales. Hence if I0 (x; y ) is Gaussian, for all (x; y ), then wj (x; y ) is also Gaussian, for all j and (x; y ). This is not the case for the pyramidal median transform. However, with this method, we have found noise suppression based on the assumption of Gaussian noise at all levels to be heuristically acceptable and successful. We now look at the cases when the noise associated with I0 is not Gaussian. If it is Poisson, and if the mean pixel value is sufficiently large, then the Anscombe (1948) transformation provides us with the image I0 to work with:
r
3 I0 = 2 I0 + 8
I0 is Gaussian with unit standard deviation.
(2)
In the case of additive Poisson and Gaussian noise, this transformation was generalized in Murtagh et al. (1995) as follows:
r
3 I0 = 2 I0 + 2 + 2 ? 8
3
(3)
where is the gain, or number of (Poisson-distributed) electrons per data number; and and are the mean and standard deviation (in electron units) of the Gaussian noise component. Murtagh et al. (1995) discuss the conditions under which these variance stabilizing transforms are valid (broadly speaking, in regimes of sufficient signal). Noise filtering is carried out as follows. Values of wj which are insignificant, hence noisy, are set to 0, and the image is reconstructed using equation (1). This filtering is seen to be adaptive, taking account of local image conditions. Starck et al. (1995) discuss relationships with work in statistics and elsewhere which uses such thresholding in wavelet space to suppress noise. A satisfactory filtering implies that the error between the image I00 obtained by reconstructing the thresholded coefficients, and the original image I0 , i.e. E = I00 ? I0, contains only noise. Such may not be quite the case in practice (due to approximate noise assumptions, to begin with). We can improve the filtering by iterating: 1. Determine the error, E , between original image and filtered image. 2. Filter E in the same manner. Then E En are the insignificant residuals.
Es ? En, where Es is the filtered result, and
=
3. Add Es to the current estimate of the filtered image, I00 . Define a new error, and reiterate. About a half-dozen iterations usually suffice. Alternatively, we can iterate until convergence of the criterion:
j ((E
?
(n 1)
)
? (E
(n)
= (E (n)) j
))
where n is the iteration number, and is the standard deviation. Noise suppression may be effectively incorporated into image compression. The following example uses the pyramidal median transform, with noise suppression based on a standard deviation of 3.84. The original image is shown in Fig. 3, and the reconstructed image (without noise) in Fig. 4. The latter, when compressed, occupies 3.9% of the storage space required by the 16-bit version of the original. Further details of this compression method may be found in Starck et al. (1995c).
4 Support The quality of operations such as filtering or restoration on astronomical images is not measured with global criteria (e.g. square error, or robust equivalents of this). Instead, local criteria are considered of greatest importance. Prime among these are astrometric and photometric invariance. Generally, given that astronomical objects do not have privileged directions of elongation, the former should be fine if the image operations are based on isotropic convolution kernels. Photometric accuracy needs to be considered in more detail. Filtering, for instance, may very well affect object pixel values at some resolution level. Thus we “protect” such pixels related to an object in the image. A 4
strategy such as one of the following may then be used: a pixel, found as significant as some level, will be protected at all levels; or in an iterative algorithm (we have seen an iterative filtering algorithm above, and will look at iterative restoration algorithms below), a pixel found once to be significant at some level will continue to be protected in subsequent iterations. “Protection” of coefficients wj , i.e. they cannot be set to zero, is implemented with a mask image corresponding to level j . This is called the support image. Visualization of the support may use the image
X p
j
2
M (j )
j =1
where M (j ) is the support at scale j .
5 Application to Image Restoration The techniques introduced above – multiresolution image transforms, noise removal, and the multiresolution support – may be introduced into iterative image restoration. The observed image (“data”), I , is related to the deconvolved image (“image”), O, given additive noise, N , through I = OP +N (4)
where is a convolution, and P is the point spread function (PSF). The Richardson-Lucy (RL) iterative solution for O involves the coupled iterative equations
O(n+1)
? O [I=I P ] (n)
(n)
I (n)
? P O
(n)
(5)
where P is the transpose of the PSF, and the iteration is indicated by n. The residual at iteration n is
R(n) = I ? P O(n) R(n) may be multiresolution-transformed.
(6)
It may therefore, at each iteration n, be noisefiltered. This is found to considerably speed up convergence. It also provides a more smoothed image – the noise peaks are not fitted with PSFs, since they are (if all works well!) removed. This approach to image restoration is a regularized one, where the regularization is local and adaptive. Although we have discussed RL restoration, the noise suppression method is equally relevant for other restoration methods – van Cittert, fixed step gradient, etc. A simulated globular cluster (see short description in Hanisch, 1993), corresponding to the pre-refurbished Hubble Space Telescope Wide Field/Planetary Camera (WF/PC-1) image, will be used as an example. The number of point sources used in the simulation was 467, their centers being located at sub-pixel locations. Fig. 5 (upper left) shows the image following addition of noise and PSF-blurring. Fig. 5 (upper right) shows the result using multiresolution-based regularization and the generalized Anscombe (Poisson and Gaussian noise) formula. The regularized result 5
based on Gaussian noise was visually much worse than the regularized result based on Poisson noise, and the latter is shown in Fig. 5 (lower left). Fig. 5 (lower right) shows the Richardson-Lucy result with noise ignored. In these figures, the contours used were 4, 12, 20 and 28 (and note that the central regions of stars have been purposely omitted). The number of iterations used in all cases was 100. For regularization, a 3 noise threshold was applied, and 4 levels were used in the a` trous wavelet transform. Initialization of the support image also took the known object central positions into account (the practical usefulness of taking central positions, only, was limited by the spread of the objects). The “truth” locations of the stars were used to carry out aperture photometry on these results. A concentric circular aperture of 3 pixels in diameter was used. The results are shown in Fig. 6. The “Noise ignored” case corresponds to the unregularized RL restoration. A robust regression (least median of squares linear regression) was also used on the results shown in Fig. 6 to confirm that the Gaussian/Poisson noise model provided photometric results of the best quality. Further discussion of these results can be found in Murtagh et al. (1995). An investigation of object detectability was also carried out, and results are summarized in Fig. 7 (models used: “P+G”: Poisson plus Gaussian; “G”: Gaussian; “P”: Poisson; “Unreg”: unregularized). The clear distinction between the “P+G” and “P” cases, on the one hand, and the “G” and “Unreg” cases, on the other, can be seen. The results of “P+G” and “P” are relatively very similar: “P+G” is found to have a higher hit rate, at low detection thresholds, leading to both a higher number of correct objects, and a higher false alarm rate. Further details are provided in Murtagh et al. (1995).
6 Application to Object Detection The multiresolution support, discussed in section 4, provides an excellent first phase of processing for astronomical object detection and characterization. The use of the Minkowski closing operator (erosion with a suitable structuring element, followed by dilation) may allow spurious “objects” to be removed: very linear detector artifacts such as bleeding or charge-overflow effects, or cosmic ray hits on CCD images. An example of this is shown in Figs. 8, 9 and 10. NGC 4636 was discussed by Kissler et al. (1993), and characterized as a rich globular cluster system in a normal elliptical galaxy. An extracted 512 512 subimage from a European Southern Observatory New Technology Telescope (NTT) image is shown in Fig. 8. A multiresolution support of this image was obtained. This multiresolution support image with values 0, 2, 4, and 8 represents, in one image, the 4 resolution levels examined. The result of booleanizing this support image – any pixel with value greater than 2 was assigned a 1-value – and then carrying out 2 openings, is shown in Fig. 10. The object “islands” are labeled, and a report file produced, with Gaussian profile fits and other features. This information can be used to discriminate between objects of interest.
6
7 Conclusions A powerful theoretical framework has been overviewed, based on wavelet and other multiresolution image transforms. A principal objective has been to separate signal from noise. A data structure which is a by-product of this processing is the multi-level support image. The unique characteristics of astronomical image data imply that assessment of quality of results is carried out on the basis of photometric, astrometric and morphological faithfulness to the “real” data.
Acknowledgements Figs. 1 and 2 draw on work with H. Bohnhardt. ¨ Figs. 3 and 4 draw on work with B. Pirenne. Figs. 8, 9 and 10 draw on work with W. Zeilinger.
References 1. F.J. Anscombe, “The transformation of Poisson, binomial and negative-binomial data”, Biometrika 15, 246–254 (1948). 2. B. Hanisch, “WFPC simulation data sets available”, in B. Hanisch, ed., Restoration (1993), 76–77. 3. M. Holdschneider, R. Kronland-Martinet, J. Morlet and Ph. Tchamitchian, “A real-time algorithm for signal analysis with the help of the wavelet transform”, in Wavelets, J.M. Combes et al., eds, Springer-Verlag, Berlin (1989), 286–297. 4. M. Kissler, T. Richtler, E.V. Held, E.K. Grebel, S. Wagner and M. Cappaccioli, “NGC 4636 – a rich globular cluster system in a normal elliptical galaxy”, The Messenger No. 73 (1993), 32–34. 5. F. Murtagh, W.W. Zeilinger, J.-L. Starck and H. Bohnhardt, ¨ “Detection of faint extended structures by multiresolution wavelet analysis”, The Messenger No. 73 (1993), 37–39. 6. F. Murtagh, J.-L. Starck and F. Murtagh, “Image restoration with noise suppression using a multiresolution support”, Astronomy and Astrophysics (1995) in press. 7. D.L. Snyder, A.M. Hammoud and R.L. White, “Image recovery from data acquired with a charge-coupled-device camera”, Journal of the Optical Society of America 10, 1014–1023 (1993). 8. J.L. Starck and F. Murtagh, “Image restoration with noise suppression using the wavelet transform”, Astronomy and Astrophysics 288, 342–348 (1994).
7
9. J.-L. Starck, F. Murtagh and A. Bijaoui, “Multiresolution and astronomical image processing”, in Astronomical Data Analysis Software and Systems III, D. Shaw, H. Payne and J. Hayes, eds., Astron. Soc. Pac., 1995a, in press. (Preprint available at: http://http.hq.eso.org/fmurtagh/wavelets.html). 10. J.-L. Starck, F. Murtagh and A. Bijaoui, “Multiresolution support applied to image filtering and restoration”, 1995b, submitted. (Preprint available at: http://http.hq.eso.org/fmurtagh/wavelets.html). 11. J.-L. Starck, F. Murtagh and M. Louys, “Astronomical image compression using the pyramidal median transform”, in Astronomical Data Analysis Software and Systems III, D. Shaw, H. Payne and J. Hayes, eds., Astron. Soc. Pac., 1995c, in press. (Preprint available at: http://http.hq.eso.org/fmurtagh/wavelets.html). 12. A.M. Tekalp and G. Pavlovi´c, “Restoration of scanned photographic images”, in Digital Image Restoration, A.K. Katsaggelos, ed., New York, Springer-Verlag, 209–238 (1991).
8
List of Figure Captions Figure 1: Comet P/Swift-Tuttle, October-November 1992. Figure 2: Wavelet transform of Fig. 1 image, enhancing coma. Figure 3: Galaxy NGC 2683, Guide Star Scan image, original. Figure 4: Galaxy NGC 2683, Guide Star Scan image, following noise removal, compression, and uncompression. Figure 5: Top left: globular cluster, blurred and noisy image. Top right: regularized RL restoration with Poisson + Gaussian noise. Bottom left: regularized RL restoration with Poisson noise assumed. Bottom right: unregularized RL restoration. Figure 6: Linearity tests, using aperture photometry, for three different noise assumptions and for the unregularized (“noise ignored”) RL restoration. Zero-point lines are shown for reference. Figure 7: Plots of false alarms/non-detections for different noise models. Figure 8: Part of ESO NTT image of NGC 4636. Figure 9: Multiresolution support corresponding to Fig. 8, booleanized. Figure 10: Booleanized multiresolution support, following 2 Minkowski closings.
9
516 5 5
516
Figure 1: Comet P/Swift-Tuttle, October-November 1992. 10
516 5 5
516
Figure 2: Wavelet transform of Figure 1 image, enhancing coma 11
171 2 2
171
Figure 3: Galaxy NGC 2683, Guide Star Scan image, original. 12
171 2 2
171
Figure 4: Galaxy NGC 2683, Guide Star Scan image, following noise removal, compression, and uncompression. 13
Figure 5: Top left: globular cluster, blurred and noisy image. Top right: regularized RL restoration with Poisson + Gaussian noise. Bottom left: regularized RL restoration with Poisson noise assumed. Bottom right: unregularized RL restoration. 14
4 0
1
2
3
2 0 -2 -4
-1
0
1
2
0
1
Input magnitudes
2
3
3
4
4
•
• •
2
4
-1
-2
Poisson noise
2
•
•
4
-2
• •
• • • • ••• • • •• ••••• •••• •• ••••• • •• ••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• •• • • • • ••• • •• •• •• •••• • ••••• • • • • • • •• • • • • • •• • • • • • • ••• • • •• • • • • •
•
-4
-4 -2
•
Gaussian noise
•• • • •• • •• •• • • • ••• •••••••••••••••••••••• •••••••••••••••••••••• ••••••••••••••••••••••••••••••••••••••• • •• • •• • • •••• ••••••• •••••••••••••••••••••• ••••••• ••• • •••• • ••••• ••• •••• ••••••••• •••• • • • •• • • • • • •• • • ••• • • •• • • • •
••
•
Input magnitudes
•
•
•
• • •• ••• •• • •• •••• ••• ••••••• •••••••• •• ••••• • •••••••••••••••••••••••••••••••••••••••••••••••••••• ••••••••••••••••••••••••••••••••••••••••••••••• ••••••••••••• • •• • • • •• • • • • ••• • •••• •• • ••• • • • • •• •• • •• •• ••• ••• • • • • •• • •• • • • • • •
Input magnitudes
•
-2
0
4
0
-1
•
-2
•
Delta magnitudes
•
Delta magnitudes
•
-2
0
2
• • • • •• • • • •• •• •• •••• •• • • • • • • •• ••••• • •• ••••••••••••••••••• • •••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• • • • •••• •••• • • • •• • • •• • • •• • •••• • • •• • •• •• • ••• • •• • •• • •• • • ••• • • • • •• • • • •
-4
Delta magnitudes
4
Noise ignored
•
-2
Delta magnitudes
15
Figure 6: Linearity tests, using aperture photometry, for three different noise assumptions and for the unregularized (“noise ignored”) RL restoration. Zero-point lines are shown for reference.
Gaussian and Poisson noise
-1
0
1
Input magnitudes
2
3
4
1.2 1.0 0.8 0.6 0.4
G
0.2
P
P+G Unreg
0.0
Number of false alarms / Number of non-detections
1.4
Object detection confusion plots
0
5
10
15
20
25
Detection threshold
Figure 7: Plots of false alarms/non-detections for different noise models.
16
30
1748 950 400
1198
Figure 8: Part of ESO NTT image of NGC 4636. 17
799 1 1
799
Figure 9: Multiresolution support corresponding to Fig. 8, booleanized. 18
799 1 1
799
Figure 10: Booleanized multiresolution support, following 2 Minkowski closings. 19
