Map Source Regions to Pixels

Getting Started

Download the sample data: 5794 (SDSS J1004+4112)

unix% download_chandra_obsid 5794 

This thread requires an input image and either a source list or a stack of region files. Below is a summary of how a user might generate such products: running celldetect on reprocessed data created by fluximage and creating regions using the roi tool.

unix% chandra_repro 5794 5794/repro
unix% cd 5794/repro
unix% fluximage acisf05794_repro_evt2.fits 5794 bin=1 psfecf=0.8
unix% celldetect 5794_broad_thresh.img expstk=5794_broad_thresh.expmap psffile=5794_broad_thresh.psfmap out=cell_src.fits
unix% roi in=cell_src.fits out='5794_%04d.reg' group=individual mode=h

Details about each command can be found in the tools help files or in the Related Links above.

Users may of course choose to run wavdetect or have their own list of regions. The working directory now looks like

unix% /bin/ls 5794*
5794_0001.reg  5794_0007.reg  5794_0013.reg  5794_0019.reg  5794_0025.reg	 5794_broad_thresh.expmap
5794_0002.reg  5794_0008.reg  5794_0014.reg  5794_0020.reg  5794_0026.reg	 5794_broad_thresh.img
5794_0003.reg  5794_0009.reg  5794_0015.reg  5794_0021.reg  5794_0027.reg
5794_0004.reg  5794_0010.reg  5794_0016.reg  5794_0022.reg  5794_0028.reg
5794_0005.reg  5794_0011.reg  5794_0017.reg  5794_0023.reg  5794_0029.reg
5794_0006.reg  5794_0012.reg  5794_0018.reg  5794_0024.reg  5794_broad_flux.img

which contains an image, 5794_broad_thresh.img, and several region files. The image with regions overlaid is shown in Figure 1

Figure 1: Image with source regions

[Version: full-size]

Figure 1: Image with source regions

Broad band (0.5-7.0keV) counts image of ObsID 5794 with celldetect source regions shown. Data are binned by one (0.492"/pixel). Only the central part of the dataset is shown.

The objective of this thread will be to assign a number, simply 1 to N (where N is the number of sources), to each of the regions and then create an image showing which pixels map to which source.

Create the Source Map

The source map will be created by filtering an image with each of the regions, individually, and assigning the pixel values inside the region to be the same as the region number. The pixel values in the input image are not important, just the dimensions and meta-data (WCS and other keywords).

This process begins by creating an image whose pixel values are all equal to 1. This is done with a dmimgcalc command.

unix% dmimgcalc 5794_broad_thresh.img none ones.fits op="imgout=1.0+img1-img1" clob+

The img1-img1 term produces an image the same size as the input but whose values are all 0. The scalar 1.0 is then added to all pixels. The result of this is an image whose pixels are all now 1

unix% dmstat ones.fits cen-
ones.fits
    min:	1 	      @:	( 3444 3190 )
    max:	1 	      @:	( 3444 3190 )
   mean:	1 
  sigma:	0 
    sum:	1453230 
   good:	1453230 
   null:	0 

but has the same image size and same meta-data (WCS and other keywords) as the original file.

NaNs

Note that if the input is a floating point/real-valued image and has values = NaN, the result of img1-img1 is NaN, not 0. This will show up as a non-zero number of null pixels in the above dmstat command.

To create the map, the process is to loop through each of the regions; multiply the image by the region number, and filter the image with the region. This will create an image whose pixels are 0 outside the region and the pixels inside the region will be the region number. There are two important things to know.

First: when an image is filtered, the default behavior is for the output image to shrink to just the bounding box around the region. This can make combining images difficult, so we will be overriding this behavior using [opt full]. (See ahelp dmopt.)

Second, when any filter is applied to a file, the DM will encode that information in the image subspace. When data are combined, these subspaces then intersect. That behavior here is unnecessary and may make the combination produce unexpected results (for example many GTI components). The subspace update can also be suppressed using [opt update=no].

The ones.fits file is now individually filtered with each region and scaled to match the region number.

unix% dmimgcalc "ones.fits[sky=region(5794_0001.reg)][opt full,update=no]" none out_1.img op=imgout="img1*1.0"
unix% dmimgcalc "ones.fits[sky=region(5794_0002.reg)][opt full,update=no]" none out_2.img op=imgout="img1*2.0"
unix% dmimgcalc "ones.fits[sky=region(5794_0003.reg)][opt full,update=no]" none out_3.img op=imgout="img1*3.0"
...
unix% dmimgcalc "ones.fits[sky=region(5794_0028.reg)][opt full,update=no]" none out_28.img op=imgout="img1*28.0"
unix% dmimgcalc "ones.fits[sky=region(5794_0029.reg)][opt full,update=no]" none out_29.img op=imgout="img1*29.0"

While shown here running one command at a time; users can easily script this step for large number of regions.

Next these individual maps need to be combined into 1 image. If it were known that the regions could not possibly overlap then these images could simply be summed together.

unix% dmimgcalc "out_*.img" none src.map op="imgout=(img1+img2+img3+...+img28+img29)"

Unfortunately this is not always the case. There frequently will be regions that overlap. The kind of map being generated can only represent one region per pixel so a choice or preference needs to be made as to which source number to use.

The easiest choice is simply to use the highest region number. This can be done using non-linear filtering tool dmimgfilt.

unix% dmimgfilt "out_*.img" src.map function=max mask="point(0,0)" clob+

This command takes the stack of per-region maps as input. The values at the same point in all of the images is extracted. The output pixel value at this location is then set to the maximum of these extracted values. The process is repeated for all the pixels in images -- basically (0,0) is replaced with each X,Y pixel pair. [This is why having all the input images retain their original size, [opt full], is important.] For most pixels there will be at most 1 pixel in the stack that is not 0 and that non-zero pixel (if any) will be represented in the output. For pixels where there are more than one non-zero input pixels, the maximum value is used.

Figure 2: Source Map

[Version: full-size]

Figure 2: Source Map

Image showing the location of which pixels belong to which source region. The dark gray pixels have the value of 0 and do not belong to any source region. Each individual source region has a unique pixel value and color coded differently.

The pixel values are color coded using the ImageJ 16_ramps contributed color look-up table

unix% ds9 src.map -cmap load $ASCDS_CONTRIB/data/16_ramps.lut 

Note

Users can pre-sort their source list before beginning this thread to optimize the max behavior. One choice might be to sort the regions based on signal-to-noise ratio in ascending order. The brightest sources will have the highest source numbers and thus if they overlap with a fainter source the higher SNR source will be identified.

unix% dmsort cell_src.fits cell_src_sort.fits snr

This process works for a reasonable number of regions and with a reasonable image size. However, for many regions and/or large images, the single dmimgfilt command may begin to consume a large amount of memory. If this is the case, then the steps above can be interleaved

unix% dmimgcalc 5794_broad_thresh.img none src.map op="imgout=img1-img1" clob+
unix% dmimgcalc src.map none ones.fits op="imgout=1.0+img1" clob+

unix% dmimgcalc "ones.fits[sky=region(5794_0001.reg)][opt full,update=no]" none out.img op=imgout="img1*1.0" clob+
unix% dmimgfilt "out.img,src.map" src.map.tmp function="max" mask="point(0,0)" clob+
unix% mv src.map.tmp src.map

unix% dmimgcalc "ones.fits[sky=region(5794_0002.reg)][opt full,update=no]" none out.img op=imgout="img1*2.0" clob+
unix% dmimgfilt "out.img,src.map" src.map.tmp function="max" mask="point(0,0)" clob+
unix% mv src.map.tmp src.map
...

The disadvantage of this method is the large amount of file I/O and the output src.map file will grow large in size as many HISTORY keywords are written to the file with each iteration. But it will work for an arbitrary number of regions and image sizes.

More complex example

The ellipses shown in Figure 1 are common and easy to visualize. However, some regions are very complex and may even contain holes in them. This same approach can be used to visualize those regions.

The Voronoi Tessellation and Percolation source detection tool, vtpdetect, is optimized to detect faint, diffuse, irregularly shaped emission. This observation has such emission that was not detected by celldetect. The regions that vtpdetect detect are complex polygons. This same technique can also be used for those regions.

First the source detection process is run. Setting the scale value helps to merge any filamentary structure with the main emission.

unix% vtpdetect 5794_broad_thresh.img 5794_broad_thresh.expmap vtp_out.fits scale=0.7 mode=h

The output from vtpdetect has the same kind of source list, SRCLIST, information as celldetect or wavdetect including elliptical shape parameters. However, it also includes the actual complex polygons it used in a second SRC_REGION extension.

unix% dmlist vtp_out.fits blocks
--------------------------------------------------------------------------------
Dataset: vtp_out.fits
--------------------------------------------------------------------------------
 
     Block Name                          Type         Dimensions
--------------------------------------------------------------------------------
Block    1: PRIMARY                        Null        
Block    2: SRCLIST                        Table        24 cols x 33       rows
Block    3: SRC_REGION                     Table         5 cols x 48       rows

Those regions, shown in Figure 3, are very complicated and actually include holes (shown as polygons with a red hash/exclude line through them).

Figure 3: VTP Source Regions

[Version: full-size]

Figure 3: VTP Source Regions

vtpdetect source regions are complex polygons. As shown they can even contain holes. These regions can be loaded directly into ds9 by specifying the extension name.

unix% ds9 5794_broad_thresh.img \
  -region "vtp_out.fits[src_region]"

The technique above can now be applied to the vtpdetect regions. The only difference here is that rather than having 1 file per region as above, this example shows how to access the individual regions in the SRC_REGION extension. As per the ASC-FITS-REGION specification the COMPONENT number in the FITS regions files is used to combine multiple shapes (in this case polygons) together. So to select all the polygons associated with a particular source, one simply needs to filter on the COMPONENT value. First, the number of regions, via the number of COMPONENT values, is determined

unix% dmstat vtp_out.fits"[src_region][cols component]" verbose=0
unix% pget dmstat out_max
33

Then the technique shown above is performed by applying a COMPONENT filter to the region extension. To reduce the number of temporary files, this thread makes use of the power of the DM on-the-fly filtering syntax, specifically that

unix% dmcopy "vtp_out.fits[src_region][component=1]" v1.reg
unix% dmimgcalc "ones.fits[sky=region(v1.reg)][opt full,update=no]" none v1.img op="imgout=img1*1"

is the exact same as

unix% dmimgcalc "ones.fits[sky=region(vtp_out.fits[src_region][component=1])][opt full,update=no]" none v1.img op="imgout=img1*1"

The region file itself can include DM filters. Using this syntax the region filtering and image scaling are applied to all 33 regions.

unix% dmimgcalc "ones.fits[sky=region(vtp_out.fits[src_region][component=1])][opt full,update=no]" none v1.img op="imgout=img1*1"
unix% dmimgcalc "ones.fits[sky=region(vtp_out.fits[src_region][component=2])][opt full,update=no]" none v2.img op="imgout=img1*2"
unix% dmimgcalc "ones.fits[sky=region(vtp_out.fits[src_region][component=3])][opt full,update=no]" none v3.img op="imgout=img1*3"
...
unix% dmimgcalc "ones.fits[sky=region(vtp_out.fits[src_region][component=33])][opt full,update=no]" none v33.img op="imgout=img1*33"

Finally, the individual per-region images are then combined just as before.

unix% dmimgfilt "v*.img" vtp.map function=max mask="point(0,0)"

The results in Figure 4 are conceptually the same as Figure 2: the pixel values identify which source region each pixel belongs to. This example shows the holes in the vtpdetect region in the extended emission more clearly than simply viewing the polygons.

Figure 4: VTP Source Map

[Version: full-size]

Figure 4: VTP Source Map

The source map for the complex vtpdetect source polygons. As shown the regions can even contain holes.

Summary

This thread showed several different techniques to create a map of which image pixels belonged to a collection of source regions. Different techniques were shown for simple and complex regions, as well as for small and large datasets.

These kinds of maps are generally useful beyond just visualization. They can be used with tools such as dmimgpick to tag individual events with their source number.

unix% dmimgpick acisf05794_repro_evt2.fits src.map mapped_evts.fits method=closest
unix% dmlist mapped_evts.fits cols
...
  17   src_map                  Real8          -Inf:+Inf             

They can be used with tools such as dmimgthresh and dmimghull to create simple convex polygons around each source.

unix% dmimgthresh vtp.map - cut=19:19 value=0 | dmimghull - 19.reg

And they can also be used with dmmaskfill to "paint" source properties onto the image

unix% dmmaskfill 'cell_src.fits[cols component,snr]' 'src.map[opt type=i4]' SNR.map

which replaces the source number with, in this example, the Signal-to-Noise ratio.

These are just some of the kinds of further processing that can be done with such maps.