Automasking Guide CASA 5: Difference between revisions

From CASA Guides
Jump to navigationJump to search
← Older edit
Alipnick (talk | contribs)
80 edits
m Estarr moved page Automasking Guide to Automasking Guide CASA 5: moving this page to an archived version for CASA 5, allowing the main URL to redirect to the newest version
 
(85 intermediate revisions by 4 users not shown)
Line 1: Line 1:
'''UNDER CONSTRUCTION'''
== Introduction ==


== 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. It can be used by setting the ''usemask'' option in {{tclean}} to ’auto-multithresh’.
Starting with CASA 5.1, a new algorithm called ‘auto-multithresh’ has been incorporated into {{tclean}} to automatically mask regions during the cleaning process. 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 [https://casa.nrao.edu/casadocs/casa-5.4.0/synthesis-imaging/masks-for-deconvolution CASA documentation] and [https://iopscience.iop.org/article/10.1088/1538-3873/ab5e14 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.


A full description of the auto-multithresh algorithm is given in the [https://casa.nrao.edu/casadocs/casa-5.1.1/synthesis-imaging/masks-for-deconvolution plone documentation]. Here we briefly describe the fundamental features of the algorithm. A key feature of the auto-multithresh algorithm is that it operates on the residual image at the beginning of every minor cycle, so that the mask updates as the clean progresses. 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.
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 ==
== Parameters ==
Most of the parameters for the auto-multithresh algorithm are defined in terms of fundamental image properties (e.g., beam size, single-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 five primary auto-multithresh parameters control what thresholds the algorithm uses and how regions are removed from the mask via pruning.  
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.
# '''''Noisethreshold''''' (type: double) -- controls 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) -- controls 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.  
There are five primary auto-multithresh parameters that control what thresholds the algorithm uses and how regions are removed from the mask via pruning.  
# '''''Minbeamfrac''''' (type: double) -- controls the minimum size a region must be to be considered real emission in terms of fractions of the beam size. Any masks smaller than this will be removed before the next round of cleaning. Note that this parameter operates on both the initial threshold mask and the low signal-to-noise mask.
# '''''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.  
# '''''Lownoisethreshold''''' (type: double) -- controls the expansion of the initial mask (which is determined by ''noisethreshold'' or ''sidelobethreshold'') into low signal-to-noise regions.
# '''''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.  
# '''''Negativethreshold''''' (type: double) -- controls 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.  
# '''''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.
# '''''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.  
# '''''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 parameters that usually do not need to be changed from their default values:
Auto-multithresh has several additional parameters that usually do not need to be changed from their default values:
<ol type="1" start=6>
<ol type="1" start=6>
<li> '''''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.  
<li> '''''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.
<li> '''''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.
<li> '''''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.
<li> '''''Smoothfactor''''' (type: double) -- controls how much to smooth the initial mask. We do not recommend changing this from its default value.
<li> '''''Smoothfactor''''' (type: double) -- controls how much to smooth the initial mask. We do not recommend changing this from its default value.
<li> '''''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.
<li> '''''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.
<li> '''''Verbose''''' (type: bool) -- turns on/off per-channel logging information on the masks (default = False)
<li> '''''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 (mask). The latter estimation may be more accurate when emission covers a significant portion of the field of view (default=True).
</ol>
</ol>


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.
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.
Here are some 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.
* 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.


== Table of Standard Values ==
== 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.
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 row are tentative and should serve as a starting point for your analysis.


{| class="wikitable"
{| class="wikitable"
Line 75: Line 74:
|}
|}


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 [https://safe.nrao.edu/wiki/bin/view/ALMA/GetBaselineStats au.getBaselineStats]. You can obtain the Analysis Utilites tasks by following these instructions: https://casaguides.nrao.edu/index.php/Analysis_Utilities
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 [https://safe.nrao.edu/wiki/bin/view/ALMA/GetBaselineStats au.getBaselineStats]. You can obtain the Analysis Utilities tasks by following these instructions: https://casaguides.nrao.edu/index.php/Analysis_Utilities


Below, we provide some examples of how changing the input parameters changes the resulting mask.
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.


== Simple Examples of the Tclean Automasking Parameters and "What Can Go Wrong" ==
== Guidelines for optimizing 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 “over-tune” a parameter 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://casaguides.nrao.edu/index.php/First_Look_at_Imaging#Getting_the_Data
* 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 this feature does not work when tclean is run in parallel on cubes.) 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. Set cycleniter=niter in tclean to get identical behavior.


We’ll be using the calibrated, continuum subtracted data, specifically the 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 [https://casaguides.nrao.edu/index.php/First_Look_at_Line_Imaging 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.  
== 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.
 
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 [https://casaguides.nrao.edu/index.php/First_Look_at_Line_Imaging 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 ===
=== 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.  
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.
 
<source lang="python">
#In CASA
listobs('sis14_twhya_calibrated_flagged.ms.contsub')
</source>
 
<pre style="background-color: #fffacd;">
Observer: cqi    Project: uid://A002/X327408/X6f 
Observation: ALMA
Data records: 53161      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      0 TW Hya                    8514  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:08:09.6 - 08:13:57.3    16      0 TW Hya                  10360  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:19:53.9 - 08:25:41.7    20      0 TW Hya                  10321  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:32:00.5 - 08:37:48.2    24      0 TW Hya                  10324  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:43:45.6 - 08:49:33.4    28      0 TW Hya                    9462  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              09:05:15.6 - 09:07:31.6    36      0 TW Hya                    4180  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
          (nRows = Total number of rows per scan)
Fields: 1
  ID  Code Name                RA              Decl          Epoch  SrcId      nRows
  0    none TW Hya              11:01:51.796000 -34.42.17.36600 J2000  0          53161
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
</pre>
 
With the information above and from what we know of TW Hydra from the [https://casaguides.nrao.edu/index.php/First_Look_at_Line_Imaging First Look at Spectral Line Imaging] CASA Guide, we get the basic information for our tclean call to create the dirty cube.  


<source lang="python">
<source lang="python">
Line 110: Line 147:
interactive=False,
interactive=False,
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=0,
threshold=’0.0Jy’,
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)
</source>
</source>


Using the viewer and selecting a region in a line free channel (I used 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''.
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 ===  
=== Base Parameters ===  
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.  
To figure out what automasking parameter values we should use for this data, we run [https://safe.nrao.edu/wiki/bin/view/ALMA/GetBaselineStats au.getBaselineStats] which will tell us the "b75" value for our data and, thus, whether our data is "long" or "short" baseline.
 
<source lang="python">
#In CASA
au.getBaselineStats('sis14_twhya_calibrated_flagged.ms.contsub')
</source>
<pre style="background-color: #fffacd;">
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
</pre>
 
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.  
 
 
[[Image:Base_parameters.png|600px|thumb|<caption>The black contours mark the final mask created using the base, 12m short baseline automasking parameters. 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.</caption>]]
 


<source lang="python">
<source lang="python">
Line 150: Line 198:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
# Automasking Parameters below this line
usemask='auto-multithresh',
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
noisethreshold=4.25,
Line 157: Line 206:
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
growiterations=75,
negativethreshold=15.0)
negativethreshold=15.0,
verbose=True)
</source>
</source>


=== Turning Down ''noisethreshold'' ===
With ''verbose=True'' set, automasking via auto-multithresh 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. Masking is then performed and then pruned, only "error" messages are reported per channel (i.e. no regions removed or no regions found). The mask is then smoothed and then in later cycles grown into lower signal-to-noise regions and pruned again. Finally, negative emission is masked if present. After this, 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 future clean cycles. After this, 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. Finally, cleaning is performed on the masked image and the progress (peak residual, number of iterations, etc.) is reported per channel.
 
<pre style="background-color: #fffacd;">
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
[twhya_base_params] Peak residual (max,min) over full image : (0.149898,-0.081815)
[twhya_base_params] Total Model Flux : 0
[twhya_base_params] Setting up an auto-mask
Generating AutoMask
SidelobeLevel = 0.219791
prune size=13.535(minbeamfrac=0.3 * beampix=45.1166)
Start thresholding: create an initial mask by threshold
End thresholding: time to create the initial threshold mask:  real 0.37s ( user 0.02s, system 0.02s)
Start pruning: the initial threshold mask
[C0] No regions are found in this plane.
[C1] No regions are found in this plane.
[C2] No regions are found in this plane.
[C3] No regions are found in this plane.
[C7] No regions are removed in pruning process.
[C8] No regions are found in this plane.
[C9] No regions are found in this plane.
[C10] No regions are found in this plane.
[C11] No regions are found in this plane.
[C12] No regions are found in this plane.
[C14] No regions are found in this plane.
End pruning: time to prune the initial threshold mask: real 0.24s (user 0.05s, system 0.02s)
Start smoothing: the initial threshold mask
End smoothing: time to create the smoothed initial threshold mask: real 2.25s (user 2.92s, system 1.39s)
Start thresholding: create a negative mask
No negative region was found by auotmask.
End thresholding: time to create the negative mask: real 1.19s (user 0.54s, system 0.26s)
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=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
========== 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        5.1095e-05  0.0195051  -0.081815  noise  0.0828968        0  0  --  --  0
[C1] F        -6.20534e-05  0.0185451  -0.0749006  noise  0.0788167        0  0  --  --  0
[C2] F        -2.53807e-06  0.0193195  0.0714881  noise  0.0821077        0  0  --  --  0
[C3] F        7.66067e-05  0.0185695  -0.0724906  noise  0.0789204        0  0  --  --  0
[C4] F        -0.000102242  0.0192787  0.0862985  noise  0.0819347        2  2  --  --  0
[C5] T        -0.000765545  0.0224421  0.133081  noise  0.0953791        2  1  --  --  0
[C6] T        -0.00150215  0.0256315  0.149898  noise  0.108934        2  1  --  --  0
[C7] T        -0.000590894  0.0212334  0.116813  noise  0.0902421        1  0  --  --  0
[C8] F        1.70437e-05  0.0176792  -0.0750749  noise  0.0751365        0  0  --  --  0
[C9] F        0.000338539  0.0175967  0.0657392  noise  0.0747862        0  0  --  --  0
[C10] F        -3.07989e-05  0.0178246  -0.0780685  noise  0.0757547        0  0  --  --  0
[C11] F        -0.000249928  0.0176007  -0.0744855  noise  0.074803        0  0  --  --  0
[C12] F        0.000142292  0.0176476  -0.0685521  noise  0.0750021        0  0  --  --  0
[C13] F        -0.000252258  0.0178448  0.0768286  noise  0.0758403        1  1  --  --  0
[C14] F        0.000200417  0.0171804  0.0673071  noise  0.0730168        0  0  --  --  0
========== END of automask summary ==========
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_base_params] Number of pixels in the clean mask : 1209 out of a total of 937500 pixels. [ 0.12896 % ]
[twhya_base_params] Peak residual (max,min) within mask : (0.149898,-0.0296383) over full image : (0.149898,-0.081815)
[twhya_base_params] Total Model Flux : 0
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
[twhya_base_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.0356, CycleNiter=100000, Gain=0.1
[twhya_base_params:C0] iters=0->0 [0], model=0->0, peakres=0.081815->0.081815, Skipped this plane. Zero mask.
[twhya_base_params:C1] iters=0->0 [0], model=0->0, peakres=0.0749006->0.0749006, Skipped this plane. Zero mask.
[twhya_base_params:C2] iters=0->0 [0], model=0->0, peakres=0.0714881->0.0714881, Skipped this plane. Zero mask.
[twhya_base_params:C3] iters=0->0 [0], model=0->0, peakres=0.0724906->0.0724906, Skipped this plane. Zero mask.
[twhya_base_params:C4] iters=0->0 [0], model=0->0, peakres=0.0862985->0.0862985, Skipped this plane. Zero mask.
[twhya_base_params:C5] iters=0->85 [85], model=0->0.513613, peakres=0.133081->0.0350699, Reached cyclethreshold.
[twhya_base_params:C6] iters=85->177 [92], model=0->0.564421, peakres=0.149898->0.0355323, Reached cyclethreshold.
[twhya_base_params:C7] iters=177->256 [79], model=0->0.450467, peakres=0.116813->0.0349898, Reached cyclethreshold.
[twhya_base_params:C8] iters=256->256 [0], model=0->0, peakres=0.0750749->0.0750749, Skipped this plane. Zero mask.
[twhya_base_params:C9] iters=256->256 [0], model=0->0, peakres=0.0657392->0.0657392, Skipped this plane. Zero mask.
[twhya_base_params:C10] iters=256->256 [0], model=0->0, peakres=0.0780685->0.0780685, Skipped this plane. Zero mask.
[twhya_base_params:C11] iters=256->256 [0], model=0->0, peakres=0.0744855->0.0744855, Skipped this plane. Zero mask.
[twhya_base_params:C12] iters=256->256 [0], model=0->0, peakres=0.0685521->0.0685521, Skipped this plane. Zero mask.
[twhya_base_params:C13] iters=256->256 [0], model=0->0, peakres=0.0768286->0.0768286, Skipped this plane. Zero mask.
[twhya_base_params:C14] iters=256->256 [0], model=0->0, peakres=0.0673071->0.0673071, Skipped this plane. Zero mask.
[twhya_base_params] Total model flux (over all planes) : 1.5285    Peak Residual (over all planes) : 0.0862985 in C4:P0
Completed 256 iterations.
</pre>
 
=== Initial mask threshold calculation (sidelobethreshold and noisethreshold parameters)===
 
The initial mask threshold is calculated by the following method:
 
<blockquote>
sidelobeThresholdValue = sidelobeThreshold * sidelobeLevel * peak in residual image <br>
noiseThresholdValue = noiseThreshold * rms in residual image
</blockquote>
where the sidelobeLevel is calculated by {{tclean}} and reported in the logger (you can see this in the example logger output below).  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 will been improved in a future version of CASA.
 
The initial threshold is max(sidelobeThresholdValue,noiseThresholdValue).
 
In practice, the ''sidlobethreshold'' often sets the limit in the first few major cycle 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.5.   
First we start by turning down the value for noisethreshold to something like 1.5.   


<figure id="small_noisethreshold.png">
 
[[Image:small_noisethreshold.png|700px|thumb|<caption>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.</caption>]]
[[Image:small_noisethreshold.png|600px|thumb|<caption>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.</caption>]]
</figure>
 


<source lang="python">
<source lang="python">
Line 190: Line 339:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
usemask='auto-multithresh',
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=1.5,
noisethreshold=1.5,
Line 197: Line 346:
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
growiterations=75,
negativethreshold=15.0)
negativethreshold=15.0,
verbose=True)
</source>
</source>


What we see is that since ''noisethreshold'' is too low, ''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.  
By reducing ''noisethreshold'', auto-multithresh is now using ''sidelobetheshold'' 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 in the "automask summary" in the "thresh_type" column.  


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
Found full images : 1
[twhya_lownoise_params] Peak residual (max,min) over full image : (0.149898,-0.081815)
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] Total Model Flux : 0
[twhya_lownoise_params] Setting up an auto-mask within PB mask limit
[twhya_lownoise_params] Setting up an auto-mask
Generating AutoMask
Generating AutoMask
SidelobeLevel = 0.219842
SidelobeLevel = 0.219791
prune size=13.5358(minbeamfrac=0.3 * beampix=45.1195)
prune size=13.535(minbeamfrac=0.3 * beampix=45.1166)
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]
Start thresholding: create an initial mask by threshold
Using sidelobe threshold for chan 0 threshold=0.0308191
End thresholding: time to create the initial threshold mask:  real 0.22s ( user 0.03s, system 0.01s)
Using sidelobe threshold for chan 1 threshold=0.0304624
Start pruning: the initial threshold mask
Using sidelobe threshold for chan 2 threshold=0.0314644
End pruning: time to prune the initial threshold mask: real 2.27s (user 0.96s, system 0.04s)
Using sidelobe threshold for chan 3 threshold=0.0307041
Start smoothing: the initial threshold mask
Using sidelobe threshold for chan 4 threshold=0.0380167
End smoothing: time to create the smoothed initial threshold mask: real 5.69s (user 10.47s, system 6.63s)
Using sidelobe threshold for chan 5 threshold=0.0585222
Start thresholding: create a negative mask
Using sidelobe threshold for chan 6 threshold=0.0658874
No negative region was found by auotmask.
Using sidelobe threshold for chan 7 threshold=0.0513938
End thresholding: time to create the negative mask: real 1.23s (user 1.15s, system 0.18s)
Using sidelobe threshold for chan 8 threshold=0.0326665
========== automask summary ==========
Using sidelobe threshold for chan 9 threshold=0.0289218
chan masking? median  RMS        peak  thresh_type  thresh_value  N_reg N_pruned N_grow N_grow_pruned N_neg_pix
Using sidelobe threshold for chan 10 threshold=0.0288309
[C0] T        5.1095e-05  0.0195051  -0.081815  sidelobe  0.030867        271  168  --  --  0
Using sidelobe threshold for chan 11 threshold=0.028183
[C1] T        -6.20534e-05  0.0185451  -0.0749006  sidelobe  0.0304619        254  153  --  --  0
Using sidelobe threshold for chan 12 threshold=0.0299507
[C2] T        -2.53807e-06  0.0193195  0.0714881  sidelobe  0.0314249        261  170  --  --  0
Using sidelobe threshold for chan 13 threshold=0.0337656
[C3] T        7.66067e-05  0.0185695  -0.0724906  sidelobe 0.0306932        267  171  --  --  0
Using sidelobe threshold for chan 14 threshold=0.0296112
[C4] T        -0.000102242  0.0192787  0.0862985  sidelobe 0.0379353        172  136  --  --  0
[C5] T        -0.000765545  0.0224421  0.133081  sidelobe 0.0585003        28  25  --  --  0
[C6] T        -0.00150215  0.0256315  0.149898  sidelobe 0.0658923        27  19  --  --  0
[C7] T        -0.000590894  0.0212334  0.116813  sidelobe 0.0513489        44  37  --  --  0
[C8] T        1.70437e-05  0.0176792  -0.0750749  sidelobe 0.0326501        204  151  --  --  0
[C9] T        0.000338539  0.0175967  0.0657392  sidelobe 0.0288978        266  177  --  --  0
[C10] T        -3.07989e-05  0.0178246  -0.0780685  sidelobe 0.0288354        278  186  --  --  0
[C11] T        -0.000249928  0.0176007  -0.0744855  sidelobe 0.0281572        281  184  --  --  0
[C12] T        0.000142292  0.0176476  -0.0685521  sidelobe 0.029941        225  155  --  --  0
[C13] T        -0.000252258  0.0178448  0.0768286  sidelobe 0.0337725        184  146  --  --  0
[C14] T        0.000200417  0.0171804  0.0673071  sidelobe 0.029587        250  182  --  --  0
========== END of automask summary ==========
</pre>
</pre>


In later cycles, we can see that some channels switch to using the ''noisethreshold''. Which threshold {{tclean}} uses can differ per 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 clearly noise spikes are being picked up in both line and line-free channels.  
In later clean 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 Figure 2, 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.  
 
 
As well as choosing between ''noisethreshold'' and ''sidelobethreshold'', automasking additionally chooses between ''lownoisethreshold'' and ''sidelobethreshold'' during the second round of masking when the mask is grown into low signal-to-noise regions. The actual logic and values that are used for each threshold are the following.
 
<blockquote>
lowNoiseThresholdValue = lowNoiseThreshold * rms in residual image <br>
sidelobeThresholdValue = sidelobeThreshold * sidelobeLevel * peak in residual image <br>
noiseThresholdValue = noiseThreshold * rms in residual image
</blockquote>


Where the sidelobeLevel is calculated by {{tclean}} and reported in the logger (you can see this in the example logger output above). The value that is used in the first round of masking is
==== 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.


<blockquote>
max(sidelobeThresholdValue,noiseThresholdValue).
</blockquote>


Then, for the second round of masking where the mask is expanded into low signal-to-noise regions the threshold that is used is determined by
[[Image:small_sidelobethreshold.png|600px|thumb|<caption>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.</caption>]]


<blockquote>
max(sidelobeThresholdValue,lowNoiseThresholdValue).
</blockquote>
In practice, the ''sidlobethreshold'' often sets the limit before ''noisethreshold'' because the sidelobe emission is often larger than the random noise in the image. Also, if the ''sidelobethreshold'' is used to set the initial mask, then the mask will not be grown at all because it will also be larger than the ''lownoisetheshold''.
=== 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">
[[Image:small_sidelobethreshold.png|700px|thumb|<caption>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.</caption>]]
</figure>


<source lang="python">
<source lang="python">
Line 285: Line 421:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
usemask='auto-multithresh',
sidelobethreshold=1.0,
sidelobethreshold=1.0,
noisethreshold=4.25,
noisethreshold=4.25,
Line 295: Line 431:
</source>
</source>


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.
In Figure 3 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 been convolved with a Gaussian to create a buffer around the mask.
 
==== Turning off pruning ====


=== Turning down ''noisethreshold'' With No Pruning (i.e. ''minbeamfrac=0.0'') ===
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.  
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.  
 
 
[[Image:small_noisethreshold_nopruning.png|600px|thumb|<caption> 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.</caption>]]


<figure id="small_noisethreshold_nopruning.png">
[[Image:small_noisethreshold_nopruning.png|700px|thumb|<caption> 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.</caption>]]
</figure>


<source lang="python">
<source lang="python">
Line 327: Line 468:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
usemask='auto-multithresh',
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=3.5,
noisethreshold=4.25,
lownoisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=0.0,
minbeamfrac=0.0,
Line 337: Line 478:
</source>
</source>


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.  
In Figure 4 we show a line free channel. We see that using the base parameters, no masks were created but with no pruning, noise spikes are picked up, kept, and then cleaned by the algorithm.


=== Base Parameters With Over Pruning ===
==== Increasing the Pruning ====
Next, we set ''minbeamfrac'' to 1.0 so that only masks smaller than the beam will get pruned.
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.


<source lang="python">
<source lang="python">
Line 365: Line 506:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
usemask='auto-multithresh',
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
noisethreshold=4.25,
lownoisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=1.0,
minbeamfrac=2.0,
growiterations=75,
growiterations=75,
negativethreshold=15.0)
negativethreshold=15.0)
</source>
</source>


Masks were found but then immediately removed and a dirty cube was generated! Whomp.  
Masks were found but then immediately removed, thus the resulting mask is empty. This can be seen by looking at the log tclean messages.


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
CASA$ WARN 3 of 15 channels had all regions removed by pruning. Try decreasing minbeamfrac to remove fewer regions.  
CASA$ WARN 5 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.  
CASA$ WARN Restoring with an empty model image. Only residuals will be processed to form the output restored image.  
</pre>
</pre>


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.  
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
<blockquote>
lowNoiseThresholdValue = lowNoiseThreshold * rms in residual image <br>
max(sidelobeThresholdValue,lowNoiseThresholdValue).
</blockquote>
 
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.


=== 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.  
Finally, we lower the ''lownoisethreshold'' value to try and grow the mask into regions of lower signal-to-noise during the binary dilation stage.  


Line 410: Line 562:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
usemask='auto-multithresh',
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
noisethreshold=4.25,
Line 420: Line 572:
</source>
</source>


Hmm... this offered absolutely no change compared to the base values 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 ''lownoisethreshold'' is set lower than ''sidelobethreshold'', then ''sidelobethreshold'' becomes the lowest cap. Therefore, if the ''sidelobeThresholdValue'' is used to set the mask and is larger than ''lowNoiseThresholdValue'', then the mask is not grown at all. 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. However, for demonstration purposes, we can force it to use a small ''lownoisethreshold'' by also setting the ''sidelobethreshold'' very small.  
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.  
 


<figure id="Small_lownoisethreshold.png">
[[Image:Small_lownoisethreshold.png|600px|thumb|<caption> 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.</caption>]]
[[Image:Small_lownoisethreshold.png|700px|thumb|<caption> 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.</caption>]]
</figure>


<source lang="python">
<source lang="python">
Line 448: Line 599:
chanchunks=-1,
chanchunks=-1,
niter=100000,
niter=100000,
threshold=’0.0356Jy’,
threshold='0.0374Jy',
usemask='auto-multithresh’,
usemask='auto-multithresh',
sidelobethreshold=0.5,
sidelobethreshold=0.5,
noisethreshold=4.25,
noisethreshold=4.25,
Line 458: Line 609:
</source>
</source>


We can see the result of this in Figure 4 which shows a much more extended mask that covers some nearby emission but extends much further into noise structures.
We can see the result of this in Figure 5 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, and cases 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 happen 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 do this, we simply use the ''restart=True'' option within {{tclean}}.
 
First, we perform a round of cleaning by hand by setting ''interactive=True'' so that {{tclean}} will let us manually create masks before each minor 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'' to be equal to ''niter''. By doing this, you are essentially setting ''cycleniter'' to such a large number that the iteration limit will not be reached during the minor cycle, instead {{tclean}} will continue with 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.
 
<source lang="python">
os.system('rm -rf twhya_advanced_use.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_advanced_use',
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=True,
chanchunks=-1,
niter=100000,
cycleniter=100000,
threshold='0.0374Jy',
usemask='user’)
</source>
 
As a time saving tip, if you already have a CASA region file (from an earlier clean or one you have created), you can point {{tclean}} to it using the mask parameter, ''mask=’filename.mask’'', and {{tclean}} will load your masks automatically. This can be useful if you know that some manual cleaning is needed before automasking will start picking up extended emission but several iterations may be needed until you discover that limit. The region file method can also be used non-interactively on it’s own.
 
Once you are satisfied that you have masked all the obvious emision 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’'', 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.
 
<source lang="python">
os.system('rm -rf twhya_advanced_use.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_advanced_use',
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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
restart=True)
</source>
 
As a word of caution, automasking will not remove any regions, it will only add to the current mask. Therefore, if you create a fluffy, extended mask by hand and then start automasking, it may not be altered by the algorithm. Therefore, during the manual masking stage, mask conservatively.

Latest revision as of 20:32, 8 March 2024

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. 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 documentation 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 parameters that 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 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.
  3. 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.
  4. 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.
  5. 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 parameters that usually do not need to be changed from their default values:

  1. 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.
  2. 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.
  3. 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.
  4. Smoothfactor (type: double) -- controls how much to smooth the initial mask. We do not recommend changing this from its default value.
  5. 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.
  6. Verbose (type: bool) -- turns on/off per-channel logging information on the masks (default = False)
  7. 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 (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 row are tentative and should serve as a starting point for your analysis.

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 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 this feature does not work when tclean is run in parallel on cubes.) 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. Set cycleniter=niter in tclean 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.

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. It's always a good idea to perform a listobs first to get the basic information of the dataset. 

#In CASA
listobs('sis14_twhya_calibrated_flagged.ms.contsub')

Observer: cqi     Project: uid://A002/X327408/X6f  
Observation: ALMA
Data records: 53161       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      0 TW Hya                    8514  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:08:09.6 - 08:13:57.3    16      0 TW Hya                   10360  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:19:53.9 - 08:25:41.7    20      0 TW Hya                   10321  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:32:00.5 - 08:37:48.2    24      0 TW Hya                   10324  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:43:45.6 - 08:49:33.4    28      0 TW Hya                    9462  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              09:05:15.6 - 09:07:31.6    36      0 TW Hya                    4180  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
           (nRows = Total number of rows per scan)
Fields: 1
  ID   Code Name                RA               Decl           Epoch   SrcId      nRows
  0    none TW Hya              11:01:51.796000 -34.42.17.36600 J2000   0          53161
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 Spectral Line Imaging CASA Guide, we get the basic information for our tclean call to create the dirty cube. 

#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=0,
threshold='0.0Jy')

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.


The black contours mark the final mask created using the base, 12m short baseline automasking parameters. 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.


#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.0374Jy',
# 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,
verbose=True)

With verbose=True set, automasking via auto-multithresh 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. Masking is then performed and then pruned, only "error" messages are reported per channel (i.e. no regions removed or no regions found). The mask is then smoothed and then in later cycles grown into lower signal-to-noise regions and pruned again. Finally, negative emission is masked if present. After this, 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 future clean cycles. After this, 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. Finally, cleaning is performed on the masked image and the progress (peak residual, number of iterations, etc.) is reported per channel. 

----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
[twhya_base_params] Peak residual (max,min) over full image : (0.149898,-0.081815)
[twhya_base_params] Total Model Flux : 0
[twhya_base_params] Setting up an auto-mask
Generating AutoMask
SidelobeLevel = 0.219791
prune size=13.535(minbeamfrac=0.3 * beampix=45.1166)
Start thresholding: create an initial mask by threshold
End thresholding: time to create the initial threshold mask:  real 0.37s ( user 0.02s, system 0.02s)
Start pruning: the initial threshold mask
[C0] No regions are found in this plane.
[C1] No regions are found in this plane.
[C2] No regions are found in this plane.
[C3] No regions are found in this plane.
[C7] No regions are removed in pruning process.
[C8] No regions are found in this plane.
[C9] No regions are found in this plane.
[C10] No regions are found in this plane.
[C11] No regions are found in this plane.
[C12] No regions are found in this plane.
[C14] No regions are found in this plane.
End pruning: time to prune the initial threshold mask: real 0.24s (user 0.05s, system 0.02s)
Start smoothing: the initial threshold mask
End smoothing: time to create the smoothed initial threshold mask: real 2.25s (user 2.92s, system 1.39s)
Start thresholding: create a negative mask
No negative region was found by auotmask.
End thresholding: time to create the negative mask: real 1.19s (user 0.54s, system 0.26s)
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=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
========== 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        5.1095e-05  0.0195051  -0.081815  noise  0.0828968        0  0  --  --  0
[C1] F        -6.20534e-05  0.0185451  -0.0749006  noise  0.0788167        0  0  --  --  0
[C2] F        -2.53807e-06  0.0193195  0.0714881  noise  0.0821077        0  0  --  --  0
[C3] F        7.66067e-05  0.0185695  -0.0724906  noise  0.0789204        0  0  --  --  0
[C4] F        -0.000102242  0.0192787  0.0862985  noise  0.0819347        2  2  --  --  0
[C5] T        -0.000765545  0.0224421  0.133081  noise  0.0953791        2  1  --  --  0
[C6] T        -0.00150215  0.0256315  0.149898  noise  0.108934        2  1  --  --  0
[C7] T        -0.000590894  0.0212334  0.116813  noise  0.0902421        1  0  --  --  0
[C8] F        1.70437e-05  0.0176792  -0.0750749  noise  0.0751365        0  0  --  --  0
[C9] F        0.000338539  0.0175967  0.0657392  noise  0.0747862        0  0  --  --  0
[C10] F        -3.07989e-05  0.0178246  -0.0780685  noise  0.0757547        0  0  --  --  0
[C11] F        -0.000249928  0.0176007  -0.0744855  noise  0.074803        0  0  --  --  0
[C12] F        0.000142292  0.0176476  -0.0685521  noise  0.0750021        0  0  --  --  0
[C13] F        -0.000252258  0.0178448  0.0768286  noise  0.0758403        1  1  --  --  0
[C14] F        0.000200417  0.0171804  0.0673071  noise  0.0730168        0  0  --  --  0
========== END of automask summary ==========
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_base_params] Number of pixels in the clean mask : 1209 out of a total of 937500 pixels. [ 0.12896 % ]
[twhya_base_params] Peak residual (max,min) within mask : (0.149898,-0.0296383) over full image : (0.149898,-0.081815)
[twhya_base_params] Total Model Flux : 0
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
[twhya_base_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.0356, CycleNiter=100000, Gain=0.1
[twhya_base_params:C0] iters=0->0 [0], model=0->0, peakres=0.081815->0.081815, Skipped this plane. Zero mask.
[twhya_base_params:C1] iters=0->0 [0], model=0->0, peakres=0.0749006->0.0749006, Skipped this plane. Zero mask.
[twhya_base_params:C2] iters=0->0 [0], model=0->0, peakres=0.0714881->0.0714881, Skipped this plane. Zero mask.
[twhya_base_params:C3] iters=0->0 [0], model=0->0, peakres=0.0724906->0.0724906, Skipped this plane. Zero mask.
[twhya_base_params:C4] iters=0->0 [0], model=0->0, peakres=0.0862985->0.0862985, Skipped this plane. Zero mask.
[twhya_base_params:C5] iters=0->85 [85], model=0->0.513613, peakres=0.133081->0.0350699, Reached cyclethreshold.
[twhya_base_params:C6] iters=85->177 [92], model=0->0.564421, peakres=0.149898->0.0355323, Reached cyclethreshold.
[twhya_base_params:C7] iters=177->256 [79], model=0->0.450467, peakres=0.116813->0.0349898, Reached cyclethreshold.
[twhya_base_params:C8] iters=256->256 [0], model=0->0, peakres=0.0750749->0.0750749, Skipped this plane. Zero mask.
[twhya_base_params:C9] iters=256->256 [0], model=0->0, peakres=0.0657392->0.0657392, Skipped this plane. Zero mask.
[twhya_base_params:C10] iters=256->256 [0], model=0->0, peakres=0.0780685->0.0780685, Skipped this plane. Zero mask.
[twhya_base_params:C11] iters=256->256 [0], model=0->0, peakres=0.0744855->0.0744855, Skipped this plane. Zero mask.
[twhya_base_params:C12] iters=256->256 [0], model=0->0, peakres=0.0685521->0.0685521, Skipped this plane. Zero mask.
[twhya_base_params:C13] iters=256->256 [0], model=0->0, peakres=0.0768286->0.0768286, Skipped this plane. Zero mask.
[twhya_base_params:C14] iters=256->256 [0], model=0->0, peakres=0.0673071->0.0673071, Skipped this plane. Zero mask.
[twhya_base_params] Total model flux (over all planes) : 1.5285     Peak Residual (over all planes) : 0.0862985 in C4:P0
Completed 256 iterations.

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 below). 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 will been improved in a future version of CASA.

The initial threshold is max(sidelobeThresholdValue,noiseThresholdValue).

In practice, the sidlobethreshold often sets the limit in the first few major cycle 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.5.


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.


#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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=2.0,
noisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
verbose=True)

By reducing noisethreshold, auto-multithresh is now using sidelobetheshold 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 in the "automask summary" in the "thresh_type" column. 

----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
[twhya_lownoise_params] Peak residual (max,min) over full image : (0.149898,-0.081815)
[twhya_lownoise_params] Total Model Flux : 0
[twhya_lownoise_params] Setting up an auto-mask
Generating AutoMask
SidelobeLevel = 0.219791
prune size=13.535(minbeamfrac=0.3 * beampix=45.1166)
Start thresholding: create an initial mask by threshold
End thresholding: time to create the initial threshold mask:  real 0.22s ( user 0.03s, system 0.01s)
Start pruning: the initial threshold mask
End pruning: time to prune the initial threshold mask: real 2.27s (user 0.96s, system 0.04s)
Start smoothing: the initial threshold mask
End smoothing: time to create the smoothed initial threshold mask: real 5.69s (user 10.47s, system 6.63s)
Start thresholding: create a negative mask
No negative region was found by auotmask.
End thresholding: time to create the negative mask: real 1.23s (user 1.15s, system 0.18s)
========== 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        5.1095e-05  0.0195051  -0.081815  sidelobe  0.030867        271  168  --  --  0
[C1] T        -6.20534e-05  0.0185451  -0.0749006  sidelobe  0.0304619        254  153  --  --  0
[C2] T        -2.53807e-06  0.0193195  0.0714881  sidelobe  0.0314249        261  170  --  --  0
[C3] T        7.66067e-05  0.0185695  -0.0724906  sidelobe  0.0306932        267  171  --  --  0
[C4] T        -0.000102242  0.0192787  0.0862985  sidelobe  0.0379353        172  136  --  --  0
[C5] T        -0.000765545  0.0224421  0.133081  sidelobe  0.0585003        28  25  --  --  0
[C6] T        -0.00150215  0.0256315  0.149898  sidelobe  0.0658923        27  19  --  --  0
[C7] T        -0.000590894  0.0212334  0.116813  sidelobe  0.0513489        44  37  --  --  0
[C8] T        1.70437e-05  0.0176792  -0.0750749  sidelobe  0.0326501        204  151  --  --  0
[C9] T        0.000338539  0.0175967  0.0657392  sidelobe  0.0288978        266  177  --  --  0
[C10] T        -3.07989e-05  0.0178246  -0.0780685  sidelobe  0.0288354        278  186  --  --  0
[C11] T        -0.000249928  0.0176007  -0.0744855  sidelobe  0.0281572        281  184  --  --  0
[C12] T        0.000142292  0.0176476  -0.0685521  sidelobe  0.029941        225  155  --  --  0
[C13] T        -0.000252258  0.0178448  0.0768286  sidelobe  0.0337725        184  146  --  --  0
[C14] T        0.000200417  0.0171804  0.0673071  sidelobe  0.029587        250  182  --  --  0
========== END of automask summary ==========

In later clean 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 Figure 2, 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.


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.


#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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=1.0,
noisethreshold=4.25,
lownoisethreshold=0.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

In Figure 3 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 been 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.


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.


#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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
minbeamfrac=0.0,
growiterations=75,
negativethreshold=15.0)

In Figure 4 we show a line free channel. We see that using the base parameters, no masks were created but with no pruning, noise spikes are picked up, kept, and then cleaned by the algorithm.

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_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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
minbeamfrac=2.0,
growiterations=75,
negativethreshold=15.0)

Masks were found but then immediately removed, thus the resulting mask is empty. This can be seen by looking at the log tclean messages. 

CASA$ WARN 5 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.

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_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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=0.5,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0)

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.


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.
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.0374Jy',
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 5 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, and cases 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 happen 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 do this, we simply use the restart=True option within tclean.

First, we perform a round of cleaning by hand by setting interactive=True so that tclean will let us manually create masks before each minor 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 to be equal to niter. By doing this, you are essentially setting cycleniter to such a large number that the iteration limit will not be reached during the minor cycle, instead tclean will continue with 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. 

os.system('rm -rf twhya_advanced_use.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_advanced_use',
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=True,
chanchunks=-1,
niter=100000,
cycleniter=100000,
threshold='0.0374Jy',
usemask='user’)

As a time saving tip, if you already have a CASA region file (from an earlier clean or one you have created), you can point tclean to it using the mask parameter, mask=’filename.mask’, and tclean will load your masks automatically. This can be useful if you know that some manual cleaning is needed before automasking will start picking up extended emission but several iterations may be needed until you discover that limit. The region file method can also be used non-interactively on it’s own.

Once you are satisfied that you have masked all the obvious emision 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’, 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. 

os.system('rm -rf twhya_advanced_use.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
imagename='twhya_advanced_use',
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.0374Jy',
usemask='auto-multithresh',
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5, 
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
restart=True)

As a word of caution, automasking will not remove any regions, it will only add to the current mask. Therefore, if you create a fluffy, extended mask by hand and then start automasking, it may not be altered by the algorithm. Therefore, during the manual masking stage, mask conservatively.

Retrieved from "https://casaguides.nrao.edu/index.php?title=Automasking_Guide_CASA_5&oldid=36727"