Source: http://cxc.harvard.edu/ciao/download/doc/detect_manual/wav_theory.html
Timestamp: 2019-04-25 16:28:43+00:00

Document:
This chapter is a draft copy of the paper "A Wavelet-Based Algorithm for the Spatial Analysis of Poisson Data'' by P.E. Freeman, V. Kashyap, R. Rosner, D.Q. Lamb, Ap. J. Suppl. Series, Vol. 138, p. 185, 2002. The published paper is available online via astro-ph: http://arXiv.org/abs/astro-ph/0108429 .
may be used as a wavelet function. is the scaling or dilation parameter and is the translation parameter.
The importance of wavelet functions only began to be recognized in early 1980s, when they were used in the analysis of seismic data (Goupillaud et al. 1984). Wavelets are now used as tools in atomic physics, the study of fractal structures, time series, sound, and image analysis, etc. Slezak, Bijaoui, & Mars (1990) were the first to apply wavelets in astronomy, applying them to galaxy catalogue data to search for statistically significant clusters and voids. The interested reader may find more information on wavelets in Daubechies (1992) or Holschneider (1995).
Wavelet transforms, convolutions of an arbitrary analyzable function with a set of suitably scaled wavelet functions, provides a superior means to detect and characterize sources in an image, compared with classical "sliding-box''13.1 and Fourier methods. Because a wavelet is a function with limited spatial extent, it can be used to determine source locations, like the sliding-box. Unlike the sliding-box, however, it can also be used to identify the dominant frequencies in the analyzed function, i.e. to characterize source shape and extent. Such characterization could also be done with Fourier analysis, but unlike Fourier methods, wavelets determine localized frequency information, as a consequence of their localization in both the spatial and Fourier domains.
The program wtransform applies wavelets to the problem of detecting sources in X-ray data, and the associated program wrecon analyzes the detected sources and constructs the source list. While others have developed codes with similar analysis goals (e.g. Vikhlinin, Forman, & Jones 1994; Rosati et al. 1995; Grebenev, Forman, & Jones 1995), wtransform and wrecon are considerably more complex, e.g. correcting for exposure variation across the field of view (FOV), and making error estimates, while not requiring that the point spread function (PSF) have a particular functional form (which is particularly important for analysis of Chandra images). Another program of similar complexity but using different methods is the palermo wdetect program of Damiani et al. (1997).
In §13.2 we describe the basic recipe one would follow to detect sources in data, without reference to any particular wavelet function. We do this is to underscore the fact that the basic details of the method are not dependent on the choice of wavelet function, so long as that function is simple, with one centralized positive amplitude mode. The fact that the PSF of X-ray detectors often has Gaussian-like shape motivates our use of the Marr wavelet, or "Mexican Hat" (MH) function, which we describe in Appendix 13.5.
The program wtransform uses the correlation of the MH function with a given image to identify putative source pixels. In §13.3 we describe in detail the steps that wtransform follows, in particular discussing how the program estimates the local background at each pixel, how it corrects for the systematic shift that a pixel's correlation value will take if the exposure varies within the limited spatial extent of the wavelet function, and how it performs error analysis. In appendices which relate to §13.3, we present derivations of analytic quantities related to the MH function (Appendix 13.6), as well as describe the simulations used to estimate the threshold correlation value for source detection in each pixel (Appendix 13.7). In §13.4 we discuss in detail how the source list is constructed by wrecon , and in an appendix we describe how wrecon computes a default background image for use in determining source properties (Appendix 13.8).
(Note that the data are assumed to have units counts, and not counts per second.) The expectation value of is zero, if there are no sources within the limited spatial extent of the wavelet function, and the background count rate is locally constant, because the normalization of is zero.
A clump of counts will manifest itself as a local maximum in a correlation image if the scale sizes are approximately the same as, or larger than, the size of the clump. One can see why this is so more easily if we write , where and denote the positive and negative amplitude portions of the wavelet function, respectively. If the clump is contained completely within , then the contribution of the positive term to outweighs that of the negative term, producing a maximum.13.3 If the scale sizes are smaller, then the wavelet function will extend over a smaller region within the clump, within which the inferred counts amplitude will tend more and more to a constant. Hence . For larger scale sizes, the correlation value tends asymptotically to a maximum, .
One continues to iterate until is determined. (The usual number of iterations depends upon many factors, but is usually 3-4.) With this final background estimate, one computes a final significance for each value (the correlation of the wavelet function with the raw image data) so that a final listing of source pixels may be made.
After this algorithm is used to determine lists of source pixels for many wavelet scale sizes, cross-identification of pixels across scales is performed to create the final source list. This cross-identification allows the weeding out of spurious sources which may be detected even if a high significance threshold is used.
is an estimate of the number of counts that would be detected if there were no exposure variations; it differs from the number of photons seen from a similarly sized region of sky by a factor related to the efficiencies of the mirror assembly and detector. Note that is set to zero in pixels with exposure less than a user-specified threshold.
We stress that one should not interpret blindly as the actual background amplitude ; it should only be viewed as a computation made in order to select source pixels. For instance, if a small wavelet scale is applied in a region with a large PSF, then around any source, because the estimate of is unavoidably biased by the source counts which lie in the region spanned by . In wrecon , we provide a default method of combining information across scales to estimate (which we discuss in Appendix 13.8); it is this quantity which is used to help determine source properties.
The wavelet is centered in pixel (e.g. = [0.5,0.5] for pixel = [1,1]), and and represent the limits of integration over pixel . The integral in Eq. (13.6) can be solved analytically if the MH function (Eq. [13.37]) is used; we present this solution in Appendix 13.6. For the MH function, summation over and within an ellipse with semi-major and minor axes of length 5 is sufficient to determine to high accuracy.
13.3.2.0.2 Calculation using the Fast Fourier Transform.
The right side of this equation may be read as follows: "the inverse FFT of the product of a normalization factor , the FFT of , and the complex conjugate of the FFT of ." To avoid problems associated with FFT aliasing, we pad the image with zeros; an axis-length in the padded image is, in general, the next power of two larger than the length of that axis in the unpadded image (e.g. for an image of size 1024 1024, the padded image has size 2048 2048, with the contents of the unpadded image in the center). We compute the FFT using realft and associated routines of Press et al. (1992). For the particular case of the MH function, we may compute the Fourier Transform analytically; we present the derivation in Appendix 13.6.
In this example, 0. Thus exposure variations can lead to the detection of spurious sources near shadows (but not within them) and near the FOV edge. Below, we present two methods for correcting the estimate .
The source counts, , may also be altered by variations in exposure if the source is extended; however, correcting requires knowledge of , which is a priori unknown. Also, our prime motivation for correcting is to reduce the number of spurious sources, and correcting does not help us in this regard.
This correction is deemed "fast" because wtransform calculates once, as opposed to once per iteration. To ensure accuracy, however, the full exposure correction is performed once the estimate of has stabilized.
In §13.2, we show how the significance is computed in each pixel. However, for reasons discussed fully in Appendix 13.7, wtransform does not compute significances directly, but rather compares each correlation value with a threshold correlation value . This threshold is found by estimating the lower bound on the integral in Eq. (13.3), given , a user-specified significance threshold (e.g. 10 ), and . If , pixel is identified as a possible source pixel.
In wtransform , the user actually specifies two threshold significance values: , for source cleansing, and . During iterations, is compared to , which may be set to be much larger than (e.g. 10 ) to help reduce the bias caused by source counts in estimates of .
If depends directly on , the covariance terms are zero, and the approximation formula simplies to the familiar formula .
The random variables have variance .
This quantity is calculated only after the data are cleansed of sources. It may be calculated analytically (Eq. [13.16]) or using FFTs (Eq. [13.17]).
13.3.3.0.2 Correlation: No Exposure Correction.
This variance is computed once, after the first iteration, either analytically (Eq. [13.18]) or through the use of FFTs (Eq. [13.19]). To speed computation when using Fourier Transforms, we use the analytically-derived quantity . Details on its derivation may be found in Appendix 13.6.
13.3.3.0.3 Correlation: Full Exposure Correction.
Initial tests indicate that this formula is accurate away from sources, but greatly overestimates the variance around sources. These tests indicate that the true variance around sources does not differ greatly from the true variance away from sources ( 50%), so that this equation can provide the analyst a good estimate of the magnitude of the variances.
and "flux images" computed within wrecon itself.
Next we discuss how the flux images are computed, then show how they are used to define the cells within the FOV that delineate the extent of a putative source and which allow us to estimate its properties.
We use as the smoothing function because it has the desirable properties of being localized, and, for the particular case of the MH function, of mimicking the shape of a canonical Gaussian PSF.
The scale pairs at which to compute the flux images are user-specified; the number of pairs can be less than the number of wavelet scale pairs examined overall. Using fewer pairs reduces computation time, but may lead to a less-precise determination of source properties. Hereafter, we assume that , forcing symmetric wavelets.
To determine source properties, it is necessary to define a "cell" on the detector, in which some portion of a source's counts are detected. If we choose the size and shape of the cell, then we need detailed knowledge of the PSF to estimate the fraction of a source's counts that are recorded in the cell. However, wrecon and the other Chandra detection algorithms assume no details about a detector's PSF except for its characteristic size , which is computed for Chandra data using the routine psf_calc_size; its functional form is unknown. But even if we did have detailed knowledge of the PSF, its use would provide a good estimate of the total source counts only if the source were point-like.
and determining those pixels for which is the local maximum in -space. Pixels for which 0 are not considered.
The use of in Eq. (13.21) makes a normalized quantity, which is beneficial since exposure variations can cause spurious flux-image maxima, which would reduce the cell size and adversely affect the determination of source properties. We stress, however, that the flux image is not used in the actual determination of source properties.
and the creation of the flux image involves smearing, or de-focusing, data, so that cell boundaries will lie outside the region where the source is evident to the eye. Hence, at least for isolated sources, nearly all source counts should have been detected within cell boundaries.
This approach is not always optimal when analyzing extended sources: such sources may have more than one local maximum in counts space, leading to multiple "sources" being listed. One way to counteract this is to set , as this will increase the probability of the extended object being associated with one maximum in -space (which will ensure a sufficiently large cell). However, it also increases the probability that point sources in the vicinity of the extended source will be "swallowed up" by the larger source cell, so care must be exercised when interpreting derived source properties such as count rate. This approach is also not optimal if PSF is bimodal, as it is far off-axis for Chandra ; in this case, wtransform will output two correlation maxima for a well-sampled source, and wrecon will characterize each maximum independently of the other. Only with the a posteriori application of a detailed PSF can the two "sources" be properly combined.
It is also possible for a source cell to contain two or more underlying sources. This is indicated when there are multiple correlation maxima in the source cell at the PSF scale. The code indicates such a condition with a flag in the source list; it does not currently attempt to "break up" the cell to refine source property estimates.
does not contain one or more sources first detected at scales smaller than the currently considered scale . This is an important check when , as sources in a crowded field merge, possibly creating "new" sources at new locations in the FOV.
After determining for the putative source its cell and properties, wrecon first scans through the lists of all correlation maxima at all scales, and marks as "examined" those which lie within the current source cell, and then verifies that the putative source appears at one or more scales . It does this by determining the scale size closest to the PSF size, (by minimizing ), then counting the number of correlation maxima that lie within the source cell for . If the number is zero, wrecon rejects the putative source.
where , and represents the source cell. We note that the estimate of location may be biased if the source cell is truncated because of the presence of a nearby source.
A preferred weighting function would be the source flux , but its use greatly complicates the computation of variances and . (For similar reasons, is not used as a weighting function.) and will give similar estimates if the source is strong or if background is small.
The missing covariance terms are positive, so Eq. (13.26) underestimates the true variance.
where and are given in eqs. (13.26) and (13.28) respectively. Equation (13.30) underestimates the true variance.
13.4.2.0.5 Background Counts in Source Cell.
13.4.2.0.6 Background Count Rate in Source Cell.
and the associated local correlation maximum is contained within a source cell.
Otherwise, = 0. 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 putative sources rejected by wrecon are also not included in the image.
The source detection algorithm we present in this paper should not depend on the details of function itself, at least for simple wavelets which have one central positive mode. In this Appendix, we describe the Marr wavelet function, or "Mexican Hat function,'' (MH) which is used by wtransform and wrecon .
The integral of over all space is, by definition, zero. This allows us to disregard the factor in Eq. (13.37) when performing correlations.
The area of this ellipse, used to compute correlation thresholds, is . The amplitude of the MH function is 2 at its center, regardless of scale; its minimum amplitude is -0.27 on the ellipse defined by axis-lengths 2 and 2 .
The MH function has several advantageous properties which motivate its use.
When rotationally symmetric, its Gaussian-like positive kernel has similar shape to a canonical PSF.
The kernel's limited extent allows only those sources with intrinsically or instrumentally broadened size to be strongly detected.
It has an analytically-derivable Fourier Transform, speeding its correlation with data via the Correlation Theorem (Appendix 13.6).
It has a limited extent in Fourier space, such that limited, discrete sampling in and (e.g. at values separated by factors of two) is sufficient to sample the entire frequency domain.
Its limited extent in both spatial and Fourier domains helps to minimize effects of aliasing.
If the width of the wavelet function is of order the size of the image pixel, then using FFTs to compute and other related quantities becomes impractible. To determine , for instance, we integrate on a pixel-by-pixel basis, as shown in Eq. (13.7).
0 if are outside an ellipse with axis-lengths 5 and 5 , limiting the range of integration.
If exposure variations and the FOV edge may be ignored in the computation of the background, then we may replace in Eq. (13.5) with a background normalization factor, which we derive here.
The determinant of the Jacobian of the transformation from is .
The wave-number equals , where is the pixel number in Fourier space, is the Nyquist wave-number, and is half the length of the relevant axis in the padded image. The fourth integral in Eq. (13.47) and the second integral in Eq. (13.48) are zero, as the integrands are products of even and odd functions.
is the inferred number of background counts within the limited spatial extent of the wavelet function, and is the probability sampling distribution for given . This sampling distribution does not depend upon the scale size of the wavelet function. We assume that the background count rate is locally constant, so that instead of , we use the expected number of background counts in the positive kernel of the wavelet, a quantity which we denote . For the MH function, .
and recording in bins of size = 0.2, for -6.9 3.1, with one bin being used for all values of -6.9.13.6 From these distributions , we can determine .
If a simulated dataset has no counts, we do not analyze it. We chose 3.25 as a upper limit because of the report by Damiani et al. (1997) that for 3, the probability sampling distribution is analytically representable as a Gaussian with width . (Note that the value of that we use in this work is larger than that used by Damiani et al. by a factor of .) We return to this point below.
By analyzing over 50,000 simulated images, we are able to determine 25 values of in each bin, for values of 10 , using the central 68% (17) of the values to estimate variances. We fit simple functions to the estimated threshold values, using the method and the estimated variances. We present our results below. These functions describe the observed threshold values well, except in the regime and , where a binned look-up table is used instead. We note that these fits allow us in principle to compute threshold correlation values for significance values below our computational lower limit 10 (such as for 10 , the significance corresponding to one false source pixel in an Chandra HRC image).
We find that if we use Eq. (13.58) in the regime 3, then the derived values of are larger than those predicted by Eq. (13.61) above. Because it is more conservative, we use Eq. (13.58) for all values 0.

References: V. 
 §13
 §13
 §13
 §13
 §13