Calculating Statistics of Images

[CIAO 3.4 Science Threads]

Overview

Last Update: 1 Dec 2006 - reviewed for CIAO 3.4: ChIPS version

Synopsis:

Although in many cases it is easier (and faster) to use the dmstat tool, the sstats.sl script described in this thread allows you to ignore regions of the image by using a mask file.

Purpose:

Use the S-Lang interpreter in ChIPS and Sherpa to calculate statistics of images.

Read this thread if:

you are interested in calculating the statistics of an ACIS or HRC image.

Related Links:

Estimating Source Counts thread: related functionality of the dmextract tool.
dmopt and dmstat ahelp pages: an alternative method of marking image regions as "uninteresting" - setting such pixels to the Null value for that image type - is discussed. Note that the sstats.sl script will not handle such cases.

Proceed to the HTML or hardcopy (PDF: A4 | letter) version of the thread.

Getting Started
Create an image of the region of interest
Statistics of the selected region
1. dmstat
2. sstats.sl
Further examples
History
Images
- Annulus in which we wish to measure the signal level
- Image mosaic

Getting Started

Sample ObsID used: 1838 (ACIS-S, G21.5-09)

It is assumed that you have completed the Single Chip Exposure Map thread.

This thread uses the sstats.sl script. The most recent version of sstats.sl is v0.4 (4 Oct 2001):

unix% grep Id $ASCDS_CONTRIB/share/slsh/local-packages/sstats.sl
% $Id: sstats.sl,v 0.4 2001/10/04 18:45:12 dburke Exp $

Note that $ASCDS_CONTRIB/share/slsh/local-packages/ is the default path in the standard CIAO scripts installation; see the Scripts page for more information. Please check that you are using the most recent version before continuing. If you do not have the script installed or need to update to a newer version, please refer to the Scripts page.

Create an image of the region of interest

For this thread, we want to look at the average signal level in an annulus surrounding the source center (this is for the 1838_img_s3.fits file created in the Single Chip Exposure Map thread). We saved the selected region to stat.reg in CIAO format; its contents are:

unix% cat stat.reg
# Region file format: CIAO version 1.0
annulus(4070.5,4248.5,130,218)

Further details are available for interactive region definition using ds9.

We create an image of the region using this spatial filter. For the mask region - which should contain pixel values greater than 0 only for those pixels within the filter - we use a filtered copy of the exposure map (we can not use img_stat.fits as the mask image, since some of the pixels within the annulus contain no counts).

unix% dmcopy "1838_img_s3.fits[sky=region(stat.reg)]" img_stat.fits
unix% dmcopy "expmap_1.7kev_7.fits[sky=region(stat.reg)]" emap_stat.fits

We show a mosaic of the original image (log scaling) together with the filtered image (sqrt scaling) and the filtered exposure map (linear scaling).

Statistics of the selected region

We show results from using the dmstat tool and the sstats.sl S-Lang script:

A. dmstat

unix% punlearn dmstat
unix% dmstat img_stat.fits centroid=no
EVENTS_IMAGE
    min:        0             @:        ( 4057.6402575 4030.9302575 )
    max:        2             @:        ( 4089.6402575 4082.9302575 )
   mean:        0.027786151623 
  sigma:        0.16761614994 
    sum:        2673 
   good:        96199 
   null:        94770

The mean and standard deviation values represent the signal level within the annulus. This could equivalently be found from the original image:

unix% dmstat "1838_img_s3.fits[sky=region(stat.reg)]" centroid=no
EVENTS_IMAGE
    min:        0             @:        ( 4057.6402575 4030.9302575 )
    max:        2             @:        ( 4089.6402575 4082.9302575 )
   mean:        0.027786151623 
  sigma:        0.16761614994 
    sum:        2673 
   good:        96199 
   null:        94770

B. sstats.sl

unix% chips --slscript sstats.sl

Welcome to ChIPS, version CIAO 3.4
Copyright (C) 1999-2003, Smithsonian Astrophysical Observatory

chips> sstats( "img_stat.fits" ) 
Statistics of img_stat.fits.
Mean =  1.3997e-02  RMS =  1.1977e-01  Total      =  2.6730e+03
Min  =  0.0000e+00  Max =  2.0000e+00  Num pixels = 190969

The statistics have been calculated over the whole image. To restrict the calculation to just those pixels within the annulus we need a mask: an image the same size as img_stat.fits in which those pixels which should be considered for the calculation have values > 0. In this case we use the filtered exposure map (emap_stat.fits):

chips> sstats( "img_stat.fits", "emap_stat.fits" )
Statistics of img_stat.fits using (emap_stat.fits > 0) as a mask.
Mean =  2.7784e-02  RMS =  1.6761e-01  Total      =  2.6730e+03
Min  =  0.0000e+00  Max =  2.0000e+00  Num pixels = 96205

It is also possible to pass the routine "virtual files", so it is not necessary to create the images beforehand:

chips> sstats( "1838_img_s3.fits[sky=region(stat.reg)]", "expmap_1.7kev_7.fits[sky=region(stat.reg)]" )
Statistics of 1838_img_s3.fits[sky=region(stat.reg)] using (expmap_1.7kev_7.fits[sky=region(stat.reg)] > 0) as a mask.
Mean =  2.7784e-02  RMS =  1.6761e-01  Total      =  2.6730e+03
Min  =  0.0000e+00  Max =  2.0000e+00  Num pixels = 96205

chips> exit

The reason for the small difference in the RMS values reported by dmstat and sstats() is because dmstat calculates

  rms = sqrt( (mean(x*x) - mean(x)*mean(x)) / N )

whereas sstats() uses

  rms = sqrt( (mean(x*x) - mean(x)*mean(x)) / (N-1) )

where x is a list of N numbers and mean(x) refers to the mean value of the list x.

Further examples

Above we show an example of looking at the signal level within the region in different energy bands, using the virtual file syntax to filter and bin the event list "on the fly". In this example, we do not apply any spatial filtering to the event list; instead we create a filtered, full-size image of the exposure map using the DM [opt full] syntax.

unix% dmcopy "expmap_1.7kev_7.fits[sky=region(stat.reg)][opt full]" emap_stat_full.fits

unix% more chip_s3.reg 
# Region file format: CIAO version 1.0
rotbox(4182.860000,4393.150000,1024,1024,93.192576)

unix% chips

Welcome to ChIPS, version CIAO 3.4
Copyright (C) 1999-2003, Smithsonian Astrophysical Observatory

chips> () = evalfile( "sstats.sl" )

chips> sstats( "acisf01838N001_evt2.fits[energy=500:2000,sky=region(chip_s3.reg)][bin sky=1]", "emap_stat_full.fits" )
Statistics of acisf01838N001_evt2.fits[energy=500:2000,sky=region(chip_s3.reg)][bin sky=1] using (emap_stat_full.fits > 0)
 as a mask.
Mean =  1.0561e-02  RMS =  1.0273e-01  Total      =  1.0160e+03
Min  =  0.0000e+00  Max =  2.0000e+00  Num pixels = 96205

chips> sstats( "acisf01838N001_evt2.fits[energy=2000:8000,sky=region(chip_s3.reg)][bin sky=1]", "emap_stat_full.fits" )
Statistics of acisf01838N001_evt2.fits[energy=2000:8000,sky=region(chip_s3.reg)][bin sky=1] using (emap_stat_full.fits > 0)
 as a mask.
Mean =  1.4583e-02  RMS =  1.2100e-01  Total      =  1.4030e+03
Min  =  0.0000e+00  Max =  2.0000e+00  Num pixels = 96205

History

14 Dec 2004	updated for CIAO 3.2: path for script
20 Dec 2005	reviewed for CIAO 3.3: no changes
01 Dec 2006	reviewed for CIAO 3.4: ChIPS version