Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://hea-www.harvard.edu/astrostat/Stat310_1112/2002ApJS..138..185F.pdf Äàòà èçìåíåíèÿ: Tue Feb 7 23:24:03 2012 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 06:10:54 2012 Êîäèðîâêà: koi8-r

THE ASTROPHYSICAL JOURNAL SUPPLEMENT SERIES, 138 : 185 õ 218, 2002 January
( 2002. The American Astronomical Society. All rights reserved. Printed in U.S.A.

A WAVELET-BASED ALGORITHM FOR THE SPATIAL ANALYSIS OF POISSON DATA P. E. FREEMAN,1 V. KASHYAP,1 R. ROSNER,2 AND D. Q. LAMB2
Received 1999 April 9 ; accepted 2001 August 27

ABSTRACT Wavelets are scalable, oscillatory functions that deviate from zero only within a limited spatial regime and have average value zero, and thus may be used to simultaneously characterize the shape, location, and strength of astronomical sources. But in addition to their use as source characterizers, wavelet functions are rapidly gaining currency within the source detection ïeld. Wavelet-based source detection involves the correlation of scaled wavelet functions with binned, two-dimensional image data. If the chosen wavelet function exhibits the property of vanishing moments, signiïcantly nonzero correlation coefficients will be observed only where there are high-order variations in the data ; e.g., they will be observed in the vicinity of sources. Source pixels are identiïed by comparing each correlation coefficient with its probability sampling distribution, which is a function of the (estimated or a priori known) background amplitude. In this paper, we describe the mission-independent, wavelet-based source detection algorithm "" WAVDETECT,îî part of the freely available Chandra Interactive Analysis of Observations (CIAO) software package. Our algorithm uses the Marr, or "" Mexican Hat îî wavelet function, but may be adapted for use with other wavelet functions. Aspects of our algorithm include : (1) the computation of local, exposure-corrected normalized (i.e., ÿat-ïelded) background maps ; (2) the correction for exposure variations within the ïeld of view (due to, e.g., telescope support ribs or the edge of the ïeld) ; (3) its applicability within the low-counts regime, as it does not require a minimum number of background counts per pixel for the accurate computation of source detection thresholds ; (4) the generation of a source list in a manner that does not depend upon a detailed knowledge of the point spread function (PSF) shape ; and (5) error analysis. These features make our algorithm considerably more general than previous methods developed for the analysis of X-ray image data, especially in the low count regime. We demonstrate the robustness of WAVDETECT by applying it to an image from an idealized detector with a spatially invariant Gaussian PSF and an exposure map similar to that of the Einstein IPC ; to Pleiades Cluster data collected by the ROSAT PSPC ; and to simulated Chandra ACIS-I image of the Lockman Hole region. Subject headings : methods : data analysis õ techniques : image processing õ X-rays : general
1

. INTRODUCTION

The detection and characterization of astronomical sources becomes increasingly difficult as we attempt to observe these sources in the EUV, X-ray, and gamma-ray spectral regimes. There are several reasons for this. First, in these high-energy regimes, source data may consist of only a few counts, so that we must rely on the Poisson distribution when making statistical inferences rather than using Gaussian statistics that are considerably easier to apply, but that are strictly applicable only in the high-counts limit. Second, spatially extended sources, such as supernova remnants and galaxy clusters, exhibit bright diuse emission at high energies which may overlap with point sources, rendering the latter more difficult both to detect and characterize. And third, the present generation of broadband high-energy telescopes, unlike optical telescopes, have spatially nonuniform point spread functions (PSFs) as an unavoidable byproduct of their design. For instance, the PSF of the Wolter I-type High Resolution Mirror Assembly (HRMA) on the Chandra X-Ray Observatory (CXO) has a 50% encircled energy radius that varies in width from B0A on-axis to .3
1 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138. 2 Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637.

Z10A near the outer edges of an ACIS-I chip (B10@ oaxis). A standard method for the analysis of Poisson count data involves the application of the so-called "" sliding cell îî (Harnden et al. 1984).3 In sliding-cell analysis, two coaligned but dierently sized square cells are placed at each image pixel, with the number of counts in the annular region between cells providing an estimate of the local background amplitude at the pixel. This amplitude is then used to compute the Poisson signiïcance of the observed number of counts in the inner cell. The usefulness of a sliding-cell algorithm is limited both when the ïeld of view (FOV) is crowded, since overlapping sources cannot be handled and/or nearby sources contributing counts to the estimated background may decrease detection sensitivity, and when sources are observed o-axis, since uncertainties in the model of the PSF, which generally increase with oaxis angle, can greatly aect source property estimates (see, e.g., Kashyap et al. 1994). The sliding cell also may not provide accurate source property estimates for extended sources. Source characterization can be improved by ïtting the detected sources using maximum likelihood methods (as in the algorithm of Hasinger, Johnston, & Verbunt 1994),

3 The method used, e.g., in the CIAO source detection routine CELLDETECT.

185


186

FREEMAN ET AL.

Vol. 138

but the accuracy of this method is still limited by uncertainties in the PSF. Within the last decade, astronomers have begun to apply wavelet functions to the problem of source detection. (For an introduction to the theory of wavelet functions, see, e.g., Mallat 1998 and references therein.) These functions are scalable, oscillatory, have a ïnite support (i.e., are nonzero within a limited spatial regime), and have an average value of zero ; they can be used to deïne a set of basis functions that act as highly localized ïlters in both the spatial and frequency domains, and thus are superior source characterizers. Certain wavelet functions also exhibit the property of vanishing moments, which is important for source detection : the integral of the product of a wavelet with N vanishing moments and a polynomial of degree ¹N [ 1 is zero. Thus, the correlation of a suitably chosen wavelet function with a photon counts image will yield correlation coefficients which are signiïcantly large only in the vicinity of localized high-order variations in the data, e.g., in the vicinity of astronomical sources, which appear as PSFbroadened bumps with inïnite power-series representation. Wavelet-based source detection thus boils down to the statistical problem of identifying image pixels with "" signiïcantly large îî correlation coefficients. Damiani et al. (1997a) were the ïrst to present a waveletbased generalized method for source detection and characterization, which they apply to ROSAT PSPC data in a subsequent work (Damiani et al. 1997b). Among its features, this method, unlike the others, uses exposure maps to handle situations in which, e.g., support rib shadows or the edge of the FOV lie within the wavelet support, allowing the analysis of the entire FOV. Their algorithm is, however, instrument-dependent in that they must make signiïcant changes to it in order to account for qualitative dierences between exposure maps of dierent detectors (e.g., between the exposure maps of the ROSAT PSPC and ROSAT HRI ; Micela et al. 1999 provide a short description of the changes that are necessary to perform HRI image analysis). In addition, while one can apply the Damiani et al. method to data from other photon-counting detectors, the PSFs must be very nearly Gaussian. Last, they publish a function for the computation of source detection thresholds which is strictly valid only when the mean background amplitude per image pixel is Z(0.1 counts pixel~1)/p2, where p is the wavelet scale size. In this paper (and in Freeman et al. 2002) we describe our own algorithm for wavelet source detection and characterization that has been developed for a generic detector and implemented in the Chandra Interactive Analysis of Observations (CIAO) routine WAVDETECT.4 In subsequent papers (e.g., Kashyap et al. 2002, in preparation) we will discuss the application of this algorithm to speciïc scientiïc problems. Our algorithm is considerably more general and ÿexible than others that have been developed, in that it can : (1) operate eectively in the low background counts regime, which is crucial because of the expected low particle and cosmic background count rates for the Chandra detectors (the overall rate being D10~6 and 10~7 counts s~1 pixel~1

for the Chandra ACIS and HRC detectors, respectively) ; and (2) operate eectively regardless of the PSF shape, also crucial because of the (non-Gaussian) nature of the o-axis Chandra PSFs.5 It also (3) corrects for the eect of exposure variations in a general, nonõdetector-speciïc manner. Thus, our algorithm may be immediately adapted for the analysis of data from virtually any other photon-counting detector. In ° 2, we provide a brief description of wavelet functions, and deïne the Marr, or "" Mexican Hat îî (MH) wavelet function which we use in our algorithm.6 We then present a simple example in which we apply the MH function to idealized data, in order to build the readerîs intuition about how to interpret the results. The MH function has been used often in astronomical wavelet analyses,7 in, for example : the analysis of 13CO spectral maps of the molecular cloud L1551 (Gill & Henrikson 1990 ; see also Langer, Wilson, & Anderson 1993, who apply Laplacian pyramid transforms to a 13CO image map of Barnard 5) ; the analysis of galaxy cluster structure in optically derived catalogs (Slezak, Bijaoui, & Mars 1990) ; the detection and localization of features in optical CCD images of galaxies (Coupinot et al. 1992) ; the analysis of substructures observed in ROSAT PSPC images of the Coma Cluster (Vikhlinin, Forman, & Jones 1994) ; the examination of Einstein HRI and ROSAT PSPC images of Abell 1367 (Grebenev et al. 1995) ; the detection of serendipitous X-ray clusters in archival ROSAT PSPC data (SHARC ; cf. Ulmer et al. 1995 ; Freeman et al. 1996 ; Nichol et al. 1997) ; and the analysis and modeling of X-ray emission in Abell clusters (Lazzati & Chincarini 1998 ; see also Lazzati et al. 1998) and stellar clusters (Damiani et al. 1997b). In ° 3 we describe the basic steps of our algorithm, in which sources are detected and characterized. The basic steps of source detection include : the correlation of the wavelet function with the data to create the correlation image ; the estimation of the local background (if necessary) in each pixel ; the computation of the source detection threshold in each pixel ; accounting for the eects of exposure variations within the FOV on the correlation value ; and the identiïcation and "" cleansing îî of source counts from the image. Source count cleansing is done iteratively, ïrst by analyzing the raw data, then the ïrst version of the cleansed data, etc., until the background is estimated. This ïnal background is used to compute source detection

4 The CIAO software http ://cxc.harvard.edu/ciao/. "" WTRANSFORM,îî a source generator ; these programs may

package may be downloaded from WAVDETECT is composed of detector, and "" WRECON,îî a source list be run separately.

5 While it can operate to a limited extent if nothing at all is known about the PSF, our algorithm is most eective if characteristic PSF sizes, e.g., the radii of circles containing 50% of the encircled energy for dierent o-axis angles, are computable. 6 We note that while our algorithm uses the MH function exclusively, one can adapt our algorithm to work with other wavelet functions. 7 However, we must note that the Mexican Hat is not the only wavelet function used by astronomers ; for instance, Rosati et al. (1995) and Vikhlinin et al. (1998) use the Morlet wavelet function to detect and analyze X-ray clusters in ROSAT PSPC images, while Slezak, Durret, & Gerbal (1994), Biviano et al. (1996), and Pierre & Starck (1998) use the so-called ` "" B-splines îî in their cluster analyses. Also, there is the a trous algorithm of Starck, Murtagh, & Bijaoui (1995), Starck & Murtagh (1998), and Starck & Pierre (1998). Methods based on this algorithm are similar to the one which we describe in this paper ; however, because their use is limited to particular problems, these methods are not sufficiently generalized to be applied to the full-FOV data of an arbitrary detector. For instance, these methods generally do not take into account exposure variations within the FOV or are applied only in limited regions not greatly aected by vignetting, and they also generally ignore local variations in the background.


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

187

thresholds, which are subsequently used to identify putative sources. Source characterization involves the combination of outputs from a number of wavelet scales, and the creation of source cells on the initial image within which source properties (positions, count rates, etc.) are estimated. In ° 4 we demonstrate the efficacy of the algorithm by applying it to a large variety of cases, such as : an idealized image containing simulated extended and point sources, with spatially invariant Gaussian PSF and an exposure map similar to that of the Einstein IPC ; a 32 ks ROSAT PSPC observation of the Pleiades cluster (cf. Micela et al. 1996) ; and a simulated 30 ks Chandra ACIS-I observation of the Lockman Hole region (T. Gaetz 1998, private communication). We then describe the dierences between our method and previously published methods, especially that of Damiani et al., in ° 5, and provide our summary and conclusions in ° 6. In future works, our algorithm will be applied to the analysis of speciïc scientiïc problems, e.g., observed spatial variations in the diuse X-ray background of the Pleiades Cluster (see Kashyap et al. 1996 and Kashyap et al. 2002, in preparation).
2

from which is created the Marr, or "" Mexican Hat îî (MH), wavelet function : 1 x2 y2 xy \ 2[ [ , W@ mp p 2np p p2 p2 xy x y xy 2)~(y2@2p2) ]e~(x2@2px y \

A

B

A

B

. WAVELETS AND SOURCE DETECTION

2.1. W avelet Properties Wavelets can be used to ïlter an image at a given length scale.8 This can be seen by considering the dierence between two smoothed functions : one formed by convolving an (arbitrary) function f with a real-valued, nonnegative, inïnitely dierentiable smoothing function / ½ C=, whose size is characterized by a scale p and satisïes the condition / / \ 1 (e.g., the Gaussian function), and one formed by the convolution of f with the same smoothing function with scale size p ] dp (Holschneider 1995). Because all structure at length scales smaller than the smoothing scales would be suppressed, the dierence of these two smoothed functions will provide information about the details of f that are introduced at the scale p itself. In the situation of interest to us, the analysis of twodimensional images, the relationship between / and the wavelet function W @ is, in the limits dp , dp ] 0, x y xy xy L L , , \[ p ]p / . (1) W@ x Lp y Lp mp p pp xy x y xy (The reason for our use of a prime symbol will become apparent below.) This is an analyzing wavelet, or mother wavelet (hence the subscript m). Other members of the same wavelet family (so-called "" atoms îî) can be generated from this mother wavelet via dilations and translations :

xy 1 , . (5) W 2np p m p p xy xy This is the wavelet function we use in our source detection algorithm. (Note that we use the function W instead of W @ , m m to be consistent with Damiani et al. 1997a.) The MH function has a positive kernel with shape similar to a canonical PSF (PW ), surrounded by a negative annulus (NW ; Fig. 1). Ellipses with semiaxes of length J2p , 2p , and 5p x y x y x y describe the boundary between PW and NW , the minimum of the MH function, and the eective support of the MH function, respectively. The MH function oers several advantages which motivate its use for source detection : (1) it has two vanishing moments (i.e., the correlation of the MH function with constant or linear functions is zero), so that it acts to suppress the contribution of the (generally spatially constant) background to the correlation coefficients ; (2) while it does not have compact support, and cannot be used to construct a set of orthogonal basis functions, a dyadic sequence of MH functions (i.e., with MH functions with scales separated by factors of 2) is sufficient to sample the entire frequency domain, because of the limited extent of an MH functionîs Fourier transform ; (3) its limited extent in both spatial and Fourier domains helps to minimize eects of aliasing ; and (4) it is analytically manipulable, so that many numerical operations may be performed using analytically derivable functions (see Appendix A), which can reduce computation time signiïcantly. 2.2. A Simple Example Before describing our algorithm in detail, we present a simpliïed example to help build the readerîs intuition. We assume that we are analyzing a subset of an evenly exposed image with the MH wavelet function, far from any edge of the FOV. Within this image, the counts in each pixel, D , i, j are sampled from a function which has a constant, relatively large (?1) amplitude, which we denote B. These assumptions allow us to build intuition with a minimum of unnecessary detail (such as correcting for exposure variations, etc.), and should not be construed as being reÿective of limitations on the applicability of our algorithm. As will be described in ° 3, the ïrst step in our algorithm is to correlate the wavelet function W and the binned image data D. We denote correlation using the symbol S ... % ... T ; in this case, we would write C \ SW % DT, or for a particular pixel, C \ SW % DT .9 Because the i, j i, j average value of the wavelet is zero and the data are sampled from a constant amplitude function, the mean correlation value will tend asymptotically to zero, with statistical sampling causing positively and negatively valued deviations from zero in individual pixels. Indeed, for our
9 This notation deviates somewhat from that of Mallat (1998), in which C would be written WD[i, j] ; however, we feel our notation makes i, j complicated expressions in the remainder of this work more easily interpretable.

A

B

(4)

A

BA

BA

B

x[a y[b 1 , . (2) W @(a, b ; p , p ; x, y) 4 W@ xy mp p pp x y xy p and p are the dilation parameters, and a and b are the x y translation parameters. There are many functions / that can be used to create a mother wavelet. One example is the two-dimensional Gaussian function : /

A

B

A

xy 1 x2 y2 , \ exp [ [ , pp 2np p 2p2 2p2 y xy xy x

B

A

B

(3)

8 For a general introduction to wavelets, see, e.g., Mallat (1998) and references therein, and Daubechies (1992).


188

FREEMAN ET AL.

Vol. 138

FIG. 1.õTwo-dimensional Marr, or Mexican Hat, wavelet function (eq. [5])

simple situation, the resulting probability sampling distribution (PSD) of correlation values, denoted p(C o B), will tend asymptotically to a zero-mean Gaussian, with width p P (p p B)1@2 (Freeman et al. 1996 ; Damiani et al. G, C xy 1996, 1997a). We now assume that there is a tightly bunched clump of counts in our otherwise source-free subset image, which could be caused by a Poisson sampling ÿuctuation or by an astronomical source. We assume the clump has Gaussian shape with amplitude A and width p , which is comparaG G ble to the size of the PSF of the instrument. The correlation of the MH function with a Gaussian yields a MH function that has its centroid at the Gaussian centroid and has centroid amplitude C (p , p ) proportional to A : max x y G A pp Gxy C\ max J(p2 ] p2)(p2 ] p2) G xG y p2 p2 G[ G . (6) ] 2[ p2 ] p2 p2 ] p2 G x G y For our simple example, C [ 0, i.e., C is always max greater than the mean of the max PSD. To determine whether we should associate the clump with an astronomical source, we would compute the integral of the PSD from C to O ; max if this quantity is smaller than a predeïned signiïcance threshold (e.g., \10~6, or C [ 4.75p ), then we associG, C ate the clump with a source. max Because p P p p , source detectability will vary as a G, C function of wavelet x y size. In the limits (p , p ) ] 0 and scale y ]O, the ratio C (p , p )(p p )~1@2 goes tox zero, i.e., the max x y x Ify we apply symmetric wavesource becomes undetectable. lets in our analysis, the maximum value of this ratio occurs at p \ p \ J3p . Hence, sources are most easily detected G x y when one analyzes the image with a wavelet function that

has a size similar to that of the source. Since such a function acts to "" ïlter,îî or selectively enhance, structures of similar scale, this behavior is to be expected. If the clump is associated with a previously known source, but is not detected, we can use equation (6) to determine the upper limit on source counts, by substituting the source detection threshold for C , and solving for A (see, max G e.g., ° 4.1). (Because the Gaussian is normalized, A is idenG the use tically the number of source counts.) We stress that of this method to place upper limits on source counts is limited to cases where the PSF shape is that of a twodimensional Gaussian function, and should not be used in the place of simulations if the PSF has arbitrary shape.
3.

ALGORITHM

A

B

3.1. Source Detection Our ïrst objective is source detection : the identiïcation of putative source pixels in binned, two-dimensional image data. This identiïcation is normally done by carrying out the steps in the algorithm described below separately with each of a number of wavelet functions (see also Fig. 2). The basic steps in this algorithm include (1) the correlation of the data with the given wavelet function (° 3.1.1) and (2) the identiïcation of image pixels with correlation values larger than predeïned thresholds for source detection (° 3.1.2). The second step requires knowledge of the local background in each pixel (due to unresolved point sources, diuse astrophysical emission, particle background, etc.). If the background has not been determined previously, then it must be estimated, e.g., via the method described in ° 3.1.3. Also, instrumental artifacts such as support rib structures, hot spots, and the edge of the FOV can adversely aect the second step, resulting in the detection of instrumental features in the image, which are astrophysically uninteresting.


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

189

FIG. 2.õFlow chart illustrating the source detection algorithm described in ° 3.1, as carried out at selected scales (p , p ). Parentheses indicate optional xy input or computations. The local background map may be estimated or provided by the user ; a ÿow chart illustrating how it can be iteratively estimated using the input data (and exposure map) is given in Fig. 3.

In ° 3.1.4, we describe a method for rejecting such "" instrumental sources.îî Finally, in ° 3.1.5 we describe how the variances of the correlation and local background amplitudes are estimated.
3.1.1. Correlation of the W avelet Function and the Data

The ïrst step in the source detection process at a given scale is to compute the correlation of the wavelet function W (p , p ; x, y) with the binned, two-dimensional image y datax D.10 The translation parameters a and b shown in equation (2) correspond to image pixel indices i and j ; for a given pixel, the correlation value is (7) C \; ; W { {D { { i~i , j~j i , j i, j i@ j@ (8) 4 SW % DT . i, j [Henceforth, if a quantity is given with subscripts, we are referring to its value at pixel (i, j) ; otherwise, we are referring to the quantityîs array of values.] The interested reader may ïnd details about how SW % DT is computed in Appendix A.
3.1.2. Computation of the Source Detection T hresholds

S is dubbed the signiïcance (or the Type I error ; see, e.g., i, j Eadie et al. 1971, pp. 215õ216), and it is the estimated probability that we would erroneously identify pixel (i, j) with a source. In ° 2.2, we indicated how S could be determined i analytically in the high-counts limit,, jwhere p(C o B ) tends i asymptotically to a zero-mean Gaussian of, j width (2np p B )1@2. In general, however, p(C o B ) must be x y i, j determinedj via simulations (see Appendix B). i,Because the number of simulations we carried out is only sufficient to directly determine signiïcances S Z 10~7, and because i, j strong sources may have much greater signiïcances (in a qualitative sense), we instead compute a source detection threshold C (S , B ) via the equation o, i, j o i, j = dCp(C o B ) , (10) S\ i, j o Co, i, j where S is the user-speciïed threshold signiïcance.11 If C [ C o , we associate the pixel (i, j) with a source. i, j o, i, j

P

3.1.3. Background Estimation

To determine whether pixel (i, j) should be associated with a source, we compute the probability of observing the correlation value C if there are only background counts i, j present within the support of W : S\ i, j

P

=
j

Ci,

dCp(C o B ) . i, j

(9)

In this section, we describe how the local background counts amplitude B is estimated if it is unknown a priori i, j (see also Fig. 3). While there is no unique way to make this estimate, we seek a method that does not depend on a detectorîs PSF, both for increased generality and computational speed, and also because we want our algorithm to be able to detect and analyze sources of arbitrary size, not just point sources. We can fulïll this condition by creating back11 One choice is S \ P~1, where P is the number of analyzed pixels in o the image ; with this choice, the average number of false detections is one per image.

10 In this section, we do not speciïcally refer to the Mexican Hat wavelet to underscore the fact that our algorithm may be adapted for use with other wavelet functions.


190

FREEMAN ET AL.

Vol. 138

FIG. 3.õFlow chart illustrating the iterative local background estimation method described in ° 3.1.3. Note that a background map is output for each selected scale pair (p , p ). Parentheses indicate optional computations. xy

ground maps at each wavelet scale,12 using a localized function that is wavelet-scale, and not PSF-scale, dependent : the wavelet negative annulus NW (p , p ) (i.e., W with positive values reset to zero). While the x y function NW is related to a wavelet function, we stress that it itself is not a wavelet function and that its use in background estimation does not constitute a transform. If the exposure within the support of NW is constant, and if p and p are sufficiently large such that the integral of the x y source counts distribution (i.e., the PSF for a point source) over the NW is insigniïcant with respect to the integrated background, then one can estimate the background using the formula exp (1) B\ o SNW i, j 4np p x y into a maps are later combined
%

where 4np p /exp (1) is the integrated volume of NW (p , p ) x y ° A3 for derivation). If on the other hand (see xy the exposure is not constant, the exposure map13 can be used as a weighting function : B \E B i, j i, j norm, i, j SNW % DT i, j (12) \E i, j SNW % ET i, j ; @ ; @ NW { { D { { ij i~i , j~j i , j , \E (13) i, j ; @ ; @ NW { { E { { ij i~i , j~j i , j where B is the normalized (i.e., ÿat-ïelded) number of norm, i, j expected background counts at pixel (i, j). We ignore the distinction between vignetted and nonvignetted components of the background (e.g., the particle background) in
13 If one does not provide an exposure map, a ÿat one is assumed, to account for the edge of the FOV.

DT o , i, j

(11)

12 These lation of source properties. See ° 3.2.1.

single map used in the calcu-


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

191

our estimate, because of degeneracy. We note that if the amplitude of the nonvignetted background component BNV is known, then one could in principle estimate the background using a variation of equation (13) : B \E i, j i, j SNW % (D [ B NV)T i, j ] BNV . SNW % ET i, j

There are three situations in which counts from sources may bias the local background estimate : (1) if the estimate is being made within diuse extended emission ; (2) if the pixel in which the estimate is being made is a source pixel, but p x and/or p is smaller than the source size s ; or (3) if sources y are located within the NW . The ïrst situation is a nonissue, because if the analysis goal is detection, say of a source within a supernova

FIG. 4.õIllustration of how source counts may bias a ïnal background estimate (° 3.1.3), causing "" bumps îî if the wavelet scale sizes (p , p characteristic PSF size at a given pixel. We use a 50 ] 50 subïeld of the ROSAT PSPC Pleiades Cluster image in which r B 3 pixels. x p To y PSF as top, but for surface plot of the background estimate for p \ p \ 1 pixel. Middle : same as top, but for p \ p \ 2J2 pixels. Bottom : same x y x y pixels.

) [ r , the PSF : image and p \p \8 x y


192

FREEMAN ET AL.

Vol. 138

remnant, then the diuse emission should be treated as a local background component : for instance, should a clump of counts be associated with a source, or with Poisson ÿuctuations in the background and the diuse emission ? In the second situation, the background map will exhibit a "" bump îî at the location of the source, whose amplitude is greatest where the number of source counts within the support of the NW is maximized, generally at the source centroid (Fig. 4). Within the bump, the source detection threshold is overestimated, and thus this perhaps otherwise detectable source may remain undetected at the given scale. However, this is not a critical problem, since if the source is detectable, it will be detected when (p , p ) Z s, a scale xy regime where the bump is minimized, and the issue becomes moot. Note that these bumps do not adversely aect source property estimation, since they are eliminated when the background map that is used for source property estimates is computed (° 3.2.1). The third situation can be more problematic for source detection. If there are sources in the FOV, then there will always be some image pixels for which the background amplitude is overestimated : for instance, for a symmetric wavelet function, these pixels will surround sources in circular "" rings îî of radius B2 p (the radius at which NW achieves its minimum value ; see Fig. 5). In these rings, C is o overestimated, so that otherwise detectable sources whose locations coincide with these rings may go undetected. Rings can appear regardless of the scale or source size, and thus can impede source detection at all scales. They can also adversely aect the computation of source properties, as a background that is overestimated in the vicinity of a source can lead to underestimated count rates, etc., in the ïnal source list. (This last point is demonstrated below in Figure 16, which shows the eect of rings on the estimation of Pleiades Cluster source properties. See ° 4.2.) One way to remove the rings is to remove source counts iteratively from the raw data, via the following algorithm : 1. Identify pixels to be "" cleansed îî using p(C o B) and the initial background map, which we will dub B ; 1 2. Mask out these pixels or replace their data with other values, creating a new image we denote D ; 3. Estimate B (D ), where n is the 2iteration number nn (n º 2) ; and 4. If the background map is to be reïned yet again, determine C (D ), identify pixels to be cleansed using p(C o B ), n n and return nto step 2. Then, when the ïnal background map B is determined, final compare the original values C with p(C o B ) to create a final ïnal list of putative source pixels. Regarding steps 1 and 4, since the goal of this iterative approach is to remove as many source counts as possible from the raw data, we advocate an aggressive approach to identifying the pixels to be cleansed : the threshold signiïcance should be set high, e.g., to S \ 10~2 (although never o higher than 0.05, which corresponds to the oft-used 95% rejection level of statistics). Regarding step 2, while masking is used by Damiani et al. (1997a) in their two-iteration approach to source detection, it would preclude us from using FFTs to calculate C(D ). Thus, we replace the data in n cleansed pixels with the inferred background amplitude. There are no rigorous quantitative rules governing how one should specify the number of iterations, as that can depend on the crowdedness of the ïeld, the source distribu-

tion, the source strengths, and the wavelet scale size.14 We do note that iterative cleansing will cease if the background map does not change from one iteration to the next, i.e., if no new pixels are marked for cleansing.
3.1.4. T reating the Eect of Exposure V ariations on Source Detection

In ° 3.1.2, we describe how we use probability sampling distributions p(C o B ) to identify sources. These distribui, j tions are derived assuming a spatially constant exposure. If, for instance, the exposure map exhibits localized high-order variations, then the list of detected sources may contain a mixture of astrophysically interesting sources and "" instrumental sources îî aligned near support rib shadows, near the edge of the FOV, at the location of hot pixels, etc. (See, e.g., Fig. 1 of Damiani et al. 1997b.) Thus, an efficacious source detection algorithm should include additional calculations that act to decrease the detectability of instrumental sources, while leaving the detection efficiency of astrophysical sources unchanged.15,16 One possibility is to construct new sampling distributions p(C o B , E) for each observation, taking into account all i, j the exposure variations that can appear within the wavelet support ; however, this is not computationally practical. Instead, we estimate the systematic eect that exposure variations have on the correlation coefficients. Assuming the null hypothesis, we may write the observed correlation coefficient as (14) C \ SW % DT \ SW % BT ] *C , i, j i, j i, j i, j where B is the estimated (noise-free) background intensity and *C i, j is the noise (and possibly source count) contribuj tion to i,C . We ignore the latter term (see the caveats i, below) andj rewrite the former so that its dependence on exposure variations is explicit : SW
%

BT \ SW % EB T i, j norm i, j \; ; W { {E { {B i~i , j~j i , j norm, i{, j{ i@ j@ \ ; ; W { {[E B i~i , j~j i, j norm, i{, j{ i@ j@ [ (E [ E { {)B ] i, j i , j norm, i{, j{ T [ SW % (dEB )T . \ E SW % B norm i, j norm i, j i, j (15)

The quantity dE { { encapsulates exposure variability within the waveleti, j i , j support, and thus the last term in equation (15) encapsulates the eect of exposure variability upon
14 For a typical, uncrowded Chandra ïeld, two iterations (i.e., one round of source count cleansing) are usually sufficient, because the high resolution of Chandra reduces source crowding relative to that observed in, e.g., ROSAT data. However, one should always verify that this is the case with oneîs speciïc image ! 15 The distributions p(C o B ) are also derived assuming a spatially coni, j stant background map, from which simulated data are sampled. Thus, an efficacious source detection algorithm should also include calculations that mitigate the eect of background variations caused by, e.g., X-ray shadows. The current algorithm does not take such variations into account ; anecdotal evidence (e.g., in ° 4.2) indicates that they have little eect, possibly because they are far less "" sharp îî than variations induced by support rib shadows, etc. Note that if the background is known a priori, one can in principle remove the eect of high-order background variations on correlation values using a transformation similar to the one described below. 16 Note that exposure corrections are not mandatoryõfor instance, the user may choose to have no corrections made if the analysis goal is scaleby-scale characterization of sources in correlation space. See the caveats below.


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

193

FIG. 5.õIllustration of how source counts may bias an initial background estimate (° 3.1.3), causing "" rings îî if sources are located within the negative annulus of the wavelet. We use a 50 ] 50 subïeld of the ROSAT PSPC Pleiades Cluster image. Top : image and surface plot of the background estimate after one iteration (i.e., no "" cleansing îî). Middle : same as top, but after two iterations. Bottom : same at top, but after three iterations.

C ; subtracting this term from C yields an "" exposurei, j i, corrected îî quantity that containsj only information of astrophysical value : \ C [ SW C cor, i, j i, j \ C [ SW i, j
% %

If B is constant (or linear) within the wavelet support, norm then equation (16) reduces to Capprox \ C [ B SW % ET . (17) cor, i, j i, j norm, i, j i, j Because SW % ET is computed only once, we dub this the "" fast îî exposure correction, as opposed to the "" full îî exposure correction of equation (16). One should not use the "" fast îî correction if nonlinear structures (caused, e.g., by X-ray shadows) are apparent in the background map.

(dEB

)T norm i, j EB T ] E SW norm i, j i, j

%

B T. norm i, j (16)

It is this quantity that is compared with the distribution p(C o B ) to determine whether (i, j) is a source pixel. i, j


194 One should keep the following caveats in mind :

FREEMAN ET AL.

Vol. 138

1. Strictly speaking, C still cannot be directly comcor, i, j pared with the probability sampling distribution p(C o B ) i, j because the noise term *C is itself uncorrected. If we i, j concentrate on the issue of false positives (i.e., assume that there is no source count contribution to *C ), the impori, j tant question is : is the asymptotic width of the distribution from which *C is sampled smaller than the width of the cor, i, j distribution from which *C is sampled ? If so, then the i, j rate of false detections will still be greater than expected. This is a problem if and only if for a given pixel, E is i, j smaller than the average exposure E over the wavelet ave support, i.e., this is only a problem within troughs or beyond the edge of the FOV. To see this, harken back to the simple example of ° 2.2 : what would happen to the width of p(C o B ) if we were to reduce the exposure ? Fewer counts i, j would be detected, so B would decrease, and the width of i, j the noise distribution, which is P (B )1@2, would also i, decrease. Thus, we suggest that one shouldj carefully scrutinize all sources detected in low-exposure regions ([0.2E ). max 2. Note the distinction between the correlation maps C and C , especially if the analysis goal is not just source cor

detection, but also image decomposition (the scale-by-scale characterization of sources in correlation space).17 The quantity C represents a pixel-by-pixel wavelet ïltering cor not of the raw data D, but of the quantity D [ EB norm ]E B (eq. [16]). Because B may be estimated i, j norm norm using the NW , which is most sensitive to low-frequency components of the data (see Fig. 6), these modiïed data may be "" contaminated îî with low-frequency information (although wavelet ïltering [eq. (16)] mitigates the eect of the contamination). 3. We note that while systematic overestimates of B norm (caused for reasons discussed in ° 3.1.3) adversely aect the computation of C , they will not lead to an increased cor number of false detections. This is because the only situation where such a systematic overestimate aects non17 Note that source characterization in WAVDETECT is done using the raw data themselves and not using correlation coefficients. Aside from source detection, the only other place where the exposure-corrected correlation map is used in WAVDETECT is in the creation of the noise-free, exposure-corrected image of detected sources (° 3.2.5). Thus, this caveat is only an issue if the user desires to analyze correlation maps outside of WAVDETECT.

FIG. 6.õSample power spectra of the Mexican Hat wavelet function (W ) and of the negative annulus of the Mexican Hat wavelet function (NW ). N is the number of pixels in the (padded) image. The NW has much larger response to low-frequency components, but this signal would be suppressed upon correlation with W (see ° 3.1.4).


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

195

source pixels in the ïnal, reïned background image is when the background is computed in regions of extended diuse emission. In this situation, the algorithm treats the sum of the real background and the diuse emission as the "" background,îî so the rate of false detections (which is independent of background amplitude) will be unchanged.
3.1.5. V ariance Estimation

We estimate the variances of B (if one does not provide i, j a background map), and of C (or C or Capprox ), using i, j cor, i, j cor, i, j the standard formula (Eadie et al. 1971, p. 23) V [Y ] \ V [; ; a X ] i, j i, j ij \ ; ; a2 V [X ] i, j i, j ij ]2; ; ; ; a a { { cov[X , X { {] (18) i, j i , j i, j i j i i@;i j j@;j where Y is quantity of interest and X are random varii, ables (either D or functions of D ). Forj instance, i, j i, j V [C ] \ V [; ; W { { D { {] i~i , j~j i , j i, j i@ j@ \ ; ; W 2 { { V [D { {] i~i , j~j i ,j i@ j@ (19) \; ; W2 { {D{ { i~i , j~j i , j @ j@ i \ SW 2 % DT . i, j Note that we make two assumptions when deriving this formula : (1) the datum D { { is sampled from a Poisson ,j distribution with variance i D { { ; and (2) each pixelîs raw i ,j datum is independently sampled (so covariance terms do not contribute to V [C ]). i, j We list variance formulae related to source detection in Table 1. The reader should keep in mind two important caveats about them : 1. These formulae ignore the contribution of the covariance terms, which are nonzero for V [C ], V [Capprox], and cor V [B] if the data are iteratively cleansed,cor if B is a not just i.e., a function of the raw data D only. We ignore these terms because even the simplest covariance computation, that for a two-iteration background map (see Appendix C), has a staggeringly high computational cost : we ïnd that the CPU time needed to compute the variance increases by a factor DO(d d p2 p2), where d and d are the x- and y-axis x y pixels, respectively. Also, additional arrays conx y lengths in x y taining information needed to compute the covariance

terms must be kept in memory, so there is a resource cost as well. We ïnd that including covariance terms increases the variance by a median value of B7%, and at most by only B30% adjacent to strong sources, although this is a sourcestrengthõ and source-geometryõdependent result that obviously cannot be blindly applied to all ïelds. Ultimately, it is up to the user to judge whether adding the computation of covariance to our base algorithm is worthwhile. 2. When computing V [B ] for the ïnal background i map, we make the simplifying, jassumption that the variance of a cleansed datum is equal to the cleansed datum itself ; for instance, for a two-iteration background map, V [B ] \ ; i, j i@ \; i@ where NW { { i~i , j~j . \E (21) i, j SNW % ET i, j We make this assumption because if the data are a mixture of raw data and background estimates, the variance estimates become increasingly complicated : for one iteration, a i, i{, j, j
{

; j@ ; j@

a2 { { V [D { {] i, i , j, j 2, i , j a2 { { D { { , i, i , j, j 2, i , j (20)

V [D ] \ D , i, j i, j for two iterations,

4D j V [D ] \ 5 ;i,; a2 D 2, i, j 6 i@ j@ i, i{, j, j{ i{, j
etc.

0 0

uncleansed pixel ,
{

cleansed pixel ,

(22)

3.2. Source Characterization Once we have identiïed putative source pixels at each of a number of wavelet scale size pairs (p , p ), our next objective is source characterization, wherein x combine informawe y tion derived at each scale pair to generate a ïnal source list and to estimate source properties (see Fig. 7). Unlike source detection, source characterization algorithms can be arbitrarily complex, depending, for instance, upon whether one wishes to use detailed PSF information. Our method is particularly simple, in that we use only the characteristic PSF size at a given pixel, r . This size may be associated with, e.g., 50% encircledPSF, i, j ; we ïnd that smaller values, energy such as 39.3% (which corresponds to the integral of a symmetric normalized two-dimensional Gaussian to radius p ), G work better than large values. While using the detailed PSF

TABLE 1 SOURCE DETECTION VARIANCE FORMULAE Property Correlation (C ) ............................................... i, j Exposure-corrected correlation (C cor, i, j ) ..................... SW 2
%

Variance DT i, j
%

V[C ] ] i, j ])T ] 2E SW 2 SW 2 % (E2V [B norm i, j i, j V[C ] ] B SW 2 i, j norm, i, j
%

(EV [B ])T ] E2 SW 2 norm i, j i, j

%

V [B

norm

])T

i, j

Approximate exposure-corrected correlation (C approx ) ...... cor, i, j Normalized background (B

E2T i, j

) ........................... SNW 2 % D T /(SNW % ET )2 norm, i, j N i, j i, j NOTES.õi and j are pixel indices, and N is the number of iterations used to compute the background map, if it is not provided.


196

FREEMAN ET AL.

Vol. 138

FIG. 7.õFlow chart illustrating the source list generation algorithm described in ° 3.2. Parentheses indicate optional input.

shape may allow for more accurate estimates of source properties, the simplicity of our scheme makes it more immediately applicable to images from virtually any counts detector (including those for which calibration is on-going or is for other reasons incomplete).
3.2.1. Corrected Background Estimate

values when extended sources are analyzed. We use equation (18) to estimate the variance of B@ , with covarinorm, i, j ance terms ignored : N v pp 2 i, j, k x, k y, k ; V [B ]. (25) norm, i, j, k ;N v pp k/1 k/1 i, j, k x, k y, k V [B ] is estimated using the approximate equation norm, i, j, k listed in Table 1.
3.2.2. Source Cells

A

B

In ° 3.1.3, we describe how we calculate a scale-dependent normalized local background estimate B , by assuming that there are no source counts in the norm negative annulus, NW . Source detection itself is not markedly aected if this assumption is violated and the background overestimated, for reasons given in ° 3.1.3, but an overestimated background will adversely aect the computation of source properties. We create a new, corrected, background estimate by combining information across scales (denoted with a subscript k), noting that the assumption that there are no source counts in NW is always violated around sources if either p or p is less than the source size : x, k y, k ;N v ppB B@ \ k/1 i, j, k x, k y, k norm, i, j, k . (23) norm, i, j ;N v pp k/1 i, j, k x, k y, k N is the number of scale pairs used, and

4 1 min (J2px , J2py) º mrPSF, i, j , 0 v \5 i, j, k 6 0 otherwise . 0

(24)

The quantity m is a multiplicative factor, set to 1 when estimating the properties of point sources, and to larger

The next step in source characterization is to determine which pixels of the original image D are to be associated with each detected source. We term contiguous pixels which are associated with a particular source a source cell, and as we describe below in ° 3.2.4, we use the raw data in a source cell to determine a sourceîs properties. We need to create source cells because, as noted above, we do not use PSF shape information in our algorithm, and because we want an algorithm that is applicable to both point and extended sources. Source cells are not computed in algorithms such as the sliding cell, where the integrated PSF volume and sum of (point) source counts for a given user-deïned cell can be used to estimate the total number of (point) source counts, etc. To create source cells, we ïrst compute source count images, smoothing the raw data with the positive kernel of the wavelet function, PW , at user-speciïed scales, and then


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

197

subtracting the corrected background image, B@ : SC i, j, k \ max

DT i, j, k [ B@ , 0 . (26) norm, i, j ET i, j, k (We denote these images SC to avoid confusion with either the signiïcance S or the correlation image C.) We use PW as the smoothing function because it has the desirable properties of being localized, and, for the particular case of a symmetric MH function, of mimicking the shape of a canonical Gaussian PSF. In regions where there are no sources, source count image values are either zero, or positive and nearly zero ; only in the vicinity of sources do the values deviate markedly from zero. Thus, the source count images appear to contain numerous "" islands îî of nonzero ÿux in a sea of zero values, with their relative size increasing with size of the smoothing function PW (see Fig. 8). Each island contains one or more peaks, and subislands may be deïned using each peak, with saddle points providing the boundaries between them. (Sub-) islands observed in selected source count images deïne the source cells. To deïne a source cell, we must select a source counts image and then determine to which (sub-) island the putative source belongs. The selection proceeds as follows. The location of a putative source in correlation-space is assumed to be the location of the correlation maximum,
% %

A

SPW SPW

B

(i , j ). At this location, we compute the PSF size, in pixels, CC. r (This introduces a bias toward point sources ; we C PSF, iC , jto this point below.) We then select the source count return image with smoothing scale "" closest îî to r by miniC , jC mizing o log p [ log r o , where pPSF, i deïned for is 2k 2 PSF, iC , jC k each scale pair : p \ exp k

C

log (p ) ] log (p ) x y. 2

D

On the selected source count image, we examine pixel (i , j ): the (sub-) island to which this pixel belongs deïnes CC the source cell. A source cell deïned in this manner has advantageous , then nearly all isolated point properties : (1) if p B r C , jC k source counts should PSF, iwithin a source cell (Fig. 9) ; (2) lie exposure variations are taken into account via the use of SPW % ET in equation (26), so that source cells are not truncated near, e.g., support rib shadows ; and (3), as noted above, saddle points in the source count images provide natural boundaries between sources in crowded ïelds (Fig. 10). We note two situations where care must be exercised when interpreting results. First, the source cell for an extended source may be too small if the smoothing scale is . The steps one must take to deal with this situBr PSF, iC , jC

FIG. 8.õSample source counts image (° 3.2.2), created by ïrst smoothing the data in a 50 ] 50 subïeld of the ROSAT PSPC Pleiades Cluster image with a PW function of size p \ p \ 2 pixels, then subtracting the estimated background. Pixels associated with the strongest local maxima comprise the source x y cells that are used for source property estimation. See Figs. 9 and 10.


198

FREEMAN ET AL.

Vol. 138

FIG. 9.õTop : counts data showing an isolated Pleiades Cluster source, as observed by the ROSAT PSPC. Middle : source counts image data, created smoothing the counts data with a PW function of size p \ p \ 2 y pixels, then subtracting the estimated background. Bottom : the x source cell deïned using the source counts image data (° 3.2.2). This cell, used in the estimation of source properties, contains nearly all, if not all, of the counts from this source.

ation will vary depending upon analysis circumstances, but one possible step is to create only one source count image, , and to use this image to deïne the cell, with p ? r PSF, iC , jC while k being careful to note whether previously detected point sources are located within it. (See, e.g., ° 4.1.) Another situation for which care must be exercised is when the PSF

FIG. 10.õTop : counts data showing two nearly overlapping Pleiades Cluster sources, as observed by the ROSAT PSPC. Middle : source counts image data, created smoothing the counts data with a PW function of size p \ p \ 2 pixels, then subtracting the estimated background. Bottom : x y the source cells deïned using the source counts image data. The saddle point seen in the middle image deïnes the boundary between the cells (° 3.2.2).

is bimodal or otherwise strangely behaved (such as the oaxis Chandra PSF) ; two or more source cells could be created for one detected source. The necessary steps to deal with this situation depend upon the details of the detector


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA
TABLE 2 SOURCE PROPERTY EXPECTATION VALUES Property X location (pixels) ........................ Y location (pixels) ........................ Bkg counts ................................ Counts ..................................... Exposure ................................... Count rate (counts per exp. unit) ....... Bkg count rate (counts per exp. unit)... NOTES.õi and j are pixel indices ; sc within the sourceîs cell ; D \ ¸ ¸ s i map. in the same units as the exposure|sc j Expectation Value

199

itself, and thus we will not discuss this particular situation further here.
3.2.3. Source Rejection

Because the same source will generally be detected at multiple scales, and to further decrease the possibility of ïnding false sources, it is necessary to reject sources from the lists of correlation maxima generated at each scale. A maximum observed at (i , j ), for the scale pair (p , p ), is CC x , k y, k rejected from further consideration if any of the following conditions are met : (1) (i , j ) lies in a previously deïned C \ C0 ; (3) the ellipse deïned by source cell ; (2) SC iC , jC , k (x [ i )2 (y [ j )2 C] C \1 2p2 2p2 x, k y, k contains one or more previously deïned sources detected at smaller scales (this can occur when p or p is Zr , k y, k PSF since previously identiïed sources willx,eventually merge if the ïeld is crowded, creating "" new îî sources at new locations) ; and (4) if, after the source cell is deïned for a particular scale, it is found not to contain any correlation maxima at that scale. This last check is aimed at rejecting small-scale Poisson ÿuctuations that may be observed in the background data.
3.2.4. Source Properties

x \ (1/D )¸ ¸ Di s s i |sc j |sc i, j y \ (1/D )¸ ¸ Dj s s i |sc j |sc i, j B \¸ ¸ E B@ s i |sc j |sc i, j norm, i, j C \D [B s s s t \ (t /D )¸ ¸ DE s o s i |sc j |sc i, j i, j R \ C /t s, C ss \ B /t B ss those image pixels which lie and t is the total exposure, o

R s, denotes D; |sc i, j

exposure variations removed : NC SD \ ; cor, i, j, k l . (27) i, j i, j, k pp k/1 x, k y, k Dividing the correlation value by p and p restores the x, y, k normalization contained in equationk (5), allowing a scaleby-scale summation. The quantity l \ 1 if (1) C [ i, j, k 0, (2) the local maximum corresponding to (i, j) cor, i, jbeen has , k identiïed as a source pixel, and (3) the associated local correlation maximum is contained within a source cell (the second condition ensures that random maxima which are not associated with a source but which happen to lie within a source cell are not included in the source image ; the last condition ensures that rejected sources are not included) ; otherwise, l \ 0. i, j, k
4

We present the formulae we use to estimate source properties and their variances in Tables 2 and 3, respectively. The summations performed when making these estimates are carried out over the pixels within the source cell. We use the raw counts data, D , as a weighting function, j instead of the source ÿuence D i,[ E B@ , because the i, j i, j the j use of the latter can greatly complicatenorm, i, estimation of variances. Using the data rather than the source ÿuence will lead to similar estimates when the background amplitude is small relative to source amplitude.
3.2.5. Noise-Free Source Image

. VERIFICATION

We can use the information present in the correlation images and the source cell image to create a "" noise-free îî rendering of the observed source data, SD, with the eect of

To verify its source detection and characterization capabilities, we apply our algorithm to : (1) 1 and 10 ks observations by an idealized detector with a spatially invariant PSF ; (2) a 32 ks ROSAT PSPC observation of the Pleiades

TABLE 3 SOURCE PROPERTY VARIANCE FORMULAE Property X location (x ) (pixels) .............................. s Y location (y ) (pixels) .............................. s Bkg counts (B ) ...................................... s Counts (C ) .......................................... s Exposure (t ) (s) ..................................... s Count rate (R s, C ) (counts per exp. unit) .......... Variance (1/D2)¸ ¸ D (i [ x )2 s i |sc j |sc i, j s (1/D2 )¸ ¸ D ( j [ y )2 s i |sc j |sc i, j s ¸ ¸ E2 V [B@ ] i |sc j |sc i, j norm, i, j D ] V [B ] s s (1/D2 )¸ ¸ D (t E [ t )2 s i |sc j |sc i, j o i, j s (1/t2 )(V [C ] ] R2 V [t ]) s s s, C s

) (counts per exp. unit) ...... (1/t2 )(V [B ] ] R2 V [t ]) s, B s s s, B s NOTES.õi and j are pixel indices ; sc denotes those image pixels which lie within the sourceîs cell ; D \ ¸ ¸ D ; and t is the total exposure, in the same units as the s j deïned exposure map. V [B@ i |sc ] is|sc i, j in eq. o(25). norm, i, j Bkg count rate (R


200

FREEMAN ET AL. closest "" standard îî scale to the assumed PSF size (if we assume scale sizes separated by factors of J2 rather than 2 for greater source detection efficiency). In both images, our algorithm detects one false source, which is consistent with the assumption of S \ 10~6 for a 512 ] 512 image. o Creating source cells by using information derived at a wavelet scale close to the PSF size is not optimal when one wishes to analyze and describe extended sources, as discussed in ° 3.2.2, and indeed by examining Figure 13 we can see that the extended source cells are undersized. There are many ways by which an analyst may wish to treat extended sources ; here, we show how one could derive an image showing the (normalized) number of counts per pixel within the extended source. We use as our example the largest extended object in Image B. First, we must expand the source cell so that it just encloses the extended source. To do this, we increase the minimum scale size at which a source counts image is to be computed (in this example, from 2J2 pixels to 8J2 pixels ; see Figs. 14aõ14b). (Note that the parameter m in eq. [24] must also be increased so that the background is not overestimated within the source.) Second, we would use the new source cell as a spatial ïlter, applying it to the original source counts image created at the PSF scale (e.g., Fig. 14c), or to the data, etc. (If one outputs source count image data, one can then, in principle, ït directly to them. For instance, the image may be of a galaxy cluster, and one may wish to assess the detectability of its constituent galaxies. However, care must be exercised since the data in contiguous pixels are not independent. Taking into account the width of the PW used to smooth the raw image data, we can state that data more than J2 p pixels apart are independent. Thus, simple statistical ïtting can be done to a sparse grid of data. This process of ïtting would be essentially equivalent to the "" decimation îî method described by Lazzati et al. 1998, except that they ït to correlation image data.) 4.2. ROSAT PSPC : T he Pleiades Cluster Next, we demonstrate that our algorithm outperforms the sliding cell in efficiently detecting sources in a crowded ïeld, by applying it to the deep (32 ks) ROSAT PSPC observation of the core of the Pleiades Cluster (RP200068, cf. Micela et al. 1996). These data were obtained in two segments separated by roughly six months, and the slight boresight oset between the two segments has been corrected using a method described by Micela et al. The data were also ïltered to exclude times of high background contamination, and to exclude pulse-heightõinvariant (PI)

Cluster ; and (3) a simulated 30 ks Chandra ACIS-I observation of the Lockman Hole region. These tests allow us to demonstrate that our algorithm can efficiently detect and accurately describe well-sampled sources in an uncrowded ïeld and can eectively analyze crowded ïelds, even in the low-background limit of the Chandra detectors. 4.1. Idealized Detector with Spatially Invariant PSF We ïrst demonstrate that our algorithm can efficiently detect, and accurately describe, well-sampled sources in an uncrowded ïeld. We apply it to two 512 ] 512 images, hereafter Images A and B, that represent 1 and 10 ks observations by an idealized detector with an eective area 1000 cm2 and a spatially invariant Gaussian PSF of width p \ 2.56 pixels (Fig. 11). The exposure map for this PSF detector is similar to that of the Einstein IPC. Within each image we randomly place 42 point sources, and four extended sources with elliptical shape. The ÿuxes of the point sources were sampled from a log N [ log S distribution with slope [1.5, above 10~14 erg cm~2 s~1. We also simulate a locally variable background (amplitude D10~5 counts s~1 pixel~1) by setting the background amplitude at ïve reference points, performing minimum-curvaturesurface interpolation, and sampling background data in each pixel. We analyze the images assuming the input parameters listed for Test 1 in Table 4. In Figure 12, we plot the source counts for detected sources, and the upper limits for undetected sources, against the number of predicted counts. Upper limits are deïned using the source detection threshold values at the correlation maxima nearest the location of the undetected sources and are computed using equation (6), with p \ p \ J3p \ 4.43 pixels. We conclude that x y our algorithm efficiently G detects and describes point sources with Z10 counts. This is not an absolute quantity : the minimum number of counts needed for source detection varies as a function of the background amplitude and source size (see ° 2.2). Hence, while we may conclude that a 1 ks observation is sufficient for the detection by Chandra of nearly on-axis point sources with ÿuxes Z10~14 erg cm~2 s~1, since it has a similar eective area as, and lower expected background count rates than, our idealized detector, it may not be sufficient far o-axis, or for detectors that have higher rates of background accumulation. In Figure 13, we show the cells for the detected sources, along with the input source locations. We use a wavelet scale of 2J2 pixels to create the source counts image that is in turn used to delineate the source cells. This scale is the

TABLE 4 NUMBER OF DETECTED SOURCES : ROSAT PSPC PLEIADES IMAGE Scale Separation 2 J2 2 J2 2 J2 2 J2 Exposure Correction Fast (eq. [17]) Fast Fast Fast Full (eq. [16]) Full Fast Fast Number of Sources 129 136 130 137 132 138 144 148

Test 1 2 3 4 5 6 7 8 ...... ...... ...... ...... ...... ...... ...... ......

Iterations 2 2 3 3 3 3 2 2

Signiïcance 10 10 10 10 10 10 10 10 ~6 ~6 ~6 ~6 ~6 ~6 ~5 ~5

NOTES.õScale sizes range from p \ 1 to p \ 16 pixels. The threshold signiïcance for source cleansing is S \ 10~2. o


FIG. 11.õTop : 512 ] 512 pixel image (Image A) showing a 1 ks observation by an idealized detector with eective area 1000 cm2, a spatially invariant Gaussian PSF of width p \ 2.56 pixels, and an exposure map similar to that of the Einstein IPC (° 4.1). Within this image were placed 42 point sources and PSF four extended sources. The background is assumed to be locally variable, with amplitude D10~5 counts s~1 pixel~1. Bottom : same as top, except the observation time is 10 ks (Image B).


202

FREEMAN ET AL.

FIG. 12.õTop : comparison of observed counts, with estimated 1 p errors, and upper limits for undetected sources, with predicted counts for the point sources of Image A. Upper limits are deïned using the source detection threshold values at the correlation maxima nearest the location of the undetected sources and are computed using eq. (6), with p \ p \ J3p \ 4.43 pixels. The source exposure time (see Table 2), not the total observation time, is used to x y compute the predicted counts. Bottom : same as top, but forGImage B. The two farthest outliers on this ïgure correspond to observed "" sources îî which are actually composed of two sources each.

channels at both the low-energy (PI \ 20, to avoid the socalled "" ghost image îî problem ; Nousek & Lesser 1993) and high-energy (PI [ 201, where no instrument map is available to determine exposure variations) ends of the spectral response. We have computed exposure maps taking into account these changes using software developed by Snowden & Kuntz (1998). In Table 4, we show how the number of detected sources varies as a function of the number of iterations, the spacing of scales, the exposure correction method, and the source detection signiïcance. Comparing Tests 1 through 6, for which S is constant, we ïnd that the number of detected o sources changes little if the number of iterations is increased beyond two, or if the full exposure correction method (eq. [16]) is used instead of the fast one (eq. [17]). However, there is an B5% increase in source yield by analyzing the image with the scale sizes spaced by factors of J2 instead of 2. The extra sources are relatively weak sources whose probability of detection is maximized around the scales J2 pixels, 2J2 pixels, etc. It is the decision of the user as to whether the increase in weak-source detection efficiency is worth nearly doubling the computation time. We compare the source detection results for our Test 8 with those shown in the WGACAT and ROSATSRC catalogs (White, Giommi, & Angelini 1994, and Voges et al. 1994, respectively),18 as well as those shown in Micela et al. and Damiani et al. (1997b). (The WGACAT and ROSATSRC catalog teams, and Micela et al., use variants of the sliding-cell algorithm.) We choose Test 8 because it most closely resembles the analysis of Damiani et al., who perform a two-iteration analysis of the Pleiades image with wavelet functions of scale sizes 1, J2, . . ., 16 pixels, and with threshold signiïcance S \ 1.33 ] 10~5 (4.2 p). Our result is o nearly equal to that of Damiani et al., as we detect 148 sources, while they detect 150. We cannot directly compare our numerical results with those of Damiani et al., who only
18 Available from HEASARC : http ://heasarc.gsfc.nasa.gov/.

publish ïgures showing the performance of their algorithm on a subset of the Pleiades ïeld. However, we can emulate Damiani et al. by showing how our results compare with those of the WGACAT and ROSATSRC catalogs (Table 5, cf. Damiani et al. Table 1). Our results are virtually identical to those of Damiani et al., in, e.g., how many sources are detected only by our algorithm, etc. We thus may conclude that our algorithm and that of Damiani et al. generate a largely similar source list. A large fraction of the sources that we detect, but that are not included in the WGACAT or ROSATSRC catalogs, are located either near the inner telescope support ring (B20@ o-axis) or near the precipitous drop-os in exposure caused by telescope vignetting near the edge of the FOV (Fig. 15 ; see also Fig. 16). A visual examination of these sources indicates that they are not spurious. We further examine in detail the two sources near the edge of the FOV that were not included in either the WGACAT or ROSATSRC catalogs but are in regions covered by other
TABLE 5 COMPARISON AMONG SOURCE DETECTION ALGORITHMS : ROSAT PSPC PLEIADES IMAGE Method Our Algorithm (Test 8) ................. Damiani et al. 1997b ................... Micela et al. 1996 ....................... WGACAT ............................... MPE ...................................... Common to : All methods ........................... Our algorithm and WGACAT ...... Our algorithm and MPE ............ WGACAT and MPE ................ Our algorithm only .................. WGACAT only ....................... MPE only ............................. Number of Detected Sources 148 150 99 127 100

81 27 12 3 28 16 4


FIG. 13.õTop : source cell image generated during the analysis of Image A, for wavelet scale size 2J2 pixels. The overlaid small circles and larger ellipses represent, respectively, the 42 point sources and four extended sources randomly placed in the image. Bottom : same as top, but for Image B.


204

FREEMAN ET AL.

Vol. 138

ROSAT PSPC pointings described by Stauer et al. (1994). We ïnd that these sources lie within r pixels of the StaufPSF fer et al. sources 159 and 197 (see their Table 2), which have reported count rates B0.023 and 0.016 counts s~1, respectively. These rates are consistent with our derived count rates. This demonstrates the robustness of our simple exposure correction method. We also note the intriguing result that the local background map computed by our algorithm indicates the presence of an X-ray shadow in the core of the cluster (Fig. 17 ; Kashyap et al. 2002, in preparation). Because the Pleiades cluster is located beyond the edge of the local bubble of hot gas (Frisch 1995), this shadow, which is most pronounced at low energies, is of more distant sources of the diuse X-ray background (DXBG), such as the extragalactic component and obscured stars in the Pleiades cluster itself. The depth of the shadow places constraints on the nature of the stellar mass-function at low masses, and in particular rules out models where the mass-function is extrapolated at a constant slope from higher masses, thus providing independent, X-ray observational support for optical observations that report drops in the mass function at low masses (M D 10 ; B see, e.g., Tinney, Mould, & Reid 1992 ; Bahcall et al. 1994). 4.3. Chandra ACIS-I : T he L ockman Hole Finally, we demonstrate that our algorithm can handle the low background amplitudes which are characteristic of Chandra observations, while continuing to outperform the sliding cell, by applying it to a simulated 30 ks Chandra ACIS-I image of the Lockman Hole (T. Gaetz 1998, private communication ; Fig. 18). The ACIS-I is comprised of four 1024 ] 1024 pixel CCDs conïgured in a 2 ] 2 square, with FOV B 17@ ] 17@ (B50 times smaller than that of the ROSAT PSPC). Within the ACIS-I ïeld are placed : (1) 12 optically identiïed ROSAT PSPC X-ray sources catalogued by Schmidt et al. (1998), including one extended cluster source ;19 (2) B6000 point sources sampled from the Hasinger et al. (1998) log N [ log S distribution between 10~17 and 5 ] 10~15 erg cm~2 s~1 ; and (3) 19500 particle background counts (corresponding to a particle background rate of 1.5 ] 10~7 counts pixel~1 s~1). If we assume the input parameters listed for Test 1 in Table 4, and do not use an exposure map, we detect 171 sources in the full ACIS-I ïeld, of which four are directly associated with the extended galactic cluster (Fig. 18). In Figure 19, we show dierential log N [ log S distributions for both detected, and all, points sources in the FOV. On the basis of this ïgure, we may conclude that our algorithm will efficiently detect sources in ACIS-I images with ÿuxes Z10~15 erg cm~2 s~1 (0.5õ2 keV) and has the ability to detect Poisson sampling ÿuctuations for sources with ÿuxes [10~16 erg cm~2 s~1. Our result compares very favorably with the performance of "" CELLDETECT,îî which detects 51 sources with S/N º 3. (Further comparison of the source detection efficiencies of WAVDETECT and CELLDETECT is provided by Kim et al. 2002, in preparation.) We conclude that the relative superiority of our wavelet detection algorithm, with respect to the sliding cell, is inversely proportional to the background amplitude. In Figure 20, we plot the osets of the locations of detected sources from their actual locations within the
19 Other detected X-ray sources that lie within the ACIS-I FOV but were not optically identiïed, have not been included.

FIG. 14.õTop : subset of the source-cell image shown at the bottom of Fig. 13, showing the cells generated in the vicinity of the largest extended source in Image B. Middle : same as top, except that now the source cells are deïned using a source counts image with smoothing size 8J2 pixels instead of 2J2 pixels. Bottom : the normalized count-rate image for the extended source (and three point sources), generated by creating a source counts image with smoothing size 2 pixels, then ïltering out all pixels not contained in the source cell shown immediately above. If this source were a galaxy cluster, the analyst could proceed to ït to these data (respecting the caveats discussed in ° 4.1).


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

205

FIG. 15.õTop left : full ROSAT PSPC image of the Pleiades Cluster. T op right : ellipses representing the 148 sources detected by our algorithm. The ellipse sizes are set by deriving the 1 p principal axes and rotation angle for each source. Middle left and right : same as top left and right, with only the central 1¡ ] 1¡ portion of the image shown. Bottom left and right : same as top left and right, with only the central 30@ ] 30@ portion of the image shown.

FOV. (These osets are caused by the asymmetry of the Chandra PSF, and their magnitude increases with o-axis angle.) We note two characteristics of these osets. First, the variation in the osets does not depend on source strength, signifying that the source-locating process is insensitive to

the strength of the source. Second, the observed osets are signiïcantly smaller than the expected mean separation between sources (B50 pixels), implying that the observed osets are due to the asymmetries inherent in the PSF and not due to source misidentiïcations. We thus ïnd that the


FIG. 16.õIllustration of how the number of iterations used to estimate the local background amplitude (° 3.1.3) aects the estimation of Pleiades Cluster source properties (° 4.2). Left : the ratio of source predicted counts C /C given backgrounds estimated using one iteration, and two iterations, as a function of 12 C . Right : same as left, but instead computed for two and three iterations. 2

FIG. 17.õCorrected background map generated during the analysis of the ROSAT PSPC image of the Pleiades Cluster (° 4.2) via the method of ° 3.2.1. The contour levels are 0, 0.9, 1.05, 1.2, 1.35, and 1.5 counts, with darker areas having more counts. Since the estimated error in each pixel is D0.01 counts, the perceived structure is real and is indicative of X-ray shadowing (Kashyap et al. 2002, in preparation).

206


FIG. 18.õTop : simulated 30 ks Chandra ACIS-I observation of the Lockman Hole. Data in all four chips are shown, rebinned by a factor of 2 for greater visual clarity. The gaps between chips are B15 pixels. Bottom : ellipses representing the 171 sources detected by our algorithm. The ellipses, whose sizes are normally set by deriving the 1 p principal axes and rotation angle for each source, have minimum axis lengths of 10 pixels for greater visual clarity.

207


208

FREEMAN ET AL.

Vol. 138

FIG. 19.õDierential log N [ log S distributions for detected point sources (dotted line) and all point sources (solid line)inthe Chandra ACIS-I Lockman Hole ïeld. We conclude that our algorithm will efficiently detect sources in ACIS-I images with ÿuxes Z10~15 erg cm~2 s~1 (0.5õ2 keV) and has the ability to detect Poisson ÿuctuations for sources with ÿuxes [ 10~16 erg cm~2 s~1.

weak sources are as "" well-behaved îî as strong sources (about whose detection and identiïcation there can be little doubt), and hence infer that even the weakest detected sources are real.
5

Grebenev et al., and Lazzati et al., with the exception that Rosati et al. use the symmetric Morlet wavelet function 2 1 W (r) \ e~r2@p2 [ e~r2@2p2 M p2 2

C

D

(28)

. COMPARISON WITH EXISTING ALGORITHMS

The source detection algorithm which we present in this work resembles algorithms published previously by Vikhlinin et al. (1994), Rosati et al. (1995), Grebenev et al. (1995), Damiani et al. (1997a), and Lazzati et al. (1998). In this section, we highlight important dierences between our algorithm and these other algorithms, all of which, having been developed for analyzing ROSAT PSPC data, suer from inherent limitations that do not allow them to be directly applied to data from, e.g., Chandra. We do not discuss how our method of source characterization diers from those described previously, because these other methods are built upon the premise that the PSF has Gaussian shape. Thus, they are simply not directly applicable in situations where the PSF has a more complex shape. 5.1. Correlation Image Our basic method for computing the correlation images is the same as that used by Vikhlinin et al., Rosati et al.,

instead of the MH function. However, these authors do not attempt to correct for exposure variations, as they focus their attention upon the center of the ROSAT PSPC FOV, and they use a fundamentally dierent method to determine source detection thresholds (see below). The method by which Damiani et al. correct for exposure and detect sources diers substantially from ours. They ïrst divide the raw data image (which they refer to as the "" photon image îî) by an exposure map on a pixel-by-pixel basis, to create the so-called "" count-rate image.îî Pixels with relative exposure less than a certain amount (e.g., 0.2) are not included, which introduces a sharp edge in the countrate image where telescope vignetting becomes important. Damiani et al. correlate the wavelet function with this image, applying an analytic correction to correlation values near the edge (as given in their eq. [12]). Because the data in the count-rate image are not Poisson-distributed, Damiani et al. must convert source detection threshold values, derived from their background map in the photon image, to


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

209

FIG. 20.õOsets of the locations of detected sources from their actual locations within the Chandra ACIS-I Lockman Hole ïeld, as a function of o-axis angle. These osets are caused by the asymmetry of the Chandra PSF, and the solid line represents its 95% encircled energy radius. Triangles represent sources with fewer than 6 counts ; crosses, sources with between 6 and 15 counts ; and stars, sources with more than 15 counts. This ïgure demonstrates that our ability to associate detected sources with actual sources does not vary as a function of source strength.

values appropriate for the count-rate image. They accomplish this by dividing the photon-image detection thresholds by an eective exposure time t , which is not a source eff exposure time that can be used to convert estimated source counts to count rates. Because Damiani et al. derive the equation with which they compute eective exposure time in the Gaussian limit (see their Appendix B), their exposure correction method is eectively limited to the high-counts regime : it cannot be applied as is to, e.g., typical Chandra data. 5.2. Background and Source Detection As described in ° 3, we estimate the local background counts amplitude in each pixel by assuming that within the negative annulus of the wavelet, NW , there are no source counts (eq. [12]). We then use that inferred amplitude to determine source detection thresholds. This is a generalization of the approach used by Vikhlinin et al., Rosati et al., Grebenev et al., and Lazzati et al., in which the correlation variance SW 2 % DT is used to determine thresholds (e.g., SW % DTSW 2 % DT1@2 º 3.5). This approach will work

if the background is locally ÿat or has a locally constant gradient, in which case the correlation of wavelet and background is zero. Such an approach is obviously insufficient for use with the whole FOV of X-ray detector, where support rib shadows and vignetting will aect the background, and/or in situations where the noninstrumental background varies markedly (such as in the Pleiades ; see Fig. 17). Damiani et al. use the same basic approach to source detection as we do, in that they compute a local background amplitude, and use it to compute source detection thresholds. Damiani et al. compute the background map by ïrst smoothing the raw data with a Gaussian with width º2p . They then interpolate background estimates PSF, i, j within support rib shadows. At each scale, they compute the median value of the smoothed background within a square region with side-length l \ 4(p2] p2 ) centered on the PSF, i, j given pixel. (We explicitly use p2 instead of r2 to highPSF light their assumption of Gaussian PSF shape.)PSF After point sources are detected (using a source cleansing signiïcance that is 0.2 p smaller than their source detection signiïcance),


210

FREEMAN ET AL.

Vol. 138

regions around them containing 95% of the counts are masked out, and a reïned background estimate is made by interpolating over the masked regions.
6

. SUMMARY

In this work, we present a generalized wavelet-based source detection algorithm that, in principle, can be applied immediately to image data collected by any photon counts detector, although it was developed speciïcally for the analysis of Chandra X-ray Observatory image data. We exclusively use the Marr, or Mexican Hat, wavelet function in this paper, but the basic details of our algorithm would be unchanged if we use other wavelet functions, such as the Haar or Morlet wavelet functions. Aspects of our algorithm include : (1) the computation of the correlation of the wavelet function and the data image using either analytic or FFT methods ; (2) the computation of a local, exposurecorrected normalized (i.e., ÿat-ïelded) background in each pixel ; (3) its applicability within the low-counts regime, as it does not require a minimum number of background counts per pixel for the accurate computation of source detection thresholds ; (4) the correction of those correlation values which are aected by large exposure variations within the wavelet support (due to, e.g., telescope support ribs or the edge of the ïeld of view), using either one of two methods given by equations (16) and (17) ; (5) the generation of a source list in a manner that does not depend upon the details of the PSF shape, including the creation of general, data-dependent source cells for the estimation of source properties and ïltering of extended source data ; and (6) full error analysis. In these respects, our algorithm is considerably more general than the similar wavelet-based methods developed by Vikhlinin et al. (1994), Rosati et al. (1995), Grebenev et al. (1995), Damiani et al. (1997a, 1997b), Starck & Pierre (1998), and Lazzati et al. (1998). Nearly all of these methods were developed speciïcally for treating data collected by the ROSAT PSPC (Starck & Pierre 1998 apply their method to data from the ROSAT HRI) ; none except Damiani et al. attempt to correct for variations in exposure within the FOV ; and all except for Damiani et al. assume a

ÿat background across the region of interest. The relatively more general method of Damiani et al. is limited to analyzing data from detectors which have PSFs with Gaussian shape, and which have high rates of background accumulation. These limitations make the Damiani et al. approach, as published, inappropriate for use with, e.g., Chandra data. In ° 4, we demonstrate the robustness of our algorithm by applying it, without algorithmic changes, to data collected by an idealized detector with a spatially invariant Gaussian PSF ; to ROSAT PSPC data of the crowded ïeld of the Pleiades Cluster ; and to a simulated Chandra ACIS-I image of the Lockman Hole region. Collectively, these test cases indicate that our algorithm : (1) eectively detects and describes point sources and can be applied to the study of extended sources ; (2) does not detect more spurious sources than expected ; (3) is more sensitive than sliding-cell methods ; and (4) has equal sensitivity to the Damiani et al. method for the speciïc case of ROSAT PSPC data. We ïnd that while we can use the algorithm as presented to analyze extended sources, such analysis requires careful monitoring on the part of the analyst. Work on generalizing this analysis is on-going and will be reported in a future work. The authors acknowledge the support of the CXC Beta Test Site grant at the University of Chicago and NASA grants NAG5-3173, NAG5-3189, NAG5-3195, NAG53196, NAG5-3831, NAG5-6755, and NAG5-7226. We are grateful to T. Gaetz and A. Vikhlinin for making available the simulated ACIS-I image of the Lockman Hole, analyzed in ° 4.3. We would like to thank R. C. Nichol, B. Holden, J. Flanagan, T. Calderwood, and A. Dobryzcki for their assistance with coding issues, E. Kolaczyk for useful discussions regarding the computation of covariance, and F. Damiani, S. Sciortino, and A. Bijaoui for useful discussions. Finally, we would like to thank the anonymous referees for their detailed comments and suggestions, which greatly helped the preparation of this manuscript. In particular, our description of the vanishing moments property of wavelets and its relationship to source detection owes much to one of these referees.

APPENDIX A DERIVATION OF FORMULAE ASSOCIATED WITH THE MEXICAN HAT FUNCTION The source detection algorithm that we have presented does not explicitly depend on the details of wavelet function itself, and thus should be applicable with other, simple, wavelet functions which have one central positive mode. Hence, we have deferred to this Appendix and the next the derivation of various formulae that we use in our algorithm which are valid speciïcally when using the Marr, or Mexican Hat (MH), wavelet function.
A1.

PIXEL-BY-PIXEL INTEGRATION OF THE MH FUNCTION

Before carrying out source detection at given scales (p , p ), we must determine the grid of values W { { (see eq. [7]). This i~i , j~j is accomplished by integrating the function W (x/p , yxp ), deïned in equations (4)õ(5), on a pixel-by-pixel basis within an /y mwithxvalues farther from the origin set to zero : y ellipse with axes with full lengths 10p and 10p , x y x2 y 2 xy dx dy W , (A1) W { {\ mp p i~i , j~j x1 y1 xy x2 y2 x2 y2 x2 y2 dx dy 2 [ [ exp [ exp [ . (A2) \ p2 p2 2p2 2p2 x1 y1 x y x y (x, y) represent coordinates relative to the center of pixel (i, j), and (x , x ) and (y , y ) denote the limits of integration over 12 12

PP A PPA

B

BA

BA

B


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

211

FIG. 21.õIllustration of variables used when computing correlation values. x and y are continuous variables describing the wavelet function centered in pixel (i, j). The correlation value for pixel (i, j) is computed by summing the product of the data in pixels (i@, j@) and the integral of the wavelet function in those pixels. The dotted line shows the boundary between the positive kernel PW and negative annulus NW of a Mexican Hat wavelet with p \ p \ 1 pixel, x y while the solid line shows the extent over which the summation is carried out ; beyond this line, the amplitude of wavelet function is B0.

pixel (i@, j@) (Fig. 21) : x \ i@ [ i [ 1 , y \ j@ [ j [ 1 , 1 2 1 2 x \ i@ [ i ] 1 , y \ j@ [ j ] 1 . 2 2 2 2 We expand equation (A2) and perform each integral separately. The ïrst integral in the expansion is computed using

P

x2 x2 t2 dx exp [ \ J2p dt exp ([ x 2p2 x1 x t1 t2 \ dt exp ([t2) [ 0 n x 2 p erf \ 2x J2p x n \ px , 2 x erf

A

B

P

P P SCA B S B P

t2) t1 dt exp ([t2) 0 x 1 [ erf J2p x

A BD

where the substitution t \ x/J2p is made, and erf(x) is the error function. Performing the same integration over y, we x determine that the ïrst term is 2(n/2)p p x y \ np p x y . x y erf erf To compute the second and third integrals, we must x y erf erf e.g., determine,

P

x

2

dx

x1

x2 x2 t2 exp [ \ 2J2p dt t2 exp ([t2) . x p2 2p2 x x t1

A

(A3)


212

FREEMAN ET AL.

Vol. 138

This integral can be solved by parts, by taking the derivative of t and the integral of t exp ([t2). We present the solution : x2 x2 x exp [ 1 [ x exp [ 2 ] 2 1 2p2 2p2 x x The ïnal solution is W i~i{, \ np p x y [ j~j x y erf erf
{

A

B

A

BS S

n p x \x ] diff 2 x erf

S

n px . 2 x erf

(A4)

S

n py x ] 2 y erf diff

A

n px [ 2 x erf

BS

n px y ] 2 x erf diff

A

S

n py . 2 y erf

B

(A5)

A2.

FOURIER TRANSFORM OF THE MH FUNCTION

Analytic computation of correlation values (eq. [7]) may be too computationally intensive if the image and/or wavelet scale sizes are too large. So in WAVDETECT, for instance, fast Fourier transforms20 (FFTs) are used if scale sizes are º2 pixels : C B FFT~1[N ] FFT(W ) ] FFT*(D)] . (A6)

To mitigate the eect of "" wrap-around,îî any image that is to be transformed is padded with zeros ; the minimum width of the padding is 10 ] max (p , p ). xy In our own FORTRAN code, we use the analytic Fourier transform of the MH function, which we denote FT(W ): dx dy W (x, y)e2ni(kx x`ky y) ~= ~= `= `= `= `= dy e2niky y dx W (x, y) cos (2nk x) ] i dy e2niky y dx W (x, y) sin (2nk x) (A7) \ x x ~= ~= ~= ~= `= `= \ dy e2niky y dx W (x, y) cos (2nk x) x ~= ~= `= `= `= dy cos (2nk y) ] i dy sin (2nk y) dx W (x, y) cos (2nk x) (A8) \ y y x ~= ~= ~= `= `= dy cos (2nk y) dx W (x, y) cos (2nk x) \ y x ~= ~= `= dy cos (2nk y) ] C . \ y W ~= The wavenumber k equals [2(i [ 1)k ]/l, where i is the pixel number in Fourier space, k \ 1 is the Nyquist wavenumber, k N k N2 and l is half the length of the relevant axis in the padded image. The fourth integral in equation (A7) and the second integral in equation (A8) are zero, as the integrands are products of even and odd functions. Solving for C , we ïnd W `= y2 x2 y2 x2 dx 2 [ [ exp [ cos (2nk x) C \ exp [ x W 2p2 p2 p2 2p2 y ~= x y x y2 y2 1 1 1 2[ ( 2nk , ,0 [ ( 2nk , . (A9) \ exp [ x J2p x J2p 2p2 p2 0, c p2 2, c y y x x x ( (p, q, j) and ( (a, p) represent two integral solutions which one may ïnd, e.g., in Gradshteyn & Ryzhik (1980 ; 0, c 2, c formulae [3.896-2] and [3.952-4], respectively) : FT(W ) \

PP P P CP P P

`= `=

P P

P

P

P

DP

P

A A

BP A BAB BCA B A B A
p2 1 cos (pj) , (p, q, j) \ Jn exp [ 0, c 4q2 q a2 2p2[ a2 . exp [ 4p2 4p5

BD

( and

A

B

(A10)

(

(a, p) \ Jn 2, c

A

B

(A11)

Substituting these quantities into equation (A9), one ïnds that

y2 y2 [ 2n2k2 p2 p J2n 4n2k2 p2] 1 [ . C \ exp [ xx x xx W 2p2 p2 y y

A

B

C

D

(A12)

20 The CIAO WAVDETECT code (written in C) makes use of the FFTW package (Frigo & Johnson 1998), while our own FORTRAN code uses the publicly available MFFT algorithm (Nobile & Roberto 1986).


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

213

Substituting equation (A12) into equation (A8) and solving, we ïnd FT(W ) \ exp ([2n2k2 p2)p J2n xxx

`= y2 y2 dy cos (2nk y) exp [ 4n 2k2 p2] 1 [ y xx 2p2 p2 ~= y y 1 1 1 \ exp ([2n2k2 p2)p J2n (4n2k2 p2] 1)( (2nk , , 0) [ ( 2nk , xxx xx 0, c y J2p 2, c y J2p p2 y y y \ (2n)3p p (k2 p2] k2 p2) exp [[2n2(k2 p2] k2 p2)] . (A13) xyxx yy xx yy By using FT(W ) instead of FFT(W ), we reduce computation time, but numerical estimates are less accurate (with respect to estimates derived analytically). We quantify the discrepancy between any two of our three methods (FT(W ), FFT(W ), or analytic) using

P C

A

BC

D

A

BD

o *C o . (A14) T an We use the analytically derived detection threshold T in the denominator because the expectation value of C is zero. It also an provides an intuitive way to describe the discrepancy at the detection threshold. We ïnd d D 10~4, while d Z 10~2, FFT, an FT, an increasing with C/T . (See Fig. 22.) The discrepancy an results from the fact that we analyze binned images. Because the data are binned, we should actually compute FT(W ) using the equation d\ dx dy W { { e2ni(kx x`ky y) , (A15) i~i , j~j ~= ~= If W (x, y) exhibits signiïcant curvature in a bin, then W { { D W (x, y) within that bin ; the eect of this discrepancy in the ïnal correlation value is proportional to the number ofi~i , j~j in the bin. Thus, d counts will increase with source strength, as FT, an shown in the right panel of Figure 22. If we wish to compute the error of the correlation value in a pixel, we calculate SW 2 % DT , saving computation time by i, j using the analytic Fourier transform FT(W 2). The derivation of this function is similar to the derivation of FT(W ). A new integral which appears is FT(W ) \ `= x4 x2 dx exp [ cos (2nk x) . (A16) x p4 p2 ~= x x One may solve this integral by parts, dierentiating the term (x3/p2) cos (2nk x) and integrating the term (x/p2) x x x

PP

`= `=

P

AB

FIG. 22.õL eft : discrepancy between the correlation values determined by the analytic method and by using the FFT, for the ROSAT PSPC data of the Pleiades cluster. In this example, p \ p \ 2 pixels. Right : same as left, but comparing the analytic method with the use of the analytic Fourier transform for x y the MH function.


214 exp [[(x2/p2)]. The integral x

FREEMAN ET AL.

Vol. 138

x2 dx x3 sin (2nk x) exp [ x p2 ~= x has solution (see, e.g., Gradshteyn & Ryzhik 1980, formula [3.952-5])

P

`=

AB

(A17)

P

The ïnal solution is

x2 3n(k /p2) [ 2n3k3 xx x exp ([n2p2 k2) . dx x3 sin (2nk x) exp [ \ Jnp7 x x xx p2 2 ~= x

`=

AB

(A18)

FT(W 2) \ [2np p ] n5p p (k2 p2] k2 p2)2] exp [[n2(k2 p2] k2 p2)] . (A19) xy xyxx yy xx yy Use of FT(W 2) instead of FFT(W 2) leads to a less accurate accounting of the errors as derived using purely analytic methods, for reasons given above. We ïnd the magnitude of the discrepancy in error estimates, as a fraction of the detection threshold, to be similar to that seen above (D10~2) ; if we use the analytic error in the denominator instead of the detection threshold, we ïnd that the discrepancy is constant as a function of C/T ,at D10~2.
A3.

INTEGRATION OF THE MH FUNCTION NEGATIVE ANNULUS

If exposure variations and the FOV edge may be ignored in the computation of the background, then we may replace SNW % ET in equation (12) with a background normalization factor, N , which we derive here. i, j B The average value of the MH function is zero. Hence, N may be derived by integrating the MH function over either PW or B NW . We choose the former : py*2~(x2@p2)+1@2 x x2 y2 x2 y2 py dy 2 [ [ exp [ exp [ . (A20) p2 p2 2p2 2p2 ~S2px ~*2~(x2@p2)+1@2 x x y x y The limits of integration reÿect that the core extends over an ellipse with axes J2p and J2p . We reparameterize the x integral using polar coordinates (r, h), after ïrst mapping the boundary ellipse to a boundary circley using the transformation y@ \ (p /p )y : xy px*2~(x2@p2)+1@2 x x2] y@2 x2] y@2 p S2px dx dy@ 2 [ exp [ (A21) N\y Bp p2 2p2 ~px*2~(x2@p2)+1@2 x x x x ~S2px S2px r2 r2 p 2n dr r 2 [ exp [ (A22) dh \y p2 2p2 p x x 0 x0 S2px r2 S2px r3 r2 p dr 2r exp [ [ dr exp [ . (A23) \ 2n y 2p2 p2 2p2 p x x x 0 x0 The determinant of the Jacobian of the transformation from (x, y@) ] (r, h)is r. We evaluate the second integral using integration by parts : N\ B S2p
x

P

dx

P

C

DA

BA

B

PP C DA PP C DAB CP A BP A A BC A BK D P
x

B

BD A

S2px r3 r2 r2 S2px r2 exp [ \[ r 2 exp [ ] dr2r exp [ . p2 2p2 2p 2p2 0 x x x0 0 x The second integral in equation (A24) cancels with the ïrst integral in equation (A23), leaving dr r2 S2p p N \ 2n y r2 exp [ B 2p p x0 x p 2p2 x \ 2n y p exp (1) x 4np p x y. \ exp (1)

P

S2px

B

(A24)

C

A

BK D

(25)

APPENDIX B COMPUTATION OF DETECTION THRESHOLDS We associate an image pixel (i, j) with a source if the signiïcance user-deïned detection threshold signiïcance S : o = dCp(C o B ) [ S\ i, j i, j Ci, j

P

S of its correlation value C is greater than a i, j i, j S, o (B1)


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

215

where B is either SNW % DT , the convolution of the wavelet negative annulus NW with (possibly cleansed) image data D, i, j i, j or the amplitude in an input background map. p(C o B ) is the probability sampling distribution (PSD) for C given B . (We i, j i, j are assuming here that the ïeld of view is evenly exposed, which is true for the simulations we describe below. In a more general expression, we would replace SNW % DT with E SNW % DT /SNW % ET .) By deïning the signiïcance in this i, j i, j i, j i, j manner, we are making the assumption that the expected background counts amplitude is constant throughout the wavelet support and that if we are computing the background, there are no source counts in the negative annulus. Since the PSD is dimensionless, it is independent of scale size (e.g., the PSD for B would be the same if we were to double p and p and i, j x y reduce the count rate by a factor of 4). We replace B with the equivalent quantity q , the expected number of background counts within the spatial region i, j i, j spanned by the positive kernel of the wavelet, PW . (q is related to B by multiplicative factors.) The PSD does not have i, j i, j analytic form when the total number of expected counts within the positive kernel of the wavelet, q , is small ([ 1) ; this is i, j demonstrated by Damiani et al. (1997a). Like Damiani et al., we use simulations to estimate signiïcances and detection thresholds ; we carry these simulations out both within the regime they examine and also for lower expected counts values. Because this is a computationally intensive problem, we cannot accurately estimate signiïcances for each pixel if S [ 10~7. i, j Instead, we compare the value of C in each pixel with the detection threshold C (S , B ), deïned by i, j i, j, o o i, j = dCp(C o B ) S\ i, j o Ci, j, o = dCp(C o q ) . \ i, j Ci, j, o Here we recast the equation using the variable q , the number of expected counts in PW , as its use clariïes our description i, j below of how we estimate C . i,),, we simulated over 50,000 1024 ] 1024 ÿat-ïeld images.21 For each image, we : jo To determine C (S , q i, j, o o i, j 1. Randomly selected values log q from the range [10 ¹ log q ¹ 3.25 (we describe why we chose this upper limit below) ; o o 2. Determined the expected background amplitude B \ q /2np2 in each image pixel (p \ 4 pixels) ; o o distribution given rate B ; 3. Sampled data in each pixel from the appropriate Poisson o 4. Computed C and q for each image pixel ; and i, j in bins of size *(log q ) \ 0.2, for [6.9 ¹ log q ¹ 3.1, with one bin being used for all values of i, j 5. Recorded C j iC i log q \ [6.9.22i,From these distributions p,(j o q ), we can determine C, j . i, j i, j i, j, o Because there is an inverse correlation between observed values of C and q , it is important to record p(C o q ) and not i i, j i, j p(C o q ). Use of the latter distribution leads to underestimated detection, jthresholds, and thus to larger numbers of false source o detections than one would expect, given S . We determined 25 values of C (S , q) inoeach log q bin, for values of S Z 10~7, using the central 68% (17 values) to estimate oo o "" 1 p îî errors on C . We then ïtted these data with simple functions, minimizing the s2 statistic. These functions we use o detection thresholds well, except in the regime log q [ 0 and log S Z [4, where we use a look-up describe the observed i o table instead. These ïts allow us in principle to compute detection threshold, jfor signiïcances below our computational lower limit S D 10~7 (such as for S D 10~9, the signiïcance corresponding to one false source pixel in an AXAF HRC image). o o In the regime log q [ 0 and log S [ [4, we compute log C (q ) using the function i, j o o i, j log C (q ) \ A (log q )2] B (log q ) ] C , (B2) o i, j lo i, j lo i, j lo where

P P

A \ 0.00462 , lo B \ 0.0661 , lo C \[0.0154(log S )2[ 0.252(log S ) [ 0.031 . lo o o In the regime log q Z 0, we compute C (q ) using the function i, j o i, j C (q ) \ A Jq ] B , o i, j hi i, j hi

(B3)

21 Images that had at least one count ; empty images were ignored. 22 This is of the order of the machine precision.


216 where

FREEMAN ET AL.

Vol. 138

4 [0.509(log So) ] 1.897 [ 0.00172(log So ] 7)3.606 , log So º [7 , 0 A \5 hi 6 [0.509(log S ) ] 1.897 , log S \ [7 , 0 o o B \[1.115(log S ) [ 1.038 . hi o This function is also used by Damiani et al. to ït detection thresholds in this counts regime, though their derived coefficients dier from ours. For values log q D 0, we ïnd that we must use the formula i, j log C (q ) \ A (log q ) ] B (B4) o i, j mid i, j mid to compute the detection threshold, with
A \ 0.00182(log S )3] 0.0279(log S )2] 0.158(log S ) ] 0.607 mid o o o and S ) ] 0.612 [ 0.00015(log S ] 7)4.75 , [2 \ log S ¹ [1 , o o o S ) ] 0.612 [ 0.00085(log S ] 7)3.5 , [7 ¹ log S ¹ [2 , (B6) o o o S ) ] 0.612 , log S \ [7 . o o For our simulations, we chose log q \ 3.25 as a upper limit because of the contention of Damiani et al. that if log q Z 3, o i, j the PSD probability sampling distribution is analytically representable as a Gaussian with width p \ (q )1@2. (Note that the i,this is the case, the j value of q that we use in this work is larger than that used by Damiani et al. by a factor of 2n.) If i, signiïcancejis given by B C2 = 1 dC exp [ i, j S\ i, j J2nq 2q i, j i, j Ci, j Ci, j 1 1 C2 \1[ [ dC exp [ i, j 2 J2nq 2q i, j i, j 0 C 1 i, j \ 1 [ erf . (B7) 2 J2q i, j We ïnd, however, that if we use this formula in the regime log q Z 3, the derived values of C are smaller than those i, j i, j, o predicted by equation (B3) above. Thus, because it is more conservative, we use equation (B3) to compute detection thresholds for all values log q Z 0. i, j (B5)

4 [0.064(log \ 5 [0.064(log mid 6 [0.064(log

0 0

P

A

P

B A

B

C

A

BD

APPENDIX C COVARIANCE ESTIMATE : TWO-ITERATION BACKGROUND If the data are cleansed (see ° 3.1.3), then an exact calculation of V [B ] will include nonzero covariance terms. In this i section, we derive these terms assuming that the data are cleansed only, jonce, i.e., we compute the variance for B ,a background estimate made by convolving the wavelet negative annulus with data D that are a mixture of raw data 2, i,and Dj 2 ïrst-iteration background estimates B : 1 SNW % D T 2 i, j \ N ; ; NW { { D@{ { , (C1) B \E i, j @ @ i~i , j~j i , j 2, i, j i, j SNW % ET i, j ij where N \ E /SNW % ET . i, j ,j i, j Using iequation (18), V [B ] \ ; ; (N NW { {)2V [D@{ {] ] 2 ; ; ; ; (N NW { {)(N NW {{ {{)cov[D@{ { , D@{{ {{] . 2, i, j i, j i~i , j~j i ,j i, j i~i , j~j i, j i~i , j~j i ,j i ,j i@ j@ i@ i@@;i@ j@ j@@;j@ To calculate the variance, we must estimate both V [D@{ {] and cov[D@{ { , D@{{ {{]. The estimate of the former is given in , ij i ,j equation (22), derived assuming that the raw datum D { { iisj sampled from ,a Poisson distribution with variance D { { and that i ,j i ,j each pixelîs raw datum is independently sampled. Estimation of cov[D@{ { , D@{{ {{] is more complicated. For the two-iteration background case, there are three possibilities : ,j i D@{ { \ D { { and D@{{ {{i \ D {{ , j ; D@{ { \ D { { and D@{{ {{ \ B {{ {{ (or vice-versa) ; or D@{ { \ B { { and D@{{ {{ \ B {{ {{ . j , j{{ j i ,j i ,that the raw datai are independently sampled, in the iïrst case the covariance isi ,zero. Ini the second ,case, we i can i ,j i ,j i ,j i ,j 1, , j 1, , j ij 1, , j Assuming


No. 1, 2002

SPATIAL ANALYSIS OF POISSON DATA

217

FIG. 23.õIllustration of the eect of including the computation of covariances in the calculation of a two-iteration background (see Appendix C for details). L eft : ratio of variances V [B ] /V [B ] , for a 50 ] 50 subïeld of the ROSAT PSPC Pleiades Cluster image, analyzed with a wavelet with covar 2 values scale sizes p \ p \ 4 pixels. Light 2regions havenocovar B1, while dark regions have values up to B1.3. Right : 50 ] 50 subïeld of the ROSAT PSPC x for which the ratio of variances shown at left was computed. Sources in this ïeld correspond with light (low-amplitude) regions in the ratio y Pleiades image image.

estimate the covariance using the approximation (Eadie et al., p. 27) cov[D i{, j{ ,B i{{, j{{

LD { { LB {{ {{ i ,j 1, i , j cov[D , D { {] ]\; ; ; ; 1, k, l k , l Lk Lk { { @ l@ klk k, l k ,l LD { { LB {{ {{ i ,j 1, i , j V [D ] \; ; k, l Lk Lk k, l k, l kl LD { { LB {{ {{ i ,j 1, i , j V [D ] \; ; k, l LD LD k, l k, l kl LB {{ {{ 1, i , j V [D { {] \ i ,j LD { { i ,j L N ; ; NW {{ {{{ {{ {{{ D {{{ {{{ \ V [D { {] i ~i , j ~j i , j i , j LD { { i{{, j{{ @@@ @@@ i ,j ij \ V [D { {]N {{ {{ NW {{ { {{ { i ,j i ,j i ~i , j ~j \ D { { N {{ {{ NW {{ { {{ { . i ,j i ,j i ~i , j ~j In the above equations, k represents the expectation value of the sampling distribution for D ; for a Poisson distribution with variance equal to the datum, k will be equal to the datum. In the third case, we ïnd (while skipping some steps) LB { { LB {{ {{ 1, i , j 1, i , j V [D ] cov[B { { , B {{ {{] \ ; ; k, l 1, i , j 1, i , j LD LD kl k, l k, l N V [D ] \ ; ; N { { NW { NW {{ k, l i ,j i ~k, j{~l i{{, j{{ i ~k, j{{~l kl \ ; ; N { { NW { N NW {{ D. i ,j i ~k, j{~l i{{, j{{ i ~k, j{{~l k, l kl To demonstrate the eect of including covariance terms, we compute V [B ] for a 50 ] 50 subïeld of the ROSAT PSPC Pleiades image, with p \ p \ 4 pixels. In Figure 23, we show the ratio V [B 2] /V [B ] . We infer the following from y 2 covar 2 nocovar this computation : (1) itxis quantitatively unimportant at the location of sources and within voids, where the change in variance is [ 1% ; (2) in the vicinity of a strong sources, the change in variance is D10% (with maximum B30% for this particular ïeld, in the vicinity of an 800 count source) and is aected by source crowding ; and (3) including covariance terms increases the CPU time needed to compute V [B ] by a factor D O(d d p2 p2), where d and d are the x- and y-axis lengths, 2 xyxy x y respectively, in pixels.

A BA B A BA B A BA B AB

A

BA

B


218

FREEMAN ET AL.

REFERENCES Bahcall, J. N., Flynn, C., Gould, A., & Kirhakos, S. 1994, ApJ, 435, L51 Lazzati, D., Campana, S., Rosati, P., Chincarini, G., & Giacconi, R. 1998, Biviano, A., Durret, F., Gerbal, D., Le Fevre, O., Lobo, C., Mazure, A., & A&A, 331, 41 Slezak, E. 1996, A&A, 311, 95 Lazzati, D., & Chincarini, G. 1998, A&A, 339, 52 Coupinot, G., Hecquet, J., Auriere, M., & Futaully, R. 1992, A&A, 259, 701 Mallat, S. 1998, A Wavelet Tour of Signal Processing (London : Academic Damiani, F., Maggio, A., Micela, G., & Sciortino, S. 1996, in Press) Rontgenstrahlung from the Universe, ed. H. U. Zimmerman, J. Micela, G., Sciortino, S., Kashyap, V., Harnden, F. R., Jr., & Rosner, R. Trumper, & H. Yorke (MPE Report 263), 679 1996, ApJS, 102, 75 õõõ. 1997a, ApJ, 483, 350 Micela, G., et al. 1999, A&A, 341, 751 õõõ. 1997b, ApJ, 483, 370 Nichol, R. C., Holden, B. P., Romer, A. K., Ulmer, M. P., Burke, D. J., & Daubechies, I. 1992, Ten Lectures on Wavelets (Philadelphia : SIAM) Collins, C. A. 1997, ApJ, 481, 644 Eadie, W. T., Drijard, D., James, F. E., Roos, M., & Sadoulet, B. 1971, Nobile, A., & Roberto, V. 1986, Comput. Phys. Commun., 42, 233 Statistical Methods in Experimental Physics (Amsterdam : NorthNousek, J. A., & Lesser, A. 1993, ROSAT Newslett., 8, 13 Holland) Pierre, M., & Starck, J.-L. 1998, A&A, 330, 801 Freeman, P. E., Kashyap, V., Rosner, R., & Lamb, D. Q. 2002, in Statistical Rosati, P., Della Cecai, R., Burg, R., Norman, C., & Giacconi, R. 1995, ApJ, Challenges in Modern Astronomy III, ed. G. J. Babu & E. D. Feigelson 445, L11 (New York : Springer), in press Schmidt, M., et al. 1998, A&A, 329, 495 Freeman, P. E., Kashyap, V., Rosner, R., Nichol, R., Holden, B., & Lamb, Slezak, E., Bijaoui, A., & Mars, G. 1990, A&A, 227, 301 D. Q. 1996, in Astronomical Data Analysis Software and Systems V, ed. Slezak, E., Durret, F., & Gerbal, D. 1994, AJ, 108, 1996 G. H. Jacoby & J. Barnes (San Francisco : ASP), 163 Snowden, S. L., & Kuntz, K. D. 1998, Cookbook for Analysis Procedures Frigo, M., & Johnson, S. G. 1998, Proc. ICASSP, 3, 1381 for ROSAT XRT Observations of Extended Objects and the Diuse Frisch, P. C. 1995, Space Sci. Rev., 72, 499 Background, Part I : Individual Observations, Technical Report, NASA/ Gill, A. G., & Henrikson, R. N. 1990, ApJ, 365, L27 GSFC Gradshteyn, I., & Ryzhik, I. 1980, Table of Integrals, Series, and Products Starck, J.-L., & Murtagh, F. 1998, PASP, 110, 193 (San Diego : Academic Press) Starck, J.-L., Murtagh, F., & Bijaoui, A. 1995, in Astronomical Data Grebenev, S. A., Forman, W., Jones, C., & Murray, S. 1995, ApJ, 445, 607 Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Harnden, F. R., Jr., Fabricant, D. G., Harris, D. E., & Schwarz, J. 1984, Hayes (San Francisco : ASP), 279 SAO Spec. Rep. 393 Starck, J.-L., & Pierre, M. 1998, A&AS, 128, 397 Hasinger, G., Burg, R., Giacconi, R., Schmidt, M., Trumper, J., & ZamoStauer, J. R., Caillault, J.-P., Gagne, M., Prosser, C. F., & Hartmann, rani, G. 1998, A&A, 329, 482 L. W. 1994, ApJS, 91, 625 Hasinger, G., Johnston, H., & Verbunt, F. 1994, A&A, 288, 466 Tinney, C. G., Mould, J. R., & Reid, I. N. 1992, ApJ, 396, 173 Holschneider, M. 1995, Wavelets : An Analysis Tool (Oxford : Oxford Ulmer, M. P., Romer, A. K., Nichol, R. C., Holden, B., Collins, C., & Burke, Univ. Press) D. 1995, BAAS, 187, 95.03 Kashyap, V., Micela, G., Sciortino, S., Harnden, F. R., Jr., & Rosner, R. Vikhlinin, A., Forman, W., & Jones, C. 1994, ApJ, 435, 162 1994, in The Soft X-Ray Cosmos, ed. E. M. Schlegel & R. Petre (New Vikhlinin, A., McNamara, B. R., Forman, W., Jones, C., Quintana, H., & York : AIP), 239 Hornstrup, A. 1998, ApJ, 502, 558 Kashyap, V., Rosner, R., Harnden, F. R., Jr., Micela, G., & Sciortino, S. Voges, W., Gruber, R., Haberl, F., Kuerster, M., Pietsch, W., & Zimmer1996, in Cool Stars, Stellar Systems, and the Sun, ed. R. Pallavicini & A. mann, U. 1994, ROSAT News 32 K. Dupree (San Francisco : ASP), 365 White, N. E., Giommi, P., & Angelini, L. 1994, IAU Circ. 6100 Langer, W. D., Wilson, R. W., & Anderson, C. H. 1993, ApJ, 408, L45