Automasking Guide CASA 6.5.4
Last checked on CASA Version 6.5.4
Introduction
This guide features CARTA, the “Cube Analysis and Rendering Tool for Astronomy,” which is the new NRAO visualization tool for images and cubes. The CASA viewer (imview) has not been maintained for a few years and will be removed from future versions of CASA. We strongly recommend using CARTA, as it provides a much more efficient, stable, and feature rich user experience. A comparison of the CASA viewer and CARTA, as well as instructions on how to use CARTA at NRAO, is provided in the CARTA section of the CASA docs.
Starting with CASA 5.1, a new algorithm has been incorporated into tclean to automatically mask regions during the cleaning process. Referred to as ‘auto-multithresh’, this algorithm is intended to mimic what an experienced user would do when manually masking images. It can be used by setting the usemask option in tclean to ’auto-multithresh’.
A full description of the auto-multithresh algorithm is given in the CASA docs and Kepley et al. 2020. Here we briefly describe the fundamental features of the algorithm. First, the algorithm identifies regions that are either above a signal-to-noise limit or a sidelobe level (as defined by the user), whichever is higher. If the mask regions are smaller than some fraction of the beam, they are removed (“pruned”) to avoid including spurious noise peaks in the masked region. If there is extended low signal to noise emission surrounding the initial mask then this emission is included in the mask by expanding the noise/sidelobe mask into the low signal-to-noise regions using a process called binary dilation. Finally, the mask is expanded to include a buffer region around the masked emission by convolving the mask by a Gaussian and then retaining only regions that are above n% of the peak in the smoothed mask. Absorption regions can also be masked using a method similar to that for the initial threshold mask, but are not pruned or extended into low signal-to noise regions. Note that for cubes, the algorithm masks each channel independently. The algorithm operates on the residual image at the beginning of every minor cycle, so that the mask updates as the clean progresses.
As of CASA 5.6, auto-multithresh now functions with polarization data. It applies the same algorithms to the Stokes QUV images as used for the Stokes I image. This means that the full masking process is applied to the positive emission (including the prune and grow steps), but that the masking of the negative emission only includes the initial threshold mask (no prune or grow).
Parameters
The parameters for the auto-multithresh algorithm are defined in terms of fundamental image properties (e.g., fraction of beam, signal-to-noise, etc) rather than exact numerical values (e.g., 3”, 6 mJy/beam). This feature allows the same set of parameters to be used for multiple images. The auto-multithresh algorithm calculates the numerical value of each parameter internally by estimating the noise in the residual image (currently done using the median absolute deviation) and the beam size and sidelobe level calculated by tclean internally. The numerical values are recalculated at the beginning of every minor cycle and may change as the clean progresses.
There are five primary auto-multithresh subparameters that control what thresholds the algorithm uses and how regions are removed from the mask via pruning.
- noisethreshold (type: double) -- sets the signal-to-noise threshold above which significant emission is masked during the initial round of mask creation. Note that either noisethreshold or sidelobethreshold is used depending on which threshold is higher.
- sidelobethreshold (type: double) -- sets a threshold based on the sidelobe level above which significant emission is masked during the initial round of mask creation. Note that either noisethreshold or sidelobethreshold is used depending on which threshold is higher.
- lownoisethreshold (type: double) -- sets the threshold into which the initial mask (which is determined by either noisethreshold or sidelobethreshold) is expanded in order to include low signal-to-noise regions in the mask.
- minbeamfrac (type: double) -- sets the minimum size a region must be to be retained in the mask. The parameter is specified as a fractional beam size. Any masks smaller than this will be removed (or pruned) from the mask. Note that this parameter is used to control the pruning for both the initial threshold mask and the low signal-to-noise mask.
- negativethreshold (type: double) -- sets the signal-to-noise threshold for absorption features to be masked. Note that any absorption features that are masked are not pruned or expanded into low signal-to-noise regions. This parameter should have a positive value, i.e., negativethreshold=5.0 will create a mask for values less than -5.0 sigma.
Auto-multithresh has several additional subparameters that usually do not need to be changed from their default values:
- pbmask (type: double) -- primary beam cutoff
- smoothfactor (type: double) -- controls how much to smooth the initial mask. We do not recommend changing this from its default value.
- cutthreshold (type: double) -- controls what regions of the smoothed mask are retained to form the final mask. We do not recommend changing this from its default value.
- growiterations (type: int) -- controls the maximum number of iterations that binary dilation performs. The binary dilation may halt earlier than the maximum number of iterations if there is no change in the mask. Lowering this value can reduce the computational time needed to expand the mask into low signal-to-noise regions at the expense of potentially not expanding the mask all the way to its edges. A value between 75 and 100 is usually adequate.
- dogrowprune (type: bool) -- allows you to turn off the pruning of the low signal-to-noise mask, which speeds up masking for images and cubes with complex low signal-to-noise emission.
- minpercentchange (type: double) -- allows you to stop masking when the mask changes by less than n% and the cyclethreshold is equal to the threshold for the previous cycle. This should be used with care and is turned off (minpercentchange=-1) by default.
- verbose (type: bool) -- turns on/off per-channel logging information on the masks (default = False)
Additionally, the following parameter should be considered when using automasking (this not a subparameter of auto-multithresh).
- fastnoise (type: bool) -- If True, estimate the noise using the median absolute deviation. If False, estimate the noise using the chauvenet (no mask) or the median absolute deviations of pixels outside the mask but inside the primary beam mask. The latter estimation may be more accurate when emission covers a significant portion of the field of view (default=True).
The default automasking parameters in tclean have been set to values appropriate for the ALMA 12m array in its more compact configurations (i.e. not long baseline data). These values depend on the PSF (i.e., uv-coverage) of the data. Data with PSFs significantly different than the fiducial ALMA 12m arrays (e.g. 7m data) will likely need the auto-multithresh parameters to be modified in order to achieve the best mask. The default pipeline parameters presented here (Table below) have been extensively tested and most likely will work for your data set, although you may wish to optimize them for your particular imaging case. Data from other instruments (e.g., VLA, ATCA, etc) has been shown to work with this algorithm in limited testing, but the parameter values may need to be optimized.
Table of Standard Values
The values presented in this table are the standard values that the pipeline uses when it performs automasking and have been extensively tested. The values listed for combined data in the last column are tentative and should serve as a starting point for your analysis. This table can be found in the Pipeline User's Guide section 9.39 (hif_makeimages).
| Parameter | 7m (continnum/line) | 12m (b75 < 300m) | 12m (300m < b75 < 400m) | 12m (b75 > 400m) | 12m + 7m combined TENTATIVE |
|---|---|---|---|---|---|
| noisethreshold | 5.0 | 4.25 | 5.0 | 5.0 | 4.25 |
| sidelobethreshold | 1.25 | 2.0 | 2.0 | 2.5 | 2.0 |
| lownoisethreshold | 2.0 | 1.5 | 1.5 | 1.5 | 1.5 |
| minbeamfrac | 0.1 | 0.3 | 0.3 | 0.3 | 0.3 |
| negativethreshold | 0.0 | 0.0 (continuum) 15.0 (line) |
0.0 (continuum) 7.0 (line) |
0.0 (continuum) 7.0 (line) |
0.0 |
| fastnoise | False | False | False | True | False |
The term “b75” refers to the 75th percentile of baselines. You can find this value in the pipeline weblog by going to the antenna configuration page and clicking on the “baselines” tab, the table lists the baseline length in increasing distance and the percentile. If you don’t have a weblog for your data, you may use the Analysis Utilities task au.getBaselineStats. You can obtain the Analysis Utilities tasks by following the Analysis Utilities CASA Guide.
The 75% baseline split is only valid for ALMA data. For data from other telescopes, evaluate the uv-coverage of the observation. If the uv-coverage of a data set is good, then the PSF will have relatively low sidelobes and the default 12m parameters are likely a good start. We have had success in limited testing using the short baseline ALMA parameters for non-snapshot VLA observations. If the uv-coverage of a data set is poor, then the PSF will have higher sidelobes and the 7m array parameters are probably a better start. An example of this case would be snapshot VLA observations or ATCA observations.
If you are using pipeline imaging tasks, these automasking parameters are also able to be defined within those task calls. See the Pipeline Tasks Reference Manual section 2.12.
Guidelines for Optimizing auto-multithresh Parameters
- Start with the pipeline parameters for the data with the most similar parameters to your data set. When in doubt, choose the default short baseline 12m parameters.
- For cubes, use a small (~10 channel) cube with representative emission as a test case to avoid long run times.
- Change one parameter at a time and keep track of your previous images so that you can directly compare different runs.
- Inspect the log output. The ‘auto-multithresh’ algorithm reports extensive information to the logger on what it is doing including what thresholds are being used, how many regions are being pruned, etc. If using CASA 5.3+, set verbose=True to get this information.
- Save the masks for each major cycle. To save the intermediate masks, type the following on the casa command line:
#In CASA
import os
os.environ['SAVE_ALL_AUTOMASKS']="true"
This feature does NOT work when tclean is run with parallel=True on cubes.
- You can also inspect the intermediate masks by setting interactive=True. However, the resulting masks will NOT necessarily be the same as masks generated with interactive=False because how often the major cycle is triggered can be different in the two cases. Set cycleniter=niter when using non-interactive mode to get identical behavior.
Exploring the auto-multithresh Parameters
This section should serve as a guide to help you gain a more intuitive understanding of how the automasking algorithm works and understand better what each parameter does during the automasking process. Additionally, we will show you what can happen when you set a parameter to a non-optimal value and end up masking noise features or not masking real emission.
Obtain the Data
To demonstrate the use of the automasking parameters, we will be using the TW Hydra dataset which is a small dataset used in the “First Look” CASA Guides. Specifically, we will use a flagged, calibrated, self-calibrated, and continuum subtracted Band 7 Measurement Set that is created as part of the First_Look_at_Line_Imaging guide. Download and untar: twhya_selfcal.ms.contsub.tar (~300 MB).
#In bash
tar xvf twhya_selfcal.ms.contsub.tar
Almost all of the values used in our following tclean commands remain constant throughout this section, and are explained in and taken from the First_Look_at_Line_Imaging guide. The main difference is that we are using the usemask=’auto-multithresh’ parameter to use the automasking feature and will be adjusting the auto-masking subparameters to highlight their use.
Make a Dirty Cube
First we make a dirty cube to figure out the value for threshold and get an idea of what the emission/noise looks like. It's always a good idea to perform a listobs first to get the basic information of the dataset.
#In CASA
listobs('twhya_selfcal.ms.contsub', listfile='twhya_selfcal.ms.contsub.listobs.txt')
Observer: cqi Project: uid://A002/X327408/X6f
Observation: ALMA
Data records: 44772 Total elapsed time = 4268.11 seconds
Observed from 19-Nov-2012/07:56:23.5 to 19-Nov-2012/09:07:31.6 (UTC)
ObservationID = 0 ArrayID = 0
Date Timerange (UTC) Scan FldId FieldName nRows SpwIds Average Interval(s) ScanIntent
19-Nov-2012/07:56:23.5 - 08:02:11.3 12 5 TW Hya 7616 [0] [6.05] [OBSERVE_TARGET#ON_SOURCE]
08:08:09.6 - 08:13:57.3 16 5 TW Hya 8442 [0] [6.05] [OBSERVE_TARGET#ON_SOURCE]
08:19:53.9 - 08:25:41.7 20 5 TW Hya 8389 [0] [6.05] [OBSERVE_TARGET#ON_SOURCE]
08:32:00.5 - 08:37:48.2 24 5 TW Hya 8409 [0] [6.05] [OBSERVE_TARGET#ON_SOURCE]
08:43:45.6 - 08:49:33.4 28 5 TW Hya 8514 [0] [6.05] [OBSERVE_TARGET#ON_SOURCE]
09:05:15.6 - 09:07:31.6 36 5 TW Hya 3402 [0] [6.05] [OBSERVE_TARGET#ON_SOURCE]
(nRows = Total number of rows per scan)
Fields: 1
ID Code Name RA Decl Epoch SrcId nRows
5 none TW Hya 11:01:51.796000 -34.42.17.36600 J2000 4 44772
Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)
SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) BBC Num Corrs
0 ALMA_RB_07#BB_2#SW-01#FULL_RES 384 TOPO 372533.086 610.352 234375.0 372649.9688 2 XX YY
With the information above and from what we know of TW Hydra from the First_Look_at_Line_Imaging guide, we get the basic information for our tclean call to create the dirty cube (set niter=0 to perform no cleaning).
#In CASA
os.system('rm -rf twhya_dirtycube.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_dirtycube',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=0)
Determine Clean Threshold
In this tutorial we will clean non-interactively, and thus need to set a threshold to avoid over-cleaning emission. This is true for all non-interactive cleaning, and also a good idea for interactive cleaning (the threshold parameter is not a subparameter of auto-multithresh). After tclean has finished, you should now open the dirty image in CARTA. If using NRAO machines, you can open a new terminal tab, cd to the working directory, then type:
#In bash
carta --no_browser
Copy the output URL into a browser to view your CARTA session. Select and load twhya_dirtycube.image. Selecting a region in a line free channel (e.g., channel 3) we find that the RMS is ~29 mJy/beam. We choose the threshold to be 3x the RMS: threshold='87mJy' (The "/beam" is implied).
Base Parameters
To figure out what automasking parameter values we should use for this data, we run au.getBaselineStats which will tell us the "b75" value for our data. If you are not using an NRAO machine, follow the Analysis_Utilities guide.
#In CASA
import analysisUtils as au
au.getBaselineStats('twhya_selfcal.ms.contsub')
Found 1 obsIDs Found 325 baselines (26 stations) Unprojected lengths: min=15.063874, max=374.719747, rms=162.602204 m number=325, min=15.06m, max=374.72m, median=139.34m, mean=142.78m, std=77.81m 20%ile=69.4m 25%ile=79.0m, 30%ile=87.6m, 75%ile=197.4m, 90%ile=250.2m
From the output above (which was printed to the terminal window), we see that b75 (indicated as 75%ile) is 197.4m, so we will use the pipeline default values in the "12m (b75 < 300m)" column of the Table of Standard Values above. This will give us a mask to compare our testing against which will be displayed as black contours in all of the following figures.
#In CASA
os.system('rm -rf twhya_base_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_base_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)
- Open the following three images in CARTA (with File > Open and ctrl+click, or one at a time with File > Append):
twhya_base_params.image
twhya_base_params.residual
twhya_base_params.mask - In the Image List, lock all 3 images in XY (spatial matching), Z (spectral matching), and R (raster scaling matching). If .image is not already set as the reference for matching, right click the image name and select "Set as ... reference." To reorder images in the list, click an index, then click-and-drag the index.
- Click the settings gear of the Image Viewer, go to the Global tab, enable multi-panel, and choose 2 columns and 1 row. Adjust the images to an appropriate position (Zoom to fit is recommended).
- Use the Animator to page through the channels. Notice that the mask image is blank for all channels except 5 and 7. For this example we choose channel 7.
- Adjust the raster scaling in Render Configuration.
- Open the Contour dialog. Unlock data source. Choose the mask. The Histogram, Generator mode, and Parameters provide different ways to calculate Levels, however for a mask contour we can ignore these. A mask image has value 1.0 where the mask exists, and 0.0 elsewhere (see this by moving your cursor around a mask image). Enter 0.5 directly into the Levels box, then click apply. This displays a contour outline of the mask on all 3 images. Go to the Styling tab to adjust color and thickness.
With verbose=True set, automasking reports its progress per-channel in the logger window. Below is the logger output for the first major clean cycle using the base parameters above.
- The minimum mask size and sidelobeLevel are calculated in order to define different thresholds based on the user input.
----------------------------------------------------------- Run Major Cycle 1 ------------------------------------- Absolute Peak residual within mask : 0.190839, over full image : 0.190839 [twhya_base_params] Setting up an auto-mask within PB mask limit ---------------------------------------------------- Run Automask --------------------------------------------- Set Deconvolution Options for [twhya_base_params] : hogbom Absolute Peak residual within mask : 0.190839, over full image : 0.190839 Processing channels in range [0, 14] Generating AutoMask [C0] Using Chauvenet algorithm for the image statistics [C1] Using Chauvenet algorithm for the image statistics [C2] Using Chauvenet algorithm for the image statistics [C3] Using Chauvenet algorithm for the image statistics [C4] Using Chauvenet algorithm for the image statistics [C5] Using Chauvenet algorithm for the image statistics [C6] Using Chauvenet algorithm for the image statistics [C7] Using Chauvenet algorithm for the image statistics [C8] Using Chauvenet algorithm for the image statistics [C9] Using Chauvenet algorithm for the image statistics [C10] Using Chauvenet algorithm for the image statistics [C11] Using Chauvenet algorithm for the image statistics [C12] Using Chauvenet algorithm for the image statistics [C13] Using Chauvenet algorithm for the image statistics [C14] Using Chauvenet algorithm for the image statistics SidelobeLevel = 0.21961 prune size=9.47266(minbeamfrac=0.3 * beampix=31.5755)
- Masking is performed and then pruned. Only "error" messages are reported per channel (i.e. no regions removed or no regions found).
- The mask is smoothed and then in later cycles grown into lower signal-to-noise regions and pruned again.
- Negative emission is masked if present.
*** Start auto-multithresh processing for Channel 0*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 1*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 2*** No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 3*** No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 4*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 5*** No regions are removed in pruning process. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 6*** No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 7*** No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 8*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 9*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 10*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 11*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 12*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 13*** No regions are found in this plane. No negative region was found by auotmask. *** Start auto-multithresh processing for Channel 14*** No regions are found in this plane. No negative region was found by auotmask. *** Timing summary for whole planes *** Total time to create the initial threshold mask: real 0.02s ( user 0.48s, system 0s) Total time to prune the initial threshold mask: real 0.03s (user 0.66s, system 0.01s) Total time to create the smoothed initial threshold mask: real 0.75s (user 2.13s, system 1.19s) Total time to grow the previous mask: real 0s (user 0s, system 0s) Total time to prune the grow mask: real 0s (user 0s, system 0s) Total time to create the smoothed grow mask: real 0s (user 0s, system 0s)
- Automasking evaluates whether or not to continue to try masking each channel. If no masks were found in a channel or all masks were removed through the pruning process, the channel is no longer considered in the upcoming minor clean cycles.
Stopping masking for chan=0 Stopping masking for chan=1 Stopping masking for chan=2 Stopping masking for chan=3 Stopping masking for chan=4 Stopping masking for chan=6 Stopping masking for chan=8 Stopping masking for chan=9 Stopping masking for chan=10 Stopping masking for chan=11 Stopping masking for chan=12 Stopping masking for chan=13 Stopping masking for chan=14
- A summary of the automasking process is given in the "Automask summary" table. Here, you can see the threshold that was used to create the mask (either "noise" or "sidelobe"), the number of regions found, number of regions pruned, etc. Our automasking results in masks for channels 5 and 7.
========== automask summary ========== chan masking? median RMS peak thresh_type thresh_value N_reg N_pruned N_grow N_grow_pruned N_neg_pix [C0] F -- -- 0.12525 noise -- 0 0 0 0 0 [C1] F -- -- -0.113465 noise -- 0 0 0 0 0 [C2] F -- -- 0.136005 noise -- 1 1 0 0 0 [C3] F -- -- 0.148636 noise -- 1 1 0 0 0 [C4] F -- -- 0.119875 noise -- 0 0 0 0 0 [C5] T 0.000376255 0.0328259 0.190839 noise 0.139886 1 0 0 0 0 [C6] F -- -- 0.152505 noise -- 1 1 0 0 0 [C7] T -3.61604e-05 0.0308196 0.166526 noise 0.130947 2 1 0 0 0 [C8] F -- -- 0.102411 noise -- 0 0 0 0 0 [C9] F -- -- 0.106988 noise -- 0 0 0 0 0 [C10] F -- -- -0.117285 noise -- 0 0 0 0 0 [C11] F -- -- 0.109285 noise -- 0 0 0 0 0 [C12] F -- -- 0.108688 noise -- 0 0 0 0 0 [C13] F -- -- -0.111106 noise -- 0 0 0 0 0 [C14] F -- -- -0.106717 noise -- 0 0 0 0 0 ========== END of automask summary ==========
- Cleaning is performed on the masked image and the progress (peak residual, number of iterations, etc.) is reported per channel. Notice how minor cycles (iterations/iters) are performed on our masked channels 5 and 7 only.
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_base_params] Number of pixels in the clean mask : 362 out of a total of 937500 pixels. [ 0.0386133 % ]
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Minor Cycle controls : {'cycleniter': 100000, 'cyclethreshold': 0.08699999749660492, 'loopgain': 0.10000000149011612, 'nsigma': 0.0, 'thresholdreached': True}
itsNsigma=0
---------------------------------------------------- Run Minor Cycle Iterations ---------------------------------------------
Set Deconvolution Options for [twhya_base_params] : hogbom
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
itsNsigma=0
[twhya_base_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.087, CycleNiter=100000, Gain=0.1
[twhya_base_params:C0] iters=0->0 [0], model=0->0, peakres=0.12525->0.12525, Skipped this plane. Zero mask.
[twhya_base_params:C1] iters=0->0 [0], model=0->0, peakres=0.113465->0.113465, Skipped this plane. Zero mask.
[twhya_base_params:C2] iters=0->0 [0], model=0->0, peakres=0.136005->0.136005, Skipped this plane. Zero mask.
[twhya_base_params:C3] iters=0->0 [0], model=0->0, peakres=0.148636->0.148636, Skipped this plane. Zero mask.
[twhya_base_params:C4] iters=0->0 [0], model=0->0, peakres=0.119875->0.119875, Skipped this plane. Zero mask.
[twhya_base_params:C5] iters=0->16 [16], model=0->0.190617, peakres=0.190839->0.085914, Reached cyclethreshold.
[twhya_base_params:C6] iters=16->16 [0], model=0->0, peakres=0.152505->0.152505, Skipped this plane. Zero mask.
[twhya_base_params:C7] iters=16->39 [23], model=0->0.258554, peakres=0.166526->0.0863384, Reached cyclethreshold.
[twhya_base_params:C8] iters=39->39 [0], model=0->0, peakres=0.102411->0.102411, Skipped this plane. Zero mask.
[twhya_base_params:C9] iters=39->39 [0], model=0->0, peakres=0.106988->0.106988, Skipped this plane. Zero mask.
[twhya_base_params:C10] iters=39->39 [0], model=0->0, peakres=0.117285->0.117285, Skipped this plane. Zero mask.
[twhya_base_params:C11] iters=39->39 [0], model=0->0, peakres=0.109285->0.109285, Skipped this plane. Zero mask.
[twhya_base_params:C12] iters=39->39 [0], model=0->0, peakres=0.108688->0.108688, Skipped this plane. Zero mask.
[twhya_base_params:C13] iters=39->39 [0], model=0->0, peakres=0.111106->0.111106, Skipped this plane. Zero mask.
[twhya_base_params:C14] iters=39->39 [0], model=0->0, peakres=0.106717->0.106717, Skipped this plane. Zero mask.
[twhya_base_params] Total model flux (over all planes) : 0.44917 Peak Residual (over all planes) : 0.0863384 in C7:P0
Completed 39 iterations.
Initial Mask Threshold Calculation
The initial mask threshold is calculated by the following method:
noiseThresholdValue = noisethreshold * rms in residual image
sidelobeThresholdValue = sidelobethreshold * sidelobeLevel * peak in residual image
where the sidelobeLevel is calculated by tclean and reported in the logger (you can see this in the example logger output above).
The initial threshold is max(noiseThresholdValue, sidelobeThresholdValue).
In practice, the sidelobethreshold often sets the limit in the first few major cycles because the sidelobe emission is often larger than the random noise in the image. Once the strongest emission has been removed from the residual image, later cycles tend to use the noisethreshold.
Turning Down noisethreshold
First we start by turning down the value for noisethreshold to something like 1.75.
#In CASA
os.system('rm -rf twhya_low_noise_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_noise_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=1.75, # lowered from base
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)
Open the following in CARTA:
twhya_low_noise_params.image
twhya_low_noise_params.residual
twhya_low_noise_params.mask
twhya_base_params.mask
Using the method described above, display the .image and .residual sibe-by-side, overlay the base mask as black contours, and overlay the lownoise mask as white contours.
By reducing noisethreshold, auto-multithresh is now using sidelobethreshold instead to draw the initial masks for most channels. You can see this in the logger. The "thresh_type" column of the "automask summary" reports which threshold is being used to find maskable emission in each channel. The output below is trimmed for brevity indicated by "*****".
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[twhya_lownoise_params] Setting up an auto-mask within PB mask limit
---------------------------------------------------- Run Automask ---------------------------------------------
Set Deconvolution Options for [twhya_lownoise_params] : hogbom
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
Generating AutoMask
*****
SidelobeLevel = 0.21961
prune size=9.47266(minbeamfrac=0.3 * beampix=31.5755)
*****
========== automask summary ==========
chan masking? median RMS peak thresh_type thresh_value N_reg N_pruned N_grow N_grow_pruned N_neg_pix
[C0] T 0.000326007 0.0306058 0.12525 sidelobe 0.0553385 206 152 0 0 0
[C1] T 0.000123019 0.0282244 -0.113465 sidelobe 0.0499591 206 152 0 0 0
[C2] T 4.33802e-05 0.03036 0.136005 sidelobe 0.0597796 150 105 0 0 0
[C3] T 0.000242864 0.0281953 0.148636 sidelobe 0.065527 95 85 0 0 0
[C4] T 0.000411201 0.0290491 0.119875 sidelobe 0.0530629 210 156 0 0 0
[C5] T 0.000376255 0.0328259 0.190839 sidelobe 0.0841967 43 40 0 0 0
[C6] T 0.000135499 0.0357602 0.152505 sidelobe 0.0671187 126 82 0 0 0
[C7] T -3.61604e-05 0.0308196 0.166526 sidelobe 0.0731053 56 40 0 0 0
[C8] T 0.000578858 0.0269325 0.102411 noise 0.0477108 233 178 0 0 0
[C9] T 5.75308e-05 0.0267095 0.106988 sidelobe 0.0470487 229 155 0 0 0
[C10] T -3.76976e-05 0.0264451 -0.117285 sidelobe 0.0514762 189 154 0 0 0
[C11] T 5.7293e-06 0.0268588 0.109285 sidelobe 0.048006 212 150 0 0 0
[C12] T -0.000308202 0.0278082 0.108688 noise 0.0483562 219 154 0 0 0
[C13] T 1.40047e-05 0.0268257 -0.111106 sidelobe 0.0488139 214 159 0 0 0
[C14] T 0.00029441 0.0257197 -0.106717 sidelobe 0.0471667 220 160 0 0 0
========== END of automask summary ==========
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_lownoise_params] Number of pixels in the clean mask : 109074 out of a total of 937500 pixels. [ 11.6346 % ]
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Minor Cycle controls : {'cycleniter': 100000, 'cyclethreshold': 0.08699999749660492, 'loopgain': 0.10000000149011612, 'nsigma': 0.0, 'thresholdreached': True}
itsNsigma=0
---------------------------------------------------- Run Minor Cycle Iterations ---------------------------------------------
Set Deconvolution Options for [twhya_lownoise_params] : hogbom
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
itsNsigma=0
[twhya_lownoise_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.087, CycleNiter=100000, Gain=0.1
[twhya_lownoise_params:C0] iters=0->30 [30], model=0->0.245732, peakres=0.12525->0.0869206, Reached cyclethreshold.
[twhya_lownoise_params:C1] iters=30->47 [17], model=0->0.0652784, peakres=0.105279->0.0865448, Reached cyclethreshold.
[twhya_lownoise_params:C2] iters=47->90 [43], model=0->0.310039, peakres=0.136005->0.0869305, Reached cyclethreshold.
[twhya_lownoise_params:C3] iters=90->109 [19], model=0->0.194673, peakres=0.148636->0.0853744, Reached cyclethreshold.
[twhya_lownoise_params:C4] iters=109->137 [28], model=0->0.162311, peakres=0.119875->0.0869172, Reached cyclethreshold.
[twhya_lownoise_params:C5] iters=137->160 [23], model=0->0.258899, peakres=0.190839->0.0869799, Reached cyclethreshold.
[twhya_lownoise_params:C6] iters=160->243 [83], model=0->0.851988, peakres=0.152505->0.0869064, Reached cyclethreshold.
[twhya_lownoise_params:C7] iters=243->288 [45], model=0->0.419446, peakres=0.166526->0.0862288, Reached cyclethreshold.
[twhya_lownoise_params:C8] iters=288->297 [9], model=0->0.00968076, peakres=0.102411->0.0869284, Reached cyclethreshold.
[twhya_lownoise_params:C9] iters=297->310 [13], model=0->0.0652534, peakres=0.106988->0.0856935, Reached cyclethreshold.
[twhya_lownoise_params:C10] iters=310->324 [14], model=0->0.0348655, peakres=0.117285->0.0865639, Reached cyclethreshold.
[twhya_lownoise_params:C11] iters=324->345 [21], model=0->0.105682, peakres=0.109285->0.086327, Reached cyclethreshold.
[twhya_lownoise_params:C12] iters=345->357 [12], model=0->0.0756125, peakres=0.108688->0.0868853, Reached cyclethreshold.
[twhya_lownoise_params:C13] iters=357->372 [15], model=0->0.10012, peakres=0.0996836->0.0860267, Reached cyclethreshold.
[twhya_lownoise_params:C14] iters=372->382 [10], model=0->0.0373442, peakres=0.106717->0.0866348, Reached cyclethreshold.
[twhya_lownoise_params] Total model flux (over all planes) : 2.93693 Peak Residual (over all planes) : 0.0869799 in C5:P0
Completed 382 iterations.
Notice that masking and cleaning occurs in all channels. Only channels 8 and 12 use thresh_type = noise. In later major cycles, we can see that more channels switch from sidelobe to noise. The threshold used by auto-multithresh are independent for each channel and can change between major cycles as clean progresses. In Figure 3, we see that fainter emission gets picked up by automasking when we reduce the noisethreshold value (as expected), but that noise spikes are clearly being picked up in both line and line-free channels. We can also see this clearly in the log as the number of masks created and reported in the automask summary is much, much higher than in the base parameters log output.
Turning Down sidelobethreshold
Now we turn sidelobethreshold down to see what impact this has on the masks that are created. Note that because of how masks are drawn during each step of the automasking process, we need to also turn down either noisethreshold or lownoisethreshold to guarantee that the sidelobethreshold is used by the algorithm during either the initial step or the secondary growing of the masks step.
#In CASA
os.system('rm -rf twhya_low_sidelobe_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_sidelobe_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=1.0, # lowered from base
lownoisethreshold=0.5, # lowered from base
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)
In Figure 4 we see that turning down sidelobethreshold makes a mask that covers a much larger region and clearly covers areas of the image that do not contain real emission.
The Effects of Pruning
Pruning is used to remove small noise-like regions from the resulting mask. Note that regions are pruned before being convolved with a Gaussian to create a buffer around the mask.
Turning Off Pruning
Now we use the base parameters but turn off pruning by setting minbeamfrac to 0.0 so that no regions are removed from the mask.
#In CASA
os.system('rm -rf twhya_low_minbeamfrac_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_minbeamfrac_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=0.0, # lowered from base
negativethreshold=15.0,
verbose=True,
fastnoise=False)
In Figure 5 we see that with no pruning, a noise spike near the edge of channel 7 was picked up, kept, and cleaned. Notice that several line-free channels include masks around noise spikes as well.
Increasing The Pruning
Next, we set minbeamfrac to 2.0 and leave the rest of parameters at their base values, so that all masks smaller than twice the beam will get pruned.
#In CASA
os.system('rm -rf twhya_high_minbeamfrac_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_high_minbeamfrac_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=2.0, # raised from base
negativethreshold=15.0,
verbose=True,
fastnoise=False)
Masks were found but then immediately removed, thus the resulting mask is empty. This can be seen by looking at the log messages:
SidelobeLevel = 0.21961 prune size=63.151(minbeamfrac=2 * beampix=31.5755) ***** No regions found by automasking 15 of 15 channels used the noise threshold to mask, but no emission was found. Try decreasing your noisethreshold parameter if you want to capture emission in these channels. 5 of 15 channels had all regions removed by pruning. Try decreasing minbeamfrac to remove fewer regions ***** [twhya_high_minbeamfrac_params] Number of pixels in the clean mask : 0 out of a total of 937500 pixels. [ 0 % ] ***** Reached global stopping criterion : zero mask ***** Restoring with an empty model image. Only residuals will be processed to form the output restored image.
In this case, no cleaning is done, a dirty cube is generated, and tclean is exited.
Note that pruning does not calculate whether the emission within a region is significant, so you can inadvertently prune point sources with high signal-to-noise emission (e.g., calibrators). In this case, reduce the minbeamfrac.
Growing The Mask Into Low Signal-to-Noise Regions
The lownoisethreshold parameter is used to define a low signal-to-noise mask. The mask from the previous cycle is then grown into the region defined by the low signal-to-noise mask via binary dilation. To protect against masking sidelobes with the low signal-to-noise mask, the threshold for this mask is determined using the following heuristic
lowNoiseThresholdValue = lownoisethreshold * rms in residual image
max(sidelobeThresholdValue, lowNoiseThresholdValue).
If the value of the lownoisethreshold is lower than the value of the sidelobethreshold, then the mask will not be expanded into low signal-to-noise regions.
Finally, we lower the lownoisethreshold value to try and grow the mask into regions of lower signal-to-noise during the binary dilation stage.
#In CASA
os.system('rm -rf twhya_low_lownoise_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_lownoise_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.0, # lowered from base
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)
This mask is the same as the base mask even if we try to make the lownoisethreshold very, very small (i.e. something like 0.0001). This is due to the fact that if the value of the lownoisethreshold is less than the value of the sidelobethreshold then the sidelobethreshold is used to set mask. This effectively attempts to act as a “fail-safe” to ensure that overmasking does not occur accidentally by setting one of the thresholds very small. For demonstration purposes, we can force it to use a small lownoisethreshold by also setting the sidelobethreshold very small.
#In CASA
os.system('rm -rf twhya_low_lownoise+sidelobe_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_lownoise+sidelobe_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=0.5, # lowered from base
lownoisethreshold=1.0, # lowered from base
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)
We can see the result of this in Figure 7 which shows a much more extended mask that covers some nearby emission but extends down into lower signal-to-noise regions.
Advanced Use Case - Merging User Masks with Automasking
In some situations, you may find that automasking is not masking obvious emission even after adjusting the thresholds. This most commonly occurs in cases:
- with extremely bright, diffuse emission throughout the entire field of view
- where there is low-level absorption that has a similar depth to any negative bowls in the same channel
- where the cube has not been continuum subtracted.
The first case is due to the difficulty of estimating noise in the presence of wide-spread signal. The last two are cases when the astronomer “knows” where the real signal is (likely through some additional information), but where the algorithm does not have enough knowledge to distinguish.
In these cases, it may be helpful to “kick start” automasking by first defining masks by hand for a round of cleaning and then having automasking pick up from there. To continue cleaning an existing image, we simply use the restart=True option within tclean (which is the default value).
First, we perform a round of cleaning by hand by setting interactive=True so that tclean will let us manually create masks before each major cycle. As mentioned above in the “Make a Dirty Cube” section, in order to get identical cleaning behavior between interactive and non-interactive modes, you must set cycleniter equal to niter. With this arbitrary very large number (cycleniter=niter=100000), the minor cycle iteration limit per major cycle will not be reached. Instead, tclean will continue cleaning until the internally computed cyclethreshold or the user specified threshold is hit. One last parameter to set is usemask='user' to indicate that you would like to manually mask emission. Below is an example tclean command using the same dataset as above.
#In CASA
os.system('rm -rf twhya_manual+automask.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_manual+automask',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=True, # manual masking in interactive GUI
niter=100000,
cycleniter=100000, # set equal to niter to have interactive mode replicate the behavior of non-interactive mode
threshold='87mJy',
usemask='user' # default value, but shown here explicitly in place of auto-multithresh
)
Once you are satisfied that you have masked all the obvious emission and cleaned to a level where automasking will start to pick up extended emission, you can exit tclean by pressing the red “X”. After tclean finishes, you can start tclean again but now in non-interactive mode (interactive=False), with automasking enabled (usemask='auto-multithresh'), and using the appropriate automasking thresholds that match your data (see the “Table of Standard Values” above). Finally, we set restart=True which will make tclean pick up where it left off using the existing .mask file, along with the existing .model, .psf, .image, etc. that share the same base imagename.
Side tip: if you already have a CASA region file from an earlier clean or one you have created, and is has a different name, you can read it into tclean using the mask parameter, mask='basename.mask'. The region file method can also be used non-interactively on it's own.
You may see the following warnings. This is a known issue, and can be ignored for the above command.
----------------------------------------------------------- Run Major Cycle 2 ------------------------------------- 2024-04-02 23:11:40 WARN task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1572) Error : Error in copying internal T/F mask : Error in deleting internal T/F mask : Region mask0 could not be removed 2024-04-02 23:11:40 WARN task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1572)+ Cannot delete the mask (used in another process)
However, this does prevent restart=True from working correctly. There is a known issue with file locks in tclean. To use restart=True, exit CASA at this step and start a new session before running the command below.
#In CASA
# os.system('rm -rf twhya_manual+automask.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_manual+automask',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False,
restart=True # default value, but shown here explicitly
)
As a word of caution, automasking will not remove any regions, it will only add to the current mask. Therefore, during the manual masking stage, mask conservatively. If you create a fluffy, extended mask by hand and then start automasking, it may not be altered by the algorithm.
Version History
The older version of this guide is here: Automasking Guide CASA 5
