Point source detection performance of Hard Xray Modulation Telescope imaging observation
Abstract
The Hard Xray Modulation Telescope (HXMT) will perform an allsky survey in hard Xray band as well as deep imaging of a series of small sky regions. We expect various compact objects to be detected in these imaging observations. Point source detection performance of HXMT imaging observation depends not only on the instrument but also on its data analysis since images are reconstructed from HXMT observed data with numeric methods. Denoising technique plays an import part in HXMT imaging data analysis pipeline alongside with demodulation and source detection. In this paper we have implemented several methods for denoising HXMT data and evaluated the point source detection performances in terms of sensitivities and location accuracies. The results show that direct demodulation with fold cross correlation should be the default reconstruction and regularization methods, although both sensitivity and location accuracy could be further imporved by selecting and tuning numerical methods in data analysis of HXMT imaging observations.
Vol.0 (200x) No.0, 000–000
1 Introduction
1.1 Background
Hard Xray Modulation Telescope (HXMT) is a planned Chinese space Xray telescope. The telescope will perform an allsky survey in hard Xray band ( – ), a series of deep imaging observations at small sky regions and pointed observations.
We expect a large number of Xray sources, e.g., AGNs, to be detected in its allsky survey. We also expect through a series of deep imaging observations at the Galactic plane Xray transients can be detected(Li, 2007; Lu, 2012). Therefore the point source detection performance is one of our concerns on HXMT data analysis.
Methods and corresponding sensitivities of pointed observation have been discussed by Jin et al. (2010). In imaging observations such as allsky survey and deep imaging at small sky regions, a variety of data analysis aspects and methods are involved. First, images are computed instead of recorded directly by the optical instrument. Mapping as well as image reconstruction methods are useful. The direct scientific data products from imaging observations of HXMT are still scientific events, more specifically, Xray photon arrival events. The attitude control system (ACS) of the spacecraft reports the attitudes periodically. We take these reported attitudes as nodes to perform certain interpolations to determine the instantaneous attitude for each scientific event. In this way, a set of parameters are assigned to each event, including the coordinates on the celestial sphere where the telescope is pointing at. Hence we call this process events mapping, where scientific events are mapped from time domain to the celestial sphere. The product of this process is refered to as the observed image, which implies the dimensionality of the data. Second, the exposure to a specific source is limited more strictly thus the signaltonoise ratio (SNR) is tightly restricted. Hence the importance of regularization methods becomes prominent. Finally, a picture is worth a thousand words. Various information can be extracted from an image by numerical methods.
In this paper we investigate the point source detection performance of the imaging and detecting system synthesised from the telescope as well as diverse combinations of data analysis methods, especially the regularization methods.
1.2 Modulation in HXMT imaging observation
The PSF of HXMT HE, which is a composite telescope consisting of 18 collimated detectors, describes the response of the telescope to a point source when the telescope is pointing at the source as well as different locations around it. In other word, the PSF is a density function of the distribution of responses occur in different observation states, which is denoted by the instantaneous attitude of the telescope as a spacecraft. So the PSF takes the attitude of the telescope as its input.
We use the proper Euler’s angles to describe the attitude of the telescope, i.e., and the longitude and latitude of the pointing, as well as denoting the rotation angle of the telescope around its own pointing axis, namely, the position angle. The modulation equation that corresponds to the imaging observation over the sphere surface is
(1) 
where is the object function (i.e., the image) defined in a compact subset of the sphere surface , is the modulation kernel function that relates the value of object function defined on a neighbourhood of the point of the subset to the instantaneous response of the telescope, while its status is .
The modulation kernel function is determined by the point spread function (PSF) of the telescope. Suppose a unit point source is located at the zenith of the celestial sphere, i.e., the point in the corresponding Cartesian coordinate system. Fix the position angle , while slew the telescope across the polar cap, and assign responses of the telescope to the unit source to pixels on the celestial sphere. Then we obtain a map represents the PSF on the celestial sphere with fixed rotation angle .
The map is then projected to a tangent plane of the celestial sphere by gnomonic mapping, i.e.,
(2) 
where and are local Cartesian coordinates on the tangent plane. Now we have
(3) 
i.e., the PSF defined on the tangent plane.
To provide for data analysis we have two discrete forms of ,
(4) 
where is the neighbourhood of the pixel of the 2D pixel grid, and the normalized form
(5) 
Given the discrete image , the detection area and the duration of exposure on each pixel , the observed data is
(6) 
where is the constant background count rate, denotes the convolution, and is the normalized PSF on the tangent plane. Eq. 6 approximates modulations in HXMT imaging observations.
Distortion occurs when projecting a set of points from (or to) a sphere surface to (or from) a plane. For example, distance between any two points, area of any continuous subset, angle between any two crossing lines (or tangents of curves) are altered nonuniformly. On the other hand, the rotation angle is not always fixed during HXMT imaging observations. Both distortions in image reconstruction from HXMT observed data and rotations during imaging observations have been discussed by Huo & Zhou (2013). Here for the sake of simplicity, we ignore them in this article.
2 Numerical methods
2.1 Single point source detection performance estimation
We estimate the single point source performance in terms of sensitivity and position accuracy through the following procedures.

Determine flux threshold for point source detection.

Simulate a frame of observed data contains only background counts.

Run denoise program on the simulated data to try to increase the signaltonoise ratio.

Demodulate the denoised data.

Run SExtractor, a source detection program, by Bertin & Arnouts (1996) on the demodulated image to detect point sources and extract their intensities, coordinates and other parameters. At this point a catalog of point sources is compiled from the simulated data. Point sources detected here, i.e., from images demodulated from background data are false detections.

Choose a cut from the histogram as the flux threshold so that a certain percent of the false detections are rejected and the rejection percentage is precise enough. The rejection percentage, e.g., or etc., reflects the significance of detections above the corresponding threshold.


Estimate detection efficiency and position accuracy of a point source of specific flux intensity.

Simulate observed data contains a single point source of given flux intensity located at in the model image.

Examine each detection is the catalog. Provided the th source in the catalog is detected at in the demodulated image, the distance between the extracted source and the true point source
(7) as well as the flux of the extracted source are investigated to determine whether the th source is a true source or not. We define the score of the current catalog in detecting the single point source as
(8) where is the index of the current catalog, and are position accuracy and flux thresholds, therefore the score reveals whether the th catalog contains the true source or not, in other word, through the previous steps (simulated observation, denoising, demodulation, source extraction and thresholding) whether we have detected the true source effectively. If we have, the outcome of these steps is counted as an effective detection of the true source, otherwise it is ineffective.

Repeat the previous steps (from 2a to 2c) times and calculate the percentage of effective detections amongst all detections, namely,
(9) which is defined as the detection efficiency.
Let be the position of the brightest source in the th catalog, the position accuracy is calculated as
(10)


Find an flux intensity so that the corresponding detection effeciency . The intensity marks the point source sensitivity of the detecting system synthesised from both the telescope in specific status and the data analysis program chain.
2.2 Imaging and mapping
2.2.1 Demodulation
Direct demodulation (DD) method (Li & Wu, 1994) is used to estimate the true image from observed data. Residual map calculated with CLEAN algorithm (Högbom, 1974) is used as lower limit constraint in DD method. The skewness of the residual map is calculated in each CLEAN iteration and its minimum absolute value serves as the main stopping criterion of iterations.
2.2.2 Crosscorrelation
In contrary to detecting and extracting point sources from demodulated images, it is also feasible to do so from crosscorrelated maps as long as the point sources are isolated with each other, compared to the FWHM of the PSF, since the position of a peak of the expected value of such a map coincides with the position of a source regardless of the PSF. Crosscorrelating the observed data and the PSF yields the correlated map
(11) 
where denotes cross correlation, and denote the observed data and the PSF respectively. The peak of that coincides with a point source of flux intensity is
(12) 
according to Eq. 6. On the other hand, the background variance of the correlated map is
(13) 
since follows Poisson distribution. Hence the significance of the peak in term of numbers of is
(14) 
where is total duration of exposure on the 2D pixel grid, is the square root of the arithme mean of over the 2D pixel grid, which is determined by the PSF as well as the range of the pixel grid only, provided the pixel grid is fine enough (see Eq. 5).
Crosscorrelation significance of isolated point source defined Eq. 14 can be evaluated directly, given only the flux intensity of the source, the background count rate, the duration of exposure, the detection area and the PSF. Hence it is determined by the object (i.e., the point source), the telescope and the status of observation thus effects from data analysis programs are minimized.
We use the crosscorrelation significance as a reference. For example, in our simulations an isolated point source of flux has significance.
2.3 Denoising
2.3.1 Linear methods
Gaussian smoothing is often used in digital image processing to suppress the noise at the cost of reduction in resolution. Tradeoff between noise suppression and resolution conservation is adjusted through the standard deviation of the Gaussian distribution serving as the smoothing kernel function. The best resolution of HXMT HE observed data is about , limited by the FWHM of its narrowfield collimator. We set to so that the FWHM of the Gaussian kernel is also .
fold cross correlation transform () can be used in DD method to regularize the illposed problem, more specifically speaking, to ensure the convergence as well as stability of the solution(Li, 2003). Here we put this technique in the denoising category. We have tested fold and fold cross correlated DD methods respectively in this article.
2.3.2 Nonlinear methods
Nonlocal means denoising (Buades et al., 2005) is an edgepreserving nonlinear denoising method. To increase its performance we have implemented this method with fast fourier transforms (FFTs). The pixelwisely evaluation of the general Euclidean distance between the th pixel and other pixels of an image is replaced by
(15) 
where is neighbourhood of the th pixel, is the pixel in the neighbourhood, is the image, is weight coefficient of the distance function. We use pixels Gaussian kernel with standard deviation (in pixels) as the weight coefficient in our simulations. We reduce the complexity of NLMeans method by computing convolutions in Eq. 15 with FFTs instead of using searching windows. (Buades et al., 2008).
Median filter is another nonlinear edgepreserving denoising method. This method is effective in removing saltandpepper noise in digital images. In HXMT observed data such noise is incurred typically by missing data or charged particles of cosmic rays. We fixed the size of the filter at (about pixels) in HXMT HE data denoising.
The last nonlinear denoising method included in this article is the adaptive wavelet thresholding with multiresolution support (Starck & Murtagh, 2006). The multiresolution support of a noisy image is a subset that contains significant coefficients only. So wavelet coefficients that dominated by noise are discarded. In this article we implemented the noniterative algorithm. The spline wavelet is used for multiresolution decomposition.
3 Simulation and results
3.1 Inorbit background simulation
HXMT HE inorbit background count rate ranges from to (Li et al., 2009). We use a constant count rate to simulate the average inorbit background of HXMT HE in this article.
3.2 Source energy spectrum and telescope detection efficiency
We use the formula proposed by Massaro et al. (2000) together with parameters fitted by Jourdain & Roques (2009)
(16) 
to model energy spectra of Crablike sources from to , where is in and the flux is in .
The detector efficiency of HXMT HE is derived from its simulated energy response, as shown in Fig. 1.
Detection efficiency of HXMT HE to a Crablike source is . Count rate of HXMT HE corresponds to the intensity is , given the detection area of HXMT HE is about .
3.3 PSF and modulation
We use the diagram shown in Fig. 2 to simulate the PSF on the tangent plane.
We use the concentric average of , namely,
(17) 
and the cumulative sum of ,
(18) 
to characterize the radial fadeout of the PSF and the concentration of the PSF respectively, as plotted in Fig. 3.
From Fig. 3 we see that the FWHM of the simulated PSF is about while responses occur in a diameter of .
Despite the fact that the direct observed data is scientific events instead of 2dimensional images, we start our simulation from simulated observed data in the form of images defined on 2dimensional cartesian pixel grid. We use a model image for simulations. Given the diameter of the PSF, the central region is fully modulated, that is, all contributions to observed data in this region are from the model image only. The surrounding region is partially modulated, i.e., part of the contributions to observed data is this region are from the model image.
The average exposure per unit solid angle is in HXMT halfyear allsky survey. The partially modulated region is discretized by pixel grid, so . The detection area of each HXMT HE detector is approximately so the total area of all the 17 HXMT HE detectors is .
3.4 Results
We have implemented several denoising methods (see Section 2.3). Fig. 4 shows denoising results by these methods.
We have simulated frames of observed data that contains the inorbit background counts only for each methods to estimate the corresponding flux thresholds by the method specified in Procedure 1 of Section 2.1. From false detections we have obtained and flux thresholds, see Table 1 for results.
Denoising  thres., in  thres., in 

w/o denoise  
1fold CCT  
2fold CCT  
Gaussian,  
NLMeans  
Median filter,  
Median filter,  
Wavelet thres. 
We have simulated frames of observd data that contains a Crablike point source of , , , , , , , , and respectively, i.e., frames of observed data in total.
With these simulated data we have estimated location accuracies as well as detection efficiencies by the method described in Procedure 2 of the following methods:

DD without denoising,

DD with fold cross correlation,

DD with fold cross correlation,

DD with Gaussian smoothing, where ,

DD with NLMeans filtering, the size of the filter is and (both parameters are in pixels),

DD with median filter,

DD with median filter, and

DD with adaptive wavelet thresholding.
Implementation details of the above methods are in Sect. 2.3.
Table 2 shows the location accuracies on simulated data of , , and point sources.
w/o denoise  

1fold CCT  
2fold CCT  
Gaussian,  
NLMeans  
Median filter,  
Median filter,  
Wavelet thres. 
Table 3 shows the detection efficiencies on simulated data of , , , and point sources.
Although the RL iteration itself we employed in DD method is totalcountsconservative, i.e., the sum of counts of each pixel is conserved after the iteration (Richardson, 1972), the regularizations including the background constrains as well as various denoising techniques are not necessarily countsconservative. As a result, the absolute flux threshold (Table 1) for rejecting false detections doesn’t reflect the sensitivity directly but the detection efficiency (Table 3) does.
Errors of flux thresholds, location accuracies and detection efficiencies in Table 1, Table 2 and Table 3 are calculated by bootstrapping. For example, following Procedure 1 a set of frames of demodulated images are obtained, from which false detections are calculated and a histogram is plotted, where both the and thresholds are determined. Now let’s generate a new set with the same volume by resampling from the original set with replacement in order to calculate both the thresholds again. We repeat this resampling process until we get enough thresholds calculated to estimate their standard deviations. In Table 1, Table 2 and Table 3 each of the errors is a standard deviation that calculated from resampled sets.
w/o denoise  

1fold CCT  
2fold CCT  
Gaussian,  
NLMeans  
Median filter,  
Median filter,  
Wavelet thres. 
A comprehensive summary of all the tested methods on all simulated data is shown in Fig. 5.
4 Conclusion
According to the results from our tests no denoising method shows significant advantage over fold cross correlated DD method in single point source detection efficiency. Therefore it’s suggested that fold cross correlation should be the default regularization method for single point source detection in HXMT imaging data analysis.
In the other hand, the location accuracy can be improved with alternative denoising methods, such as median filter, wavelet thresholding, Gaussian smoothing with little kernel, or without denoising, according to the results in this work.
This article is focused on the single point source detection of HXMT imaging data analysis, where other interesting topics can not be all covered. Although the alternative denoising methods have been outperformed more or less by the default fold cross correlation in their contributions to the detection efficiency, they have shown certain advantages in location accuracy, and these features are promising for locating bright transients, multiple sources resolving, and so on.
Acknowledgements
In this work we made use of SciPy(Jones et al., 2001–), PyFITS and AIRE. PyFITS is a product of the Space Telescope Science Institute, which is operated by AURA for NASA. AIRE is a set of computing facilities initiated by Tsinghua Centre for Astrophysics. This work was supported by National Natrual Science Foundation of China (NSFC) under grants No. 11373025, No. 11173038 and No. 11403014 as well as Tsinghua University Initiative Scientific Research Program under grant No. 20111081102.
References
 Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, Astronomy and Astrophysics Supplement Series, 117, 393
 Buades et al. (2005) Buades, A., Coll, B., & Morel, J.M. 2005, in Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, vol. 2, 60–65 (IEEE)
 Buades et al. (2008) Buades, A., Coll, B., & Morel, J.M. 2008, International Journal of Computer Vision, 76, 123
 Högbom (1974) Högbom, J. A. 1974, \aaps, 15, 417
 Huo & Zhou (2013) Huo, Z.X., & Zhou, J.F. 2013, Research in Astronomy and Astrophysics, 13, 991
 Jin et al. (2010) Jin, J., Chen, Y., Zhang, S.N., et al. 2010, Chinese Physics C, 34, 66
 Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python
 Jourdain & Roques (2009) Jourdain, E., & Roques, J. 2009, The Astrophysical Journal, 704, 17
 Li et al. (2009) Li, G., Wu, M., Zhang, S., & Jin, Y.K. 2009, Chinese Astronomy and Astrophysics, 33, 333
 Li (2003) Li, T.P. 2003, in High Energy Processes and Phenomena in Astrophysics, IAU Symposium, vol. 214, edited by X. D. Li, V. Trimble, & Z. R. Wang, 70
 Li (2007) Li, T.P. 2007, Nuclear Physics B  Proceedings Supplements, 166, 131, proceedings of the Third International Conference on Particle and Fundamental Physics in Space
 Li & Wu (1994) Li, T.P., & Wu, M. 1994, Astrophysics and Space Science, 215, 213
 Lu (2012) Lu, F. 2012, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, vol. 8443
 Massaro et al. (2000) Massaro, E., Cusumano, G., Litterio, M., & Mineo, T. 2000, \aap, 361, 695
 Richardson (1972) Richardson, W. H. 1972, JOSA, 62, 55
 Starck & Murtagh (2006) Starck, J.L., & Murtagh, F. 2006, Astronomical Image and Data Analysis (Berlin: SpringerVerlag)