Using Pixel Mask Filters
![[CXC Logo]](../../imgs/cxc-logo.gif){kind=logo}
CIAO 4.17 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:
Related Links:
Last Update: 25 Jan 2022 - Review for CIAO 4.14. Updated for Repro5 and CALDB 4.9.6.
Contents
- Getting Started
- Traditional approach: using CIAO regions
- Alternative approach: using pixel masks
- Summary
- History
-
Images
- Figure 1: Broad band image of Abell 2626
- Figure 2: Broad band image of Abell 2626 with wavdetect sources
- Figure 3: Broad band image with point source excluded
- Figure 4: Radial profile from broad band image with excluded sources using regions
- Figure 5: Unit mask filtered with source list
- Figure 6: Inverted unit mask filtered with source list
- Figure 7: Pixel mask with excluded sources and detector edges
- Figure 8: Visualize mask with ds9
- Figure 9: Visualize mask with matplotlib
- Figure 10: Visualizing mask with dmimg2jpg
- Figure 11: Radial profile of broad-band image with sources excluded using a pixel mask
- Figure 12: Mask created from FOV file
- Figure 13: Result of FOV mask intersected with source-free mask
- Figure 14: Difference in FOV vs. exposure map based masks
- Figure 15: Radial profile comparison: image vs events
Getting Started
Download the sample data: 16136 (Abell 2626)
unix% download_chandra_obsid 16136
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
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]](ds9_03.png)
[Version: full-size]
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]](../../imgs/important.png)
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]](ds9_01.png)
[Version: full-size]
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(4111.96,2756,27.5003,24.0959,115.223)&!Ellipse(3951.3,3719.56,6.42102,3.40998,126.953) &!Ellipse(4403.63,3773.54,6.51762,4.9292,164.678)&!Ellipse(3790.33,3882.92,11.4101,7.56801,81.8236) # POS(X,Y) SHAPE R[2] ROTANG[2] COMPONENT NaN NaN Field NaN NaN NaN NaN 1 4111.9611650485 2756.0 !Ellipse 27.5003337860 24.0958766937 115.2232360840 NaN 1 3951.3040935673 3719.5614035088 !Ellipse 6.4210228920 3.4099793434 126.9527511597 NaN 1 4403.6296296296 3773.5432098765 !Ellipse 6.5176234245 4.9291977882 164.6782531738 NaN 1 3790.3269230769 3882.9230769231 !Ellipse 11.4101200104 7.5680093765 81.8236083984 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]](ds9_02.png)
[Version: full-size]
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(3088.5,1706.5,5188.5,6031.5) &!Ellipse(4111.96,2756,27.5003,24.0959,115.223)&!Ellipse(3951.3,3719.56,6.42102,3.40998,126.953) &!Ellipse(4403.63,3773.54,6.51762,4.9292,164.678)&!Ellipse(3790.33,3882.92,11.4101,7.56801,81.8236) &!Ellipse(4130.26,3905.47,3.45655,2.50034,101.599)&!Ellipse(3732.98,3912.9,11.2455,6.21996,60.9637) &!Ellipse(3850.06,3917.7,9.23236,7.19514,13.5605)&!Ellipse(4398.22,3940.63,8.22292,5.29508,148.626) &!Ellipse(3582.72,4135.44,9.6198,6.23803,76.9147)&!Ellipse(3822.38,4173.19,7.54269,5.33875,110.211) &!Ellipse(3967.88,4194.16,5.7766,5.44011,43.025)&!Ellipse(4047.87,4216.18,5.5322,5.06669,49.6533) &!Ellipse(4215.2,4217.3,6.75448,5.53103,15.551)&!Ellipse(4465.41,4231.8,8.63881,3.58085,179.454) &!Ellipse(4202.73,4247.17,4.61904,4.43388,64.3435)&!Ellipse(3803.28,4311.35,6.40808,5.75475,55.1055) &!Ellipse(4268.75,4346.63,10.2877,7.38126,90.4414)&!Ellipse(4208.88,4396.98,4.67631,4.39647,68.2531) &!Ellipse(3610.93,4532.4,13.6672,8.49983,62.2819)...[truncated] 1 sky Real8 Field area = 9.0825e+06 Region area = 9.03182e+06 1 sky [ 1] x 3088.50: 5188.50 1 sky [ 2] y 1706.50: 6031.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 ...
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]](../../imgs/note.png)
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(3088.5,1706.5,5188.5,6031.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]](../../imgs/tip.png)
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]](chips_01.png)
[Version: full-size]
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")
plt.show()
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 @: ( 3089 1707 ) max: 1 @: ( 3089 1707 ) mean: 1 sum: 9082500 good: 9082500 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]](ds9_10.png)
[Version: full-size]
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]](ds9_11.png)
[Version: full-size]
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.
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]](ds9_12.png)
[Version: full-size]
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]](ds9_13.png)
[Version: full-size]
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
Figure 9: Visualize mask with matplotlib
![[Thumbnail image: using python to view mask]](chips_40.png)
[Version: full-size]
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)
plt.show()
Figure 10: Visualizing mask with dmimg2jpg
![[Thumbnail image: image with mask made with dmimg2jpg]](dmimg2jpg_40.jpg)
[Version: full-size]
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]](chips_21.png)
[Version: full-size]
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")
plt.show()
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] ---------------------------------------------------------------------- Compare Tables ---------------------------------------------------------------------- Compare Table Structure: Block name: HISTOGRAM Compare Column Details: Compare Virtual Column Details: Compare Column Data: Column: Row: Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ CEL_FLUX 8 Values are not equal 1.54796289051984e-07 1.54947194223538e-07 +1.50905e-10 (+0.0975%) CEL_FLUX 19 Values are not equal 4.87414599423975e-08 4.87452272412216e-08 +3.7673e-12 (+0.00773%) CEL_FLUX 34 Values are not equal 1.85854066666275e-08 1.85860940866614e-08 +6.8742e-13 (+0.0037%)
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]](ds9_30.png)
[Version: full-size]
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]](ds9_31.png)
[Version: full-size]
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]](ds9_32.png)
[Version: full-size]
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] ---------------------------------------------------------------------- Compare Tables ---------------------------------------------------------------------- Compare Table Structure: Block name: HISTOGRAM Compare Column Details: Compare Virtual Column Details: Compare Column Data: Column: Row: Message: Value(s): Diff: ---------------- -------------- -------------------------------------- ---------------------------------- ------------------------ CEL_FLUX 8 Values are not equal 1.54796289051984e-07 1.54947194223538e-07 +1.50905e-10 (+0.0975%) CEL_FLUX 19 Values are not equal 4.87414599423975e-08 4.87452272412216e-08 +3.7673e-12 (+0.00773%) CEL_FLUX 34 Values are not equal 1.85854066666275e-08 1.85860940866614e-08 +6.8742e-13 (+0.0037%)
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]](chips_11.png)
[Version: full-size]
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")
plt.show()
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. |
| 25 Jan 2022 | Review for CIAO 4.14. Updated for Repro5 and CALDB 4.9.6. |