3 Detection procedures

Without attempting to provide a review of the available techniques in the literature, we briefly describe here the procedures we have tested. They are summarized in Table 3.

 

 

Table 3: Designation, implementation and short description of the procedures or method used for detection and analysis

Procedure	Implementation	Version	Method
EMLDETECT	XMM-SAS v5.0	3.7.2	Cell detection + Maximum likelihood
VTPDETECT	Chandra CIAO	2.0.2	Voronoi Tessellation and percolation
WAVDETECT	Chandra CIAO	2.0.2	Wavelet
EWAVELET	XMM-SAS v5.0	2.4	Wavelet
G+SE	Gauss + SExtractor	2.1.6	Mixed - Gauss convolution followed by SExtractor detection
MR/1+SE	MR/1 + SExtractor	2.1.6	Mixed - multi-resolution filtering followed by SExtractor detection

3.1 Sliding cell detection and maximum likelihood (ML) method

Historically, the sliding cell detection method was first used for Einstein Observatory observations (e.g. EMSS - Gioia et al. 1990). It is included in ROSAT, Chandra and XMM-Newton data analysis tools and a good description can be found in the specific documentation for each of those missions.

The X-ray image is scanned by a sliding square box and if the signal-to-noise of the source centered in the box is greater than the specified threshold value it is marked as an object. The signal is derived from the pixel values inside the cell and noise is estimated from the neighboring pixels. Secondly, the objects and some zone around them are removed from the image forming the so-called "cheese'' image which is interpolated later by a suitable function (generally a spline) to create a smooth background image. The original image is scanned again but this time using a threshold from the estimated background inside the running cell to give the map detection object list.

The procedure is fast and robust and does not rely on a priori assumptions. However it has difficulties, especially in detecting extended features, close objects and sources near the detection limit. Many refinements are now implemented improving the sliding cell method: (1) consecutive runs with increasing cell size, (2) matched filter detection cell where the cell size depends on the off-axis angle. However, the most important improvement was the addition of the maximum likelihood (ML) technique to further analyze the detected sources.

The ML technique was first applied to analyze ROSAT  observations (Cruddace et al. 1988, 1991; Hasinger et al. 1993). It was used to produce all general X-ray surveys from ROSAT mission (e.g. RASS - Voges et al. 1999, WARPS survey - Ebeling et al. 2000). The two lists from local and map detection passes can be merged to form the input objects list for the ML pass. It is useful to feed the ML procedure with as many candidate objects as possible, having in mind that large numbers of objects could be very CPU-expensive. The spatial distribution of an input source is compared to the PSF model - the likelihood that both distributions are the same - is calculated with varying the input source parameters (position, extent, counts) and the corresponding confidence limits can be naturally computed. A multi-PSF fit is also implemented which helps in deblending and reconstructing the parameters of close sources. In the output list, only sources with a likelihood above a threshold are kept.

The ML method performs well and has many valuable features, however, it has some drawbacks - it needs a PSF model to perform the likelihood calculation and thus favours point-like source analysis, extent likelihood could be reliably taken only for bright sources, it cannot detect objects which are not already presented in the input list (e.g. missing detections in the local or map passes).

Here we have used EMLDETECT - an implementation of the method specifically adapted for XMM-SAS (Brunner 1996). In the map mode sliding cell pass we used a low signal-to-noise ratio ($\sim$$3\sigma$) above the background in order to have as many as possible input objects for the ML pass. The likelihood limit (given by $L = -\ln P$, where P is the probability of finding an excess above the background) was taken to be 10, which corresponds roughly to $4\sigma$detection. A multi-PSF fitting mode with the maximum of 6 simultaneous PSF profile fits was used.

3.2 VTP

VTP - the Voronoi Tessellation and Percolation method (Ebeling & Wiedenmann 1993; Ebeling 1993) is a general method for detecting structures in a distribution of points (photons in our case) by choosing regions with enhanced surface density with respect to an underlying distribution (Poissonian in X-ray images). It treats the raw photon distribution directly without any recourse to a PSF model or a geometrical shape of the objects it finds. Each photon defines a centre of a polygon in the Voronoi tessellation image and the surface brightness is simply the inverse area of the polygon (assuming one single photon per cell). The distribution function of the inverse areas of all photons is compared to that expected from a Poisson distribution and all the cells above a given threshold are flagged and percolated, i.e. connected to form an object. This method was successfully used with ROSAT data (Scharf et al. 1997) and is currently incorporated in the Chandra DETECT package (Dobrzycki et al. 1999).

Apart from these advantages, VTP has some drawbacks which are especially important for XMM-Newton observations: (1) because of the telescope's high sensitivity and rather large PSF with strong tails, the percolation procedure tends to link nearby objects; (2) excessive CPU time for images with relatively large number of photons; (3) there is no simple way to estimate the extension of objects.

3.3 Wavelet detection procedures

In the past few years a new approach has been extensively used: the wavelet technique (WT). This method consists in convolving an image with a wavelet function:

$$w_a (x,y) = I(x,y) \otimes \psi\left(\frac{x}{a},\frac{y}{a}\right),$$ (1)

where w_a are the wavelet coefficients corresponding to a scale a, I(x,y) is the input image and $\psi$ is the wavelet function. The wavelet function must have zero normalization and satisfy a simple scaling relation

$$\psi_a(x,y) = \frac{1}{a}\psi_1\left(\frac{x}{a},\frac{y}{a}\right).$$ (2)

The choice of $\psi$ is dictated by the nature of the problem but most often the second derivative of the Gaussian function or the so called "Mexican hat'' is used.

The WT procedure consists of decomposing the original image into a given number of wavelet coefficient images, w_a, within the chosen set of scales a. In each wavelet image, features with characteristic sizes close to the scale are magnified and the problem is to mark the significant ones, i.e. those which are not due to noise. In most cases, this selection of significant wavelet coefficients cannot be performed analytically because of the redundancy of the WT introducing cross-correlation between pixels. For Gaussian white noise, w_a are distributed normally, allowing easy thresholding. This is not the case for X-ray images which are in the Poissonian photon noise regime.

Various techniques were developed for selecting the significant wavelet coefficients in X-ray images. In Vikhlinin et al. (1997) a local Gaussian noise was assumed; Slezak et al. (1994) used the Ascombe transformation to transform an image with Poissonian noise into an image with Gaussian noise; in Slezak et al. (1993), Starck & Pierre (1998) a histogram of the wavelet function is used. In recent years a technique based on Monte Carlo simulations is used successfully (e.g. Grebenev et al. 1995; Damiani et al. 1997; Lazzati et al. 1999).

Once the significant coefficients at each scale are chosen, the local maxima at all scales are collected and cross-identified to define objects. Different characteristics, such as centroids, light distribution etc., can be computed, as well as an indication of the source size at the scale where the object wavelet coefficient is maximal.

WT has many advantages - the multiresolution approach is well suited both for point-like and extended sources but favours circularly symmetric ones. Because of the properties of the wavelet function a smoothly varying background is automatically removed. Extensive description of wavelet transform and its different applications can be found in Starck et al. (1998).

In this work we have tested two WT procedures:

WAVDETECT:

one of the Chandra wavelet-based detection techniques (Freeman et al. 1996; Dobrzycki et al. 1999). It uses a "Mexican hat'' as wavelet function and identifies the significant coefficients by Monte Carlo simulations. The background level needed to empirically estimate the significance is taken directly from the negative annulus of the wavelet function. Source properties are computed inside the detection cell defined by minimizing the function $\vert\log r_{\rm PSF} - \log \sigma_{\rm F}\vert$, where $r_{\rm PSF}$ is the size of the PSF encircling a given fraction of the total PSF flux and $\sigma_{\rm F}$ is the size of the object at the scale closest to the PSF size. No detailed PSF shape information is needed to perform this minimization - just $r_{\rm PSF}$ as a function of the off-axis angle. The control parameters which need attention are the significance threshold and the set of scales. We have used a significance threshold of 10^-4 corresponding to 1 false event in 10000 ($\sim$$4\sigma$ in Gaussian case) and "$\sqrt{2}$ sequence'' for the scales, where $a=1,\sqrt{2},2,2\sqrt{2},4,3\sqrt{2}\dots 16$;

EWAVELET:

XMM-SAS package, based on a wavelet analysis. It is very similar to WAVDETECT but implements some new ideas. The sources are assumed to have a Gaussian shape in order to analytically derive their counts and extent. Currently, the PSF information is ignored, which can be regarded as a serious drawback. We have used the significance threshold of 10^-4 ($\sim$$4\sigma$) and wavelet scales 1, 2, 4, 8 and 16.

3.4 Mixed approach

This method combines a source detection by an elaborated procedure over filtered/smoothed raw photon images.

The use of such a mixed approach is motivated by the fact that procedures for source detection in astronomical images have been developed for many years and the steps and problems of deblending, photometry, classification of objects are now quite well understood. The raw photon image manipulations can be performed with very simple smoothing procedures (for example a Gaussian convolution) or with more sophisticated methods like the "matching filter'' technique, adaptive smoothing or multiresolution (wavelet) filtering.

We have used two different types of raw image filtering:

(1)   Gaussian convolution. For our simulated images we applied two convolutions with $FWHM=12\hbox {$^{\prime \prime }$ }$ and $20\hbox{$^{\prime\prime}$ }$ constrained by the image characteristics (see Sect. 2);

(2)   Multiresolution iterative threshold filtering (MR/1 package, Starck et al. 1998) by means of "à trous'' (with holes) wavelet method and Poissonian noise model (Slezak et al. 1993; Starck & Pierre 1998) to flag the significant wavelet coefficients (also known as auto-convolution or wavelet function histogram method). The control parameters are the significance threshold, wavelet scales and the number of iterations for the image reconstruction. We took 10^-4 significance threshold (corresponding to $\sim$$4\sigma$ for Gaussian distribution), 6 scales for the "à trous'' method allowing analysis of structures with characteristic sizes up to 2⁶ pixels ($\sim$ $4\hbox{$^\prime$ }$). We have used 25 iterations or a reconstruction error less than 10^-5. More iterations can improve the photometry but can also lead to the appearance of very faint real or sometimes false features. Optionally the last-scale wavelet smoothed image in the iterative restoration process is used for analysis of large-scale background variations or very extended features.

Source detection over the filtered images was performed with SExtractor (Bertin & Arnouts 1996) - one of the most widely used procedures for object detection in astronomy with a simple interface and very fast execution time. Originally, SExtractor was developed for optical images but it can be applied very successfully on X-ray images after they are properly filtered. The flexibility of SExtractor is ensured by a large number of control parameters requiring dedicated adjustments. In the following we will briefly outline those having a significant influence on our results:

• Detection and analysis threshold - important parameters for Gaussian convolved images. It is preferable to use a value greater than $4\sigma$, even sacrificing very faint sources but reducing the number of false detections. For MR/1 filtering this is irrelevant because the features in the images are already at $4\sigma$ significance;
• Background mesh size - influences object detection as well as photometry. Using too coarse a background grid (less grid points) will lead to a smoother background image but some faint objects will be missed in the detection. While too detailed a map (more grid points) will lead to many false detections and worse photometry, depending on the local noise properties. We have found after many experiments that local background meshes with 32 to 64 nodes give best results balancing both effects;
• Minimum detection area - in some cases this parameter can help to avoid spurious detections especially when one uses wavelet scales below the characteristic size of objects;
• Deblending parameters - the number of the deblending levels and the minimum contrast. We used the extreme values of 64 levels and zero contrast so that any saddle point will lead to object splitting;
• Photometry - we have decided to use the automatic aperture photometry procedure implemented in SExtractor following Kron (1980) and Infante (1987) which gives better results than the fixed aperture or isophotal methods;
• Classification - we discuss the classification in the next subsection.

   
3.5 Classification

The problem of classifying sources as resolved or unresolved (or stars/galaxies, extended/point-like) has a long history and discussions can be found in widespread detection packages like FOCAS (Valdes et al. 1990), INVENTORY (Kruszewski 1988), SExtractor (Bertin & Arnouts 1996), Neural Extractor (Andreon et al. 2000). This task becomes even more difficult for X-ray images - we have already mentioned various X-ray telescope effects that blur and distort the images as a function of the off-axis angle.

EMLDETECT

- the extension likelihood is calculated by the deviation of the object fitted from the PSF model when varying the size (extension) of the object. It was used with ROSAT observations to flag possible extended objects and in most of the cases the results were positive. The extension likelihood depends on the chosen PSF model and for faint objects is not reliable (De Grandi et al. 1997);

WAVDETECT

- the ratio between object size and the PSF size ( $R_{\rm PSF}$) is used for classification purposes. It depends on the adopted definition for object size and can be misleading, especially for faint objects;

SExtractor

- one parameter for classification is the stellarity index. It is based on neural-network training performed over set of some 10⁶ simulated optical images. There are 10 parameters included in the neural network training - 8 isophotal areas, the maximal intensity and the "seeing''. The only controllable parameter is the "seeing'' which, however, has no meaning for space observations. It can be tuned to the mean PSF size.
The half-light radius R₅₀ is another classification parameter. It can be used as a robust discriminator because it depends only on the photon distribution inside the object ("luminosity profile''), which is assumed to be different for resolved and unresolved objects.

Additional indications might be taken into account when the classification is ambiguous - e.g. the wavelet scale at which object's wavelet coefficient is maximal or even spectral signatures.