Last modified: 3 Apr 2019

URL: https://cxc.cfa.harvard.edu/ciao/threads/masks/

Using Pixel Mask Filters

CIAO 4.12 Science Threads


Overview

Synopsis:

A pixel mask is a 2D image whose pixel values are used to filter a dataset (both images and tables). Masks can be used as an alternative to classic regions shapes (circle, box, polygon, etc) or can be used in conjunction with them to create complex filters.

This thread will provide an introduction to pixel masks in the context of the task of creating a radial profile for extended emission while excluding overlapping point sources.

This technique mimics steps taken in the Obtain and Fit a Radial Profile thread.

Purpose:

Users should read this thread if they want to learn how to use pixel masks to filter their data.

Related Links:

Last Update: 3 Apr 2019 - Updated to use matplotlib for plotting.


Contents


Getting Started

Download the sample data: 16136 (Abell 2626)

It is assumed that users have recalibrated their dataset using chandra_repro before continuing.

Pixel masks can be used with both tables (event files) and images. Conceptually working with images is a little bit easier so we begin by creating a broad-band (0.5-7.0keV) image of the Abell 2626 cluster in OBS_ID 16136 along with the corresponding exposure map.

unix% fluximage acisf16136_repro_evt2.fits out=abell2626 bin=1 clob+ psfecf=0.9
Running fluximage
Version: 28 November 2018

Using CSC ACIS broad science energy band.
Aspect solution pcadf498642854N002_asol1.fits found.
Bad-pixel file acisf16136_repro_bpix1.fits found.
Mask file acisf16136_001N002_msk1.fits found.

The output images will have 2092 by 4308 pixels, pixel size of 0.492 arcsec,
    and cover x=3115.5:5207.5:1,y=1762.5:6070.5:1.

Running tasks in parallel with 4 processors.
Creating aspect histograms for obsid 16136
Creating 4 instrument maps for obsid 16136
Creating 4 exposure maps for obsid 16136
Combining 4 exposure maps for obsid 16136
Thresholding data for obsid 16136
Exposure-correcting image for obsid 16136
Creating PSF map for obsid 16136

The following files were created:

 The clipped counts image is:
     abell2626_broad_thresh.img

 The clipped exposure map is:
     abell2626_broad_thresh.expmap

 The PSF map is:
     abell2626_broad_thresh.psfmap

 The exposure-corrected image is:
     abell2626_broad_flux.img

The output from fluximage is shown in Figure 1

Figure 1: Broad band image of Abell 2626

[Thumbnail image: broad band image of Abell 2626 cluster shown with radial profile annuli]

[Version: full-size]

[Print media version: broad band image of Abell 2626 cluster shown with radial profile annuli]

Figure 1: Broad band image of Abell 2626

(Left): The broad band (0.5-7.0keV) image of Abell 2626. (Right): The same image with a series of annuli centered at (x,y)=(4077.3,4269.7) with radii going from 0 to 600 pixels in 10-pixel increments.

As can be seen there are several point-like sources superimposed on the the extended emission. Those sources need to be excluded from the dataset. Also, the outer annuli extend beyond of the edge of the detector so the field-of-view boundary will need to be included to compute the area correctly.

The point sources can be detected with wavdetect and the output source list can be used to filter the dataset.

unix% punlearn wavdetect wrecon wtransform
unix% pset wavdetect \
  infile=abell2626_broad_thresh.img \
  expfile=abell2626_broad_thresh.expmap \
  psffile=abell2626_broad_thresh.psfmap \
  outfile=abell2626.srcs \
  scellfile=abell2626.cells \
  imagefile=abell2626.recon \
  defnbkgfile=abell2626.nbkg \
  scales="1.4 2.0 4.0 8.0 12.0 16.0" \
  interdir=`pwd` clobber=yes

unix% wavdetect mode=h

The choice of wavelet scales was selected to only identify the point-like objects while leaving the structure in the extended emission unresolved. The output from wavdetect is shown in Figure 2.

[IMPORTANT]
Important

Users should always carefully scrutinize the output from any source detection algorithm.

Figure 2: Broad band image of Abell 2626 with wavdetect sources

[Thumbnail image: Broad band image of Abell 2626 with wavdetect sources]

[Version: full-size]

[Print media version: Broad band image of Abell 2626 with wavdetect sources]

Figure 2: Broad band image of Abell 2626 with wavdetect sources

The broad band image of Abell 2626 with the wavdetect detected sources shown.

The wavdetect source list is sufficient for the purposes of this thread.

Now that the source have been identified the objective is to exclude those source from the dataset when creating a radial profile using the annuli shown in Figure 1 (right).



Traditional approach: using CIAO regions

We start by computing the radial profile using the traditional CIAO region files so that we can use it for comparison later on. CIAO filtering supports excluding regions explicitly. This technique is used in the An Image of Diffuse Emission thread as well as the Obtain and Fit a Radial Profile thread.

To be more efficient, the we follow the advice in this FAQ and will manually invert the region before using it.

unix% python -c 'from region import *;src=CXCRegion("abell2626.srcs");(field()-src).write("excluded_srcs.reg",fits=True);'
unix% dmlist excluded_srcs.reg"[#row=1:5]" data,clean
Region Block: Field()&!Ellipse(3976.07,3767.58,10.5009,6.77157,131.978)&!Ellipse(4428.43,3822.05,13.0252,8.25924,73.7475)&!Ellipse(3815.19,3930.45,9.5287,6.53082,99.2819)&!Ellipse(3757.5,3960.5,11.1467,6.0769,61.0829)
#  POS(X,Y)                                 SHAPE              R[2]                                     ROTANG[2]                                COMPONENT
                                  NaN NaN Field                                               NaN NaN                                  NaN NaN          1
      3976.0663265306      3767.5765306122 !Ellipse                  10.5008506775         6.7715659142                       131.9777679443 NaN          1
      4428.4257425743      3822.0495049505 !Ellipse                  13.0252161026         8.2592449188                        73.7475128174 NaN          1
      3815.1914893617      3930.4468085106 !Ellipse                   9.5287017822         6.5308151245                        99.2818756104 NaN          1
      3757.4956521739      3960.5043478261 !Ellipse                  11.1467132568         6.0769000053                        61.0829200745 NaN          1

The source list has been converted into a CIAO region file where each source is now explicitly excluded from the entire field().

The broad band image is then filtered with this region using dmcopy

unix% dmcopy "abell2626_broad_thresh.img[sky=region(excluded_srcs.reg)]" image_point_srcs_removed_regions clob+

The resulting image is shown in Figure 3 with the source regions drawn for reference. The point sources have now been removed from the dataset.

Figure 3: Broad band image with point source excluded

[Thumbnail image: Broad band image with point source excluded]

[Version: full-size]

[Print media version: Broad band image with point source excluded]

Figure 3: Broad band image with point source excluded

The broad band image after the point sources have been excluded.

CIAO retains the information about all filters that have been applied to a dataset in the file's subspace. This can be displayed using dmlist

unix% dmlist image_point_srcs_removed_regions subspace

--------------------------------------------------------------------------------
Data subspace for block EVENTS_IMAGE: Components: 3 Descriptors: 17 
--------------------------------------------------------------------------------
 
 --- Component 1 --- 
   1 sky                  Real8               Rectangle(3115.5,1762.5,5207.5,6070.5)
&!Ellipse(3976.07,3767.58,10.5009,6.77157,131.978)&!Ellipse(4428.43,3822.05,13.0252,8.25924,73.7475)
&!Ellipse(3815.19,3930.45,9.5287,6.53082,99.2819)&!Ellipse(3757.5,3960.5,11.1467,6.0769,61.0829)
&!Ellipse(3874.58,3965.54,9.35435,6.95154,2.23643)&!Ellipse(3607.18,4183.19,9.9418,7.06736,73.3857)
&!Ellipse(3846.73,4220.56,4.81187,3.51667,89.5624)&!Ellipse(3992.24,4241.9,3.70426,3.54155,10.6357)
&!Ellipse(4071.33,4263.32,5.58344,5.06664,32.5627)&!Ellipse(4239.77,4264.7,3.91227,3.60294,131.364)
&!Ellipse(4490,4279.52,9.35693,3.81444,2.35721)&!Ellipse(4226.91,4295.07,11.3243,7.12405,97.8399)
&!Ellipse(3827.62,4358.72,15.8993,10.2424,175.903)&!Ellipse(4293.27,4394.1,10.3766,7.30043,88.0386)
&!Ellipse(4233.58,4444.42,5.35299,4.67135,106.273)&!Ellipse(3635.71,4580.23,13.7746,8.69293,61.3198)
&!Ellipse(4331.15,4605.21,18.3689,10.2362,28.6763)&!Ellipse(4135.83,4737.6,20.031,14.7404,32.8505)
&!Ellipse(4333.77,4837.71,15.2453,10.6939,30.6961)...[truncated]
   1 sky                  Real8               Field area = 9.01234e+06 Region area = 8.96331e+06

   1 sky                  [ 1] x                   3115.50:     5207.50 
   1 sky                  [ 2] y                   1762.50:     6070.50 
   2 expno                Int4                3:35762 
   3 ccd_id               Int2                5:5,7:7 
   4 node_id              Int2                0:3 
   5 chipx                Int2                1:1024 
   6 chipy                Int2                1:1024 
   7 tdetx                Int2                1:8192 
   8 tdety                Int2                1:8192 
   9 detx                 Real4                   0.50:     8192.50 
  10 dety                 Real4                   0.50:     8192.50 
  11 phas                 Int2                -4096:4095 
  12 pha                  Int4                0:36855 
  13 pha_ro               Int4                0:36855 
  14 energy               Real4                     500.0:     7000.0 
  15 pi                   Int4                1:1024 
  16 fltgrade             Int2                0:255 
  17 grade                Int2                0:0,2:2,3:3,4:4,6:6 
 --- Component 2 --- 
 ...

Here we can see each of the Ellipses that have been excluded from the dataset. The information in the subspace is used by dmextract when creating the radial profile.

[NOTE]
Note

The truncated message only means that dmlist has truncated the value being displayed; all of the excluded ellipses are stored in the subspace and are used by CIAO tools.

As is highlighted, the sky subspace starts with Rectangle(3115.5,1762.5,5207.5,6070.5). This represents the entire image, not just the part of the image that has data (ie has non-zero exposure). That is, the boundary of the detector is not included in the dataset. The Field Of View (fov) boundary will need to be added when computing the radial profile since the outer annuli extend past the edge of the detect. This will be done with the fov file.

The radial profile is then extracted with dmextract

unix% dmextract image_point_srcs_removed_regions"[sky=region(acisf16136_repro_fov1.fits)][bin sky=annulus(4077.3,4269.7,0:600:10)]" \
  exp="abell2626_broad_thresh.expmap" \
  out=radial_profile_with_regions_from_image mode=h op=generic clob+
[TIP]
Background

For purposes of this thread background has been ignored. Users would typically need to include background when working with regions that cover large areas.

Figure 4: Radial profile from broad band image with excluded sources using regions

[Thumbnail image: Radial Profile from broad band image with excluded sources using regions]

[Version: full-size]

[Print media version: Radial Profile from broad band image with excluded sources using regions]

Figure 4: Radial profile from broad band image with excluded sources using regions

Radial profile of Abell 2626 broad-band (0.5-7.0keV) image. Contaminating point sources have been removed using CIAO regions and the edge of the detector has been taken into account by using the FOV file.

The plot is created in python using the virtual columns CEL_R and CEL_FLUX.

from pycrates import *
import matplotlib.pylab as plt

tab=read_file("radial_profile_with_regions_from_image")
x = tab.get_column("cel_rmid").values
y = tab.get_column("cel_flux").values
e = tab.get_column("cel_flux_err").values

plt.errorbar(x,y,yerr=e, drawstyle="steps-mid", ecolor="gray")
plt.yscale("log")
plt.xlabel("CEL_R (arcsec)")
plt.ylabel("CEL_FLUX (photon/cm**2/arcsec**2/s)")
plt.title("Abell 2626")

The output from this command is shown in Figure 4. This command takes approximately 60 seconds to complete; most of the time is spent dealing with the large number of ellipses in the region subspace intersecting with the field-of-view polygons, and with each of the annuli.

This is the final product. Later sections will compare their results to this output.



Alternative approach: using pixel masks

In this section the radial profile will be computed in the same way as before, except pixel masks will be used to exclude the point sources and to account for the edge of the detector.

What is a pixel mask?

A pixel mask is a two dimensional image whose pixel values are used to filter a dataset (images or tables). An image pixel value equal to zero indicates that the data are to be filtered out (removed, excluded); any finite, non-zero value indicates that the data are allowed to pass the filter (remain included).

Fast facts:

  • A pixel mask may have a world coordinate system attached to it. If so, it must match the WCS of the dataset being filtered.

  • Pixel masks may be any image data type (integer or real valued).

  • Special image value including IEEE NaN, +/-Inf, and integer BLANK values are all treated as "0".

The help file for dmmasks contains all the details about properties of pixel masks.


How to create a pixel mask

Since a pixel mask is just a 2D image, any CIAO tool or applications that outputs a 2D image can be used.

One way to create is mask to start with a unit image (all pixel values equal to 1). Then use regions to filter the image into the desired mask.

When working with images, the world-coordinate system in the mask and the dataset being filtered must match, so it is useful to begin with the image itself and use it as the starting point for creating the mask.

First, dmimgcalc is used to create the unit image (all pixel values equal to 1)

unix% dmimgcalc abell2626_broad_thresh.img none ones.fits op="imgout=(1+(img1-img1))" clob+

dmgimgcalc was used with a single input image. The operation, op, subtracts the image from itself so that all pixel values are zero, and then adds 1 to all the pixels. We can check that all pixels are equal to 1 using dmstat

unix% dmstat ones.fits sig- cen- med-
ones.fits
    min:	1 	      @:	( 3116 1763 )
    max:	1 	      @:	( 3116 1763 )
   mean:	1 
    sum:	9012336 
   good:	9012336 
   null:	0 

Since min equals max, and since there are no null values, we can be sure that all the image pixel values are equal to 1.

Next, the unit image is filtered with the original wavdetect source list:

unix% dmcopy "ones.fits[sky=region(abell2626.srcs)][opt full]" srcs.mask clob+

The output is shown in Figure 5. The opt full option is used to retain the original image size.

Figure 5: Unit mask filtered with source list

[Thumbnail image: Image showing sources]

[Version: full-size]

[Print media version: Image showing sources]

Figure 5: Unit mask filtered with source list

The unit image has been filtered with the wavdetect source list to produce a map showing where the point sources are located. Pixel values equal to 1 indicate source pixels, pixels values equal to 0 are outside the source regions.

The mask now has pixel values equal to 1 in all pixels that are included in any source region. However, what is needed in this exercise is to exclude those pixel in the source region. Users can use an [exclude ] filter with pixel masks, but it is also possible to invert the mask using generic datamodel image manipulation tools. For example

unix% dmimgcalc srcs.mask none no_srcs.mask op="imgout=(1-img1)" clob+

The mask has been inverted by subtracting it from "1". The inverted mask is shown in Figure 6.

Figure 6: Inverted unit mask filtered with source list

[Thumbnail image: Image showing sources]

[Version: full-size]

[Print media version: Image showing sources]

Figure 6: Inverted unit mask filtered with source list

The pixel mask shown in Figure 5 has been inverted using dmimgcalc.

Now, all the pixels equal to 1 (white) will be included in the radial profile, whereas the pixels equal to 0 (black) will be excluded.

Looking at this image it is clear that the edge of the detector has not yet be taken into account (the field of pixel values equal to 1 extends to the left and right edges of the image, but based on Figure 1 we know the detector does not extend that far.

The edge of the detector can be folded into the mask in several ways. One way is to use the exposure map. We can use dmimgthresh apply a threshold to the exposure map and use that to filter the inverted source mask.

unix% dmimgthresh no_srcs.mask out=no_srcs_with_fov.mask exp=abell2626_broad_thresh.expmap cut=1 value=0 clob+

Here dmimgthresh applies a threshold to the exposure map file, abell2626_broad_thresh.expmap. The input image is copied to the output image. Then pixels in output that correspond to pixels in the exposure map which are less than 1.0 are set to 0. The output is shown in Figure 7.

[TIP]
Tip

In this example cut=1 was used because typical exposure map pixel vales are ~100 cm2*ONTIME, thus typically in the ~1.0e6 range. Since the threshold is applied to values strictly below the cut value, a cut=0 would not produce the desired result. Therefore a small number, much less than the typical values is used.

Figure 7: Pixel mask with excluded sources and detector edges

[Thumbnail image: Pixel mask with excluded sources and detector edges]

[Version: full-size]

[Print media version: Pixel mask with excluded sources and detector edges]

Figure 7: Pixel mask with excluded sources and detector edges

A threshold was applied to the exposure map to filter the mask file. The edges of the detectors are now properly accounted for.

The mask shown in Figure 7 is the equivalent mask to the source regions and field-of-view filters used above.


Visualize mask

There are several ways to visualize a mask on top of another dataset. Figure 8 shows how to use ds9, Figure 9 shows how to use matplotlib, and Figure 10 shows how to use dmimg2jpg

Figure 8: Visualize mask with ds9

[Thumbnail image: using ds9 to view mask]

[Version: full-size]

[Print media version: using ds9 to view mask]

Figure 8: Visualize mask with ds9

The ds9 -mask option can be used to display a mask on top of another image. Both images must be the same size with the same WCS.

unix% ds9 abell2626_broad_thresh.img -scale log \
  -block 4 -scale log -pan to 4077.3 4269.7 physical \
  -mask color green -mask transparency 60 -mask no_srcs_with_fov.mask

The -mask option does not work very well when displaying event files.

Figure 9: Visualize mask with matplotlib

[Thumbnail image: using python to view mask]

[Version: full-size]

[Print media version: using python to view mask]

Figure 9: Visualize mask with matplotlib

The mask can be displayed by using overlapping images with matplotlib with adjustable transparency, alpha.

import matplotlib.pylab as plt
from pycrates import read_file
from crates_contrib.utils import scale_image_crate

cr = read_file("abell2626_broad_thresh.img")
scale_image_crate(cr,"asinh")
pixels = cr.get_image().values
msk = read_file("no_srcs_with_fov.mask").get_image().values

plt.imshow(pixels, cmap="gray", origin="lower", aspect="equal")
plt.imshow(msk, cmap="Greens", alpha=0.4, origin="lower", aspect="equal")
plt.xlim(400,1400)
plt.ylim(2000,3000)

Figure 10: Visualizing mask with dmimg2jpg

[Thumbnail image: image with mask made with dmimg2jpg]

[Version: full-size]

[Print media version: image with mask made with dmimg2jpg]

Figure 10: Visualizing mask with dmimg2jpg

dmimg2jpg can also be used to visualize the mask. As with other regions, it draws the outline around the mask rather than shading in the entire masked region.

unix% dmimg2jpg infile="abell2626_broad_thresh.img[bin sky=4]" outfile=dmimg2jpg_40.jpg \
  regionfile="mask(no_srcs_with_fov.mask[bin sky=4][opt full])" regioncolor=")colors.green" mode=h clob+

The [opt full] is needed here to retain the original mask image size. Without it, the mask would shrink to the bounding box around the excluded sources.


Filtering with pixel mask file

Now that the mask has been created it will be used to filter the broad band image.

The syntax is similar to the region filter, but using the mask() token keyword

unix% dmcopy "abell2626_broad_thresh.img[sky=mask(no_srcs_with_fov.mask)][opt full]" image_with_mask clob+

Again, the opt full option is used to retain the original image size. The output image looks the same as what is shown in Figure 3.

Similar to that file, the subspace can be displayed:

unix% dmlist image_with_mask subspace

--------------------------------------------------------------------------------
Data subspace for block EVENTS_IMAGE: Components: 3 Descriptors: 17 
--------------------------------------------------------------------------------
 
 --- Component 1 --- 
   1 sky                  Real8               TABLE MASK
                                              MASK(MASK)
                                              Field area = 9.01234e+06 Region area = 4.35209e+06

                          [ 1] x              
                                                   3115.50:     5207.50
                          [ 2] y              
                                                   1762.50:     6070.50
   2 expno                Int4                3:35762 
   3 ccd_id               Int2                5:5,7:7 
   4 node_id              Int2                0:3 
   5 chipx                Int2                1:1024 
   6 chipy                Int2                1:1024 
   7 tdetx                Int2                1:8192 
   8 tdety                Int2                1:8192 
   9 detx                 Real4                   0.50:     8192.50 
  10 dety                 Real4                   0.50:     8192.50 
  11 phas                 Int2                -4096:4095 
  12 pha                  Int4                0:36855 
  13 pha_ro               Int4                0:36855 
  14 energy               Real4                     500.0:     7000.0 
  15 pi                   Int4                1:1024 
  16 fltgrade             Int2                0:255 
  17 grade                Int2                0:0,2:2,3:3,4:4,6:6 
 --- Component 2 --- 
 ...
 

which shows a MASK has now been introduced into the subspace. . The (MASK) identifies another block within the dataset where the mask information is stored. We can display this when we list the blocks in the dataset:

unix% dmlist image_with_mask blocks

--------------------------------------------------------------------------------
Dataset: image_with_mask
--------------------------------------------------------------------------------
 
     Block Name                          Type         Dimensions
--------------------------------------------------------------------------------
Block    1: EVENTS_IMAGE                   Image      Int4(2092x4308)
Block    2: MASK                           Image      Byte(2092x4308)

The MASK block is the same size and dimension as the image, and has been stored as a 1-byte image. This can be visualized using one of the the above techniques:

unix% ds9 image_with_mask -scale log -block 4 -scale log -pan to 4077.3 4269.7 physical \
  -mask color green -mask transparency 60 -mask "image_with_mask[mask]"

This produces the same image as in Figure 8.


Use masked image to create radial profile

The masked image can then be used to create a radial profile

unix% dmextract "image_with_mask[bin sky=annulus(4077.3,4269.7,0:600:10)]" \
  exp="abell2626_broad_thresh.expmap" \
  out=radial_profile_with_mask_from_image mode=h op=generic clob+

As before, the annulus are intersected with the spatial filters encoded in the input file's subspace to compute the correct AREA and thus flux values. The output is shown in Figure 11 and is compared with the output using regions.

Figure 11: Radial profile of broad-band image with sources excluded using a pixel mask

[Thumbnail image: Radial profile of broad-band image with sources excluded using a pixel mask]

[Version: full-size]

[Print media version: Radial profile of broad-band image with sources excluded using a pixel mask]

Figure 11: Radial profile of broad-band image with sources excluded using a pixel mask

The radial profile created by excluding point sources using pixel masks (black) and using regions (red) are shown. As expected, they are nearly identical. The difference between the two is shown in the bottom panel.

from pycrates import read_file
import matplotlib.pylab as plt
import numpy as np

a = read_file("radial_profile_with_mask_from_image")
ax = a.get_column("cel_rmid").values
ay = a.get_column("cel_flux").values
ae = a.get_column("cel_flux_err").values

b = read_file("radial_profile_with_regions_from_image")
bx = b.get_column("cel_rmid").values
by = b.get_column("cel_flux").values
be = b.get_column("cel_flux_err").values

grid = plt.GridSpec(4, 1, wspace=0.0, hspace=0.0)
plt.subplot(grid[0:3,0])
plt.errorbar(ax,ay,yerr=ae, drawstyle="steps-mid", ecolor="gray")
plt.errorbar(bx,by,yerr=be, drawstyle="steps-mid", ecolor="gray", color="red")
plt.yscale("log")
plt.ylabel("CEL_FLUX (photon/cm**2/arcsec**2/s)")
plt.title("Abell 2626")

plt.subplot(grid[3,0])
delta = ay - by
quad = np.sqrt(ae*ae + be*be)
plt.errorbar(ax,delta,yerr=quad, drawstyle="steps-mid", ecolor="gray")
plt.xlabel("CEL_R (arcsec)")
plt.ylabel(r"$\Delta$ CEL_FLUX")

The difference between these two methods comes down to how the edge of the image is treated. The field-of-view region creates a loose boundary around the edge of the detector, whereas the exposure map provides an accurate footprint based on the detector geometry and aspect solution. The effect is only seen when the radii intersect the edge of the detector (large radii). We can see this if we diff the outputs

unix% dmdiff radial_profile_with_regions_from_image"[cols cel_flux]" radial_profile_with_mask_from_image"[cols cel_flux]" sub- key- 
Infile 1:  radial_profile_with_regions_from_image[cols cel_flux]
Infile 2:  radial_profile_with_mask_from_image[cols cel_flux]


TABLE NAME: HISTOGRAM

-----------------------
TABLE VALUE DIFFERENCES
-----------------------

Message:                             Row:        Column:              Value(s):                    Diff:
--------                             ----        -------              ---------                    -----
Values are not equal                  42        CEL_FLUX 1.22601768094434e-08  1.22607280200636e-08 +5.51211e-13 (+0.0045 %)
Values are not equal                  44        CEL_FLUX 1.14021673758562e-08  1.14026103782445e-08 +4.43002e-13 (+0.00389 %)
Values are not equal                  51        CEL_FLUX 8.54549095788864e-09  8.54572658892842e-09 +2.35631e-13 (+0.00276 %)
Values are not equal                  53        CEL_FLUX 8.14868724708034e-09  8.14868724708034e-09 +1.65436e-24 (+2.03e-14 %)
Values are not equal                  55        CEL_FLUX 7.85983505037149e-09  7.86015345515113e-09 +3.18405e-13 (+0.00405 %)
Values are not equal                  57        CEL_FLUX 7.68046043214958e-09  7.68070245333169e-09 +2.42021e-13 (+0.00315 %)
Values are not equal                  59        CEL_FLUX 7.08735096784064e-09  7.08735096784064e-09 -8.27181e-25 (-1.17e-14 %)
Values are not equal                  60        CEL_FLUX 6.91583037576813e-09  6.91583037576813e-09 +8.27181e-25 (+1.2e-14 %)

The differences are very small.


Alternative edge treatment: intersecting masks

As mentioned above there are several different ways to encode the edge of the detector into the mask. Another alternative approach is to use the same FOV region file to create its own mask file similar to how the source-free mask was created, and then intersect those two masks together.

This approach starts with the same unit image created above, but now filters it with the FOV file

unix% dmcopy "ones.fits[sky=region(acisf16136_repro_fov1.fits)][opt full]" fov.mask clob+

Figure 12: Mask created from FOV file

[Thumbnail image: mask image created from fov file]

[Version: full-size]

[Print media version: mask image created from fov file]

Figure 12: Mask created from FOV file

The mask image created by filtering the unit image with the field of view file similar to how the source-mask was created.

The mask in Figure 12 now just needs to be intersected with the mask in Figure 6. In terms of logical operations, intersection is accomplished by multiplication, so the two mask files just need to be multiplied together

unix% dmimgcalc no_srcs.mask fov.mask no_srcs_fov.mask op=mul clob+

The intersected mask is shown in Figure 13 and the difference in the masks is shown in Figure 14.

Figure 13: Result of FOV mask intersected with source-free mask

[Thumbnail image: Result of FOV mask intersected with source-free mask]

[Version: full-size]

[Print media version: Result of FOV mask intersected with source-free mask]

Figure 13: Result of FOV mask intersected with source-free mask

The result of intersecting the FOV mask with the source free mask looks very similar to the mask created in Figure 7 where the edge was created from the exposure map.

Figure 14: Difference in FOV vs. exposure map based masks

[Thumbnail image: Difference in FOV vs. exposure map based masks]

[Version: full-size]

[Print media version: Difference in FOV vs. exposure map based masks]

Figure 14: Difference in FOV vs. exposure map based masks

The difference in the mask created using the FOV region file and by applying a threshold to the exposure map. Red pixels included in the FOV region based mask, blue pixels (none visible) are only included in the exposure map based mask. Grey pixels are the same in both.

unix% dmimgcalc no_srcs_fov.mask no_srcs_with_fov.mask diff_fov.mask op=sub clob+

The radial profile can be extracted using this mask in the same way

unix% dmextract "abell2626_broad_thresh.img[sky=mask(no_srcs_fov.mask)][bin sky=annulus(4077.3,4269.7,0:600:10)]" \
  exp="abell2626_broad_thresh.expmap" \
  out=radial_profile_with_mask_from_image2 mode=h op=generic clob+
# dmextract (CIAO): WARNING: Input file, "abell2626_broad_thresh.img[sky=mask(no_srcs_fov.mask)]", has no rows in it.

The WARNING here can be ignored; it because in this example the mask filter is applied on-the-fly rather than creating a separate input file.

This output is then nearly identical to the profile created from regions

unix% dmdiff radial_profile_with_regions_from_image"[cols cel_flux]" radial_profile_with_mask_from_image2"[cols cel_flux]"  key- sub-
Infile 1:  radial_profile_with_regions_from_image[cols cel_flux]
Infile 2:  radial_profile_with_mask_from_image2[cols cel_flux]


TABLE NAME: HISTOGRAM

-----------------------
TABLE VALUE DIFFERENCES
-----------------------

Message:                             Row:        Column:              Value(s):                    Diff:
--------                             ----        -------              ---------                    -----
Values are not equal                  42        CEL_FLUX 1.22601768094434e-08  1.22607280200636e-08 +5.51211e-13 (+0.0045 %)
Values are not equal                  44        CEL_FLUX 1.14021673758562e-08  1.14026103782445e-08 +4.43002e-13 (+0.00389 %)
Values are not equal                  51        CEL_FLUX 8.54549095788864e-09  8.54572658892842e-09 +2.35631e-13 (+0.00276 %)
Values are not equal                  55        CEL_FLUX 7.85983505037149e-09  7.86015345515113e-09 +3.18405e-13 (+0.00405 %)
Values are not equal                  57        CEL_FLUX 7.68046043214958e-09  7.68070245333169e-09 +2.42021e-13 (+0.00315 %)

The differences are due to small changes in how the pixelated area is computed for these complex regions.


Why use pixel masks?

As the results between the region based filtering and the pixel mask based filtering are equivalent, users may question why use one approach or the other?

One important difference not shown here directly is that the mask-filtered dmextract command only took 30 seconds to complete, compared to the traditional region-based filtering which took 60 seconds. While this may seem like a modest increase, consider that the mask based filter is fixed -- based only on the number of pixels in the mask -- whereas the regions based filter is based on the number of sources being excluded. Therefore, if this example had 50 or 500 sources, the mask based filter would continue to run in the same, fixed amount of time. In contrast it is expected that the region based filtering would at best scale linearly, O(N), with the number of sources. If the region were not inverted manually, the run-time would grow exponentially.

Another advantage to using masks is how easily they are to manipulate. For example to invert the mask, it was a simple dmimgcalc command. To invert the region required creating intermediate files and somewhat complicated (at best arcane) set of commands.


Tables: Pixel masks work the same way with event files

It is easy to visualize the use of pixel masks with images where there must be a pixel-by-pixel match between the image being filtered and the image being used as a mask. Pixel masks work the same way with event files too. The area/extent of the pixel is used to filter the rows in the event file the same way as pixels are filtered in images. The same mask() syntax is used

unix% dmcopy "acisf16136_repro_evt2.fits[sky=mask(no_srcs_with_fov.mask),energy=500:7000]" events_point_srcs_removed_mask clob+

The mask is stored in the subspace in the same way:

unix% dmlist events_point_srcs_removed_mask subspace
 
--------------------------------------------------------------------------------
Data subspace for block EVENTS: Components: 4 Descriptors: 16 
--------------------------------------------------------------------------------
 
 --- Component 1 --- 
   1 time                 Real8               TABLE GTI7
                                              
                                              498642950.2601068020:498755270.7022138238
   2 expno                Int4                3:35762 
   3 ccd_id               Int2                7:7 
   4 node_id              Int2                0:3 
   5 chip                 [ 1] chipx          1:1024 
   5 chip                 [ 2] chipy          1:1024 
   6 tdet                 [ 1] tdetx          1:8192 
   6 tdet                 [ 2] tdety          1:8192 
   7 det                  [ 1] detx               0.50:     8192.50 
   7 det                  [ 2] dety               0.50:     8192.50 
   8 sky                  Real4               TABLE MASK
                                              MASK(MASK)
                                              Field area = 6.71089e+07 Region area = 4.35209e+06

                          [ 1] x              
                                                   3115.50:     5207.50
                          [ 2] y              
                                                   1762.50:     6070.50
   9 phas                 Int2                -4096:4095 
  10 pha                  Int4                0:36855 
  11 pha_ro               Int4                0:36855 
  12 energy               Real4                     500.0:     7000.0 
  13 pi                   Int4                1:1024 
  14 fltgrade             Int2                0:255 
  15 grade                Int2                0:0,2:2,3:3,4:4,6:6 
  16 status               Bit                 
  ...

The mask is stored in a separate extension

unix% dmlist events_point_srcs_removed_mask blocks

--------------------------------------------------------------------------------
Dataset: events_point_srcs_removed_mask
--------------------------------------------------------------------------------
 
     Block Name                          Type         Dimensions
--------------------------------------------------------------------------------
Block    1: PRIMARY                        Null        
Block    2: EVENTS                         Table        16 cols x 326105   rows
Block    3: GTI7                           Table         2 cols x 1        rows
Block    4: GTI5                           Table         2 cols x 1        rows
Block    5: GTI6                           Table         2 cols x 5        rows
Block    6: GTI8                           Table         2 cols x 5        rows
Block    7: MASK                           Image      Byte(2092x4308)

The radial profile is then also extracted in the same way

unix% dmextract events_point_srcs_removed_mask"[bin sky=annulus(4077.3,4269.7,0:600:10)]" \
  exp="abell2626_broad_thresh.expmap" \
  out=radial_profile_with_mask mode=h op=generic clob+

This is where things start to differ. Since event locations have continuous, real-values, the individual events which are included or excluded using the filters is different compared to using image pixel. This is because for images only the discrete, center pixel location is tested for being inside the region. Therefore we expect to see some differences in the radial profile extracted from the image compared to the radial profile extracted from an event file. This is illustrated in Figure 15.

Figure 15: Radial profile comparison: image vs events

[Thumbnail image: diff in radial profiles for images vs events]

[Version: full-size]

[Print media version: diff in radial profiles for images vs events]

Figure 15: Radial profile comparison: image vs events

A comparison of radial profile for the broad-band Abell 2626 cluster with sources excluded using pixels masks. The black curve is the radial profile extracted from the binned image. The red curve is the radial profile extracted from the events file.

The bottom panel shows the absolute difference in the histograms. The plot was made with these commands:

from pycrates import read_file
import matplotlib.pylab as plt
import numpy as np

a = read_file("radial_profile_with_mask_from_image")
ax = a.get_column("cel_rmid").values
ay = a.get_column("cel_flux").values
ae = a.get_column("cel_flux_err").values

b = read_file("radial_profile_with_mask")
bx = b.get_column("cel_rmid").values
by = b.get_column("cel_flux").values
be = b.get_column("cel_flux_err").values

grid = plt.GridSpec(4, 1, wspace=0.0, hspace=0.0)
plt.subplot(grid[0:3,0])
plt.errorbar(ax,ay,yerr=ae, drawstyle="steps-mid", ecolor="gray")
plt.errorbar(bx,by,yerr=be, drawstyle="steps-mid", ecolor="gray", color="red")
plt.yscale("log")
plt.ylabel("CEL_FLUX (photon/cm**2/arcsec**2/s)")
plt.title("Abell 2626")

plt.subplot(grid[3,0])
delta = ay - by
quad = np.sqrt(ae*ae + be*be)
plt.errorbar(ax,delta,yerr=quad, drawstyle="steps-mid", ecolor="gray")
plt.xlabel("CEL_R (arcsec)")
plt.ylabel(r"$\Delta$ CEL_FLUX")

The differences shown in Figure 15 are not in any way due to pixel masks but is purely the result of comparing binned image data to event list data.



Summary

In this thread users were introduced to pixel masks filters. Users learned

  • What a pixel mask is.
  • How to create a pixel mask.
  • How to manipulate and edit a pixel mask
    • Using region filter
    • By applying thresholds
    • Inverting a mask
    • Intersecting masks
  • How to visualize pixel masks
  • How to use a pixel mask to filter a dataset (both images and event files)
  • How pixel masks are stored in the datasets subspace
  • Why pixel masks are useful

There are countless other uses of pixel masks filters beyond just extracting radial profiles and equally countless other ways to edit and manipulate pixel masks beyond those methods shown in this introductory thread.


History

01 Apr 2018 Initial version.
10 Dec 2018 Updated for CIAO 4.11. fluximage now creates the PSF map.
03 Apr 2019 Updated to use matplotlib for plotting.