Документ взят из кэша поисковой машины. Адрес оригинального документа : http://xmm.vilspa.esa.es/sas/5.4.1/doc/ewavelet.ps.gz Дата изменения: Fri Jan 10 22:11:11 2003 Дата индексирования: Tue Oct 2 04:32:19 2012 Кодировка: Поисковые слова: guide 8.0

XMM-Newton Science Analysis System Page: 1
ewavelet
January 10, 2003
Abstract
The task ewavelet detects sources in an image using the mexican hat wavelet
algorithm. The algorithm detects both point sources and extended sources provided
that the source extent is not too large (< 1/8th the image size). The nal source list
will contain source positions, count rates, plus an indication of the source extent.
1 Instruments/Modes
Instrument Mode
EPIC MOS: IMAGING
EPIC PN: IMAGING
2 Use
pipeline processing yes
interactive analysis yes
3 Description
The ewavelet tasks detects sources using wavelet transforms. For this task and its description we have
made use of the Chandra Detect User Guide by Dobrzycki et al. ([1]) and the articles by Damiani et al.
([2],[3]), but the implementation will be new.
Wavelet transformations are an extension of Fourier transformations, but unlike the Fourier transfor-
mation functions (sines and cosines) wavelet functions have a well dened location in space. A wavelet
function should have a zero normalization and satisfy the property:
W a; (x) = 1
 W
 x a


; (1)
 is called the scaling or dilation parameter and a is the translation parameter. A good source about
wavelets is Holschneider ([4]). For this task we have used the two-dimensional Mexican Hat (MH) function
W (; x; y), which can be derived from the two-dimensional Gaussian function
(x; y; ) =
1
2 2
exp
 x 2 + y 2
2 2

(2)
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 2
0
0.05
0.1
0.15
­6 ­4 ­2 0 2 4 6
Z
X
Mexican Hat
Gaussian
Figure 1: Solid line: The mexican hat wavelet
function. Dashed line: a gaussian function
with the same parameter .
using the following relationship:
W (; x; y) = [(x; y)r + 2](x; y) = 1
2 2
h
2 x 2 + y 2
 2
i
exp
 x 2 + y 2
2 2

: (3)
The central part of the MH is positive and resembles a gaussian function (see Fig. 1). Outside the circle
with axis length
p
2 the MH is negative. Since the image is convolved with the MH function the MH acts
like a sort of sliding cell: the positive part being the source cell and the negative part being background
area. From this it is clear that the wavelet scale should be smaller than the image itself in order to get a
well dened image convolution. As a rule of thumb, the scale size should not exceed 1/8th of the image
size.
Sources are detected by convolving the image with the MH function for a given scale parameter .
Local maxima in the convolved image correspond to sources. The image is convolved using several scale
parameters. In this implementation, like in most implementations, the scale size is increased for each
convolution with a factor
p
2. Unlike the Chandra software we use circular symmetric wavelet functions.
The reason is that allowing  x 6=  y , increases the number of convolutions, which would make the task
considerably slower. Moreover, it would only increase the sensitivity for elongations in two very distinct
directions, ignoring source position angles in between 0 ф and 90 ф .
If we assume that the source shape has a gaussian form we can analytically derive the value of the
maximum correlation at a given wavelet scale . So the source is described by:
S(r) = N src
2a 2
exp( r 2
2a 2
) + b; (4)
where we have used the symbol a instead of  for the gaussian parameter, in order to set it apart from the
wavelet scale parameter . We have here included a constant background term, b, which is the number
of background photons per pixel. N src is the total number of photons, or another suitable normalization
(e.g. total countrate). The maximum of the convolved source occurs when the source position matches
the wavelet position, this gives for the maximum of the convolution:
Cmax = [W  S] max = 1
2 2
N src
2a 2
Z 1
0
2 r 2
 2


exp r 2
2 2


exp r 2
2a 2


2rdr =
N src
2a 2
2
1 + (=a) 2

1 1
1 + (=a) 2

(5)
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 3
Cmax peaks for a =  with Cmax = N src =4a 2 . Note that the background does not appear in eq. (5, as
the integral of the wavelet function amounts, per denition, to zero. This suggests that we can estimate
source counts and countrate by comparing C i;j for several wavelet scales and picking the wavelet scale,
max for which Cmax peaks. This gives an estimate for the source size a = max and consequently we
can estimate the source normalization with N src = Cmax 4 2
max . This aspect of the wavelet analysis
(namely source extent estimation) is neglected by the Chandra software. In section 3.1 it will be shown
how one can improve this estimate.
3.1 Countrate estimation
The wavelet detection routine is most and for all a detection routine. Countrates can, however, be
estimated in a simple way, but with less accuracy than for example a maximum likelihood tting method.
Countrates are estimated using a countrate image. This countrate image consists of the input image
divided by the exposure map. Both the input image and the countrate image are convolved with the
Mexican Hat wavelet. If a maximum is found in the convolved input image, a corresponding maximum
in the convolved countrate image will be searched for. Multiplying this maximum correlation with 4 2
gives an estimate for the source countrate.
In a separate routine the best correlation for a given source for all dierent wavelet scales will be compared
and the most appropriate wavelet scale will be used, in order to obtain the source countrate. In fact, we
can do a little better than that: correlation values at two dierent scales can be used to infer a better
estimate for the scale at which the correlation peaks. From eq. (5) it can be derived that the best
estimate for the scale,  0 , at which the wavelet correlation peaks is:
 0 =
v u u u t
1
2  2
2
q
C1
C2  2
1
q
C1
C2
1
2
; (6)
 1 ,  2 are here the two wavelet scales at which the source was detected, and C 1 and C 2 are the corre-
sponding wavelet correlation values.
Note that this algorithm is more like described by Damiani et al. ([2]), than the algorithm in the Chandra
wavelet code.
To estimate the error in the number of counts we have to calculate the variance, and for the error in
eq. (6) also the covariance, related to eq. (5). We assume here that the image to be analysed is a photon
count image, i.e. the errors per pixel follow a poisson distribution. For a gaussian source shape we obtain
for the covariance between wavelet scales  1 and  2 :
[(W 1 W 2 )  S] max =
N src
(2) 3  2
1
 2
2
a 2
Z 1
0

2 r 2
 2
1

2 r 2
 2
2

exp r 2
2 2
1
r 2
2 2
2


N src exp r 2
2a 2


+ 2a 2 b

2rdr =
N src
2 2 (a 2  2
1
+ a 2  2
2
+  2
1
 2
2
)

2 s 2
 2
1
s 2
 2
2
+ s 4
 2
1
 2
2

+ 4 2
1
 2
2
b
( 2
1
+  2
2
) 3
; (7)
where s 2 = 2a 2  2
1
 2
2
=( 2
1
 2
2
+  2
2
a 2 + a 2  2
1
).
To calculate the variance we can simply set  1 =  2 = :
[W 2
 S] max = 1
2 2
1
2 2
1
2a 2
Z 1
0

2 r 2
 2
 2
exp r 2
 2


N src exp r 2
2a 2


+ 2a 2 b

2rdr =
1
2 2
N src
2a 2
1
1 + 1
2
(=a) 2

2
2
1 + 1
2
(=a) 2
+
1
(1 + 1
2
(=a) 2 ) 2

+ b
2 2
: (8)
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 4
Figure 2: Histograms of wavelet convolved images. The wavelet scale was 2 p
2 and the background value
were respectively b = 0:01/q ' 0:5, b = 0:05/q ' 2:5, b = 0:1/q ' 5:0 and b = 1/q ' 50 (cf. Damiani et
al. [2]).
The error in the correlation is the square root of this expression.
3.2 The detection threshold
In order to determine the signicance of a detection, we have to calculate what the chance is that we nd
a peak in the wavelet convolved image in the absence of real sources. The uctuations in the wavelet
map is given by Eq. (8). For large number of (background) photons we may expect that the uctuations
in W approach a gaussian distribution, with standard deviation
q
b
2 2
:
P (Cmax > C 0 ) = 
p
b
Z 1
C0
exp
 1
2
2 2 C 2
b

dC: (9)
Because we used a dierent normalization constant in eq. (5), this equation is dierent from Damiani
et al. ([2]) and Dobrzycki et al. ([1]). However, using the transformation q = 2 2 b and C 0 = 2 2 we
obtain:
P (C 0
max > C 0
0
) = 1
p
2q
Z 1
C 0
0
exp
 1
2
C 02
q

dC 0 = 1
2
h
1 erf
 C 0
p
2q
i
: (10)
The transformation is useful, as both Dobrzycki et al. ([1]) and Damiani et al. ([2]) have made Monte
Carlo simulations to approximate the statistical behaviour in case the gaussian approximation is not valid.
They observed that this is the case for log q < 3:25. To determine signicant thresholds for the case the
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 5
Figure 3: The response of a gaussian source (dotted line) to the convolution of the image with the
bakcground smoothing function. The dashed line shows the function itself and the solid line shows the
response to a gaussian ( bgd = 2 gauss ). The gaussian has been multiplied by a factor 1/5.
gaussian approximation is not valid Damiani et al. ([2]) tted the results of Monte Carlo simulations to
the following formula:
C 0
0
(k; q) = C 0 2 2 = k p
q + (c 1 + c 2 k + c 3 k 2 ); (11)
where C 0
0
and C 0 are the detection threshold for MH wavelets with and without normalization constants
respectively, and k is the number of gaussian sigma's. The numerical constants were found to be: c 1 =
0:2336, c 2 = 0:0354 and c 3 = 0:1830. Clearly, this t has the correct assymptotic behaviour.
Note that all the above apply to correlation values of individual pixels. In reality the probabilities to
detect spurious sources above the background is a little bit smaller, as detections are dened as local
maxima above the threshold. Furthermore, since we are only interested in maxima not in minima eq 11
refers to the \single tail" probability of a gaussian distribution. As an example, in a 1024  1024 image,
1363 = 0:0013  1024 2 pixels are expected above the 3 threshold.
3.3 The convolution
From version 2.0 on ewavelet uses the fast fourier transform (FFT) to convolve images with the wavelet
lter and the background images with other smoothing lters. For small wavelet scales direct convolutions
are still used. The FFT routines are those of the external tw package. Care is taken that the input
images are only once transformed and stored, in order to avoid unnecessary calls to the FFT routines. The
Fourier Transforms of the mexican hat function and the background smoothing functions are calculated
directly using the analytical expressions for those functions (rather than using a FFT).
The direct convolution (small wavelet scales only) uses a window size of approximately 4:5. A smaller
window size decreases the computational eort at the expense of the convolution accuracy. (Note that
for a two dimensional gaussian function 3 does not correspond to a fraction of 99.73%, but to 98.9%.)
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 6
3.4 Background estimation
The wavelet detection algorithm itself does not use any separate background estimate. However, for
the determing the signicance of a candidate source detection the background value at a given pixel
coordinate has to be known.
At the moment the background is calculated by smearing the image with a function of the form
r 2 = 2 exp( 0:5r 2 = 2 ) .
This function is related to the mexican hat function through its Fourier Transform. It has the nice
property that it weighs more the outside ring, so that possible sources do not give too much of a strong
signal at the source position. For the  parameter I choose to use a value twice the size of the wavelet
scale under investigation (so for every scale a dierent background map is generated).
After one iteration the input for the background estimate is the original image minus the sources, detected
at scales smaller or equal to the wavelet scale under investigation. In this way sources with a large extent
are modelled in the background map. This prevents spurious detections of small sources within the
envelope of extended sources.
3.5 Avoiding spurious detections near edges in the image
The algorithm used to prevent the detection of spurious sources near sharp gradients in the exposure
(e.g. chip edges) has been changed as of version ewavelet v3.3 .
The background map is convolved with the current wavelet lter. This convolved image should show
similar features as the convolved input image, if no real sources were present. We only retain sources
whose value of Cmax exceeds the corresponding pixel value in the convolved background map by at least
a factor of edgethreshold.
3.6 Creating reconstruction image and background image from an input
source list
If the parameter withinputlist is set to \true", the reconstruction and/or background image will be
created using source parameters from an input source list. The lename of the input list is specied by
the parameter inputlist. Use parameters makerecon and savebkgrnd to specify which of the output
images you want to create. Only sources originally detected at scales up to maxscale will be used for to
calculate the images. The value of maxscale will also control the scale of the kernel used to lter the
background image.
3.7 Examples
1. Run ewavelet on MOS1 image, saving reconstruction image and background image. Use all permitted
convolution scales and a detection threshold of 4.0.
ewavelet imageset=M1_IMAGE_8000.FIT expmapset=M1_EXPMAP8000.FIT
srclistset=ewavlist_m1.fits minscale=1 maxscale=32 threshold=4.0
makerecon=yes recimageset=m1_rec.fits
savebkgrnd=yes bkgrndset=m1_bkg.fits
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 7
2. Use ewavelet to create a reconstruction image and background image based on the existing ewavelet
source list \my ewavlist m1.ts". Use maxscale=32 in order to include sources detected at all scales. In
this mode no source detection is performed.
ewavelet imageset=M1_IMAGE_8000.FIT expmapset=M1_EXPMAP8000.FIT
withinputlist=yes inputlist=my_ewavlist_m1.fits
makerecon=yes recimageset=my_m1_rec.fits
savebkgrnd=yes bkgrndset=my_m1_bkg.fits
maxscale=32
3.8 Implementation
This task is written in C++ and uses an image class containing the image data and the FITS parameters
associated with the world coordinates. FFT routines are taken from the FFTW package, which is fast.
4 Parameters
This section documents the parameters recognized by this task (if any).
Parameter Mand Type Default Constraints
imageset yes string image.fits
Name of the input ts image.
srclistset no string ewavlist.fits
Output data set containing the wavelet source list.
expmapset no string expmap.fits
Name of the input ts image containing the exposure map.
useexpmap no boolean true
Use an external exposure map or not (in case no proper exposure map is available).
minscale no real 2.0 1:0{16:0
Minimum wavelet scale.
maxscale no real 8.0 2:0{32:0
Maximum wavelet scale.
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 8
threshold no real 5 3{7
Detection threshold in , only sources with a signicance larger than this threshold will be stored. The
choice of the threshold will also in uence the parameters of sources above the threshold, since the back-
ground estimate depends on the number of detected sources.
edgethreshold no real 0.1 0.001{1.0
Controls the supression of spurious detection near sharp edges in the exposure mask. The exposure
mask is convolved with the MH function and only regions, where the resulting image value is less than
edgethreshold, are used for detection. Small values lead to better supression of spurious sources, but
also real sources near edges might be missed. With edgethreshold=1.0, no sources will be rejected.
niter no integer 2 2{4
Number of iterations to be used.
makerecon no boolean false
Boolean controlling whether an image with detected sources represented by gaussians will be saved as a
ts image.
recimageset no string reconstruction.ts
Name of output ts image containing the reconstructed input image. See also parameter makerecon.
savebkgrnd no boolean false
Boolean controlling whether the nal background image will be saved as a ts image.
bkgrndset no string background.ts
Name of the ts output image containing the nal background image.
withinputlist no boolean false
Boolean controlling whether reconstruction and/or background image will be derived from an input source
list (parameter inputlist) . If set to \true", no source detection will be performed.
inputlist no string ewavlist.ts
If withinputlist=\true", this le will be used as input source list.
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 9
5 Errors
This section documents warnings and errors generated by this task (if any). Note that warnings and
errors can also be generated in the SAS infrastructure libraries, in which case they would not be docu-
mented here. Refer to the index of all errors and warnings available in the HTML version of the SAS
documentation.
paramError (error)
A minimum wavelet scale was requested, larger than the maxium wavelet scale
inputImError (error)
The task encountered a serious error when trying to read the input image.
inputMapError (error)
The task encountered a serious error when trying to read the input exposure map.
reconError (error)
The task encountered problems making a model image, containing gaussians for all detected
sources. This map is used to subtract sources from the input image in order to generate a
better background image.
bgdCalc (error)
Error generating background map.
wavCalc (error)
Wavelet analysis error. Should be accompanied by warnings, which give more specic infor-
mation.
srclistError (error)
Error generating nal source list. No sources may have been found.
unsupportedDataType (warning)
corrective action: The task uses single oating point arrays to store images, doubles are
automatically changed to single oating point values, but this may cause errors.
noExpMap (warning)
In some cases an exposure map may not be available, this is, however, not recommended
and will result in this warning.
corrective action: A at exposure map is made, using the EXPOSURE keyword of the input
image
noMatch (warning)
One of the following keywords in the input image/exposure map did not match: OBS ID,EXP ID
and/or INSTRUMENT.
corrective action:
reconError (warning)
This will useally mean that the source list is empty, so that a "reconstructed image" could
not be made
corrective action: No reconstructed image will be written/or in some cases a proper back-
ground map cannot be made, which generates an error.
diCoord (warning)
World coordinates of two images (e.g. exposure map and input image) are dierent
corrective action:
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 10
noSources (warning)
No source are present in the source list.
corrective action:
coordError (warning)
Signies an error in the image world coordinates system.
corrective action: The columns RA and DEC will remain empty in the nal source list
6 Input Files
1. PPPs product (from task evselect): EPIC FITS images (one broad band image)
2. PPS product (from task eexpmap): Exposure image (for the broad band image)
7 Output Files
1. EPIC ewavelet source list. The source list will contain the following columns: SRC NUM,
IMAX, JMAX, X IMA, Y IMA, X IMA ERR, Y IMA ERR, RA,DEC, RADEC ERR, SCNTS, SCNTS ERR,
WSCALE, WCORR, XYCORR, RATE, EXTENT, EXT ERR, EXP MAP and BG MAP.
Most of them can also be found in the tasks eboxdetect and emldetect. New are:
(a) IMAX, JMAX: pixels (integer) where the maximum correlation value was found (X IMA,
Y IMA give a rened meausurement of the position)
(b) WSCALE: the wavelet scale at which the source was found. This gives an indication of
the source extent (in gaussian ). Note that if the value is comparable to the point
spread function one is probably dealing with a point source.
(c) WCORR: the maximum correlation value (C i;j in this description).
(d) EXTENT and EXT ERR: gives the source extent and its error (not correcting for the spread
function). It is dened as  0 p
8 ln 2 (see eq. 6). This corresponds to the FWHM for a
gaussian.
(e) XYCORR: Correlation value between X and Y . Large values are an indication of an
elongated source with a position angle around 45 ф .
2. FITS image containing the nal background map generated and used by ewavelet.
3. A ts image with the \reconstructed image": the detected sources are represented by gaus-
sians with  =wavelet scale and normalization the number of counts found by the algorithm.
This image will make it easy to compare results of the task with the input image.
8 Algorithm
The wavelet detection involves the following steps:
 Read the input image and the exposure map.
 Generate the countrate image from the input image and exposure map.
 Generate the initial background map by scaling the exposure map in such a way that the
background map contains as many \counts" as the input image.
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 11
next wavelet scale
Next iteration
Cross correlate sources between
wavelet scales
and update the source list
WT computation image and
countrate image
Find candidate sources
and store their properties
Save final source list
Read image and exposure map
Generate countrate image
set scale = min scale
Generate background map
(use source list if present)
Figure 4: Block diagram of the wavelet source
detection algorithm.
 Do the following for several wavelet scales:
1 generate a background map for this particular scale
2 convolve image with the MH wavelet function
3 nd local maxima in the convolved image
4 Check if the local maximum is signicant
5 If the maximum is signicant, add the point to the source list and calculate the number
of counts from the wavelet correlation value (C i;j )
 Determine if some sources were detected at several wavelets scales. If so, choose the wavelet
scale with the highest local correlation value.
 Use the detections of the same source at several wavelet scale to get some information on
the source properties.
 Make a new iteration or generate the nal source list.
xmmsas 20030110 1802-5.4.1

XMM-Newton Science Analysis System Page: 12
9 Comments
10 Future developments
 FFT routines may still need some performance optimization, as not all the functionality of
the FFTW package is currently used.
 Galactic coordinates are not yet calculated, but this can change in the near future.
 Currently preference is given to small sources over extended sources if they are present near
the same coordinates. It would be nice to have a good selection criterium to include also
extended sources in the nal source list if it is clear that the signal is not caused by the
small scale source.
 More work has to be done on the error estimates of source parameters.
 Countrates are estimated assuming that the source shape is a gaussian. Unresolved sources
have the shape of the PSF. This means that for point-like source we can in future calibrate
the true countrates versus the wavelet countrates, for small wavelet scales. However, this
involves a lot of calibration work and/or numerical simulations.
References
[1] K. Glotfelty et al. A. Dobrzycki, H. Ebeling. CHANDRA detect 1.0 User Guide. Technical Report
Version 1.0, Chandra Science Center, February 1999.
[2] G. Micela S. Sciortino F. Damiani, A. Maggio. A method based on wavelet transforms for source
detectionin photon-counting detector images. I Theory and general properties. ApJ., 483:350{369,
1997.
[3] G. Micela S. Sciortino F. Damiani, A. Maggio. A method based on wavelet transforms for source
detectionin photon-counting detector images. II Application to ROSAT PSPC images. ApJ., 483:370{
389, 1997.
[4] M. Holschneider. Wavelets: An Analysis Tool. Oxford University Press, 1995.
xmmsas 20030110 1802-5.4.1