Automasking Guide CASA 5

UNDER CONSTRUCTION

Introduction

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 and 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 documentation. 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.

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 using the noise derived from the residual image itself (via the median absolute deviation) and the beam size and sidelobe level supplied by tclean internally. The numerical values are recalculated at the beginning of every minor cycle and may change as the clean progresses.

The five primary auto-multithresh parameters control what thresholds the algorithm uses and how regions are removed from the mask via pruning.

1. 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.
2. Sidelobethreshold (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.
3. Minbeamfrac (type: double) -- sets the minimum size of a region must have 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.
4. Lownoisethreshold (type: double) -- sets the threshold into which the initial mask (which is determined by noisethreshold or sidelobethreshold) is expanded in order to include low signal-to-noise regions in the mask.
5. Negativethreshold (type: double) -- set 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 will create a mask for values less than -5 sigma.

Auto-multithresh has several additional parameters that usually do not need to be changed from their default values:

6. Growiterations (type: int) -- controls the number of iterations that binary dilation performs. 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.
7. Smoothfactor (type: double) -- controls how much to smooth the initial mask. We do not recommend changing this from its default value.
8. 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.

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.

While all cycle 5 pipeline products will use the CASA 5.1 version of auto-multithresh, auto-multithresh has been updated in CASA 5.3 with significant speed improvements as well as changes to several parameters. The growiterations parameter is CASA 5.3 now indicates the maximum number of iterations that binary dilation forms. The binary dilation may halt earlier than the maximum number of iterations if there is no change in the mask. The dogrowprune parameter 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. The minpercentchange parameter 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. The latter parameter is should be used with care and is turned off (minpercentchange=-1) by default. Finally, verbose=True turned on per-channel information on the masks.

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. As for combined data (the last row of the table), these values are listed as tentative and should serve as a starting point for your analysis and could be dependent on the type of emission you have.

Array	sidelobethreshold	noisethreshold	minbeamfrac	lownoisethreshold	negativethreshold
12m (short) b75<300m	2.0	4.25	0.3	1.5	0.0 (continuum)/15.0 (line)
12m (long) b75>300m	3.0	5.0	0.3	1.5	0.0 (continuum)/7.0 (line)
7m (continnum/line)	1.25	5.0	0.1	2.0	0.0
12m + 7m combined TENTATIVE	2.0	4.25	0.3	1.5	0.0

The term “b75” refers to the 75th percentile of baselines. You can find this value in the 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 these instructions: https://casaguides.nrao.edu/index.php/Analysis_Utilities

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 good 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.

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: os.environ['SAVE_ALL_AUTOMASKS']="true". Note that 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 minor cycle is triggered can be different in the two cases.

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 not optimal value and end up masking noise features or not masking real emission.

To demonstrate the use of the automasking parameters, we will be using the TW Hydra dataset which is a small dataset from the “First Look” CASA Guides and is available here: https://bulk.cv.nrao.edu/almadata/sciver/TWHya/sis14_twhya_calibrated_flagged.ms.contsub.tgz

This file is the calibrated, continuum subtracted Band 7 data. Untar the file and you should see a file named sis14_twhya_calibrated_flagged.ms.contsub. Almost all of the values used in our tclean commands remain constant throughout this section and are explained and taken from the First Look at Spectral Line Imaging CASA 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.

#In CASA
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_dirtycube',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0Jy’,
# Automasking Parameters below this line
usemask='auto-multithresh’,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

In the following tutorial, we will run tclean non-interactively, and thus need to set a *clean* threshold to avoid over-cleaning emission. The clean threshold parameter is not related to the auto-multithresh parameters. Using the viewer and selecting a region in a line free channel (e.g., channel 3) we find that the RMS in a line free channel is 1.78e-2 Jy/beam. To set our threshold for cleaning, we set it to twice the RMS so threshold= 2*RMS = 0.0374 Jy. (The "/beam" is implied here.)

You do not need to set a clean threshold for auto-masking to work. The auto-multithresh algorithm works in both interactive and non-interactive mode. You can even modify the masks that auto-multithresh produces in interactive mode. However, the masks produced interactively may differ from the masks produced in non-interactive mode because of differences in how often major cycles are triggered. To get identical behavior, set cycleniter=niter in tclean.

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 and, thus, whether our data is "long" or "short" baseline.

#In CASA
au.getBaselineStats('sis14_twhya_calibrated_flagged.ms.contsub')

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 in the log output) is 197.4m, which puts us solidly in the "short" baseline category. Therefore, here we use the “standard”, base automasking values that the pipeline uses for 12m short baseline data from the table above. This will give us a mask to compare our testing against and will be displayed as black contours in all of the following figures.

#In CASA
os.system('rm -rf twhya_base_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_base_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5, 
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

Initial mask threshold calculation (sidelobethreshold and noisethreshold parameters)

The initial mask threshold is calculated by the following method:

sidelobeThresholdValue = sidelobeThreshold * sidelobeLevel * peak in residual image
noiseThresholdValue = noiseThreshold * rms 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 rms in the residual image is calculated using the median absolute deviation in CASA 5.1, which can yield poor estimates in the case where emission covers a large fraction of the field. The noise calculation has been improved in CASA 5.3.

The initial threshold is max(sidelobeThresholdValue,noiseThresholdValue).

In practice, the sidlobethreshold often sets the limit before noisethreshold because the sidelobe emission is often larger than the random noise in the image.

Turning Down noisethreshold

First we start by turning down the value for noisethreshold to something like 1.5.

<figure id="small_noisethreshold.png">

On the left is the image and on the right is the residual of channel 6 of 15 of the cube created with the base parameters. The black contours represent the mask created when using the base parameters and the red contours represent the masks created using a small noisethreshold value.

</figure>

#In CASA
os.system('rm -rf twhya_lownoise_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_lownoise_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=2.0,
noisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

By reducing noisethreshold, auto-multithresh is now using sidelobetheshold is used instead to draw the initial masks. You can see this in the logger while the automasking is working, tclean reports which threshold is being used to find maskable emission in the image.

----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
Found full images : 1
Found part images : 0
[twhya_lownoise_params] Peak residual (max,min) over full image : (0.149852,-0.0818088)
[twhya_lownoise_params] Total Model Flux : 0
[twhya_lownoise_params] Setting up an auto-mask within PB mask limit
Generating AutoMask
SidelobeLevel = 0.219842
prune size=13.5358(minbeamfrac=0.3 * beampix=45.1195)
rms from MAD (mads*1.4826)= [0.0194964, 0.018565, 0.0193204, 0.0185711, 0.0192612, 0.0224527, 0.0256428, 0.0212457, 0.0176657, 0.0175954, 0.0178157, 0.0176049, 0.0176475, 0.0178578, 0.0171668]
Using sidelobe threshold for chan 0 threshold=0.0308191
Using sidelobe threshold for chan 1 threshold=0.0304624
Using sidelobe threshold for chan 2 threshold=0.0314644
Using sidelobe threshold for chan 3 threshold=0.0307041
Using sidelobe threshold for chan 4 threshold=0.0380167
Using sidelobe threshold for chan 5 threshold=0.0585222
Using sidelobe threshold for chan 6 threshold=0.0658874
Using sidelobe threshold for chan 7 threshold=0.0513938
Using sidelobe threshold for chan 8 threshold=0.0326665
Using sidelobe threshold for chan 9 threshold=0.0289218
Using sidelobe threshold for chan 10 threshold=0.0288309
Using sidelobe threshold for chan 11 threshold=0.028183
Using sidelobe threshold for chan 12 threshold=0.0299507
Using sidelobe threshold for chan 13 threshold=0.0337656
Using sidelobe threshold for chan 14 threshold=0.0296112

In later cycles, we can see that some channels switch to using the noisethreshold. The threshold used by auto-multithresh are independent for each channel and can change between major cycles as clean progresses. In the image below 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.

Turning Down sidelobethreshold

Now we turn sidelobethreshold down to see what impact this has on the masks that are created. Note that because of what we found out above about 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.

<figure id="small_sidelobethreshold.png">

On the left is the image and on the right is the residual of channel 6 of 15 of the cube created with the base parameters. The black contours represent the mask created when using the base parameters and the red contours represent the mask created using a small sidelobethreshold value.

</figure>

#In CASA
os.system('rm -rf twhya_lowsidelobe_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_lowsidelobe_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=1.0,
noisethreshold=4.25,
lownoisethreshold=0.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

In Figure 2 we see that by 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.

Turning down noisethreshold With No Pruning (i.e. minbeamfrac=0.0)

Now we turn down noisethreshold but also turn off pruning by setting minbeamfrac to 0.0 so that no regions are removed from the mask.

<figure id="small_noisethreshold_nopruning.png">

On the left is the image and on the right is the residual of channel 4 of 15 of the cube created with the base parameters. The red contours represent the mask created using a small noisethreshold value and no pruning.

</figure>

#In CASA
os.system('rm -rf twhya_lownoise_noprune_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_lownoise_noprune_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=2.0,
noisethreshold=3.5,
lownoisethreshold=1.5,
minbeamfrac=0.0,
growiterations=75,
negativethreshold=15.0)

In Figure 3 we show a line free channel. We see that using the base parameters, no masks were created but with a small noisethreshold value and no pruning, noise spikes are picked up and masked by the algorithm.

Base Parameters With Over Pruning

Next, we set minbeamfrac to 1.0 so that only masks smaller than the beam will get pruned.

#In CASA
os.system('rm -rf twhya_highminbeamfrac_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_highminbeamfrac_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
minbeamfrac=1.0,
growiterations=75,
negativethreshold=15.0)

Masks were found but then immediately removed and a dirty cube was generated! Whomp.

CASA$ WARN 3 of 15 channels had all regions removed by pruning. Try decreasing minbeamfrac to remove fewer regions. 
CASA$ WARN Restoring with an empty model image. Only residuals will be processed to form the output restored image. 

This is because automasking initially finds regions of pixels with high signal-to-noise which are often not large patches of emission and smaller than the beam size. Thus with minbeamfrac set to 1.0, all the masks are pruned.

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 by binary dilating the mask.

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 lownoisethreshold is lower than the sidelobethreshold, then the mask will not be expanded into low signal-to-noise regions.

Turning Down lownoisethreshold

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_lowlownoise_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_lowlownoise_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=0.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

<figure id="Small_lownoisethreshold.png">

On the left is the image and on the right is the residual of channel 6 of 15 of the cube created with the base parameters. The black contours represent the mask created when using the base parameters, the red contours represent the mask created using a small lownoisethreshold value.

</figure>

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, as we saw above, 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.

os.system('rm -rf twhya_low_lownoise+sidelobe_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_low_lownoise+sidelobe_params',
field='0',
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.08arcsec',
weighting='briggs',
robust=0.5,
restoringbeam='common',
interactive=False,
chanchunks=-1,
niter=100000,
threshold=’0.0356Jy’,
usemask='auto-multithresh’,
sidelobethreshold=0.5,
noisethreshold=4.25,
lownoisethreshold=1.0,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

We can see the result of this in Figure 4 which shows a much more extended mask that covers some nearby emission but extends down into lower signal-to-noise regions.