Automasking Guide CASA 6.5.4: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Created page with "== 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 [https://casa.nrao..."
 
 
(103 intermediate revisions by the same user not shown)
Line 1: Line 1:
{{checked_6.5.4}}
== 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’.  
{{CARTA_6.5.4}}
 
Starting with CASA 5.1, a new algorithm has been incorporated into {{tclean_6.5.4}} 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_6.5.4}} 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://casadocs.readthedocs.io/en/v6.5.4/notebooks/synthesis_imaging.html#Automasking CASA docs] 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.


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).
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).
Line 9: Line 13:
== Parameters ==
== 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.
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_6.5.4}} 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.  
There are five primary auto-multithresh subparameters that control what thresholds the algorithm uses and how regions are removed from the mask via pruning.
# '''''Noisethreshold''''' (type: double) -- sets the signal-to-noise threshold above which significant emission is masked during the initial round of mask creation. Note that either ''noisethreshold'' or ''sidelobethreshold'' is used depending on which threshold is higher.  
# '''''noisethreshold''''' (type: double) -- sets the signal-to-noise threshold above which significant emission is masked during the initial round of mask creation. Note that either ''noisethreshold'' or ''sidelobethreshold'' is used depending on which threshold is higher.
# '''''Sidelobethreshold''''' (type: double) -- sets a threshold based on the sidelobe level above which significant emission is masked during the initial round of mask creation. Note that either ''noisethreshold'' or ''sidelobethreshold'' is used depending on which threshold is higher.  
# '''''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.
# '''''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.
# '''''Lownoisethreshold''''' (type: double) -- sets the threshold into which the initial mask (which is determined by either ''noisethreshold'' or ''sidelobethreshold'') is expanded in order to include low signal-to-noise regions in the mask.
# '''''minbeamfrac''''' (type: double) -- sets the minimum size a region must be to be retained in the mask. The parameter is specified as a fractional beam size. Any masks smaller than this will be removed (or pruned) from the mask. Note that this parameter is used to control the pruning for both the initial threshold mask and the low signal-to-noise mask.
# '''''Negativethreshold''''' (type: double) -- sets the signal-to-noise threshold for absorption features to be masked. Note that any absorption features that are masked are not pruned or expanded into low signal-to-noise regions. This parameter should have a positive value, i.e., ''negativethreshold''=5.0 will create a mask for values less than -5.0 sigma.
# '''''negativethreshold''''' (type: double) -- sets the signal-to-noise threshold for absorption features to be masked. Note that any absorption features that are masked are not pruned or expanded into low signal-to-noise regions. This parameter should have a positive value, i.e., ''negativethreshold''=5.0 will create a mask for values less than -5.0 sigma.


Auto-multithresh has several additional subparameters that usually do not need to be changed from their default values:
<ol type="1" start=6>
<li> '''''pbmask''''' (type: double) -- primary beam cutoff
<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> '''''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> '''''verbose''''' (type: bool) -- turns on/off per-channel logging information on the masks (default = False)
</ol>


Auto-multithresh has several additional parameters that usually do not need to be changed from their default values:
Additionally, the following parameter should be considered when using automasking (this not a subparameter of auto-multithresh).
<ol type="1" start=6>
<ol type="1" start=13>
<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> '''''fastnoise''''' (type: bool) -- If True, estimate the noise using the median absolute deviation. If False, estimate the noise using the chauvenet (no mask) or the median absolute deviations of pixels outside the mask but inside the primary beam mask. The latter estimation may be more accurate when emission covers a significant portion of the field of view (default=True).
<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> '''''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_6.5.4}} 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 ==
== 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.
The values presented in this table are the standard values that the pipeline uses when it performs automasking and have been extensively tested.  The values  listed for combined data in the last column are tentative and should serve as a starting point for your analysis. This table can be found in the [https://almascience.nrao.edu/documents-and-tools/cycle10/alma_pipeline_users_guide_2023#subsection.9.39 Pipeline User's Guide section 9.39 (hif_makeimages)].


{| class="wikitable"
{| class="wikitable"
! Parameter
! 7m (continnum/line)
! 12m (b75 < 300m)
! 12m (300m < b75 < 400m)
! 12m (b75 > 400m)
! 12m + 7m combined ''TENTATIVE''
|-
| style="text-align:center;" | ''noisethreshold''
| 5.0
| '''4.25'''
| 5.0
| 5.0
| 4.25
|-
|-
| Array
| style="text-align:center;" | ''sidelobethreshold''
|style="text-align:center;"| ''sidelobethreshold''
| 1.25
| ''noisethreshold''
| '''2.0'''
| ''minbeamfrac''
| 2.0
| ''lownoisethreshold''
| 2.5
| ''negativethreshold''
| 2.0
|-
|-
| 12m (short) b75<300m
| style="text-align:center;" | ''lownoisethreshold''
| 2.0
| 2.0
| 4.25
| '''1.5'''
| 0.3
| 1.5
| 1.5
| 1.5
| 1.5
| 0.0 (continuum)/15.0 (line)
|-
|-
| 12m (long) b75>300m
| style="text-align:center;" | ''minbeamfrac''
| 3.0
| 0.1
| 5.0
| '''0.3'''
| 0.3
| 0.3
| 0.3
| 0.3
| 1.5
| 0.0 (continuum)/7.0 (line)
|-
|-
| 7m (continnum/line)
| style="text-align:center;" | ''negativethreshold''
| 1.25
| 0.0
| 5.0
| '''0.0 (continuum)<br>15.0 (line)'''
| 0.1
| 0.0 (continuum)<br>7.0 (line)
| 2.0
| 0.0 (continuum)<br>7.0 (line)
| 0.0
| 0.0
|-
|-
| 12m + 7m combined '''TENTATIVE'''
| style="text-align:center;" | ''fastnoise''
| 2.0
| False
| 4.25
| '''False'''
| 0.3
| False
| 1.5
| True
| 0.0
| False
|}
|}


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
The term “b75” refers to the 75th percentile of baselines. You can find this value in the pipeline weblog by going to the antenna configuration page and clicking on the “baselines” tab, the table lists the baseline length in increasing distance and the percentile. If you don’t have a weblog for your data, you may use the Analysis Utilities task [https://safe.nrao.edu/wiki/bin/view/ALMA/GetBaselineStats au.getBaselineStats]. You can obtain the Analysis Utilities tasks by following the [https://casaguides.nrao.edu/index.php/Analysis_Utilities Analysis Utilities CASA Guide].


The 75% baseline split is only valid for ALMA data. For data from other telescopes, evaluate the uv-coverage  of the observation. If the uv-coverage of a data set is good, then the PSF will have relatively low sidelobes and the default 12m parameters are likely a good start. We have had success in limited testing using the short baseline ALMA parameters for non-snapshot VLA observations. If the uv-coverage of a data set is poor, then the PSF will have higher sidelobes and the 7m array parameters are probably a better start. An example of this case would be snapshot VLA observations or ATCA observations.
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 ==
If you are using pipeline imaging tasks, these automasking parameters are also able to be defined within those task calls. See the [https://almascience.nrao.edu/documents-and-tools/cycle10/reference-manual-2023#section.2.12 Pipeline Tasks Reference Manual section 2.12].
 
== Guidelines for Optimizing auto-multithresh Parameters ==


* Start with the pipeline parameters for the data with the most similar parameters to your data set. When in doubt, choose the default short baseline 12m parameters.
* 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.
Line 84: Line 107:
* Change one parameter at a time and keep track of your previous images so that you can directly compare different runs.
* 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.
* 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.
* Save the masks for each major cycle. To save the intermediate masks, type the following on the casa command line:
 
<source lang="python">
#In CASA
import os
os.environ['SAVE_ALL_AUTOMASKS']="true"
</source>
 
This feature does NOT work when {{tclean_6.5.4}} is run with ''parallel=True'' on cubes.  
* You can also inspect the intermediate masks by setting ''interactive=True''. However, the resulting masks will NOT necessarily be the same as masks generated with ''interactive=False'' because how often the major cycle is triggered can be different in the two cases. Set ''cycleniter=niter'' when using non-interactive mode to get identical behavior.


== Exploring the auto-multithresh parameters ==
== 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.
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
=== Obtain the Data ===
 
To demonstrate the use of the automasking parameters, we will be using the TW Hydra dataset which is a small dataset used in the “First Look” CASA Guides. Specifically, we will use a flagged, calibrated, self-calibrated, and continuum subtracted Band 7 Measurement Set that is created as part of the [[First_Look_at_Line_Imaging]] guide. Download and untar: [https://bulk.cv.nrao.edu/almadata/public/ALMA_firstlooks/twhya_selfcal.ms.contsub.tar twhya_selfcal.ms.contsub.tar] (~300 MB).
 
<source lang="bash">
#In bash
tar xvf twhya_selfcal.ms.contsub.tar
</source>


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.  
Almost all of the values used in our following {{tclean_6.5.4}} commands remain constant throughout this section, and are explained in and taken from the [[First_Look_at_Line_Imaging]] guide. The main difference is that we are using the ''usemask=’auto-multithresh’'' parameter to use the automasking feature and will be adjusting the auto-masking subparameters to highlight their use.


=== Make a Dirty Cube ===
=== 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.  
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_6.5.4}} first to get the basic information of the dataset.  


<source lang="python">
<source lang="python">
#In CASA
#In CASA
listobs('sis14_twhya_calibrated_flagged.ms.contsub')
listobs('twhya_selfcal.ms.contsub', listfile='twhya_selfcal.ms.contsub.listobs.txt')
</source>
</source>


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Observer: cqi    Project: uid://A002/X327408/X6f   
  Observer: cqi    Project: uid://A002/X327408/X6f   
Observation: ALMA
Observation: ALMA
Data records: 53161       Total elapsed time = 4268.11 seconds
Data records: 44772       Total elapsed time = 4268.11 seconds
   Observed from  19-Nov-2012/07:56:23.5  to  19-Nov-2012/09:07:31.6 (UTC)
   Observed from  19-Nov-2012/07:56:23.5  to  19-Nov-2012/09:07:31.6 (UTC)


   ObservationID = 0        ArrayID = 0
   ObservationID = 0        ArrayID = 0
   Date        Timerange (UTC)          Scan  FldId FieldName            nRows    SpwIds  Average Interval(s)    ScanIntent
   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]
   19-Nov-2012/07:56:23.5 - 08:02:11.3    12      5 TW Hya                    7616 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:08:09.6 - 08:13:57.3    16      0 TW Hya                   10360 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:08:09.6 - 08:13:57.3    16      5 TW Hya                   8442 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:19:53.9 - 08:25:41.7    20      0 TW Hya                   10321 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:19:53.9 - 08:25:41.7    20      5 TW Hya                   8389 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:32:00.5 - 08:37:48.2    24      0 TW Hya                   10324 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:32:00.5 - 08:37:48.2    24      5 TW Hya                   8409 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:43:45.6 - 08:49:33.4    28      0 TW Hya                    9462 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               08:43:45.6 - 08:49:33.4    28      5 TW Hya                    8514 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               09:05:15.6 - 09:07:31.6    36      0 TW Hya                    4180 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
               09:05:15.6 - 09:07:31.6    36      5 TW Hya                    3402 [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
           (nRows = Total number of rows per scan)
           (nRows = Total number of rows per scan)  
Fields: 1
Fields: 1
   ID  Code Name                RA              Decl          Epoch  SrcId      nRows
   ID  Code Name                RA              Decl          Epoch  SrcId      nRows
   0   none TW Hya              11:01:51.796000 -34.42.17.36600 J2000  0         53161
   5   none TW Hya              11:01:51.796000 -34.42.17.36600 J2000  4         44772
Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
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   
   SpwID  Name                          #Chans  Frame  Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz) BBC Num  Corrs   
Line 124: Line 163:
</pre>
</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.  
With the information above and from what we know of TW Hydra from the [[First_Look_at_Line_Imaging]] guide, we get the basic information for our {{tclean_6.5.4}} call to create the dirty cube (set ''niter=0'' to perform no cleaning).
 
[[Image:twhya_dirtycube.image.rms_CASA_6.5.4.png|300px|thumb|<caption>Figure 1: In CARTA, use a region in a line-free channel to determine RMS, then use a multiple of this (~3x) for the threshold parameter.</caption>]]


<source lang="python">
<source lang="python">
#In CASA
#In CASA
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
os.system('rm -rf twhya_dirtycube.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_dirtycube',
imagename='twhya_dirtycube',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 141: Line 183:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=0)
niter=0,
</source>
threshold='0.0Jy')
 
=== Determine Clean Threshold ===
In this tutorial we will clean non-interactively, and thus need to set a threshold to avoid over-cleaning emission. This is true for all non-interactive cleaning, and also a good idea for interactive cleaning (the ''threshold'' parameter is not a subparameter of ''auto-multithresh''). After {{tclean_6.5.4}} has finished, you should now open the dirty image in CARTA. If using NRAO machines, you can open a new terminal tab, cd to the working directory, then type:
 
<source lang="bash">
#In bash
carta --no_browser
</source>
</source>


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.)
Copy the output URL into a browser to view your CARTA session. Select and load '''twhya_dirtycube.image'''. Selecting a region in a line free channel (e.g., channel 3) we find that the RMS is ~29 mJy/beam. We choose the threshold to be 3x the RMS: ''threshold='87mJy' ''(The "/beam" is implied).


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.
<br clear="all">


=== Base Parameters ===  
=== Base Parameters ===  
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.  
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. If you are not using an NRAO machine, follow the [[Analysis_Utilities]] guide.


<source lang="python">
<source lang="python">
#In CASA
#In CASA
au.getBaselineStats('sis14_twhya_calibrated_flagged.ms.contsub')
import analysisUtils as au
au.getBaselineStats('twhya_selfcal.ms.contsub')
</source>
</source>
<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Found 1 obsIDs
Found 325 baselines (26 stations)
Found 325 baselines (26 stations)
Unprojected lengths:  min=15.063874, max=374.719747, rms=162.602204 m
Unprojected lengths:  min=15.063874, max=374.719747, rms=162.602204 m
Line 169: Line 220:
</pre>
</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.  
From the output above (which was printed to the terminal window), we see that b75 (indicated as 75%ile) is 197.4m, so we will use the pipeline default values in the '''"12m (b75 < 300m)"''' column of the [[#Table_of_Standard_Values]] above. This will give us a mask to compare our testing against which will be displayed as black contours in all of the following figures.  
 
 
[[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>]]


[[Image:automask_base_parameters_CASA_6.5.4.png|800px|thumb|<caption>Figure 2: In CARTA, we display channel 7 of the cube created with the base parameters for '''12m (b75 < 300m)'''. On the left is the image, and on the right is the residual. The black contours mark the final mask created for this channel.</caption>]]


<source lang="python">
<source lang="python">
#In CASA
#In CASA
os.system('rm -rf twhya_base_params.*')
os.system('rm -rf twhya_base_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_base_params',
imagename='twhya_base_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 191: Line 240:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking Parameters below this line
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,  
lownoisethreshold=1.5,  
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0,
verbose=True)
verbose=True,
fastnoise=False)
</source>
</source>


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.  
# Open the following three images in CARTA (with File > Open and ctrl+click, or one at a time with File > Append): </br> ''twhya_base_params.image'' </br> ''twhya_base_params.residual'' </br> ''twhya_base_params.mask'' </br>
# In the Image List, lock all 3 images in XY (spatial matching), Z (spectral matching), and R (raster scaling matching). If '''.image''' is not already set as the reference for matching, right click the image name and select "Set as ... reference." To reorder images in the list, click an index, then click-and-drag the index.
# Click the settings gear of the Image Viewer, go to the Global tab, enable multi-panel, and choose 2 columns and 1 row. Adjust the images to an appropriate position (Zoom to fit is recommended).
# Use the Animator to page through the channels. Notice that the mask image is blank for all channels except 5 and 7. For this example we choose channel 7.
# Adjust the raster scaling in Render Configuration.
# Open the Contour dialog. Unlock data source. Choose the mask. The Histogram, Generator mode, and Parameters provide different ways to calculate Levels, however for a mask contour we can ignore these. A mask image has value 1.0 where the mask exists, and 0.0 elsewhere (see this by moving your cursor around a mask image). Enter 0.5 directly into the Levels box, then click apply. This displays a contour outline of the mask on all 3 images. Go to the Styling tab to adjust color and thickness.
 
With ''verbose=True'' set, automasking reports its progress per-channel in the logger window. Below is the logger output for the first major clean cycle using the base parameters above.
 
* The minimum mask size and sidelobeLevel are calculated in order to define different thresholds based on the user input.


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
[twhya_base_params] Peak residual (max,min) over full image : (0.149898,-0.081815)
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[twhya_base_params] Total Model Flux : 0
[twhya_base_params] Setting up an auto-mask within PB mask limit
[twhya_base_params] Setting up an auto-mask
---------------------------------------------------- Run Automask  ---------------------------------------------
Set Deconvolution Options for [twhya_base_params] : hogbom
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
Generating AutoMask
Generating AutoMask
SidelobeLevel = 0.219791
[C0] Using Chauvenet algorithm for the image statistics
prune size=13.535(minbeamfrac=0.3 * beampix=45.1166)
[C1] Using Chauvenet algorithm for the image statistics
Start thresholding: create an initial mask by threshold
[C2] Using Chauvenet algorithm for the image statistics
End thresholding: time to create the initial threshold mask:  real 0.37s ( user 0.02s, system 0.02s)
[C3] Using Chauvenet algorithm for the image statistics
Start pruning: the initial threshold mask
[C4] Using Chauvenet algorithm for the image statistics
[C0] No regions are found in this plane.
[C5] Using Chauvenet algorithm for the image statistics
[C1] No regions are found in this plane.
[C6] Using Chauvenet algorithm for the image statistics
[C2] No regions are found in this plane.
[C7] Using Chauvenet algorithm for the image statistics
[C3] No regions are found in this plane.
[C8] Using Chauvenet algorithm for the image statistics
[C7] No regions are removed in pruning process.
[C9] Using Chauvenet algorithm for the image statistics
[C8] No regions are found in this plane.
[C10] Using Chauvenet algorithm for the image statistics
[C9] No regions are found in this plane.
[C11] Using Chauvenet algorithm for the image statistics
[C10] No regions are found in this plane.
[C12] Using Chauvenet algorithm for the image statistics
[C11] No regions are found in this plane.
[C13] Using Chauvenet algorithm for the image statistics
[C12] No regions are found in this plane.
[C14] Using Chauvenet algorithm for the image statistics
[C14] No regions are found in this plane.
SidelobeLevel = 0.21961
End pruning: time to prune the initial threshold mask: real 0.24s (user 0.05s, system 0.02s)
prune size=9.47266(minbeamfrac=0.3 * beampix=31.5755)
Start smoothing: the initial threshold mask
</pre>
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
* Masking is performed and then pruned. Only "error" messages are reported per channel (i.e. no regions removed or no regions found).
* The mask is smoothed and then in later cycles grown into lower signal-to-noise regions and pruned again.
* Negative emission is masked if present.
 
<pre style="background-color: #fffacd;">
*** Start auto-multithresh processing for Channel 0***
No regions are found in this plane.
No negative region was found by auotmask.
No negative region was found by auotmask.
End thresholding: time to create the negative mask: real 1.19s (user 0.54s, system 0.26s)
*** Start auto-multithresh processing for Channel 1***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 2***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 3***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 4***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 5***
No regions are removed in pruning process.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 6***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 7***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 8***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 9***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 10***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 11***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 12***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 13***
No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 14***
No regions are found in this plane.
No negative region was found by auotmask.
*** Timing summary for whole planes ***
Total time to create the initial threshold mask: real 0.02s ( user 0.48s, system 0s)
Total time to prune the initial threshold mask: real 0.03s (user 0.66s, system 0.01s)
Total time to create the smoothed initial threshold mask: real 0.75s (user 2.13s, system 1.19s)
Total time to grow the previous mask: real 0s (user 0s, system 0s)
Total time to prune the grow mask: real 0s (user 0s, system 0s)
Total time to create the smoothed grow mask: real 0s (user 0s, system 0s)
</pre>
 
* Automasking evaluates whether or not to continue to try masking each channel. If no masks were found in a channel or all masks were removed through the pruning process, the channel is no longer considered in the upcoming minor clean cycles.
 
<pre style="background-color: #fffacd;">
Stopping masking for chan=0
Stopping masking for chan=0
Stopping masking for chan=1
Stopping masking for chan=1
Line 245: Line 360:
Stopping masking for chan=3
Stopping masking for chan=3
Stopping masking for chan=4
Stopping masking for chan=4
Stopping masking for chan=6
Stopping masking for chan=8
Stopping masking for chan=8
Stopping masking for chan=9
Stopping masking for chan=9
Line 252: Line 368:
Stopping masking for chan=13
Stopping masking for chan=13
Stopping masking for chan=14
Stopping masking for chan=14
</pre>
* A summary of the automasking process is given in the "Automask summary" table. Here, you can see the threshold that was used to create the mask (either "noise" or "sidelobe"), the number of regions found, number of regions pruned, etc. Our automasking results in masks for channels 5 and 7.
<pre style="background-color: #fffacd;">
========== automask summary ==========
========== automask summary ==========
chan masking? median  RMS        peak  thresh_type  thresh_value  N_reg N_pruned N_grow N_grow_pruned N_neg_pix
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
[C0] F        -- -- 0.12525 noise  --       0  0  0 0 0
[C1] F        -6.20534e-05 0.0185451 -0.0749006 noise  0.0788167       0  0  -- -- 0
[C1] F        --  -- -0.113465 noise  --       0  0  0 0 0
[C2] F        -2.53807e-06 0.0193195 0.0714881 noise  0.0821077       0  0 --  -- 0
[C2] F        --  -- 0.136005 noise  --       1  1  0  0  0
[C3] F        7.66067e-05 0.0185695 -0.0724906 noise  0.0789204       0  0 --  -- 0
[C3] F        -- -- 0.148636 noise  --       1  1  0  0  0
[C4] F        -0.000102242 0.0192787 0.0862985 noise  0.0819347        2 2 -- -- 0
[C4] F        -- -- 0.119875 noise  --        0  0 0 0 0
[C5] T        -0.000765545 0.0224421 0.133081 noise  0.0953791       2 1 -- -- 0
[C5] T        0.000376255 0.0328259 0.190839 noise  0.139886       1 0 0 0 0
[C6] T       -0.00150215 0.0256315 0.149898 noise  0.108934       2 1  -- -- 0
[C6] F       -- -- 0.152505 noise  --       1 1  0 0 0
[C7] T        -0.000590894 0.0212334 0.116813 noise  0.0902421       1  0  --  -- 0
[C7] T        -3.61604e-05 0.0308196 0.166526 noise  0.130947       2  1  0  0 0
[C8] F        1.70437e-05 0.0176792 -0.0750749 noise  0.0751365       0  0  -- -- 0
[C8] F        -- -- 0.102411 noise  --       0  0  0 0 0
[C9] F        0.000338539 0.0175967 0.0657392 noise  0.0747862       0  0  -- -- 0
[C9] F        -- -- 0.106988 noise  --       0  0  0 0 0
[C10] F        -3.07989e-05 0.0178246 -0.0780685 noise  0.0757547       0  0  -- -- 0
[C10] F        --  -- -0.117285 noise  --       0  0  0 0 0
[C11] F        -0.000249928 0.0176007 -0.0744855 noise  0.074803       0  0  -- -- 0
[C11] F        -- -- 0.109285 noise  --       0  0  0 0 0
[C12] F        0.000142292 0.0176476 -0.0685521 noise  0.0750021       0  0  -- -- 0
[C12] F        -- -- 0.108688 noise  --       0  0  0 0 0
[C13] F        -0.000252258 0.0178448 0.0768286 noise  0.0758403        1 1 -- -- 0
[C13] F        -- -- -0.111106 noise  --        0  0 0 0 0
[C14] F        0.000200417 0.0171804 0.0673071 noise  0.0730168       0  0  -- -- 0
[C14] F        -- -- -0.106717 noise  --       0  0  0 0 0
========== END of automask summary ==========
========== END of automask summary ==========
</pre>
* Cleaning is performed on the masked image and the progress (peak residual, number of iterations, etc.) is reported per channel. Notice how minor cycles (iterations/iters) are performed on our masked channels 5 and 7 only.
<pre style="background-color: #fffacd;">
----------------------------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_base_params] Number of pixels in the clean mask : 1209 out of a total of 937500 pixels. [ 0.12896 % ]
[twhya_base_params] Number of pixels in the clean mask : 362 out of a total of 937500 pixels. [ 0.0386133 % ]
[twhya_base_params] Peak residual (max,min) within mask : (0.149898,-0.0296383) over full image : (0.149898,-0.081815)
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[twhya_base_params] Total Model Flux : 0
Minor Cycle controls : {'cycleniter': 100000, 'cyclethreshold': 0.08699999749660492, 'loopgain': 0.10000000149011612, 'nsigma': 0.0, 'thresholdreached': True}
    itsNsigma=0
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
[twhya_base_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.0356, CycleNiter=100000, Gain=0.1
Set Deconvolution Options for [twhya_base_params] : hogbom
[twhya_base_params:C0] iters=0->0 [0], model=0->0, peakres=0.081815->0.081815, Skipped this plane. Zero mask.
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[twhya_base_params:C1] iters=0->0 [0], model=0->0, peakres=0.0749006->0.0749006, Skipped this plane. Zero mask.
Processing channels in range [0, 14]
[twhya_base_params:C2] iters=0->0 [0], model=0->0, peakres=0.0714881->0.0714881, Skipped this plane. Zero mask.
    itsNsigma=0
[twhya_base_params:C3] iters=0->0 [0], model=0->0, peakres=0.0724906->0.0724906, Skipped this plane. Zero mask.
[twhya_base_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.087, CycleNiter=100000, Gain=0.1
[twhya_base_params:C4] iters=0->0 [0], model=0->0, peakres=0.0862985->0.0862985, Skipped this plane. Zero mask.
[twhya_base_params:C0] iters=0->0 [0], model=0->0, peakres=0.12525->0.12525, 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:C1] iters=0->0 [0], model=0->0, peakres=0.113465->0.113465, Skipped this plane. Zero mask.
[twhya_base_params:C6] iters=85->177 [92], model=0->0.564421, peakres=0.149898->0.0355323, Reached cyclethreshold.
[twhya_base_params:C2] iters=0->0 [0], model=0->0, peakres=0.136005->0.136005, Skipped this plane. Zero mask.
[twhya_base_params:C7] iters=177->256 [79], model=0->0.450467, peakres=0.116813->0.0349898, Reached cyclethreshold.
[twhya_base_params:C3] iters=0->0 [0], model=0->0, peakres=0.148636->0.148636, Skipped this plane. Zero mask.
[twhya_base_params:C8] iters=256->256 [0], model=0->0, peakres=0.0750749->0.0750749, Skipped this plane. Zero mask.
[twhya_base_params:C4] iters=0->0 [0], model=0->0, peakres=0.119875->0.119875, Skipped this plane. Zero mask.
[twhya_base_params:C9] iters=256->256 [0], model=0->0, peakres=0.0657392->0.0657392, Skipped this plane. Zero mask.
[twhya_base_params:C5] iters=0->16 [16], model=0->0.190617, peakres=0.190839->0.085914, Reached cyclethreshold.
[twhya_base_params:C10] iters=256->256 [0], model=0->0, peakres=0.0780685->0.0780685, Skipped this plane. Zero mask.
[twhya_base_params:C6] iters=16->16 [0], model=0->0, peakres=0.152505->0.152505, 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:C7] iters=16->39 [23], model=0->0.258554, peakres=0.166526->0.0863384, Reached cyclethreshold.
[twhya_base_params:C12] iters=256->256 [0], model=0->0, peakres=0.0685521->0.0685521, Skipped this plane. Zero mask.
[twhya_base_params:C8] iters=39->39 [0], model=0->0, peakres=0.102411->0.102411, 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:C9] iters=39->39 [0], model=0->0, peakres=0.106988->0.106988, 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:C10] iters=39->39 [0], model=0->0, peakres=0.117285->0.117285, 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
[twhya_base_params:C11] iters=39->39 [0], model=0->0, peakres=0.109285->0.109285, Skipped this plane. Zero mask.
Completed 256 iterations.
[twhya_base_params:C12] iters=39->39 [0], model=0->0, peakres=0.108688->0.108688, Skipped this plane. Zero mask.
[twhya_base_params:C13] iters=39->39 [0], model=0->0, peakres=0.111106->0.111106, Skipped this plane. Zero mask.
[twhya_base_params:C14] iters=39->39 [0], model=0->0, peakres=0.106717->0.106717, Skipped this plane. Zero mask.
[twhya_base_params] Total model flux (over all planes) : 0.44917     Peak Residual (over all planes) : 0.0863384 in C7:P0
Completed 39 iterations.
</pre>
</pre>


=== Initial mask threshold calculation (sidelobethreshold and noisethreshold parameters)===
=== Initial Mask Threshold Calculation ===


The initial mask threshold is calculated by the following method:
The initial mask threshold is calculated by the following method:


<blockquote>
<blockquote>
sidelobeThresholdValue = sidelobeThreshold * sidelobeLevel * peak in residual image <br>
noiseThresholdValue = ''noisethreshold'' * rms in residual image <br>
noiseThresholdValue = noiseThreshold * rms in residual image
sidelobeThresholdValue = ''sidelobethreshold'' * sidelobeLevel * peak in residual image
</blockquote>
</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.
where the sidelobeLevel is calculated by {{tclean_6.5.4}} and reported in the logger (you can see this in the example logger output above).


The initial threshold is max(sidelobeThresholdValue,noiseThresholdValue).  
The initial threshold is max(noiseThresholdValue, sidelobeThresholdValue).  


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.
In practice, the ''sidelobethreshold'' often sets the limit in the first few major cycles because the sidelobe emission is often larger than the random noise in the image. Once the strongest emission has been removed from the residual image, later cycles tend to use the noisethreshold.


==== Turning Down ''noisethreshold'' ====
==== 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.75.   
 
 
[[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>]]


[[Image:small_noisethreshold_CASA_6.5.4.png|600px|thumb|<caption>Figure 3: On the left is the image and on the right is the residual of channel 7 of 15 of the cube created with these parameters. The black contours represent the mask created when using the base parameters, and the white contours represent the masks created by setting a small ''noisethreshold'' value, which favors the use of ''sidelobethreshold''.</caption>]]


<source lang="python">
<source lang="python">
#In CASA
#In CASA
os.system('rm -rf twhya_lownoise_params.*')
os.system('rm -rf twhya_low_noise_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_lownoise_params',
imagename='twhya_low_noise_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 332: Line 461:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
noisethreshold=1.75, # lowered from base
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=1.5,
lownoisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0,
verbose=True)
verbose=True,
fastnoise=False)
</source>
</source>


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.  
Open the following in CARTA: </br> ''twhya_low_noise_params.image'' </br> ''twhya_low_noise_params.residual'' </br> ''twhya_low_noise_params.mask'' </br> ''twhya_base_params.mask''
 
Using the method described above, display the '''.image''' and '''.residual''' sibe-by-side, overlay the base mask as black contours, and overlay the lownoise mask as white contours.
 
By reducing ''noisethreshold'', auto-multithresh is now using ''sidelobethreshold'' instead to draw the initial masks for most channels. You can see this in the logger. The "thresh_type" column of the "automask summary" reports which threshold is being used to find maskable emission in each channel. The output below is trimmed for brevity indicated by "*****".


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
[twhya_lownoise_params] Peak residual (max,min) over full image : (0.149898,-0.081815)
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[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
---------------------------------------------------- Run Automask  ---------------------------------------------
Set Deconvolution Options for [twhya_lownoise_params] : hogbom
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
Generating AutoMask
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)
SidelobeLevel = 0.21961
Start pruning: the initial threshold mask
prune size=9.47266(minbeamfrac=0.3 * beampix=31.5755)
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 ==========
========== automask summary ==========
chan masking? median  RMS        peak  thresh_type  thresh_value  N_reg N_pruned N_grow N_grow_pruned N_neg_pix
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
[C0] T        0.000326007 0.0306058 0.12525 sidelobe  0.0553385       206 152 0 0 0
[C1] T        -6.20534e-05 0.0185451 -0.0749006 sidelobe  0.0304619       254 153 -- -- 0
[C1] T        0.000123019 0.0282244 -0.113465 sidelobe  0.0499591       206 152 0 0 0
[C2] T        -2.53807e-06 0.0193195 0.0714881 sidelobe  0.0314249       261 170 -- -- 0
[C2] T        4.33802e-05 0.03036 0.136005 sidelobe  0.0597796       150 105 0 0 0
[C3] T        7.66067e-05 0.0185695 -0.0724906 sidelobe  0.0306932       267 171 -- -- 0
[C3] T        0.000242864 0.0281953 0.148636 sidelobe  0.065527       95 85 0 0 0
[C4] T        -0.000102242 0.0192787 0.0862985 sidelobe  0.0379353       172 136 -- -- 0
[C4] T        0.000411201 0.0290491 0.119875 sidelobe  0.0530629       210 156 0 0 0
[C5] T        -0.000765545 0.0224421 0.133081 sidelobe  0.0585003       28 25 -- -- 0
[C5] T        0.000376255 0.0328259 0.190839 sidelobe  0.0841967       43 40 0 0 0
[C6] T        -0.00150215 0.0256315 0.149898 sidelobe  0.0658923       27 19 -- -- 0
[C6] T        0.000135499 0.0357602 0.152505 sidelobe  0.0671187       126 82 0 0 0
[C7] T        -0.000590894 0.0212334 0.116813 sidelobe  0.0513489       44 37 -- -- 0
[C7] T        -3.61604e-05 0.0308196 0.166526 sidelobe  0.0731053       56 40 0 0 0
[C8] T        1.70437e-05 0.0176792 -0.0750749 sidelobe 0.0326501       204 151 -- -- 0
[C8] T        0.000578858 0.0269325 0.102411 noise 0.0477108       233 178 0 0 0
[C9] T        0.000338539 0.0175967 0.0657392 sidelobe  0.0288978       266 177 -- -- 0
[C9] T        5.75308e-05 0.0267095 0.106988 sidelobe  0.0470487       229 155 0 0 0
[C10] T        -3.07989e-05  0.0178246 -0.0780685 sidelobe  0.0288354       278 186 -- -- 0
[C10] T        -3.76976e-05  0.0264451 -0.117285 sidelobe  0.0514762       189 154 0 0 0
[C11] T        -0.000249928 0.0176007 -0.0744855 sidelobe  0.0281572       281 184 -- -- 0
[C11] T        5.7293e-06 0.0268588 0.109285 sidelobe  0.048006       212 150 0 0 0
[C12] T        0.000142292 0.0176476 -0.0685521 sidelobe 0.029941       225 155 -- -- 0
[C12] T        -0.000308202 0.0278082 0.108688 noise 0.0483562       219 154 0 0 0
[C13] T        -0.000252258 0.0178448 0.0768286 sidelobe  0.0337725       184 146 -- -- 0
[C13] T        1.40047e-05 0.0268257 -0.111106 sidelobe  0.0488139       214 159 0 0 0
[C14] T        0.000200417 0.0171804 0.0673071 sidelobe  0.029587       250 182 -- -- 0
[C14] T        0.00029441 0.0257197 -0.106717 sidelobe  0.0471667       220 160 0 0 0
========== END of automask summary ==========
========== END of automask summary ==========
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_lownoise_params] Number of pixels in the clean mask : 109074 out of a total of 937500 pixels. [ 11.6346 % ]
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Minor Cycle controls : {'cycleniter': 100000, 'cyclethreshold': 0.08699999749660492, 'loopgain': 0.10000000149011612, 'nsigma': 0.0, 'thresholdreached': True}
    itsNsigma=0
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
Set Deconvolution Options for [twhya_lownoise_params] : hogbom
Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
    itsNsigma=0
[twhya_lownoise_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.087, CycleNiter=100000, Gain=0.1
[twhya_lownoise_params:C0] iters=0->30 [30], model=0->0.245732, peakres=0.12525->0.0869206, Reached cyclethreshold.
[twhya_lownoise_params:C1] iters=30->47 [17], model=0->0.0652784, peakres=0.105279->0.0865448, Reached cyclethreshold.
[twhya_lownoise_params:C2] iters=47->90 [43], model=0->0.310039, peakres=0.136005->0.0869305, Reached cyclethreshold.
[twhya_lownoise_params:C3] iters=90->109 [19], model=0->0.194673, peakres=0.148636->0.0853744, Reached cyclethreshold.
[twhya_lownoise_params:C4] iters=109->137 [28], model=0->0.162311, peakres=0.119875->0.0869172, Reached cyclethreshold.
[twhya_lownoise_params:C5] iters=137->160 [23], model=0->0.258899, peakres=0.190839->0.0869799, Reached cyclethreshold.
[twhya_lownoise_params:C6] iters=160->243 [83], model=0->0.851988, peakres=0.152505->0.0869064, Reached cyclethreshold.
[twhya_lownoise_params:C7] iters=243->288 [45], model=0->0.419446, peakres=0.166526->0.0862288, Reached cyclethreshold.
[twhya_lownoise_params:C8] iters=288->297 [9], model=0->0.00968076, peakres=0.102411->0.0869284, Reached cyclethreshold.
[twhya_lownoise_params:C9] iters=297->310 [13], model=0->0.0652534, peakres=0.106988->0.0856935, Reached cyclethreshold.
[twhya_lownoise_params:C10] iters=310->324 [14], model=0->0.0348655, peakres=0.117285->0.0865639, Reached cyclethreshold.
[twhya_lownoise_params:C11] iters=324->345 [21], model=0->0.105682, peakres=0.109285->0.086327, Reached cyclethreshold.
[twhya_lownoise_params:C12] iters=345->357 [12], model=0->0.0756125, peakres=0.108688->0.0868853, Reached cyclethreshold.
[twhya_lownoise_params:C13] iters=357->372 [15], model=0->0.10012, peakres=0.0996836->0.0860267, Reached cyclethreshold.
[twhya_lownoise_params:C14] iters=372->382 [10], model=0->0.0373442, peakres=0.106717->0.0866348, Reached cyclethreshold.
[twhya_lownoise_params] Total model flux (over all planes) : 2.93693    Peak Residual (over all planes) : 0.0869799 in C5:P0
Completed 382 iterations.
</pre>
</pre>


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.  
Notice that masking and cleaning occurs in all channels. Only channels 8 and 12 use thresh_type = noise. In later major cycles, we can see that more channels switch from sidelobe to noise. The threshold used by auto-multithresh are independent for each channel and can change between major cycles as clean progresses. In Figure 3, we see that fainter emission gets picked up by automasking when we reduce the ''noisethreshold'' value (as expected), but that noise spikes are clearly being picked up in both line and line-free channels. We can also see this clearly in the log as the number of masks created and reported in the automask summary is much, much higher than in the base parameters log output.


==== Turning Down ''sidelobethreshold'' ====
==== 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.  
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.  


 
[[Image:small_sidelobethreshold_CASA_6.5.4.png|600px|thumb|<caption>Figure 4: On the left is the image and on the right is the residual of channel 7 of 15 of the cube created with these parameters. The black contours represent the mask created when using the base parameters and the white contours represent the mask created using a small ''sidelobethreshold'' value.</caption>]]
[[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>]]
 


<source lang="python">
<source lang="python">
#In CASA
#In CASA
os.system('rm -rf twhya_lowsidelobe_params.*')
os.system('rm -rf twhya_low_sidelobe_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_lowsidelobe_params',
imagename='twhya_low_sidelobe_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 414: Line 573:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
sidelobethreshold=1.0,
noisethreshold=4.25,
noisethreshold=4.25,
lownoisethreshold=0.5,
sidelobethreshold=1.0, # lowered from base
lownoisethreshold=0.5, # lowered from base
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0)
verbose=True,
fastnoise=False)
</source>
</source>


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.
In Figure 4 we see that turning down ''sidelobethreshold'' makes a mask that covers a much larger region and clearly covers areas of the image that do not contain real emission.


=== The Effects of Pruning ===
=== 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.
Pruning is used to remove small noise-like regions from the resulting mask. Note that regions are pruned before being convolved with a Gaussian to create a buffer around the mask.
 
==== Turning off pruning ====


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


 
[[Image:small_minbeamfrac_CASA_6.5.4.png|600px|thumb|<caption>Figure 5: 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 white contours represent the mask created with no pruning.</caption>]]
[[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>]]
 


<source lang="python">
<source lang="python">
#In CASA
#In CASA
os.system('rm -rf twhya_lownoise_noprune_params.*')
os.system('rm -rf twhya_low_minbeamfrac_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_lownoise_noprune_params',
imagename='twhya_low_minbeamfrac_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 461: Line 618:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=0.0,
minbeamfrac=0.0, # lowered from base
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0)
verbose=True,
fastnoise=False)
</source>
</source>


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.
In Figure 5 we see that with no pruning, a noise spike near the edge of channel 7 was picked up, kept, and cleaned. Notice that several line-free channels include masks around noise spikes as well.


==== Increasing the Pruning  ====
==== 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.
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">
#In CASA
#In CASA
os.system('rm -rf twhya_highminbeamfrac_params.*')
os.system('rm -rf twhya_high_minbeamfrac_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_highminbeamfrac_params',
imagename='twhya_high_minbeamfrac_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 499: Line 657:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,
lownoisethreshold=1.5,
minbeamfrac=2.0,
minbeamfrac=2.0, # raised from base
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0)
verbose=True,
fastnoise=False)
</source>
</source>


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


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
CASA$ WARN 5 of 15 channels had all regions removed by pruning. Try decreasing minbeamfrac to remove fewer regions.  
SidelobeLevel = 0.21961
CASA$ WARN Restoring with an empty model image. Only residuals will be processed to form the output restored image.  
prune size=63.151(minbeamfrac=2 * beampix=31.5755)
*****
No regions found by automasking
15 of 15 channels used the noise threshold to mask, but no emission was found. Try decreasing your noisethreshold parameter if you want to capture emission in these channels.
5 of 15 channels had all regions removed by pruning. Try decreasing minbeamfrac to remove fewer regions
*****
[twhya_high_minbeamfrac_params] Number of pixels in the clean mask : 0 out of a total of 937500 pixels. [ 0 % ]
*****
Reached global stopping criterion : zero mask
*****
Restoring with an empty model image. Only residuals will be processed to form the output restored image.
</pre>
</pre>


In this case, no cleaning is done, a dirty cube is generated, and tclean is exited.
In this case, no cleaning is done, a dirty cube is generated, and {{tclean_6.5.4}} 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.
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 ===  
=== 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
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>
<blockquote>
lowNoiseThresholdValue = lowNoiseThreshold * rms in residual image <br>
lowNoiseThresholdValue = ''lownoisethreshold'' * rms in residual image <br>
max(sidelobeThresholdValue,lowNoiseThresholdValue).
max(sidelobeThresholdValue, lowNoiseThresholdValue).
</blockquote>
</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.
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.


Finally, we lower the ''lownoisethreshold'' value to try and grow the mask into regions of lower signal-to-noise during the binary dilation stage.  
[[Image:Small_lownoisethreshold_sidelobe_limited_CASA_6.5.4.png|600px|thumb|<caption>Figure 6: On the left is the channel 7 mask created with base parameters, and on the right is the channel 7 mask created by lowering ''lownoisethreshold'' only. There is no difference, because the default value of ''sidelobethreshold'' takes over.</caption>]]


<source lang="python">
<source lang="python">
#In CASA
#In CASA
os.system('rm -rf twhya_lowlownoise_params.*')
os.system('rm -rf twhya_low_lownoise_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_lowlownoise_params',
imagename='twhya_low_lownoise_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 555: Line 726:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.0, # lowered from base
lownoisethreshold=0.5,
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0)
verbose=True,
fastnoise=False)
</source>
</source>


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


 
[[Image:Small_lownoisethreshold_CASA_6.5.4.png|600px|thumb|<caption>Figure 7: On the left is the image and on the right is the residual of channel 7 of 15 of the cube created with these parameters. The black contours represent the mask created when using the base parameters, and the white contours represent the mask created using a small ''lownoisethreshold'' value (with an even lower ''sidelobethreshold'' so that it doesn't become the limiter).</caption>]]
[[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>]]


<source lang="python">
<source lang="python">
#In CASA
os.system('rm -rf twhya_low_lownoise+sidelobe_params.*')
os.system('rm -rf twhya_low_lownoise+sidelobe_params.*')
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_lownoise+sidelobe_params',
imagename='twhya_low_lownoise+sidelobe_params',
field='0',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 592: Line 764:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
sidelobethreshold=0.5,
noisethreshold=4.25,
noisethreshold=4.25,
lownoisethreshold=1.0,
sidelobethreshold=0.5, # lowered from base
lownoisethreshold=1.0, # lowered from base
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0)
verbose=True,
fastnoise=False)
</source>
</source>


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


== Advanced Use Case - Merging User Masks with Automasking ==
== 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 some situations, you may find that automasking is not masking obvious emission even after adjusting the thresholds. This most commonly occurs in cases:
* with extremely bright, diffuse emission throughout the entire field of view
* where there is low-level absorption that has a similar depth to any negative bowls in the same channel
* where the cube has not been continuum subtracted.
 
The first case is due to the difficulty of estimating noise in the presence of wide-spread signal. The last two are cases when the astronomer “knows” where the real signal is (likely through some additional information), but where the algorithm does not have enough knowledge to distinguish.
 
In these cases, it may be helpful to “kick start” automasking by first defining masks by hand for a round of cleaning and then having automasking pick up from there. To continue cleaning an existing image, we simply use the ''restart=True'' option within {{tclean_6.5.4}} (which is the default value).


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_6.5.4}} will let us manually create masks before each major cycle. As mentioned above in the “Make a Dirty Cube” section, in order to get identical cleaning behavior between interactive and non-interactive modes, you must set ''cycleniter'' equal to ''niter''. With this arbitrary very large number (''cycleniter=niter=100000''), the minor cycle iteration limit per major cycle will not be reached. Instead, {{tclean_6.5.4}} will continue cleaning until the internally computed ''cyclethreshold'' or the user specified ''threshold'' is hit. One last parameter to set is ''usemask='user''' to indicate that you would like to manually mask emission. Below is an example {{tclean_6.5.4}} command using the same dataset as above.


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.  
[[Image:twhya_manual_mask_CASA_6.5.4.png|400px|thumb|<caption>Figure 8: The tclean interactive GUI, with a manually drawn mask for channel 7 shown in white.</caption>]]


<source lang="python">
<source lang="python">
os.system('rm -rf twhya_advanced_use.*')
#In CASA
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
os.system('rm -rf twhya_manual+automask.*')
imagename='twhya_advanced_use',
tclean(vis='twhya_selfcal.ms.contsub',
field='0',
imagename='twhya_manual+automask',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 634: Line 815:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=True,
interactive=True, # manual masking in interactive GUI
chanchunks=-1,
niter=100000,
niter=100000,
cycleniter=100000,
cycleniter=100000, # set equal to niter to have interactive mode replicate the behavior of non-interactive mode
threshold='0.0374Jy',
threshold='87mJy',
usemask='user’)
usemask='user' # default value, but shown here explicitly in place of auto-multithresh
)
</source>
</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 emission and cleaned to a level where automasking will start to pick up extended emission, you can exit {{tclean_6.5.4}} by pressing the red “X”. After {{tclean_6.5.4}} finishes, you can start {{tclean_6.5.4}} again but now in non-interactive mode (''interactive=False''), with automasking enabled (''usemask='auto-multithresh'''), and using the appropriate automasking thresholds that match your data (see the “Table of Standard Values” above). Finally, we set ''restart=True'' which will make {{tclean_6.5.4}} pick up where it left off using the existing .mask file, along with the existing .model, .psf, .image, etc. that share the same base ''imagename''.
 
Side tip: if you already have a CASA region file from an earlier clean or one you have created, and is has a different name, you can read it into {{tclean_6.5.4}} using the mask parameter, ''mask='basename.mask'''. 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.  
You may see the following warnings. This is a known issue, and can be ignored for the above command.
<pre style="background-color: #fffacd;">
----------------------------------------------------------- Run Major Cycle 2 -------------------------------------
2024-04-02 23:11:40 WARN task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1572) Error : Error in copying internal T/F mask : Error in deleting internal T/F mask : Region mask0 could not be removed
2024-04-02 23:11:40 WARN task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1572)+ Cannot delete the mask (used in another process)
</pre>
 
<font color="red">
However, this does prevent ''restart=True'' from working correctly. There is a known issue with file locks in {{tclean_6.5.4}}. To use ''restart=True'', exit CASA at this step and start a new session before running the command below.
</font>
 
[[Image:twhya_manual+automask_CASA_6.5.4.png|500px|thumb|<caption>Figure 9: On the left is the image and on the right is the residual of channel 7 of the cube created with manual + automasking. The black contours represent the mask created when using automasking only with the base parameters, and the white contours represent the mask created by first doing an interactive clean and then an automasking non-interactive clean.</caption>]]


<source lang="python">
<source lang="python">
os.system('rm -rf twhya_advanced_use.*')
#In CASA
tclean(vis='sis14_twhya_calibrated_flagged.ms.contsub',
# os.system('rm -rf twhya_manual+automask.*')
imagename='twhya_advanced_use',
tclean(vis='twhya_selfcal.ms.contsub',
field='0',
imagename='twhya_manual+automask',
field='5',
spw='0',
spw='0',
specmode='cube',
specmode='cube',
Line 665: Line 860:
gridder='standard',
gridder='standard',
imsize=[250,250],
imsize=[250,250],
cell='0.08arcsec',
cell='0.1arcsec',
weighting='briggs',
weighting='briggsbwtaper',
robust=0.5,
robust=0.5,
restoringbeam='common',
restoringbeam='common',
interactive=False,
interactive=False,
chanchunks=-1,
niter=100000,
niter=100000,
threshold='0.0374Jy',
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
sidelobethreshold=2.0,
noisethreshold=4.25,
lownoisethreshold=1.5,  
lownoisethreshold=1.5,  
minbeamfrac=0.3,
minbeamfrac=0.3,
growiterations=75,
negativethreshold=15.0,
negativethreshold=15.0,
restart=True)
verbose=True,
fastnoise=False,
restart=True # default value, but shown here explicitly
)
</source>
</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.
As a word of caution, automasking will not remove any regions, it will only add to the current mask. Therefore, during the manual masking stage, mask conservatively.
If you create a fluffy, extended mask by hand and then start automasking, it may not be altered by the algorithm.
 
== Version History ==
 
The older version of this guide is here: [[Automasking Guide CASA 5]]

Latest revision as of 19:10, 6 August 2024

Last checked on CASA Version 6.5.4

Introduction

This guide features CARTA, the “Cube Analysis and Rendering Tool for Astronomy,” which is the new NRAO visualization tool for images and cubes. The CASA viewer (imview) has not been maintained for a few years and will be removed from future versions of CASA. We strongly recommend using CARTA, as it provides a much more efficient, stable, and feature rich user experience. A comparison of the CASA viewer and CARTA, as well as instructions on how to use CARTA at NRAO, is provided in the CARTA section of the CASA docs.

Starting with CASA 5.1, a new algorithm has been incorporated into tclean to automatically mask regions during the cleaning process. Referred to as ‘auto-multithresh’, this algorithm is intended to mimic what an experienced user would do when manually masking images. It can be used by setting the usemask option in tclean to ’auto-multithresh’.

A full description of the auto-multithresh algorithm is given in the CASA docs and Kepley et al. 2020. Here we briefly describe the fundamental features of the algorithm. First, the algorithm identifies regions that are either above a signal-to-noise limit or a sidelobe level (as defined by the user), whichever is higher. If the mask regions are smaller than some fraction of the beam, they are removed (“pruned”) to avoid including spurious noise peaks in the masked region. If there is extended low signal to noise emission surrounding the initial mask then this emission is included in the mask by expanding the noise/sidelobe mask into the low signal-to-noise regions using a process called binary dilation. Finally, the mask is expanded to include a buffer region around the masked emission by convolving the mask by a Gaussian and then retaining only regions that are above n% of the peak in the smoothed mask. Absorption regions can also be masked using a method similar to that for the initial threshold mask, but are not pruned or extended into low signal-to noise regions. Note that for cubes, the algorithm masks each channel independently. The algorithm operates on the residual image at the beginning of every minor cycle, so that the mask updates as the clean progresses.

As of CASA 5.6, auto-multithresh now functions with polarization data. It applies the same algorithms to the Stokes QUV images as used for the Stokes I image. This means that the full masking process is applied to the positive emission (including the prune and grow steps), but that the masking of the negative emission only includes the initial threshold mask (no prune or grow).

Parameters

The parameters for the auto-multithresh algorithm are defined in terms of fundamental image properties (e.g., fraction of beam, signal-to-noise, etc) rather than exact numerical values (e.g., 3”, 6 mJy/beam). This feature allows the same set of parameters to be used for multiple images. The auto-multithresh algorithm calculates the numerical value of each parameter internally by estimating the noise in the residual image (currently done using the median absolute deviation) and the beam size and sidelobe level calculated by tclean internally. The numerical values are recalculated at the beginning of every minor cycle and may change as the clean progresses.

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

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

  1. pbmask (type: double) -- primary beam cutoff
  2. smoothfactor (type: double) -- controls how much to smooth the initial mask. We do not recommend changing this from its default value.
  3. 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.
  4. 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.
  5. 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.
  6. 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.
  7. verbose (type: bool) -- turns on/off per-channel logging information on the masks (default = False)

Additionally, the following parameter should be considered when using automasking (this not a subparameter of auto-multithresh).

  1. fastnoise (type: bool) -- If True, estimate the noise using the median absolute deviation. If False, estimate the noise using the chauvenet (no mask) or the median absolute deviations of pixels outside the mask but inside the primary beam mask. The latter estimation may be more accurate when emission covers a significant portion of the field of view (default=True).

The default automasking parameters in tclean have been set to values appropriate for the ALMA 12m array in its more compact configurations (i.e. not long baseline data). These values depend on the PSF (i.e., uv-coverage) of the data. Data with PSFs significantly different than the fiducial ALMA 12m arrays (e.g. 7m data) will likely need the auto-multithresh parameters to be modified in order to achieve the best mask. The default pipeline parameters presented here (Table below) have been extensively tested and most likely will work for your data set, although you may wish to optimize them for your particular imaging case. Data from other instruments (e.g., VLA, ATCA, etc) has been shown to work with this algorithm in limited testing, but the parameter values may need to be optimized.

Table of Standard Values

The values presented in this table are the standard values that the pipeline uses when it performs automasking and have been extensively tested. The values listed for combined data in the last column are tentative and should serve as a starting point for your analysis. This table can be found in the Pipeline User's Guide section 9.39 (hif_makeimages).

Parameter 7m (continnum/line) 12m (b75 < 300m) 12m (300m < b75 < 400m) 12m (b75 > 400m) 12m + 7m combined TENTATIVE
noisethreshold 5.0 4.25 5.0 5.0 4.25
sidelobethreshold 1.25 2.0 2.0 2.5 2.0
lownoisethreshold 2.0 1.5 1.5 1.5 1.5
minbeamfrac 0.1 0.3 0.3 0.3 0.3
negativethreshold 0.0 0.0 (continuum)
15.0 (line) 		0.0 (continuum)
7.0 (line) 		0.0 (continuum)
7.0 (line) 		0.0
fastnoise False False False True False

The term “b75” refers to the 75th percentile of baselines. You can find this value in the pipeline weblog by going to the antenna configuration page and clicking on the “baselines” tab, the table lists the baseline length in increasing distance and the percentile. If you don’t have a weblog for your data, you may use the Analysis Utilities task au.getBaselineStats. You can obtain the Analysis Utilities tasks by following the Analysis Utilities CASA Guide.

The 75% baseline split is only valid for ALMA data. For data from other telescopes, evaluate the uv-coverage of the observation. If the uv-coverage of a data set is good, then the PSF will have relatively low sidelobes and the default 12m parameters are likely a good start. We have had success in limited testing using the short baseline ALMA parameters for non-snapshot VLA observations. If the uv-coverage of a data set is poor, then the PSF will have higher sidelobes and the 7m array parameters are probably a better start. An example of this case would be snapshot VLA observations or ATCA observations.

If you are using pipeline imaging tasks, these automasking parameters are also able to be defined within those task calls. See the Pipeline Tasks Reference Manual section 2.12.

Guidelines for Optimizing auto-multithresh Parameters

  • Start with the pipeline parameters for the data with the most similar parameters to your data set. When in doubt, choose the default short baseline 12m parameters.
  • For cubes, use a small (~10 channel) cube with representative emission as a test case to avoid long run times.
  • Change one parameter at a time and keep track of your previous images so that you can directly compare different runs.
  • Inspect the log output. The ‘auto-multithresh’ algorithm reports extensive information to the logger on what it is doing including what thresholds are being used, how many regions are being pruned, etc. If using CASA 5.3+, set verbose=True to get this information.
  • Save the masks for each major cycle. To save the intermediate masks, type the following on the casa command line:
#In CASA
import os
os.environ['SAVE_ALL_AUTOMASKS']="true"

This feature does NOT work when tclean is run with parallel=True on cubes.

  • You can also inspect the intermediate masks by setting interactive=True. However, the resulting masks will NOT necessarily be the same as masks generated with interactive=False because how often the major cycle is triggered can be different in the two cases. Set cycleniter=niter when using non-interactive mode to get identical behavior.

Exploring the auto-multithresh Parameters

This section should serve as a guide to help you gain a more intuitive understanding of how the automasking algorithm works and understand better what each parameter does during the automasking process. Additionally, we will show you what can happen when you set a parameter to a non-optimal value and end up masking noise features or not masking real emission.

Obtain the Data

To demonstrate the use of the automasking parameters, we will be using the TW Hydra dataset which is a small dataset used in the “First Look” CASA Guides. Specifically, we will use a flagged, calibrated, self-calibrated, and continuum subtracted Band 7 Measurement Set that is created as part of the First_Look_at_Line_Imaging guide. Download and untar: twhya_selfcal.ms.contsub.tar (~300 MB). 

#In bash
tar xvf twhya_selfcal.ms.contsub.tar

Almost all of the values used in our following tclean commands remain constant throughout this section, and are explained in and taken from the First_Look_at_Line_Imaging guide. The main difference is that we are using the usemask=’auto-multithresh’ parameter to use the automasking feature and will be adjusting the auto-masking subparameters to highlight their use.

Make a Dirty Cube

First we make a dirty cube to figure out the value for threshold and get an idea of what the emission/noise looks like. It's always a good idea to perform a listobs first to get the basic information of the dataset. 

#In CASA
listobs('twhya_selfcal.ms.contsub', listfile='twhya_selfcal.ms.contsub.listobs.txt')

   Observer: cqi     Project: uid://A002/X327408/X6f  
Observation: ALMA
Data records: 44772       Total elapsed time = 4268.11 seconds
   Observed from   19-Nov-2012/07:56:23.5   to   19-Nov-2012/09:07:31.6 (UTC)

   ObservationID = 0         ArrayID = 0
  Date        Timerange (UTC)          Scan  FldId FieldName             nRows     SpwIds   Average Interval(s)    ScanIntent
  19-Nov-2012/07:56:23.5 - 08:02:11.3    12      5 TW Hya                    7616  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:08:09.6 - 08:13:57.3    16      5 TW Hya                    8442  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:19:53.9 - 08:25:41.7    20      5 TW Hya                    8389  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:32:00.5 - 08:37:48.2    24      5 TW Hya                    8409  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              08:43:45.6 - 08:49:33.4    28      5 TW Hya                    8514  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
              09:05:15.6 - 09:07:31.6    36      5 TW Hya                    3402  [0]  [6.05] [OBSERVE_TARGET#ON_SOURCE]
           (nRows = Total number of rows per scan) 
Fields: 1
  ID   Code Name                RA               Decl           Epoch   SrcId      nRows
  5    none TW Hya              11:01:51.796000 -34.42.17.36600 J2000   4          44772
Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)
  SpwID  Name                           #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz) BBC Num  Corrs  
  0      ALMA_RB_07#BB_2#SW-01#FULL_RES    384   TOPO  372533.086       610.352    234375.0 372649.9688        2  XX  YY

With the information above and from what we know of TW Hydra from the First_Look_at_Line_Imaging guide, we get the basic information for our tclean call to create the dirty cube (set niter=0 to perform no cleaning).

Figure 1: In CARTA, use a region in a line-free channel to determine RMS, then use a multiple of this (~3x) for the threshold parameter.
#In CASA
os.system('rm -rf twhya_dirtycube.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_dirtycube',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=0)

Determine Clean Threshold

In this tutorial we will clean non-interactively, and thus need to set a threshold to avoid over-cleaning emission. This is true for all non-interactive cleaning, and also a good idea for interactive cleaning (the threshold parameter is not a subparameter of auto-multithresh). After tclean has finished, you should now open the dirty image in CARTA. If using NRAO machines, you can open a new terminal tab, cd to the working directory, then type: 

#In bash
carta --no_browser

Copy the output URL into a browser to view your CARTA session. Select and load twhya_dirtycube.image. Selecting a region in a line free channel (e.g., channel 3) we find that the RMS is ~29 mJy/beam. We choose the threshold to be 3x the RMS: threshold='87mJy' (The "/beam" is implied).


Base Parameters

To figure out what automasking parameter values we should use for this data, we run au.getBaselineStats which will tell us the "b75" value for our data. If you are not using an NRAO machine, follow the Analysis_Utilities guide. 

#In CASA
import analysisUtils as au
au.getBaselineStats('twhya_selfcal.ms.contsub')

Found 1 obsIDs
Found 325 baselines (26 stations)
Unprojected lengths:  min=15.063874, max=374.719747, rms=162.602204 m
number=325, min=15.06m, max=374.72m, median=139.34m, mean=142.78m, std=77.81m
20%ile=69.4m 25%ile=79.0m, 30%ile=87.6m, 75%ile=197.4m, 90%ile=250.2m

From the output above (which was printed to the terminal window), we see that b75 (indicated as 75%ile) is 197.4m, so we will use the pipeline default values in the "12m (b75 < 300m)" column of the #Table_of_Standard_Values above. This will give us a mask to compare our testing against which will be displayed as black contours in all of the following figures.

Figure 2: In CARTA, we display channel 7 of the cube created with the base parameters for 12m (b75 < 300m). On the left is the image, and on the right is the residual. The black contours mark the final mask created for this channel.
#In CASA
os.system('rm -rf twhya_base_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_base_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5, 
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)
  1. Open the following three images in CARTA (with File > Open and ctrl+click, or one at a time with File > Append):
    twhya_base_params.image
    twhya_base_params.residual
    twhya_base_params.mask
  2. In the Image List, lock all 3 images in XY (spatial matching), Z (spectral matching), and R (raster scaling matching). If .image is not already set as the reference for matching, right click the image name and select "Set as ... reference." To reorder images in the list, click an index, then click-and-drag the index.
  3. Click the settings gear of the Image Viewer, go to the Global tab, enable multi-panel, and choose 2 columns and 1 row. Adjust the images to an appropriate position (Zoom to fit is recommended).
  4. Use the Animator to page through the channels. Notice that the mask image is blank for all channels except 5 and 7. For this example we choose channel 7.
  5. Adjust the raster scaling in Render Configuration.
  6. Open the Contour dialog. Unlock data source. Choose the mask. The Histogram, Generator mode, and Parameters provide different ways to calculate Levels, however for a mask contour we can ignore these. A mask image has value 1.0 where the mask exists, and 0.0 elsewhere (see this by moving your cursor around a mask image). Enter 0.5 directly into the Levels box, then click apply. This displays a contour outline of the mask on all 3 images. Go to the Styling tab to adjust color and thickness.

With verbose=True set, automasking reports its progress per-channel in the logger window. Below is the logger output for the first major clean cycle using the base parameters above.

  • The minimum mask size and sidelobeLevel are calculated in order to define different thresholds based on the user input.
----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[twhya_base_params] Setting up an auto-mask within PB mask limit 
---------------------------------------------------- Run Automask  ---------------------------------------------
Set Deconvolution Options for [twhya_base_params] : hogbom
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
Generating AutoMask
[C0] Using Chauvenet algorithm for the image statistics 
[C1] Using Chauvenet algorithm for the image statistics 
[C2] Using Chauvenet algorithm for the image statistics 
[C3] Using Chauvenet algorithm for the image statistics 
[C4] Using Chauvenet algorithm for the image statistics 
[C5] Using Chauvenet algorithm for the image statistics 
[C6] Using Chauvenet algorithm for the image statistics 
[C7] Using Chauvenet algorithm for the image statistics 
[C8] Using Chauvenet algorithm for the image statistics 
[C9] Using Chauvenet algorithm for the image statistics 
[C10] Using Chauvenet algorithm for the image statistics 
[C11] Using Chauvenet algorithm for the image statistics 
[C12] Using Chauvenet algorithm for the image statistics 
[C13] Using Chauvenet algorithm for the image statistics 
[C14] Using Chauvenet algorithm for the image statistics 
SidelobeLevel = 0.21961
prune size=9.47266(minbeamfrac=0.3 * beampix=31.5755)
  • Masking is performed and then pruned. Only "error" messages are reported per channel (i.e. no regions removed or no regions found).
  • The mask is smoothed and then in later cycles grown into lower signal-to-noise regions and pruned again.
  • Negative emission is masked if present.
*** Start auto-multithresh processing for Channel 0***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 1***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 2***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 3***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 4***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 5***
 No regions are removed in pruning process.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 6***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 7***
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 8***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 9***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 10***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 11***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 12***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 13***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Start auto-multithresh processing for Channel 14***
 No regions are found in this plane.
No negative region was found by auotmask.
*** Timing summary for whole planes ***
Total time to create the initial threshold mask:  real 0.02s ( user 0.48s, system 0s)
Total time to prune the initial threshold mask: real 0.03s (user 0.66s, system 0.01s)
Total time to create the smoothed initial threshold mask: real 0.75s (user 2.13s, system 1.19s)
Total time to grow the previous mask: real 0s (user 0s, system 0s)
Total time to prune the grow mask: real 0s (user 0s, system 0s)
Total time to create the smoothed grow mask: real 0s (user 0s, system 0s)
  • Automasking evaluates whether or not to continue to try masking each channel. If no masks were found in a channel or all masks were removed through the pruning process, the channel is no longer considered in the upcoming minor clean cycles.
Stopping masking for chan=0
Stopping masking for chan=1
Stopping masking for chan=2
Stopping masking for chan=3
Stopping masking for chan=4
Stopping masking for chan=6
Stopping masking for chan=8
Stopping masking for chan=9
Stopping masking for chan=10
Stopping masking for chan=11
Stopping masking for chan=12
Stopping masking for chan=13
Stopping masking for chan=14
  • A summary of the automasking process is given in the "Automask summary" table. Here, you can see the threshold that was used to create the mask (either "noise" or "sidelobe"), the number of regions found, number of regions pruned, etc. Our automasking results in masks for channels 5 and 7.
========== automask summary ==========
chan masking? median   RMS         peak   thresh_type   thresh_value   N_reg N_pruned N_grow N_grow_pruned N_neg_pix
[C0] F        --  --  0.12525  noise  --        0  0  0  0  0
[C1] F        --  --  -0.113465  noise  --        0  0  0  0  0
[C2] F        --  --  0.136005  noise  --        1  1  0  0  0
[C3] F        --  --  0.148636  noise  --        1  1  0  0  0
[C4] F        --  --  0.119875  noise  --        0  0  0  0  0
[C5] T        0.000376255  0.0328259  0.190839  noise  0.139886        1  0  0  0  0
[C6] F        --  --  0.152505  noise  --        1  1  0  0  0
[C7] T        -3.61604e-05  0.0308196  0.166526  noise  0.130947        2  1  0  0  0
[C8] F        --  --  0.102411  noise  --        0  0  0  0  0
[C9] F        --  --  0.106988  noise  --        0  0  0  0  0
[C10] F        --  --  -0.117285  noise  --        0  0  0  0  0
[C11] F        --  --  0.109285  noise  --        0  0  0  0  0
[C12] F        --  --  0.108688  noise  --        0  0  0  0  0
[C13] F        --  --  -0.111106  noise  --        0  0  0  0  0
[C14] F        --  --  -0.106717  noise  --        0  0  0  0  0
========== END of automask summary ==========
  • Cleaning is performed on the masked image and the progress (peak residual, number of iterations, etc.) is reported per channel. Notice how minor cycles (iterations/iters) are performed on our masked channels 5 and 7 only.
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_base_params] Number of pixels in the clean mask : 362 out of a total of 937500 pixels. [ 0.0386133 % ]
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Minor Cycle controls : {'cycleniter': 100000, 'cyclethreshold': 0.08699999749660492, 'loopgain': 0.10000000149011612, 'nsigma': 0.0, 'thresholdreached': True}
    itsNsigma=0
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
Set Deconvolution Options for [twhya_base_params] : hogbom
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
    itsNsigma=0
[twhya_base_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.087, CycleNiter=100000, Gain=0.1
[twhya_base_params:C0] iters=0->0 [0], model=0->0, peakres=0.12525->0.12525, Skipped this plane. Zero mask.
[twhya_base_params:C1] iters=0->0 [0], model=0->0, peakres=0.113465->0.113465, Skipped this plane. Zero mask.
[twhya_base_params:C2] iters=0->0 [0], model=0->0, peakres=0.136005->0.136005, Skipped this plane. Zero mask.
[twhya_base_params:C3] iters=0->0 [0], model=0->0, peakres=0.148636->0.148636, Skipped this plane. Zero mask.
[twhya_base_params:C4] iters=0->0 [0], model=0->0, peakres=0.119875->0.119875, Skipped this plane. Zero mask.
[twhya_base_params:C5] iters=0->16 [16], model=0->0.190617, peakres=0.190839->0.085914, Reached cyclethreshold.
[twhya_base_params:C6] iters=16->16 [0], model=0->0, peakres=0.152505->0.152505, Skipped this plane. Zero mask.
[twhya_base_params:C7] iters=16->39 [23], model=0->0.258554, peakres=0.166526->0.0863384, Reached cyclethreshold.
[twhya_base_params:C8] iters=39->39 [0], model=0->0, peakres=0.102411->0.102411, Skipped this plane. Zero mask.
[twhya_base_params:C9] iters=39->39 [0], model=0->0, peakres=0.106988->0.106988, Skipped this plane. Zero mask.
[twhya_base_params:C10] iters=39->39 [0], model=0->0, peakres=0.117285->0.117285, Skipped this plane. Zero mask.
[twhya_base_params:C11] iters=39->39 [0], model=0->0, peakres=0.109285->0.109285, Skipped this plane. Zero mask.
[twhya_base_params:C12] iters=39->39 [0], model=0->0, peakres=0.108688->0.108688, Skipped this plane. Zero mask.
[twhya_base_params:C13] iters=39->39 [0], model=0->0, peakres=0.111106->0.111106, Skipped this plane. Zero mask.
[twhya_base_params:C14] iters=39->39 [0], model=0->0, peakres=0.106717->0.106717, Skipped this plane. Zero mask.
[twhya_base_params] Total model flux (over all planes) : 0.44917     Peak Residual (over all planes) : 0.0863384 in C7:P0
Completed 39 iterations.

Initial Mask Threshold Calculation

The initial mask threshold is calculated by the following method:

noiseThresholdValue = noisethreshold * rms in residual image
sidelobeThresholdValue = sidelobethreshold * sidelobeLevel * peak in residual image

where the sidelobeLevel is calculated by tclean and reported in the logger (you can see this in the example logger output above).

The initial threshold is max(noiseThresholdValue, sidelobeThresholdValue).

In practice, the sidelobethreshold often sets the limit in the first few major cycles because the sidelobe emission is often larger than the random noise in the image. Once the strongest emission has been removed from the residual image, later cycles tend to use the noisethreshold.

Turning Down noisethreshold

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

Figure 3: On the left is the image and on the right is the residual of channel 7 of 15 of the cube created with these parameters. The black contours represent the mask created when using the base parameters, and the white contours represent the masks created by setting a small noisethreshold value, which favors the use of sidelobethreshold.
#In CASA
os.system('rm -rf twhya_low_noise_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_noise_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=1.75, # lowered from base
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)

Open the following in CARTA:
twhya_low_noise_params.image
twhya_low_noise_params.residual
twhya_low_noise_params.mask
twhya_base_params.mask

Using the method described above, display the .image and .residual sibe-by-side, overlay the base mask as black contours, and overlay the lownoise mask as white contours.

By reducing noisethreshold, auto-multithresh is now using sidelobethreshold instead to draw the initial masks for most channels. You can see this in the logger. The "thresh_type" column of the "automask summary" reports which threshold is being used to find maskable emission in each channel. The output below is trimmed for brevity indicated by "*****". 

----------------------------------------------------------- Run Major Cycle 1 -------------------------------------
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
[twhya_lownoise_params] Setting up an auto-mask within PB mask limit 
---------------------------------------------------- Run Automask  ---------------------------------------------
Set Deconvolution Options for [twhya_lownoise_params] : hogbom
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
Generating AutoMask

*****

SidelobeLevel = 0.21961
prune size=9.47266(minbeamfrac=0.3 * beampix=31.5755)

*****

========== automask summary ==========
chan masking? median   RMS         peak   thresh_type   thresh_value   N_reg N_pruned N_grow N_grow_pruned N_neg_pix
[C0] T        0.000326007  0.0306058  0.12525  sidelobe  0.0553385        206  152  0  0  0
[C1] T        0.000123019  0.0282244  -0.113465  sidelobe  0.0499591        206  152  0  0  0
[C2] T        4.33802e-05  0.03036  0.136005  sidelobe  0.0597796        150  105  0  0  0
[C3] T        0.000242864  0.0281953  0.148636  sidelobe  0.065527        95  85  0  0  0
[C4] T        0.000411201  0.0290491  0.119875  sidelobe  0.0530629        210  156  0  0  0
[C5] T        0.000376255  0.0328259  0.190839  sidelobe  0.0841967        43  40  0  0  0
[C6] T        0.000135499  0.0357602  0.152505  sidelobe  0.0671187        126  82  0  0  0
[C7] T        -3.61604e-05  0.0308196  0.166526  sidelobe  0.0731053        56  40  0  0  0
[C8] T        0.000578858  0.0269325  0.102411  noise  0.0477108        233  178  0  0  0
[C9] T        5.75308e-05  0.0267095  0.106988  sidelobe  0.0470487        229  155  0  0  0
[C10] T        -3.76976e-05  0.0264451  -0.117285  sidelobe  0.0514762        189  154  0  0  0
[C11] T        5.7293e-06  0.0268588  0.109285  sidelobe  0.048006        212  150  0  0  0
[C12] T        -0.000308202  0.0278082  0.108688  noise  0.0483562        219  154  0  0  0
[C13] T        1.40047e-05  0.0268257  -0.111106  sidelobe  0.0488139        214  159  0  0  0
[C14] T        0.00029441  0.0257197  -0.106717  sidelobe  0.0471667        220  160  0  0  0
========== END of automask summary ==========
----------------------------------------------------------------------------------------------------------------------------------------
[twhya_lownoise_params] Number of pixels in the clean mask : 109074 out of a total of 937500 pixels. [ 11.6346 % ]
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Minor Cycle controls : {'cycleniter': 100000, 'cyclethreshold': 0.08699999749660492, 'loopgain': 0.10000000149011612, 'nsigma': 0.0, 'thresholdreached': True}
    itsNsigma=0
---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------
Set Deconvolution Options for [twhya_lownoise_params] : hogbom
 Absolute Peak residual within mask : 0.190839, over full image : 0.190839
Processing channels in range [0, 14]
    itsNsigma=0
[twhya_lownoise_params] Run Hogbom minor-cycle on 15 chans | CycleThreshold=0.087, CycleNiter=100000, Gain=0.1
[twhya_lownoise_params:C0] iters=0->30 [30], model=0->0.245732, peakres=0.12525->0.0869206, Reached cyclethreshold.
[twhya_lownoise_params:C1] iters=30->47 [17], model=0->0.0652784, peakres=0.105279->0.0865448, Reached cyclethreshold.
[twhya_lownoise_params:C2] iters=47->90 [43], model=0->0.310039, peakres=0.136005->0.0869305, Reached cyclethreshold.
[twhya_lownoise_params:C3] iters=90->109 [19], model=0->0.194673, peakres=0.148636->0.0853744, Reached cyclethreshold.
[twhya_lownoise_params:C4] iters=109->137 [28], model=0->0.162311, peakres=0.119875->0.0869172, Reached cyclethreshold.
[twhya_lownoise_params:C5] iters=137->160 [23], model=0->0.258899, peakres=0.190839->0.0869799, Reached cyclethreshold.
[twhya_lownoise_params:C6] iters=160->243 [83], model=0->0.851988, peakres=0.152505->0.0869064, Reached cyclethreshold.
[twhya_lownoise_params:C7] iters=243->288 [45], model=0->0.419446, peakres=0.166526->0.0862288, Reached cyclethreshold.
[twhya_lownoise_params:C8] iters=288->297 [9], model=0->0.00968076, peakres=0.102411->0.0869284, Reached cyclethreshold.
[twhya_lownoise_params:C9] iters=297->310 [13], model=0->0.0652534, peakres=0.106988->0.0856935, Reached cyclethreshold.
[twhya_lownoise_params:C10] iters=310->324 [14], model=0->0.0348655, peakres=0.117285->0.0865639, Reached cyclethreshold.
[twhya_lownoise_params:C11] iters=324->345 [21], model=0->0.105682, peakres=0.109285->0.086327, Reached cyclethreshold.
[twhya_lownoise_params:C12] iters=345->357 [12], model=0->0.0756125, peakres=0.108688->0.0868853, Reached cyclethreshold.
[twhya_lownoise_params:C13] iters=357->372 [15], model=0->0.10012, peakres=0.0996836->0.0860267, Reached cyclethreshold.
[twhya_lownoise_params:C14] iters=372->382 [10], model=0->0.0373442, peakres=0.106717->0.0866348, Reached cyclethreshold.
[twhya_lownoise_params] Total model flux (over all planes) : 2.93693     Peak Residual (over all planes) : 0.0869799 in C5:P0
Completed 382 iterations.

Notice that masking and cleaning occurs in all channels. Only channels 8 and 12 use thresh_type = noise. In later major cycles, we can see that more channels switch from sidelobe to noise. The threshold used by auto-multithresh are independent for each channel and can change between major cycles as clean progresses. In Figure 3, we see that fainter emission gets picked up by automasking when we reduce the noisethreshold value (as expected), but that noise spikes are clearly being picked up in both line and line-free channels. We can also see this clearly in the log as the number of masks created and reported in the automask summary is much, much higher than in the base parameters log output.

Turning Down sidelobethreshold

Now we turn sidelobethreshold down to see what impact this has on the masks that are created. Note that because of how masks are drawn during each step of the automasking process, we need to also turn down either noisethreshold or lownoisethreshold to guarantee that the sidelobethreshold is used by the algorithm during either the initial step or the secondary growing of the masks step.

Figure 4: On the left is the image and on the right is the residual of channel 7 of 15 of the cube created with these parameters. The black contours represent the mask created when using the base parameters and the white contours represent the mask created using a small sidelobethreshold value.
#In CASA
os.system('rm -rf twhya_low_sidelobe_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_sidelobe_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=1.0, # lowered from base
lownoisethreshold=0.5, # lowered from base
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)

In Figure 4 we see that turning down sidelobethreshold makes a mask that covers a much larger region and clearly covers areas of the image that do not contain real emission.

The Effects of Pruning

Pruning is used to remove small noise-like regions from the resulting mask. Note that regions are pruned before being convolved with a Gaussian to create a buffer around the mask.

Turning Off Pruning

Now we use the base parameters but turn off pruning by setting minbeamfrac to 0.0 so that no regions are removed from the mask.

Figure 5: 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 white contours represent the mask created with no pruning.
#In CASA
os.system('rm -rf twhya_low_minbeamfrac_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_minbeamfrac_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=0.0, # lowered from base
negativethreshold=15.0,
verbose=True,
fastnoise=False)

In Figure 5 we see that with no pruning, a noise spike near the edge of channel 7 was picked up, kept, and cleaned. Notice that several line-free channels include masks around noise spikes as well.

Increasing The Pruning

Next, we set minbeamfrac to 2.0 and leave the rest of parameters at their base values, so that all masks smaller than twice the beam will get pruned. 

#In CASA
os.system('rm -rf twhya_high_minbeamfrac_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_high_minbeamfrac_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5,
minbeamfrac=2.0, # raised from base
negativethreshold=15.0,
verbose=True,
fastnoise=False)

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

SidelobeLevel = 0.21961
prune size=63.151(minbeamfrac=2 * beampix=31.5755)
*****
No regions found by automasking
15 of 15 channels used the noise threshold to mask, but no emission was found. Try decreasing your noisethreshold parameter if you want to capture emission in these channels.
5 of 15 channels had all regions removed by pruning. Try decreasing minbeamfrac to remove fewer regions
*****
[twhya_high_minbeamfrac_params] Number of pixels in the clean mask : 0 out of a total of 937500 pixels. [ 0 % ]
*****
Reached global stopping criterion : zero mask
*****
Restoring with an empty model image. Only residuals will be processed to form the output restored image.

In this case, no cleaning is done, a dirty cube is generated, and tclean is exited.

Note that pruning does not calculate whether the emission within a region is significant, so you can inadvertently prune point sources with high signal-to-noise emission (e.g., calibrators). In this case, reduce the minbeamfrac.

Growing The Mask Into Low Signal-to-Noise Regions

The lownoisethreshold parameter is used to define a low signal-to-noise mask. The mask from the previous cycle is then grown into the region defined by the low signal-to-noise mask via binary dilation. To protect against masking sidelobes with the low signal-to-noise mask, the threshold for this mask is determined using the following heuristic

lowNoiseThresholdValue = lownoisethreshold * rms in residual image
max(sidelobeThresholdValue, lowNoiseThresholdValue).

If the value of the lownoisethreshold is lower than the value of the sidelobethreshold, then the mask will not be expanded into low signal-to-noise regions.

Finally, we lower the lownoisethreshold value to try and grow the mask into regions of lower signal-to-noise during the binary dilation stage.

Figure 6: On the left is the channel 7 mask created with base parameters, and on the right is the channel 7 mask created by lowering lownoisethreshold only. There is no difference, because the default value of sidelobethreshold takes over.
#In CASA
os.system('rm -rf twhya_low_lownoise_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_lownoise_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.0, # lowered from base
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)

This mask is the same as the base mask even if we try to make the lownoisethreshold very, very small (i.e. something like 0.0001). This is due to the fact that if the value of the lownoisethreshold is less than the value of the sidelobethreshold then the sidelobethreshold is used to set mask. This effectively attempts to act as a “fail-safe” to ensure that overmasking does not occur accidentally by setting one of the thresholds very small. For demonstration purposes, we can force it to use a small lownoisethreshold by also setting the sidelobethreshold very small.

Figure 7: On the left is the image and on the right is the residual of channel 7 of 15 of the cube created with these parameters. The black contours represent the mask created when using the base parameters, and the white contours represent the mask created using a small lownoisethreshold value (with an even lower sidelobethreshold so that it doesn't become the limiter).
#In CASA
os.system('rm -rf twhya_low_lownoise+sidelobe_params.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_low_lownoise+sidelobe_params',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=0.5, # lowered from base
lownoisethreshold=1.0, # lowered from base
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False)

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

Advanced Use Case - Merging User Masks with Automasking

In some situations, you may find that automasking is not masking obvious emission even after adjusting the thresholds. This most commonly occurs in cases:

  • with extremely bright, diffuse emission throughout the entire field of view
  • where there is low-level absorption that has a similar depth to any negative bowls in the same channel
  • where the cube has not been continuum subtracted.

The first case is due to the difficulty of estimating noise in the presence of wide-spread signal. The last two are cases when the astronomer “knows” where the real signal is (likely through some additional information), but where the algorithm does not have enough knowledge to distinguish.

In these cases, it may be helpful to “kick start” automasking by first defining masks by hand for a round of cleaning and then having automasking pick up from there. To continue cleaning an existing image, we simply use the restart=True option within tclean (which is the default value).

First, we perform a round of cleaning by hand by setting interactive=True so that tclean will let us manually create masks before each major cycle. As mentioned above in the “Make a Dirty Cube” section, in order to get identical cleaning behavior between interactive and non-interactive modes, you must set cycleniter equal to niter. With this arbitrary very large number (cycleniter=niter=100000), the minor cycle iteration limit per major cycle will not be reached. Instead, tclean will continue cleaning until the internally computed cyclethreshold or the user specified threshold is hit. One last parameter to set is usemask='user' to indicate that you would like to manually mask emission. Below is an example tclean command using the same dataset as above.

Figure 8: The tclean interactive GUI, with a manually drawn mask for channel 7 shown in white.
#In CASA
os.system('rm -rf twhya_manual+automask.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_manual+automask',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=True, # manual masking in interactive GUI
niter=100000,
cycleniter=100000, # set equal to niter to have interactive mode replicate the behavior of non-interactive mode
threshold='87mJy',
usemask='user' # default value, but shown here explicitly in place of auto-multithresh
)

Once you are satisfied that you have masked all the obvious emission and cleaned to a level where automasking will start to pick up extended emission, you can exit tclean by pressing the red “X”. After tclean finishes, you can start tclean again but now in non-interactive mode (interactive=False), with automasking enabled (usemask='auto-multithresh'), and using the appropriate automasking thresholds that match your data (see the “Table of Standard Values” above). Finally, we set restart=True which will make tclean pick up where it left off using the existing .mask file, along with the existing .model, .psf, .image, etc. that share the same base imagename.

Side tip: if you already have a CASA region file from an earlier clean or one you have created, and is has a different name, you can read it into tclean using the mask parameter, mask='basename.mask'. The region file method can also be used non-interactively on it's own.

You may see the following warnings. This is a known issue, and can be ignored for the above command. 

 	----------------------------------------------------------- Run Major Cycle 2 -------------------------------------
2024-04-02 23:11:40	WARN	task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1572)	Error : Error in copying internal T/F mask : Error in deleting internal T/F mask : Region mask0 could not be removed
2024-04-02 23:11:40	WARN	task_tclean::SynthesisImagerVi2::runCubeGridding (file src/code/synthesis/ImagerObjects/SynthesisImagerVi2.cc, line 1572)+	Cannot delete the mask (used in another process)

However, this does prevent restart=True from working correctly. There is a known issue with file locks in tclean. To use restart=True, exit CASA at this step and start a new session before running the command below.

Figure 9: On the left is the image and on the right is the residual of channel 7 of the cube created with manual + automasking. The black contours represent the mask created when using automasking only with the base parameters, and the white contours represent the mask created by first doing an interactive clean and then an automasking non-interactive clean.
#In CASA
# os.system('rm -rf twhya_manual+automask.*')
tclean(vis='twhya_selfcal.ms.contsub',
imagename='twhya_manual+automask',
field='5',
spw='0',
specmode='cube',
nchan=15,
start='0.0km/s',
width='0.5km/s',
outframe='LSRK',
restfreq='372.67249GHz',
deconvolver='hogbom',
gridder='standard',
imsize=[250,250],
cell='0.1arcsec',
weighting='briggsbwtaper',
robust=0.5,
restoringbeam='common',
interactive=False,
niter=100000,
threshold='87mJy',
# Automasking parameters below this line
usemask='auto-multithresh',
noisethreshold=4.25,
sidelobethreshold=2.0,
lownoisethreshold=1.5, 
minbeamfrac=0.3,
negativethreshold=15.0,
verbose=True,
fastnoise=False,
restart=True # default value, but shown here explicitly
)

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

Version History

The older version of this guide is here: Automasking Guide CASA 5

Retrieved from "https://casaguides.nrao.edu/index.php?title=Automasking_Guide_CASA_6.5.4&oldid=37429"