Map Source Regions to Pixels
CIAO 4.14 Science Threads
- Running celldetect
- Using Contrib Color Look-Up Tables
- Single Chip ACIS Exposure Map and Exposure-corrected Image
- ASC-FITS-REGION specification, Rots and McDowell.
Last Update: 15 Feb 2022 - Review for CIAO 4.14. Updated for Repro-5.
- Getting Started
- Create the Source Map
- More complex example
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
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.
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
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
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
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
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.
|30 Jan 2014||Initial version.|
|07 Feb 2014||Minor edits.|
|23 Dec 2014||Review for CIAO 4.7; no changes.|
|07 Dec 2018||Updated for CIAO 4.11 celldetect and fluximage changes.|
|15 Feb 2022||Review for CIAO 4.14. Updated for Repro-5.|